#ifndef Gruen_h
#define Gruen_h
/**
* @file
* $Revision: 1.4 $
* $Date: 2010/06/14 20:45:52 $
*
* 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 & 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 "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;
/**
* @brief Gruen pattern matching algorithm
*
* This class provides adaptive image (chip) registration using
* the AutoReg factory architecture. This algorithm uses an
* Affine transform to iteratively adjust the search chip at
* each iteration. Each iteration solves for new adjustments to
* the Affine transform until the 6 affine parameters fall below the
* tolerances as specified in AffineThreshHold1 and AffineThreshHold2.
*
* This class minimizes the 6 specifiable Affine transform components of a 3x3
* matrix. The three Affine components for X (sample) and Y (line) are scale,
* shear and translation. Gruen provides control over the maximum values
* these three components should attain in order for the registration to
* converge to a successful match. These limits are specified
* by AffineScaleTolerance, AffineShearTolerance and
* AffineTranslationTolerance. AffineShearTolerance is
* optional and if not specified, it defaults to the value of
* AffineScaleTolerance. These tolerances specify the maximum amount
* of translation pixels can be shifted between one Gruen iteration and
* another. For example, AffineTranslationTolerance = 0.2 means that a
* subsearch chip cannot move in sample or line direction more than 0.2 pixels
* in order to satisfy convergence. AffineScaleTolerance constrains
* the sample and line scale elements of the Affine transformation. And
* AffineShearTolerance constrains the sample and line shear elements
* of the Affine transformation. The scale and shear parameters are scaled by
* size of the Chip. The sample scale and shear Affine components maximum
* limit is computed as AffineScaleTolerance/((#Samples-1)/2) and
* AffineShearTolerance/((#Samples-1)/2). Likewise, the line scale and
* shear maximums are computed using #Lines in the preceding equation.
*
* @ingroup PatternMatching
*
* @see MinimumGruen AdaptiveGruen AutoReg
*
* @author 2009-09-09 Kris Becker
*
* @internal
* @history 2009-09-10 Kris Becker Changed AffineThreshHold1 to
* AffineTranslationTolerance and AffineThreshHold2 to
* AffineScaleTolerance. Added AffineShearTolerance. Added logging
* of Gruen parameters.
* @history 2009-09-11 Kris Becker Minor error in computation of radiometric
* shift
* @history 2010-06-14 Jeannie Walldren - Capitalized getTransform and
* setTransform to match changes in Chip class
* @history 2011-05-23 Kris Becker - Reworked major portions of
* implementation for a more modular support.
*/
class Gruen : public AutoReg {
public:
/** Default constructor */
Gruen();
Gruen(Pvl &pvl);
//! Destructor
virtual ~Gruen() { }
/** Returns the current call count */
BigInt CallCount() const { return (m_callCount); }
void WriteSubsearchChips(const QString &pattern = "SubChip");
AffineTolerance getAffineTolerance() const;
/** Returns the SPICE tolerance constraint as read from config file */
inline double getSpiceConstraint() const { return (m_spiceTol); }
/** Returns the Affine tolerance constraint as read from config file */
inline double getAffineConstraint() const { return (m_affineTol); }
void setAffineRadio(const AffineRadio &affrad);
void setAffineRadio();
/** Returns default settings for Affine/Radiometric parameters */
const AffineRadio &getDefaultAffineRadio() const { return (m_defAffine);}
/** Return current state of Affine/Radio state */
const AffineRadio &getAffineRadio() const { return (m_affine); }
/**
* @brief Returns the register state of the last successful Gruen match
*
* This method returns the full match condition of the last call to Gruen
* Register function that was successful.
*
* BEWARE: This is only valid if Register returns successfully! This is
* due to AutoReg returning conditions that occur prior to the actual
* Gruen algorithm being called.
*
* @author Kris Becker - 6/4/2011
*
* @return MatchPoint
*/
MatchPoint getLastMatch() const { return (m_point); }
protected:
/** Returns the default name of the algorithm as Gruen */
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);
/** Returns the ideal fit for a perfect Gruen result */
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:
/** Error enumeration values */
enum ErrorTypes { NotEnoughPoints = 1, CholeskyFailed = 2,
EigenSolutionFailed = 3, AffineNotInvertable = 4,
MaxIterationsExceeded = 5, RadShiftExceeded = 6,
RadGainExceeded = 7, MaxEigenExceeded = 8,
AffineDistExceeded = 9
};
/** Struct that maintains error counts */
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;
};
/// Declaration of error count list
typedef CollectorMap 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();
/**
* @brief Helper method to initialize parameters
*
* This method will check the existance of a keyword and extract the value
* if it exists to the passed parameter (type). If it doesn't exist, the
* default values is returned.
*
* @param T Templated variable type
* @param conf Parameter profile container
* @param keyname Name of keyword to get a value from
* @param defval Default value it keyword/value doesn't exist
* @param index Optional index of the value for keyword arrays
*
* @return T Return type
*/
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);
};
/**
* @brief Helper method to initialize parameters
*
* This method will check the existance of a keyword and extract the value
* if it exists to the passed parameter (type). If it doesn't exist, the
* default values is returned.
*
* @param T Templated variable type
* @param conf Parameter profile container
* @param keyname Name of keyword to get a value from
* @param defval Default value it keyword/value doesn't exist
* @param index Optional index of the value for keyword arrays
*
* @return T Return type
*/
template
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);
};
/**
* @brief Keyword formatter for Gruen parameters
*
* Constructs a keyword with actual user/programmer values if provided,
* otherwise sets the value to "Unbounded".
*
* @param T Type of value to record
* @param keyname Name of keyword to create
* @param value Value to set keyword to if in profile
* @param unit Optional unit value
*
* @return PvlKeyword Constructed keyword for key/value
*/
template
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"));
};
/**
* @brief Checks value of key, produces appropriate value
*
* This function checks the value of the keyword for specialness
* and will create the appropriate keyword if it is special.
*
* @param keyname Name of keyword to create
* @param value Keyword value
* @param unit Optional unit qualifer with value
*
* @return PvlKeyword Returns newly created keyword/value
*/
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);
/** Returns the default radiometric gain value */
inline Radiometric getDefaultRadio() const {
return (Radiometric(m_defShift, m_defGain));
}
/** Returns number of degrees of freedom of points */
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