Isis 3 Programmer Reference
VimsSkyMap.cpp
Go to the documentation of this file.
1 
20 #include "VimsSkyMap.h"
21 
22 #include <iostream>
23 #include <iomanip>
24 
25 #include <QDebug>
26 
27 #include "Camera.h"
28 #include "FileName.h"
29 #include "IException.h"
30 #include "IString.h"
31 #include "iTime.h"
32 #include "LeastSquares.h"
33 #include "PolynomialBivariate.h"
34 #include "Preference.h"
35 #include "SpecialPixel.h"
36 
37 using namespace std;
38 
39 namespace Isis {
48  VimsSkyMap::VimsSkyMap(Camera *parent, Pvl &lab) :
49  CameraSkyMap(parent) {
50  }
51 
52 
69  void VimsSkyMap::Init(Pvl &lab) {
70  PvlGroup inst = lab.findGroup("Instrument", Pvl::Traverse);
71 
72  // Vis or IR
73  p_channel = (QString) inst ["Channel"];
74  // Get the start time in et
75  QString stime = (QString) inst ["NativeStartTime"];
76  QString intTime = stime.split(".").first();
77  stime = stime.split(".").last();
78 
79  p_etStart = p_camera->getClockTime(intTime).Et();
80  p_etStart += toDouble(stime) / 15959.0;
81 
82  //----------------------------------------------------------------------
83  // Because of inaccuracy with the 15 Mhz clock, the IR exposure and
84  // interline delay need to be adjusted.
85  //----------------------------------------------------------------------
86  p_irExp = (toDouble(inst ["ExposureDuration"][0]) * 1.01725) / 1000.;
87  p_visExp = toDouble(inst ["ExposureDuration"][1]) / 1000.;
88  p_interlineDelay = (toDouble(inst ["InterlineDelayDuration"]) * 1.01725) / 1000.;
89 
90  // Get summation mode
91  QString sampMode = QString((QString)inst ["SamplingMode"]).toUpper();
92 
93  // Get sample/line offsets
94  int sampOffset = inst ["XOffset"];
95  int lineOffset = inst ["ZOffset"];
96 
97  // Get swath width/length which will be image size unless occultation image
98  p_swathWidth = inst ["SwathWidth"];
99  p_swathLength = inst ["SwathLength"];
100 
101  if (p_channel == "VIS") {
102  if (sampMode == "NORMAL") {
103  p_xPixSize = p_yPixSize = 0.00051;
104  p_xBore = p_yBore = 31;
105  p_camSampOffset = sampOffset - 1;
106  p_camLineOffset = lineOffset - 1;
107 
108  }
109  else if (sampMode == "HI-RES") {
110  p_xPixSize = p_yPixSize = 0.00051 / 3.0;
111  p_xBore = p_yBore = 94;
112  //New as of 2009-08-04 per Dyer Lytle's email
113  // Old Values:p_camSampOffset = 3 * (sampOffset - 1) + p_swathWidth;
114  // p_camLineOffset = 3 * (lineOffset - 1) + p_swathLength;
115  p_camSampOffset = (3 * (sampOffset + p_swathWidth / 2)) - p_swathWidth / 2;
116  p_camLineOffset = (3 * (lineOffset + p_swathLength / 2)) - p_swathLength / 2;
117  }
118  else {
119  string msg = "Unsupported SamplingMode [" + IString(sampMode) + "]";
121  }
122  }
123  else if (p_channel == "IR") {
124  if (sampMode == "NORMAL") {
125  p_xPixSize = p_yPixSize = 0.000495;
126  p_xBore = p_yBore = 31;
127  p_camSampOffset = sampOffset - 1;
128  p_camLineOffset = lineOffset - 1;
129  }
130  else if (sampMode == "HI-RES") {
131  p_xPixSize = 0.000495 / 2.0;
132  p_yPixSize = 0.000495;
133  p_xBore = 62.5;
134  p_yBore = 31;
135  p_camSampOffset = 2 * ((sampOffset - 1) + ((p_swathWidth - 1) / 4));
136  p_camLineOffset = lineOffset - 1;
137  }
138  else {
139  string msg = "Unsupported SamplingMode [" + IString(sampMode) + "]";
141  }
142  }
143 
144  // Calculate ra/dec maps
145  for(int line = 0; line < p_camera->ParentLines(); line++) {
146  for(int samp = 0; samp < p_camera->ParentSamples(); samp++) {
147  p_raMap[line][samp] = Isis::NULL8;
148  p_decMap[line][samp] = Isis::NULL8;
149  }
150  }
151 
152  //---------------------------------------------------------------------
153  // Loop for each pixel in cube, get pointing information and calculate
154  // control point (line,sample,lat,lon) for later use in latlon_to_linesamp.
155  //---------------------------------------------------------------------
156  p_minRa = 99999.;
157  p_minDec = 99999.;
158  p_maxRa = -99999.;
159  p_maxDec = -99999.;
160 
161  p_camera->IgnoreProjection(true);
162  for(int line = 0; line < p_camera->ParentLines(); line++) {
163 
164  // VIS exposure is for a single line. According to SIS,
165  // NATIVE_START_TIME is for the first pixel of the IR exposure.
166  // "The offset from IR start to VIS start is calculated by
167  // IrExposMsec - VisExposMsec)/2".
168  // This needs to be moved to forward routine
169 
170  if(p_channel == "VIS") {
171  double et = ((double)p_etStart + (((p_irExp * p_swathWidth) - p_visExp) / 2.)) +
172  ((line + 0.5) * p_visExp);
173  p_camera->setTime(et);
174  }
175 
176  for(int samp = 0; samp < p_camera->ParentSamples(); samp++) {
177  if(p_channel == "IR") {
178  double et = (double)p_etStart +
179  (line * p_camera->ParentSamples() * p_irExp) +
180  (line * p_interlineDelay) + ((samp + 0.5) * p_irExp);
181  p_camera->setTime(et);
182  }
183 
184  p_camera->SetImage((double) samp + 1, (double)line + 1);
185  double ra = p_camera->RightAscension();
186  double dec = p_camera->Declination();
187  if(ra < p_minRa) p_minRa = ra;
188  if(ra > p_maxRa) p_maxRa = ra;
189  if(dec < p_minDec) p_minDec = dec;
190  if(dec > p_maxDec) p_maxDec = dec;
191  p_raMap[line][samp] = ra;
192  p_decMap[line][samp] = dec;
193  }
194  }
195  p_camera->IgnoreProjection(false);
196 
197  }
198 
199 
212  bool VimsSkyMap::SetFocalPlane(const double ux, const double uy,
213  const double uz) {
214  p_ux = ux;
215  p_uy = uy;
216  p_uz = uz;
217 
218  double imgSamp = ux;
219  double imgLine = uy;
220 
221  if((imgLine < 0.5) || (imgLine > p_camera->ParentLines() + 0.5) ||
222  (imgSamp < 0.5) || (imgSamp > p_camera->ParentSamples() + 0.5)) {
223  return false;
224  }
225  imgLine--;
226  imgSamp--;
227 
228  // does interline_delay & exposure-duration account for summing modes?
229  // if not, won't use p_parentLine/p_parentSample
230  double et = 0.;
231  if(p_channel == "VIS") {
232  et = (p_etStart + ((p_irExp * p_swathWidth) - p_visExp) / 2.) +
233  ((imgLine + 0.5) * p_visExp);
234  }
235  else if(p_channel == "IR") {
236  et = (double)p_etStart +
237  (imgLine * p_camera->ParentSamples() * p_irExp) +
238  (imgLine * p_interlineDelay) + ((imgSamp + 0.5) * p_irExp);
239  }
240  p_camera->setTime(et);
241 
242  SpiceDouble lookC[3];
243  LookDirection(lookC);
244 
245  SpiceDouble unitLookC[3];
246  vhat_c(lookC, unitLookC);
247  return p_camera->SetLookDirection(unitLookC);
248  }
249 
258  bool VimsSkyMap::SetSky(const double ra, const double dec) {
259  if(ra < p_minRa || ra > p_maxRa ||
260  dec < p_minDec || dec > p_maxDec) {
261  return false;
262  }
263  // Find closest points ??? what tolerance ???
264  double minDist = 9999.;
265  int minSamp = -1;
266  int minLine = -1;
267 
268  for(int line = 0; line < p_camera->ParentLines(); line++) {
269 
270  for(int samp = 0; samp < p_camera->ParentSamples(); samp++) {
271  double mapRa = p_raMap[line][samp];
272  if(mapRa == Isis::NULL8) continue;
273  double mapDec = p_decMap[line][samp];
274  if(mapDec == Isis::NULL8) continue;
275  // If on boundary convert lons. If trying to find 360, convert
276  // lons on other side of meridian to values greater than 360. If
277  // trying to find 1.0, convert lons on other side to negative numbers.
278  if(abs(mapRa - ra) > 180) {
279  if((ra - mapRa) > 0) {
280  mapRa = 360. + mapRa;
281  }
282  else if((ra - mapRa) < 0) {
283  mapRa = mapRa - 360.;
284  }
285  }
286  double dist = ((ra - mapRa) * (ra - mapRa)) +
287  ((dec - mapDec) * (dec - mapDec));
288  if(dist < minDist) {
289  minDist = dist;
290  minSamp = samp;
291  minLine = line;
292  }
293  }
294  }
295 
296  //-----------------------------------------------------------------
297  // If dist is less than some ??? tolerance ??? this is the
298  // closest point. Use this point and surrounding 8 pts as
299  // control pts.
300  //----------------------------------------------------------------
301  if(minDist >= 9999.) return false;
302 
303  //-------------------------------------------------------------
304  // Set-up for LU decomposition (least2 fit).
305  // Assume we will have 9 control points, this may not be true
306  // and will need to be adjusted before the final solution.
307  //-------------------------------------------------------------
308  PolynomialBivariate sampBasis(1);
309  PolynomialBivariate lineBasis(1);
310  LeastSquares sampLsq(sampBasis);
311  LeastSquares lineLsq(lineBasis);
312  vector<double> known(2);
313 
314  for(int line = minLine - 1; line < minLine + 2; line++) {
315  if(line < 0 || line > p_camera->ParentLines() - 1) continue;
316  for(int samp = minSamp - 1; samp < minSamp + 2; samp++) {
317  // Check for edges
318  if(samp < 0 || samp > p_camera->ParentSamples() - 1) continue;
319 
320  double mapRa = p_raMap[line][samp];
321  double mapDec = p_decMap[line][samp];
322  if((mapRa == Isis::NULL8) || (mapDec == Isis::NULL8)) continue;
323 
324  // If on boundary convert lons. If trying to find 360, convert
325  // lons on other side of meridian to values greater than 360. If
326  // trying to find 1.0, convert lons on other side to negative numbers.
327  if(abs(mapRa - ra) > 180) {
328  if((ra - mapRa) > 0) {
329  mapRa = 360. + mapRa;
330  }
331  else if((ra - mapRa) < 0) {
332  mapRa = mapRa - 360.;
333  }
334  }
335 
336  known[0] = mapDec;
337  known[1] = mapRa;
338  sampLsq.AddKnown(known, samp + 1);
339  lineLsq.AddKnown(known, line + 1);
340  }
341  }
342  if(sampLsq.Knowns() < 3) return false;
343 
344  sampLsq.Solve();
345  lineLsq.Solve();
346 
347  // Solve for new sample position
348  known[0] = dec;
349  known[1] = ra;
350  double inSamp = sampLsq.Evaluate(known);
351  double inLine = lineLsq.Evaluate(known);
352 
353  if(inSamp < 0 || inSamp > p_camera->ParentSamples() + 0.5 ||
354  inLine < 0 || inLine > p_camera->ParentLines() + 0.5) {
355  return false;
356  }
357 
358  p_camera->IgnoreProjection(true);
359  p_camera->SetImage(inSamp, inLine);
360  p_camera->IgnoreProjection(false);
361  p_focalPlaneX = inSamp;
362  p_focalPlaneY = inLine;
363 
364  return true;
365  }
366 
367 
383  void VimsSkyMap::LookDirection(double v[3]) {
384  double x = p_ux - 1. + p_camSampOffset;
385  double y = p_uy - 1. + p_camLineOffset;
386 
387  // Compute pointing angles based on pixel size separation
388  double theta = Isis::HALFPI - (y - p_yBore) * p_yPixSize;
389  double phi = (-1. * Isis::HALFPI) + (x - p_xBore) * p_xPixSize;
390 
391  v[0] = sin(theta) * cos(phi);
392  v[1] = cos(theta);
393  v[2] = sin(theta) * sin(phi) / -1.;
394  }
395 }
double RightAscension()
Returns the right ascension angle (sky longitude).
Definition: Sensor.cpp:576
void LookDirection(double v[3])
Determines the look direction in the camera coordinate system.
Definition: VimsSkyMap.cpp:383
double p_decMap[64][64]
Declination map.
Definition: VimsSkyMap.h:108
double p_focalPlaneY
Undistorted y value for the focal plane.
Definition: CameraSkyMap.h:77
double p_minRa
Minimum right ascension.
Definition: VimsSkyMap.h:103
double p_maxDec
Maximum declination.
Definition: VimsSkyMap.h:106
PvlGroupIterator findGroup(const QString &name, PvlGroupIterator beg, PvlGroupIterator end)
Find a group with the specified name, within these indexes.
Definition: PvlObject.h:141
virtual bool SetFocalPlane(const double ux, const double uy, const double uz)
Compute ground position from focal plane coordinate.
Definition: VimsSkyMap.cpp:212
SpiceDouble p_etStart
Start ephemeris time.
Definition: VimsSkyMap.h:79
int Knowns() const
The number of knowns (or times AddKnown was invoked) linear combination of the variables.
Definition: LeastSquares.h:145
const double HALFPI
The mathematical constant PI/2.
Definition: Constants.h:57
Nth degree Polynomial with two variables.
double p_uy
Distorted focal plane y, in millimeters.
Definition: VimsSkyMap.h:84
Namespace for the standard library.
double p_yPixSize
Y pixel size.
Definition: VimsSkyMap.h:88
double p_interlineDelay
InterlineDelayDuration keyword value from the instrument group of the labels, divided by 1000...
Definition: VimsSkyMap.h:81
virtual bool SetSky(const double ra, const double dec)
Sets the sky position to the given ra and dec.
Definition: VimsSkyMap.cpp:258
Search child objects.
Definition: PvlObject.h:170
bool SetLookDirection(const double v[3])
Sets the look direction of the spacecraft.
Definition: Sensor.cpp:156
double toDouble(const QString &string)
Global function to convert from a string to a double.
Definition: IString.cpp:164
This error is for when a programmer made an API call that was illegal.
Definition: IException.h:162
int ParentLines() const
Returns the number of lines in the parent alphacube.
Definition: Camera.cpp:2818
int p_swathLength
SwathLength keyword value from the instrument group of the labels.
Definition: VimsSkyMap.h:98
bool SetImage(const double sample, const double line)
Sets the sample/line values of the image to get the lat/lon values.
Definition: Camera.cpp:170
double Evaluate(const std::vector< double > &input)
Invokes the BasisFunction Evaluate method.
int p_swathWidth
SwathWidth keyword value from the instrument group of the labels.
Definition: VimsSkyMap.h:96
double p_maxRa
Maximum right ascension.
Definition: VimsSkyMap.h:104
double p_yBore
Y boresight.
Definition: VimsSkyMap.h:90
int p_camSampOffset
Sample offset.
Definition: VimsSkyMap.h:100
Contains multiple PvlContainers.
Definition: PvlGroup.h:57
double p_raMap[64][64]
Right ascension map.
Definition: VimsSkyMap.h:107
#define _FILEINFO_
Macro for the filename and line number.
Definition: IException.h:40
iTime getClockTime(QString clockValue, int sclkCode=-1, bool clockTicks=false)
This converts the spacecraft clock ticks value (clockValue) to an iTime.
Definition: Spice.cpp:977
double p_xBore
X boresight.
Definition: VimsSkyMap.h:89
int ParentSamples() const
Returns the number of samples in the parent alphacube.
Definition: Camera.cpp:2828
void setTime(const iTime &time)
By setting the time you essential set the position of the spacecraft and body as indicated in the cla...
Definition: Sensor.cpp:112
double p_visExp
VIS exposure duration, divided by 1000.
Definition: VimsSkyMap.h:94
Generic least square fitting class.
Definition: LeastSquares.h:115
double p_ux
Distorted focal plane x, in millimeters.
Definition: VimsSkyMap.h:83
Container for cube-like labels.
Definition: Pvl.h:135
void Init(Pvl &lab)
Initialize vims sky model.
Definition: VimsSkyMap.cpp:69
void AddKnown(const std::vector< double > &input, double expected, double weight=1.0)
Invoke this method for each set of knowns.
QString p_channel
Channel keyword value from the instrument group of the labels.
Definition: VimsSkyMap.h:92
double p_focalPlaneX
Undistorted x value for the focal plane.
Definition: CameraSkyMap.h:76
double p_minDec
Minimum declination.
Definition: VimsSkyMap.h:105
Isis exception class.
Definition: IException.h:107
Adds specific functionality to C++ strings.
Definition: IString.h:181
Namespace for ISIS/Bullet specific routines.
Definition: Apollo.h:31
Convert between undistorted focal plane and ra/dec coordinates.
Definition: CameraSkyMap.h:48
Camera * p_camera
The main camera to calculate distortions on.
Definition: CameraSkyMap.h:72
void IgnoreProjection(bool ignore)
Set whether or not the camera should ignore the Projection.
Definition: Camera.cpp:2937
double p_xPixSize
X pixel size.
Definition: VimsSkyMap.h:87
double Et() const
Returns the ephemeris time (TDB) representation of the time as a double.
Definition: iTime.h:139
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...
double p_irExp
IR exposure duration, divided by 1000.
Definition: VimsSkyMap.h:95
double p_uz
Distorted focal plane z, in millimeters.
Definition: VimsSkyMap.h:85
int p_camLineOffset
Line offset.
Definition: VimsSkyMap.h:101
double Declination()
Returns the declination angle (sky latitude).
Definition: Sensor.cpp:587