Gruen.h

#ifndef Gruen_h
#define Gruen_h
 
/* SPDX-License-Identifier: CC0-1.0 */
 
 
#include "AutoReg.h"
#include "tnt/tnt_array2d.h"
#include "tnt/tnt_array1d.h"
#include "Affine.h"
#include "DbProfile.h"
#include "GruenTypes.h"
#include "CollectorMap.h"
#include "Statistics.h"
#include "Pvl.h"
#include "PvlGroup.h"
 
namespace Isis {
  class Chip;
 
  class Gruen : public AutoReg {
    public:
 
      Gruen();
      Gruen(Pvl &pvl);
 
      virtual ~Gruen() { }
 
      BigInt CallCount() const { return (m_callCount); }
 
      void WriteSubsearchChips(const QString &pattern = "SubChip");
 
      AffineTolerance getAffineTolerance() const;
 
      inline double getSpiceConstraint() const { return (m_spiceTol); }
 
      inline double getAffineConstraint() const { return (m_affineTol); }
 
      void setAffineRadio(const AffineRadio &affrad);
      void setAffineRadio();
 
      const AffineRadio &getDefaultAffineRadio() const { return (m_defAffine);}
 
      const AffineRadio &getAffineRadio() const { return (m_affine);  }
 
      MatchPoint getLastMatch() const { return (m_point);   }
 
    protected:
      virtual QString AlgorithmName() const {
        return ("Gruen");
      }
 
      int algorithm(Chip &pattern, Chip &subsearch, const Radiometric &radio,
                    BigInt &ptsUsed, double &resid, GMatrix &atai,
                    AffineRadio &affrad);
 
      Analysis errorAnalysis(const BigInt &npts, const double &resid,
                             const GMatrix &atai);
 
      //  These methods are for AutoReg non-adaptive requirements
      virtual double MatchAlgorithm(Chip &pattern, Chip &subsearch);
      virtual bool CompareFits(double fit1, double fit2);
 
      virtual double IdealFit() const {  return (0.0); }
 
      // AutoReg adaptive method
      virtual AutoReg::RegisterStatus Registration(Chip &sChip,
          Chip &pChip,
          Chip &fChip,
          int startSamp,
          int startLine,
          int endSamp,
          int endLine,
          int bestSamp,
          int bestLine);
 
      virtual Pvl AlgorithmStatistics(Pvl &pvl);
 
    private:
      enum ErrorTypes {  NotEnoughPoints = 1, CholeskyFailed = 2,
                         EigenSolutionFailed = 3, AffineNotInvertable = 4,
                         MaxIterationsExceeded = 5, RadShiftExceeded = 6,
                         RadGainExceeded = 7, MaxEigenExceeded = 8,
                         AffineDistExceeded = 9
                       };
 
 
      struct ErrorCounter {
        ErrorCounter() : m_gerrno(0), m_keyname("Unknown"), m_count(0) { }
        ErrorCounter(int gerrno, const QString &keyname) : m_gerrno(gerrno),
          m_keyname(keyname), m_count(0) { }
 
        inline int  Errno() const {  return (m_gerrno); }
        inline BigInt Count() const { return (m_count); }
        inline void BumpIt() { m_count++; }
        PvlKeyword LogIt() const { return (PvlKeyword(m_keyname, toString(m_count))); }
 
        int         m_gerrno;
        QString m_keyname;
        BigInt      m_count;
      };
 
      typedef CollectorMap<int, ErrorCounter> ErrorList;
 
      // Iteration loop variables
      BigInt       m_callCount;
      QString  m_filePattern;
 
      ErrorList    m_errors;            // Error logger
      BigInt       m_unclassified;      // Unclassified errors
 
      // Tolerance and count parameters
      int          m_maxIters;          // Maximum iterations
      int          m_nIters;            // Number iterations
      BigInt       m_totalIterations;   // Running total iteration count
 
