Isis 3 Programmer Reference
NewHorizonsMvicFrameCameraDistortionMap.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 "Constants.h"
16#include "FunctionTools.h"
17#include "IString.h"
18#include "NewHorizonsMvicFrameCameraDistortionMap.h"
19
20#include "CameraFocalPlaneMap.h"
21
22
23#include <QDebug>
24
25
26using namespace boost::math;
27using namespace std;
28using namespace Isis;
29
30namespace Isis {
44 vector<double> xDistortionCoeffs,
45 vector<double> yDistortionCoeffs) :
46 CameraDistortionMap(parent, 1.0) {
47
48 m_xDistortionCoeffs = xDistortionCoeffs;
49 m_yDistortionCoeffs = yDistortionCoeffs;
50
51 double pixelPitch = p_camera->PixelPitch();
52
53 m_focalPlaneHalf_x = 0.5 * p_camera->Samples() * pixelPitch; // 32.5 mm
54 m_focalPlaneHalf_y = 0.5 * p_camera->Lines() * pixelPitch; // 0.832 mm
55 }
56
57
62
63
64// /**
65// * Testing method to output corrections in x and y at pixel centers for entire focal plane.
66// * Output in csv format for viewing/plotting in Excel.
67// */
68//bool NewHorizonsMvicFrameCameraDistortionMap::outputDeltas() {
69//
70// QString ofname("mvic_frame_deltas.csv");
71// std::ofstream fp_out(ofname.toLatin1().data(), std::ios::out);
72// if (!fp_out)
73// return false;
74//
75// char buf[1056];
76//
77// double deltax, deltay;
78//
79// for (double line = 0.5; line <= 128.5; line += 1.0) { // loop in y direction
80// for (double sample=0.5; sample <= 5000.5; sample += 1.0) { // loop in x direction
81//
82// p_camera->FocalPlaneMap()->SetDetector(sample,line);
83//
84// double fplanex = p_camera->FocalPlaneMap()->FocalPlaneX();
85// double fplaney = p_camera->FocalPlaneMap()->FocalPlaneY();
86//
87// SetFocalPlane(fplanex,fplaney);
88//
89// deltax = fplanex - p_undistortedFocalPlaneX;
90// deltay = fplaney - p_undistortedFocalPlaneY;
91//
92// sprintf(buf, "%lf,%lf,%lf,%lf\n", sample, deltax/0.013, line, deltay/0.013);
93//
94// fp_out << buf;
95// }
96// }
97//
98// fp_out.close();
99//
100// return true;
101//}
102
103
114 bool NewHorizonsMvicFrameCameraDistortionMap::SetFocalPlane(const double dx, const double dy) {
115
116 p_focalPlaneX = dx;
117 p_focalPlaneY = dy;
118
119 // in case of failures, initialize undistorted focal plane coordinates to the distorted
120 // coordinate values
123
124 // if x and/or y lie outside of the detector, do NOT apply distortion
125 // set undistorted focal plane values to be identical to raw values
126 if ((fabs(dx) > m_focalPlaneHalf_x) || (fabs(dy) > m_focalPlaneHalf_y)) {
127 return true;
128 }
129
130 // shift from ISIS MVIC FT image coordinate system with +x to the left and +y down to
131 // the desired system of +x to the right and +y up
132 // e.g. negate x and y
133
134 // scale x and y to lie in the range -1.0 to +1.0
135 // this is requirement for Legendre Polynomials, man
136 double xscaled = -dx/m_focalPlaneHalf_x;
137 double yscaled = -dy/m_focalPlaneHalf_y;
138
139 // compute distortion corrections in x and y using Legendre Polynomials
140 // these corrections are also in the -1.0 to +1.0 range
141 // if Legendre computations fail, we set undistorted focal plane x and y
142 // to be identical to distorted x and y and return true
143 double deltax, deltay;
144 if (!computeDistortionCorrections(xscaled, yscaled, deltax, deltay)) {
145 return true;
146 }
147
148 // apply the corrections to original x,y
149 xscaled += deltax;
150 yscaled += deltay;
151
152 // scale back from range of '-1.0 to +1.0' to the detector, '-32.5 to 32.5 mm'
154 p_undistortedFocalPlaneY = -yscaled * m_focalPlaneHalf_y;
155
156 return true;
157 }
158
159
160// bool NewHorizonsMvicFrameCameraDistortionMap::SetFocalPlane(const double dx, const double dy) {
161
162// p_focalPlaneX = dx;
163// p_focalPlaneY = dy;
164
165// // if x and/or y lie outside of the detector, do NOT apply distortion
166// // set undistorted focal plane values to be identical to raw values
167// if ((fabs(dx) > m_focalPlaneHalf_x) || (fabs(dy) > m_focalPlaneHalf_y)) {
168// p_undistortedFocalPlaneX = dx;
169// p_undistortedFocalPlaneY = dy;
170
171// return true;
172// }
173
174// // scale x and y to lie in the range -1.0 to +1.0
175// // this is requirement for Legendre Polynomials, man
176// double xscaled = dx/m_focalPlaneHalf_x;
177// double yscaled = dy/m_focalPlaneHalf_y;
178
179// // compute distortion corrections in x and y using Legendre Polynomials
180// // these corrections are also in the -1.0 to +1.0 range
181// double deltax, deltay;
182// computeDistortionCorrections(xscaled, yscaled, deltax, deltay);
183
184// // apply the corrections
185// xscaled += deltax;
186// yscaled += deltay;
187
188// // scale back from range of '-1.0 to +1.0' to the detector, '-32.656 to 32.656 mm'
189// p_undistortedFocalPlaneX = xscaled * m_focalPlaneHalf_x;
190// p_undistortedFocalPlaneY = yscaled * m_focalPlaneHalf_y;
191
192// return true;
193// }
194
195
209
210 // image coordinates prior to introducing distortion
213
214 double xScaledDistortion, yScaledDistortion;
215
216 // scale undistorted coordinates to range of -1.0 to +1.0
217 double xtScaled = -ux/m_focalPlaneHalf_x;
218 double ytScaled = -uy/m_focalPlaneHalf_y;
219
220 double uxScaled = xtScaled;
221 double uyScaled = ytScaled;
222
223 double xScaledPrevious = 1000000.0;
224 double yScaledPrevious = 1000000.0;
225
226 double tolerance = 0.000001;
227
228 bool bConverged = false;
229
230 // iterating to introduce distortion...
231 // we stop when the difference between distorted coordinates
232 // in successive iterations is at or below the given tolerance
233 for( int i = 0; i < 50; i++ ) {
234
235 // compute distortion in x and y (scaled to -1.0 - +1.0) using Legendre Polynomials
236 computeDistortionCorrections(xtScaled, ytScaled, xScaledDistortion, yScaledDistortion);
237
238 // update scaled image coordinates
239 xtScaled = uxScaled - xScaledDistortion;
240 ytScaled = uyScaled - yScaledDistortion;
241
242 // check for convergence
243 if((fabs(xtScaled - xScaledPrevious) <= tolerance) && (fabs(ytScaled - yScaledPrevious) <= tolerance)) {
244 bConverged = true;
245 break;
246 }
247
248 xScaledPrevious = xtScaled;
249 yScaledPrevious = ytScaled;
250 }
251
252 if (bConverged) {
253 // scale coordinates back to detector (-32.5 to +32.5)
254 xtScaled *= -m_focalPlaneHalf_x;
255 ytScaled *= -m_focalPlaneHalf_y;
256
257 // set distorted coordinates
258 p_focalPlaneX = xtScaled;
259 p_focalPlaneY = ytScaled;
260 }
261
262 return bConverged;
263 }
264// bool NewHorizonsMvicFrameCameraDistortionMap::SetUndistortedFocalPlane(const double ux, const double uy) {
265
266// // image coordinates prior to introducing distortion
267// p_undistortedFocalPlaneX = ux;
268// p_undistortedFocalPlaneY = uy;
269
270// double xScaledDistortion, yScaledDistortion;
271
272// // scale undistorted coordinates to range of -1.0 to +1.0
273// double xtScaled = ux/m_focalPlaneHalf_x;
274// double ytScaled = uy/m_focalPlaneHalf_y;
275
276// double uxScaled = xtScaled;
277// double uyScaled = ytScaled;
278
279// double xScaledPrevious = 1000000.0;
280// double yScaledPrevious = 1000000.0;
281
282// double tolerance = 0.000001;
283
284// bool bConverged = false;
285
286// // iterating to introduce distortion...
287// // we stop when the difference between distorted coordinates
288// // in successive iterations is at or below the given tolerance
289// for( int i = 0; i < 50; i++ ) {
290
291// // compute distortion in x and y (scaled to -1.0 - +1.0) using Legendre Polynomials
292// computeDistortionCorrections(xtScaled, ytScaled, xScaledDistortion, yScaledDistortion);
293
294// // update scaled image coordinates
295// xtScaled = uxScaled - xScaledDistortion;
296// ytScaled = uyScaled - yScaledDistortion;
297
298// // check for convergence
299// if((fabs(xtScaled - xScaledPrevious) <= tolerance) && (fabs(ytScaled - yScaledPrevious) <= tolerance)) {
300// bConverged = true;
301// break;
302// }
303
304// xScaledPrevious = xtScaled;
305// yScaledPrevious = ytScaled;
306// }
307
308// if (bConverged) {
309// // scale coordinates back to detector (-32.656 to +32.656)
310// xtScaled *= m_focalPlaneHalf_x;
311// ytScaled *= m_focalPlaneHalf_y;
312
313// // set distorted coordinates
314// p_focalPlaneX = xtScaled;
315// p_focalPlaneY = ytScaled;
316// }
317
318// return bConverged;
319// }
320
321
337 const double yscaled,
338 double &deltax, double &deltay) {
339
340 double lpx0, lpx1, lpx2, lpx3, lpx4, lpx5;
341 double lpy0, lpy1, lpy2, lpy3, lpy4, lpy5;
342
343 // Legendre polynomials
344 // boost library method legendre_p will generate an exception if xscaled or yscaled do not lie
345 // between -1 to 1 (inclusive). In this event we return false.
346 try {
347 lpx0 = legendre_p(0,xscaled);
348 lpx1 = legendre_p(1,xscaled);
349 lpx2 = legendre_p(2,xscaled);
350 lpx3 = legendre_p(3,xscaled);
351 lpx4 = legendre_p(4,xscaled);
352 lpx5 = legendre_p(5,xscaled);
353 lpy0 = legendre_p(0,yscaled);
354 lpy1 = legendre_p(1,yscaled);
355 lpy2 = legendre_p(2,yscaled);
356 lpy3 = legendre_p(3,yscaled);
357 lpy4 = legendre_p(4,yscaled);
358 lpy5 = legendre_p(5,yscaled);
359 }
360 // TESTING NOTE: Could not find a way to cause this error. If one is found a test should be added
361 catch (const std::exception& e) {
362 return false;
363 }
364
365 deltax =
366 m_xDistortionCoeffs[0] * lpx0 * lpy1 +
367 m_xDistortionCoeffs[1] * lpx1 * lpy0 +
368 m_xDistortionCoeffs[2] * lpx0 * lpy2 +
369 m_xDistortionCoeffs[3] * lpx1 * lpy1 +
370 m_xDistortionCoeffs[4] * lpx2 * lpy0 +
371 m_xDistortionCoeffs[5] * lpx0 * lpy3 +
372 m_xDistortionCoeffs[6] * lpx1 * lpy2 +
373 m_xDistortionCoeffs[7] * lpx2 * lpy1 +
374 m_xDistortionCoeffs[8] * lpx3 * lpy0 +
375 m_xDistortionCoeffs[9] * lpx0 * lpy4 +
376 m_xDistortionCoeffs[10] * lpx1 * lpy3 +
377 m_xDistortionCoeffs[11] * lpx2 * lpy2 +
378 m_xDistortionCoeffs[12] * lpx3 * lpy1 +
379 m_xDistortionCoeffs[13] * lpx4 * lpy0 +
380 m_xDistortionCoeffs[14] * lpx0 * lpy5 +
381 m_xDistortionCoeffs[15] * lpx1 * lpy4 +
382 m_xDistortionCoeffs[16] * lpx2 * lpy3 +
383 m_xDistortionCoeffs[17] * lpx3 * lpy2 +
384 m_xDistortionCoeffs[18] * lpx4 * lpy1 +
385 m_xDistortionCoeffs[19] * lpx5 * lpy0;
386
387 deltay =
388 m_yDistortionCoeffs[0] * lpx0 * lpy1 +
389 m_yDistortionCoeffs[1] * lpx1 * lpy0 +
390 m_yDistortionCoeffs[2] * lpx0 * lpy2 +
391 m_yDistortionCoeffs[3] * lpx1 * lpy1 +
392 m_yDistortionCoeffs[4] * lpx2 * lpy0 +
393 m_yDistortionCoeffs[5] * lpx0 * lpy3 +
394 m_yDistortionCoeffs[6] * lpx1 * lpy2 +
395 m_yDistortionCoeffs[7] * lpx2 * lpy1 +
396 m_yDistortionCoeffs[8] * lpx3 * lpy0 +
397 m_yDistortionCoeffs[9] * lpx0 * lpy4 +
398 m_yDistortionCoeffs[10] * lpx1 * lpy3 +
399 m_yDistortionCoeffs[11] * lpx2 * lpy2 +
400 m_yDistortionCoeffs[12] * lpx3 * lpy1 +
401 m_yDistortionCoeffs[13] * lpx4 * lpy0 +
402 m_yDistortionCoeffs[14] * lpx0 * lpy5 +
403 m_yDistortionCoeffs[15] * lpx1 * lpy4 +
404 m_yDistortionCoeffs[16] * lpx2 * lpy3 +
405 m_yDistortionCoeffs[17] * lpx3 * lpy2 +
406 m_yDistortionCoeffs[18] * lpx4 * lpy1 +
407 m_yDistortionCoeffs[19] * lpx5 * lpy0;
408
409 return true;
410 }
411}
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 Lines() const
Returns the number of lines in the image.
Definition Camera.cpp:2816
int Samples() const
Returns the number of samples in the image.
Definition Camera.cpp:2806
std::vector< double > m_xDistortionCoeffs
distortion coefficients in x and y as determined
NewHorizonsMvicFrameCameraDistortionMap(Camera *parent, vector< double > xDistortionCoeffs, vector< double > yDistortionCoeffs)
Camera distortion map constructor.
virtual bool SetFocalPlane(const double dx, const double dy)
Compute undistorted focal plane x/y.
double m_focalPlaneHalf_x
half of focal plane x and y dimensions in mm
virtual bool SetUndistortedFocalPlane(const double ux, const double uy)
Compute distorted focal plane x/y.
std::vector< double > m_yDistortionCoeffs
by Keith Harrison (Interface Control Document section 10.3.1.2)
bool computeDistortionCorrections(const double xscaled, const double yscaled, double &deltax, double &deltay)
Compute distortion corrections in x and y direction.
This is free and unencumbered software released into the public domain.
Definition Apollo.h:16
Namespace for the standard library.