diff --git a/ale/base/type_distortion.py b/ale/base/type_distortion.py index f8421711a6a21923103ba4246817d0b44889f972..0649b1c1da0624359977fa113dff50311ea0fc03 100644 --- a/ale/base/type_distortion.py +++ b/ale/base/type_distortion.py @@ -1,3 +1,5 @@ +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] + } + } diff --git a/ale/base/type_sensor.py b/ale/base/type_sensor.py index f49e6c7ccffc5ac68cf1ee64ecf8350408d918e0..6c70eec580b1bcb452a609440c52b98300d78946 100644 --- a/ale/base/type_sensor.py +++ b/ale/base/type_sensor.py @@ -61,7 +61,7 @@ class LineScanner(): ephemeris times split based on image lines """ if not hasattr(self, "_ephemeris_time"): - self._ephemeris_time = np.linspace(self.ephemeris_start_time, self.ephemeris_stop_time, self.image_lines + 1) + self._ephemeris_time = np.linspace(self.ephemeris_start_time, self.ephemeris_stop_time, self.image_lines + 1) return self._ephemeris_time @property diff --git a/ale/drivers/msl_drivers.py b/ale/drivers/msl_drivers.py index e826e560ff3cf5628da351f7c6d22ab96c9ad342..7af3d1ef3546a90c8c28ab4b25e2bd0a726784e0 100644 --- a/ale/drivers/msl_drivers.py +++ b/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 @@ -106,7 +109,23 @@ class MslMastcamPds3NaifSpiceDriver(Cahvor, Framer, Pds3Label, NaifSpice, NoDist focal plane to detector samples """ 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): """ diff --git a/include/ale/Distortion.h b/include/ale/Distortion.h index dd51f8646313d01595707a0a991860aed5f98bb3..07ad43d4889a3ffc651af0557c2210c5b131b755 100644 --- a/include/ale/Distortion.h +++ b/include/ale/Distortion.h @@ -7,7 +7,8 @@ namespace ale { RADIAL, KAGUYALISM, DAWNFC, - LROLROCNAC + LROLROCNAC, + CAHVOR }; } diff --git a/src/Util.cpp b/src/Util.cpp index d62a5ae5311977c4e0d9ef136c6884ed9fbf98a8..297c46c1cb8ffd580cfe1147e4c480098ff30f92 100644 --- a/src/Util.cpp +++ b/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."); diff --git a/tests/pytests/test_cahvor_mixin.py b/tests/pytests/test_cahvor_mixin.py index 0a2600975a2e025deddcc4d37805978b85aeb1ed..26eee41f9aad1ec1838e47686f1e9f7133cb077e 100644 --- a/tests/pytests/test_cahvor_mixin.py +++ b/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()) diff --git a/tests/pytests/test_distortion.py b/tests/pytests/test_distortion.py index 3e8faeca09080f155e42015f942552f820104082..d560c7444fed5317b8aacaac9958d9f8bb6e72e1 100644 --- a/tests/pytests/test_distortion.py +++ b/tests/pytests/test_distortion.py @@ -1,9 +1,9 @@ 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]