      DbProfile    m_prof;              //  Parameter profile
      double       m_transTol;          //  Affine translation tolerance
      double       m_scaleTol;          //  Affine scale tolerance
      double       m_shearTol;          //  Affine shear tolerance
      double       m_spiceTol;          // SPICE tolerance
      double       m_affineTol;         // Affine Tolerance
 
      double       m_shiftTol;          // Radiometric shift tolerance
      double       m_rgainMinTol;       // Radiometric Gain minimum
      double       m_rgainMaxTol;       // Radiometric Gain maximum
 
      double       m_defGain;           // Default Radiometric gain
      double       m_defShift;          // Default Radiometric shift
 
      // These are for recomputing SMTK points
      AffineRadio  m_defAffine;         // Default affine/radiometric settings
      AffineRadio  m_affine;            // Incoming affine setting
      MatchPoint   m_point;             // Last resuting registration result
 
      //  Statistics gathered during processing
      Statistics   m_eigenStat;
      Statistics   m_iterStat;
      Statistics   m_shiftStat;
      Statistics   m_gainStat;
 
 
      //  Static init of default PVL parameters
      static Pvl &getDefaultParameters();
 
      QString ConfKey(const DbProfile &conf, const QString &keyname,
                const QString &defval, int index = 0) const {
        if(!conf.exists(keyname)) {
          return (defval);
        }
        if(conf.count(keyname) < index) {
          return (defval);
        }
        QString iValue(conf.value(keyname, index));
        IString tmp(iValue);
        QString value = tmp.ToQt();  // This makes it work with a string?
        return (value);
      };
 
      template <typename T>
      T ConfKey(const DbProfile &conf, const QString &keyname,
                const T &defval, int index = 0) const {
        if(!conf.exists(keyname)) {
          return (defval);
        }
        if(conf.count(keyname) < index) {
          return (defval);
        }
        QString iValue(conf.value(keyname, index));
        IString tmp(iValue);
        T value = tmp;  // This makes it work with a string?
        return (value);
      };
 
 
      template <typename T>
      PvlKeyword ParameterKey(const QString &keyname,
                              const T &value,
                              const QString &unit = "") const {
        if(m_prof.exists(keyname)) {
          return(ValidateKey(keyname, value, unit));
        }
        return (PvlKeyword(keyname, "Unbounded"));
      };
 
 
      inline PvlKeyword ValidateKey(const QString keyname,
                                    const double &value,
                                    const QString &unit = "") const {
        if(IsSpecial(value)) {
          return (PvlKeyword(keyname, "NULL"));
        }
        else {
          return (PvlKeyword(keyname, toString(value), unit));
        }
      }
 
      ErrorList initErrorList() const;
      int logError(int gerrno, const QString &gerrmsg);
      PvlGroup StatsLog() const;
      PvlGroup ParameterLog() const;
 
      void init(Pvl &pvl);
 
      inline Radiometric getDefaultRadio() const {
        return (Radiometric(m_defShift, m_defGain));
      }
 
      inline double DegreesOfFreedom(const int npts) const {
        return ((double) (npts - NCONSTR));
      }
 
      void resetStats();
 
      GMatrix Identity(const int &ndiag = 3) const;
      GMatrix Choldc(const GMatrix &a, GVector &p) const;
      GMatrix Cholsl(const GMatrix &a, const GVector &p,
                     const GMatrix &b, const GMatrix &x) const;
      int Jacobi(const GMatrix &a, GVector &evals, GMatrix &evecs,
                  const int &MaxIters = 50) const;
      void EigenSort(GVector &evals, GMatrix &evecs) const;
      BigInt MinValidPoints(BigInt totalPoints) const;
      bool ValidPoints(BigInt totalPoints, BigInt nPoints) const;
 
      int CheckConstraints(MatchPoint &point);
      Coordinate getChipUpdate(Chip &sChip, MatchPoint &point) const;
      AutoReg::RegisterStatus Status(const MatchPoint &result);
  };
 
}  // namespace Isis
#endif