From f098d9a721c2d0708a60a46016717177b67aaee7 Mon Sep 17 00:00:00 2001
From: Jesse Mapel <jmapel@usgs.gov>
Date: Wed, 22 Jun 2022 10:24:53 -0700
Subject: [PATCH] Initial frame jitter (#399)

---
 include/usgscsm/UsgsAstroFrameSensorModel.h |   5 +-
 include/usgscsm/Utilities.h                 |  13 ++
 src/UsgsAstroFrameSensorModel.cpp           |  57 ++++++--
 src/Utilities.cpp                           |  82 ++++++++++++
 tests/Fixtures.h                            |  32 ++++-
 tests/FrameCameraTests.cpp                  |  29 ++++
 tests/UtilitiesTests.cpp                    | 139 ++++++++++++++++++++
 tests/data/jitterFramer.json                |  97 ++++++++++++++
 8 files changed, 441 insertions(+), 13 deletions(-)
 create mode 100644 tests/data/jitterFramer.json

diff --git a/include/usgscsm/UsgsAstroFrameSensorModel.h b/include/usgscsm/UsgsAstroFrameSensorModel.h
index b24d4ea..9ee36fa 100644
--- a/include/usgscsm/UsgsAstroFrameSensorModel.h
+++ b/include/usgscsm/UsgsAstroFrameSensorModel.h
@@ -38,7 +38,7 @@ class UsgsAstroFrameSensorModel : public csm::RasterGM,
 
   std::string constructStateFromIsd(const std::string &jsonIsd,
                                     csm::WarningList *warnings);
-  
+
   // Apply a rotation and translation to a state string. The effect is
   // to transform the position and orientation of the camera in ECEF
   // coordinates.
@@ -365,6 +365,9 @@ class UsgsAstroFrameSensorModel : public csm::RasterGM,
   std::vector<double> m_iTransS;
   std::vector<double> m_iTransL;
   std::vector<double> m_boresight;
+  std::vector<double> m_lineJitter;
+  std::vector<double> m_sampleJitter;
+  std::vector<double> m_lineTimes;
   double m_majorAxis;
   double m_minorAxis;
   double m_focalLength;
diff --git a/include/usgscsm/Utilities.h b/include/usgscsm/Utilities.h
index 3636762..1378d99 100644
--- a/include/usgscsm/Utilities.h
+++ b/include/usgscsm/Utilities.h
@@ -24,6 +24,19 @@ void computeDistortedFocalPlaneCoordinates(
     const double &startingLine, const double iTransS[], const double iTransL[],
     double &distortedX, double &distortedY);
 
+void removeJitter(
+    const double &line, const double &sample,
+    const std::vector<double> lineJitterCoeffs, const std::vector<double> sampleJitterCoeffs,
+    const std::vector<double> lineTimes,
+    double &dejitteredLine, double &dejitteredSample);
+
+void addJitter(
+    const double &line, const double &sample,
+    const double &tolerance, const int &maxIts,
+    const std::vector<double> lineJitterCoeffs, const std::vector<double> sampleJitterCoeffs,
+    const std::vector<double> lineTimes,
+    double &jitteredLine, double &jitteredSample);
+
 void computePixel(const double &distortedX, const double &distortedY,
                   const double &sampleOrigin, const double &lineOrigin,
                   const double &sampleSumming, const double &lineSumming,
diff --git a/src/UsgsAstroFrameSensorModel.cpp b/src/UsgsAstroFrameSensorModel.cpp
index 763af33..3c4e36b 100644
--- a/src/UsgsAstroFrameSensorModel.cpp
+++ b/src/UsgsAstroFrameSensorModel.cpp
@@ -81,6 +81,9 @@ void UsgsAstroFrameSensorModel::reset() {
   m_iTransS = std::vector<double>(3, 0.0);
   m_iTransL = std::vector<double>(3, 0.0);
   m_boresight = std::vector<double>(3, 0.0);
+  m_lineJitter.clear();
+  m_sampleJitter.clear();
+  m_lineTimes.clear();
   m_parameterType =
       std::vector<csm::param::Type>(NUM_PARAMETERS, csm::param::REAL);
   m_referencePointXyz.x = 0;
@@ -163,6 +166,18 @@ csm::ImageCoord UsgsAstroFrameSensorModel::groundToImage(
                m_startingDetectorSample, m_startingDetectorLine, &m_iTransS[0],
                &m_iTransL[0], line, sample);
 
+  // Optionally apply rolling shutter jitter correction
+  if (m_lineTimes.size()) {
+    double jitteredLine, jitteredSample;
+    addJitter(
+        line, sample,
+        desired_precision, 20,
+        m_lineJitter, m_sampleJitter, m_lineTimes,
+        jitteredLine, jitteredSample);
+    line = jitteredLine;
+    sample = jitteredSample;
+  }
+
   MESSAGE_LOG("Computed groundToImage for {}, {}, {} as line, sample: {}, {}",
               groundPt.x, groundPt.y, groundPt.z, line, sample);
 
@@ -211,8 +226,13 @@ csm::EcefCoord UsgsAstroFrameSensorModel::imageToGround(
       m[0][0], m[0][1], m[0][2], m[1][0], m[1][1], m[1][2], m[2][0], m[2][1],
       m[2][2]);
 
-  // Apply the principal point offset, assuming the pp is given in pixels
-  double xl, yl, zl;
+  // Optionally apply rolling shutter jitter correction
+  if (m_lineTimes.size()) {
+    double dejitteredLine, dejitteredSample;
+    removeJitter(line, sample, m_lineJitter, m_sampleJitter, m_lineTimes, dejitteredLine, dejitteredSample);
+    line = dejitteredLine;
+    sample = dejitteredSample;
+  }
 
   // Convert from the pixel space into the metric space
   double x_camera, y_camera;
@@ -229,6 +249,7 @@ csm::EcefCoord UsgsAstroFrameSensorModel::imageToGround(
               undistortedY);
 
   // Now back from distorted mm to pixels
+  double xl, yl, zl;
   xl = m[0][0] * undistortedX + m[0][1] * undistortedY -
        m[0][2] * -m_focalLength;
   yl = m[1][0] * undistortedX + m[1][1] * undistortedY -
@@ -737,6 +758,9 @@ std::string UsgsAstroFrameSensorModel::getModelState() const {
       {"m_transY", {m_transY[0], m_transY[1], m_transY[2]}},
       {"m_iTransS", {m_iTransS[0], m_iTransS[1], m_iTransS[2]}},
       {"m_iTransL", {m_iTransL[0], m_iTransL[1], m_iTransL[2]}},
+      {"m_lineJitter", m_lineJitter},
+      {"m_sampleJitter", m_sampleJitter},
+      {"m_lineTimes", m_lineTimes},
       {"m_majorAxis", m_majorAxis},
       {"m_minorAxis", m_minorAxis},
       {"m_spacecraftVelocity",
@@ -783,7 +807,7 @@ void UsgsAstroFrameSensorModel::applyTransformToState(ale::Rotation const& r, al
   // The vector having the rotation and translation from camera to
   // ECEF
   std::vector<double> currentParameterValue = j["m_currentParameterValue"];
-  
+
   // Extract the quaternions from the frame camera, apply the
   // rotation, and put them back.
   std::vector<double> quaternions;
@@ -792,7 +816,7 @@ void UsgsAstroFrameSensorModel::applyTransformToState(ale::Rotation const& r, al
   applyRotationToQuatVec(r, quaternions);
   for (size_t it = 0; it < 4; it++)
     currentParameterValue[it + 3] = quaternions[it];
-  
+
   // Extract the positions from the frame camera, apply to them the
   // rotation and translation, and put them back.
   std::vector<double> positions;
@@ -949,17 +973,23 @@ void UsgsAstroFrameSensorModel::replaceModelState(
         state.at("m_currentParameterCovariance").get<std::vector<double>>();
 
     // These are optional and may not exist in all state files
-    if (state.find("m_maxElevation") != state.end()) 
+    if (state.find("m_maxElevation") != state.end())
       m_maxElevation = state.at("m_maxElevation").get<double>();
-    if (state.find("m_minElevation") != state.end()) 
+    if (state.find("m_minElevation") != state.end())
       m_minElevation = state.at("m_minElevation").get<double>();
-    if (state.find("m_originalHalfLines") != state.end()) 
+    if (state.find("m_originalHalfLines") != state.end())
       m_originalHalfLines = state.at("m_originalHalfLines").get<double>();
-    if (state.find("m_originalHalfSamples") != state.end()) 
+    if (state.find("m_originalHalfSamples") != state.end())
       m_originalHalfSamples = state.at("m_originalHalfSamples").get<double>();
-    if (state.find("m_pixelPitch") != state.end()) 
+    if (state.find("m_pixelPitch") != state.end())
       m_pixelPitch = state.at("m_pixelPitch").get<double>();
-    
+    if (state.find("m_lineJitter") != state.end())
+      m_lineJitter = state.at("m_lineJitter").get<std::vector<double>>();
+    if (state.find("m_sampleJitter") != state.end())
+      m_sampleJitter = state.at("m_sampleJitter").get<std::vector<double>>();
+    if (state.find("m_lineTimes") != state.end())
+      m_lineTimes = state.at("m_lineTimes").get<std::vector<double>>();
+
   } catch (std::out_of_range &e) {
     MESSAGE_LOG("State keywords required to generate sensor model missing: " +
                 std::string(e.what()) + "\nUsing model string: " + stringState +
@@ -1150,6 +1180,13 @@ std::string UsgsAstroFrameSensorModel::constructStateFromIsd(
   state["m_iTransL"] = ale::getFocal2PixelLines(parsedIsd);
   state["m_iTransS"] = ale::getFocal2PixelSamples(parsedIsd);
 
+  // optional rolling shutter jitter
+  if (parsedIsd.find("jitter") != parsedIsd.end()) {
+    state["m_lineJitter"] = parsedIsd.at("jitter").at("lineJitterCoefficients");
+    state["m_sampleJitter"] = parsedIsd.at("jitter").at("sampleJitterCoefficients");
+    state["m_lineTimes"] = parsedIsd.at("jitter").at("lineExposureTimes");
+  }
+
   // We don't pass the pixel to focal plane transformation so invert the
   // focal plane to pixel transformation
   try {
diff --git a/src/Utilities.cpp b/src/Utilities.cpp
index e0c552e..4640229 100644
--- a/src/Utilities.cpp
+++ b/src/Utilities.cpp
@@ -1,5 +1,6 @@
 #include "Utilities.h"
 
+
 #include <Error.h>
 #include <cmath>
 #include <stack>
@@ -88,6 +89,87 @@ void computeDistortedFocalPlaneCoordinates(
   distortedY = p21 * t1 + p22 * t2;
 }
 
+
+// Compute the de-jittered pixel coordinate given a jittered image coordinate
+// a set of jitter coefficients and line exposure times for rolling shutter.
+// Jitter coefficients are in largest power first order. There is no constant
+// coefficient. For example {1, 2, 3} would correspond to 1*t^3 + 2*t&2 + 3*t.
+void removeJitter(
+    const double &line, const double &sample,
+    const std::vector<double> lineJitterCoeffs, const std::vector<double> sampleJitterCoeffs,
+    const std::vector<double> lineTimes,
+    double &dejitteredLine, double &dejitteredSample) {
+  // Check input
+  if (lineJitterCoeffs.size() != sampleJitterCoeffs.size() ||
+      lineTimes.size() == 0) {
+    throw csm::Error(
+        csm::Error::INDEX_OUT_OF_RANGE,
+        "Jitter coefficient vectors must be the same size.",
+        "removeJitter");
+  }
+  if (lineTimes.size() == 0) {
+    throw csm::Error(
+        csm::Error::INDEX_OUT_OF_RANGE,
+        "Line exposure times must be non-empty.",
+        "removeJitter");
+  }
+  double lineJitter = 0;
+  double sampleJitter = 0;
+  // Bound line index to the vector of line exposure times;
+  double time = lineTimes[std::max(std::min((int)std::round(line), (int)lineTimes.size()), 1) - 1];
+  for (unsigned int n = 0; n < lineJitterCoeffs.size(); n++) {
+    double timeTerm = pow(time, lineJitterCoeffs.size() - n);
+    lineJitter += lineJitterCoeffs[n] * timeTerm;
+    sampleJitter += sampleJitterCoeffs[n] * timeTerm;
+  }
+  dejitteredLine = line - lineJitter;
+  dejitteredSample = sample - sampleJitter;
+
+  return;
+}
+
+
+// Compute the jittered pixel coordinate given a de-jittered image coordinate
+// a set of jitter coefficients and line exposure times for rolling shutter.
+// Jitter coefficients are in largest power first order. There is no constant
+// coefficient. For example {1, 2, 3} would correspond to 1*t^3 + 2*t&2 + 3*t.
+// This uses an iterative method so a tolerance and maximum number of iteration
+// are required to determine when to stop iterating.
+void addJitter(
+    const double &line, const double &sample,
+    const double &tolerance, const int &maxIts,
+    const std::vector<double> lineJitterCoeffs, const std::vector<double> sampleJitterCoeffs,
+    const std::vector<double> lineTimes,
+    double &jitteredLine, double &jitteredSample) {
+  int iteration = 0;
+  double dejitteredLine = line - 1;
+  double dejitteredSample = sample - 1;
+  double currentLine = line;
+  double currentSample =  sample;
+
+  while (iteration < maxIts) {
+    removeJitter(
+        currentLine, currentSample,
+        lineJitterCoeffs, sampleJitterCoeffs, lineTimes,
+        dejitteredLine, dejitteredSample);
+
+    if (fabs(dejitteredLine - line) < tolerance &&
+        fabs(dejitteredSample - sample) < tolerance) {
+      break;
+    }
+
+    currentLine = line + currentLine - dejitteredLine;
+    currentSample = sample + currentSample - dejitteredSample;
+    iteration++;
+  }
+
+  jitteredLine = currentLine;
+  jitteredSample = currentSample;
+
+  return;
+}
+
+
 // Compute the image pixel for a distorted focal plane coordinate
 // in - line
 // in - sample
diff --git a/tests/Fixtures.h b/tests/Fixtures.h
index 35a2251..350bf8f 100644
--- a/tests/Fixtures.h
+++ b/tests/Fixtures.h
@@ -131,6 +131,34 @@ class OrbitalFrameSensorModel : public ::testing::Test {
   }
 };
 
+class JitterFrameSensorModel : public ::testing::Test {
+ protected:
+  csm::Isd jitterIsd;
+  csm::Isd isd;
+  std::shared_ptr<csm::RasterGM> jitterModel;
+  std::shared_ptr<csm::RasterGM> model;
+
+  void SetUp() override {
+    csm::Model *genericModel = nullptr;
+
+    UsgsAstroPlugin cameraPlugin;
+
+    isd.setFilename("data/orbitalFramer.img");
+
+    genericModel = cameraPlugin.constructModelFromISD(
+      isd, UsgsAstroFrameSensorModel::_SENSOR_MODEL_NAME);
+    model = std::shared_ptr<csm::RasterGM>(dynamic_cast<csm::RasterGM *>(genericModel));
+    ASSERT_NE(model, nullptr);
+
+    jitterIsd.setFilename("data/jitterFramer.img");
+
+    genericModel = cameraPlugin.constructModelFromISD(
+      jitterIsd, UsgsAstroFrameSensorModel::_SENSOR_MODEL_NAME);
+    jitterModel = std::shared_ptr<csm::RasterGM>(dynamic_cast<csm::RasterGM *>(genericModel));
+    ASSERT_NE(jitterModel, nullptr);
+  }
+};
+
 class FrameIsdTest : public ::testing::Test {
  protected:
   csm::Isd isd;
@@ -255,7 +283,7 @@ class OrbitalLineScanSensorModel : public ::testing::Test {
     model = std::shared_ptr<csm::Model>(cameraPlugin.constructModelFromISD(
       isd, UsgsAstroLsSensorModel::_SENSOR_MODEL_NAME));
     sensorModel = dynamic_cast<UsgsAstroLsSensorModel *>(model.get());
-    
+
     ASSERT_NE(sensorModel, nullptr);
   }
 
@@ -369,7 +397,7 @@ class OrbitalPushFrameSensorModel : public ::testing::Test {
   void SetUp() override {
     isd.setFilename("data/orbitalPushFrame.img");
     UsgsAstroPlugin cameraPlugin;
-    
+
     csm::Model * model = cameraPlugin.constructModelFromISD(
        isd, UsgsAstroPushFrameSensorModel::_SENSOR_MODEL_NAME);
 
diff --git a/tests/FrameCameraTests.cpp b/tests/FrameCameraTests.cpp
index 535ff30..48208c1 100644
--- a/tests/FrameCameraTests.cpp
+++ b/tests/FrameCameraTests.cpp
@@ -189,6 +189,35 @@ TEST_F(OrbitalFrameSensorModel, ImageToProximateImagingLocus) {
   EXPECT_TRUE(warnings.empty());
 }
 
+// ****************************************************************************
+//   Jitter sensor model tests
+// ****************************************************************************
+
+TEST_F(JitterFrameSensorModel, Center) {
+  // Test jitter is y = -1/256 t^2 + 1/16 t which is 0.25 at t=8
+  csm::ImageCoord imagePt(7.75, 8.25);
+  csm::EcefCoord groundPt = model->imageToGround(imagePt, 0.0);
+
+  csm::ImageCoord jitterImagePt(8.0, 8.0);
+  csm::EcefCoord jitterGroundPt = jitterModel->imageToGround(jitterImagePt, 0.0);
+
+  EXPECT_DOUBLE_EQ(groundPt.x, jitterGroundPt.x);
+  EXPECT_DOUBLE_EQ(groundPt.y, jitterGroundPt.y);
+  EXPECT_DOUBLE_EQ(groundPt.z, jitterGroundPt.z);
+}
+
+TEST_F(JitterFrameSensorModel, Inversion) {
+  for (int line = 0; line <= 16; line++) {
+    for (int sample = 0; sample <= 16; sample++) {
+      csm::ImageCoord imagePt1(line, sample);
+      csm::EcefCoord groundPt = jitterModel->imageToGround(imagePt1, 0.001);
+      csm::ImageCoord imagePt2 = jitterModel->groundToImage(groundPt);
+      EXPECT_NEAR(imagePt1.line, imagePt2.line, 0.001);
+      EXPECT_NEAR(imagePt1.samp, imagePt2.samp, 0.001);
+    }
+  }
+}
+
 // ****************************************************************************
 //   Modified sensor model tests
 // ****************************************************************************
diff --git a/tests/UtilitiesTests.cpp b/tests/UtilitiesTests.cpp
index 86be3f4..b5c1eed 100644
--- a/tests/UtilitiesTests.cpp
+++ b/tests/UtilitiesTests.cpp
@@ -96,6 +96,145 @@ TEST(UtilitiesTests, computeDistortedFocalPlaneCoordinatesStart) {
   EXPECT_DOUBLE_EQ(undistortedFocalPlaneY, -0.2);
 }
 
+TEST(UtilitiesTests, removeJitter) {
+  std::vector<double> lineJit = {2.0, 3.0, 4.0};
+  std::vector<double> sampleJit = {-10.0, -20.0, -30.0};
+  std::vector<double> lineTimes = {0.0, 1.0, 2.0, 3.0, 4.0};
+
+  double testline, testSample;
+  double dejitteredLine, dejitteredSample;
+  double lineJitter, sampleJitter;
+
+  // Test at t=0, this should do nothing
+  testline = 1;
+  testSample = 0;
+  dejitteredLine = 10;
+  dejitteredSample = 10;
+  removeJitter(testline, testSample, lineJit, sampleJit, lineTimes, dejitteredLine, dejitteredSample);
+  EXPECT_DOUBLE_EQ(dejitteredLine, testline);
+  EXPECT_DOUBLE_EQ(dejitteredSample, testSample);
+
+  // Test at t=1 to check proper coefficients
+  testline = 2;
+  testSample = 5.0;
+  dejitteredLine = 10;
+  dejitteredSample = 10;
+  removeJitter(testline, testSample, lineJit, sampleJit, lineTimes, dejitteredLine, dejitteredSample);
+  EXPECT_DOUBLE_EQ(dejitteredLine, testline - (2.0 + 3.0 + 4.0));
+  EXPECT_DOUBLE_EQ(dejitteredSample, testSample + (10.0 + 20.0 + 30.0));
+
+  // Test at t=2 to check exponents
+  testline = 3;
+  testSample = 5.0;
+  dejitteredLine = 10;
+  dejitteredSample = 10;
+  removeJitter(testline, testSample, lineJit, sampleJit, lineTimes, dejitteredLine, dejitteredSample);
+  EXPECT_DOUBLE_EQ(dejitteredLine, testline - (8 * 2.0 + 4 * 3.0 + 2 * 4.0));
+  EXPECT_DOUBLE_EQ(dejitteredSample, testSample + (8 * 10.0 + 4 * 20.0 + 2 * 30.0));
+
+  // Test at extreme negative line to ensure it bounds to the image
+  testline = 1;
+  testSample = 0;
+  dejitteredLine = 10;
+  dejitteredSample = 10;
+  removeJitter(testline, testSample, lineJit, sampleJit, lineTimes, dejitteredLine, dejitteredSample);
+  lineJitter = testline - dejitteredLine;
+  sampleJitter = testSample - dejitteredSample;
+  testline = -100;
+  dejitteredLine = 10;
+  dejitteredSample = 10;
+  removeJitter(testline, testSample, lineJit, sampleJit, lineTimes, dejitteredLine, dejitteredSample);
+  EXPECT_DOUBLE_EQ(testline - dejitteredLine, lineJitter);
+  EXPECT_DOUBLE_EQ(testSample - dejitteredSample, sampleJitter);
+
+  // Test at too large of a line to ensure it bounds to the image
+  testline = 5.0;
+  testSample = 0;
+  dejitteredLine = 10;
+  dejitteredSample = 10;
+  removeJitter(testline, testSample, lineJit, sampleJit, lineTimes, dejitteredLine, dejitteredSample);
+  lineJitter = testline - dejitteredLine;
+  sampleJitter = testSample - dejitteredSample;
+  testline = 100;
+  dejitteredLine = 10;
+  dejitteredSample = 10;
+  removeJitter(testline, testSample, lineJit, sampleJit, lineTimes, dejitteredLine, dejitteredSample);
+  EXPECT_DOUBLE_EQ(testline - dejitteredLine, lineJitter);
+  EXPECT_DOUBLE_EQ(testSample - dejitteredSample, sampleJitter);
+
+  // Test at fractional lines to check rounding
+  testline = 1.0;
+  testSample = 5.0;
+  dejitteredLine = 10;
+  dejitteredSample = 10;
+  removeJitter(testline, testSample, lineJit, sampleJit, lineTimes, dejitteredLine, dejitteredSample);
+  lineJitter = testline - dejitteredLine;
+  sampleJitter = testSample - dejitteredSample;
+  testline = 0.75;
+  dejitteredLine = 10;
+  dejitteredSample = 10;
+  removeJitter(testline, testSample, lineJit, sampleJit, lineTimes, dejitteredLine, dejitteredSample);
+  EXPECT_DOUBLE_EQ(testline - dejitteredLine, lineJitter);
+  EXPECT_DOUBLE_EQ(testSample - dejitteredSample, sampleJitter);
+  testline = 1.35;
+  dejitteredLine = 10;
+  dejitteredSample = 10;
+  removeJitter(testline, testSample, lineJit, sampleJit, lineTimes, dejitteredLine, dejitteredSample);
+  EXPECT_DOUBLE_EQ(testline - dejitteredLine, lineJitter);
+  EXPECT_DOUBLE_EQ(testSample - dejitteredSample, sampleJitter);
+}
+
+TEST(UtilitiesTests, removeJitterErrors) {
+  std::vector<double> lineJit = {2.0, 3.0, 4.0};
+  std::vector<double> shortLineJit = {2.0, 3.0};
+  std::vector<double> sampleJit = {-10.0, -20.0, -30.0};
+  std::vector<double> longSampleJit = {-10.0, -20.0, -30.0, -40.0};
+  std::vector<double> lineTimes = {0.0, 1.0, 2.0, 3.0, 4.0};
+  std::vector<double> emptyLineTimes;
+  double testLine = 4;
+  double testSample = 15;
+  double dejitteredLine, dejitteredSample;
+
+  EXPECT_THROW(
+      removeJitter(testLine, testSample, shortLineJit, sampleJit, lineTimes, dejitteredLine, dejitteredSample),
+      csm::Error);
+  EXPECT_THROW(
+      removeJitter(testLine, testSample, lineJit, longSampleJit, lineTimes, dejitteredLine, dejitteredSample),
+      csm::Error);
+  EXPECT_THROW(
+      removeJitter(testLine, testSample, lineJit, sampleJit, emptyLineTimes, dejitteredLine, dejitteredSample),
+      csm::Error);
+}
+
+TEST(UtilitiesTests, addJitter) {
+  // The line jitter must be sufficiently small to make a unique solution.
+  // Extremely high jitter can result in multiple solutions to the inversion.
+  std::vector<double> lineJit = {0.01, 0.02, 0.03};
+  std::vector<double> sampleJit = {-0.01, -0.02, -0.03};
+  std::vector<double> lineTimes = {0.0, 1.0, 2.0, 3.0, 4.0};
+
+  double testline = 4;
+  double testSample = 15.0;
+  double dejitteredLine = 10;
+  double dejitteredSample = 10;
+  removeJitter(
+      testline, testSample,
+      lineJit, sampleJit, lineTimes,
+      dejitteredLine, dejitteredSample);
+  double jitteredLine = 10;
+  double jitteredSample = 10;
+  double tolerance = 1e-7;
+  int maxIts = 30;
+  addJitter(
+    dejitteredLine, dejitteredSample,
+    tolerance, maxIts,
+    lineJit, sampleJit, lineTimes,
+    jitteredLine, jitteredSample);
+
+  EXPECT_LE(fabs(jitteredLine - testline), tolerance);
+  EXPECT_LE(fabs(jitteredSample - testSample), tolerance);
+}
+
 TEST(UtilitiesTests, computePixel) {
   double iTransS[] = {0.0, 0.0, 10.0};
   double iTransL[] = {0.0, 10.0, 0.0};
diff --git a/tests/data/jitterFramer.json b/tests/data/jitterFramer.json
new file mode 100644
index 0000000..133fd95
--- /dev/null
+++ b/tests/data/jitterFramer.json
@@ -0,0 +1,97 @@
+{
+  "name_model" : "USGS_ASTRO_FRAME_SENSOR_MODEL",
+  "name_platform" : "TEST_PLATFORM",
+  "name_sensor" : "TEST_SENSOR",
+  "center_ephemeris_time": 1000.0,
+  "detector_sample_summing": 1,
+  "detector_line_summing": 1,
+  "focal2pixel_lines": [0.0, 10.0, 0.0],
+  "focal2pixel_samples": [0.0, 0.0, 10.0],
+  "focal_length_model": {
+    "focal_length": 50
+  },
+  "image_lines": 16,
+  "image_samples": 16,
+  "detector_center" : {
+    "line" : 8.0,
+    "sample" : 8.0
+  },
+  "interpolation_method": "lagrange",
+  "optical_distortion": {
+    "radial": {
+      "coefficients": [0.0, 0.0, 0.0]
+    }
+  },
+  "radii": {
+    "semimajor": 1000,
+    "semiminor": 1000,
+    "unit":"km"
+  },
+  "jitter": {
+    "lineJitterCoefficients": [-0.00390625, 0.0625],
+    "sampleJitterCoefficients": [0.00390625, -0.0625],
+    "lineExposureTimes": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
+  },
+  "reference_height": {
+    "maxheight": 1000,
+    "minheight": -1000,
+    "unit": "m"
+  },
+  "instrument_position": {
+    "unit": "km",
+    "positions": [
+      [1050.0, 0.0, 0.0]
+    ],
+    "velocities": [
+      [ 0.0, 0.0, -1.0]
+    ],
+    "ephemeris_times": [
+      50
+    ],
+    "reference_frame": 1
+  },
+  "instrument_pointing": {
+    "quaternions": [
+      [-0.707106781186547, 0.0, -0.707106781186547, 0]
+    ],
+    "angular_velocities": [
+      [0, 0, 0]
+    ],
+    "ephemeris_times": [
+      50
+    ],
+    "reference_frame": 1
+  },
+  "body_rotation": {
+    "quaternions": [
+      [1, 0, 0, 0]
+    ],
+    "ephemeris_times": [
+      50
+    ],
+    "angular_velocities": [
+      [0.0, 0.0, 0.0]
+    ],
+    "reference_frame": 1,
+    "constant_frames": [1, 2, 3],
+    "time_dependent_frames": [4, 5, 6]
+  },
+  "starting_detector_line": 0,
+  "starting_detector_sample": 0,
+  "detector_sample_summing": 1.0,
+  "detector_line_summing": 1.0,
+  "starting_ephemeris_time": 0,
+  "sun_position": {
+    "positions": [
+      [100000.0, 100000.0, 0.0]
+    ],
+    "velocities": [
+      [0.0, 0.0, 0.0]
+    ],
+    "ephemeris_times": [
+      50
+    ],
+    "unit": "m",
+    "reference_frame": 1
+  }
+}
-- 
GitLab