diff --git a/isis/src/hayabusa2/apps/hyb2onc2isis/hyb2onc2isis.cpp b/isis/src/hayabusa2/apps/hyb2onc2isis/hyb2onc2isis.cpp index 987e6ae55d..222d27ae4c 100644 --- a/isis/src/hayabusa2/apps/hyb2onc2isis/hyb2onc2isis.cpp +++ b/isis/src/hayabusa2/apps/hyb2onc2isis/hyb2onc2isis.cpp @@ -53,6 +53,7 @@ void hyb2onc2isis(UserInterface &ui) { Pvl hyb2onc2isis(QString fitsFileName, QString outputCubeFileName, CubeAttributeOutput att, QString target) { + Pvl finalPvl; ProcessImportFits importFits; importFits.setFitsFile(FileName(fitsFileName)); importFits.setProcessFileStructure(0); @@ -138,7 +139,7 @@ Pvl hyb2onc2isis(QString fitsFileName, QString outputCubeFileName, CubeAttribute // Translate the Instrument group - FileName transFile(transDir + "hyb2oncInstrument1.trn"); + FileName transFile(transDir + "hyb2oncInstrument.trn"); if (updatedKeywords) { transFile = transDir+"hyb2oncInstrumentUpdated.trn"; @@ -255,14 +256,16 @@ Pvl hyb2onc2isis(QString fitsFileName, QString outputCubeFileName, CubeAttribute // Convert the image data importFits.Progress()->SetText("Importing Hayabusa2 image"); - Pvl finalLabel = *(outputCube->label() ); + importFits.StartProcess(); importFits.Finalize(); - return finalLabel; - + + return outputLabel; + } -} + + diff --git a/isis/src/hayabusa2/apps/hyb2onccal/Hyb2OncCalUtils.h b/isis/src/hayabusa2/apps/hyb2onccal/Hyb2OncCalUtils.h index 967cf07920..dfefe3e3fe 100644 --- a/isis/src/hayabusa2/apps/hyb2onccal/Hyb2OncCalUtils.h +++ b/isis/src/hayabusa2/apps/hyb2onccal/Hyb2OncCalUtils.h @@ -120,14 +120,12 @@ static double g_compfactor(1.0); // Default if OutputMode = LOSS-LESS; 16.0 for static QString g_iofCorrection("IOF"); //!< Is I/F correction to be applied? // I/F variables -static double g_solarDist(1.0); /**< Distance from the Sun to the target body -(used to calculate g_iof) */ -static double g_iof(1.0); //!< I/F conversion value + static double g_iofScale(1.0); static double g_solarFlux(1.0); //!< The solar flux (used to calculate g_iof). -// TODO: we do not have this conversion factor for Hayabusa 2 ONC. -static double g_v_standard(1.0); -// static double g_v_standard(3.42E-3);//!< Base conversion for all filters (Tbl. 9) +static double g_sensitivity(1.0); +static double g_effectiveBandwidth(1.0); + namespace Isis { @@ -186,12 +184,74 @@ bool newton_rapheson(double Iobs,double x0, double g[3],double &result, double e /** -* @brief Apply radiometric correction to each line of an AMICA image. + * @brief linearFun: The linear correction function (used by the newton_rapheson method) + * @author 2019-02-12 Tyler Wilson + * @param Iobs: The observed intensity + * @param x: The ideal intensity. + * @param g: The vector of empirically derived coefficients for the third-order polynomial + * modelling the linear correction (for DN values < 3400 DN) + * @return The value of the function at the point x. + */ +double linearFun(double Iobs,double x, double g[3]) { + return Iobs - g[0]*x -g[1]*pow(x,2.0) -g[2]*pow(x,3.0); + +} + +/** + * @brief dFun: The first-order derivative of linearFun + * @author 2019-02-12 Tyler Wilson + * @param x: The ideal intensity. + * @param g: The vector of empirically derived coefficients for the third-order polynomial + * modelling the linear correction (for DN values < 3400 DN) + * @return + */ +double dFun(double x, double g[3]) { + return -g[0] - 2*g[1]*x -3*g[2]*pow(x,2.0); + +} + +/** + * @brief newton_rapheson + * @author 2019-02-12 Tyler Wilson + * @param Iobs: The observed DN intensity + * @param x0: The starting value for the Newton-Rapheson method + * @param g: A vector of the coefficients for the linearity function. It is a third-order + * polynomial. + * @param result: The final approximation of the root of the equation. + * @param epsilon: The tolerance on the final solution. + * @return A root of the linearity correction function, centered near the origin. + */ +bool newton_rapheson(double Iobs,double x0, double g[3],double &result, double epsilon=1e-6 ) { + + double x[2]; + double dx = 1.0; + int iter = 0; + int maxIterations=500; + x[0] = x0; + while (dx > epsilon) { + + x[1]=x[0] - linearFun(Iobs,x[0],g)/dFun(x[0],g); + dx = fabs(x[1]-x[0]) ; + x[0]=x[1]; + iter++; + if (iter > maxIterations) { + + return false; + } + } + result = x[1]; + return true; +} + + +/** +* @brief Apply radiometric correction to each line of a Hayabusa2 image. * @author 2016-03-30 Kris Becker * @param in Raw image and flat field * @param out Radometrically corrected image * @internal * @history 2017-07-2017 Ian Humphrey & Kaj Williams - Adapted from amicacal. +* @history 2019-02-12 Tyler Wilson - Modified to support new calibration settings/formulas. */ void Calibrate(vector& in, vector& out) { @@ -218,8 +278,6 @@ void Calibrate(vector& in, vector& out) { for (int i = 0; i < imageIn.size(); i++) { imageOut[i] = imageIn[i]*pow(2.0,12-g_bitDepth); - - // Check for special pixel in input image and pass through if ( IsSpecial(imageOut[i]) ) { imageOut[i] = imageIn[i]; @@ -227,12 +285,10 @@ void Calibrate(vector& in, vector& out) { } // Apply compression factor here to raise LOSSY dns to proper response - imageOut[i] *= g_compfactor; // 1) BIAS Removal - Only needed if not on-board corrected - if ( !g_onBoardSmearCorrection ) { @@ -241,31 +297,24 @@ void Calibrate(vector& in, vector& out) { continue; } else { + imageOut[i] = imageOut[i] - g_bias; } } - double dn = imageOut[i]; - double linearCorrection; + + double dn = imageOut[i]; double result = 1.0; double x0 = 1.0; - newton_rapheson(imageOut[i],x0, g_L,result ); + newton_rapheson(dn,x0, g_L,result ); imageOut[i] = result; - //qDebug() << dn << ","<< result; - // DARK Current - imageOut[i] = imageOut[i] - g_darkCurrent; - - + imageOut[i] = imageOut[i] - g_darkCurrent; - // READOUT Smear Removal - Not needed if on-board corrected. Binning is - // accounted for in computation of c1 before loop. - // if (nsubImages <= 1) { - // imageOut[i] = c1*(imageOut[i] - smear); - // } + //Smear correction if (!g_onBoardSmearCorrection) { double smear = 0; @@ -277,10 +326,6 @@ void Calibrate(vector& in, vector& out) { } - //Linearity Correction - //In the SIS this adjustment is made just after the bias, but - //in the Calibration paper it happens just before the flat field correction. - // FLATFIELD correction // Check for any special pixels in the flat field (unlikely) @@ -294,16 +339,11 @@ void Calibrate(vector& in, vector& out) { } else { if (flatField[i] != 0) { - imageOut[i] /= flatField[i]; + imageOut[i] /= (flatField[i]*g_sensitivity*g_texp); } } } - // TODO: once the radiance values are known for each band, we can correctly compute I/F. - // For now, g_iof is 1, so output will be in DNs. - // 7) I/F or Radiance Conversion (or g_iof might = 1, in which case the output will be in DNs) - imageOut[i] *= g_iof; - } @@ -312,107 +352,6 @@ void Calibrate(vector& in, vector& out) { -/** - * @brief Load required NAIF kernels required for timing needs. - * - * This method maintains the loading of kernels for HAYABUSA timing and - * planetary body ephemerides to support time and relative positions of planet - * bodies. - */ -/* Helper function for sunDistanceAu, don't need this until we have radiance calibration - parameters for Hayabusa2 ONC-T filters to calculate radiance and I/F -static void loadNaifTiming() { - static bool naifLoaded = false; - if (!naifLoaded) { - -// Load the NAIF kernels to determine timing data - Isis::FileName leapseconds("$base/kernels/lsk/naif????.tls"); - leapseconds = leapseconds.highestVersion(); - Isis::FileName sclk("$hayabusa2/kernels/sclk/hyb2_20141203-20161231_v01.tsc"); - Isis::FileName pck1("$hayabusa2/kernels/tspk/de430.bsp"); - Isis::FileName pck2("$hayabusa2/kernels/tspk/jup329.bsp"); - Isis::FileName pck3("$hayabusa2/kernels/tspk/sat375.bsp"); - Isis::FileName pck4("$hayabusa2/kernels/spk/hyb2_20141203-20141214_0001m_final_ver1.oem.bsp"); - Isis::FileName pck5("$hayabusa2/kernels/spk/hyb2_20141203-20151231_0001h_final_ver1.oem.bsp"); - Isis::FileName pck6("$hayabusa2/kernels/spk/hyb2_20151123-20151213_0001m_final_ver1.oem.bsp"); - -// Load the kernels - QString leapsecondsName(leapseconds.expanded()); - QString sclkName(sclk.expanded()); - - QString pckName1(pck1.expanded()); - QString pckName2(pck2.expanded()); - QString pckName3(pck3.expanded()); - QString pckName4(pck4.expanded()); - QString pckName5(pck5.expanded()); - QString pckName6(pck6.expanded()); - - furnsh_c(leapsecondsName.toLatin1().data()); - furnsh_c(sclkName.toLatin1().data()); - - furnsh_c(pckName1.toLatin1().data()); - furnsh_c(pckName2.toLatin1().data()); - furnsh_c(pckName3.toLatin1().data()); - furnsh_c(pckName4.toLatin1().data()); - furnsh_c(pckName5.toLatin1().data()); - furnsh_c(pckName6.toLatin1().data()); - - -// Ensure it is loaded only once - naifLoaded = true; - } - return; -} -*/ - - -/** - * @brief Computes the distance from the Sun to the observed body. - * - * This method requires the appropriate NAIK kernels to be loaded that - * provides instrument time support, leap seconds and planet body ephemeris. - * - * @return @b double Distance in AU between Sun and observed body. - */ -/* commented out until we have radiance values (RAD/IOF group in calibration trn) for Hayabusa2. - static bool sunDistanceAU(const QString &scStartTime, - const QString &target, - double &sunDist) { - - // Ensure NAIF kernels are loaded - loadNaifTiming(); - sunDist = 1.0; - - // Determine if the target is a valid NAIF target - SpiceInt tcode; - SpiceBoolean found; - bodn2c_c(target.toLatin1().data(), &tcode, &found); - - if (!found) return (false); - - // Convert starttime to et - double obsStartTime; - scs2e_c(-37, scStartTime.toLatin1().data(), &obsStartTime); - - // Get the vector from target to sun and determine its length - double sunv[3]; - double lt; - spkpos_c(target.toLatin1().data(), obsStartTime, "J2000", "LT+S", "sun", - sunv, <); - NaifStatus::CheckErrors(); - - double sunkm = vnorm_c(sunv); - - - // Return in AU units - sunDist = sunkm / 1.49597870691E8; - - //cout << "sunDist = " << sunDist << endl; - return (true); -} -*/ - - /** * @brief Translates a 1-banded Isis::Cube to an OpenMat object * @@ -538,13 +477,16 @@ void translate(Cube *flatField,double *transform, QString fname) { * @return FileName Path and name of flat file file * @internal * @history 2017-07-27 Ian Humphrey & Kaj Williams - Adapted from amicacal. +* @history 2019-02-12 Tyler Wilson - Modified to support new calibration settings/formulas. */ FileName DetermineFlatFieldFile(const QString &filter) { QString fileName = "$hayabusa2/calibration/flatfield/"; // FileName consists of binned/notbinned, camera, and filter - fileName += "flat_" + filter.toLower() + ".cub"; + fileName += "flat_" + filter.toLower() + "_norm.cub"; + FileName final(fileName); + return final; } @@ -574,6 +516,8 @@ QString loadCalibrationVariables(const QString &config) { PvlGroup &DarkCurrent = g_configFile.findGroup("DarkCurrent"); PvlGroup &Smear = g_configFile.findGroup("SmearRemoval"); PvlGroup &solar = g_configFile.findGroup("SOLARFLUX"); + PvlGroup &sensitivity = g_configFile.findGroup("SENSITIVITYFACTOR"); + PvlGroup & effectiveBW = g_configFile.findGroup("EFFECTIVEBW"); PvlGroup &linearity = g_configFile.findGroup("Linearity"); // PvlGroup &iof = g_configFile.findGroup("RAD"); @@ -614,11 +558,13 @@ QString loadCalibrationVariables(const QString &config) { // Compute BIAS correction factor (it's a constant so do it once!) g_bias = g_b0+g_b1*g_CCD_T_temperature+g_b2*g_ECT_T_temperature; - g_bias *= (g_bae0 + g_bae1*g_AEtemperature); //correction factor - qDebug() << "Bias: " << g_bias; + //double correction_factor = (g_bae0 + g_bae1*g_AEtemperature); + g_bias *= (g_bae0 + g_bae1*g_AEtemperature); //bias correction factor // Load the Solar Flux for the specific filter g_solarFlux = solar[g_filter.toLower()]; + g_sensitivity = sensitivity[g_filter.toLower()]; + g_effectiveBandwidth = effectiveBW[g_filter.toLower()]; //Load the linearity variables g_L[0] = linearity["L"][0].toDouble(); @@ -626,16 +572,10 @@ QString loadCalibrationVariables(const QString &config) { g_L[2] = linearity["L"][2].toDouble(); - - // radiance = g_v_standard * g_iofScale - // iof = radiance * pi *dist_au^2 - // g_iofScale = iof[g_filter]; - return ( calibFile.original() ); } - } #endif diff --git a/isis/src/hayabusa2/apps/hyb2onccal/hyb2onccal.cpp b/isis/src/hayabusa2/apps/hyb2onccal/hyb2onccal.cpp index 145fa318df..a6a6778e84 100644 --- a/isis/src/hayabusa2/apps/hyb2onccal/hyb2onccal.cpp +++ b/isis/src/hayabusa2/apps/hyb2onccal/hyb2onccal.cpp @@ -39,9 +39,6 @@ using namespace std; -//void Calibrate(vector& in, vector& out); - -//QString loadCalibrationVariables(const QString &config); @@ -126,7 +123,6 @@ void hyb2onccal(UserInterface &ui) { catch (IException &e) { QString msg = "Unable to read [BitDepth] keyword in the Instrument group " "from input file [" + icube->fileName() + "]"; - //qDebug() << msg; g_bitDepth = 12; } @@ -151,7 +147,7 @@ void hyb2onccal(UserInterface &ui) { try { g_AEtemperature = inst["ONCAETemperature"][0].toDouble(); - qDebug() << "g_AEtemperature: " << g_AEtemperature; + } catch(IException &e) { QString msg = "Unable to read [ONCAETemperature] keyword in the Instrument group " @@ -162,7 +158,7 @@ void hyb2onccal(UserInterface &ui) { try { g_CCD_T_temperature = inst["ONCTCCDTemperature"][0].toDouble(); - qDebug() << "g_CCD_T_Temperature: " << g_CCD_T_temperature; + } catch(IException &e) { QString msg = "Unable to read [ONCTCCDTemperature] keyword in the Instrument group " @@ -172,7 +168,7 @@ void hyb2onccal(UserInterface &ui) { try { g_ECT_T_temperature = inst["ONCTElectricCircuitTemperature"][0].toDouble(); - qDebug() << "g_ECT_T_temperature: " << g_ECT_T_temperature; + } catch(IException &e) { QString msg = "Unable to read [ONCTElectricCircuitTemperature] keyword in the Instrument group " @@ -208,11 +204,16 @@ void hyb2onccal(UserInterface &ui) { if (smearCorrection=="ONBOARD") { g_onBoardSmearCorrection=true; } + else { + qDebug() << icube->fileName(); + + } QString compmode = inst["Compression"]; // TODO: verify that the compression factor/scale is actually 16 for compressed Hayabusa2 images. g_compfactor = ( "lossy" == compmode.toLower() ) ? 16.0 : 1.0; + QString target = inst["TargetName"]; g_target = target; @@ -273,14 +274,12 @@ void hyb2onccal(UserInterface &ui) { Cube *flatOriginal = new Cube(flatfile.expanded() ); - //int transform[5] = {binning,startSample,startLine,lastSample,lastLine}; double alphaStartSample = alphaCube["AlphaStartingSample"][0].toDouble(); double alphaStartLine = alphaCube["AlphaStartingLine"][0].toDouble(); double alphaEndSample = alphaCube["AlphaEndingSample"][0].toDouble(); double alphaEndLine = alphaCube["AlphaEndingLine"][0].toDouble(); - //int transform[5] = {binning,startSample,startLine,lastSample,lastLine}; double transform[5] = {(double)binning,alphaStartSample,alphaStartLine,alphaEndSample,alphaEndLine}; // Translates and scales the flatfield image. Scaling @@ -302,13 +301,14 @@ void hyb2onccal(UserInterface &ui) { Cube *ocube = p.SetOutputCube("TO"); - //QString fname = ocube->fileName(); + QString calfile = loadCalibrationVariables(ui.GetAsString("CONFIG")); + + g_timeRatio = g_Tvct/(g_texp + g_Tvct); - g_iof = 1.0; // Units of DN QString g_units = "DN"; // if ( "radiance" == g_iofCorrection.toLower() ) { @@ -350,7 +350,10 @@ void hyb2onccal(UserInterface &ui) { calibrationLog.addKeyword(PvlKeyword("CalibrationFile", calfile)); calibrationLog.addKeyword(PvlKeyword("FlatFieldFile", flatfile.originalPath() + "/" + flatfile.name())); - calibrationLog.addKeyword(PvlKeyword("CompressionFactor", toString(g_compfactor, 2))); + + PvlKeyword sensitivityFactor("SensitivityFactor"); + sensitivityFactor.addValue(toString(g_sensitivity,16)); + calibrationLog.addKeyword(sensitivityFactor); // Parameters PvlKeyword bn("Bias_Bn"); @@ -362,18 +365,7 @@ void hyb2onccal(UserInterface &ui) { PvlKeyword bnae("Bias_AECorrection"); bnae.addValue(toString(g_bae0, 8)); bnae.addValue(toString(g_bae1, 8)); - calibrationLog.addKeyword(bnae); - -#if 0 - PvlKeyword linearityCoefs("Linearity"); - linearityCoefs.addValue(toString(g_L0,8)); - linearityCoefs.addValue(toString(g_L1,8)); - linearityCoefs.addValue(toString(g_L2,8)); - calibrationLog.addKeyword(linearityCoefs); - -#endif - - + calibrationLog.addKeyword(bnae); calibrationLog.addKeyword(PvlKeyword("Bias_AETemp", toString(g_AEtemperature, 16))); @@ -400,13 +392,27 @@ void hyb2onccal(UserInterface &ui) { calibrationLog.addKeyword(PvlKeyword("Smear_texp", toString(g_texp, 16))); calibrationLog.addKeyword(PvlKeyword("CalibrationUnits", g_iofCorrection)); - calibrationLog.addKeyword(PvlKeyword("RadianceStandard", toString(g_v_standard, 16))); + calibrationLog.addKeyword(PvlKeyword("RadianceScaleFactor", toString(g_iofScale, 16))); - calibrationLog.addKeyword(PvlKeyword("SolarDistance", toString(g_solarDist, 16), "AU")); - calibrationLog.addKeyword(PvlKeyword("SolarFlux", toString(g_solarFlux, 16))); - calibrationLog.addKeyword(PvlKeyword("IOFFactor", toString(g_iof, 16))); + calibrationLog.addKeyword(PvlKeyword("SolarFlux", toString(g_solarFlux, 16))); calibrationLog.addKeyword(PvlKeyword("Units", g_units)); + PvlKeyword linearityCoefs("LinearityCoefficients"); + linearityCoefs.addValue(toString(g_L[0],16)); + linearityCoefs.addValue(toString(g_L[1],16)); + linearityCoefs.addValue(toString(g_L[2],16)); + calibrationLog.addKeyword(linearityCoefs); + + PvlKeyword darkCurrentCoefs("DarkCurrentCoefficients"); + darkCurrentCoefs.addValue(toString(g_d0,16)); + darkCurrentCoefs.addValue(toString(g_d1,16)); + calibrationLog.addKeyword(darkCurrentCoefs); + + PvlKeyword darkCurrent("DarkCurrent"); + darkCurrent.addValue(toString(g_darkCurrent,16)); + calibrationLog.addKeyword(darkCurrent); + + // Write Calibration group to output file ocube->putGroup(calibrationLog); Application::Log(calibrationLog); diff --git a/isis/src/hayabusa2/apps/hyb2onccal/hyb2onccal.xml b/isis/src/hayabusa2/apps/hyb2onccal/hyb2onccal.xml index 0a60b43555..0767c19d0b 100644 --- a/isis/src/hayabusa2/apps/hyb2onccal/hyb2onccal.xml +++ b/isis/src/hayabusa2/apps/hyb2onccal/hyb2onccal.xml @@ -1,7 +1,7 @@ - Calibrates Hayabusa 2 ONC-T images. Calibration includes bias and dark current correction as @@ -66,91 +66,6 @@ xsi:noNamespaceSchemaLocation="http://isis.astrogeology.usgs.gov/Schemas/Applica

- - -

Step 2: Flat-field correction:

Performs a correction for pixel-to-pixel variation in CCD response and vignetting @@ -162,93 +77,24 @@ xsi:noNamespaceSchemaLocation="http://isis.astrogeology.usgs.gov/Schemas/Applica pixel in the flat-field image.

- - -

-

Notes

-
    -
  1. - Smear correction is not currently provided, as we do not have the readout time - (charge-transfer period) for the instrument. -
    2. -
  2. - Similarly, since no radiance values could be found for the instrument filters, the output - cannot be correctly converted to radiance or I/F. -
    3. -
  3. -
    4. -
-

References:

  1. S. Kameda et al. "Preflight Calibration Test Results for Optical Navigation Camera Telescope (ONC-T) Onboard the Hayabusa2 Spacecraft". - Space Sci Rev (2017) 208:17-31.
    2. + Space Sci Rev (2017) 208:17-31. + +
  2. E. Tatsumi et al. "Updated Inflight Calibration of Hayabusa2’s Optical Navigation Camera (ONC) for + Scientific Observations during the Cruise Phase". + https://arxiv.org/abs/1810.11065 [astro-ph.IM] 25 Oct, 2018. +
    3. +
  3. H. Suzuki et al. "Initial inflight calibration for Hayabusa2 optical navigation camera + (ONC) for science observations of asteroid Ryugu". + Icarus. 300. 10.1016/j.icarus.2017.09.011 +

@@ -270,6 +116,9 @@ values for each band for Hayabusa2 ONC. Original version. Adapted from amicacal. + + Modified version to handle the updated calibration values/formulas. + @@ -315,11 +164,14 @@ values for each band for Hayabusa2 ONC. the parameters as needed. *.trn - + /scratch/isis3data/hayabusa2/calibration/onc/hyb2oncCalibration????.trn --> + /usgs/cpkgs/isis3/datalocal/hayabusa2/calibration/onc/hyb2oncCalibration????.trn + +