Added lronacpho application for photometric calibration of Lunar Reconnaissance Orbiter NAC images.

#include <cmath>

#include <cstdlib>

#include <iomanip>

#include <iostream>

#include <memory>

#include <sstream>

#include <SpiceUsr.h>

#include <SpiceZfc.h>

#include <SpiceZmc.h>

#include "Camera.h"

#include "DbProfile.h"

#include "LROCEmpirical.h"

#include "PvlObject.h"

using namespace std;

namespace Isis {

 /**

  * Create an LROCEmpirical photometric object

  * 

  * @author 2016-08-15 Victor Silva

  *

  * @internal

  *   @history 2016-08-15 - Code adapted from lrowacpho written by Kris Becker

  */

  LROCEmpirical::LROCEmpirical(PvlObject &pvl, Cube &cube, bool useCamera) :

                 PhotometricFunction(pvl, cube, useCamera) {

    init(pvl, cube);

  }

  /**

   * Destructor

   */

  LROCEmpirical::~LROCEmpirical() {};

  /**

   * @brief Initialize class from input PVL and Cube files

   *

   * This method is typically called at the class instantiation

   * time but is reentrant. It reads the parameter PVL file and 

   * extracts Photometric model and Normalization models from it.

   * The cube is needed to match all potential profiles for each

   * band.

   *

   * @param pvl Input PVL parameter files

   * @param cube Input cube file to correct

   * 

   * @throws IException::User "Errors in the input PVL file."

   * @throws IException::User "Band with wavelength Center does not have

   *                           PhotometricModel Algorithm group/profile."

   *

   * @author Victor Silva 2016-08-15

   *

   * @internal

   *   @history 2016-08-15 Victor Silva - Adapted from the lrowacpho application

   *                           written by Kris Becker.

   */

  void LROCEmpirical::init( PvlObject &pvl, Cube &cube ) {

    //  Make it reentrant

    m_profiles.clear();

    m_bandpho.clear();

    //  Interate over all Photometric groups

    m_normProf = DbProfile(pvl.findObject("NormalizationModel").findGroup("Algorithm", Pvl::Traverse));

    m_iRef = toDouble(ConfKey(m_normProf, "IncRef", toString(30.0)));

    m_eRef = toDouble(ConfKey(m_normProf, "EmaRef", toString(0.0)));

    m_gRef = toDouble(ConfKey(m_normProf, "PhaRef", toString(m_iRef)));

    PvlObject &phoObj = pvl.findObject("PhotometricModel");

    DbProfile phoProf = DbProfile(phoObj);

    PvlObject::PvlGroupIterator algo = phoObj.beginGroup();

    while (algo != phoObj.endGroup()) {

      if (algo->name().toLower() == "algorithm") {

        m_profiles.push_back(DbProfile(phoProf, DbProfile(*algo)));

      }

      ++algo;

    }

    Pvl *label = cube.label();

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

    QString errs("");

    for (int i = 0; i < cube.bandCount(); i++) {

      Parameters parms = findParameters(toDouble(center[i]));

      if (parms.IsValid()) {

        parms.band = i + 1;

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

        m_bandpho.push_back(parms);

      }

      else { // Appropriate photometric parameters not found

        ostringstream mess;

        mess << "Band [" << i + 1 << "] with wavelength Center = [" << center[i] << 

        "] does not have PhotometricModel Algorithm group/profile";

        IException e(IException::User, mess.str(), _FILEINFO_);

        errs += e.toString() + "\n";

      }

    }

    // Check for errors and throw them all at the same time

    if (!errs.isEmpty()) {

      errs += " --> Errors in the input PVL file [" + pvl.fileName() + "]";

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

    }

    return;

  }

