Official websites use .gov
A .gov website belongs to an official government organization in the United States.

Secure .gov websites use HTTPS
A lock ( ) or https:// means you’ve safely connected to the .gov website. Share sensitive information only on official, secure websites.

Isis 3 Programmer Reference
ReseauDistortionMap.cpp
1
5
6/* SPDX-License-Identifier: CC0-1.0 */
7#include "ReseauDistortionMap.h"
8#include "LeastSquares.h"
9#include "BasisFunction.h"
10#include "CameraFocalPlaneMap.h"
11#include "Pvl.h"
12#include "Statistics.h"
13#include <float.h>
14
15using namespace std;
16namespace Isis {
28 ReseauDistortionMap::ReseauDistortionMap(Camera *parent, Pvl &labels, const QString &fname) :
29 CameraDistortionMap(parent, 1.0) {
30 // Set up distortion coefficients
31 Pvl mast(fname);
32 PvlGroup dim = mast.findGroup("Dimensions");
33 p_distortedLines = dim.findKeyword("DistortedLines");
34 p_distortedSamps = dim.findKeyword("DistortedSamples");
35 p_undistortedLines = dim.findKeyword("UndistortedLines");
36 p_undistortedSamps = dim.findKeyword("UndistortedSamples");
37 PvlGroup mastRes = mast.findGroup("MasterReseaus");
38 PvlKeyword mline = mastRes.findKeyword("Line");
39 PvlKeyword msamp = mastRes.findKeyword("Sample");
40 p_numRes = mline.size();
41 if(mline.size() != msamp.size()) {
42 string msg = "The number of lines and samples for the master reseaus are";
43 msg += "not equal, the data file may be bad";
44 throw IException(IException::User, msg, _FILEINFO_);
45 }
46 for(int i = 0; i < p_numRes; i++) {
47 p_mlines.push_back(toDouble(mline[i]));
48 p_msamps.push_back(toDouble(msamp[i]));
49 }
50 p_pixelPitch = parent->PixelPitch();
51
52 PvlGroup refRes = labels.findGroup("Reseaus", Pvl::Traverse);
53 PvlKeyword rline = refRes.findKeyword("Line");
54 PvlKeyword rsamp = refRes.findKeyword("Sample");
55 if(rline.size() != rsamp.size()) {
56 string msg = "The number of lines and samples for the refined reseaus are";
57 msg += "not equal, the data file may be bad";
58 throw IException(IException::User, msg, _FILEINFO_);
59 }
60 for(int i = 0; i < p_numRes; i++) {
61 p_rlines.push_back(toDouble(rline[i]));
62 p_rsamps.push_back(toDouble(rsamp[i]));
63 }
64 if(p_mlines.size() != p_rlines.size()) {
65 string msg = "The number of master reseaus and refined reseaus";
66 msg += "do not appear to be equal";
67 throw IException(IException::User, msg, _FILEINFO_);
68 }
69 }
70
81 const double dy) {
82 p_focalPlaneX = dx;
83 p_focalPlaneY = dy;
84
85 // Convert distorted x,y position to a sample, line position
86 double focalSamp = dx / p_pixelPitch +
87 p_camera->FocalPlaneMap()->DetectorSampleOrigin();
88 double focalLine = dy / p_pixelPitch +
89 p_camera->FocalPlaneMap()->DetectorLineOrigin();
90
91 // Find distance from input point to all nominal reseaus
92 double distances[p_numRes], wt[5];
93 int closepts[5];
94 double ldiffsq, sdiffsq;
95 for(int i = 0; i < p_numRes; i++) {
96 sdiffsq = (focalSamp - p_rsamps[i]) * (focalSamp - p_rsamps[i]);
97 ldiffsq = (focalLine - p_rlines[i]) * (focalLine - p_rlines[i]);
98 distances[i] = ldiffsq + sdiffsq;
99 }
100
101 // Get 5 closest reseaus to the input point
102 for(int ifpt = 0; ifpt < 5; ifpt++) {
103 int imin = 0;
104 for(int i = 0; i < p_numRes; i++) {
105 if(distances[i] < distances[imin]) imin = i;
106 }
107 closepts[ifpt] = imin;
108 wt[ifpt] = distances[imin];
109 distances[imin] = DBL_MAX; /* Arbitrary big number */
110 }
111
112 // Make sure 5 closest points are not colinear
113 Statistics lstats, sstats;
114 double line[5], samp[5];
115 for(int k = 0; k < 5; k++) {
116 line[k] = p_rlines[closepts[k]];
117 samp[k] = p_rsamps[closepts[k]];
118 }
119 lstats.AddData(line, 5);
120 sstats.AddData(samp, 5);
121
122 // If they are colinear, return false
123 if(lstats.StandardDeviation() < 10.0 || sstats.StandardDeviation() < 10.0) {
124 return false;
125 }
126
127 // Weight each of the 5 closest points, and solve the system of equations
128 if(wt[0] > 0.0) {
129 double scale = wt[0];
130 double rfitlines[5], rfitsamps[5], mfitlines[5], mfitsamps[5];
131 for(int ifpt = 0; ifpt < 5; ifpt++) {
132 int index = closepts[ifpt];
133 rfitlines[ifpt] = p_rlines[index];
134 rfitsamps[ifpt] = p_rsamps[index];
135 mfitlines[ifpt] = p_mlines[index];
136 mfitsamps[ifpt] = p_msamps[index];
137 wt[ifpt] = scale / wt[ifpt];
138 }
139 BasisFunction bsX("bilinearInterpX", 3, 3);
140 BasisFunction bsY("bilinearInterpY", 3, 3);
141 LeastSquares lsqX(bsX);
142 LeastSquares lsqY(bsY);
143
144 vector<double> known;
145 known.resize(3);
146 for(int i = 0; i < 5; i++) {
147 known[0] = 1.0;
148 known[1] = rfitsamps[i];
149 known[2] = rfitlines[i];
150 lsqX.AddKnown(known, mfitsamps[i], wt[i]);
151 lsqY.AddKnown(known, mfitlines[i], wt[i]);
152 }
153 lsqX.Solve();
154 lsqY.Solve();
155
156 known[1] = focalSamp;
157 known[2] = focalLine;
158
159 // Test to make sure the point is inside of the image
160 double undistortedFocalSamp = lsqX.Evaluate(known);
161 double undistortedFocalLine = lsqY.Evaluate(known);
162 if(undistortedFocalSamp < 0.5) return false;
163 if(undistortedFocalLine < 0.5) return false;
164 if(undistortedFocalSamp > p_undistortedSamps + 0.5) return false;
165 if(undistortedFocalLine > p_undistortedLines + 0.5) return false;
166
167 // Convert undistorted sample, line position to an x,y position
168 p_undistortedFocalPlaneX = (undistortedFocalSamp - p_undistortedSamps
169 / 2.0) * p_pixelPitch;
170 p_undistortedFocalPlaneY = (undistortedFocalLine - p_undistortedLines
171 / 2.0) * p_pixelPitch;
172 }
173 else { // If the point passed in is a reseau...
174 int index = closepts[0];
176 * p_pixelPitch;
177 p_undistortedFocalPlaneY = (p_mlines[index] - p_undistortedLines / 2.0)
178 * p_pixelPitch;
179 }
180 return true;
181 }
182
193 const double uy) {
196
197 // Convert undistorted values to sample, line positions
198 double undistortedFocalSamp = ux / p_pixelPitch + p_undistortedSamps / 2.0;
199 double undistortedFocalLine = uy / p_pixelPitch + p_undistortedLines / 2.0;
200
201 // Find distance from input point to all nominal reseaus
202 double distances[p_numRes], wt[5];
203 int closepts[5];
204 double ldiffsq, sdiffsq;
205 for(int i = 0; i < p_numRes; i++) {
206 sdiffsq = (undistortedFocalSamp - p_msamps[i]) *
207 (undistortedFocalSamp - p_msamps[i]);
208 ldiffsq = (undistortedFocalLine - p_mlines[i]) *
209 (undistortedFocalLine - p_mlines[i]);
210 distances[i] = ldiffsq + sdiffsq;
211 }
212
213 // Get 5 closest reseaus to the input point
214 for(int ifpt = 0; ifpt < 5; ifpt++) {
215 int imin = 0;
216 for(int i = 0; i < p_numRes; i++) {
217 if(distances[i] < distances[imin]) imin = i;
218 }
219 closepts[ifpt] = imin;
220 wt[ifpt] = distances[imin];
221 distances[imin] = DBL_MAX; /* Arbitrary big number */
222 }
223
224 // Make sure 5 closest points are not colinear
225 Statistics lstats, sstats;
226 double line[5], samp[5];
227 for(int k = 0; k < 5; k++) {
228 line[k] = p_rlines[closepts[k]];
229 samp[k] = p_rsamps[closepts[k]];
230 }
231 lstats.AddData(line, 5);
232 sstats.AddData(samp, 5);
233 // If they are colinear, return false
234 if(lstats.StandardDeviation() < 10.0 || sstats.StandardDeviation() < 10.0) {
235 return false;
236 }
237
238 // Weight each of the 5 closest points and solve the system of equations
239 if(wt[0] > 0.0) {
240 double scale = wt[0];
241 double rfitlines[5], rfitsamps[5], mfitlines[5], mfitsamps[5];
242 for(int ifpt = 0; ifpt < 5; ifpt++) {
243 int index = closepts[ifpt];
244 mfitlines[ifpt] = p_mlines[index];
245 mfitsamps[ifpt] = p_msamps[index];
246 rfitlines[ifpt] = p_rlines[index];
247 rfitsamps[ifpt] = p_rsamps[index];
248 wt[ifpt] = scale / wt[ifpt];
249 }
250 BasisFunction bsX("bilinearInterpX", 3, 3);
251 BasisFunction bsY("bilinearInterpY", 3, 3);
252 LeastSquares lsqX(bsX);
253 LeastSquares lsqY(bsY);
254
255 vector<double> known;
256 known.resize(3);
257 for(int i = 0; i < 5; i++) {
258 known[0] = 1.0;
259 known[1] = mfitsamps[i];
260 known[2] = mfitlines[i];
261 lsqX.AddKnown(known, rfitsamps[i], wt[i]);
262 lsqY.AddKnown(known, rfitlines[i], wt[i]);
263 }
264 lsqX.Solve();
265 lsqY.Solve();
266
267 known[1] = undistortedFocalSamp;
268 known[2] = undistortedFocalLine;
269
270 // Test points to make sure they are in the image
271 double distortedFocalSamp = lsqX.Evaluate(known);
272 double distortedFocalLine = lsqY.Evaluate(known);
273 if(distortedFocalSamp < 0.5) return false;
274 if(distortedFocalLine < 0.5) return false;
275 if(distortedFocalSamp > p_undistortedSamps + 0.5) return false;
276 if(distortedFocalLine > p_undistortedLines + 0.5) return false;
277
278 // Convert distorted sample, line position back to an x,y position
279 p_focalPlaneX = (distortedFocalSamp -
280 p_camera->FocalPlaneMap()->DetectorSampleOrigin()) * p_pixelPitch;
281 p_focalPlaneY = (distortedFocalLine -
282 p_camera->FocalPlaneMap()->DetectorLineOrigin()) * p_pixelPitch;
283 }
284
285 else { // If the point passed in is a reseau...
286 int index = closepts[0];
287 p_focalPlaneX = (p_rsamps[index] -
288 p_camera->FocalPlaneMap()->DetectorSampleOrigin()) * p_pixelPitch;
289 p_focalPlaneY = (p_rlines[index] -
290 p_camera->FocalPlaneMap()->DetectorLineOrigin()) * p_pixelPitch;
291 }
292 return true;
293 }
294
295} // end namespace isis
Generic linear equation class.
CameraDistortionMap(Camera *parent, double zDirection=1.0)
Camera distortion map constructor.
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
Isis exception class.
Definition IException.h:91
@ User
A type of error that could only have occurred due to a mistake on the user's part (e....
Definition IException.h:126
Generic least square fitting class.
double Evaluate(const std::vector< double > &input)
Invokes the BasisFunction Evaluate method.
int Solve(Isis::LeastSquares::SolveMethod method=SVD)
After all the data has been registered through AddKnown, invoke this method to solve the system of eq...
void AddKnown(const std::vector< double > &input, double expected, double weight=1.0)
Invoke this method for each set of knowns.
PvlKeyword & findKeyword(const QString &name)
Find a keyword with a specified name.
Contains multiple PvlContainers.
Definition PvlGroup.h:41
Container for cube-like labels.
Definition Pvl.h:119
A single keyword-value pair.
Definition PvlKeyword.h:87
int size() const
Returns the number of values stored in this keyword.
Definition PvlKeyword.h:133
@ Traverse
Search child objects.
Definition PvlObject.h:158
PvlGroupIterator findGroup(const QString &name, PvlGroupIterator beg, PvlGroupIterator end)
Find a group with the specified name, within these indexes.
Definition PvlObject.h:129
std::vector< double > p_msamps
Master Reseau Locations.
std::vector< double > p_rsamps
Refined Reseau Locations.
double p_undistortedSamps
Dimensions of undistorted cube.
virtual bool SetFocalPlane(const double dx, const double dy)
Finds the undistorted x/y position of the given distorted point.
double p_pixelPitch
Pixel Pitch of parent Camera.
int p_numRes
Number of Reseaus.
ReseauDistortionMap(Camera *parent, Pvl &labels, const QString &fname)
Creates a ReseauDistortionMap object.
virtual bool SetUndistortedFocalPlane(const double ux, const double uy)
Finds the distorted x/y position of the given undistorted point.
double p_distortedSamps
Dimensions of distorted cube.
This class is used to accumulate statistics on double arrays.
Definition Statistics.h:93
void AddData(const double *data, const unsigned int count)
Add an array of doubles to the accumulators and counters.
double StandardDeviation() const
Computes and returns the standard deviation.
This is free and unencumbered software released into the public domain.
Definition Apollo.h:16
double toDouble(const QString &string)
Global function to convert from a string to a double.
Definition IString.cpp:149
Namespace for the standard library.