Isis 3 Programmer Reference
LineScanCameraGroundMap.cpp
1
6/* SPDX-License-Identifier: CC0-1.0 */
7
8#include "LineScanCameraGroundMap.h"
9
10#include <iostream>
11#include <iomanip>
12
13#include <QTime>
14#include <QList>
15#include <QFile>
16#include <QTextStream>
17
18#include "IException.h"
19#include "IString.h"
20#include "Camera.h"
21#include "CameraDistortionMap.h"
22#include "CameraFocalPlaneMap.h"
23#include "Distance.h"
24#include "LineScanCameraDetectorMap.h"
25#include "iTime.h"
26#include "Latitude.h"
27#include "Longitude.h"
28#include "Statistics.h"
29#include "SurfacePoint.h"
30#include "FunctionTools.h"
31
32
33using namespace std;
34using namespace Isis;
35
36bool ptXLessThan(const QList<double> l1, const QList<double> l2);
37
44 public std::function<double(double)> {
45 public:
46
48 m_camera = camera;
49 m_surfacePoint = surPt;
50 }
51
52
54
55
64 double operator()(double et) {
65 double lookC[3] = {0.0, 0.0, 0.0};
66 double ux = 0.0;
67 double uy = 0.0;
68 double dx = 0.0;
69 double dy = 0.0;
70
71 // Verify the time is with the cache bounds
72 double startTime = m_camera->cacheStartTime().Et();
73 double endTime = m_camera->cacheEndTime().Et();
74 if (et < startTime || et > endTime) {
75 IString msg = "Ephemeris time passed to LineOffsetFunctor is not within the image "
76 "cache bounds";
77 throw IException(IException::Programmer, msg, _FILEINFO_);
78 }
79
80 m_camera->Sensor::setTime(et);
81
82 // Set ground
83 if (!m_camera->Sensor::SetGround(m_surfacePoint, false)) {
84 IString msg = "Sensor::SetGround failed for surface point in LineScanCameraGroundMap.cpp"
85 " LineOffsetFunctor";
86 throw IException(IException::Programmer, msg, _FILEINFO_);
87 }
88
89 // Calculate the undistorted focal plane coordinates
90 m_camera->Sensor::LookDirection(lookC);
91 ux = m_camera->FocalLength() * lookC[0] / lookC[2];
92 uy = m_camera->FocalLength() * lookC[1] / lookC[2];
93
94
95 // This was replaced with the code below to get Chandrayaan M3 to work.
96 // SetUndistortedFocalPlane was failing a majority of the time, causing most SetGround calls
97 // to fail. Even when it did succeed, it was producing non-continous return values.
98 // Set the undistorted focal plane coordinates
99 // if (!m_camera->DistortionMap()->SetUndistortedFocalPlane(ux, uy)) {
100 // IString msg = "DistortionMap::SetUndistoredFocalPlane failed for surface point in "
101 // "LineScanCameraGroundMap.cpp LineOffsetFunctor";
102 // throw IException(IException::Programmer, msg, _FILEINFO_);
103 // }
104
105 // Get the natural (distorted focal plane coordinates)
106 // dx = m_camera->DistortionMap()->FocalPlaneX();
107 // dy = m_camera->DistortionMap()->FocalPlaneY();
108 // std::cout << "use dist" << std::endl;
109 //}
110
111
112 // Try to use SetUndistortedFocalPlane, if that does not work use the distorted x,y
113 // under the assumption (bad|good) that extrapolating the distortion
114 // is causing the distorted x,y to be way off the sensor, and thus not very good anyway.
115 if (m_camera->DistortionMap()->SetUndistortedFocalPlane(ux, uy)) {
116 // Get the natural (distorted focal plane coordinates)
117 dx = m_camera->DistortionMap()->FocalPlaneX();
118 dy = m_camera->DistortionMap()->FocalPlaneY();
119 }
120 else {
121 dx = ux;
122 dy = uy;
123 }
124
125 if (!m_camera->FocalPlaneMap()->SetFocalPlane(dx, dy)) {
126 IString msg = "FocalPlaneMap::SetFocalPlane failed for surface point in "
127 "LineScanCameraGroundMap.cpp LineOffsetFunctor";
128 throw IException(IException::Programmer, msg, _FILEINFO_);
129 }
130
131 // Return the offset
132 return (m_camera->FocalPlaneMap()->DetectorLineOffset() -
133 m_camera->FocalPlaneMap()->DetectorLine());
134 }
135
136
137 private:
138 SurfacePoint m_surfacePoint;
139 Camera* m_camera;
140};
141
142
149 public std::function<double(double)> {
150
151 public:
153 m_camera = camera;
154 surfacePoint = surPt;
155 }
156
157
159
160
161 double operator()(double et) {
162 double s[3], p[3];
163
164 //verify the time is with the cache bounds
165 double startTime = m_camera->cacheStartTime().Et();
166 double endTime = m_camera->cacheEndTime().Et();
167
168 if (et < startTime || et > endTime) {
169 IString msg = "Ephemeris time passed to SensorSurfacePointDistanceFunctor is not within the image "
170 "cache bounds";
171 throw IException(IException::Programmer, msg, _FILEINFO_);
172 }
173 m_camera->Sensor::setTime(et);
174 if (!m_camera->Sensor::SetGround(surfacePoint, false)) {
175 IString msg = "Sensor::SetGround failed for surface point in LineScanCameraGroundMap.cpp"
176 "SensorSurfacePointDistanceFunctor";
177 }
178 m_camera->instrumentPosition(s);
179 m_camera->Coordinate(p);
180 return sqrt((s[0] - p[0]) * (s[0] - p[0]) +
181 (s[1] - p[1]) * (s[1] - p[1]) +
182 (s[2] - p[2]) * (s[2] - p[2]) ); //distance
183 }
184
185 private:
186 SurfacePoint surfacePoint;
187 Camera* m_camera;
188};
189
190
191namespace Isis {
192
198
199
204
213 const Longitude &lon) {
214 Distance radius(p_camera->LocalRadius(lat, lon));
215
216 if (radius.isValid()) {
217 return SetGround(SurfacePoint(lat, lon, radius));
218 }
219 else {
220 return false;
221 }
222 }
223
224
233 bool LineScanCameraGroundMap::SetGround(const SurfacePoint &surfacePoint, const int &approxLine) {
234 FindFocalPlaneStatus status = FindFocalPlane(approxLine, surfacePoint);
235 if (status == Success) return true;
236 //if(status == Failure) return false;
237 return false;
238 }
239
240
248 FindFocalPlaneStatus status = FindFocalPlane(-1, surfacePoint);
249
250 if (status == Success) return true;
251
252 return false;
253 }
254
255
256 double LineScanCameraGroundMap::FindSpacecraftDistance(int line,
257 const SurfacePoint &surfacePoint) {
258
259 CameraDetectorMap *detectorMap = p_camera->DetectorMap();
260 detectorMap->SetParent(p_camera->ParentSamples() / 2, line);
261 if (!p_camera->Sensor::SetGround(surfacePoint, false)) {
262 return DBL_MAX;
263 }
264
265 return p_camera->SlantDistance();
266 }
267
268
269 LineScanCameraGroundMap::FindFocalPlaneStatus
270 LineScanCameraGroundMap::FindFocalPlane(const int &approxLine,
271 const SurfacePoint &surfacePoint) {
272
273 //CameraDistortionMap *distortionMap = p_camera->DistortionMap();
274 //CameraFocalPlaneMap *focalMap = p_camera->FocalPlaneMap();
275
276 double approxTime = 0;
277 double approxOffset = 0;
278 double lookC[3] = {0.0, 0.0, 0.0};
279 double ux = 0.0;
280 double uy = 0.0;
281 //double dx = 0.0, dy = 0.0;
282 //double s[3], p[3];
283 const double cacheStart = p_camera->Spice::cacheStartTime().Et();
284 const double cacheEnd = p_camera->Spice::cacheEndTime().Et();
285
286 double lineRate = ((LineScanCameraDetectorMap *)p_camera->DetectorMap())->LineRate(); //line rate
287
288 if (lineRate == 0.0) return Failure;
289
290 LineOffsetFunctor offsetFunc(p_camera,surfacePoint);
291 SensorSurfacePointDistanceFunctor distanceFunc(p_camera,surfacePoint);
292
293 // METHOD #1
294 // Use the line given as a start point for the secant method root search.
295 if (approxLine >= 0.5) {
296
297 // convert the approxLine to an approximate time and offset
298 p_camera->DetectorMap()->SetParent(p_camera->ParentSamples() / 2.0, approxLine);
299 approxTime = p_camera->time().Et();
300
301 approxOffset = offsetFunc(approxTime);
302
303 // Check to see if there is no need to improve this root, it's good enough
304 if (fabs(approxOffset) < 1e-2) {
305 p_camera->Sensor::setTime(approxTime);
306 // check to make sure the point isn't behind the planet
307 if (!p_camera->Sensor::SetGround(surfacePoint, true)) {
308 return Failure;
309 }
310 p_camera->Sensor::LookDirection(lookC);
311 ux = p_camera->FocalLength() * lookC[0] / lookC[2];
312 uy = p_camera->FocalLength() * lookC[1] / lookC[2];
313
314 p_focalPlaneX = ux;
315 p_focalPlaneY = uy;
316
317 return Success;
318 }
319
320 double fl, fh, xl, xh;
321
322 // starting times for the secant method, kept within the domain of the cache
323 xh = approxTime;
324 if (xh + lineRate < cacheEnd) {
325 xl = xh + lineRate;
326 }
327 else {
328 xl = xh - lineRate;
329 }
330
331 // starting offsets
332 fh = approxOffset; //the first is already calculated
333 fl = offsetFunc(xl);
334
335 // Iterate to refine the given approximate time that the instrument imaged the ground point
336 for (int j=0; j < 10; j++) {
337 if (fl-fh == 0.0) {
338 return Failure;
339 }
340 double etGuess = xl + (xh - xl) * fl / (fl - fh);
341
342 if (etGuess < cacheStart) etGuess = cacheStart;
343 if (etGuess > cacheEnd) etGuess = cacheEnd;
344
345 double f = offsetFunc(etGuess);
346
347
348 // elliminate the node farthest away from the current best guess
349 if (fabs( xl- etGuess) > fabs( xh - etGuess)) {
350 xl = etGuess;
351 fl = f;
352 }
353 else {
354 xh = etGuess;
355 fh = f;
356 }
357
358 // See if we converged on the point so set up the undistorted focal plane values and return
359 if (fabs(f) < 1e-2) {
360 p_camera->Sensor::setTime(approxTime);
361 // check to make sure the point isn't behind the planet
362 if (!p_camera->Sensor::SetGround(surfacePoint, true)) {
363 return Failure;
364 }
365 p_camera->Sensor::LookDirection(lookC);
366 ux = p_camera->FocalLength() * lookC[0] / lookC[2];
367 uy = p_camera->FocalLength() * lookC[1] / lookC[2];
368
369 p_focalPlaneX = ux;
370 p_focalPlaneY = uy;
371
372 return Success;
373 }
374 } // End itteration using a guess
375 // return Failure; // Removed to let the lagrange method try to find the line if secant fails
376 } // End use a guess
377
378
379 // METHOD #2
380 // The guess or middle line did not work so try estimating with a quadratic
381 // The offsets are typically quadratic, so three points will be used to approximate a quadratic
382 // as a first order attempt to find the root location(s)
383
384 // The three nodes to be used to approximate the quadratic
385 double offsetNodes[3];
386 double timeNodes[3];
387 double timeAverage;
388 double scale;
389 QList<double> root;
390 QList<double> offset;
391 QList<double> dist;
392
393 timeNodes[0] = cacheStart;
394 timeNodes[2] = cacheEnd;
395 timeNodes[1] = (cacheStart+cacheEnd) / 2.0; // middle time
396
397 double quadPoly[3];
398 double temp;
399
400 for (int i=0; i<3; i++) {
401 offsetNodes[i] = offsetFunc(timeNodes[i]);
402 }
403
404 // centralize and normalize the data for stability in root finding
405 timeAverage = (timeNodes[0] + timeNodes[1] + timeNodes[2]) / 3.0;
406 timeNodes[0] -= timeAverage;
407 timeNodes[1] -= timeAverage;
408 timeNodes[2] -= timeAverage;
409
410 scale = 1.0 / sqrt((timeNodes[0] - timeNodes[2]) *
411 (timeNodes[0] - timeNodes[2]) +
412 (offsetNodes[0] - offsetNodes[2]) *
413 (offsetNodes[0] - offsetNodes[2]));
414
415 timeNodes[0] *= scale;
416 timeNodes[1] *= scale;
417 timeNodes[2] *= scale;
418
419 offsetNodes[0] *= scale;
420 offsetNodes[1] *= scale;
421 offsetNodes[2] *= scale;
422
423 // Use lagrange interpolating polynomials to find the coefficients of the quadratic,
424 // there are many ways to do this; I chose to do it this way because it is pretty straight
425 // forward and cheap
426 quadPoly[0] = quadPoly[1] = quadPoly[2] = 0.0;
427
428 temp = offsetNodes[0] / ((timeNodes[0] - timeNodes[1]) * (timeNodes[0] - timeNodes[2]));
429 quadPoly[0] += temp;
430 quadPoly[1] += temp * (-timeNodes[1] - timeNodes[2]);
431 quadPoly[2] += temp * timeNodes[1] * timeNodes[2];
432
433 temp = offsetNodes[1] / ((timeNodes[1] - timeNodes[0]) * (timeNodes[1] - timeNodes[2]));
434 quadPoly[0] += temp;
435 quadPoly[1] += temp * (-timeNodes[0] - timeNodes[2]);
436 quadPoly[2] += temp * timeNodes[0] * timeNodes[2];
437
438 temp = offsetNodes[2] / ((timeNodes[2] - timeNodes[0]) * (timeNodes[2] - timeNodes[1]));
439 quadPoly[0] += temp;
440 quadPoly[1] += temp * (-timeNodes[0] - timeNodes[1]);
441 quadPoly[2] += temp * timeNodes[0] * timeNodes[1];
442
443 // Now that we have the coefficients of the quadratic look for roots
444 // (see Numerical Recipes Third Edition page 227)
445 temp = quadPoly[1] * quadPoly[1] - 4.0 * quadPoly[0] * quadPoly[2]; //discriminant
446
447 // THIS IS A PREMATURE FAILURE RETURN. IT SHOULD TRY THE NEXT METHON BEFORE FAILING
448 if (temp < 0.0) {
449 return Failure; // there are apparently not any real roots on this image
450 }
451
452 if (quadPoly[1] >= 0.0) {
453 temp = -0.5 * (quadPoly[1] + sqrt(temp));
454 }
455 else {
456 temp = -0.5 * (quadPoly[1] - sqrt(temp));
457 }
458
459 if (quadPoly[0] != 0.0) {
460 root.push_back(temp/quadPoly[0]);
461 }
462
463 if (quadPoly[2] != 0.0) {
464 root.push_back(quadPoly[2]/temp);
465 }
466
467 // check to see if the roots are in the time interval of the cache
468 for (int i=root.size()-1; i>=0; i--) {
469 if ( root[i] < timeNodes[0] || root[i] > timeNodes[2] ) {
470 root.removeAt(i);
471 }
472 }
473
474 // return the calculated roots to the original system
475 for (int i=0; i<root.size(); i++) {
476 root[i] = root[i]/scale + timeAverage;
477 }
478
479 // THIS IS A PREMATURE FAILURE RETURN. IT SHOULD TRY THE NEXT METHON BEFORE FAILING
480 if (root.size() == 0) {
481 return Failure; // there are apparently not any roots on this image
482 }
483
484 // At the time of this writing ISIS made no attempt to support any sensors that were not "1 to 1".
485 // Meaning they imaged the same point on the ground in multiple lines of the image
486 // therefore we must somehow reduce multiple possible roots to a single one, the legacy
487 // code (replaced with this code) did this based on distance from the sensor to the target
488 // the shortest distance being the winner. For legacy consistency I have used the same logic below.
489
490 for (int i=0; i<root.size(); i++) { // Offset/dist calculation loop
491 dist << distanceFunc(root[i]);
492 offset << offsetFunc(root[i]);
493 }
494
495 // Save the root with the smallest dist
496 {
497 int j=0;
498 for (int i=1; i<root.size(); i++) {
499 if (dist[i] < dist[j]) j=i;
500 }
501
502 approxTime = root[j]; // Now we have our time
503 approxOffset = offset[j]; // The offsets are saved to avoid recalculating it later
504 }
505
506 if (fabs(approxOffset) < 1.0e-2) { // No need to iteratively improve this root, it's good enough
507 p_camera->Sensor::setTime(approxTime);
508
509 // Check to make sure the point isn't behind the planet
510 if (!p_camera->Sensor::SetGround(surfacePoint, true)) {
511 return Failure;
512 }
513
514 p_camera->Sensor::LookDirection(lookC);
515 ux = p_camera->FocalLength() * lookC[0] / lookC[2];
516 uy = p_camera->FocalLength() * lookC[1] / lookC[2];
517
518 p_focalPlaneX = ux;
519 p_focalPlaneY = uy;
520
521 return Success;
522 }
523
524
525 // METHOD #3
526 // Estimated line and quadratic approximation insufficient, try Brent's method
527 // The offsets are typically quadratic, so three points will be used to approximate a quadratic
528 // as a first order attempt to find the root location(s)
529
530 // The above sections are sufficient for finding the correct times for the vast majority of
531 // back projection solutions. The following section exists for those few particularly
532 // stuborn problems. They are typically characterised by being significantly non-quadratic.
533 // Examples mostly include images with very long exposure times.
534 //
535 // Further, while the preceeding sections is intended to be fast, this section is intended to be
536 // thurough. Brents method (Numerical Recipies 454 - 456) will be used to find all the roots,
537 // that are bracketed by the five points defined in the quadratic solution method above.
538 // The root with the shortest distance to the camera will be returned
539
540 // Get everything ordered for iteration combine and sort the five already defined points
541 QList <QList <double> > pts;
542
543 for (int i=0; i<3; i++) {
544 QList <double> pt;
545 pt << timeNodes[i] / scale + timeAverage;
546 pt << offsetNodes[i] / scale;
547 pts << pt;
548 }
549
550 for (int i=0; i<root.size(); i++) {
551 QList <double> pt;
552 pt << root[i];
553 pt << offset[i];
554 pts << pt;
555 }
556
557 qSort(pts.begin(), pts.end(), ptXLessThan);
558
559 root.clear();
560 for (int i=1; i<pts.size(); i++) {
561 // If the signs of the two offsets are not the same they bracket at least one root
562 if ( (pts[i-1][1] > 0) - (pts[i-1][1] < 0) != (pts[i][1] > 0) - (pts[i][1] < 0) ) {
563 double temp;
564 if (FunctionTools::brentsRootFinder <LineOffsetFunctor> (offsetFunc, pts[i-1], pts[i],
565 1.0e-3, 200, temp)) {
566 root << temp;
567 }
568 }
569 }
570
571 // Discard any roots that are looking through the planet
572 for (int i = root.size()-1; i>=0; i--) {
573 p_camera->Sensor::setTime(root[i]);
574 //check to make sure the point isn't behind the planet
575 if (!p_camera->Sensor::SetGround(surfacePoint, true)) {
576 root.removeAt(i);
577 }
578 }
579
580 // If none of the roots remain...
581 if (root.size() == 0) {
582 return Failure;
583 }
584
585 // Choose from the remaining roots, the solution with the smallest distance to target
586 dist.clear();
587 offset.clear();
588 for (int i=0; i<root.size(); i++) { // Offset/dist calculation loop
589 dist << distanceFunc(root[i]);
590 offset << offsetFunc(root[i]);
591 }
592
593 // Save the root with the smallest dist
594 {
595 int j=0;
596 for (int i=1; i<root.size(); i++) {
597 if (dist[i] < dist[j]) j=i;
598 }
599
600 p_camera->Sensor::setTime(root[j]);
601 }
602
603 // No need to make sure the point isn't behind the planet, it was done above
604 p_camera->Sensor::LookDirection(lookC);
605 ux = p_camera->FocalLength() * lookC[0] / lookC[2];
606 uy = p_camera->FocalLength() * lookC[1] / lookC[2];
607
608 p_focalPlaneX = ux;
609 p_focalPlaneY = uy;
610
611 return Success;
612 }
613}
614
615
616bool ptXLessThan(const QList<double> l1, const QList<double> l2) {
617 return l1[0] < l2[0];
618}
Convert between parent image coordinates and detector coordinates.
virtual bool SetParent(const double sample, const double line)
Compute detector position from a parent image coordinate.
virtual bool SetUndistortedFocalPlane(double ux, double uy)
Compute distorted focal plane x/y.
double FocalPlaneX() const
Gets the x-value in the focal plane coordinate system.
double FocalPlaneY() const
Gets the y-value in the focal plane coordinate system.
virtual bool SetFocalPlane(const double dx, const double dy)
Compute detector position (sample,line) from focal plane coordinates.
Convert between undistorted focal plane and ground coordinates.
double p_focalPlaneX
Camera's x focal plane coordinate.
Camera * p_camera
Camera.
double p_focalPlaneY
Camera's y focal plane coordinate.
CameraDetectorMap * DetectorMap()
Returns a pointer to the CameraDetectorMap object.
Definition Camera.cpp:2876
int ParentSamples() const
Returns the number of samples in the parent alphacube.
Definition Camera.cpp:2846
double FocalLength() const
Returns the focal length.
Definition Camera.cpp:2762
CameraDistortionMap * DistortionMap()
Returns a pointer to the CameraDistortionMap object.
Definition Camera.cpp:2856
CameraFocalPlaneMap * FocalPlaneMap()
Returns a pointer to the CameraFocalPlaneMap object.
Definition Camera.cpp:2866
Distance measurement, usually in meters.
Definition Distance.h:34
Isis exception class.
Definition IException.h:91
@ Programmer
This error is for when a programmer made an API call that was illegal.
Definition IException.h:146
Adds specific functionality to C++ strings.
Definition IString.h:165
This class is designed to encapsulate the concept of a Latitude.
Definition Latitude.h:51
Convert between parent image coordinates and detector coordinates.
LineScanCameraGroundMap(Camera *cam)
Constructor.
virtual bool SetGround(const Latitude &lat, const Longitude &lon)
Compute undistorted focal plane coordinate from ground position.
This class is designed to encapsulate the concept of a Longitude.
Definition Longitude.h:40
virtual double SlantDistance() const
Return the distance between the spacecraft and surface point in kmv.
Definition Sensor.cpp:646
void Coordinate(double p[3]) const
Returns the x,y,z of the surface intersection in BodyFixed km.
Definition Sensor.cpp:198
Distance LocalRadius() const
Returns the local radius at the intersection point.
Definition Sensor.cpp:269
virtual iTime cacheEndTime() const
Accessor method for the cache end time.
Definition Spice.cpp:765
void instrumentPosition(double p[3]) const
Returns the spacecraft position in body-fixed frame km units.
Definition Spice.cpp:829
virtual iTime time() const
Returns the ephemeris time in seconds which was used to obtain the spacecraft and sun positions.
Definition Spice.cpp:891
virtual iTime cacheStartTime() const
Accessor method for the cache start time.
Definition Spice.cpp:750
This class defines a body-fixed surface point.
double operator()(double et)
Compute the number of lines between the current line (i.e., the line imaged at the et as set in the c...
This is free and unencumbered software released into the public domain.
Definition Apollo.h:16
Namespace for the standard library.