Unverified Commit c9f5f8d6 authored by Victor Silva's avatar Victor Silva Committed by GitHub
Browse files

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

* updated to use new lroc empirical model and parameters and ability to get from calibration file but still backwards compatible if old parameters file exists will use old model

* refactored app to be callable and created gtests

* removed unnecessary new lines and added items to lronacpho data directory for gtest

* cropped adding test cub

* fixed issues from code review and added history comments to functional test
parent e5a0a9d8

CHANGELOG.md

+1 −1
Original line number Diff line number Diff line
@@ -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

+62 −42
Original line number Diff line number Diff line
@@ -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,6 +216,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

   */

  void LROCEmpirical::report( PvlContainer &pvl ) {
@@ -219,16 +225,18 @@ namespace Isis { 
    pvl.addComment(" where:");

    pvl.addComment("  mu0 = cos(incidence)");

    pvl.addComment("  mu = cos(emission)");



    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,10 +246,13 @@ 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];
@@ -250,10 +261,10 @@ namespace Isis { 
        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));

        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;
@@ -261,15 +272,15 @@ namespace Isis { 
    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

   *
@@ -277,7 +288,7 @@ namespace Isis { 
   * 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

+19 −11
Original line number Diff line number Diff line
@@ -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

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

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

0 → 100644
+222 −0
Original line number Diff line number Diff line 
#include "ProcessByLine.h"

#include "LROCEmpirical.h"

#include "PhotometricFunction.h"

#include "Buffer.h"

#include "Camera.h"

#include "iTime.h"

#include "SpecialPixel.h"

#include "Pvl.h"

#include "PvlGroup.h"

#include "UserInterface.h"

#include "lronacpho.h"



using namespace std;



namespace Isis{



  // globals

  static bool g_useDEM;

  PhotometricFunction *g_phoFunction;



  // forward declare helper functions

  void phoCal (Buffer &in, Buffer &out);

  void phoCalWithBackplane (vector<Isis::Buffer *> &in, vector<Isis::Buffer *> &out );



  /**

    * @brief Photometric application for the LRO NAC cameras

    *

    * This application provides features that allow multiband cubes for LRO NAC cameras

    * to be photometrically correctedConstructor

    *

    * @param ui The user interfact to parse the parameters from. 

    */

  void lronacpho(UserInterface &ui, Pvl *log){

    Cube iCube(ui.GetCubeName("FROM"));

    lronacpho(&iCube, ui, log);

  }



  /**

    * @brief Photometric application for the LRO NAC cameras

    *

    * This application provides features that allow multiband cubes for LRO NAC cameras

    * to be photometrically corrected

    * @author 2016-09-16 Victor Silva

    *

    * @internal

    *  @history 2016-09-19 Victor Silva - Adapted from lrowacpho written by Kris Becker

    *	 @history 2021-03-12 Victor Silva - Updates"PostCall" include ability to run with default values

    * 																			Added new values for 2019 version of LROC Empirical function.

    *  @history 2022-04-18 Victor Silva - Refactored to make callable for GTest framework                                   

    *

    * @param iCube The input cube to be photometrically corrected.

    * @param ui The user interfact to parse the parameters from. 

    */

  void lronacpho(Cube *iCube, UserInterface &ui, Pvl *log){

    ProcessByLine p;

    p.SetInputCube(iCube);



    Cube *oCube = p.SetOutputCube(ui.GetCubeName("TO"), ui.GetOutputAttribute("TO"));//, 

    

    // Backplane option

    bool useBackplane = false;

    if (ui.WasEntered("BACKPLANE")) {

      CubeAttributeInput backplaneCai = ui.GetInputAttribute("BACKPLANE");



      Cube bpCube;

      bpCube.open(ui.GetFileName("BACKPLANE"));

      int bpBands = bpCube.bandCount();

      bpCube.close();



      int bpCaiBands = backplaneCai.bands().size();

      if (bpBands < 3 || (bpCaiBands != 3 && bpBands > 3)) {

        string msg = "Invalid Backplane: The backplane must be exactly 3 bands";

        throw IException(IException::User, msg, _FILEINFO_);

      }



      if (iCube->bandCount() != 1) {

        string msg = "Invalid Image: The backplane option can only be used with a single image band at a time.";

        throw IException(IException::User, msg, _FILEINFO_);

      }



      CubeAttributeInput cai;

      bpCaiBands == 3 ? cai.setAttributes("+" + backplaneCai.bands()[0]) : cai.setAttributes("+1" ) ;

      p.SetInputCube(ui.GetFileName("BACKPLANE"), cai);

      bpCaiBands == 3 ? cai.setAttributes("+" + backplaneCai.bands()[1]) : cai.setAttributes("+2" ) ;

      p.SetInputCube(ui.GetFileName("BACKPLANE"), cai);

      bpCaiBands == 3 ? cai.setAttributes("+" + backplaneCai.bands()[2]) : cai.setAttributes("+3" ) ;

      p.SetInputCube(ui.GetFileName("BACKPLANE"), cai);



      useBackplane = true;

    }

    // Get name of parameters file

    QString algoName = "";

    QString algoFile = ui.GetAsString("PHOPAR");



    FileName algoFileName(algoFile);



    if(algoFileName.isVersioned())

      algoFileName = algoFileName.highestVersion();

    

    if(!algoFileName.fileExists()) {

      QString msg = algoFile + " does not exist.";

      throw IException(IException::User, msg, _FILEINFO_);

    }



    Pvl params(algoFileName.expanded());



    algoName = PhotometricFunction::algorithmName(params);

    algoName = algoName.toUpper();



    // Set NAC algorithm

    if (algoName == "LROC_EMPIRICAL") {

        g_phoFunction = new LROCEmpirical(params, *iCube, !useBackplane);

    }

    else {

      QString msg = " Algorithm Name [" + algoName + "] not recognized. ";

      msg += "Compatible Algorithms are:\n LROC_Empirical\n";

      throw IException(IException::User, msg, _FILEINFO_);

    }



    // Set user selected max and mins

    g_phoFunction->setMinimumPhaseAngle(ui.GetDouble("MINPHASE"));

    g_phoFunction->setMaximumPhaseAngle(ui.GetDouble("MAXPHASE"));

    g_phoFunction->setMinimumEmissionAngle(ui.GetDouble("MINEMISSION"));

    g_phoFunction->setMaximumEmissionAngle(ui.GetDouble("MAXEMISSION"));

    g_phoFunction->setMinimumIncidenceAngle(ui.GetDouble("MININCIDENCE"));

    g_phoFunction->setMaximumIncidenceAngle(ui.GetDouble("MAXINCIDENCE"));



    // Set use of DEM to calculate photometric angles

    g_useDEM = ui.GetBoolean("USEDEM");



    // Begin processing by line

    if(useBackplane) {

      p.StartProcess(phoCalWithBackplane);

    }

    else {

      p.StartProcess(phoCal);

    }

    // Start all the PVL

    PvlGroup photo("Photometry");

    g_phoFunction->report(photo);

        

    oCube->putGroup(photo);

    

    if(log){

      log->addGroup(photo);

    }

    //Application::Log(photo);



    p.EndProcess();

    p.ClearInputCubes();

    delete g_phoFunction;

  }

  

  /**

    * @brief Apply LROC Empirical photometric correction

    *

    * Short function dispatched for each line to apply the LROC Empirical photometric

    * correction function.

    *

    * @author 2016-09-19 Victor Silva

    *

    * @internal

    *   @history 2016-09-19 Victor silva - Adapted from lrowacpho written by Kris Becker

    *   @history 2022-04-18 Victor Silva - Refactored to make callable for GTest framework

    *

    * @param in Buffer containing input data

    * @param out Buffer of photometrically corrected data

    */

  void phoCal(Buffer &in, Buffer &out) {

    // Iterate through pixels

    for(int i = 0; i < in.size(); i++){

      // Ignore special pixels

      if(IsSpecial(in[i])){

        out[i] = in[i];

      }

      else{

        // Get correction and test for validity

        double ph = g_phoFunction->compute(in.Line(i), in.Sample(i), in.Band(i), g_useDEM);

        out[i] = ( IsSpecial(ph) ? Null : (in[i] * ph) );

      }

    }

    return;

  }//end phoCal



  /**

    * @brief Apply LROC Empirical photometric correction with backplane

    *

    * Short function dispatched for each line to apply the LROC Empirical photometrc

    * correction function.

    *

    * @author 2016-09-19 Victor Silva

    *

    * @internal

    *   @history 2016-09-19 Victor silva - Adapted from lrowacpho written by Kris Becker

    *   @history 2022-04-18 Victor Silva - Refactored to make callable for GTest framework

    *

    * @param in Buffer containing input data

    * @param out Buffer of photometrically corrected data

    */

  void phoCalWithBackplane( std::vector<Isis::Buffer *> &in, std::vector<Isis::Buffer *> &out ){



    Buffer &image = *in[0];

    Buffer &phase = *in[1];

    Buffer &emission = *in[2];

    Buffer &incidence = *in[3];

    Buffer &calibrated = *out[0];



    for (int i = 0; i < image.size(); i++) {

      //  Don't correct special pixels

      if (IsSpecial(image[i])) {

        calibrated[i] = image[i];

      }

      else {

      // Get correction and test for validity

        double ph = g_phoFunction->photometry(incidence[i], emission[i], phase[i], image.Band(i));

        calibrated[i] = (IsSpecial(ph) ? Null : image[i] * ph);

      }

    }

    return;

  }

}

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

0 → 100644
+13 −0
Original line number Diff line number Diff line 
#ifndef lronacpho_h

#define lronacpho_h



#include "Cube.h"

#include "lronacpho.h"

#include "UserInterface.h"



namespace Isis{

  extern void lronacpho(Cube *iCube, UserInterface &ui, Pvl *log=nullptr);

  extern void lronacpho(UserInterface &ui, Pvl *log=nullptr);

}



#endif
 
 No newline at end of file