Enable and fix the TGO CaSSIS driver

CHANGELOG.md

@@ -36,6 +36,7 @@ release.
 ## [Unreleased]
 
 ### Added
+- Re-enabled and fixed the TGO CaSSIS driver, which now emits the CaSSIS rational distortion. Validated against ISIS to within ~0.013 pixel across 130 framelets of two real stereo pairs. [#720](https://github.com/DOI-USGS/ale/pull/720)
 - `MroHiRisePds3LabelNaifSpiceDriver`, a PDS3 EDR label driver for HiRISE that generates an ISD directly from a raw EDR label without requiring an ISIS cube, paralleling the existing CTX PDS3 driver. [#702](https://github.com/DOI-USGS/ale/pull/702)
 - Added a catch to try correcting paths in metakernels (using spice_root) if they have been left as default. [#703](https://github.com/DOI-USGS/ale/pull/703)
 

ale/base/type_distortion.py

@@ -64,6 +64,34 @@ def usgscsm_distortion_model(self):
         """
         return {"radial": {"coefficients": [0.0, 0.0, 0.0]}}
 
+
+class CassisDistortion():
+    """
+    Mix-in for the TGO CaSSIS rational (ratio-of-quadratics) distortion model
+    by Tulyakov/Ivanov (EPFL), the same model ISIS uses in TgoCassisDistortionMap.
+    With chi = [x^2, x*y, y^2, x, y, 1] the two directions are
+        corrected = (A1_corr.chi)/(A3_corr.chi), ...    (distorted -> undistorted)
+        distorted = (A1_dist.chi)/(A3_dist.chi), ...    (undistorted -> distorted)
+    each using three 6-vectors. The 36 coefficients live in the IK/IAK as
+    INS<ikid>_OD_A{1,2,3}_{CORR,DIST}.
+    """
+
+    @property
+    def usgscsm_distortion_model(self):
+        """
+        Returns
+        -------
+        : dict
+          Dictionary with the 36 CaSSIS distortion coefficients, packed in the
+          order A1_corr, A2_corr, A3_corr, A1_dist, A2_dist, A3_dist (6 each),
+          matching the unpacking in USGSCSM Distortion.cpp.
+        """
+        coefficients = []
+        for vec in ["A1_CORR", "A2_CORR", "A3_CORR", "A1_DIST", "A2_DIST", "A3_DIST"]:
+            coefficients.extend(self.naif_keywords["INS{}_OD_{}".format(self.ikid, vec)])
+        return {"cassis": {"coefficients": coefficients}}
+
+
 class KaguyaSeleneDistortion():
     """
     Mix-in for sensors on the Kaguya/Selene mission. 

ale/drivers/__init__.py

@@ -25,7 +25,7 @@
 
 # Explicit list of disabled drivers
 
-__disabled_drivers__ = ["tgo_drivers", "osirisrex_drivers"]
+__disabled_drivers__ = ["osirisrex_drivers"]
 
 # dynamically load drivers
 __all__ = [os.path.splitext(os.path.basename(d))[0] for d in glob(os.path.join(os.path.dirname(__file__), '*_drivers.py'))]

ale/drivers/tgo_drivers.py

@@ -1,12 +1,15 @@
+import numpy as np
+import scipy.constants
+import spiceypy as spice
 from pyspiceql import pyspiceql
 
 from ale.base import Driver, WrongInstrumentException
 from ale.base.data_naif import NaifSpice
 from ale.base.label_isis import IsisLabel
 from ale.base.type_sensor import Framer
-from ale.base.type_distortion import NoDistortion
+from ale.base.type_distortion import CassisDistortion
 
-class TGOCassisIsisLabelNaifSpiceDriver(Framer, IsisLabel, NaifSpice, NoDistortion, Driver):
+class TGOCassisIsisLabelNaifSpiceDriver(Framer, IsisLabel, NaifSpice, CassisDistortion, Driver):
     """
     Driver for reading TGO Cassis ISIS3 Labels. These are Labels that have been ingested
     into ISIS from PDS EDR images but have not been spiceinit'd yet.
@@ -48,7 +51,7 @@ def ephemeris_start_time(self):
           ephemeris start time of the image.
         """
         if not hasattr(self, "_ephemeris_start_time"):
-            self._ephemeris_start_time = pyspiceql.utcToEt(utc=self.utc_start_time.strftime("%Y-%m-%d %H:%M:%S.%f"), searchKernels=self.search_kernels, useWeb=self.use_web)
+            self._ephemeris_start_time = pyspiceql.utcToEt(utc=self.utc_start_time.strftime("%Y-%m-%d %H:%M:%S.%f"), searchKernels=self.search_kernels, useWeb=self.use_web)[0]
         return self._ephemeris_start_time
 
     @property
@@ -68,3 +71,153 @@ def sensor_model_version(self):
     @property
     def sensor_name(self):
         return self.label['IsisCube']['Instrument']['SpacecraftName']
+
+    @property
+    def sample_summing(self):
+        """
+        CaSSIS stores SummingMode as an enum (0 = 1x1, 1 = 2x2, 2 = 4x4), not as
+        the summing factor itself. ISIS converts it as summing = sumMode * 2, then
+        falls back to 1 when that is 0 (see TgoCassisCamera). Replicate that here,
+        otherwise the CSM detector summing becomes 0 and groundToImage diverges.
+        """
+        sum_mode = self.label['IsisCube']['Instrument']['SummingMode']
+        summing = sum_mode * 2
+        if summing <= 0:
+            summing = 1
+        return summing
+
+    @property
+    def line_summing(self):
+        return self.sample_summing
+
+    @property
+    def detector_center_sample(self):
+        """
+        ISIS uses 0.5-based CCD coordinates (pixel centers at half integers),
+        so convert the IK boresight sample to the CSM 0-based convention by
+        subtracting 0.5, as the LRO, MRO, Dawn, MESSENGER, MEX, Kaguya and KPLO
+        drivers do. Without this the CSM look is offset from ISIS by half a pixel
+        in sample (and half in line), i.e. sqrt(0.5^2+0.5^2) ~ 0.707 px.
+        """
+        return super().detector_center_sample - 0.5
+
+    @property
+    def detector_center_line(self):
+        """
+        ISIS uses 0.5-based CCD coordinates; convert to the CSM 0-based
+        convention by subtracting 0.5 (see detector_center_sample).
+        """
+        return super().detector_center_line - 0.5
+
+    @property
+    def sensor_position(self):
+        """
+        CaSSIS sets LT_SURFACE_CORRECT with LIGHTTIME_CORRECTION=LT+S, so ISIS
+        applies the surface light-time correction. The shared
+        NaifSpice.sensor_position samples the target body at the raw ephemeris
+        time in that branch, but ISIS samples it at the surface-light-time
+        adjusted time (ephem - obs_tar_lt + radius_lt); the body moves along its
+        orbit during that interval, which otherwise leaves a constant
+        tens-of-meters camera-center bias versus ISIS. This override applies that
+        for CaSSIS only, so the shared path is unchanged for every other sensor.
+        It is the shared surface-light-time branch with the single change that the
+        body is sampled at adjusted_time. CaSSIS is a single-record framer, so
+        sampling the body at adjusted_time[0]..[-1] is exact.
+        """
+        if not (self.correct_lt_to_surface
+                and self.light_time_correction.upper() == 'LT+S'):
+            return super().sensor_position
+
+        if not hasattr(self, '_position'):
+            ephem = self.ephemeris_time
+            pos = []
+            vel = []
+
+            target = self.spacecraft_name
+            observer = self.target_name
+            if self.swap_observer_target:
+                target = self.target_name
+                observer = self.spacecraft_name
+
+            ephem_kwargs = {"startEt": ephem[0],
+                            "stopEt": ephem[-1],
+                            "numRecords": len(ephem),
+                            "ckQualities": ["reconstructed"],
+                            "spkQualities": ["reconstructed"],
+                            "searchKernels": self.search_kernels,
+                            "useWeb": self.use_web}
+
+            obs_tars_kwargs = {**ephem_kwargs,
+                               "target": target,
+                               "observer": observer,
+                               "frame": "J2000",
+                               "abcorr": self.light_time_correction,
+                               "mission": self.spiceql_mission}
+            ssb_obs_kwargs = {**ephem_kwargs,
+                              "target": observer,
+                              "observer": "SSB",
+                              "frame": "J2000",
+                              "abcorr": "NONE",
+                              "mission": self.spiceql_mission}
+
+            obs_tars = pyspiceql.getTargetStatesRanged(**obs_tars_kwargs)[0]
+            ssb_obs = pyspiceql.getTargetStatesRanged(**ssb_obs_kwargs)[0]
+
+            obs_tar_lts = np.array(obs_tars)[:, -1]
+            ssb_obs_states = np.array(ssb_obs)[:, 0:6]
+
+            radius_lt = (self.target_body_radii[2] + self.target_body_radii[0]) / 2 \
+                / (scipy.constants.c / 1000.0)
+            adjusted_time = ephem - obs_tar_lts + radius_lt
+
+            # The only change from the shared method: sample the target body at
+            # the surface-light-time adjusted time rather than the raw ephem time.
+            ssb_tars_kwargs = {**ephem_kwargs,
+                               "target": target,
+                               "observer": "SSB",
+                               "frame": "J2000",
+                               "abcorr": "NONE",
+                               "mission": self.spiceql_mission}
+            ssb_tars_kwargs["startEt"] = adjusted_time[0]
+            ssb_tars_kwargs["stopEt"] = adjusted_time[-1]
+            ssb_tars = pyspiceql.getTargetStatesRanged(**ssb_tars_kwargs)[0]
+            ssb_tar_states = np.array(ssb_tars)[:, 0:6]
+
+            _states = ssb_tar_states - ssb_obs_states
+
+            reference_frame_id = pyspiceql.translateNameToCode(frame=self.reference_frame,
+                                                               mission=self.spiceql_mission,
+                                                               searchKernels=self.search_kernels,
+                                                               useWeb=self.use_web)[0]
+
+            function_args = {**ephem_kwargs,
+                             "toFrame": reference_frame_id,
+                             "refFrame": 1,
+                             "mission": self.spiceql_mission}
+            function_args.pop("spkQualities")
+            rotations = pyspiceql.getTargetOrientationsRanged(**function_args)[0]
+
+            states = []
+            for i, rotation in enumerate(rotations):
+                quaternion = rotation[:4]
+                av = [0, 0, 0]
+                if len(rotation) > 4:
+                    av = rotation[4:]
+                rotation_matrix = spice.q2m(quaternion)
+                matrix = spice.rav2xf(rotation_matrix, av)
+                rotated_state = spice.mxvg(matrix, _states[i])
+                states.append(rotated_state)
+
+            for state in states:
+                if self.swap_observer_target:
+                    pos.append(-state[:3])
+                    vel.append(-state[3:])
+                else:
+                    pos.append(state[:3])
+                    vel.append(state[3:])
+
+            # SPICE works in km, so convert to m.
+            self._position = 1000 * np.asarray(pos)
+            self._velocity = 1000 * np.asarray(vel)
+            self._ephem = ephem
+        return self._position, self._velocity, self._ephem

environment.yml

@@ -13,7 +13,7 @@ dependencies:
   - networkx
   - nlohmann_json
   - numpy
-  - spiceql>=1.3.0
+  - spiceql>=1.3.0,<1.4
   - pvl
   - pytest
   - pytest-cov

include/ale/Distortion.h

@@ -11,7 +11,8 @@ namespace ale {
     CAHVOR,
     LUNARORBITER,
     RADTAN,
-    KPLOSHADOWCAM
+    KPLOSHADOWCAM,
+    CASSIS
   };
 }
 

src/Util.cpp

@@ -366,6 +366,8 @@ DistortionType getDistortionModel(json isd) {
       return DistortionType::RADTAN;
     } else if (distortion.compare("kplo_shadowcam") == 0) {
       return DistortionType::KPLOSHADOWCAM;
+    } else if (distortion.compare("cassis") == 0) {
+      return DistortionType::CASSIS;
     }
   } catch (...) {
     throw std::runtime_error("Could not parse the distortion model.");
@@ -559,6 +561,19 @@ std::vector<double> getDistortionCoeffs(json isd) {
         coefficients = std::vector<double>(1, 0.0);
       }
     } break;
+    case DistortionType::CASSIS: {
+      try {
+        coefficients = isd.at("optical_distortion")
+                          .at("cassis")
+                          .at("coefficients")
+                          .get<std::vector<double>>();
+        return coefficients;
+      } catch (...) {
+        throw std::runtime_error(
+            "Could not parse the cassis distortion model coefficients.");
+        coefficients = std::vector<double>(36, 0.0);
+      }
+    } break;
   }
   throw std::runtime_error(
       "Could not parse the distortion model coefficients.");