Unverified Commit cae8ec2d authored by Cordell Michaud's avatar Cordell Michaud Committed by GitHub
Browse files

Re-add and refactor lrowacphomap (#4961) 



* Add radiance units label in lrowaccal output cube

* Revert "Add radiance units label in lrowaccal output cube"

* Add refactored lrowacphomap app

* Add lrowacphomap functional tests and test data

* Add entry to CHANGELOG.md for the re-addition and refactoring of lrowacphomap

* Change ShapeModel keyword in Kernels group of lrowacphomap test cube label to relative path

* Add new test cube corresponding to the lrowacphomap test input cube's updated ShapeModel value

* Update expected stats result values for the LrowacphomapNoBack functional test

Co-authored-by: default avatarVictor Silva <victor.silva@asu.edu>
parent db39f2e2

CHANGELOG.md

+1 −0
Original line number Diff line number Diff line
@@ -43,6 +43,7 @@ release. 


### Added

- Improved functionality of msi2isis and MsiCamera model to support new Eros dataset, including support for Gaskell's SUMSPICE files that adjust timing, pointing and spacecraft position ephemeris. [#4886](https://github.com/USGS-Astrogeology/ISIS3/issues/4886)

- Re-added and refactored the LRO photometry application lrowacphomap to be callable for testing purposes. Issue: [#4960](https://github.com/USGS-Astrogeology/ISIS3/issues/4960), PR: [#4961](https://github.com/USGS-Astrogeology/ISIS3/pull/4961)



### Deprecated

isis/src/lro/apps/lrowacphomap/HapkeLRO.cpp

0 → 100644
+352 −0
Original line number Diff line number Diff line 
/**

 * @file

 * $Revision$

 * $Date$

 *

 *   Unless noted otherwise, the portions of Isis written by the USGS are

 *   public domain. See individual third-party library and package descriptions

 *   for intellectual property information, user agreements, and related

 *   information.

 *

 *   Although Isis has been used by the USGS, no warranty, expressed or

 *   implied, is made by the USGS as to the accuracy and functioning of such

 *   software and related material nor shall the fact of distribution

 *   constitute any such warranty, and no responsibility is assumed by the

 *   USGS in connection therewith.

 *

 *   For additional information, launch

 *   $ISISROOT/doc//documents/Disclaimers/Disclaimers.html

 *   in a browser or see the Privacy &amp; Disclaimers page on the Isis website,

 *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on

 *   http://www.usgs.gov/privacy.html.

 */

#include <cmath>

#include <cstdlib>

#include <iomanip>

#include <memory>

#include <iostream>

#include <sstream>



#include <SpiceUsr.h>

#include <SpiceZfc.h>

#include <SpiceZmc.h>



#include "Projection.h"

#include "Camera.h"

#include "HapkeLRO.h"

#include "DbProfile.h"

#include "PvlObject.h"

#include "Brick.h"



using namespace std;



namespace Isis {



  /**

   * Create an HapkeLRO photometric object

   * 

   * @author 2021-08-02 Cordell Michaud

   * 

   * @internal

   *   @history 2021-07-19 Cordell Michaud - Code adapted from PhotometricFunction written by Kris Becker

   * 

   * @param pvl       Photometric parameter files

   * @param cube      Input cube file

   * @param useCamera Enable using camera from input cube

   * @param paramMap  Parameter cube

   */

  HapkeLRO::HapkeLRO(PvlObject &pvl, Cube &cube, bool useCamera, Cube *paramMap)

    : PhotometricFunction(pvl, cube, useCamera) {

    

    m_paramMap = paramMap;

    m_paramProj = m_paramMap->projection();

    m_currentMapSample = 0;

    m_currentMapLine = 0;

    m_currentMapIndex = 0;



    PvlObject phoModel = pvl.findObject("PhotometricModel");

    m_iRef = phoModel.findKeyword("Incref")[0].toDouble(); //  Incidence refernce angle

    m_eRef = phoModel.findKeyword("EmiRef")[0].toDouble(); //  Emission  reference angle

    m_gRef = phoModel.findKeyword("Pharef")[0].toDouble(); //  Phase     reference angle

    if (QString(phoModel.findKeyword("Units")[0]).toUpper() == "DEGREES") {

      m_isDegrees = true;

    }

    else {

      m_isDegrees = false;

    }



    PvlKeyword center = cube.label()->findGroup("BandBin", Pvl::Traverse)["Center"];



    PvlKeyword paramBandNames = m_paramMap->label()->findGroup("BandBin", Pvl::Traverse).findKeyword("Name");



    // Now go through the Groups looking for all the band parameters

    for (int i = 0; i < phoModel.groups(); i++) {

      if (phoModel.group(i).isNamed("Parameters")) {

        PvlGroup paramGroup = phoModel.group(i);

        HapkeLRO::Parameters parms;



        parms.bandBinCenter = paramGroup.findKeyword("BandBinCenter")[0].toDouble();

        for (int j = 0; j < center.size(); j++) {

          if (center[j] == paramGroup.findKeyword("BandBinCenter")[0]) {

            parms.band = j + 1;

          }

        }



        PvlKeyword bands = paramGroup.findKeyword("Bands");

        for (int j = 0; j < bands.size(); j++) {

          parms.mapBands.push_back(bands[j].toInt() - 1);

          parms.names.push_back(QString(paramBandNames[bands[j].toInt() - 1]).toUpper());

          parms.values.push_back(0.0);

        }



        m_bandParameters.push_back(parms);

      }

    }



    m_paramBrickCount = paramMap->lineCount() * paramMap->sampleCount();

    m_paramBricks = new Brick*[m_paramBrickCount];

    memset(m_paramBricks, 0, m_paramBrickCount * sizeof(Brick*));



    m_normalized = true;



    m_hfunc = "HG";

    if (phoModel.hasKeyword("Hfunc")) {

      m_hfunc = phoModel.findKeyword("Hfunc")[0].toUpper();

    }



    if (m_hfunc != "HG") {

      QString msg = "Invalid HFunction: " + m_hfunc;

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

    }

  }



  /**

   * Destructor

   */

  HapkeLRO::~HapkeLRO() {

    for (int i = 0; i < m_paramBrickCount; i++) {

      delete m_paramBricks[i];

    }



    delete [] m_paramBricks;

  }





  bool HapkeLRO::normalized() const {

    return m_normalized;

  }





  /**

   * @brief Return photometric property given angles

   * 

   * This method computes the photometric property at the given 

   * cube location after ensuring a proper parameter container is 

   * found for the specified band.

   * 

   * @author 2021-08-18 Cordell Michaud

   * 

   * @internal

   *   @history 2021-07-19 Cordell Michaud - Adapted code from PhotometricFunction written by Kris Becker

   * 

   * @param i     Incidence angle at cube location

   * @param e     Emission angle at cube location

   * @param g     Phase angle at cube location

   * @param band  Band number in cube (actually is band index) for lookup purposes

   * 

   * @return double Photometric correction using parameters

   */

  double HapkeLRO::photometry(double i, double e, double g, int band) const {

    double lat = m_camera->UniversalLatitude();

    double lon = m_camera->UniversalLongitude();

    

    return photometry(i, e, g, lat, lon, band);

  }





  /**

   * @brief Return photometric property given angles

   * 

   * This method computes the photometric property at the given 

   * cube location after ensuring a proepr parameter container is 

   * found for the specified band.

   * 

   * @author 2021-08-18 Cordell Michaud

   * 

   * @internal

   *   @history 2021-08-03 Cordell Michaud - Adapted code from PhotometricFunction written by Kris Becker

   * 

   * @param i     Incidence angle at cube location

   * @param e     Emission angle at cube location

   * @param g     Phase angle at cube location

   * @param lat   Latitude at cube location

   * @param lon   Longitude at cube location

   * @param band  Band number in cube (actually is band index) for lookup purposes

   * 

   * @return double Photometric correction using parameters

   */

  double HapkeLRO::photometry(double i, double e, double g, double lat, double lon, int band) const {

    m_paramProj->SetUniversalGround(lat, lon);

    int intSamp = (int) (m_paramProj->WorldX() + 0.5);

    int intLine = (int) (m_paramProj->WorldY() + 0.5);



    if (intSamp <= 0 || intSamp > m_paramMap->sampleCount()

        || intLine <= 0 || intLine > m_paramMap->lineCount()) {

      

      return 0.0;

    }



    // If the line and sample of the parameter map are different

    //   we need to update the parameters

    if (m_currentMapSample != intSamp || m_currentMapLine != intLine) {

      m_currentMapSample = intSamp;

      m_currentMapLine = intLine;

      m_currentMapIndex = 0;

      // read the data for that line/sample into a brick buffer

      int paramBrickIndex = (intLine - 1) * m_paramMap->sampleCount() + intSamp;



      // if this brick doesn't exist, create it

      if (m_paramBricks[paramBrickIndex] == NULL) {

        Brick *newBrick = new Brick(1, 1, m_paramMap->bandCount(), m_paramMap->pixelType());

        newBrick->SetBasePosition(intSamp, intLine, 1);

        m_paramMap->read(*newBrick);



        // store the Brick, so we don't have to read it from the file again

        m_paramBricks[paramBrickIndex] = newBrick;

      }

      Brick &b = *m_paramBricks[paramBrickIndex];



      // And for each parameter, copy the appropriate value

      for (unsigned int p = 0; p < m_bandParameters.size(); p++) {

        HapkeLRO::Parameters &parms = m_bandParameters[p];



        if (parms.band == band) {

          m_currentMapIndex = p;

          for (unsigned int v = 0; v < parms.values.size(); v++) {

            parms.values[v] = b[parms.mapBands[v]];

            if (IsSpecial(parms.values[v])) {

              return Isis::Null;

            }

          }

          parms.phoStd = photometry(parms, m_iRef, m_eRef, m_gRef);

        }

      }

    }



    double ph = photometry(m_bandParameters[m_currentMapIndex], i, e, g);

    if (!normalized()) {

      return ph;

    }



    return m_bandParameters[m_currentMapIndex].phoStd / ph;

  }





  /**

   * @brief Performs actual photometric correction calculations

   *

   * This routine computes photometric correction using parameters

   * for the Exponential-Buratti-Hill equation.

   *

   * @author 2021-07-29 Cordell Michaud

   * 

   * @internal

   *  @history 2021-07-19 Cordell Michaud - Code adapted from HapkeLRO written by Kris Becker

   *

   * @param parms Container of band-specific HapkeLRO parameters

   * @param i     Incidence angle in degrees

   * @param e     Emission angle in degrees

   * @param g     Phase angle in degrees

   *

   * @return double Photometric correction parameter

   */

  double HapkeLRO::photometry(HapkeLRO::Parameters &parms, double i, double e, double g) const {

    double w = parms["W"];

    double bc0 = parms["BCO"];

    double hc = parms["HC"];

    double xi = parms["XI"];



    //  Ensure problematic values are adjusted

    if (i == 0.0) {

      i = 10.E-12;

    }

    if (e == 0.0) {

      e = 10.E-12;

    }



    // Convert to radians

    i *= rpd_c();

    e *= rpd_c();



    double cosg = 0.0;

    double tang2hc = 0.0;

    if (!m_isDegrees) {

      g *= rpd_c();

      cosg = cos(g);

      tang2hc = tan(g / 2) / hc;

    }

    else {

      cosg = cos(rpd_c() * g);

      tang2hc = tan(rpd_c() * g / 2) / hc;

    }



    // Compute Lommel-Seeliger components

    double mu = cos(e);

    double mu0 = cos(i);



    double p = (1 - pow(xi, 2)) / pow(1 - 2 * xi * cosg + pow(xi, 2), 1.5);

    double Hmu = (1 + 2 * mu) / (1 + 2 * mu * pow(1 - w, 0.5));

    double Hmu0 = (1 + 2 * mu0) / (1 + 2 * mu0 * pow(1 - w, 0.5));

    double Bc = (1 + (1 - exp(-1 * tang2hc)) / (tang2hc)) / (2 * pow(1 + tang2hc, 2));



    double rcal = w / 4 * (p + Hmu0 * Hmu - 1) * (1 + bc0 * Bc);

    

    return rcal * mu0 / (mu0 + mu);

  }







  /**

   * @brief Return parameters used for all bands

   * 

   * This method creates keyword vectors of band-specific parameters used 

   * in the photometric correction.

   * 

   * @author 2021-08-03 Cordell Michaud

   * 

   * @internal

   *   @history 2021-07-19 Cordell Michaud - Adapted code from HapkeLRO written by Kris Becker

   *   @history 2021-07-19 Cordell Michaud - Adapted code from PhotometricFunction written by Kris Becker

   * 

   * @param pvl Output PVL container write keywords

   */

  void HapkeLRO::report(PvlContainer &pvl) {

    pvl.addComment(QString("IoF/LS = w/4 * (p(g) + H(mu0,w)*H(mu,w)-1) * (1+Bc0*Bc(g,h))"));

    pvl.addComment(QString("  where:"));

    pvl.addComment(QString("    p(g) = (1-xi^2)/(1-2*xi*cos(g) + xi^2)^(3/2)"));

    pvl.addComment(QString("    H(x,w) = (1+2*x)/(1+2*x*sqrt(1-w))"));

    pvl.addComment(QString("    Bc(g,h) = (1 + (1-exp(-tan(g/2)/h))/(tan(g/2)/h))/(2*(1+tan(g/2)/h)^2)"));



    pvl += PvlKeyword("Algorithm", "HapkeLRO");

    pvl += PvlKeyword("ParameterMapCube", m_paramMap->fileName());

    pvl += PvlKeyword("IncRef", toString(m_iRef), "degrees");

    pvl += PvlKeyword("EmiRef", toString(m_eRef), "degrees");

    pvl += PvlKeyword("PhaRef", toString(m_gRef), "degrees");

    PvlKeyword units("Units");

    

    if (m_isDegrees) {

      units += "Degrees";

    }

    else {

      units += "Radians";

    }

    

    return;

  }





  void HapkeLRO::setNormalized(bool normalized) {

    m_normalized = normalized;

  }



} // namespace Isis

isis/src/lro/apps/lrowacphomap/HapkeLRO.h

0 → 100644
+123 −0
Original line number Diff line number Diff line 
#ifndef hapkelro_h

#define hapkelro_h

/**

 * @file

 * $Revision$

 * $Date$

 *

 *   Unless noted otherwise, the portions of Isis written by the USGS are

 *   public domain. See individual third-party library and package descriptions

 *   for intellectual property information, user agreements, and related

 *   information.

 *

 *   Although Isis has been used by the USGS, no warranty, expressed or

 *   implied, is made by the USGS as to the accuracy and functioning of such

 *   software and related material nor shall the fact of distribution

 *   constitute any such warranty, and no responsibility is assumed by the

 *   USGS in connection therewith.

 *

 *   For additional information, launch

 *   $ISISROOT/doc//documents/Disclaimers/Disclaimers.html

 *   in a browser or see the Privacy &amp; Disclaimers page on the Isis website,

 *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on

 *   http://www.usgs.gov/privacy.html.

 */



#include <iostream>

#include <sstream>

#include <iomanip>

#include <vector>



#include "PhotometricFunction.h"

#include "Camera.h"

#include "SpecialPixel.h"



namespace Isis {



  class PvlObject;

  class Camera;

  class Brick;

  class Projection;





  /**

   * @brief An implementation of the HapkeLRO photometric function

   *

   * This class implements the Exponential-Buratti-Hill photometric

   * equation as outline in thier paper "Multispectral Photometry

   * of the Moon and Absolute Calibration of the Clementine UV/VIS

   * Camera", published in Icaris v141, pg. 205-255 (1999).

   *

   * @author  2021-08-18 Cordell Michaud

   * 

   * @internal

   *   @history 2021-07-19 Cordell Michaud - Code adapted from PhotometricFunction and HapkeLRO written by Kris Becker

   *

   */

  class HapkeLRO : public PhotometricFunction {

    public:

      HapkeLRO(PvlObject &pvl, Cube &cube, bool useCamera, Cube *paramMap);

      ~HapkeLRO();

      bool normalized() const;

      double photometry(double i, double e, double g, int band = 1) const override;

      double photometry(double i, double e, double g, double lat, double lon, int band = 1) const;

      void report(PvlContainer &pvl);

      void setNormalized(bool normalized);



    private:

      QString m_hfunc;



      bool m_isDegrees;



      bool m_normalized;



      Cube *m_paramMap;

      mutable Projection *m_paramProj;

      mutable int m_currentMapSample, m_currentMapLine, m_currentMapIndex;



      mutable Brick **m_paramBricks;

      int m_paramBrickCount;





      /**

       * Container for band photometric correction parameters

       * 

       * @author 2021-07-30 Cordell Michaud

       * 

       * @internal

       *   @history 2021-07-19 Cordell Michaud - Code adapted from PhotoemtricFunction written by Kris Becker

       */

      class Parameters

      {

        public:

          Parameters() 

            : band(1), bandBinCenter(0.0), mapBands(), 

            names(), phoStd(0.0), values() {}



          ~Parameters() {}



          double operator[] (const QString &name) const {

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

              if (names[i] == name) {

                return values[i];

              }

            }



            QString msg = "No parameter band named '" + name + "'.";

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

          }



          int band; // band number

          double bandBinCenter; // center wavelength

          std::vector<int> mapBands; // parameter bands

          std::vector<QString> names; // parameter names

          double phoStd; // computed photometric std

          std::vector<double> values; // parameter values

      };



      mutable std::vector<Parameters> m_bandParameters; // Parameters for each band

      double photometry(Parameters &parms, double i, double e, double g) const;

  };

};



#endif
 
 No newline at end of file

isis/src/lro/apps/lrowacphomap/HapkeLROC.cpp

0 → 100644
+638 −0

File added.

Preview size limit exceeded, changes collapsed.

isis/src/lro/apps/lrowacphomap/HapkeLROC.h

0 → 100644
+131 −0
Original line number Diff line number Diff line 
#ifndef hapkelroc_h

#define hapkelroc_h

/**

 * @file

 * $Revision$

 * $Date$

 *

 *   Unless noted otherwise, the portions of Isis written by the USGS are

 *   public domain. See individual third-party library and package descriptions

 *   for intellectual property information, user agreements, and related

 *   information.

 *

 *   Although Isis has been used by the USGS, no warranty, expressed or

 *   implied, is made by the USGS as to the accuracy and functioning of such

 *   software and related material nor shall the fact of distribution

 *   constitute any such warranty, and no responsibility is assumed by the

 *   USGS in connection therewith.

 *

 *   For additional information, launch

 *   $ISISROOT/doc//documents/Disclaimers/Disclaimers.html

 *   in a browser or see the Privacy &amp; Disclaimers page on the Isis website,

 *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on

 *   http://www.usgs.gov/privacy.html.

 */



#include "PhotometricFunction.h"



namespace Isis {

  class PvlObject;

  class Camera;

  class Brick;

  class Projection;





  /**

   * @brief An implementation of the bidirectional reflectance model by Hapke (2012).

   *

   * This class is based on the Hapke PhotoModel from Isis base code. Photometry method

   * was updated to match full Hapke (2012) model and some variable names were changed

   * to match reference formula.

   * 

   * Reference:

   *    Hapke, B. (2012) Theory of Reflectance and Emittance Spectroscopy, Cambridge Univ. Press.

   *

   * @author 2021-08-18 Cordell Michaud

   * @internal

   *   @history 2021-07-19 Cordell Michaud - Code adapted from HapkeLROC written by Dan Clarke

   *   @history 2021-07-19 Cordell Michaud - Code adapted from PhotometricFunction written by Kris Becker

   *

   */

  class HapkeLROC : public PhotometricFunction {

    public:

      HapkeLROC(PvlObject &pvl, Cube &cube, bool useCamera, Cube *paramMap);

      ~HapkeLROC();



      bool normalized() const;

      double photometry(double i, double e, double g, int band = 1) const override;

      double photometry(double i, double e, double g, double lat, double lon, int band = 1) const;

      void report(PvlContainer &pvl);

      void setNormalized(bool normalized);



    private:

      QString m_hfunc;



      bool m_isDegrees;



      bool m_normalized;



      Cube *m_paramMap;

      mutable Projection *m_paramProj;

      mutable int m_currentMapSample, m_currentMapLine, m_currentMapIndex;



      mutable Brick **m_paramBricks;

      int m_paramBrickCount;



      mutable double m_chiThetaP;

      mutable double m_invChiThetaP;

      mutable double m_oldPhase;

      mutable double m_oldIncidence;

      mutable double m_oldEmission;

      mutable double m_photoB0save;

      mutable double m_photoCot2t;

      mutable double m_photoCott;

      mutable double m_photoTant;

      mutable double m_photoThetaold;



      mutable double m_result;





      /**

        * Container for band photometric correction parameters

        * 

        * @author 2021-08-02 Cordell Michaud

        * 

        * @internal

        *   @history 2021-07-19 Cordell Michaud - Code adapted from PhotometricFunction written by Kris Becker

        */

      class Parameters {

        public:

          Parameters() 

            : band(1), bandBinCenter(0.0), mapBands(), 

            names(), phoStd(0.0), values() {}



          ~Parameters() {}



          double operator[] (const QString &name) const {

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

              if (names[i] == name) {

                return values[i];

              }

            }



            QString msg = "No parameter band named '" + name + "'.";

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

          }



          int band; // band number

          double bandBinCenter; // center wavelength

          std::vector<int> mapBands; // parameter bands

          std::vector<QString> names; // parameter names

          double phoStd; // computed photometric std

          std::vector<double> values; // parameter values

      };



      mutable std::vector<Parameters> m_bandParameters; // Parameters for each band

      double photometry(const Parameters &parms, double i, double e, double g) const;

  };



};



#endif
 
 No newline at end of file