USGS

Isis 3.0 Application Source Code Reference

Home

hicalbeta.cpp

Go to the documentation of this file.
00001 //  $Id: hicalbeta.cpp,v 1.14 2009/09/15 21:56:44 kbecker Exp $
00002 #include "Isis.h"
00003 
00004 #include <cstdio>
00005 #include <string>
00006 #include <vector>
00007 #include <algorithm>
00008 #include <sstream>
00009 #include <iostream>
00010 
00011 #include "FileName.h"
00012 #include "ProcessByLine.h"
00013 #include "UserInterface.h"
00014 #include "Pvl.h"
00015 #include "PvlGroup.h"
00016 #include "IString.h"
00017 #include "HiCalConf.h"
00018 #include "CollectorMap.h"
00019 #include "HiCalTypes.h"
00020 #include "HiCalUtil.h"
00021 #include "HiCalData.h"
00022 #include "Statistics.h"
00023 #include "SplineFill.h"              // SpineFillComp.h
00024 #include "LowPassFilter.h"           // LowPassFilterComp.h
00025 #include "ZeroBufferSmooth.h"        // DriftBuffer.h  (Zf)
00026 #include "ZeroBufferFit.h"           // DriftCorrect.h (Zd)
00027 #include "ZeroReverse.h"             // OffsetCorrect.h (Zz)
00028 #include "ZeroDark.h"                // DarkSubtractComp.h (Zb)
00029 #include "GainLineDrift.h"           // GainVLineComp.h (Zg)
00030 #include "GainNonLinearity.h"        // Non-linear gain (new)
00031 #include "GainChannelNormalize.h"    // ZggModule.h (Zgg)
00032 #include "GainFlatField.h"           // FlatFieldComp.h (Za)
00033 #include "GainTemperature.h"         // TempGainCorrect.h (Zt)
00034 #include "GainUnitConversion.h"      // ZiofModule.h (Ziof)
00035 
00036 
00037 using namespace Isis;
00038 using namespace std;
00039 
00040 //!< Define the matrix container for systematic processing
00041 typedef CollectorMap<IString, HiVector, NoCaseStringCompare> MatrixList;
00042 
00043 //  Calibration parameter container
00044 MatrixList *calVars = 0;
00045 
00046 
00047 /**
00048  * @brief Apply calibration to each HiRISE image line
00049  *
00050  * This function applies the calbration equation to each input image line.  It
00051  * gets matrices and constants from the \b calVars container that is established
00052  * in the main with some user input via the configuration (CONF) parameter.
00053  *
00054  * @param in  Input raw image line buffer
00055  * @param out Output calibrated image line buffer
00056  */
00057 void calibrate(Buffer &in, Buffer &out) {
00058   const HiVector &ZBF    = calVars->get("ZeroBufferFit");
00059   const HiVector &ZRev   = calVars->get("ZeroReverse");
00060   const HiVector &ZD     = calVars->get("ZeroDark");
00061   const HiVector &GLD    = calVars->get("GainLineDrift");
00062   const HiVector &GCN    = calVars->get("GainChannelNormalize");
00063   const double GNL       = calVars->get("GainNonLinearity")[0];
00064   const HiVector &GFF    = calVars->get("GainFlatField");
00065   const HiVector &GT     = calVars->get("GainTemperature");
00066   const double GUC       = calVars->get("GainUnitConversion")[0];
00067 
00068   //  Set current line (index)
00069   int line(in.Line()-1);
00070   if (calVars->exists("LastGoodLine")) {
00071     int lastline = ((int) (calVars->get("LastGoodLine"))[0]) - 1;
00072     if ( line > lastline ) { line = lastline; }
00073   }
00074 
00075   //  Apply correction to point of non-linearity accumulating average
00076   vector<double> data;
00077   for (int i = 0 ; i < in.size() ; i++) {
00078     if (IsSpecial(in[i])) {
00079       out[i] = in[i];
00080     }
00081     else {
00082       double hdn;
00083       hdn = (in[i] - ZBF[line] - ZRev[i] - ZD[i]); // Drift, Reverse, Dark
00084       hdn = hdn / GLD[line];  // GainLineDrift
00085       data.push_back(hdn);   // Accumulate average for non-linearity
00086       out[i] = hdn;
00087     }
00088   }
00089 
00090   // Second loop require to apply non-linearity gain correction from average
00091   if ( data.size() > 0 ) {  // No valid pixels means just that - done
00092     // See HiCalUtil.h for the function that returns this stat.
00093     double NLGain = 1.0 - (GNL * GainLineStat(data));
00094     for (int i = 0 ; i < out.size() ; i++) {
00095       if (!IsSpecial(out[i])) {
00096         double hdn = out[i];
00097         hdn = hdn * GCN[i] * NLGain * GFF[i] * GT[i];  // Gain, Non-linearity
00098                                                        // gain,  FlatField,
00099                                                        // TempGain
00100         out[i] = hdn / GUC;                   // I/F or DN or DN/US
00101       }
00102     }
00103   }
00104   return;
00105 }
00106 
00107 
00108 void IsisMain(){
00109 
00110   const QString hical_program = "hicalbeta";
00111   const QString hical_version = "5.0";
00112   const QString hical_revision = "$Revision: 1.15 $";
00113   const QString hical_runtime = Application::DateTime();
00114 
00115   UserInterface &ui = Application::GetUserInterface();
00116 
00117   QString procStep("prepping phase");
00118   try {
00119 //  The output from the last processing is the input into subsequent processing
00120     ProcessByLine p;
00121 
00122     Cube *hifrom = p.SetInputCube("FROM");
00123     int nsamps = hifrom->sampleCount();
00124     int nlines = hifrom->lineCount();
00125 
00126 //  Initialize the configuration file
00127     QString conf(ui.GetAsString("CONF"));
00128     HiCalConf hiconf(*(hifrom->label()), conf);
00129     DbProfile hiprof = hiconf.getMatrixProfile();
00130 
00131 // Check for label propagation and set the output cube
00132     Cube *ocube = p.SetOutputCube("TO");
00133     if ( !IsTrueValue(hiprof,"PropagateTables", "TRUE") ) {
00134       RemoveHiBlobs(*(ocube->label()));
00135     }
00136 
00137 //  Set specified profile if entered by user
00138     if (ui.WasEntered("PROFILE")) {
00139       hiconf.selectProfile(ui.GetAsString("PROFILE"));
00140     }
00141 
00142 
00143 //  Add OPATH parameter to profiles
00144     if (ui.WasEntered("OPATH")) {
00145       hiconf.add("OPATH",ui.GetAsString("OPATH"));
00146     }
00147     else {
00148       //  Set default to output directory
00149       hiconf.add("OPATH", FileName(ocube->fileName()).path());
00150     }
00151 
00152 //  Do I/F output DN conversions
00153     QString units = ui.GetString("UNITS");
00154 
00155     //  Allocate the calibration list
00156     calVars = new MatrixList;
00157 
00158 //  Set up access to HiRISE ancillary data (tables, blobs) here.  Note it they
00159 //  are gone, this will error out. See PropagateTables in conf file.
00160     HiCalData caldata(*hifrom);
00161 
00162 ////////////////////////////////////////////////////////////////////////////
00163 //  Drift Correction (Zf) using buffer pixels
00164 //    Extracts specified regions of the calibration buffer pixels and runs
00165 //    series of lowpass filters.  Apply spline fit if any missing data
00166 //    remains.  Config file contains parameters for this operation.
00167     procStep = "ZeroBufferSmooth module";
00168     hiconf.selectProfile("ZeroBufferSmooth");
00169     hiprof = hiconf.getMatrixProfile();
00170     HiHistory ZbsHist;
00171     ZbsHist.add("Profile["+ hiprof.Name()+"]");
00172     if ( !SkipModule(hiprof) ) {
00173       ZeroBufferSmooth zbs(caldata, hiconf);
00174       calVars->add("ZeroBufferSmooth", zbs.ref());
00175       ZbsHist = zbs.History();
00176       if ( hiprof.exists("DumpModuleFile") ) {
00177         zbs.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
00178       }
00179     }
00180     else {
00181       //  NOT RECOMMENDED!  This is required for the next step!
00182       //  SURELY must be skipped with ZeroBufferSmooth step as well!
00183       calVars->add("ZeroBufferSmooth", HiVector(nlines, 0.0));
00184       ZbsHist.add("Debug::SkipModule invoked!");
00185     }
00186 
00187 /////////////////////////////////////////////////////////////////////
00188 // ZeroBufferFit
00189 //  Compute second level of drift correction.  The high level noise
00190 //  is removed from a modeled non-linear fit.
00191 //
00192     procStep = "ZeroBufferFit module";
00193     HiHistory ZbfHist;
00194     hiconf.selectProfile("ZeroBufferFit");
00195     hiprof = hiconf.getMatrixProfile();
00196     ZbfHist.add("Profile["+ hiprof.Name()+"]");
00197     if (!SkipModule(hiprof) ) {
00198       ZeroBufferFit zbf(hiconf);
00199 
00200       calVars->add(hiconf.getProfileName(),
00201                    zbf.Normalize(zbf.Solve(calVars->get("ZeroBufferSmooth"))));
00202       ZbfHist = zbf.History();
00203       if ( hiprof.exists("DumpModuleFile") ) {
00204         zbf.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
00205       }
00206     }
00207     else {
00208       calVars->add(hiconf.getProfileName(), HiVector(nlines, 0.0));
00209       ZbfHist.add("Debug::SkipModule invoked!");
00210     }
00211 
00212 
00213  ////////////////////////////////////////////////////////////////////
00214  //  ZeroReverse
00215     procStep = "ZeroReverse module";
00216     hiconf.selectProfile("ZeroReverse");
00217     hiprof = hiconf.getMatrixProfile();
00218     HiHistory ZrHist;
00219     ZrHist.add("Profile["+ hiprof.Name()+"]");
00220     if ( !SkipModule(hiprof) ) {
00221       ZeroReverse zr(caldata, hiconf);
00222       calVars->add(hiconf.getProfileName(), zr.ref());
00223       ZrHist = zr.History();
00224       if ( hiprof.exists("DumpModuleFile") ) {
00225         zr.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
00226       }
00227     }
00228     else {
00229       calVars->add(hiconf.getProfileName(), HiVector(nsamps, 0.0));
00230       ZrHist.add("Debug::SkipModule invoked!");
00231     }
00232 
00233 /////////////////////////////////////////////////////////////////
00234 // ZeroDark removes dark current
00235 //
00236     procStep = "ZeroDark module";
00237     hiconf.selectProfile("ZeroDark");
00238     hiprof =  hiconf.getMatrixProfile();
00239     HiHistory ZdHist;
00240     ZdHist.add("Profile["+ hiprof.Name()+"]");
00241     if ( !SkipModule(hiprof) ) {
00242       ZeroDark zd(hiconf);
00243       calVars->add(hiconf.getProfileName(), zd.ref());
00244       ZdHist = zd.History();
00245       if ( hiprof.exists("DumpModuleFile") ) {
00246         zd.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
00247       }
00248     }
00249     else {
00250       calVars->add(hiconf.getProfileName(), HiVector(nsamps, 0.0));
00251       ZdHist.add("Debug::SkipModule invoked!");
00252     }
00253 
00254 ////////////////////////////////////////////////////////////////////
00255 // GainLineDrift correct for gain-based drift
00256 //
00257     procStep = "GainLineDrift module";
00258     hiconf.selectProfile("GainLineDrift");
00259     hiprof = hiconf.getMatrixProfile();
00260     HiHistory GldHist;
00261     GldHist.add("Profile["+ hiprof.Name()+"]");
00262     if ( !SkipModule(hiprof) ) {
00263       GainLineDrift gld(hiconf);
00264       calVars->add(hiconf.getProfileName(), gld.ref());
00265       GldHist = gld.History();
00266       if ( hiprof.exists("DumpModuleFile") ) {
00267         gld.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
00268       }
00269     }
00270     else {
00271       calVars->add(hiconf.getProfileName(), HiVector(nlines, 1.0));
00272       GldHist.add("Debug::SkipModule invoked!");
00273     }
00274 
00275 ////////////////////////////////////////////////////////////////////
00276 //  GainNonLinearity  Correct for non-linear gain
00277     procStep = "GainNonLinearity module";
00278     hiconf.selectProfile("GainNonLinearity");
00279     hiprof =  hiconf.getMatrixProfile();
00280     HiHistory GnlHist;
00281     GnlHist.add("Profile["+ hiprof.Name()+"]");
00282     if ( !SkipModule(hiprof) ) {
00283       GainNonLinearity gnl(hiconf);
00284       calVars->add(hiconf.getProfileName(), gnl.ref());
00285       GnlHist = gnl.History();
00286       if ( hiprof.exists("DumpModuleFile") ) {
00287         gnl.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
00288       }
00289     }
00290     else {
00291       calVars->add(hiconf.getProfileName(), HiVector(1, 0.0));
00292       GnlHist.add("Debug::SkipModule invoked!");
00293     }
00294 
00295 ////////////////////////////////////////////////////////////////////
00296 //  GainChannelNormalize  Correct for sample gain with the G matrix
00297     procStep = "GainChannelNormalize module";
00298     hiconf.selectProfile("GainChannelNormalize");
00299     hiprof =  hiconf.getMatrixProfile();
00300     HiHistory GcnHist;
00301     GcnHist.add("Profile["+ hiprof.Name()+"]");
00302     if ( !SkipModule(hiprof) ) {
00303       GainChannelNormalize gcn(hiconf);
00304       calVars->add(hiconf.getProfileName(), gcn.ref());
00305       GcnHist = gcn.History();
00306       if ( hiprof.exists("DumpModuleFile") ) {
00307         gcn.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
00308       }
00309     }
00310     else {
00311       calVars->add(hiconf.getProfileName(), HiVector(nsamps, 1.0));
00312       GcnHist.add("Debug::SkipModule invoked!");
00313     }
00314 
00315 ////////////////////////////////////////////////////////////////////
00316 //  GainFlatField  Flat field correction with A matrix
00317     procStep = "GainFlatField module";
00318     hiconf.selectProfile("GainFlatField");
00319     hiprof =  hiconf.getMatrixProfile();
00320     HiHistory GffHist;
00321     GffHist.add("Profile["+ hiprof.Name()+"]");
00322     if ( !SkipModule(hiprof) ) {
00323       GainFlatField gff(hiconf);
00324       calVars->add(hiconf.getProfileName(), gff.ref());
00325       GffHist = gff.History();
00326       if ( hiprof.exists("DumpModuleFile") ) {
00327         gff.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
00328       }
00329     }
00330     else {
00331       calVars->add(hiconf.getProfileName(), HiVector(nsamps, 1.0));
00332       GffHist.add("Debug::SkipModule invoked!");
00333     }
00334 
00335 ////////////////////////////////////////////////////////////////////
00336 // GainTemperature -  Temperature-dependant gain correction
00337     procStep = "GainTemperature module";
00338     hiconf.selectProfile("GainTemperature");
00339     hiprof =  hiconf.getMatrixProfile();
00340     HiHistory GtHist;
00341     GtHist.add("Profile["+ hiprof.Name()+"]");
00342     if ( !SkipModule(hiprof) ) {
00343       GainTemperature gt(hiconf);
00344       calVars->add(hiconf.getProfileName(), gt.ref());
00345       GtHist = gt.History();
00346       if ( hiprof.exists("DumpModuleFile") ) {
00347         gt.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
00348       }
00349     }
00350     else {
00351       calVars->add(hiconf.getProfileName(), HiVector(nsamps, 1.0));
00352       GtHist.add("Debug::SkipModule invoked!");
00353     }
00354 
00355 ////////////////////////////////////////////////////////////////////
00356 //  GainUnitConversion converts to requested units
00357 //
00358     procStep = "GainUnitConversion module";
00359     hiconf.selectProfile("GainUnitConversion");
00360     hiprof = hiconf.getMatrixProfile();
00361     HiHistory GucHist;
00362     GucHist.add("Profile["+ hiprof.Name()+"]");
00363     if ( !SkipModule(hiprof) ) {
00364       GainUnitConversion guc(hiconf, units);
00365       calVars->add(hiconf.getProfileName(), guc.ref());
00366       GucHist = guc.History();
00367       if ( hiprof.exists("DumpModuleFile") ) {
00368         guc.Dump(hiconf.getMatrixSource("DumpModuleFile",hiprof));
00369       }
00370     }
00371     else {
00372       calVars->add(hiconf.getProfileName(), HiVector(1,1.0));
00373       GucHist.add("Debug::SkipModule invoked!");
00374       GucHist.add("Units[Unknown]");
00375     }
00376 
00377     //  Reset the profile selection to default
00378     hiconf.selectProfile();
00379 
00380 //----------------------------------------------------------------------
00381 //
00382 /////////////////////////////////////////////////////////////////////////
00383 //  Call the processing function
00384     procStep = "calibration phase";
00385     p.StartProcess(calibrate);
00386 
00387     // Get the default profile for logging purposes
00388     hiprof = hiconf.getMatrixProfile();
00389     const QString conf_file = hiconf.filepath(conf);
00390 
00391     // Quitely dumps parameter history to alternative format file.  This
00392     // is completely controlled by the configuration file
00393     if ( hiprof.exists("DumpHistoryFile") ) {
00394       procStep = "logging/reporting phase";
00395       FileName hdump(hiconf.getMatrixSource("DumpHistoryFile",hiprof));
00396       QString hdumpFile = hdump.expanded();
00397       ofstream ofile(hdumpFile.toAscii().data(), ios::out);
00398       if (!ofile) {
00399         QString mess = "Unable to open/create history dump file " +
00400                       hdump.expanded();
00401         IException(IException::User, mess, _FILEINFO_).print();
00402       }
00403       else {
00404         ofile << "Program:  " << hical_program << endl;
00405         ofile << "RunTime:  " << hical_runtime << endl;
00406         ofile << "Version:  " << hical_version << endl;
00407         ofile << "Revision: " << hical_revision << endl << endl;
00408 
00409         ofile << "FROM:     " << hifrom->fileName() << endl;
00410         ofile << "TO:       " << ocube->fileName()  << endl;
00411         ofile << "CONF:     " << conf_file  << endl << endl;
00412 
00413         ofile << "/* " << hical_program << " application equation */\n"
00414               << "/* hdn = (idn - ZeroBufferFit(ZeroBufferSmooth) - ZeroReverse - ZeroDark) */\n"
00415               << "/* odn = hdn / GainLineDrift * GainNonLinearity * GainChannelNormalize */\n"
00416               << "/*           * GainFlatField  * GainTemperature / GainUnitConversion */\n\n";
00417 
00418         ofile << "****** PARAMETER GENERATION HISTORY *******" << endl;
00419         ofile << "\nZeroBufferSmooth   = " << ZbsHist << endl;
00420         ofile << "\nZeroBufferFit   = " << ZbfHist << endl;
00421         ofile << "\nZeroReverse   = " << ZrHist << endl;
00422         ofile << "\nZeroDark   = " << ZdHist << endl;
00423         ofile << "\nGainLineDrift   = " << GldHist << endl;
00424         ofile << "\nGainNonLinearity   = " << GnlHist << endl;
00425         ofile << "\nGainChannelNormalize = " << GcnHist << endl;
00426         ofile << "\nGainFlatField   = " << GffHist << endl;
00427         ofile << "\nGainTemperature   = " << GtHist << endl;
00428         ofile << "\nGainUnitConversion = " << GucHist << endl;
00429 
00430         ofile.close();
00431       }
00432     }
00433 
00434 //  Ensure the RadiometricCalibration group is out there
00435     const QString rcalGroup("RadiometricCalibration");
00436     if (!ocube->hasGroup(rcalGroup)) {
00437       PvlGroup temp(rcalGroup);
00438       ocube->putGroup(temp);
00439     }
00440 
00441     PvlGroup &rcal = ocube->group(rcalGroup);
00442     rcal += PvlKeyword("Program", hical_program);
00443     rcal += PvlKeyword("RunTime", hical_runtime);
00444     rcal += PvlKeyword("Version",hical_version);
00445     rcal += PvlKeyword("Revision",hical_revision);
00446 
00447     PvlKeyword key("Conf", conf_file);
00448     key.addCommentWrapped("/* " + hical_program + " application equation */");
00449     key.addComment("/* hdn = idn - ZeroBufferFit(ZeroBufferSmooth) */");
00450     key.addComment("/*           - ZeroReverse - ZeroDark */");
00451     key.addComment("/* odn = hdn / GainLineDrift * GainNonLinearity */");
00452     key.addComment("/*           * GainChannelNormalize * GainFlatField */");
00453     key.addComment("/*           * GainTemperature / GainUnitConversion */");
00454     rcal += key;
00455 
00456     //  Record parameter generation history.  Controllable in configuration
00457     //  file.  Note this is optional because of a BUG!! in the ISIS label
00458     //  writer as this application was initially developed
00459     if ( IsEqual(ConfKey(hiprof,"LogParameterHistory",QString("TRUE")),"TRUE")) {
00460       rcal += ZbsHist.makekey("ZeroBufferSmooth");
00461       rcal += ZbfHist.makekey("ZeroBufferFit");
00462       rcal += ZrHist.makekey("ZeroReverse");
00463       rcal += ZdHist.makekey("ZeroDark");
00464       rcal += GldHist.makekey("GainLineDrift");
00465       rcal += GnlHist.makekey("GainNonLinearity");
00466       rcal += GcnHist.makekey("GainChannelNormalize");
00467       rcal += GffHist.makekey("GainFlatField");
00468       rcal += GtHist.makekey("GainTemperature");
00469       rcal += GucHist.makekey("GainUnitConversion");
00470     }
00471 
00472     p.EndProcess();
00473   }
00474   catch (IException &ie) {
00475     delete calVars;
00476     calVars = 0;
00477     QString mess = "Failed in " + procStep;
00478     throw IException(ie, IException::User, mess, _FILEINFO_);
00479   }
00480 
00481 // Clean up parameters
00482   delete calVars;
00483   calVars = 0;
00484 }
00485