OverlapNormalization.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#include "OverlapNormalization.h"
 
#include <iomanip>
 
#include "BasisFunction.h"
#include "IException.h"
#include "LeastSquares.h"
#include "Statistics.h"
 
using namespace std;
 
namespace Isis {
 
  OverlapNormalization::OverlapNormalization(std::vector<Statistics *> statsList) {
    m_gainFunction = NULL;
    m_gainLsq = NULL;
    m_offsetFunction = NULL;
    m_offsetLsq = NULL;
 
    m_statsList = statsList;
    m_gainFunction = new BasisFunction("BasisFunction",
                                       statsList.size(), statsList.size());
    m_gainLsq = new LeastSquares(*m_gainFunction);
    m_offsetFunction = new BasisFunction("BasisFunction",
                                         statsList.size(), statsList.size());
    m_offsetLsq = new LeastSquares(*m_offsetFunction);
 
    m_gains.resize(statsList.size());
    m_offsets.resize(statsList.size());
    for (unsigned int i = 0; i < statsList.size(); i++) {
      m_gains[i] = 1.0;
      m_offsets[i] = 0.0;
    }
    m_solved = false;
  }
 
  OverlapNormalization::~OverlapNormalization() {
    if (m_gainFunction != NULL) delete m_gainFunction;
    if (m_offsetFunction != NULL) delete m_offsetFunction;
    if (m_gainLsq != NULL) delete m_gainLsq;
    if (m_offsetLsq != NULL) delete m_offsetLsq;
    for (unsigned int i = 0; i < m_statsList.size(); i++) delete m_statsList[i];
  };
 
