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

CHANGELOG.md

@@ -34,7 +34,7 @@ release.
 -->
 
 ## [Unreleased]
-
+- Updated the LRO photometry application Lronacpho, to use by default the current 2019 photometric model (LROC_Empirical). The model's coefficients are found in the PVL file that exists in the LRO data/mission/calibration directory. If the old parameter file is provided, the old algorithm(2014) will be used. This functionality is desired for calculation comparisons. Issue: [#4512](https://github.com/USGS-Astrogeology/ISIS3/issues/4512), PR: [#4519](https://github.com/USGS-Astrogeology/ISIS3/pull/4519)
 ### Changed
 
 ### Added

isis/src/lro/apps/lronacpho/LROCEmpirical.cpp

@@ -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>
@@ -38,13 +37,11 @@ namespace Isis {
     init(pvl, cube);
   }
 
-
   /**
    * Destructor
    */
   LROCEmpirical::~LROCEmpirical() {};
 
-
   /**
    * @brief Initialize class from input PVL and Cube files
    *
@@ -120,7 +117,6 @@ namespace Isis {
     return;
   }
 
-
   /**
    * @brief Method to get photometric property given angles
    *
@@ -154,7 +150,6 @@ namespace Isis {
     return (m_bandpho[band - 1].phoStd / ph);
   }
 
-
   /**
    * @brief Performs actual photometric correction calculations
    *
@@ -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
@@ -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) + "]. ";
+      throw IException(IException::Programmer, mess, _FILEINFO_);
+    }
 
     return (rcal);
   }
 
-
   /**
    * @brief Return parameters used for all bands
    *
@@ -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");
@@ -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]));
     }
 
     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.
    *
@@ -316,8 +327,6 @@ namespace Isis {
     return (Parameters());
   }
 
-
-
   /*
    * @brief Extracts necessary LROCEmprical parameters from profile
    *
@@ -333,18 +342,29 @@ 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
+   *
+   *   @history 2022-04-26 Victor Silva - Changed strings casted using toString to 
+   *                                      string literal "0.0"
+   *
    */
   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), "0.0")));
+    for (int i=0; i<7; i++)
+        pars.bTerms.push_back(toDouble(ConfKey(profile, "B" + toString(i), "0.0")));
+
     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", "0"));
 
     return (pars);
   }

isis/src/lro/apps/lronacpho/LROCEmpirical.h

@@ -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 {
@@ -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);