Isis 3 Programmer Reference
NewHorizonsMvicTdiCameraDistortionMap.cpp
1
7/* SPDX-License-Identifier: CC0-1.0 */
8
9#include <cmath>
10#include <iostream>
11
12#include <boost/math/special_functions/legendre.hpp>
13
14#include "Camera.h"
15#include "CameraFocalPlaneMap.h"
16#include "Constants.h"
17#include "FunctionTools.h"
18#include "IString.h"
19#include "LineScanCameraDetectorMap.h"
20#include "NewHorizonsMvicTdiCameraDistortionMap.h"
21
22#include <QDebug>
23
24
25using namespace boost::math;
26using namespace std;
27using namespace Isis;
28
29namespace Isis {
48 Camera *parent,
49 vector<double> xDistortionCoeffs,
50 vector<double> yDistortionCoeffs,
51 vector<double> residualColDistCoeffs,
52 vector<double> residualRowDistCoeffs) :
53 CameraDistortionMap(parent, 1.0) {
54
55 m_xDistortionCoeffs = xDistortionCoeffs;
56 m_yDistortionCoeffs = yDistortionCoeffs;
57
58 m_residualColDistCoeffs = residualColDistCoeffs;
59 m_residualRowDistCoeffs = residualRowDistCoeffs;
60
61 // half of detector in x is 32.5 mm
63 }
64
65
70
71
88 // TODO: Don't really like this - we always return true, even in event of failures
89 bool NewHorizonsMvicTdiCameraDistortionMap::SetFocalPlane(const double dx, const double dy) {
90
91 p_focalPlaneX = dx;
92 p_focalPlaneY = dy;
93
94 // initialize undistorted focal plane coordinates to the distorted
95 // coordinate values
98
99 // if x lies outside of the detector, do NOT apply distortion
100 // set undistorted focal plane values to be identical to raw values
101 if ((fabs(dx) > m_focalPlaneHalf_x)) {
102 return true;
103 }
104
105 // scale x to lie in the range -1.0 to +1.0
106 // this is required for Legendre Polynomials, man
107 double xscaled = -dx / m_focalPlaneHalf_x;
108
109 // compute distortion corrections in x using Legendre Polynomials
110 // these corrections are also in the -1.0 to +1.0 range
111 double deltax1 = 0.0;
112 if (!computeDistortionCorrections(xscaled, 0.0, deltax1)) {
113 return true;
114 }
115
116 // now compute residual distortion corrections (per Jason Cook)
117 // TODO: implementation not complete
118 double deltax2 = 0.0;
119 double deltay2 = 0.0;
120 computeResidualDistortionCorrections(dx, deltax2, deltay2);
121
122 // apply the corrections
123 xscaled += deltax1;
124
125 // scale back from range of '-1.0 to +1.0' to the detector, '-32.5 to 32.5 mm'
126// p_undistortedFocalPlaneX = -xscaled * m_focalPlaneHalf_x;
127// p_undistortedFocalPlaneY = p_focalPlaneY;
128 p_undistortedFocalPlaneX = -xscaled * m_focalPlaneHalf_x + deltax2;
130
131 return true;
132 }
133
134
149
150 // image coordinates prior to introducing distortion
153
154 double xt = ux;
155 double yt = uy;
156
157// p_focalPlaneX = p_undistortedFocalPlaneX;
158// p_focalPlaneY = p_undistortedFocalPlaneY;
159// return true;
160
161 // scale undistorted coordinates to range of -1.0 to +1.0
162 double uxScaled = -ux/m_focalPlaneHalf_x;
163 double xtScaled;
164 double xprevious, yprevious;
165 double scaledDeltax;
166 double deltax2 = 0.0;
167 double deltay2 = 0.0;
168
169 xprevious = 1000000.0;
170 yprevious = 1000000.0;
171
172 double tolerance = 0.000001;
173
174 bool bConverged = false;
175
176 // iterating to introduce distortion...
177 // we stop when the difference between distorted coordinates
178 // in successive iterations is at or below the given tolerance
179 for ( int i = 0; i < 50; i++ ) {
180
181 xtScaled = -xt / m_focalPlaneHalf_x;
182
183 // Changed this failure mode from a throw to return false. These functions should never
184 // throw, the callers are not catching, so bad things happen. The callers should be checking
185 // the return status
186 if (fabs(xtScaled) > 1.0) {
187 return bConverged;
188 }
189
190 // compute scaled distortion in x (scaled to -1.0 - +1.0) using Legendre Polynomials
191 computeDistortionCorrections(xtScaled, 0.0, scaledDeltax);
192
193 // compute residual distortion in unscaled focal plane coordinates
194 computeResidualDistortionCorrections(xt, deltax2, deltay2);
195
196 // update unscaled coordinates
197 xt = -(uxScaled-scaledDeltax) * m_focalPlaneHalf_x - deltax2;
198 yt = uy - deltay2;
199
200 // check for convergenceyScaledDistortion
201 if ((fabs(xt - xprevious) <= tolerance) && (fabs(yt - yprevious) <= tolerance)) {
202 bConverged = true;
203 break;
204 }
205
206 xprevious = xt;
207 yprevious = yt;
208 }
209
210 if (bConverged) {
211 // scale x coord inate back to detector (-32.5 to +32.5)
212 // xtScaled *= -m_focalPlaneHalf_x;
213
214 // set distorted coordinates
215 p_focalPlaneX = xt;
216 p_focalPlaneY = yt;
217 }
218
219 return bConverged;
220 }
221
222
238 const double xscaled,
239 const double yscaled,
240 double &deltax) {
241
242 double lpx0, lpx1, lpx2, lpx3, lpx4, lpx5;
243 double lpy0, lpy1, lpy2, lpy3, lpy4, lpy5;
244
245 // Legendre polynomials
246 try {
247 lpx0 = legendre_p(0,xscaled);
248 lpx1 = legendre_p(1,xscaled);
249 lpx2 = legendre_p(2,xscaled);
250 lpx3 = legendre_p(3,xscaled);
251 lpx4 = legendre_p(4,xscaled);
252 lpx5 = legendre_p(5,xscaled);
253 lpy0 = legendre_p(0,yscaled);
254 lpy1 = legendre_p(1,yscaled);
255 lpy2 = legendre_p(2,yscaled);
256 lpy3 = legendre_p(3,yscaled);
257 lpy4 = legendre_p(4,yscaled);
258 lpy5 = legendre_p(5,yscaled);
259 }
260 catch (const std::exception& e) {
261 return false;
262 }
263
264 deltax =
265 m_xDistortionCoeffs[0] * lpx0 * lpy1 +
266 m_xDistortionCoeffs[1] * lpx1 * lpy0 +
267 m_xDistortionCoeffs[2] * lpx0 * lpy2 +
268 m_xDistortionCoeffs[3] * lpx1 * lpy1 +
269 m_xDistortionCoeffs[4] * lpx2 * lpy0 +
270 m_xDistortionCoeffs[5] * lpx0 * lpy3 +
271 m_xDistortionCoeffs[6] * lpx1 * lpy2 +
272 m_xDistortionCoeffs[7] * lpx2 * lpy1 +
273 m_xDistortionCoeffs[8] * lpx3 * lpy0 +
274 m_xDistortionCoeffs[9] * lpx0 * lpy4 +
275 m_xDistortionCoeffs[10] * lpx1 * lpy3 +
276 m_xDistortionCoeffs[11] * lpx2 * lpy2 +
277 m_xDistortionCoeffs[12] * lpx3 * lpy1 +
278 m_xDistortionCoeffs[13] * lpx4 * lpy0 +
279 m_xDistortionCoeffs[14] * lpx0 * lpy5 +
280 m_xDistortionCoeffs[15] * lpx1 * lpy4 +
281 m_xDistortionCoeffs[16] * lpx2 * lpy3 +
282 m_xDistortionCoeffs[17] * lpx3 * lpy2 +
283 m_xDistortionCoeffs[18] * lpx4 * lpy1 +
284 m_xDistortionCoeffs[19] * lpx5 * lpy0;
285
286 return true;
287 }
288
289
301 double &residualDeltax,
302 double &residualDeltay) {
303
304 double s = 2500.5 - dx/0.013;
305
306 residualDeltax = s * m_residualColDistCoeffs[1] +
307 pow(s,2) * m_residualColDistCoeffs[2] +
308 pow(s,3) * m_residualColDistCoeffs[3] +
309 pow(s,4) * m_residualColDistCoeffs[4] +
310 pow(s,5) * m_residualColDistCoeffs[5];
311
312 residualDeltay = s * m_residualRowDistCoeffs[1] +
313 pow(s,2) * m_residualRowDistCoeffs[2] +
314 pow(s,3) * m_residualRowDistCoeffs[3] +
315 pow(s,4) * m_residualRowDistCoeffs[4] +
316 pow(s,5) * m_residualRowDistCoeffs[5];
317
318 // convert to mm (correction for negative x to the left)
319 residualDeltax *= -p_camera->PixelPitch();
320 residualDeltay *= p_camera->PixelPitch();
321 }
322
323
325// * Testing method to output corrections in x and y at pixel centers for entire focal plane.
326// * Output in csv format for viewing/plotting in Excel.
327// */
328//bool NewHorizonsMvicTdiCameraDistortionMap::outputResidualDeltas() {
329//
330// QString ofname("mvic_tdi_residual_deltas.csv");
331// std::ofstream fp_out(ofname.toLatin1().data(), std::ios::out);
332// if (!fp_out)
333// return false;
334//
335// char buf[1056];
336//
337// double residualColDelta;
338//
339// for (int s=0; s <= 5000; s++) { // loop in sample direction
340// residualColDelta = m_residualColDistCoeffs[0] +
341// s * m_residualColDistCoeffs[1] +
342// pow(double(s),2) * m_residualColDistCoeffs[2] +
343// pow(double(s),3) * m_residualColDistCoeffs[3] +
344// pow(double(s),4) * m_residualColDistCoeffs[4] +
345// pow(double(s),5) * m_residualColDistCoeffs[5];
346// sprintf(buf, "%d,%lf\n", s, residualColDelta);
347// fp_out << buf;
348// }
349//
350// fp_out.close();
351//
352// return true;
353//}
354
355}
Distort/undistort focal plane coordinates.
double p_focalPlaneX
Distorted focal plane x.
double p_undistortedFocalPlaneX
Undistorted focal plane x.
double p_undistortedFocalPlaneY
Undistorted focal plane y.
Camera * p_camera
The camera to distort/undistort.
double p_focalPlaneY
Distorted focal plane y.
double PixelPitch() const
Returns the pixel pitch.
Definition Camera.cpp:2772
int Samples() const
Returns the number of samples in the image.
Definition Camera.cpp:2806
NewHorizonsMvicTdiCameraDistortionMap(Camera *parent, vector< double > xDistortionCoeffs, vector< double > yDistortionCoeffs, vector< double > residualColDistCoeffs, vector< double > residualRowDistCoeffs)
Camera distortion map constructor.
void computeResidualDistortionCorrections(const double dx, double &residualDeltax, double &residualDeltay)
Compute residual distortion corrections in row and column direction TODO: Implementati plete.
double m_focalPlaneHalf_x
half of focal plane x and y dimensions in mm
virtual bool SetFocalPlane(const double dx, const double dy)
Compute undistorted focal plane x/y.
std::vector< double > m_yDistortionCoeffs
by Keith Harrison (Interface Control Document section 10.3.1.2)
virtual bool SetUndistortedFocalPlane(const double ux, const double uy)
Compute distorted focal plane x/y.
std::vector< double > m_xDistortionCoeffs
distortion coefficients in x and y as determined
vector< double > m_residualColDistCoeffs
residual distortion coefficients as determined
bool computeDistortionCorrections(const double xscaled, const double yscaled, double &deltax)
Compute distortion corrections in x and y direction.
vector< double > m_residualRowDistCoeffs
by Jason Cook, SWRI (MVIC Distortion)
This is free and unencumbered software released into the public domain.
Definition Apollo.h:16
Namespace for the standard library.