USGS

Isis 3.0 Application Source Code Reference

Home

hicalbeta.cpp
Go to the documentation of this file.
1 // $Id: hicalbeta.cpp,v 1.14 2009/09/15 21:56:44 kbecker Exp $
2 #include "Isis.h"
3 
4 #include <cstdio>
5 #include <string>
6 #include <vector>
7 #include <algorithm>
8 #include <sstream>
9 #include <iostream>
10 
11 #include "FileName.h"
12 #include "ProcessByLine.h"
13 #include "UserInterface.h"
14 #include "Pvl.h"
15 #include "PvlGroup.h"
16 #include "IString.h"
17 #include "HiCalConf.h"
18 #include "CollectorMap.h"
19 #include "HiCalTypes.h"
20 #include "HiCalUtil.h"
21 #include "HiCalData.h"
22 #include "Statistics.h"
23 #include "SplineFill.h" // SpineFillComp.h
24 #include "LowPassFilter.h" // LowPassFilterComp.h
25 #include "ZeroBufferSmooth.h" // DriftBuffer.h (Zf)
26 #include "ZeroBufferFit.h" // DriftCorrect.h (Zd)
27 #include "ZeroReverse.h" // OffsetCorrect.h (Zz)
28 #include "ZeroDark.h" // DarkSubtractComp.h (Zb)
29 #include "GainLineDrift.h" // GainVLineComp.h (Zg)
30 #include "GainNonLinearity.h" // Non-linear gain (new)
31 #include "GainChannelNormalize.h" // ZggModule.h (Zgg)
32 #include "GainFlatField.h" // FlatFieldComp.h (Za)
33 #include "GainTemperature.h" // TempGainCorrect.h (Zt)
34 #include "GainUnitConversion.h" // ZiofModule.h (Ziof)
35 
36 
37 using namespace Isis;
38 using namespace std;
39 
40 //!< Define the matrix container for systematic processing
41 typedef CollectorMap<IString, HiVector, NoCaseStringCompare> MatrixList;
42 
43 // Calibration parameter container
45 
46 
47 /**
48  * @brief Apply calibration to each HiRISE image line
49  *
50  * This function applies the calbration equation to each input image line. It
51  * gets matrices and constants from the \b calVars container that is established
52  * in the main with some user input via the configuration (CONF) parameter.
53  *
54  * @param in Input raw image line buffer
55  * @param out Output calibrated image line buffer
56  */
57 void calibrate(Buffer &in, Buffer &out) {
58  const HiVector &ZBF = calVars->get("ZeroBufferFit");
59  const HiVector &ZRev = calVars->get("ZeroReverse");
60  const HiVector &ZD = calVars->get("ZeroDark");
61  const HiVector &GLD = calVars->get("GainLineDrift");
62  const HiVector &GCN = calVars->get("GainChannelNormalize");
63  const double GNL = calVars->get("GainNonLinearity")[0];
64  const HiVector &GFF = calVars->get("GainFlatField");
65  const HiVector &GT = calVars->get("GainTemperature");
66  const double GUC = calVars->get("GainUnitConversion")[0];
67 
68  // Set current line (index)
69  int line(in.Line()-1);
70  if (calVars->exists("LastGoodLine")) {
71  int lastline = ((int) (calVars->get("LastGoodLine"))[0]) - 1;
72  if ( line > lastline ) { line = lastline; }
73  }
74 
75  // Apply correction to point of non-linearity accumulating average
76  vector<double> data;
77  for (int i = 0 ; i < in.size() ; i++) {
78  if (IsSpecial(in[i])) {
79  out[i] = in[i];
80  }
81  else {
82  double hdn;
83  hdn = (in[i] - ZBF[line] - ZRev[i] - ZD[i]); // Drift, Reverse, Dark
84  hdn = hdn / GLD[line]; // GainLineDrift
85  data.push_back(hdn); // Accumulate average for non-linearity
86  out[i] = hdn;
87  }
88  }
89 
90  // Second loop require to apply non-linearity gain correction from average
91  if ( data.size() > 0 ) { // No valid pixels means just that - done
92  // See HiCalUtil.h for the function that returns this stat.
93  double NLGain = 1.0 - (GNL * GainLineStat(data));
94  for (int i = 0 ; i < out.size() ; i++) {
95  if (!IsSpecial(out[i])) {
96  double hdn = out[i];
97  hdn = hdn * GCN[i] * NLGain * GFF[i] * GT[i]; // Gain, Non-linearity
98  // gain, FlatField,
99  // TempGain
100  out[i] = hdn / GUC; // I/F or DN or DN/US
101  }
102  }
103  }
104  return;
105 }
106 
107 
108 void IsisMain(){
109 
110  const QString hical_program = "hicalbeta";
111  const QString hical_version = "5.0";
112  const QString hical_revision = "$Revision: 1.15 $";
113  const QString hical_runtime = Application::DateTime();
114 
115  UserInterface &ui = Application::GetUserInterface();
116 
117  QString procStep("prepping phase");
118  try {
119 // The output from the last processing is the input into subsequent processing
120  ProcessByLine p;
121 
122  Cube *hifrom = p.SetInputCube("FROM");
123  int nsamps = hifrom->sampleCount();
124  int nlines = hifrom->lineCount();
125 
126 // Initialize the configuration file
127  QString conf(ui.GetAsString("CONF"));
128  HiCalConf hiconf(*(hifrom->label()), conf);
129  DbProfile hiprof = hiconf.getMatrixProfile();
130 
131 // Check for label propagation and set the output cube
132  Cube *ocube = p.SetOutputCube("TO");
133  if ( !IsTrueValue(hiprof,"PropagateTables", "TRUE") ) {
134  RemoveHiBlobs(*(ocube->label()));
135  }
136 
137 // Set specified profile if entered by user
138  if (ui.WasEntered("PROFILE")) {
139  hiconf.selectProfile(ui.GetAsString("PROFILE"));
140  }
141 
142 
143 // Add OPATH parameter to profiles
144  if (ui.WasEntered("OPATH")) {
145  hiconf.add("OPATH",ui.GetAsString("OPATH"));
146  }
147  else {
148  // Set default to output directory
149  hiconf.add("OPATH", FileName(ocube->fileName()).path());
150  }
151 
152 // Do I/F output DN conversions
153  QString units = ui.GetString("UNITS");
154 
155  // Allocate the calibration list
156  calVars = new MatrixList;
157 
158 // Set up access to HiRISE ancillary data (tables, blobs) here. Note it they
159 // are gone, this will error out. See PropagateTables in conf file.
160  HiCalData caldata(*hifrom);
161 
162 ////////////////////////////////////////////////////////////////////////////
163 // Drift Correction (Zf) using buffer pixels
164 // Extracts specified regions of the calibration buffer pixels and runs
165 // series of lowpass filters. Apply spline fit if any missing data
166 // remains. Config file contains parameters for this operation.
167  procStep = "ZeroBufferSmooth module";
168  hiconf.selectProfile("ZeroBufferSmooth");
169  hiprof = hiconf.getMatrixProfile();
170  HiHistory ZbsHist;
171  ZbsHist.add("Profile["+ hiprof.Name()+"]");
172  if ( !SkipModule(hiprof) ) {
173  ZeroBufferSmooth zbs(caldata, hiconf);
174  calVars->add("ZeroBufferSmooth", zbs.ref());
175  ZbsHist = zbs.History();
176  if ( hiprof.exists("DumpModuleFile") ) {
177  zbs.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
178  }
179  }
180  else {
181  // NOT RECOMMENDED! This is required for the next step!
182  // SURELY must be skipped with ZeroBufferSmooth step as well!
183  calVars->add("ZeroBufferSmooth", HiVector(nlines, 0.0));
184  ZbsHist.add("Debug::SkipModule invoked!");
185  }
186 
187 /////////////////////////////////////////////////////////////////////
188 // ZeroBufferFit
189 // Compute second level of drift correction. The high level noise
190 // is removed from a modeled non-linear fit.
191 //
192  procStep = "ZeroBufferFit module";
193  HiHistory ZbfHist;
194  hiconf.selectProfile("ZeroBufferFit");
195  hiprof = hiconf.getMatrixProfile();
196  ZbfHist.add("Profile["+ hiprof.Name()+"]");
197  if (!SkipModule(hiprof) ) {
198  ZeroBufferFit zbf(hiconf);
199 
200  calVars->add(hiconf.getProfileName(),
201  zbf.Normalize(zbf.Solve(calVars->get("ZeroBufferSmooth"))));
202  ZbfHist = zbf.History();
203  if ( hiprof.exists("DumpModuleFile") ) {
204  zbf.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
205  }
206  }
207  else {
208  calVars->add(hiconf.getProfileName(), HiVector(nlines, 0.0));
209  ZbfHist.add("Debug::SkipModule invoked!");
210  }
211 
212 
213  ////////////////////////////////////////////////////////////////////
214  // ZeroReverse
215  procStep = "ZeroReverse module";
216  hiconf.selectProfile("ZeroReverse");
217  hiprof = hiconf.getMatrixProfile();
218  HiHistory ZrHist;
219  ZrHist.add("Profile["+ hiprof.Name()+"]");
220  if ( !SkipModule(hiprof) ) {
221  ZeroReverse zr(caldata, hiconf);
222  calVars->add(hiconf.getProfileName(), zr.ref());
223  ZrHist = zr.History();
224  if ( hiprof.exists("DumpModuleFile") ) {
225  zr.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
226  }
227  }
228  else {
229  calVars->add(hiconf.getProfileName(), HiVector(nsamps, 0.0));
230  ZrHist.add("Debug::SkipModule invoked!");
231  }
232 
233 /////////////////////////////////////////////////////////////////
234 // ZeroDark removes dark current
235 //
236  procStep = "ZeroDark module";
237  hiconf.selectProfile("ZeroDark");
238  hiprof = hiconf.getMatrixProfile();
239  HiHistory ZdHist;
240  ZdHist.add("Profile["+ hiprof.Name()+"]");
241  if ( !SkipModule(hiprof) ) {
242  ZeroDark zd(hiconf);
243  calVars->add(hiconf.getProfileName(), zd.ref());
244  ZdHist = zd.History();
245  if ( hiprof.exists("DumpModuleFile") ) {
246  zd.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
247  }
248  }
249  else {
250  calVars->add(hiconf.getProfileName(), HiVector(nsamps, 0.0));
251  ZdHist.add("Debug::SkipModule invoked!");
252  }
253 
254 ////////////////////////////////////////////////////////////////////
255 // GainLineDrift correct for gain-based drift
256 //
257  procStep = "GainLineDrift module";
258  hiconf.selectProfile("GainLineDrift");
259  hiprof = hiconf.getMatrixProfile();
260  HiHistory GldHist;
261  GldHist.add("Profile["+ hiprof.Name()+"]");
262  if ( !SkipModule(hiprof) ) {
263  GainLineDrift gld(hiconf);
264  calVars->add(hiconf.getProfileName(), gld.ref());
265  GldHist = gld.History();
266  if ( hiprof.exists("DumpModuleFile") ) {
267  gld.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
268  }
269  }
270  else {
271  calVars->add(hiconf.getProfileName(), HiVector(nlines, 1.0));
272  GldHist.add("Debug::SkipModule invoked!");
273  }
274 
275 ////////////////////////////////////////////////////////////////////
276 // GainNonLinearity Correct for non-linear gain
277  procStep = "GainNonLinearity module";
278  hiconf.selectProfile("GainNonLinearity");
279  hiprof = hiconf.getMatrixProfile();
280  HiHistory GnlHist;
281  GnlHist.add("Profile["+ hiprof.Name()+"]");
282  if ( !SkipModule(hiprof) ) {
283  GainNonLinearity gnl(hiconf);
284  calVars->add(hiconf.getProfileName(), gnl.ref());
285  GnlHist = gnl.History();
286  if ( hiprof.exists("DumpModuleFile") ) {
287  gnl.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
288  }
289  }
290  else {
291  calVars->add(hiconf.getProfileName(), HiVector(1, 0.0));
292  GnlHist.add("Debug::SkipModule invoked!");
293  }
294 
295 ////////////////////////////////////////////////////////////////////
296 // GainChannelNormalize Correct for sample gain with the G matrix
297  procStep = "GainChannelNormalize module";
298  hiconf.selectProfile("GainChannelNormalize");
299  hiprof = hiconf.getMatrixProfile();
300  HiHistory GcnHist;
301  GcnHist.add("Profile["+ hiprof.Name()+"]");
302  if ( !SkipModule(hiprof) ) {
303  GainChannelNormalize gcn(hiconf);
304  calVars->add(hiconf.getProfileName(), gcn.ref());
305  GcnHist = gcn.History();
306  if ( hiprof.exists("DumpModuleFile") ) {
307  gcn.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
308  }
309  }
310  else {
311  calVars->add(hiconf.getProfileName(), HiVector(nsamps, 1.0));
312  GcnHist.add("Debug::SkipModule invoked!");
313  }
314 
315 ////////////////////////////////////////////////////////////////////
316 // GainFlatField Flat field correction with A matrix
317  procStep = "GainFlatField module";
318  hiconf.selectProfile("GainFlatField");
319  hiprof = hiconf.getMatrixProfile();
320  HiHistory GffHist;
321  GffHist.add("Profile["+ hiprof.Name()+"]");
322  if ( !SkipModule(hiprof) ) {
323  GainFlatField gff(hiconf);
324  calVars->add(hiconf.getProfileName(), gff.ref());
325  GffHist = gff.History();
326  if ( hiprof.exists("DumpModuleFile") ) {
327  gff.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
328  }
329  }
330  else {
331  calVars->add(hiconf.getProfileName(), HiVector(nsamps, 1.0));
332  GffHist.add("Debug::SkipModule invoked!");
333  }
334 
335 ////////////////////////////////////////////////////////////////////
336 // GainTemperature - Temperature-dependant gain correction
337  procStep = "GainTemperature module";
338  hiconf.selectProfile("GainTemperature");
339  hiprof = hiconf.getMatrixProfile();
340  HiHistory GtHist;
341  GtHist.add("Profile["+ hiprof.Name()+"]");
342  if ( !SkipModule(hiprof) ) {
343  GainTemperature gt(hiconf);
344  calVars->add(hiconf.getProfileName(), gt.ref());
345  GtHist = gt.History();
346  if ( hiprof.exists("DumpModuleFile") ) {
347  gt.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
348  }
349  }
350  else {
351  calVars->add(hiconf.getProfileName(), HiVector(nsamps, 1.0));
352  GtHist.add("Debug::SkipModule invoked!");
353  }
354 
355 ////////////////////////////////////////////////////////////////////
356 // GainUnitConversion converts to requested units
357 //
358  procStep = "GainUnitConversion module";
359  hiconf.selectProfile("GainUnitConversion");
360  hiprof = hiconf.getMatrixProfile();
361  HiHistory GucHist;
362  GucHist.add("Profile["+ hiprof.Name()+"]");
363  if ( !SkipModule(hiprof) ) {
364  GainUnitConversion guc(hiconf, units);
365  calVars->add(hiconf.getProfileName(), guc.ref());
366  GucHist = guc.History();
367  if ( hiprof.exists("DumpModuleFile") ) {
368  guc.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
369  }
370  }
371  else {
372  calVars->add(hiconf.getProfileName(), HiVector(1,1.0));
373  GucHist.add("Debug::SkipModule invoked!");
374  GucHist.add("Units[Unknown]");
375  }
376 
377  // Reset the profile selection to default
378  hiconf.selectProfile();
379 
380 //----------------------------------------------------------------------
381 //
382 /////////////////////////////////////////////////////////////////////////
383 // Call the processing function
384  procStep = "calibration phase";
385  p.StartProcess(calibrate);
386 
387  // Get the default profile for logging purposes
388  hiprof = hiconf.getMatrixProfile();
389  const QString conf_file = hiconf.filepath(conf);
390 
391  // Quitely dumps parameter history to alternative format file. This
392  // is completely controlled by the configuration file
393  if ( hiprof.exists("DumpHistoryFile") ) {
394  procStep = "logging/reporting phase";
395  FileName hdump(hiconf.getMatrixSource("DumpHistoryFile",hiprof));
396  QString hdumpFile = hdump.expanded();
397  ofstream ofile(hdumpFile.toLatin1().data(), ios::out);
398  if (!ofile) {
399  QString mess = "Unable to open/create history dump file " +
400  hdump.expanded();
401  IException(IException::User, mess, _FILEINFO_).print();
402  }
403  else {
404  ofile << "Program: " << hical_program << endl;
405  ofile << "RunTime: " << hical_runtime << endl;
406  ofile << "Version: " << hical_version << endl;
407  ofile << "Revision: " << hical_revision << endl << endl;
408 
409  ofile << "FROM: " << hifrom->fileName() << endl;
410  ofile << "TO: " << ocube->fileName() << endl;
411  ofile << "CONF: " << conf_file << endl << endl;
412 
413  ofile << "/* " << hical_program << " application equation */\n"
414  << "/* hdn = (idn - ZeroBufferFit(ZeroBufferSmooth) - ZeroReverse - ZeroDark) */\n"
415  << "/* odn = hdn / GainLineDrift * GainNonLinearity * GainChannelNormalize */\n"
416  << "/* * GainFlatField * GainTemperature / GainUnitConversion */\n\n";
417 
418  ofile << "****** PARAMETER GENERATION HISTORY *******" << endl;
419  ofile << "\nZeroBufferSmooth = " << ZbsHist << endl;
420  ofile << "\nZeroBufferFit = " << ZbfHist << endl;
421  ofile << "\nZeroReverse = " << ZrHist << endl;
422  ofile << "\nZeroDark = " << ZdHist << endl;
423  ofile << "\nGainLineDrift = " << GldHist << endl;
424  ofile << "\nGainNonLinearity = " << GnlHist << endl;
425  ofile << "\nGainChannelNormalize = " << GcnHist << endl;
426  ofile << "\nGainFlatField = " << GffHist << endl;
427  ofile << "\nGainTemperature = " << GtHist << endl;
428  ofile << "\nGainUnitConversion = " << GucHist << endl;
429 
430  ofile.close();
431  }
432  }
433 
434 // Ensure the RadiometricCalibration group is out there
435  const QString rcalGroup("RadiometricCalibration");
436  if (!ocube->hasGroup(rcalGroup)) {
437  PvlGroup temp(rcalGroup);
438  ocube->putGroup(temp);
439  }
440 
441  PvlGroup &rcal = ocube->group(rcalGroup);
442  rcal += PvlKeyword("Program", hical_program);
443  rcal += PvlKeyword("RunTime", hical_runtime);
444  rcal += PvlKeyword("Version",hical_version);
445  rcal += PvlKeyword("Revision",hical_revision);
446 
447  PvlKeyword key("Conf", conf_file);
448  key.addCommentWrapped("/* " + hical_program + " application equation */");
449  key.addComment("/* hdn = idn - ZeroBufferFit(ZeroBufferSmooth) */");
450  key.addComment("/* - ZeroReverse - ZeroDark */");
451  key.addComment("/* odn = hdn / GainLineDrift * GainNonLinearity */");
452  key.addComment("/* * GainChannelNormalize * GainFlatField */");
453  key.addComment("/* * GainTemperature / GainUnitConversion */");
454  rcal += key;
455 
456  // Record parameter generation history. Controllable in configuration
457  // file. Note this is optional because of a BUG!! in the ISIS label
458  // writer as this application was initially developed
459  if ( IsEqual(ConfKey(hiprof,"LogParameterHistory",QString("TRUE")),"TRUE")) {
460  rcal += ZbsHist.makekey("ZeroBufferSmooth");
461  rcal += ZbfHist.makekey("ZeroBufferFit");
462  rcal += ZrHist.makekey("ZeroReverse");
463  rcal += ZdHist.makekey("ZeroDark");
464  rcal += GldHist.makekey("GainLineDrift");
465  rcal += GnlHist.makekey("GainNonLinearity");
466  rcal += GcnHist.makekey("GainChannelNormalize");
467  rcal += GffHist.makekey("GainFlatField");
468  rcal += GtHist.makekey("GainTemperature");
469  rcal += GucHist.makekey("GainUnitConversion");
470  }
471 
472  p.EndProcess();
473  }
474  catch (IException &ie) {
475  delete calVars;
476  calVars = 0;
477  QString mess = "Failed in " + procStep;
478  throw IException(ie, IException::User, mess, _FILEINFO_);
479  }
480 
481 // Clean up parameters
482  delete calVars;
483  calVars = 0;
484 }
485 
HiVector Solve(const HiVector &d)
Compute non-linear fit to (typically) ZeroBufferSmooth module.
GainNonLinearity Module - Applies non-linear, line-dependant gain.
void Dump(const QString &fname) const
Dumps the component to a specified file.
Definition: Module.h:146
void IsisMain()
Definition: apollo2isis.cpp:42
void RemoveHiBlobs(Pvl &label)
Deletes HiRISE specific BLOBS from cube file.
Definition: HiCalUtil.h:463
LineManager * in
Definition: crop.cpp:24
bool SkipModule(const DbProfile &prof)
Checks profile flag to skip the current Module.
Definition: HiCalUtil.h:309
Cube * ocube
Definition: demprep.cpp:20
TNT::Array1D< double > HiVector
1-D Buffer
Definition: HiCalTypes.h:38
CollectorMap< IString, HiVector, NoCaseStringCompare > MatrixList
&lt; Define the matrix container for systematic processing
Definition: hical.cpp:41
T ConfKey(const DbProfile &conf, const QString &keyname, const T &defval, int index=0)
Find a keyword in a profile using default for non-existant keywords.
Definition: HiCalUtil.h:219
const HiHistory & History() const
Return recorded history of events.
Definition: Module.h:131
int line
Definition: blend.cpp:120
Computes non-linear lsq fit of HiRISE Drift (Zd module)
Definition: ZeroBufferFit.h:56
double GainLineStat(std::vector< double > &data)
Return the statistics for a vector of data.
Definition: HiCalUtil.h:492
Computes a complex dark subtraction component (ZeroDark module)
Definition: ZeroDark.h:62
bool IsTrueValue(const DbProfile &prof, const QString &key, const QString &value="TRUE")
Determines if the keyword value is the expected value.
Definition: HiCalUtil.h:290
GaingTemperature Module - Applies temperature-dependant gain correction (column)
MatrixList * calVars
Definition: hical.cpp:44
HiVector Normalize(const HiVector &v)
Compute normalized solution vector from result.
bool IsEqual(const QString &v1, const QString &v2="TRUE")
Shortened string equality test.
Definition: HiCalUtil.h:272
Processes Reverse Clock calibration data (ZeroReverse Module)
Definition: ZeroReverse.h:60
GainFlatField Module - Computes flat field correction for sample.
Definition: GainFlatField.h:57
Processes Buffer calibration data (ZeroBufferSmooth Module)
Computes units parameters for HiRISE data calibration (Ziof Module)
const HiVector & ref() const
Return data via a const reference.
Definition: Module.h:126
void calibrate(vector< Buffer * > &in, vector< Buffer * > &out)
This applies the calculated calibration coefficients to the file.
Definition: vimscal.cpp:240
Computes a gain correction for each sample GainChannelNormalize.
Cube * out
Definition: lrowac2pds.cpp:50
Computes a gain correction for each line (Zg Module)
Definition: GainLineDrift.h:55
void add(const QString &event)
Definition: HiCalTypes.h:71
Container for HiRISE calibration data.
Definition: HiCalData.h:50
PvlKeyword makekey(const QString &name="History") const
Definition: HiCalTypes.h:83