  /**

   * @brief Method to get photometric property given angles

   *

   * This method computes photometric property at the given cube location after 

   * proper parameter container is found for the specific band.

   * 

   * @param i    Incidence angle in degrees

   * @param e    Emission angle in degrees

   * @param g    Phase angle in degrees

   * @param band cube band to be tested

   *

   * 

   * @throws IException::Programmer "Provided band out of range."

   *

   * @return @b double The photometric property

   *

   * @author 2016-08-15 Victor Silva

   *

   * @internal

   *   @history 2016-08-15 Victor Silva - Adapted code from lrowacpho application

   *                         written by Kris Becker

   */

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

    if ((band <= 0) || (band > (int) m_bandpho.size())) {

      std::string mess = "Provided band " + IString(band) + " out of range.";

      throw IException(IException::Programmer, mess, _FILEINFO_);

    }

    double ph = photometry(m_bandpho[band - 1], i, e, g);

    return (m_bandpho[band - 1].phoStd / ph);

  }

  /**

   * @brief Performs actual photometric correction calculations

   *

   * This routine computes photometric correction using parameters for the 

   * LROC Emperical equation.

   *

   * @param parms Container of band-specific exponential parameters

   * @param i     Incidence angle in degrees

   * @param e     Emission angle in degrees

   * @param g     Phase angle in degrees

   *

   * @return @b double Photometric correction parameter

   *

   * @author 2016-08-15 Victor Silva

   *

   * @internal

   *   @history 2016-08-15 Victor Silva - Adapted code from lrowacpho application

   *                         written by Kris Becker

   *

   */

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

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

    g *= parms.phaUnit; //  Apply unit normalizer

    // Compute Lommel-Seeliger components

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

    return (rcal);

  }

  /**

   * @brief Return parameters used for all bands

   * 

   * Method creates keyword vector of band-specific parameters

   * used in the photometric correction.

   *

   * @param pvl Output PVL container for keywords written

   *

   * @author 2016-08-15 Victor Silva 

   *

   * @internal

   *   @history 2016-08-15 Victor Silva - Adapted code from lrowacpho application

   *                         written by Kris Becker

   */

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

    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("IncRef", toString(m_iRef), "degrees");

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

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

    PvlKeyword units("FunctionUnits");

    PvlKeyword phostd("PhotometricStandard");

    PvlKeyword bbc("BandBinCenter");

    PvlKeyword bbct("BandBinCenterTolerance");

    PvlKeyword bbn("BandNumber");

    PvlKeyword a0("A0");

    PvlKeyword a1("A1");

    PvlKeyword a2("A2");

    PvlKeyword a3("A3");

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

    }

    pvl += units;

    pvl += phostd;

    pvl += bbc;

    pvl += bbct;

    pvl += bbn;

    pvl += a0;

    pvl += a1;

    pvl += a2;

    pvl += a3;

    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 

   * parameter and if it is less than or equal to it, a Parameter 

   * container is returned. 

   *

   * @param wavelength The wavelength to find parameters for

   *

   * @return @b LROEmpirical::parameters The parameters used for the specified

   *                                     wavelength. If a BandBinCenter is not

   *                                     found for the wavelength then a

   *                                     default parameters object is returned.

   *

   * @author 2016-08-15 Victor Silva 

   * 

   * @internal

   *   @history 2016-08-15 Victor Silva - Adapted from the lrowacpho application

   *                           written by Kris Becker.

   */

  LROCEmpirical::Parameters LROCEmpirical::findParameters( const double wavelength ) const {

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

      const DbProfile &profile = m_profiles[i];

      if (profile.exists("BandBinCenter")) {

        double p_center = toDouble(ConfKey(profile, "BandBinCenter", toString(Null)));

        double tolerance = toDouble(ConfKey(profile, "BandBinCenterTolerance", toString(1.0E-6)));

        if (fabs(wavelength - p_center) <= fabs(tolerance)) {

          Parameters pars = extract(profile);

          pars.iProfile = i;

          pars.wavelength = wavelength;

          pars.tolerance = tolerance;

          return (pars);

        }

      }

    }

    // Not found if we reach here

    return (Parameters());

  }

  /*

   * @brief Extracts necessary LROCEmprical parameters from profile

   *

   * Given a profile read from the input PVL file, this method extracts

   * needed parameters (from keywords) in the PVL profile and creates

   * a container of the converted values

   *

   * @param profile The profile to extract parameters for

   *

   * @return @b LROCEmpirical::parameters

   *

   * @author Victor Silva 2016-08-15

   *

   * @internal

   *   @history 2016-08-15 Victor Silva - Adapted from the lrowacpho application

                              written by Kris Becker.

   */

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

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

    return (pars);

  }

}

