Add a CaSSIS rational distortion type

.gitmodules

@@ -6,7 +6,7 @@
 	url = https://github.com/USGS-Astrogeology/csm
 [submodule "ale"]
 	path = ale
-	url = git@github.com:DOI-USGS/ale.git
+	url = https://github.com/oleg-alexandrov/ale.git
 [submodule "json"]
 	path = json
 	url = https://github.com/nlohmann/json.git

CHANGELOG.md

@@ -35,6 +35,9 @@ release.
 
 ## [Unreleased]
 
+### Added
+- A `CASSIS` distortion type implementing the TGO CaSSIS rational ratio-of-quadratics distortion model, matching ISIS `TgoCassisDistortionMap`, with the json-name, ALE integer-enum, and coefficient-extraction dispatch wired in. Pairs with the ALE TGO CaSSIS driver. [#512](https://github.com/DOI-USGS/usgscsm/pull/512)
+
 ## [2.1.0] - 2026-06-09
 
 ### Added

include/usgscsm/Distortion.h

@@ -8,7 +8,7 @@
 
 // This should be synched with the enum in ale/Distortion.h
 enum DistortionType { RADIAL, TRANSVERSE, KAGUYALISM, DAWNFC, LROLROCNAC, CAHVOR,
-                     LUNARORBITER, RADTAN, KPLOSHADOWCAM };
+                     LUNARORBITER, RADTAN, KPLOSHADOWCAM, CASSIS };
 
 // Transverse distortion Jacobian
 void transverseDistortionJacobian(double x, double y, double *jacobian,

src/Distortion.cpp

@@ -4,6 +4,20 @@
 #include <Error.h>
 #include <string>
 
+// TGO CaSSIS rational distortion helper. The model (Tulyakov/Ivanov, EPFL) is a
+// ratio of quadratics with basis chi = [x^2, x*y, y^2, x, y, 1]. Each direction
+// uses three 6-vectors A1, A2, A3, so x_out = (A1.chi)/(A3.chi), y_out =
+// (A2.chi)/(A3.chi). Ported from ISIS TgoCassisDistortionMap. A points at the
+// first of 6 contiguous coefficients.
+static double cassisChiDotA(double x, double y, const double *A) {
+  return A[0] * x * x + A[1] * x * y + A[2] * y * y + A[3] * x + A[4] * y + A[5];
+}
+
+// CaSSIS CCD is 2048 px at a 0.01 mm pitch. The rational model is only valid on
+// the CCD; outside the box (and where the A3 denominator has roots) ISIS returns
+// the input unchanged. Match that guard here. Bound = 0.5*pitch*width + 0.2 mm.
+static const double cassisFocalBound = 0.5 * 0.01 * 2048.0 + 0.2;
+
 /**
  * @description Calculates the Jacobian matrix of the distortion function at a specific point (x, y). 
  * The distortion model is based on the given optical distortion coefficients.
@@ -420,6 +434,35 @@ void removeDistortion(double dx, double dy, double &ux, double &uy,
       return;
     } break;
 
+    // TGO CaSSIS rational model, distorted -> undistorted (CORR direction).
+    // Ported from ISIS TgoCassisDistortionMap::SetFocalPlane. Coefficients are
+    // 36 doubles: A1_corr, A2_corr, A3_corr, A1_dist, A2_dist, A3_dist (6 each).
+    // Here we use the CORR triple (indices 0..17).
+    case CASSIS: {
+      if (opticalDistCoeffs.size() != 36) {
+        csm::Error::ErrorType errorType = csm::Error::INDEX_OUT_OF_RANGE;
+        std::string message =
+            "Distortion coefficients for CaSSIS must be of size 36, "
+            "current size: " +
+            std::to_string(opticalDistCoeffs.size());
+        std::string function = "removeDistortion";
+        throw csm::Error(errorType, message, function);
+      }
+      // Outside the CCD the rational model is invalid; return input unchanged.
+      if (fabs(dx) > cassisFocalBound || fabs(dy) > cassisFocalBound) {
+        ux = dx;
+        uy = dy;
+        return;
+      }
+      const double *A1_corr = &opticalDistCoeffs[0];
+      const double *A2_corr = &opticalDistCoeffs[6];
+      const double *A3_corr = &opticalDistCoeffs[12];
+      double divider = cassisChiDotA(dx, dy, A3_corr);
+      ux = cassisChiDotA(dx, dy, A1_corr) / divider;
+      uy = cassisChiDotA(dx, dy, A2_corr) / divider;
+      return;
+    } break;
+
   }
 }
 
@@ -705,6 +748,33 @@ void applyDistortion(double ux, double uy, double &dx, double &dy,
       return;
     } break;
 
+    // TGO CaSSIS rational model, undistorted -> distorted (DIST direction).
+    // Ported from ISIS TgoCassisDistortionMap::SetUndistortedFocalPlane. Uses the
+    // DIST triple (indices 18..35) of the 36 coefficients. Closed form, no solver.
+    case CASSIS: {
+      if (opticalDistCoeffs.size() != 36) {
+        csm::Error::ErrorType errorType = csm::Error::INDEX_OUT_OF_RANGE;
+        std::string message =
+            "Distortion coefficients for CaSSIS must be of size 36, "
+            "current size: " +
+            std::to_string(opticalDistCoeffs.size());
+        std::string function = "applyDistortion";
+        throw csm::Error(errorType, message, function);
+      }
+      if (fabs(ux) > cassisFocalBound || fabs(uy) > cassisFocalBound) {
+        dx = ux;
+        dy = uy;
+        return;
+      }
+      const double *A1_dist = &opticalDistCoeffs[18];
+      const double *A2_dist = &opticalDistCoeffs[24];
+      const double *A3_dist = &opticalDistCoeffs[30];
+      double divider = cassisChiDotA(ux, uy, A3_dist);
+      dx = cassisChiDotA(ux, uy, A1_dist) / divider;
+      dy = cassisChiDotA(ux, uy, A2_dist) / divider;
+      return;
+    } break;
+
     // Compute undistorted focal plane coordinate given a distorted
     // focal plane coordinate. This case works by iteratively adding distortion
     // until the new distorted point, r, undistorts to within a tolerance of the

src/Utilities.cpp

@@ -1834,6 +1834,8 @@ DistortionType getDistortionModel(json isd, csm::WarningList *list) {
       return DistortionType::RADTAN;
     } else if (distortion.compare("kplo_shadowcam") == 0) {
       return DistortionType::KPLOSHADOWCAM;
+    } else if (distortion.compare("cassis") == 0) {
+      return DistortionType::CASSIS;
     }
   } catch (...) {
     if (list) {
@@ -1886,6 +1888,8 @@ DistortionType getDistortionModel(int aleDistortionModel,
       return DistortionType::RADTAN;
     } else if (aleDistortionType == ale::DistortionType::KPLOSHADOWCAM) {
       return DistortionType::KPLOSHADOWCAM;
+    } else if (aleDistortionType == ale::DistortionType::CASSIS) {
+      return DistortionType::CASSIS;
     }
   } catch (...) {
     if (list) {
@@ -2099,8 +2103,26 @@ std::vector<double> getDistortionCoeffs(json isd, csm::WarningList *list) {
         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 (...) {
+        if (list) {
+          list->push_back(csm::Warning(
+              csm::Warning::DATA_NOT_AVAILABLE,
+              "Could not parse the cassis distortion model coefficients.",
+              "Utilities::getDistortion()"));
+        }
+        coefficients = std::vector<double>(36, 0.0);
+      }
+    } break;
   }
-  
+
   if (list) {
     list->push_back(
         csm::Warning(csm::Warning::DATA_NOT_AVAILABLE,

tests/DistortionTests.cpp

@@ -486,3 +486,86 @@ TEST(KploShadowCamMapping, intToEnum) {
       static_cast<int>(ale::DistortionType::KPLOSHADOWCAM));
   EXPECT_EQ(dt, DistortionType::KPLOSHADOWCAM);
 }
+
+// TGO CaSSIS rational ratio-of-quadratics distortion (Tulyakov/Ivanov, EPFL),
+// matching ISIS TgoCassisDistortionMap. With chi = [x^2, x*y, y^2, x, y, 1],
+// each direction is x_out = (A1.chi)/(A3.chi), y_out = (A2.chi)/(A3.chi). The 36
+// coefficients are A1_corr, A2_corr, A3_corr (distorted -> undistorted, used by
+// removeDistortion) then A1_dist, A2_dist, A3_dist (undistorted -> distorted,
+// used by applyDistortion), 6 each. Values are the real coefficients from the
+// ISIS addendum tgoCassisAddendum007.ti. Expected outputs below were computed
+// independently from the rational formula.
+static const std::vector<double> cassisCoeffs = {
+    0.00376130530948266,  -0.0134154156065812,   -1.86749521007237e-05,
+    1.00021352681836,     -0.000432362371703953, -0.000948065735350123,
+    9.9842559363676e-05,   0.00373543707958162,  -0.0133299918873929,
+   -0.000215311328389359,  0.995296015537294,    -0.0183542717710778,
+   -3.13320167004204e-05, -7.35655125749807e-06, -1.57664245066771e-05,
+    0.00373549465439151,  -0.0141671946930935,    1.0,
+    0.00213658795560622,  -0.00711785765064197,   1.10355974742147e-05,
+    0.573607182625377,     0.000250884350194894,  0.000550623913037132,
+   -5.69725741015406e-05,  0.00215155905679149,  -0.00716392991767185,
+    0.000124152787728634,  0.576459544392426,     0.010576940564854,
+    1.78250771483506e-05,  4.24592743471094e-06,  9.51220699036653e-06,
+    0.00215158425420738,  -0.0066835595774833,    0.573741540971609};
+
+TEST(Cassis, removeDistortion) {
+  double ux, uy;
+  removeDistortion(2.0, -8.0, ux, uy, cassisCoeffs, 874.9, DistortionType::CASSIS);
+  EXPECT_NEAR(ux, 1.9927225905155812, 1e-9);
+  EXPECT_NEAR(uy, -7.942225998826645, 1e-9);
+}
+
+TEST(Cassis, applyDistortion) {
+  double dx, dy;
+  applyDistortion(2.0, -8.0, dx, dy, cassisCoeffs, 874.9, DistortionType::CASSIS);
+  EXPECT_NEAR(dx, 2.007349177383467, 1e-9);
+  EXPECT_NEAR(dy, -8.058521300253895, 1e-9);
+}
+
+TEST(Cassis, roundTripApplyRemove) {
+  // CORR and DIST are independent EPFL fits, not exact inverses, so the round
+  // trip recovers the input to the fit residual (~5e-6 mm), as in ISIS.
+  for (auto uv : std::vector<std::pair<double, double>>{
+           {1.5, -7.0}, {-3.0, -9.0}, {0.5, -8.5}}) {
+    double dx, dy, ux, uy;
+    applyDistortion(uv.first, uv.second, dx, dy, cassisCoeffs, 874.9,
+                    DistortionType::CASSIS);
+    removeDistortion(dx, dy, ux, uy, cassisCoeffs, 874.9, DistortionType::CASSIS);
+    EXPECT_NEAR(ux, uv.first, 1e-4);
+    EXPECT_NEAR(uy, uv.second, 1e-4);
+  }
+}
+
+TEST(Cassis, offCcdIdentity) {
+  // Outside the CCD (|x| or |y| > ~10.44 mm) the model returns the input
+  // unchanged, matching ISIS.
+  double ux, uy, dx, dy;
+  removeDistortion(15.0, 0.0, ux, uy, cassisCoeffs, 874.9, DistortionType::CASSIS);
+  EXPECT_EQ(ux, 15.0);
+  EXPECT_EQ(uy, 0.0);
+  applyDistortion(0.0, -12.0, dx, dy, cassisCoeffs, 874.9, DistortionType::CASSIS);
+  EXPECT_EQ(dx, 0.0);
+  EXPECT_EQ(dy, -12.0);
+}
+
+TEST(Cassis, wrongCoefficientCount) {
+  double ux, uy;
+  std::vector<double> badCoeffs(10, 0.0);
+  EXPECT_ANY_THROW(removeDistortion(1.0, 1.0, ux, uy, badCoeffs, 874.9,
+                                    DistortionType::CASSIS));
+}
+
+TEST(CassisMapping, stringToEnum) {
+  nlohmann::json isd;
+  isd["optical_distortion"]["cassis"]["coefficients"] = cassisCoeffs;
+  EXPECT_EQ(getDistortionModel(isd), DistortionType::CASSIS);
+}
+
+TEST(CassisMapping, intToEnum) {
+  // ALE and USGSCSM keep independent DistortionType enums; this pins the
+  // cross-library integer mapping for CASSIS.
+  DistortionType dt =
+      getDistortionModel(static_cast<int>(ale::DistortionType::CASSIS));
+  EXPECT_EQ(dt, DistortionType::CASSIS);
+}