Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update lronacpho to use LROC's newer photometric model and parameters file #4519

Merged
merged 6 commits into from
Apr 28, 2022
Merged
Show file tree
Hide file tree
Changes from 5 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
101 changes: 59 additions & 42 deletions isis/src/lro/apps/lronacpho/LROCEmpirical.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@ find files of those names at the top level of this repository. **/
#include <iostream>
#include <memory>
#include <sstream>

#include <SpiceUsr.h>
#include <SpiceZfc.h>
#include <SpiceZmc.h>
Expand All @@ -38,13 +37,11 @@ namespace Isis {
init(pvl, cube);
}


/**
* Destructor
*/
LROCEmpirical::~LROCEmpirical() {};


/**
* @brief Initialize class from input PVL and Cube files
*
Expand Down Expand Up @@ -120,7 +117,6 @@ namespace Isis {
return;
}


/**
* @brief Method to get photometric property given angles
*
Expand Down Expand Up @@ -154,7 +150,6 @@ namespace Isis {
return (m_bandpho[band - 1].phoStd / ph);
}


/**
* @brief Performs actual photometric correction calculations
*
Expand All @@ -173,7 +168,8 @@ namespace Isis {
* @internal
* @history 2016-08-15 Victor Silva - Adapted code from lrowacpho application
* written by Kris Becker
*
* @history 2021-03-12 Victor Silva - Added b parameters for 2019 version of
* LROC Empirical algorithm
*/
double LROCEmpirical::photometry( const Parameters &parms, double i, double e, double g ) const {
// Ensure problematic values are adjusted
Expand All @@ -193,12 +189,20 @@ namespace Isis {
double mu = cos(e);
double mu0 = cos(i);
double alpha = g;
double rcal = exp(parms.a0 + parms.a1 * alpha + parms.a2 * mu + parms.a3 * mu0);
double rcal;

if (parms.algoVersion == 2014 || parms.algoVersion == 0)
rcal = exp(parms.aTerms[0] + parms.aTerms[1] * alpha + parms.aTerms[2] * mu + parms.aTerms[3] * mu0);
else if (parms.algoVersion == 2019)
rcal = mu0 / (mu + mu0) * exp(parms.bTerms[0] + parms.bTerms[1] * (alpha * alpha) + parms.bTerms[2] * alpha + parms.bTerms[3] * sqrt(alpha) + parms.bTerms[4] * mu + parms.bTerms[5] * mu0 + parms.bTerms[6] * (mu0 * mu0) );
else {
std::string mess = "Algorithm version in PVL file not recognized [" + IString(parms.algoVersion) + "]. ";
jessemapel marked this conversation as resolved.
Show resolved Hide resolved
throw IException(IException::Programmer, mess, _FILEINFO_);
}

return (rcal);
}


/**
* @brief Return parameters used for all bands
*
Expand All @@ -212,23 +216,27 @@ namespace Isis {
* @internal
* @history 2016-08-15 Victor Silva - Adapted code from lrowacpho application
* written by Kris Becker
* @history 2021-03-12 Victor Silva - Added b parameters for 2019 version of
* LROC Empirical algorithm
*/
void LROCEmpirical::report( PvlContainer &pvl ) {

pvl.addComment("I/F = F(mu, mu0,phase)");
pvl.addComment(" where:");
pvl.addComment(" mu0 = cos(incidence)");
pvl.addComment(" mu = cos(emission)");
pvl.addComment(" F(mu, mu0, phase) = exp(A0 + A1 * phase + A2 * mu + A3 * mu0 ");

if (m_bandpho[0].algoVersion == 2019 )
pvl.addComment(" F(mu, mu0, phase) = mu0 / (mu + mu0) * exp(B0 + B1 * (alpha * alpha) + B2 * alpha + B3 * sqrt(alpha) + B4 * mu + B5 * mu0 + B6 * (mu0 * mu0) )");
else if (m_bandpho[0].algoVersion == 2014 || m_bandpho[0].algoVersion == 0)
pvl.addComment(" F(mu, mu0, phase) = exp (A0 + A1 * phase + A2 * mu + A3 * mu0 ");
else {
std::string mess = "Could not file the correction algorithm name.";
throw IException(IException::Programmer, mess, _FILEINFO_);
}

pvl += PvlKeyword("Algorithm", "LROC_Empirical");
pvl += PvlKeyword("IncRef", toString(m_iRef), "degrees");
pvl += PvlKeyword("EmaRef", toString(m_eRef), "degrees");
pvl += PvlKeyword("Algorithm", "LROC_Empirical");
pvl += PvlKeyword("IncRef", toString(m_iRef), "degrees");
pvl += PvlKeyword("EmaRef", toString(m_eRef), "degrees");
pvl += PvlKeyword("EmaRef", toString(m_eRef), "degrees");
pvl += PvlKeyword("Algorithm", "LROC_Empirical");
pvl += PvlKeyword("AlgorithmVersion", toString(m_bandpho[0].algoVersion), "" );
pvl += PvlKeyword("IncRef", toString(m_iRef), "degrees");
pvl += PvlKeyword("EmaRef", toString(m_eRef), "degrees");
pvl += PvlKeyword("PhaRef", toString(m_gRef), "degrees");
Expand All @@ -238,46 +246,49 @@ namespace Isis {
PvlKeyword bbc("BandBinCenter");
PvlKeyword bbct("BandBinCenterTolerance");
PvlKeyword bbn("BandNumber");
PvlKeyword a0("A0");
PvlKeyword a1("A1");
PvlKeyword a2("A2");
PvlKeyword a3("A3");

std::vector<PvlKeyword> aTermKeywords;
std::vector<PvlKeyword> bTermKeywords;
for (unsigned int i = 0; i < m_bandpho[0].aTerms.size(); i++)
aTermKeywords.push_back(PvlKeyword("A" + toString((int) i)));
for (unsigned int i = 0; i < m_bandpho[0].bTerms.size(); i++)
bTermKeywords.push_back(PvlKeyword("B" + toString((int) i)));

for (unsigned int i = 0; i < m_bandpho.size(); i++) {
Parameters &p = m_bandpho[i];
units.addValue(p.units);
phostd.addValue(toString(p.phoStd));
bbc.addValue(toString(p.wavelength));
bbct.addValue(toString(p.tolerance));
bbn.addValue(toString(p.band));
a0.addValue(toString(p.a0));
a1.addValue(toString(p.a1));
a2.addValue(toString(p.a2));
a3.addValue(toString(p.a3));
Parameters &p = m_bandpho[i];
units.addValue(p.units);
phostd.addValue(toString(p.phoStd));
bbc.addValue(toString(p.wavelength));
bbct.addValue(toString(p.tolerance));
bbn.addValue(toString(p.band));
for (unsigned int j = 0; j < aTermKeywords.size(); j++)
aTermKeywords[j].addValue(toString(p.aTerms[j]));
for (unsigned int j = 0; j < bTermKeywords.size(); j++)
bTermKeywords[j].addValue(toString(p.bTerms[j]));
jessemapel marked this conversation as resolved.
Show resolved Hide resolved
}

pvl += units;
pvl += phostd;
pvl += bbc;
pvl += bbct;
pvl += bbn;
pvl += a0;
pvl += a1;
pvl += a2;
pvl += a3;

for (unsigned int i = 0; i < aTermKeywords.size(); i++)
pvl += aTermKeywords[i];
for (unsigned int i = 0; i < bTermKeywords.size(); i++)
pvl += bTermKeywords[i];

return;
}


/**
* @brief Determine LROC Empirical parameters given a wavewlength
*
* This method determines the set of LROCEmpirical parameters to
* use for a given wavelength. It iterates through all band profiles
* as read from the PVL file and computes the difference between
* wavelength parameter and the BandBinCenter keyword. The absolute
* value of this value is checke against the BandBinCenterTolerance
* value of this value is check against the BandBinCenterTolerance
* parameter and if it is less than or equal to it, a Parameter
* container is returned.
*
Expand Down Expand Up @@ -316,8 +327,6 @@ namespace Isis {
return (Parameters());
}



/*
* @brief Extracts necessary LROCEmprical parameters from profile
*
Expand All @@ -333,18 +342,26 @@ namespace Isis {
*
* @internal
* @history 2016-08-15 Victor Silva - Adapted from the lrowacpho application
written by Kris Becker.
* written by Kris Becker.
*
* @history 2021-03-12 Victor Silva - Added b parameters for 2019 version of
* LROC Empirical algorithm
*
*/
LROCEmpirical::Parameters LROCEmpirical::extract( const DbProfile &profile) const {
Parameters pars;
pars.a1 = toDouble(ConfKey(profile, "A1", toString(0.0)));
pars.a2 = toDouble(ConfKey(profile, "A2", toString(0.0)));
pars.a3 = toDouble(ConfKey(profile, "A3", toString(0.0)));

for (int i=0; i<4; i++)
pars.aTerms.push_back(toDouble(ConfKey(profile, "A" + toString(i), toString(0.0))));
for (int i=0; i<7; i++)
pars.bTerms.push_back(toDouble(ConfKey(profile, "B" + toString(i), toString(0.0))));
jessemapel marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I noticed Jesse suggested changing the toString call with the string literal "0.0", and that you made the change, was this reverted somehow? Or do you need to explicitly cast to a QString?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This must have gotten reverted when I corrected the merge issue. Thanks. Fixed in code.


pars.wavelength = toDouble(ConfKey(profile, "BandBinCenter", toString(Null)));
pars.tolerance = toDouble(ConfKey(profile, "BandBinCenterTolerance", toString(Null)));
// Determine equation units - defaults to Radians
pars.units = ConfKey(profile, "Units", QString("Radians"));
pars.phaUnit = (pars.units.toLower() == "degrees") ? 1.0 : rpd_c();
pars.algoVersion = toInt(ConfKey(profile, "AlgorithmVersion", QString("0")));
jessemapel marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same string question here.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This must have gotten reverted when I corrected the merge issue. Thanks. Fixed in code.


return (pars);
}
Expand Down
30 changes: 19 additions & 11 deletions isis/src/lro/apps/lronacpho/LROCEmpirical.h
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,8 @@ namespace Isis {
*
* @internal
* @history 2016-08-15 - Code adapted from lrowacpho written by Kris Becker
* @history 2021-03-12 Victor Silva - Added b parameters for 2019 version of
* LROC Empirical algorithm
*
*/
class LROCEmpirical : public PhotometricFunction {
Expand All @@ -53,23 +55,29 @@ namespace Isis {
*
* @internal
* @history 2016-08-05 - Code adapted from lrowacpho written by Kris Becker
* @history 2021-03-12 Victor Silva - Added b parameters for 2019 version of
* LROC Empirical algorithm
*/
struct Parameters {
Parameters() : a0(0.0), a1(0.0), a2(0.0), a3(0.0),
Parameters() : aTerms(), bTerms(),
wavelength(0.0), tolerance(0.0),
units("Degrees"), phaUnit(1.0), band(0), phoStd(0.0),
algoVersion(2019),
iProfile(-1) { }
~Parameters() { }
bool IsValid() const { return (iProfile != -1);}
double a0, a1, a2, a3; //<! Equation parameters
double wavelength; //<! Wavelength for correction
double tolerance; //<! Wavelength Range/Tolerance
QString units; //<! Phase units of equation
double phaUnit; //<! 1 for degrees, Pi/180 for radians
int band; //<! Cube band parameters
double phoStd; //<! Computed photometric std.
int iProfile; //<! Profile index of this data

bool IsValid() const {
return (iProfile != -1);
}
std::vector<double> aTerms; //<! a-terms for 2014 algorithm
std::vector<double> bTerms; //<! b-terms for 2019 algorithm
double wavelength; //<! Wavelength for correction
double tolerance; //<! Wavelength Range/Tolerance
QString units; //<! Phase units of equation
double phaUnit; //<! 1 for degrees, Pi/180 for radians
int band; //<! Cube band parameters
double phoStd; //<! Computed photometric std.
int algoVersion; //<! Algorithm version
int iProfile; //<! Profile index of this data
};

void init(PvlObject &pvl, Cube &cube);
Expand Down
Loading