Cahvor Distortion (#503)

* Adds CahvorDistortion mixin for MSL

* Small CahvorDistortion fix and adds test for CahvorDistortion

* Fixed CahvorDistortion to return all zero coefficients and offsets if no OR(E) in model

* Changed cahvor distortion test

* Reduced number of lines in line scan cameras if we don't care about the exact ck times

* More cahvor stuff

* Addressed PR feedback

* Removed extra line reduction code

* Updated math to compute x, y distortion offsets

* Changed assert to assert_allclose

ale/base/type_distortion.py

+import numpy as np
+
 class LegendreDistortion():
     """
     Mix-in for sensors that use a legendre distortion model.
@@ -106,3 +108,46 @@ class KaguyaSeleneDistortion():
                 "boresight_y" : self.boresight_y
             }
         }
+
+class CahvorDistortion():
+    """
+    Mix-in for sensors and data sets that have a CAHVOR distortion model.
+
+    This model is based on info from https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2003JE002199
+    Expects that cahvor_camera_dict and focal_length to be defined. This should be
+    a dictionary defining the CAHVOR distortion parameters.
+    """
+
+    @property
+    def usgscsm_distortion_model(self):
+        """
+        Dictionary containing the usgscsm specification for CAHVOR distortion.
+        This will be a list of coeffs for radial distortion (x0, x1, x2)
+        followed by the optical center (x, y)
+
+        Returns
+        -------
+        : dict
+        """
+        R = [0, 0, 0]
+        x = 0
+        y = 0
+        # If our model contains OR in the CAHVOR model
+        # then compute the distortion coeffs/offset
+        if (len(self.cahvor_camera_dict.keys()) >= 6):
+            A = self.cahvor_camera_dict.get("A", [0, 0, 0])
+            H = self.cahvor_camera_dict.get("H", [0, 0, 0])
+            V = self.cahvor_camera_dict.get("V", [0, 0, 0])
+            O = self.cahvor_camera_dict.get("O", [0, 0, 0])
+            i = np.dot(O, H) / np.dot(O, A)
+            j = np.dot(O, V) / np.dot(O, A)
+            x = self.pixel_size * (i - self.compute_h_c())
+            y = self.pixel_size * (self.compute_v_c() - j)
+            R = self.cahvor_camera_dict.get("R", [0, 0, 0])
+            R[1] /= self.focal_length**2
+            R[2] /= self.focal_length**4
+        return {
+            "cahvor": {
+                "coefficients": [*R, x, y]
+            }
+        }

ale/drivers/msl_drivers.py

@@ -4,11 +4,13 @@ import spiceypy as spice
 from ale.base.data_naif import NaifSpice
 from ale.base.label_pds3 import Pds3Label
 from ale.base.type_sensor import Framer
-from ale.base.type_distortion import NoDistortion
+from ale.base.type_distortion import CahvorDistortion
 from ale.base.type_sensor import Cahvor
 from ale.base.base import Driver
 
-class MslMastcamPds3NaifSpiceDriver(Cahvor, Framer, Pds3Label, NaifSpice, NoDistortion, Driver):
+class MslMastcamPds3NaifSpiceDriver(Cahvor, Framer, Pds3Label, NaifSpice, CahvorDistortion, Driver):
+    """
+    """
     @property
     def spacecraft_name(self):
         """
@@ -73,6 +75,7 @@ class MslMastcamPds3NaifSpiceDriver(Cahvor, Framer, Pds3Label, NaifSpice, NoDist
         Returns the Naif ID code for the site reference frame
         Expects REFERENCE_COORD_SYSTEM_INDEX to be defined in the camera
         PVL group. 
+
         Returns
         -------
         : int
@@ -107,6 +110,22 @@ class MslMastcamPds3NaifSpiceDriver(Cahvor, Framer, Pds3Label, NaifSpice, NoDist
         """
         return [0, 0, 1/self.pixel_size]
 
+    @property
+    def detector_center_line(self):
+        return self.label["INSTRUMENT_STATE_PARMS"]["DETECTOR_LINES"]/2
+
+    @property
+    def detector_center_sample(self):
+        return self.label["INSTRUMENT_STATE_PARMS"]["MSL:DETECTOR_SAMPLES"]/2
+
+    @property
+    def detector_start_line(self):
+        return self.label["IMAGE_REQUEST_PARMS"]["FIRST_LINE"]
+
+    @property
+    def detector_start_sample(self):
+        return self.label["IMAGE_REQUEST_PARMS"]["FIRST_LINE_SAMPLE"]
+
     @property
     def sensor_model_version(self):
         """

include/ale/Distortion.h

@@ -7,7 +7,8 @@ namespace ale {
     RADIAL,
     KAGUYALISM,
     DAWNFC,
-    LROLROCNAC
+    LROLROCNAC,
+    CAHVOR
   };
 }
 

src/Util.cpp

@@ -333,6 +333,8 @@ DistortionType getDistortionModel(json isd) {
       return DistortionType::DAWNFC;
     } else if (distortion.compare("lrolrocnac") == 0) {
       return DistortionType::LROLROCNAC;
+    } else if (distortion.compare("cahvor") == 0) {
+      return DistortionType::CAHVOR;
     }
   } catch (...) {
     throw std::runtime_error("Could not parse the distortion model.");
@@ -450,6 +452,25 @@ std::vector<double> getDistortionCoeffs(json isd) {
       coefficients = std::vector<double>(1, 0.0);
     }
   } break;
+  case DistortionType::CAHVOR:
+  {
+    try
+    {
+      coefficients = isd.at("optical_distortion")
+                         .at("cahvor")
+                         .at("coefficients")
+                         .get<std::vector<double>>();
+
+      return coefficients;
+    }
+    catch (...)
+    {
+      throw std::runtime_error(
+          "Could not parse the cahvor distortion model coefficients.");
+      coefficients = std::vector<double>(5, 0.0);
+    }
+  }
+  break;
   }
   throw std::runtime_error(
       "Could not parse the distortion model coefficients.");

