From c892d5dc2e88c4e9499a634ef00fd78b6a051960 Mon Sep 17 00:00:00 2001
From: acpaquette <acp263@nau.edu>
Date: Thu, 21 Feb 2019 15:28:31 -0700
Subject: [PATCH] Distortion Integration (#184)
* Moved transverse distortion from line scanner into new file
* Extracted radial destortion, and built associated tests
* Reconstructed distortion functions return results
* Corrected spelling, and removed commented out functions
* Renamed all dpoint variables to distortionPoint
* Removed computeUndistortedFocalPlaneCoordinates function
* Renamed distortionPoint to undistortedPoint and distortedPoint where applicable
* Removed doc string parameters
* Added combined new distortion model
* Half way through redefining parameters
* Major update to distortion and coefficient storage
* More debugging
* More debugging
* More debugging
* debugging
* debugging
* Removed debugging messages and set the output value for remove distortion
* Brought changes inline with dev
* Updated isd parsing to handle either radial or transverse
* Added Distortion Parse function
* Split apply and remove distortion, and added the distortion type to each model state
* Removed cmath
* Removing prints, and other fixes
* Updated how coefficients are parsed
* Reverted notebook
* Function namespace and parameter clean up
* More function and parameter clean up
* Reverted old code and updated distortion defaults
* Small changes for distortion
* Updated utilities and fixed failing test
---
CMakeLists.txt | 4 +-
include/usgscsm/Distortion.h | 31 +-
include/usgscsm/UsgsAstroFrameSensorModel.h | 5 +-
include/usgscsm/UsgsAstroLsSensorModel.h | 4 +-
include/usgscsm/Utilities.h | 6 +-
src/Distortion.cpp | 329 ++++++++++----------
src/UsgsAstroFrameSensorModel.cpp | 64 ++--
src/UsgsAstroLsSensorModel.cpp | 56 ++--
src/UsgsAstroPlugin.cpp | 1 -
src/Utilities.cpp | 115 ++++---
tests/DistortionTests.cpp | 123 ++++----
tests/ISDParsingTests.cpp | 22 +-
tests/PluginTests.cpp | 13 +-
tests/data/empty.json | 1 +
14 files changed, 393 insertions(+), 381 deletions(-)
create mode 100644 tests/data/empty.json
diff --git a/CMakeLists.txt b/CMakeLists.txt
index c987943..625de6a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,6 +1,6 @@
cmake_minimum_required(VERSION 3.10)
project(usgscsm VERSION 0.0.1 DESCRIPTION "usgscsm library")
-
+
include(GoogleTest)
include(cmake/gtest.cmake)
include(GNUInstallDirs)
@@ -58,7 +58,7 @@ install(DIRECTORY ${USGSCSM_INCLUDE_DIRS} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR
# Optional build or link against CSM
option (BUILD_TESTS "Build tests" ON)
if(BUILD_TESTS)
-
+
# Setup for GoogleTest
find_package (Threads)
diff --git a/include/usgscsm/Distortion.h b/include/usgscsm/Distortion.h
index 6ae0d04..574a970 100644
--- a/include/usgscsm/Distortion.h
+++ b/include/usgscsm/Distortion.h
@@ -4,22 +4,27 @@
#include <vector>
#include <math.h>
#include <tuple>
+#include <iostream>
-// Transverse Distortion
-std::tuple<double, double> removeDistortion(double dx, double dy,
- const std::vector<double> &odtX, const std::vector<double> &odtY);
-
-std::vector<std::vector<double>> distortionJacobian(double x, double y,
- const std::vector<double> &odtX, const std::vector<double> &odtY);
+enum DistortionType {
+ RADIAL,
+ TRANSVERSE
+};
-std::tuple<double, double> distortionFunction(double ux, double uy,
- const std::vector<double> &odtX, const std::vector<double> &odtY);
+// Transverse Distortion
+void distortionJacobian(double x, double y, double *jacobian,
+ const std::vector<double> opticalDistCoeffs);
-// Radial Distortion
-std::tuple<double, double> removeDistortion(double inFocalPlaneX, double inFocalPlaneY,
- const double opticalDistCoef[3], double tolerance = 1.0E-6);
+void computeTransverseDistortion(double ux, double uy, double &dx, double &dy,
+ const std::vector<double> opticalDistCoeffs);
-std::tuple<double, double> invertDistortion(double inFocalPlaneX, double inFocalPlaneY,
- const double opticalDistCoef[3], double desiredPrecision, double tolerance = 1.0E-6);
+void removeDistortion(double dx, double dy, double &ux, double &uy,
+ const std::vector<double> opticalDistCoeffs,
+ DistortionType distortionType,
+ const double tolerance = 1.0E-6);
+void applyDistortion(double ux, double uy, double &dx, double &dy,
+ const std::vector<double> opticalDistCoeffs,
+ DistortionType distortionType,
+ const double desiredPrecision = 0, const double tolerance = 1.0E-6);
#endif
diff --git a/include/usgscsm/UsgsAstroFrameSensorModel.h b/include/usgscsm/UsgsAstroFrameSensorModel.h
index 3a823e7..71c033f 100644
--- a/include/usgscsm/UsgsAstroFrameSensorModel.h
+++ b/include/usgscsm/UsgsAstroFrameSensorModel.h
@@ -9,6 +9,7 @@
#include "RasterGM.h"
#include "CorrelationModel.h"
#include "Distortion.h"
+#include "Utilities.h"
#include <json.hpp>
using json = nlohmann::json;
@@ -330,8 +331,8 @@ protected:
std::vector<double> m_currentParameterCovariance;
std::vector<csm::param::Type> m_parameterType;
std::vector<double> m_noAdjustments;
- std::vector<double> m_odtX;
- std::vector<double> m_odtY;
+ DistortionType m_distortionType;
+ std::vector<double> m_opticalDistCoeffs;
std::vector<double> m_transX;
std::vector<double> m_transY;
std::vector<double> m_spacecraftVelocity;
diff --git a/include/usgscsm/UsgsAstroLsSensorModel.h b/include/usgscsm/UsgsAstroLsSensorModel.h
index ceb9a15..5329ede 100644
--- a/include/usgscsm/UsgsAstroLsSensorModel.h
+++ b/include/usgscsm/UsgsAstroLsSensorModel.h
@@ -30,6 +30,7 @@
#include <RasterGM.h>
#include <SettableEllipsoid.h>
#include <CorrelationModel.h>
+#include "Distortion.h"
class UsgsAstroLsSensorModel : public csm::RasterGM, virtual public csm::SettableEllipsoid
{
@@ -78,7 +79,8 @@ public:
int m_ikCode;
double m_focalLength;
double m_zDirection;
- double m_opticalDistCoef[3];
+ DistortionType m_distortionType;
+ std::vector<double> m_opticalDistCoeffs;
double m_iTransS[3];
double m_iTransL[3];
double m_detectorSampleOrigin;
diff --git a/include/usgscsm/Utilities.h b/include/usgscsm/Utilities.h
index 49e84d5..719a724 100644
--- a/include/usgscsm/Utilities.h
+++ b/include/usgscsm/Utilities.h
@@ -1,6 +1,8 @@
#ifndef Utilities_h
#define Utilities_h
+#include "Distortion.h"
+
#include <vector>
#include <math.h>
#include <tuple>
@@ -75,8 +77,8 @@ double getMinHeight(nlohmann::json isd, csm::WarningList *list=nullptr);
double getMaxHeight(nlohmann::json isd, csm::WarningList *list=nullptr);
double getSemiMajorRadius(nlohmann::json isd, csm::WarningList *list=nullptr);
double getSemiMinorRadius(nlohmann::json isd, csm::WarningList *list=nullptr);
-std::vector<double> getTransverseDistortionX(nlohmann::json isd, csm::WarningList *list=nullptr);
-std::vector<double> getTransverseDistortionY(nlohmann::json isd, csm::WarningList *list=nullptr);
+DistortionType getDistortionModel(nlohmann::json isd, csm::WarningList *list=nullptr);
+std::vector<double> getDistortionCoeffs(nlohmann::json isd, csm::WarningList *list=nullptr);
std::vector<double> getRadialDistortion(nlohmann::json isd, csm::WarningList *list=nullptr);
std::vector<double> getSunPositions(nlohmann::json isd, csm::WarningList *list=nullptr);
std::vector<double> getSensorPositions(nlohmann::json isd, csm::WarningList *list=nullptr);
diff --git a/src/Distortion.cpp b/src/Distortion.cpp
index 103aa5c..7a43ceb 100644
--- a/src/Distortion.cpp
+++ b/src/Distortion.cpp
@@ -1,82 +1,5 @@
#include "Distortion.h"
-/**
- * @brief Compute undistorted focal plane x/y.
- *
- * Computes undistorted focal plane (x,y) coordinates given a distorted focal plane (x,y)
- * coordinate. The undistorted coordinates are solved for using the Newton-Raphson
- * method for root-finding if the distortionFunction method is invoked.
- *
- * @param dx distorted focal plane x in millimeters
- * @param dy distorted focal plane y in millimeters
- * @param undistortedX The undistorted x coordinate, in millimeters.
- * @param undistortedY The undistorted y coordinate, in millimeters.
- *
- * @return if the conversion was successful
- * @todo Review the tolerance and maximum iterations of the root-
- * finding algorithm.
- * @todo Review the handling of non-convergence of the root-finding
- * algorithm.
- * @todo Add error handling for near-zero determinant.
-*/
-std::tuple<double, double> removeDistortion(double dx, double dy,
- const std::vector<double> &odtX, const std::vector<double> &odtY) {
- // Solve the distortion equation using the Newton-Raphson method.
- // Set the error tolerance to about one millionth of a NAC pixel.
- const double tol = 1.4E-5;
-
- // The maximum number of iterations of the Newton-Raphson method.
- const int maxTries = 60;
-
- double x;
- double y;
- std::tuple<double, double> undistortedPoint(dx, dy);
- std::tuple<double, double> distortedPoint;
-
- // Initial guess at the root
- x = dx;
- y = dy;
-
- distortedPoint = distortionFunction(x, y, odtX, odtY);
-
- for (int count = 1; ((fabs(std::get<0>(distortedPoint)) +fabs(std::get<1>(distortedPoint))) > tol) && (count < maxTries); count++) {
-
- distortedPoint = distortionFunction(x, y, odtX, odtY);
-
- // fx = dx - fx;
- // fy = dy - fy;
- distortedPoint = std::make_tuple(dx - std::get<0>(distortedPoint), dy - std::get<1>(distortedPoint));
-
- std::vector<std::vector<double>> jacobian;
-
- jacobian = distortionJacobian(x, y, odtX, odtY);
-
- // Jxx * Jyy - Jxy * Jyx
- double determinant = jacobian[0][0] * jacobian[1][1] - jacobian[0][1] * jacobian[1][0];
- if (fabs(determinant) < 1E-6) {
- undistortedPoint = std::make_tuple(x, y);
- //
- // Near-zero determinant. Add error handling here.
- //
- //-- Just break out and return with no convergence
- return undistortedPoint;
- }
-
- //x = x + (Jyy * fx - Jxy * fy)
- x = x + (jacobian[1][1] * std::get<0>(distortedPoint) - jacobian[0][1] * std::get<1>(distortedPoint)) / determinant;
- // y = y + (Jxx * fy - Jyx * fx)
- y = y + (jacobian[0][0] * std::get<1>(distortedPoint) - jacobian[1][0] * std::get<0>(distortedPoint)) / determinant;
- }
-
- if ( (fabs(std::get<0>(distortedPoint)) + fabs(std::get<1>(distortedPoint))) <= tol) {
- // The method converged to a root.
- undistortedPoint = std::make_tuple(x, y);
- }
- // Otherwise method did not converge to a root within the maximum
- // number of iterations. Return with no distortion.
- return undistortedPoint;
-}
-
/**
* @description Jacobian of the distortion function. The Jacobian was computed
* algebraically from the function described in the distortionFunction
@@ -92,8 +15,8 @@ std::tuple<double, double> removeDistortion(double dx, double dy,
[1][0]: yx, [1][1]: yy
*/
-std::vector<std::vector<double>> distortionJacobian(double x, double y,
- const std::vector<double> &odtX, const std::vector<double> &odtY) {
+void distortionJacobian(double x, double y, double *jacobian,
+ const std::vector<double> opticalDistCoeffs) {
double d_dx[10];
d_dx[0] = 0;
@@ -118,21 +41,20 @@ std::vector<std::vector<double>> distortionJacobian(double x, double y,
d_dy[8] = 2 * x * y;
d_dy[9] = 3 * y * y;
- std::vector<std::vector<double>> jacobian(2, std::vector<double>(2));
+ jacobian[0] = 0; // xx
+ jacobian[1] = 0; // xy
+ jacobian[2] = 0; // yx
+ jacobian[3] = 0; // yy
- jacobian[0][0] = 0;
- jacobian[0][1] = 0;
- jacobian[1][0] = 0;
- jacobian[1][1] = 0;
+ int xPointer = 0;
+ int yPointer = opticalDistCoeffs.size() / 2;
- for (int i = 0; i < 10; i++) {
- jacobian[0][0] = jacobian[0][0] + d_dx[i] * odtX[i];
- jacobian[0][1] = jacobian[0][1] + d_dy[i] * odtX[i];
- jacobian[1][0] = jacobian[1][0] + d_dx[i] * odtY[i];
- jacobian[1][1] = jacobian[1][1] + d_dy[i] * odtY[i];
+ for (int i = 0; i < 10; xPointer++, yPointer++, i++) {
+ jacobian[0] = jacobian[0] + d_dx[i] * opticalDistCoeffs[xPointer];
+ jacobian[1] = jacobian[1] + d_dy[i] * opticalDistCoeffs[xPointer];
+ jacobian[2] = jacobian[2] + d_dx[i] * opticalDistCoeffs[yPointer];
+ jacobian[3] = jacobian[3] + d_dy[i] * opticalDistCoeffs[yPointer];
}
-
- return jacobian;
}
/**
@@ -142,13 +64,12 @@ std::vector<std::vector<double>> distortionJacobian(double x, double y,
*
* @param ux Undistored x
* @param uy Undistored y
- * @param odtX opticalDistCoef In X
- * @param odtY opticalDistCoef In Y
+ * @param opticalDistCoeffs For both X and Y coefficients
*
* @returns distortedPoint Newly adjusted focal plane coordinates as an x, y tuple
*/
-std::tuple<double, double> distortionFunction(double ux, double uy,
- const std::vector<double> &odtX, const std::vector<double> &odtY) {
+void computeTransverseDistortion(double ux, double uy, double &dx, double &dy,
+ const std::vector<double> opticalDistCoeffs) {
double f[10];
f[0] = 1;
@@ -162,93 +83,155 @@ std::tuple<double, double> distortionFunction(double ux, double uy,
f[8] = ux * uy * uy;
f[9] = uy * uy * uy;
- std::tuple<double, double> distortedPoint(0.0, 0.0);
- for (int i = 0; i < 10; i++) {
- distortedPoint = std::make_tuple(std::get<0>(distortedPoint) + f[i] * odtX[i],
- std::get<1>(distortedPoint) + f[i] * odtY[i]);
- }
+ int xPointer = 0;
+ int yPointer = opticalDistCoeffs.size() / 2;
- return distortedPoint;
-}
+ dx = 0.0;
+ dy = 0.0;
-/**
- * @description Compute undistorted focal plane coordinate given a distorted
- * coordinate set and the distortion coefficients
- *
- * @param inFocalPlaneX Distorted x
- * @param inFocalPlaneY Distorted y
- * @param opticalDistCoef distortion coefficients
- *
- * @returns undistortedPoint Newly adjusted focal plane coordinates as an x, y tuple
- */
-std::tuple<double, double> removeDistortion(double inFocalPlaneX, double inFocalPlaneY,
- const double opticalDistCoef[3], double tolerance) {
- double rr = inFocalPlaneX * inFocalPlaneX + inFocalPlaneY * inFocalPlaneY;
- std::tuple<double, double> undistortedPoint(inFocalPlaneX, inFocalPlaneY);
-
- if (rr > tolerance)
- {
- double dr = opticalDistCoef[0] + (rr * (opticalDistCoef[1] + rr * opticalDistCoef[2]));
- undistortedPoint = std::make_tuple(inFocalPlaneX * (1.0 - dr), inFocalPlaneY * (1.0 - dr));
+ for (int i = 0; i < 10; xPointer++, yPointer++, i++) {
+ dx = dx + f[i] * opticalDistCoeffs[xPointer];
+ dy = dy + f[i] * opticalDistCoeffs[yPointer];
}
-
- return undistortedPoint;
}
-/**
- * @description Compute undistorted focal plane coordinate given a distorted
- * focal plane coordinate. This method works by iteratively adding distortion
- * until the new distorted point, r, undistorts to within a tolerance of the
- * original point, rp.
- *
- * @param inFocalPlaneX Distorted x
- * @param inFocalPlaneY Distorted y
- * @param opticalDistCoef Distortion coefficients
- * @param desiredPrecision Convergence precision
- * @param tolerance Tolerance of r^2
- *
- * @returns undistortedPoint Newly adjusted focal plane coordinates as an x, y tuple
- */
-std::tuple<double, double> invertDistortion(double inFocalPlaneX, double inFocalPlaneY,
- const double opticalDistCoef[3], double desiredPrecision, double tolerance) {
- double rp2 = (inFocalPlaneX * inFocalPlaneX) +
- (inFocalPlaneY * inFocalPlaneY);
- std::tuple<double, double> undistortedPoint;
-
- if (rp2 > tolerance) {
- double rp = sqrt(rp2);
- // Compute first fractional distortion using rp
- double drOverR = opticalDistCoef[0]
- + (rp2 * (opticalDistCoef[1] + (rp2 * opticalDistCoef[2])));
- // Compute first distorted point estimate, r
- double r = rp + (drOverR * rp);
- double r_prev, r2_prev;
- int iteration = 0;
- do {
- // Don't get in an end-less loop. This algorithm should
- // converge quickly. If not then we are probably way outside
- // of the focal plane. Just set the distorted position to the
- // undistorted position. Also, make sure the focal plane is less
- // than 1km, it is unreasonable for it to grow larger than that.
- if (iteration >= 15 || r > 1E9) {
- drOverR = 0.0;
- break;
+void removeDistortion(double dx, double dy, double &ux, double &uy,
+ const std::vector<double> opticalDistCoeffs,
+ DistortionType distortionType,
+ const double tolerance) {
+ ux = dx;
+ uy = dy;
+
+ switch (distortionType) {
+ // Compute undistorted focal plane coordinate given a distorted
+ // coordinate set and the distortion coefficients
+ case RADIAL: {
+ double rr = dx * dx + dy * dy;
+
+ if (rr > tolerance)
+ {
+ double dr = opticalDistCoeffs[0] + (rr * (opticalDistCoeffs[1] + rr * opticalDistCoeffs[2]));
+ ux = dx * (1.0 - dr);
+ uy = dy * (1.0 - dr);
+ }
+ }
+ break;
+ // Computes undistorted focal plane (x,y) coordinates given a distorted focal plane (x,y)
+ // coordinate. The undistorted coordinates are solved for using the Newton-Raphson
+ // method for root-finding if the distortionFunction method is invoked.
+ case TRANSVERSE: {
+ // Solve the distortion equation using the Newton-Raphson method.
+ // Set the error tolerance to about one millionth of a NAC pixel.
+ // The maximum number of iterations of the Newton-Raphson method.
+ const int maxTries = 20;
+
+ double x;
+ double y;
+ double fx;
+ double fy;
+ double jacobian[4];
+
+ // Initial guess at the root
+ x = dx;
+ y = dy;
+
+ computeTransverseDistortion(x, y, fx, fy, opticalDistCoeffs);
+
+ for (int count = 1; ((fabs(fx) +fabs(fy)) > tolerance) && (count < maxTries); count++) {
+
+ computeTransverseDistortion(x, y, fx, fy, opticalDistCoeffs);
+
+ fx = dx - fx;
+ fy = dy - fy;
+
+ distortionJacobian(x, y, jacobian, opticalDistCoeffs);
+
+ // Jxx * Jyy - Jxy * Jyx
+ double determinant = jacobian[0] * jacobian[3] - jacobian[1] * jacobian[2];
+ if (fabs(determinant) < 1E-6) {
+ ux = x;
+ uy = y;
+ //
+ // Near-zero determinant. Add error handling here.
+ //
+ //-- Just break out and return with no convergence
+ return;
+ }
+
+ x = x + (jacobian[3] * fx - jacobian[1] * fy) / determinant;
+ y = y + (jacobian[0] * fy - jacobian[2] * fx) / determinant;
}
- r_prev = r;
- r2_prev = r * r;
+ if ((fabs(fx) + fabs(fy)) <= tolerance) {
+ // The method converged to a root.
+ ux = x;
+ uy = y;
- // Compute new fractional distortion:
- drOverR = opticalDistCoef[0]
- + (r2_prev * (opticalDistCoef[1] + (r2_prev * opticalDistCoef[2])));
+ return;
+ }
+ // Otherwise method did not converge to a root within the maximum
+ // number of iterations
+ }
+ break;
+ }
+}
- // Compute new estimate of r
- r = rp + (drOverR * r_prev);
- iteration++;
+void applyDistortion(double ux, double uy, double &dx, double &dy,
+ const std::vector<double> opticalDistCoeffs,
+ DistortionType distortionType,
+ const double desiredPrecision, const double tolerance)
+{
+ dx = ux;
+ dy = uy;
+
+ switch (distortionType) {
+ // 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
+ // original point, rp.
+ case RADIAL: {
+ double rp2 = (ux * ux) + (uy * uy);
+
+ if (rp2 > tolerance) {
+ double rp = sqrt(rp2);
+ // Compute first fractional distortion using rp
+ double drOverR = opticalDistCoeffs[0]
+ + (rp2 * (opticalDistCoeffs[1] + (rp2 * opticalDistCoeffs[2])));
+ // Compute first distorted point estimate, r
+ double r = rp + (drOverR * rp);
+ double r_prev, r2_prev;
+ int iteration = 0;
+ do {
+ // Don't get in an end-less loop. This algorithm should
+ // converge quickly. If not then we are probably way outside
+ // of the focal plane. Just set the distorted position to the
+ // undistorted position. Also, make sure the focal plane is less
+ // than 1km, it is unreasonable for it to grow larger than that.
+ if (iteration >= 15 || r > 1E9) {
+ drOverR = 0.0;
+ break;
+ }
+
+ r_prev = r;
+ r2_prev = r * r;
+
+ // Compute new fractional distortion:
+ drOverR = opticalDistCoeffs[0]
+ + (r2_prev * (opticalDistCoeffs[1] + (r2_prev * opticalDistCoeffs[2])));
+
+ // Compute new estimate of r
+ r = rp + (drOverR * r_prev);
+ iteration++;
+ }
+ while (fabs(r * (1 - drOverR) - rp) > desiredPrecision);
+ dx = ux / (1.0 - drOverR);
+ dy = uy / (1.0 - drOverR);
+ }
+ }
+ break;
+ case TRANSVERSE: {
+ computeTransverseDistortion(ux, uy, dx, dy, opticalDistCoeffs);
}
- while (fabs(r * (1 - drOverR) - rp) > desiredPrecision);
- undistortedPoint = std::make_tuple(inFocalPlaneX / (1.0 - drOverR),
- inFocalPlaneY / (1.0 - drOverR));
+ break;
}
- return undistortedPoint;
}
diff --git a/src/UsgsAstroFrameSensorModel.cpp b/src/UsgsAstroFrameSensorModel.cpp
index 39a1523..843a236 100644
--- a/src/UsgsAstroFrameSensorModel.cpp
+++ b/src/UsgsAstroFrameSensorModel.cpp
@@ -9,9 +9,6 @@
#include <Error.h>
#include <Version.h>
-#include "Distortion.h"
-#include "Utilities.h"
-
using json = nlohmann::json;
using namespace std;
@@ -64,8 +61,8 @@ void UsgsAstroFrameSensorModel::reset() {
m_ccdCenter = std::vector<double>(2, 0.0);
m_spacecraftVelocity = std::vector<double>(3, 0.0);
m_sunPosition = std::vector<double>(3, 0.0);
- m_odtX = std::vector<double>(10, 0.0);
- m_odtY = std::vector<double>(10, 0.0);
+ m_distortionType = DistortionType::TRANSVERSE;
+ m_opticalDistCoeffs = std::vector<double>(20, 0.0);
m_transX = std::vector<double>(3, 0.0);
m_transY = std::vector<double>(3, 0.0);
m_iTransS = std::vector<double>(3, 0.0);
@@ -128,7 +125,7 @@ csm::ImageCoord UsgsAstroFrameSensorModel::groundToImage(
// Camera rotation matrix
double m[3][3];
- calcRotationMatrix(m, adjustments);
+ calcRotationMatrix(m,adjustments);
// Sensor position
double undistortedx, undistortedy, denom;
@@ -137,14 +134,15 @@ csm::ImageCoord UsgsAstroFrameSensorModel::groundToImage(
undistortedy = (f * (m[0][1] * xo + m[1][1] * yo + m[2][1] * zo)/denom) + m_linePp;
// Apply the distortion to the line/sample location and then convert back to line/sample
- std::tuple<double, double> distortedPoint;
- distortedPoint = distortionFunction(undistortedx, undistortedy, m_odtX, m_odtY);
+ double distortedX, distortedY;
+ applyDistortion(undistortedx, undistortedy, distortedX, distortedY,
+ m_opticalDistCoeffs, m_distortionType);
// Convert distorted mm into line/sample
double sample, line;
- sample = m_iTransS[0] + m_iTransS[1] * std::get<0>(distortedPoint) + m_iTransS[2] * std::get<1>(distortedPoint) + m_ccdCenter[1];
- line = m_iTransL[0] + m_iTransL[1] * std::get<0>(distortedPoint) + m_iTransL[2] * std::get<1>(distortedPoint) + m_ccdCenter[0];
+ sample = m_iTransS[0] + m_iTransS[1] * distortedX + m_iTransS[2] * distortedY + m_ccdCenter[1];
+ line = m_iTransL[0] + m_iTransL[1] * distortedX + m_iTransL[2] * distortedY + m_ccdCenter[0];
return csm::ImageCoord(line, sample);
}
@@ -195,13 +193,15 @@ csm::EcefCoord UsgsAstroFrameSensorModel::imageToGround(const csm::ImageCoord &i
x_camera = m_transX[0] + m_transX[1] * (lo - line_center) + m_transX[2] * (so - sample_center);
// Apply the distortion model (remove distortion)
- std::tuple<double, double> undistortedPoint;
- undistortedPoint = removeDistortion(x_camera, y_camera, m_odtX, m_odtY);
+ double undistortedX, undistortedY;
+ removeDistortion(x_camera, y_camera, undistortedX, undistortedY,
+ m_opticalDistCoeffs,
+ m_distortionType);
// Now back from distorted mm to pixels
- xl = m[0][0] * std::get<0>(undistortedPoint) + m[0][1] * std::get<1>(undistortedPoint) - m[0][2] * - m_focalLength;
- yl = m[1][0] * std::get<0>(undistortedPoint) + m[1][1] * std::get<1>(undistortedPoint) - m[1][2] * - m_focalLength;
- zl = m[2][0] * std::get<0>(undistortedPoint) + m[2][1] * std::get<1>(undistortedPoint) - m[2][2] * - m_focalLength;
+ xl = m[0][0] * undistortedX + m[0][1] * undistortedY - m[0][2] * - m_focalLength;
+ yl = m[1][0] * undistortedX + m[1][1] * undistortedY - m[1][2] * - m_focalLength;
+ zl = m[2][0] * undistortedX + m[2][1] * undistortedY - m[2][2] * - m_focalLength;
double x, y, z;
double xc, yc, zc;
@@ -249,13 +249,16 @@ csm::EcefLocus UsgsAstroFrameSensorModel::imageToRemoteImagingLocus(const csm::I
double focalPlaneY = m_transY[0] + m_transY[1] * row + m_transY[2] * col;
// Distort
- std::tuple<double, double> undistortedPoint;
- undistortedPoint = removeDistortion(focalPlaneX, focalPlaneY, m_odtX, m_odtY);
+ double undistortedFocalPlaneX, undistortedFocalPlaneY;
+ removeDistortion(focalPlaneX, focalPlaneY,
+ undistortedFocalPlaneX, undistortedFocalPlaneY,
+ m_opticalDistCoeffs,
+ m_distortionType);
// Get rotation matrix and transform to a body-fixed frame
double m[3][3];
calcRotationMatrix(m);
- std::vector<double> lookC { std::get<0>(undistortedPoint), std::get<1>(undistortedPoint), m_focalLength };
+ std::vector<double> lookC { undistortedFocalPlaneX, undistortedFocalPlaneY, m_focalLength };
std::vector<double> lookB {
m[0][0] * lookC[0] + m[0][1] * lookC[1] + m[0][2] * lookC[2],
m[1][0] * lookC[0] + m[1][1] * lookC[1] + m[1][2] * lookC[2],
@@ -636,8 +639,8 @@ std::string UsgsAstroFrameSensorModel::getModelState() const {
{"m_samplePp", m_samplePp},
{"m_minElevation", m_minElevation},
{"m_maxElevation", m_maxElevation},
- {"m_odtX", m_odtX},
- {"m_odtY", m_odtY},
+ {"m_distortionType", m_distortionType},
+ {"m_opticalDistCoeffs", m_opticalDistCoeffs},
{"m_originalHalfLines", m_originalHalfLines},
{"m_originalHalfSamples", m_originalHalfSamples},
{"m_spacecraftName", m_spacecraftName},
@@ -671,8 +674,8 @@ bool UsgsAstroFrameSensorModel::isValidModelState(const std::string& stringState
"m_ccdCenter",
"m_spacecraftVelocity",
"m_sunPosition",
- "m_odtX",
- "m_odtY",
+ "m_distortionType",
+ "m_opticalDistCoeffs",
"m_transX",
"m_transY",
"m_iTransS",
@@ -745,8 +748,8 @@ void UsgsAstroFrameSensorModel::replaceModelState(const std::string& stringState
m_ccdCenter = state.at("m_ccdCenter").get<std::vector<double>>();
m_spacecraftVelocity = state.at("m_spacecraftVelocity").get<std::vector<double>>();
m_sunPosition = state.at("m_sunPosition").get<std::vector<double>>();
- m_odtX = state.at("m_odtX").get<std::vector<double>>();
- m_odtY = state.at("m_odtY").get<std::vector<double>>();
+ m_distortionType = (DistortionType)state.at("m_distortionType").get<int>();
+ m_opticalDistCoeffs = state.at("m_opticalDistCoeffs").get<std::vector<double>>();
m_transX = state.at("m_transX").get<std::vector<double>>();
m_transY = state.at("m_transY").get<std::vector<double>>();
m_iTransS = state.at("m_iTransS").get<std::vector<double>>();
@@ -793,7 +796,6 @@ void UsgsAstroFrameSensorModel::replaceModelState(const std::string& stringState
std::string UsgsAstroFrameSensorModel::constructStateFromIsd(const std::string& jsonIsd, csm::WarningList* warnings) {
-
json isd = json::parse(jsonIsd);
json state = {};
@@ -807,7 +809,6 @@ std::string UsgsAstroFrameSensorModel::constructStateFromIsd(const std::string&
state["m_startingDetectorSample"] = getDetectorStartingSample(isd, parsingWarnings);
state["m_startingDetectorLine"] = getDetectorStartingLine(isd, parsingWarnings);
-
// get focal length
state["m_focalLength"] = getFocalLength(isd, parsingWarnings);
state["m_focalLengthEpsilon"] = getFocalLengthEpsilon(isd, parsingWarnings);
@@ -849,10 +850,9 @@ std::string UsgsAstroFrameSensorModel::constructStateFromIsd(const std::string&
// sun position is not strictly necessary, but is required for getIlluminationDirection.
state["m_sunPosition"] = getSunPositions(isd);
-
// get sensor_orientation quaternion
std::vector<double> quaternion = getSensorOrientations(isd, parsingWarnings);
- if (!quaternion.empty() && quaternion.size() != 4) {
+ if (quaternion.size() != 4) {
parsingWarnings->push_back(
csm::Warning(
csm::Warning::DATA_NOT_AVAILABLE,
@@ -866,11 +866,9 @@ std::string UsgsAstroFrameSensorModel::constructStateFromIsd(const std::string&
state["m_currentParameterValue"][6] = quaternion[3];
}
-
// get optical_distortion
- state["m_odtX"] = getTransverseDistortionX(isd, parsingWarnings);
- state["m_odtY"] = getTransverseDistortionY(isd, parsingWarnings);
-
+ state["m_distortionType"] = getDistortionModel(isd);
+ state["m_opticalDistCoeffs"] = getDistortionCoeffs(isd);
// get detector_center
state["m_ccdCenter"][0] = getDetectorCenterLine(isd, parsingWarnings);
@@ -939,10 +937,10 @@ std::string UsgsAstroFrameSensorModel::constructStateFromIsd(const std::string&
parsingWarnings = nullptr;
return state.dump();
-
}
+
csm::EcefCoord UsgsAstroFrameSensorModel::getReferencePoint() const {
return m_referencePointXyz;
}
diff --git a/src/UsgsAstroLsSensorModel.cpp b/src/UsgsAstroLsSensorModel.cpp
index c998736..ec3d030 100644
--- a/src/UsgsAstroLsSensorModel.cpp
+++ b/src/UsgsAstroLsSensorModel.cpp
@@ -77,7 +77,8 @@ const std::string UsgsAstroLsSensorModel::_STATE_KEYWORD[] =
"m_ikCode",
"m_focalLength",
"m_zDirection",
- "m_opticalDistCoef",
+ "m_distortionType",
+ "m_opticalDistCoeffs",
"m_iTransS",
"m_iTransL",
"m_detectorSampleOrigin",
@@ -151,9 +152,9 @@ void UsgsAstroLsSensorModel::replaceModelState(const std::string &stateString )
m_focalLength = j["m_focalLength"];
m_zDirection = j["m_zDirection"];
-
+ m_distortionType = (DistortionType)j["m_distortionType"].get<int>();
+ m_opticalDistCoeffs = j["m_opticalDistCoeffs"].get<std::vector<double>>();
for (int i = 0; i < 3; i++) {
- m_opticalDistCoef[i] = j["m_opticalDistCoef"][i];
m_iTransS[i] = j["m_iTransS"][i];
m_iTransL[i] = j["m_iTransL"][i];
}
@@ -258,7 +259,8 @@ std::string UsgsAstroLsSensorModel::getModelState() const {
state["m_ikCode"] = m_ikCode;
state["m_focalLength"] = m_focalLength;
state["m_zDirection"] = m_zDirection;
- state["m_opticalDistCoef"] = std::vector<double>(m_opticalDistCoef, m_opticalDistCoef+3);
+ state["m_distortionType"] = m_distortionType;
+ state["m_opticalDistCoeffs"] = m_opticalDistCoeffs;
state["m_iTransS"] = std::vector<double>(m_iTransS, m_iTransS+3);
state["m_iTransL"] = std::vector<double>(m_iTransL, m_iTransL+3);
state["m_detectorSampleOrigin"] = m_detectorSampleOrigin;
@@ -323,9 +325,8 @@ void UsgsAstroLsSensorModel::reset()
m_ikCode = -85600; // 17
m_focalLength = 350.0; // 18
m_zDirection = 1.0; // 19
- m_opticalDistCoef[0] = 0.0; // 20
- m_opticalDistCoef[1] = 0.0; // 20
- m_opticalDistCoef[2] = 0.0; // 20
+ m_distortionType = DistortionType::RADIAL;
+ m_opticalDistCoeffs = std::vector<double>(3, 0.0);
m_iTransS[0] = 0.0; // 21
m_iTransS[1] = 0.0; // 21
m_iTransS[2] = 150.0; // 21
@@ -353,13 +354,13 @@ void UsgsAstroLsSensorModel::reset()
m_currentParameterValue.assign(NUM_PARAMETERS,0.0);
m_parameterType.assign(NUM_PARAMETERS,csm::param::REAL);
- m_referencePointXyz.x = 0.0; // 47
- m_referencePointXyz.y = 0.0; // 47
- m_referencePointXyz.z = 0.0; // 47
- m_gsd = 1.0; // 48
- m_flyingHeight = 1000.0; // 49
- m_halfSwath = 1000.0; // 50
- m_halfTime = 10.0; // 51
+ m_referencePointXyz.x = 0.0;
+ m_referencePointXyz.y = 0.0;
+ m_referencePointXyz.z = 0.0;
+ m_gsd = 1.0;
+ m_flyingHeight = 1000.0;
+ m_halfSwath = 1000.0;
+ m_halfTime = 10.0;
m_covariance = std::vector<double>(NUM_PARAMETERS * NUM_PARAMETERS,0.0); // 52
}
@@ -1660,12 +1661,15 @@ void UsgsAstroLsSensorModel::losToEcf(
computeDistortedFocalPlaneCoordinates(fractionalLine, sampleUSGSFull, m_detectorSampleOrigin, m_detectorLineOrigin, m_detectorSampleSumming, m_startingSample, m_iTransS, m_iTransL, natFocalPlane);
// Remove lens distortion
- std::tuple<double, double> undistortedPoint;
- undistortedPoint = removeDistortion(std::get<0>(natFocalPlane), std::get<1>(natFocalPlane), m_opticalDistCoef);
+ double undistortedFocalPlaneX, undistortedFocalPlaneY;
+ removeDistortion(std::get<0>(natFocalPlane), std::get<1>(natFocalPlane),
+ undistortedFocalPlaneX, undistortedFocalPlaneY,
+ m_opticalDistCoeffs,
+ DistortionType::RADIAL);
// Define imaging ray (look vector) in camera space
double cameraLook[3];
- createCameraLookVector(std::get<0>(undistortedPoint), std::get<1>(undistortedPoint), m_zDirection, m_focalLength, cameraLook);
+ createCameraLookVector(undistortedFocalPlaneX, undistortedFocalPlaneY, m_zDirection, m_focalLength, cameraLook);
// Apply attitude correction
double attCorr[9];
@@ -2249,18 +2253,19 @@ csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel(
double lookScale = m_focalLength / adjustedLookZ;
double focalX = adjustedLookX * lookScale;
double focalY = adjustedLookY * lookScale;
- std::tuple<double, double> undistortedPoint;
+ double distortedFocalX, distortedFocalY;
// Invert distortion
- undistortedPoint = invertDistortion(focalX, focalY, m_opticalDistCoef, desiredPrecision);
+ applyDistortion(focalX, focalY, distortedFocalX, distortedFocalY,
+ m_opticalDistCoeffs, m_distortionType, desiredPrecision);
// Convert to detector line and sample
double detectorLine = m_iTransL[0]
- + m_iTransL[1] * std::get<0>(undistortedPoint)
- + m_iTransL[2] * std::get<1>(undistortedPoint);
+ + m_iTransL[1] * distortedFocalX
+ + m_iTransL[2] * distortedFocalY;
double detectorSample = m_iTransS[0]
- + m_iTransS[1] * std::get<0>(undistortedPoint)
- + m_iTransS[2] * std::get<1>(undistortedPoint);
+ + m_iTransS[1] * distortedFocalX
+ + m_iTransS[2] * distortedFocalY;
// Convert to image sample line
double line = detectorLine + m_detectorLineOrigin;
@@ -2442,7 +2447,6 @@ double UsgsAstroLsSensorModel::determinant3x3(double mat[9]) const
std::string UsgsAstroLsSensorModel::constructStateFromIsd(const std::string imageSupportData, csm::WarningList *warnings) const
{
-
// Instantiate UsgsAstroLineScanner sensor model
json isd = json::parse(imageSupportData);
json state = {};
@@ -2468,6 +2472,9 @@ std::string UsgsAstroLsSensorModel::constructStateFromIsd(const std::string imag
state["m_ikCode"] = 0;
state["m_zDirection"] = 1;
+ state["m_distortionType"] = getDistortionModel(isd, parsingWarnings);
+ state["m_opticalDistCoeffs"] = getDistortionCoeffs(isd, parsingWarnings);
+
// Zero computed state values
state["m_referencePointXyz"] = std::vector<double>(3, 0.0);
@@ -2487,7 +2494,6 @@ std::string UsgsAstroLsSensorModel::constructStateFromIsd(const std::string imag
state["m_detectorSampleSumming"] = getSampleSumming(isd, parsingWarnings);
state["m_startingDetectorSample"] = getDetectorStartingSample(isd, parsingWarnings);
state["m_startingDetectorLine"] = getDetectorStartingLine(isd, parsingWarnings);
- state["m_opticalDistCoef"] = getRadialDistortion(isd, parsingWarnings);
state["m_detectorSampleOrigin"] = getDetectorCenterSample(isd, parsingWarnings);
state["m_detectorLineOrigin"] = getDetectorCenterLine(isd, parsingWarnings);
diff --git a/src/UsgsAstroPlugin.cpp b/src/UsgsAstroPlugin.cpp
index d731f6c..cf6e017 100644
--- a/src/UsgsAstroPlugin.cpp
+++ b/src/UsgsAstroPlugin.cpp
@@ -217,7 +217,6 @@ std::string UsgsAstroPlugin::convertISDToModelState(const csm::Isd &imageSupport
csm::Model *UsgsAstroPlugin::constructModelFromISD(const csm::Isd &imageSupportDataOriginal,
const std::string &modelName,
csm::WarningList *warnings) const {
-
std::string stringIsd = loadImageSupportData(imageSupportDataOriginal);
if (modelName == UsgsAstroFrameSensorModel::_SENSOR_MODEL_NAME) {
diff --git a/src/Utilities.cpp b/src/Utilities.cpp
index a8fc1cf..1812d66 100644
--- a/src/Utilities.cpp
+++ b/src/Utilities.cpp
@@ -551,56 +551,89 @@ double getSemiMinorRadius(json isd, csm::WarningList *list) {
return radius;
}
-std::vector<double> getTransverseDistortionX(json isd, csm::WarningList *list) {
- std::vector<double> coefficients;
- try {
- coefficients = isd.at("optical_distortion").at("transverse").at("x").get<std::vector<double>>();
- coefficients.resize(10, 0.0);
+DistortionType getDistortionModel(json isd, csm::WarningList *list) {
+ json distoriton_subset = isd.at("optical_distortion");
+
+ json::iterator it = distoriton_subset.begin();
+
+ std::string distortion = (std::string)it.key();
+
+ if (distortion.compare("transverse") == 0) {
+ return DistortionType::TRANSVERSE;
}
- catch (...) {
- if (list) {
- list->push_back(
- csm::Warning(
- csm::Warning::DATA_NOT_AVAILABLE,
- "Could not parse the transverse distortion model X coefficients.",
- "Utilities::getTransverseDistortionX()"));
- }
+ else if (distortion.compare("radial") == 0) {
+ return DistortionType::RADIAL;
}
- return coefficients;
-}
-std::vector<double> getTransverseDistortionY(json isd, csm::WarningList *list) {
- std::vector<double> coefficients;
- try {
- coefficients = isd.at("optical_distortion").at("transverse").at("y").get<std::vector<double>>();
- coefficients.resize(10, 0.0);
- }
- catch (...) {
- if (list) {
- list->push_back(
- csm::Warning(
- csm::Warning::DATA_NOT_AVAILABLE,
- "Could not parse the transverse distortion model Y coefficients.",
- "Utilities::getTransverseDistortionY()"));
- }
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the distortion model.",
+ "Utilities::getDistortionModel()"));
}
- return coefficients;
+ return DistortionType::TRANSVERSE;
}
-std::vector<double> getRadialDistortion(json isd, csm::WarningList *list) {
+std::vector<double> getDistortionCoeffs(json isd, csm::WarningList *list) {
std::vector<double> coefficients;
- try {
- coefficients = isd.at("optical_distortion").at("radial").at("coefficients").get<std::vector<double>>();
- }
- catch (...) {
- if (list) {
- list->push_back(
- csm::Warning(
- csm::Warning::DATA_NOT_AVAILABLE,
- "Could not parse the radial distortion model coefficients.",
- "Utilities::getRadialDistortion()"));
+
+ DistortionType distortion = getDistortionModel(isd);
+
+ switch (distortion) {
+ case DistortionType::TRANSVERSE: {
+ try {
+ std::vector<double> coefficientsX, coefficientsY;
+
+ coefficientsX = isd.at("optical_distortion").at("transverse").at("x").get<std::vector<double>>();
+ coefficientsX.resize(10, 0.0);
+
+ coefficientsY = isd.at("optical_distortion").at("transverse").at("y").get<std::vector<double>>();
+ coefficientsY.resize(10, 0.0);
+
+ coefficients = coefficientsX;
+
+ coefficients.insert(coefficients.end(), coefficientsY.begin(), coefficientsY.end());
+ return coefficients;
+ }
+ catch (...) {
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse a set of transverse distortion model coefficients.",
+ "Utilities::getDistortion()"));
+ }
+ coefficients = std::vector<double>(20, 0.0);
+ }
+ }
+ break;
+ case DistortionType::RADIAL: {
+ try {
+ coefficients = isd.at("optical_distortion").at("radial").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 radial distortion model coefficients.",
+ "Utilities::getDistortion()"));
+ }
+ coefficients = std::vector<double>(3, 0.0);
+ }
}
}
+ if (list) {
+ list->push_back(
+ csm::Warning(
+ csm::Warning::DATA_NOT_AVAILABLE,
+ "Could not parse the distortion model coefficients.",
+ "Utilities::getDistortion()"));
+ }
+
return coefficients;
}
diff --git a/tests/DistortionTests.cpp b/tests/DistortionTests.cpp
index b21971e..26715ef 100644
--- a/tests/DistortionTests.cpp
+++ b/tests/DistortionTests.cpp
@@ -10,117 +10,114 @@ INSTANTIATE_TEST_CASE_P(JacobianTest,ImageCoordParameterizedTest,
::testing::Values(csm::ImageCoord(2.5,2.5),csm::ImageCoord(7.5,7.5)));
TEST_P(ImageCoordParameterizedTest, JacobianTest) {
- std::vector<double> odtX = {1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0};
- std::vector<double> odtY = {0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0};
+ std::vector<double> transverseDistortionCoeffs = {1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0,
+ 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0};
double Jxx,Jxy,Jyx,Jyy;
csm::ImageCoord imagePt = GetParam();
- std::vector<std::vector<double>> jacobian;
- jacobian = distortionJacobian(imagePt.samp, imagePt.line, odtX, odtY);
+ double jacobian[4];
+ distortionJacobian(imagePt.samp, imagePt.line, jacobian, transverseDistortionCoeffs);
// Jxx * Jyy - Jxy * Jyx
- double determinant = fabs(jacobian[0][0] * jacobian[1][1] - jacobian[0][1] * jacobian[1][0]);
+ double determinant = fabs(jacobian[0] * jacobian[3] - jacobian[1] * jacobian[2]);
EXPECT_GT(determinant,1e-3);
}
TEST(Transverse, Jacobian1) {
csm::ImageCoord imagePt(7.5, 7.5);
- std::vector<double> odtX = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0};
- std::vector<double> odtY = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0};
+ std::vector<double> transverseDistortionCoeffs = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,
+ 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0};
- std::vector<std::vector<double>> jacobian;
- jacobian = distortionJacobian(imagePt.samp, imagePt.line, odtX, odtY);
+ double jacobian[4];
+ distortionJacobian(imagePt.samp, imagePt.line, jacobian, transverseDistortionCoeffs);
- EXPECT_NEAR(jacobian[0][0],56.25,1e-8 );
- EXPECT_NEAR(jacobian[0][1],112.5,1e-8);
- EXPECT_NEAR(jacobian[1][0],56.25,1e-8);
- EXPECT_NEAR(jacobian[1][1],281.25,1e-8);
+ EXPECT_NEAR(jacobian[0],56.25,1e-8 );
+ EXPECT_NEAR(jacobian[1],112.5,1e-8);
+ EXPECT_NEAR(jacobian[2],56.25,1e-8);
+ EXPECT_NEAR(jacobian[3],281.25,1e-8);
}
TEST(Transverse, distortMe_AlternatingOnes) {
csm::ImageCoord imagePt(7.5, 7.5);
- std::vector<double> odtX = {1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0};
- std::vector<double> odtY = {0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0};
+ std::vector<double> transverseDistortionCoeffs = {1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,
+ 0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0};
- double dx,dy;
- std::tuple<double, double> distortedPoint;
- distortedPoint = distortionFunction(imagePt.samp, imagePt.line, odtX, odtY);
+ double dx, dy;
+ computeTransverseDistortion(imagePt.samp, imagePt.line, dx, dy, transverseDistortionCoeffs);
- EXPECT_NEAR(std::get<0>(distortedPoint),908.5,1e-8 );
- EXPECT_NEAR(std::get<1>(distortedPoint),963.75,1e-8);
+ EXPECT_NEAR(dx,908.5,1e-8 );
+ EXPECT_NEAR(dy,963.75,1e-8);
}
TEST(Transverse, distortMe_AllCoefficientsOne) {
csm::ImageCoord imagePt(7.5, 7.5);
- std::vector<double> odtX = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
- std::vector<double> odtY = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
+ std::vector<double> transverseDistortionCoeffs = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
+ 1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
- double dx,dy;
- std::tuple<double, double> distortedPoint;
- distortedPoint = distortionFunction(imagePt.samp, imagePt.line, odtX, odtY);
+ double dx, dy;
+ computeTransverseDistortion(imagePt.samp, imagePt.line, dx, dy, transverseDistortionCoeffs);
- EXPECT_NEAR(std::get<0>(distortedPoint),1872.25,1e-8 );
- EXPECT_NEAR(std::get<1>(distortedPoint),1872.25,1e-8);
+ EXPECT_NEAR(dx,1872.25,1e-8 );
+ EXPECT_NEAR(dy,1872.25,1e-8);
}
TEST(transverse, removeDistortion1) {
csm::ImageCoord imagePt(7.5, 7.5);
+ double ux, uy;
- std::vector<double> odtX = {0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
- std::vector<double> odtY = {0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
+ std::vector<double> transverseDistortionCoeffs = {0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
+ 0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
- std::tuple<double, double> undistortedPoint;
- undistortedPoint = removeDistortion(imagePt.samp, imagePt.line, odtX, odtY);
+ removeDistortion(imagePt.samp, imagePt.line, ux, uy, transverseDistortionCoeffs, DistortionType::TRANSVERSE);
- EXPECT_NEAR(imagePt.samp,7.5,1e-8 );
- EXPECT_NEAR(imagePt.line,7.5,1e-8);
- EXPECT_NEAR(imagePt.line,std::get<0>(undistortedPoint),1e-8);
- EXPECT_NEAR(imagePt.samp,std::get<1>(undistortedPoint),1e-8);
+ EXPECT_NEAR(imagePt.samp, 7.5, 1e-8 );
+ EXPECT_NEAR(imagePt.line, 7.5, 1e-8);
+ EXPECT_NEAR(imagePt.line, ux, 1e-8);
+ EXPECT_NEAR(imagePt.samp, uy, 1e-8);
}
TEST(transverse, removeDistortion_AllCoefficientsOne) {
csm::ImageCoord imagePt(1872.25, 1872.25);
+ double ux, uy;
- std::vector<double> odtX = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
- std::vector<double> odtY = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
+ std::vector<double> transverseDistortionCoeffs = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
+ 1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
- std::tuple<double, double> undistortedPoint;
- undistortedPoint = removeDistortion(imagePt.samp, imagePt.line, odtX, odtY);
+ removeDistortion(imagePt.samp, imagePt.line, ux, uy, transverseDistortionCoeffs, DistortionType::TRANSVERSE);
// The Jacobian is singular, so the setFocalPlane should break out of it's iteration and
// returns the same distorted coordinates that were passed in.
- EXPECT_NEAR(std::get<0>(undistortedPoint),imagePt.samp,1e-8 );
- EXPECT_NEAR(std::get<1>(undistortedPoint),imagePt.line,1e-8);
+ EXPECT_NEAR(ux, imagePt.samp,1e-8 );
+ EXPECT_NEAR(uy, imagePt.line,1e-8);
}
TEST(transverse, removeDistortion_AlternatingOnes) {
csm::ImageCoord imagePt(963.75, 908.5);
+ double ux, uy;
- std::vector<double> odtX = {1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0};
- std::vector<double> odtY = {0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0};
+ std::vector<double> transverseDistortionCoeffs = {1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,
+ 0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0,0.0,1.0};
- std::tuple<double, double> undistortedPoint;
- undistortedPoint = removeDistortion(imagePt.samp, imagePt.line, odtX, odtY);
+ removeDistortion(imagePt.samp, imagePt.line, ux, uy, transverseDistortionCoeffs, DistortionType::TRANSVERSE);
- EXPECT_NEAR(std::get<0>(undistortedPoint),7.5,1e-8 );
- EXPECT_NEAR(std::get<1>(undistortedPoint),7.5,1e-8);
+ EXPECT_NEAR(ux, 7.5, 1e-8 );
+ EXPECT_NEAR(uy, 7.5, 1e-8);
}
TEST(Radial, testRemoveDistortion) {
csm::ImageCoord imagePt(0.0, 4.0);
- double dx, dy;
- double coeffs[3] = {0, 0, 0};
- std::tuple<double, double> undistortedPoint;
+ double ux, uy;
+ std::vector<double> coeffs = {0, 0, 0};
- undistortedPoint = removeDistortion(imagePt.samp, imagePt.line, coeffs);
+ removeDistortion(imagePt.samp, imagePt.line, ux, uy, coeffs, DistortionType::RADIAL);
- EXPECT_NEAR(std::get<0>(undistortedPoint),4,1e-8);
- EXPECT_NEAR(std::get<1>(undistortedPoint),0,1e-8);
+ EXPECT_NEAR(ux, 4, 1e-8);
+ EXPECT_NEAR(uy, 0, 1e-8);
}
// If coeffs are 0 then this will have the same result as removeDistortion
@@ -130,13 +127,13 @@ TEST(Radial, testInverseDistortion){
double dx, dy;
double desiredPrecision = 0.01;
- double coeffs[3] = {0, 0, 0};
- std::tuple<double, double> undistortedPoint;
+ std::vector<double> coeffs = {0, 0, 0};
- undistortedPoint = invertDistortion(imagePt.samp, imagePt.line, coeffs, desiredPrecision);
+ applyDistortion(imagePt.samp, imagePt.line, dx, dy, coeffs,
+ DistortionType::RADIAL, desiredPrecision);
- EXPECT_NEAR(std::get<0>(undistortedPoint),4,1e-8);
- EXPECT_NEAR(std::get<1>(undistortedPoint),0,1e-8);
+ EXPECT_NEAR(dx,4,1e-8);
+ EXPECT_NEAR(dy,0,1e-8);
}
TEST(Radial, testInverseOnesCoeffs){
@@ -144,11 +141,11 @@ TEST(Radial, testInverseOnesCoeffs){
double dx, dy;
double desiredPrecision = 0.01;
- double coeffs[3] = {1, 1, 1};
- std::tuple<double, double> undistortedPoint;
+ std::vector<double> coeffs = {1, 1, 1};
- undistortedPoint = invertDistortion(imagePt.samp, imagePt.line, coeffs, desiredPrecision);
+ applyDistortion(imagePt.samp, imagePt.line, dx, dy, coeffs,
+ DistortionType::RADIAL, desiredPrecision);
- EXPECT_NEAR(std::get<0>(undistortedPoint),4,1e-8);
- EXPECT_NEAR(std::get<1>(undistortedPoint),0,1e-8);
+ EXPECT_NEAR(dx,4,1e-8);
+ EXPECT_NEAR(dy,0,1e-8);
}
diff --git a/tests/ISDParsingTests.cpp b/tests/ISDParsingTests.cpp
index c2976e3..3c2c062 100644
--- a/tests/ISDParsingTests.cpp
+++ b/tests/ISDParsingTests.cpp
@@ -231,28 +231,18 @@ TEST(ISDParsing, SemiMinor) {
EXPECT_EQ(1000, getSemiMinorRadius(isd));
}
-TEST(ISDParsing, TransverseX) {
+TEST(ISDParsing, TransverseDistortion) {
json isd = {
{"optical_distortion", {
{"transverse", {
+ {"y", {-11, 21, 24}},
{"x", {-1, 2, 4}}}}
}
}
};
- std::vector<double> coefficients = {-1, 2, 4, 0, 0, 0, 0, 0, 0, 0};
- EXPECT_EQ(coefficients, getTransverseDistortionX(isd));
-}
-
-TEST(ISDParsing, TransverseY) {
- json isd = {
- {"optical_distortion", {
- {"transverse", {
- {"y", {-11, 21, 24, 16, 20}}}}
- }
- }
- };
- std::vector<double> coefficients = {-11, 21, 24, 16, 20, 0, 0, 0, 0, 0};
- EXPECT_EQ(coefficients, getTransverseDistortionY(isd));
+ std::vector<double> coefficients = {-1, 2, 4, 0, 0, 0, 0, 0, 0, 0,
+ -11, 21, 24, 0, 0, 0, 0, 0, 0, 0};
+ EXPECT_EQ(coefficients, getDistortionCoeffs(isd));
}
TEST(ISDParsing, Radial) {
@@ -264,7 +254,7 @@ TEST(ISDParsing, Radial) {
}
};
std::vector<double> coefficients = {0, 1, 2};
- EXPECT_EQ(coefficients, getRadialDistortion(isd));
+ EXPECT_EQ(coefficients, getDistortionCoeffs(isd));
}
TEST(ISDParsing, SunPosition) {
diff --git a/tests/PluginTests.cpp b/tests/PluginTests.cpp
index 795201a..f034b41 100644
--- a/tests/PluginTests.cpp
+++ b/tests/PluginTests.cpp
@@ -118,19 +118,17 @@ TEST_F(FrameIsdTest, ConstructValidCamera) {
TEST_F(FrameIsdTest, ConstructInValidCamera) {
UsgsAstroPlugin testPlugin;
- isd.setFilename("data/constVelocityLineScan.img");
+ isd.setFilename("data/empty.img");
csm::Model *cameraModel = NULL;
- csm::WarningList *warnings = new csm::WarningList;
try {
testPlugin.constructModelFromISD(
isd,
"USGS_ASTRO_FRAME_SENSOR_MODEL",
- warnings);
+ nullptr);
FAIL() << "Expected an error";
}
catch(csm::Error &e) {
EXPECT_EQ(e.getError(), csm::Error::SENSOR_MODEL_NOT_CONSTRUCTIBLE);
- EXPECT_FALSE(warnings->empty());
}
catch(...) {
FAIL() << "Expected csm SENSOR_MODEL_NOT_CONSTRUCTIBLE error";
@@ -138,9 +136,6 @@ TEST_F(FrameIsdTest, ConstructInValidCamera) {
if (cameraModel) {
delete cameraModel;
}
- if (warnings) {
- delete warnings;
- }
}
TEST_F(ConstVelLineScanIsdTest, Constructible) {
@@ -184,13 +179,13 @@ TEST_F(ConstVelLineScanIsdTest, ConstructValidCamera) {
TEST_F(ConstVelLineScanIsdTest, ConstructInValidCamera) {
UsgsAstroPlugin testPlugin;
- isd.setFilename("data/simpleFramerISD.img");
+ isd.setFilename("data/empty.img");
csm::Model *cameraModel = NULL;
try {
testPlugin.constructModelFromISD(
isd,
"USGS_ASTRO_LINE_SCANNER_SENSOR_MODEL",
- NULL);
+ nullptr);
FAIL() << "Expected an error";
}
diff --git a/tests/data/empty.json b/tests/data/empty.json
new file mode 100644
index 0000000..0967ef4
--- /dev/null
+++ b/tests/data/empty.json
@@ -0,0 +1 @@
+{}
--
GitLab