Isis 3 Programmer Reference
KaguyaMiCameraDistortionMap.cpp
Go to the documentation of this file.
1 
20 #include <cmath>
21 #include <iostream>
22 
23 #include "IString.h"
24 #include "Constants.h"
25 #include "FunctionTools.h"
27 
28 using namespace std;
29 using namespace Isis;
30 
31 namespace Isis {
44  KaguyaMiCameraDistortionMap::KaguyaMiCameraDistortionMap(Camera *parent) : CameraDistortionMap(parent, 1) {
45  }
46 
50  void KaguyaMiCameraDistortionMap::SetDistortion(const int naifIkCode) {
51  //determine if this is the VIS or the NIR sensor, by loooking at the pixel pitch
52  if (p_camera->PixelPitch() == 0.013) m_numDistCoef = 3; //VIS camera has 3 distortion coefs
53  else m_numDistCoef = 4; //NIR camera has 4 distortion coefs
54 
55  //read the distortion coefs from the NAIF Kernels
56  QString naifXKey = "INS" + toString(naifIkCode) + "_DISTORTION_COEF_X";
57  QString naifYKey = "INS" + toString(naifIkCode) + "_DISTORTION_COEF_Y";
58  for (int i=0; i < m_numDistCoef; i++) {
59  m_distCoefX[i] = p_camera->getDouble(naifXKey,i);
60  m_distCoefY[i] = p_camera->getDouble(naifYKey,i);
61  }
62 
63  //now read the boresights, or what I would typicall call the principal point offsets
64  naifXKey = "INS" + toString(naifIkCode) + "_BORESIGHT";
65  m_boreX = p_camera->getDouble(naifXKey, 0);
66  m_boreY = p_camera->getDouble(naifXKey, 1);
67  }
68 
79  bool KaguyaMiCameraDistortionMap::SetFocalPlane(const double dx, const double dy) {
80  p_focalPlaneX = dx;
81  p_focalPlaneY = dy;
82  //cout << "focal plane: " << dx << " " << dy << "\n";
83  //NOTE: the IK/FK kernel does not include the " + dx" as I do below. They also define the
84  // radial distance only in terms of Y. Erroneously (I believe) they use only the
85  // DISTORTION_COEF_X's in their model definition. Finally, they provide different distortion
86  // coefficients for each line of the CCD--despite them going throught the same optical path.
87  // From this I conclude that this distorion model is only valid if x is very near zero. Which
88  // is exactly the situation we are shooting for when modeling a line scanner (x is the along
89  // path direction for this sensor). However, we can not just arbitarily zero, or almost zero,
90  // any along path offset calculated by the back projections. Those offsets are exactly the cost
91  // being zeroed in the iterative LineScanCameraGroundMap routines to find the time that a point
92  // on the ground was imaged. Therefore it must be maintained--with the knowledge that the
93  // small adjustmenst being provided by the distortion model are only relevant as the offsets (x)
94  // approach zero.
95  if (m_numDistCoef == 3) { //VIS camera
97  m_boreX + m_distCoefX[0] + dy*(m_distCoefX[1] + dy*m_distCoefX[2]) + dx;
99  m_boreY + m_distCoefY[0] + dy*(m_distCoefY[1] + dy*m_distCoefY[2]) + dy;
100  }
101  else { //NIR camera
102  p_undistortedFocalPlaneX = m_boreX +
103  m_distCoefX[0] + dy*(m_distCoefX[1] + dy*(m_distCoefX[2] + dy*m_distCoefX[3])) + dx;
104  p_undistortedFocalPlaneY = m_boreY +
105  m_distCoefY[0] + dy*(m_distCoefY[1] + dy*(m_distCoefY[2] + dy*m_distCoefY[3])) + dy;
106  }
107  //cout << "undistorted: " << p_undistortedFocalPlaneX << " " << p_undistortedFocalPlaneY << "\n";
108  return true;
109  }
110 
121  bool KaguyaMiCameraDistortionMap::SetUndistortedFocalPlane(const double ux, const double uy) {\
122  // image coordinates prior to introducing distortion
125 
126  //std::cout << "Distortionless : " << p_undistortedFocalPlaneX << " " << p_undistortedFocalPlaneY << "\n";
127 
128  if (m_numDistCoef == 3) { //quadratic distortion model
129  //use the quadratic equation to find the distorted Y value
130  double A,B,C; //coefficients of the quadric equation
131 
132  A = m_distCoefY[2];
133  B = 1.0 + m_distCoefY[1];
134  C = m_distCoefY[0] +m_boreY - uy;
135 
137 
138  if (roots.size() == 0) return false;
139  else if (roots.size() == 1)
140  p_focalPlaneY = roots[0];
141  else //two roots to choose between--choose the one closest to uy
142  p_focalPlaneY = fabs(uy - roots[0]) < fabs(uy - roots[1]) ? roots[0] : roots[1];
143 
144  //now that we know the distortedY we can directly calculate the X
145  p_focalPlaneX = ux -
146  (m_boreX + m_distCoefX[0] + p_focalPlaneY*(m_distCoefX[1] + p_focalPlaneY*m_distCoefX[2]));
147 
148  return true;
149  }
150  else { //cubic distortion model
151  //use the cubic equation to find the distorted Y value
152  double a,b,c; //coefficients of the cubic (coef x^3 ==1)
153 
154  a = m_distCoefY[2] / m_distCoefY[3];
155  b = (1.0 + m_distCoefY[1]) / m_distCoefY[3];
156  c = (m_distCoefY[0] + m_boreY - uy) / m_distCoefY[3];
157 
158  QList<double> roots = FunctionTools::realCubicRoots(1.0,a,b,c);
159 
160  if (roots.size() == 1) { //one root
161  p_focalPlaneY = roots[0];
162  }
163  else { //pick the root closest to the orginal value
164  QList<double> delta;
165 
166  //store all the distance from undistored to focal plane
167  for (int i=0;i<roots.size();i++)
168  delta << fabs(roots[i]-uy);
169 
170  if ( roots.size() == 3) //three roots to choose among
171  p_focalPlaneY = delta[0] < delta[1] ?
172  (delta[0] < delta[2] ? roots[0]:roots[2]) :
173  (delta[1] < delta[2] ? roots[1]:roots[2]) ;
174  else //two roots to choose between
175  p_focalPlaneY = delta[0] < delta[1] ? roots[0]:roots[1] ;
176  }
177 
178  //now that we know the distortedY we can directly calculate the X
179  p_focalPlaneX = ux - (m_boreX +m_distCoefX[0] +
180  p_focalPlaneY*(m_distCoefX[1] +
181  p_focalPlaneY*(m_distCoefX[2] +
182  p_focalPlaneY* m_distCoefX[3])));
183  return true;
184  }
185  }
186 }
double p_focalPlaneX
Distorted focal plane x.
static QList< double > realCubicRoots(double coeffCubicTerm, double coeffQuadTerm, double coeffLinearTerm, double coeffConstTerm)
Find the real roots of a cubic (1, 2, or 3) (see numerical recipies 3rd edtion page 227) Form: coeffC...
virtual bool SetUndistortedFocalPlane(const double ux, const double uy)
Compute distorted focal plane x/y.
Namespace for the standard library.
QString toString(bool boolToConvert)
Global function to convert a boolean to a string.
Definition: IString.cpp:226
double p_undistortedFocalPlaneX
Undistorted focal plane x.
Camera * p_camera
The camera to distort/undistort.
Unless noted otherwise, the portions of Isis written by the USGS are public domain.
static QList< double > realQuadraticRoots(double coeffQuadTerm, double coeffLinearTerm, double coeffConstTerm)
The correct way to find the real roots of a quadratic (0, 1, or 2) (according to numerical recipies 3...
Definition: FunctionTools.h:87
Distort/undistort focal plane coordinates.
double PixelPitch() const
Returns the pixel pitch.
Definition: Camera.cpp:2754
virtual bool SetFocalPlane(const double dx, const double dy)
Compute undistorted focal plane x/y.
double p_focalPlaneY
Distorted focal plane y.
double p_undistortedFocalPlaneY
Undistorted focal plane y.
Namespace for ISIS/Bullet specific routines.
Definition: Apollo.h:31
SpiceDouble getDouble(const QString &key, int index=0)
This returns a value from the NAIF text pool.
Definition: Spice.cpp:963