  OverlapNormalization::AddStatus OverlapNormalization::AddOverlap(
    const Statistics &area1, const unsigned index1,
    const Statistics &area2, const unsigned index2,
    double weight) {
    if (index1 >= m_statsList.size()) {
      string msg = "The index 1 is outside the bounds of the list.";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    if (index2 >=  m_statsList.size()) {
      string msg = "The index 2 is outside the bounds of the list.";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    // If there is no overlapping area, then the overlap is invalid
    if (area1.ValidPixels() == 0 || area2.ValidPixels() == 0) {
      return NoOverlap;
    }
 
    // The weight must be a positive real number
    if (weight <= 0.0) {
      string msg = "All weights must be positive real numbers.";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    OverlapNormalization::Overlap o;
    o.area1 = area1;
    o.area2 = area2;
    o.index1 = index1;
    o.index2 = index2;
 
    double avg1 = area1.Average();
    double avg2 = area2.Average();
 
    // Averages must not be 0 to avoid messing up least squares
    if (avg1 == 0 || avg2 == 0) return NoContrast;
 
    m_overlapList.push_back(o);
    m_deltas.push_back(avg2 - avg1);
    m_weights.push_back(weight);
    m_solved = false;
    return Success;
  }
 
  void OverlapNormalization::Solve(SolutionType type,
                                   LeastSquares::SolveMethod method) {
    // Make sure that there is at least one overlap
    if (m_overlapList.size() == 0) {
      string msg = "None of the input images overlap";
      throw IException(IException::User, msg, _FILEINFO_);
    }
 
    // Make sure the number of valid overlaps + hold images is greater than the
    // number of input images (otherwise the least squares equation will be
    // unsolvable due to having more unknowns than knowns)
    if (m_overlapList.size() + m_idHoldList.size() < m_statsList.size()) {
      string msg = "Unable to normalize overlaps. The number of overlaps and "
                   "holds must be greater than the number of input images";
      throw IException(IException::User, msg, _FILEINFO_);
    }
    
    if ( method == LeastSquares::SPARSE ) {
      m_offsetLsq = NULL;
      m_gainLsq = NULL;
      delete m_offsetLsq;
      delete m_gainLsq;
      int sparseMatrixRows = m_overlapList.size() + m_idHoldList.size();
      int sparseMatrixCols = m_offsetFunction->Coefficients();
      m_offsetLsq = new LeastSquares(*m_offsetFunction, true, sparseMatrixRows, sparseMatrixCols, true);
      sparseMatrixCols = m_gainFunction->Coefficients();
      m_gainLsq = new LeastSquares(*m_gainFunction, true, sparseMatrixRows, sparseMatrixCols, true);
      const std::vector<double> alphaWeight(sparseMatrixCols, 1/1000.0);
      m_offsetLsq->SetParameterWeights( alphaWeight );
      m_gainLsq->SetParameterWeights( alphaWeight );
    }
 
    // Calculate offsets
    if (type != Gains && type != GainsWithoutNormalization) {
      // Add knowns to least squares for each overlap
      for (int overlap = 0; overlap < (int)m_overlapList.size(); overlap++) {
        Overlap curOverlap = m_overlapList[overlap];
        int id1 = curOverlap.index1;
        int id2 = curOverlap.index2;
 
        vector<double> input;
        input.resize(m_statsList.size());
        for (int i = 0; i < (int)input.size(); i++) input[i] = 0.0;
        input[id1] = 1.0;
        input[id2] = -1.0;
 
        m_offsetLsq->AddKnown(input, m_deltas[overlap], m_weights[overlap]);
      }
 
      // Add a known to the least squares for each hold image
      for (int h = 0; h < (int)m_idHoldList.size(); h++) {
        int hold = m_idHoldList[h];
 
        vector<double> input;
        input.resize(m_statsList.size());
        for (int i = 0; i < (int)input.size(); i++) input[i] = 0.0;
        input[hold] = 1.0;
        m_offsetLsq->AddKnown(input, 0.0, 1e30);
      }
 
      // Solve the least squares and get the offset coefficients to apply to the
      // images
      m_offsets.resize(m_statsList.size());
      m_offsetLsq->Solve(method);
      for (int i = 0; i < m_offsetFunction->Coefficients(); i++) {
        m_offsets[i] = m_offsetFunction->Coefficient(i);
      }
    }
 
    // Calculate Gains
    if (type != Offsets) {
      // Add knowns to least squares for each overlap
      for (int overlap = 0; overlap < (int)m_overlapList.size(); overlap++) {
        Overlap curOverlap = m_overlapList[overlap];
        int id1 = curOverlap.index1;
        int id2 = curOverlap.index2;
 
        vector<double> input;
        input.resize(m_statsList.size());
        for (int i = 0; i < (int)input.size(); i++) input[i] = 0.0;
        input[id1] = 1.0;
        input[id2] = -1.0;
 
        double tanp;
 
        if (type != GainsWithoutNormalization) {
          if (curOverlap.area1.StandardDeviation() == 0.0) {
            tanp = 0.0;    // Set gain to 1.0
          }
          else {
            tanp = curOverlap.area2.StandardDeviation()
                   / curOverlap.area1.StandardDeviation();
          }
        }
        else {
          if (curOverlap.area1.Average() == 0.0) {
            tanp = 0.0;
          }
          else {
            tanp = curOverlap.area2.Average() / curOverlap.area1.Average();
          }
        }
 
        if (tanp > 0.0) {
          m_gainLsq->AddKnown(input, log(tanp), m_weights[overlap]);
        }
        else {
          m_gainLsq->AddKnown(input, 0.0, 1e10); // Set gain to 1.0
        }
      }
 
      // Add a known to the least squares for each hold image
      for (int h = 0; h < (int)m_idHoldList.size(); h++) {
        int hold = m_idHoldList[h];
 
        vector<double> input;
        input.resize(m_statsList.size());
        for (int i = 0; i < (int)input.size(); i++) input[i] = 0.0;
        input[hold] = 1.0;
        m_gainLsq->AddKnown(input, 0.0, 1e10);
      }
 
      // Solve the least squares and get the gain coefficients to apply to the
      // images
      m_gains.resize(m_statsList.size());
      m_gainLsq->Solve(method);
      for (int i = 0; i < m_gainFunction->Coefficients(); i++) {
        m_gains[i] = exp(m_gainFunction->Coefficient(i));
      }
    }
 
    m_solved = true;
  }
 
  double OverlapNormalization::Average(const unsigned index) const {
    if (index >= m_statsList.size()) {
      string msg = "The index was out of bounds for the list of statistics.";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    return m_statsList[index]->Average();
  }
 
  double OverlapNormalization::Gain(const unsigned index) const {
    if (index >= m_statsList.size()) {
      string msg = "The index was out of bounds for the list of statistics.";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    return m_gains[index];
  }
 
  double OverlapNormalization::Offset(const unsigned index) const {
    if (index >= m_statsList.size()) {
      string msg = "The index was out of bounds for the list of statistics.";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    return m_offsets[index];
  }
 
  double OverlapNormalization::Evaluate(double dn, unsigned index) const {
    if (!m_solved) {
      string msg = "The least squares equation has not been successfully ";
      msg += "solved yet.";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    if (Isis::IsSpecial(dn)) return dn;
    return (dn - Average(index)) * Gain(index) + Average(index) + Offset(index);
  }
}