tests/pytests/test_cahvor_mixin.py

@@ -38,7 +38,6 @@ class test_cahvor_sensor(unittest.TestCase):
     @patch("ale.base.type_sensor.Cahvor.cahvor_camera_dict", new_callable=PropertyMock, return_value=cahvor_camera_dict())
     def test_cahvor_model_elements(self, cahvor_camera_dict):
         cahvor_matrix = self.driver.cahvor_rotation_matrix
-        print(cahvor_matrix)
         np.testing.assert_allclose(cahvor_matrix, [[-0.82067034, -0.57129702, 0.01095033],
                                                    [-0.43920248, 0.6184257, -0.65165238],
                                                    [ 0.3655151, -0.5396012, -0.7584387 ]])
@@ -51,7 +50,6 @@ class test_cahvor_sensor(unittest.TestCase):
       assert -76000 in frame_chain.nodes()
       assert -76562 in frame_chain.nodes()
       from_spice.assert_called_with(center_ephemeris_time=0, ephemeris_times=[0], sensor_frame=-76000, target_frame=10014, nadir=False, exact_ck_times=False)
-      print(frame_chain[-76562][-76000]['rotation'].quat)
       np.testing.assert_allclose(frame_chain[-76562][-76000]['rotation'].quat, [-0.28255205, 0.8940826, -0.33309383, 0.09914206])
 
     @patch("ale.base.type_sensor.Cahvor.cahvor_camera_dict", new_callable=PropertyMock, return_value=cahvor_camera_dict())

tests/pytests/test_distortion.py

 import pytest
-import pvl
+from unittest.mock import MagicMock
 
-import ale
-from ale import base
-from ale.base.type_distortion import RadialDistortion, KaguyaSeleneDistortion
+import numpy as np
+
+from ale.base.type_distortion import CahvorDistortion, RadialDistortion, KaguyaSeleneDistortion
 
 def test_radial_distortion():
     radial_distortion = RadialDistortion()
@@ -21,3 +21,33 @@ def test_kaguyaselene_distortion():
     assert kaguyaselene_distortion.usgscsm_distortion_model['kaguyalism']['y'] == [1,2,3,4]
     assert kaguyaselene_distortion.usgscsm_distortion_model['kaguyalism']['boresight_x'] == 0.1
     assert kaguyaselene_distortion.usgscsm_distortion_model['kaguyalism']['boresight_y'] == -0.1
+
+@pytest.fixture()
+def cahvor_camera_dict():
+    camera_dict = {}
+    camera_dict['C'] = np.array([0.9050933, 0.475724, -1.972196])
+    camera_dict['A'] = np.array([0.3791357, 0.9249998, 0.02481062])
+    camera_dict['H'] = np.array([-4099.206, 2247.006, 34.41227])
+    camera_dict['V'] = np.array([59.37479, 85.92623, 4648.751])
+    camera_dict['O'] = np.array([0.377433, 0.9257466, 0.02284064])
+    camera_dict['R'] = np.array([-1.510000e-04, -1.391890e-01, -1.250336e+00])
+    return camera_dict
+
+def test_cahvor_distortion(cahvor_camera_dict):
+    cahvor_distortion = CahvorDistortion()
+    cahvor_distortion.pixel_size = 0.007319440022615588
+    cahvor_distortion.focal_length = 34.0
+
+    cahvor_distortion.cahvor_camera_dict = cahvor_camera_dict
+    h_c = np.dot(cahvor_distortion.cahvor_camera_dict['A'], cahvor_distortion.cahvor_camera_dict['H'])
+    v_c = np.dot(cahvor_distortion.cahvor_camera_dict['A'], cahvor_distortion.cahvor_camera_dict['V'])
+    cahvor_distortion.compute_h_c = MagicMock(return_value=h_c)
+    cahvor_distortion.compute_v_c = MagicMock(return_value=v_c)
+    coefficients = cahvor_distortion.usgscsm_distortion_model['cahvor']['coefficients']
+    np.testing.assert_allclose(coefficients, [-0.000151, -0.00012040570934256056, -9.35644927622993e-07, 0.06295036132043122, 0.06727152372705038])
+
+    cahvor_distortion.cahvor_camera_dict = cahvor_camera_dict
+    cahvor_distortion.cahvor_camera_dict.pop('O')
+    cahvor_distortion.cahvor_camera_dict.pop('R')
+    coefficients = cahvor_distortion.usgscsm_distortion_model['cahvor']['coefficients']
+    assert coefficients == [0, 0, 0, 0, 0]