#ifndef LROCEmpirical_h

#define LROCEmpirical_h

#include <iostream>

#include <iomanip>

#include <sstream>

#include <iomanip>

#include "Camera.h"

#include "DbProfile.h"

#include "IString.h"

#include "PhotometricFunction.h"

#include "SpecialPixel.h"

namespace Isis {

  class PvlObject;

  class Camera;

  /**

   * @brief An implementation of photometric equation used for LRO NAC Cameras

   *

   * This class implements the LROC Empirical photometric equation used for the LRO NAC

   * Cameras.

   *

   * @author 2016-08-15 Victor Silva

   *

   * @internal

   *   @history 2016-08-15 - Code adapted from lrowacpho written by Kris Becker

   *

   */

  class LROCEmpirical : public PhotometricFunction {

    public:

      LROCEmpirical (PvlObject &pvl, Cube &cube, bool useCamera);

      virtual ~LROCEmpirical () ;

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

      void report( PvlContainer &pvl );

    private:

      /**

       * Container for band photometric correction parameters

       *

       * @author 2016-08-15 Victor Silva

       * 

       * @internal

       *   @history 2016-08-05 - Code adapted from lrowacpho written by Kris Becker

       */

      struct Parameters {

        Parameters() : a0(0.0), a1(0.0), a2(0.0), a3(0.0),

                       wavelength(0.0), tolerance(0.0),

                       units("Degrees"), phaUnit(1.0), band(0), phoStd(0.0),

                       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

      };

      void init(PvlObject &pvl, Cube &cube);

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

      Parameters findParameters(const double wavelength) const;

      Parameters extract( const DbProfile &profile) const;

      std::vector<DbProfile> m_profiles; //<! Profiles for all possible wavelengths

      std::vector<Parameters> m_bandpho; //<! Parameters for each band

  };

};

#endif

ifeq ($(ISISROOT), $(BLANK))

.SILENT:

error:

	echo "Please set ISISROOT";

else

	include $(ISISROOT)/make/isismake.apps

endif

 No newline at end of file

#include <SpiceUsr.h>

#include <SpiceZfc.h>

#include <SpiceZmc.h>

#include "Angle.h"

#include "Camera.h"

#include "DbProfile.h"

#include "PhotometricFunction.h"

#include "PvlObject.h"

using namespace std;

namespace Isis {

 /**

  * Construct Photometric function from Pvl and Cube file

  * 

  * @param pvl photometric parameter files

  * @param cube Input cube file

  *

  * @author 2016-08-15 Victor Silva

  *

  * @internal

  *   @history 2016-08-15 Code adapted from lrowacpho written by Kris Becker

  *

  **/

  PhotometricFunction::PhotometricFunction( PvlObject &pvl, Cube &cube , bool useCamera ) {

    if (useCamera) {

      m_camera = cube.camera();

    }

  }

  /**

   * Destructor

   */

  PhotometricFunction::~PhotometricFunction() {}

  /**

   * Set the camera used to compute photometric angles.

   * 

   * @param cam A pointer to the camera to be used

   */

  void PhotometricFunction::setCamera(Camera *cam) {

    m_camera = cam;

  }

  /**

   * Finds the name of the algorithm defined in a PVL object.

   * 

   * @param pvl The pvl to find the algorithm name in.

   * 

   * @return @b QString The algorithm name from the PVL object.

   */

