Skip to content

Commit

Permalink
Added I/F calibration step to hyb2onccal (DOI-USGS#3132)
Browse files Browse the repository at this point in the history 
* Added AlphaCube group to outputlabel in the event the image is cropped.

* Fix libtiff dependency (DOI-USGS#636) (DOI-USGS#644)

* Upgrade libtiff to 4.0.10 (DOI-USGS#636)

* Switch libtiff to 4.0.9 or higher to remove geotiff conflict

* Moved ISIS3 conda-build recipe from ISIS3_deps repository (DOI-USGS#650)

* Fixed warning in Pixel unit tests

* Made a tweak to the dimensions of the alpha cube for cropped images to give it the correct dimensions.

* Changed the transform function so it could process an array of doubles instead of an array of ints.

* Smear correction was being incorrectly applied to onboard smear corrected images and screwing things up.  Fixed it.

* Added gtest capability to hyb2onc2isis

* Removing build numbers from external libraries (DOI-USGS#660)

* Moved ISIS3 conda-build recipe from ISIS3_deps repository

* Un-pinned non-astro build numbers

* Removing build numbers from external libraries in the environment and meta.yeml files

* Final merging

* Added pixel type attribute to the output image of program shadow. Fixes DOI-USGS#5187 (DOI-USGS#659)

* Removed bolding of some text to decrease distraction.

* Fixed some typos.

* Reworded documentation.

* Added section for Environment and PreferemcesSetup in the Getting Started Section. (DOI-USGS#663)

* Updated .gitmodules to use https rather than ssh (DOI-USGS#673)

* Added build type release to conda recipe (DOI-USGS#676)

* Modified isis/tests/CMakeLists.txt as well as hyb2onc2isis to fix unresolved symbol errors (CMakeLists.txt) and to fix some bugs in hyb2onc2isis involved with turning it into a callable function.

* Added w1 test to hyb2onc2isisTests.cpp

* Modified hyb2onccal to be a callable function.

* Updates README with Discourse (DOI-USGS#690)

* Updates README with Discourse

* Update README.md

* Update README.md

* Changed the Pvl label output for Pvl hyb2onc2isis(QString fitsFileName, QString outputCubeFileName, CubeAttributeOutput att, QString target)

* Added Newton-Rapheson method to Hyb2OncCalUtils.h to solve the linearity equation.

* Final push of hayabusa2 changes made to the ingestion/calibration applications.

* Added I/F calibration step.  Cleaned up code base and worked on the documentation.  The documentation will need to be corrected before this is merged into dev.

* Fixed a missing </li> tag in hyb2onccal.xml
  • Loading branch information
@acpaquette
Tyler Wilson authored and acpaquette committed Apr 6, 2021
1 parent e3016f1 commit 9dccfa0
Show file tree
Hide file tree
Showing 5 changed files with 171 additions and 350 deletions.
13 changes: 8 additions & 5 deletions isis/src/hayabusa2/apps/hyb2onc2isis/hyb2onc2isis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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";
Expand Down Expand Up @@ -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;


}

}



228 changes: 84 additions & 144 deletions isis/src/hayabusa2/apps/hyb2onccal/Hyb2OncCalUtils.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 {

Expand Down Expand Up @@ -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<Buffer *>& in, vector<Buffer *>& out) {

Expand All @@ -218,21 +278,17 @@ void Calibrate(vector<Buffer *>& in, vector<Buffer *>& 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];
continue;
}

// 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 ) {


Expand All @@ -241,31 +297,24 @@ void Calibrate(vector<Buffer *>& in, vector<Buffer *>& 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;
Expand All @@ -277,10 +326,6 @@ void Calibrate(vector<Buffer *>& in, vector<Buffer *>& 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)
Expand All @@ -294,16 +339,11 @@ void Calibrate(vector<Buffer *>& in, vector<Buffer *>& 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;

}


Expand All @@ -312,107 +352,6 @@ void Calibrate(vector<Buffer *>& in, vector<Buffer *>& 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, &lt);
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
*
Expand Down Expand Up @@ -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;
}

Expand Down Expand Up @@ -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");

Expand Down Expand Up @@ -614,28 +558,24 @@ 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();
g_L[1] = linearity["L"][1].toDouble();
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
Expand Down
Loading

0 comments on commit 9dccfa0

Please sign in to comment.