USGS

Isis 3.0 Object Programmers' Reference

 Home
DawnVirCamera.cpp
1 #include <cctype>
2 #include <iostream>
3 #include <iomanip>
4 #include <sstream>
5 #include <algorithm>
6 
7 #include <QRegExp>
8 
9 #include "DawnVirCamera.h"
10 #include "Camera.h"
11 #include "NaifStatus.h"
12 #include "IString.h"
13 #include "IException.h"
14 #include "iTime.h"
15 #include "LineScanCameraDetectorMap.h"
16 #include "CameraFocalPlaneMap.h"
17 #include "LineScanCameraGroundMap.h"
18 #include "LineScanCameraSkyMap.h"
19 #include "NumericalApproximation.h"
20 #include "Kernels.h"
21 
22 #include "tnt/tnt_array2d_utils.h"
23 
24 // #define DUMP_INFO 1
25 
26 using namespace std;
27 namespace Isis {
28  // constructors
29  DawnVirCamera::DawnVirCamera(Cube &cube) : LineScanCamera(cube) {
30 
31  // cout << "Testing DawnVirCamera...\n";
32 
33  Pvl &lab = *cube.label();
34  PvlGroup &archive = lab.findGroup("Archive", Isis::Pvl::Traverse);
35  int procLevel = archive["ProcessingLevelId"];
36  m_is1BCalibrated = (procLevel > 2) ? true : false;
37 
38  // Get the start time from labels
39  PvlGroup &inst = lab.findGroup("Instrument", Isis::Pvl::Traverse);
40  QString channelId = inst["ChannelId"];
41 
42  QString instMode = inst["InstrumentModeId"];
43  m_slitMode = instMode[14].toAscii(); // "F" for full slit, Q for quarter slit
44 
45  // Check for presence of articulation kernel
46  bool hasArtCK = hasArticulationKernel(lab);
47 
48  // Set proper end frame
49  int virFrame(0);
50  if (channelId == "VIS") {
51  // Frame DAWN_VIR_VIS : DAWN_VIR_VIS_ZERO
52  virFrame = (hasArtCK) ? -203211 : -203221;
53  }
54  else { // (channelId == "IR)
55  // Frame DAWN_VIR_IR : DAWN_VIR_IR_ZERO
56  virFrame = (hasArtCK) ? -203213 : -203223;
57  }
58 
59  instrumentRotation()->SetFrame(virFrame);
60 
61  // We do not want to downsize the cache
62  instrumentRotation()->MinimizeCache(SpiceRotation::No);
63 
64  // Set up the camera info from ik/iak kernels
65  SetFocalLength();
66  SetPixelPitch();
67 
68  // Get other info from labels
69  PvlKeyword &frameParam = inst["FrameParameter"];
70  m_exposureTime = toDouble(frameParam[0]);
71  m_summing = toDouble(frameParam[1]);
72  m_scanRate = toDouble(frameParam[2]);
73 
74  // Setup detector map
75  // Get the line scan rates/times
76  readHouseKeeping(lab.fileName(), m_scanRate);
77  new VariableLineScanCameraDetectorMap(this, m_lineRates);
78  DetectorMap()->SetDetectorSampleSumming(m_summing);
79 
80  // Setup focal plane map
81  new CameraFocalPlaneMap(this, naifIkCode());
82 
83  // Retrieve boresight location from instrument kernel (IK) (addendum?)
84  QString ikernKey = "INS" + toString(naifIkCode()) + "_BORESIGHT_SAMPLE";
85  double sampleBoreSight = getDouble(ikernKey);
86 
87  ikernKey = "INS" + toString(naifIkCode()) + "_BORESIGHT_LINE";
88  double lineBoreSight = getDouble(ikernKey);
89 
90  FocalPlaneMap()->SetDetectorOrigin(sampleBoreSight, lineBoreSight);
91 
92  // Setup distortion map
93  new CameraDistortionMap(this);
94 
95  // Setup the ground and sky map
96  new LineScanCameraGroundMap(this);
97  new LineScanCameraSkyMap(this);
98 
99  // Set initial start time always (label start time is inaccurate)
100  setTime(iTime(startTime())); // Isis3nightly
101 // SetEphemerisTime(startTime()); // Isis3.2.1
102 
103  // Now check to determine if we have a cache already. If we have a cache
104  // table, we are beyond spiceinit and have already computed the proper
105  // point table from the housekeeping data or articulation kernel.
106  if (!instrumentRotation()->IsCached() && !hasArtCK) {
107 
108  // Create new table here prior to creating normal caches
109  Table quats = getPointingTable(channelId, virFrame);
110 
111  // Create all system tables - all kernels closed after this
112  LoadCache();
113  instrumentRotation()->LoadCache(quats);
114  }
115  else {
116  LoadCache();
117  }
118 
119 #if defined(DUMP_INFO)
120  Table cache = instrumentRotation()->Cache("Loaded");
121  cout << "Total Records: " << cache.Records() << "\n";
122 
123  for (int i = 0 ; i < cache.Records() ; i++) {
124  TableRecord rec = cache[i];
125  string separator("");
126  for (int f = 0 ; f < rec.Fields() ; f++) {
127  cout << separator << (double) rec[f];
128  separator = ", ";
129  }
130  cout << "\n";
131  }
132 #endif
133  }
134 
136  int DawnVirCamera::CkFrameId() const { return (-203000); }
138  int DawnVirCamera::CkReferenceId() const { return (1); }
140  int DawnVirCamera::SpkReferenceId() const { return (1); }
141 
142 
156  QString DawnVirCamera::scrub(const QString &text) const {
157  QString ostr;
158  for (int i = 0 ; i < text.size() ; i++) {
159  if ((text[i] > ' ') && (text[i] <= 'z')) ostr += text[i];
160  }
161  return (ostr);
162  }
163 
165  int DawnVirCamera::pixelSumming() const {
166  return (m_summing);
167  }
168 
170  double DawnVirCamera::exposureTime() const {
171  return (m_exposureTime);
172  }
173 
175  double DawnVirCamera::scanLineTime() const {
176  return (m_scanRate);
177  }
178 
180  double DawnVirCamera::lineStartTime(const double midExpTime) const {
181  return (midExpTime-(exposureTime()/2.0));
182  }
183 
185  double DawnVirCamera::lineEndTime(const double midExpTime) const {
186  return (midExpTime+(exposureTime()/2.0));
187  }
188 
190  double DawnVirCamera::startTime() const {
191  return (lineStartTime(m_mirrorData[0].m_scanLineEt));
192  }
193 
194 
196  double DawnVirCamera::endTime() const {
197  return (lineEndTime(m_mirrorData[hkLineCount()-1].m_scanLineEt));
198  }
199 
200 
202  int DawnVirCamera::hkLineCount() const {
203  return (m_mirrorData.size());
204  }
205 
217  void DawnVirCamera::readHouseKeeping(const QString &filename,
218  double lineRate) {
219  // Open the ISIS table object
220  Table hktable("VIRHouseKeeping", filename);
221 
222  m_lineRates.clear();
223  int lineno(1);
224  NumericalApproximation angFit;
225  for (int i = 0; i < hktable.Records(); i++) {
226  TableRecord &trec = hktable[i];
227  QString scet = scrub(trec["ScetTimeClock"]);
228  QString shutterMode = scrub(trec["ShutterStatus"]);
229 
230  // Compute the optical mirror angle
231  double mirrorSin = trec["MirrorSin"];
232  double mirrorCos = trec["MirrorCos"];
233  double scanElecDeg = atan(mirrorSin/mirrorCos) * dpr_c();
234  double optAng = ((scanElecDeg - 3.7996979) * 0.25/0.257812);
235  optAng /= 1000.0;
236 
237 
238  ScanMirrorInfo smInfo;
239  double lineMidTime;
240  // scs2e_c(naifSpkCode(), scet.c_str(), &lineMidTime);
241  lineMidTime = getClockTime(scet, naifSpkCode()).Et();
242  bool isDark = shutterMode.toLower() == "closed";
243 
244  // Add fit data for all open angles
245  if ( ! isDark ) { angFit.AddData(lineno, optAng); }
246 
247 #if defined(DUMP_INFO)
248  cout << "Line(" << ((isDark) ? "C): " : "O): ") << i
249  << ", OptAng(D): " << setprecision(12) << optAng * dpr_c()
250  << ", MidExpTime(ET): " << lineMidTime
251  << "\n";
252 #endif
253 
254  // Store line,
255  smInfo.m_lineNum = lineno;
256  smInfo.m_scanLineEt = lineMidTime;
257  smInfo.m_mirrorSin = mirrorSin;
258  smInfo.m_mirrorCos = mirrorCos;
259  smInfo.m_opticalAngle = optAng;
260  smInfo.m_isDarkCurrent = isDark;
261 
262  if ((!m_is1BCalibrated) || (!(m_is1BCalibrated && isDark))) {
263  m_lineRates.push_back(LineRateChange(lineno,
264  lineStartTime(lineMidTime),
265  lineRate));
266  m_mirrorData.push_back(smInfo);
267  lineno++;
268  }
269  }
270 
271  // Adjust the last time
272  LineRateChange lastR = m_lineRates.back();
273  m_lineRates.back() = LineRateChange(lastR.GetStartLine(),
274  lastR.GetStartEt(),
275  exposureTime());
276 
277  // Run through replacing all closed optical angles with fitted data.
278  // These are mostly first/last lines so must set proper extrapolation
279  // option.
280  for (unsigned int a = 0 ; a < m_mirrorData.size() ; a++) {
281  if (m_mirrorData[a].m_isDarkCurrent) {
282  m_mirrorData[a].m_opticalAngle = angFit.Evaluate(a+1,
283  NumericalApproximation::NearestEndpoint);
284  }
285  }
286 
287  // Gut check on housekeeping contents and cube lines
288  if ((int) m_lineRates.size() != Lines()) {
289  ostringstream mess;
290  mess << "Number housekeeping lines determined (" << m_lineRates.size()
291  << ") is not equal to image lines(" << Lines() << ")";
292  throw IException(IException::Programmer, mess.str(), _FILEINFO_);
293  }
294  }
295 
305  Table DawnVirCamera::getPointingTable(const QString &virChannel,
306  const int zeroFrame) {
307 
308  // Create Spice Pointing table
309  TableField q0("J2000Q0", TableField::Double);
310  TableField q1("J2000Q1", TableField::Double);
311  TableField q2("J2000Q2", TableField::Double);
312  TableField q3("J2000Q3", TableField::Double);
313  TableField av1("AV1", TableField::Double);
314  TableField av2("AV2", TableField::Double);
315  TableField av3("AV3", TableField::Double);
316  TableField t("ET", TableField::Double);
317 
318  TableRecord record;
319  record += q0;
320  record += q1;
321  record += q2;
322  record += q3;
323  record += av1;
324  record += av2;
325  record += av3;
326  record += t;
327 
328  // Get pointing table
329  Table quats("SpiceRotation", record);
330  int nfields = record.Fields();
331 
332  QString virId = "DAWN_VIR_" + virChannel;
333  QString virZero = virId + "_ZERO";
334 
335  // Allocate output arrays
336  int nvals = nfields - 1;
337  int nlines = m_lineRates.size();
338 
339  SpiceDouble eulang[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
340  SpiceDouble xform[6][6], xform2[6][6];
341  SpiceDouble m[3][3];
342  SpiceDouble q_av[7], *av(&q_av[4]);
343 
344  for (int i = 0 ; i < nlines ; i++) {
345  int index = min(i, nlines-1);
346  double etTime = m_mirrorData[index].m_scanLineEt; // mid exposure ET
347  double optAng = m_mirrorData[index].m_opticalAngle;
348  try {
349  // J2000 -> DAWN_VIR_{channel}_ZERO
350  SMatrix state = getStateRotation("J2000", virZero, etTime);
351 
352  // Set rotation of optical scan mirror (in radians)
353  eulang[1] = -optAng;
354  eul2xf_c(eulang, 1, 2, 3, xform);
355  mxmg_c(xform, &state[0][0], 6, 6, 6, xform2);
356 
357  // Transform to output format
358  xf2rav_c(xform2, m, av); // Transfers AV to output q_av via pointer
359  m2q_c(m, q_av); // Transfers quaternion
360 
361  // Now populate the table record with the line pointing
362  for (int k = 0 ; k < nvals ; k++) {
363  record[k] = q_av[k];
364  }
365 
366  // Add time to record; record to table
367  record[nvals] = etTime;
368  quats += record;
369  }
370  catch (IException &ie) {
371  ostringstream mess;
372  mess << "Failed to get point state for line " << i+1;
373  throw IException(ie, IException::User, mess.str(), _FILEINFO_);
374  }
375  }
376 
377  // Add some necessary keywords
378  quats.Label() += PvlKeyword("CkTableStartTime", toString(startTime()));
379  quats.Label() += PvlKeyword("CkTableEndTime", toString(endTime()));
380  quats.Label() += PvlKeyword("CkTableOriginalSize", toString(quats.Records()));
381 
382  // Create the time dependant frames keyword
383  int virZeroId = getInteger("FRAME_" + virZero);
384  PvlKeyword tdf("TimeDependentFrames", toString(virZeroId)); // DAWN_VIR_{ID}_ZERO
385  tdf.addValue("-203200"); // DAWN_VIR
386  tdf.addValue("-203000"); // DAWN_SPACECRAFT
387  tdf.addValue("1"); // J2000
388  quats.Label() += tdf;
389 
390  // Create constant rotation frames
391  PvlKeyword cf("ConstantFrames", toString(virZeroId));
392  cf.addValue(toString(virZeroId));
393  quats.Label() += cf;
394 
395  SpiceDouble identity[3][3];
396  ident_c(identity);
397 
398  // Store DAWN_VIR_{ID}_ZERO -> DAWN_VIR_{ID}_ZERO identity rotation
399  PvlKeyword crot("ConstantRotation");
400  for (int i = 0 ; i < 3 ; i++) {
401  for (int j = 0 ; j < 3 ; j++) {
402  crot.addValue(toString(identity[i][j]));
403  }
404  }
405 
406  quats.Label() += crot;
407 
408  return (quats);
409  }
410 
427  DawnVirCamera::SMatrix DawnVirCamera::getStateRotation(const QString &frame1,
428  const QString &frame2,
429  const double &etTime)
430  const {
431  SMatrix state(6,6);
432  NaifStatus::CheckErrors();
433  try {
434  // Get pointing w/AVs
435  sxform_c(frame1.toAscii().data(), frame2.toAscii().data(), etTime,
436  (SpiceDouble (*)[6]) state[0]);
437  NaifStatus::CheckErrors();
438  }
439  catch (IException &) {
440  try {
441  SMatrix rot(3,3);
442  pxform_c(frame1.toAscii().data(), frame2.toAscii().data(), etTime,
443  (SpiceDouble (*)[3]) rot[0]);
444  NaifStatus::CheckErrors();
445  SpiceDouble av[3] = {0.0, 0.0, 0.0 };
446  rav2xf_c((SpiceDouble (*)[3]) rot[0], av,
447  (SpiceDouble (*)[6]) state[0]);
448  }
449  catch (IException &ie2) {
450  ostringstream mess;
451  mess << "Could not get state rotation for Frame1 (" << frame1
452  << ") to Frame2 (" << frame2 << ") at time " << etTime;
453  throw IException(ie2, IException::User, mess.str(), _FILEINFO_);
454  }
455  }
456  return (state);
457  }
458 
473  bool DawnVirCamera::hasArticulationKernel(Pvl &label) const {
474  Kernels kerns(label);
475  QStringList cks = kerns.getKernelList("CK");
476  QRegExp virCk("*dawn_vir_?????????_?.bc");
477  virCk.setPatternSyntax(QRegExp::Wildcard);
478  for (int i = 0 ; i < cks.size() ; i++) {
479  if ( virCk.exactMatch(cks[i]) ) return (true);
480  }
481  return (false);
482  }
483 
484 }
485 
487 extern "C" Isis::Camera *DawnVirCameraPlugin(Isis::Cube &cube) {
488  return new Isis::DawnVirCamera(cube);
489 }