  QString PhotometricFunction::algorithmName( const PvlObject &pvl ) {

    return pvl.findObject("PhotometricModel").findGroup("Algorithm", Pvl::Traverse).findKeyword("Name")[0];

  }

 /**

  * Computes Photometric function from cube attributes

  * 

  * @param line line number in cube

  * @param sample sample number in cube

  * @param band band number in cube

  * @param useDem boolean to use provided Dem

  * 

  * @return @b double photometry calculations

  *

  * @author 2016-08-15 Victor Silva

  *

  * @internal

  *   @history 2016-08-15 Victor Silva - Code adapted from lrowacpho written by Kris Becker

  *

  **/

  double PhotometricFunction::compute( const double &line, const double &sample, int band, bool useDem) {

    // Update band if necessary

    if (m_camera->Band() != band) {

       m_camera->SetBand(band);

    }

    // Return null if not able to set image

    if (!m_camera->SetImage(sample, line)) {

      return (Null);

    }

    // calculate photometric angles

    double i = m_camera->IncidenceAngle();

    double e = m_camera->EmissionAngle();

    double g = m_camera->PhaseAngle();

    bool success = true;

    if (useDem) {

      Angle phase, incidence, emission;

      m_camera->LocalPhotometricAngles(phase, incidence, emission, success);

      if (success) {

        g = phase.degrees();

        i = incidence.degrees();

        e = emission.degrees();

      }

    }

    if ( !success || i < minimumIncidenceAngle() || i > maximumIncidenceAngle() || e < minimumEmissionAngle() || e

                        > maximumEmissionAngle() || g < minimumPhaseAngle() || g > maximumPhaseAngle()) {

      return (Null);

    }

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

  }

 /**

  * Mutator function to set minimum incidence angle

  *

  * @param double The new minimum incidence angle

  *

  * @author 2016-08-15 Victor Silva

  *

  * @internal

  *   @history 2016-08-15 Code adapted from lrowacpho written by Kris Becker

  *

  **/

  void PhotometricFunction::setMinimumIncidenceAngle(double angle) {

    m_minimumIncidenceAngle = angle;

  }

 /**

  * Mutator function to set maximum incidence angle

  *

  * @param double the new maximum incidence angle

  *

  * @author 2016-08-15 Victor Silva

  *

  * @internal

  *   @history 2016-08-15 Code adapted from lrowacpho written by Kris Becker

  *

  **/

  void PhotometricFunction::setMaximumIncidenceAngle(double angle) {

    m_maximumIncidenceAngle = angle;

  }

 /**

  * Mutator function to set minimum emission angle

  *

  * @param double The new minimum emission angle

  *

  * @author 2016-08-15 Victor Silva

  *

  * @internal

  *   @history 2016-08-15 Code adapted from lrowacpho written by Kris Becker

  *

  **/

  void PhotometricFunction::setMinimumEmissionAngle(double angle) {

    m_minimumEmissionAngle = angle;

  }

 /**

  * Mutator function to set maximum emission angle

  *

  * @param double The new maximum emission angle

  *

  * @author 2016-08-15 Victor Silva

  *

  * @internal

  *   @history 2016-08-15 Code adapted from lrowacpho written by Kris Becker

  *

  **/

  void PhotometricFunction::setMaximumEmissionAngle(double angle) {

    m_maximumEmissionAngle = angle;

  }

 /**

  * Mutator function to set minimum phase angle

  *

  * @param double The new minimum phase angle

  *

  * @author 2016-08-15 Victor Silva

  *

  * @internal

  *   @history 2016-08-15 Code adapted from lrowacpho written by Kris Becker

  *

  **/

  void PhotometricFunction::setMinimumPhaseAngle(double angle) {

    m_minimumPhaseAngle = angle;

  }

 /**

  * Mutator function to set maximum phase angle

  *

  * @param double The new maximum phase angle

  *

  * @author 2016-08-15 Victor Silva

  *