MixtureModelLinearStrategy.cpp

Go to the documentation of this file.

/*  Copyright (C) 2008 National Institute For Space Research (INPE) - Brazil.

    This file is part of the TerraLib - a Framework for building GIS enabled applications.

    TerraLib is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License,
    or (at your option) any later version.

    TerraLib is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with TerraLib. See COPYING. If not, write to
    TerraLib Team at <terralib-team@terralib.org>.
 */

/*!
  \file terralib/rp/MixtureModelLinearStrategy.cpp

  \brief Raster linear strategy for mixture model classification.
*/

// TerraLib
#include "../common/MatrixUtils.h"
#include "../common/progress/TaskProgress.h"
#include "../raster/Grid.h"
#include "../raster/Utils.h"
#include "Functions.h"
#include "Macros.h"
#include "MixtureModelLinearStrategy.h"
#include "Functions.h"

// Boost
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix.hpp>

namespace
{
  static te::rp::MixtureModelLinearStrategyFactory mixtureModelLinearStrategyFactoryInstance;
}

te::rp::MixtureModelLinearStrategy::Parameters::Parameters()
{
  reset();
}

te::rp::MixtureModelLinearStrategy::Parameters::~Parameters()
{
}

const te::rp::MixtureModelLinearStrategy::Parameters& te::rp::MixtureModelLinearStrategy::Parameters::operator=(const te::rp::MixtureModelLinearStrategy::Parameters& rhs)
{
  reset();

  return *this;
}

void te::rp::MixtureModelLinearStrategy::Parameters::reset() throw(te::rp::Exception)
{
}

te::common::AbstractParameters* te::rp::MixtureModelLinearStrategy::Parameters::clone() const
{
  return new te::rp::MixtureModelLinearStrategy::Parameters(*this);
}

te::rp::MixtureModelLinearStrategy::MixtureModelLinearStrategy()
{
  m_isInitialized = false;
}

te::rp::MixtureModelLinearStrategy::~MixtureModelLinearStrategy()
{
}

bool te::rp::MixtureModelLinearStrategy::initialize(te::rp::StrategyParameters const* const strategyParams) throw(te::rp::Exception)
{
  m_isInitialized = false;

  te::rp::MixtureModelLinearStrategy::Parameters const* paramsPtr = dynamic_cast<te::rp::MixtureModelLinearStrategy::Parameters const*>(strategyParams);

  if(!paramsPtr)
    return false;

  m_parameters = *(paramsPtr);

  m_isInitialized = true;

  return true;
}

bool te::rp::MixtureModelLinearStrategy::execute(const te::rst::Raster& inputRaster, const std::vector<unsigned int>& inputRasterBands,
                                                 const std::vector<std::string>& inputSensorBands, const std::map<std::string, std::vector<double> >& components,
                                                 std::vector<te::rst::Raster*>& outputRaster, const bool normalize, const bool enableProgressInterface) throw(te::rp::Exception)
{
  TERP_TRUE_OR_RETURN_FALSE(m_isInitialized, "Instance not initialized")

  const unsigned int nComponents = (unsigned int)components.size();
  const unsigned int nBands = (unsigned int)inputRasterBands.size();

  double dummy = outputRaster[0]->getBand(0)->getProperty()->m_noDataValue;

  std::vector<double> Qmax;
  std::vector<double> Qmin;
  for (unsigned i = 0; i < inputSensorBands.size(); i++)
  {
    Qmax.push_back(GetDigitalNumberBandMax(inputSensorBands[i]));
    Qmin.push_back(GetDigitalNumberBandMin(inputSensorBands[i]));
  }

  te::rst::Grid* tmpgrid = new te::rst::Grid(*inputRaster.getGrid());

  std::vector<te::rst::BandProperty*> tmpbandsProperties;
  for (unsigned i = 0; i < outputRaster[0]->getNumberOfBands(); i++)
  {
    te::rst::BandProperty* bprop = new te::rst::BandProperty(*(outputRaster[0]->getBand(i)->getProperty()));
    bprop->m_colorInterp = te::rst::GrayIdxCInt;
    bprop->m_type = te::dt::DOUBLE_TYPE;
    tmpbandsProperties.push_back(bprop);
  }
  std::unique_ptr<te::rst::Raster> tmpRaster (te::rst::RasterFactory::make("MEM", tmpgrid, tmpbandsProperties, std::map<std::string, std::string>()));

  // used to calculate minimum e maximum values of outputraster using statitical minmax algorithm
  std::vector<double> tmp_min;
  std::vector<double> tmp_max;
  for (unsigned i = 0; i < tmpRaster->getNumberOfBands(); i++)
  {
    tmp_min.push_back(std::numeric_limits< double >::max());
    tmp_max.push_back(-std::numeric_limits< double >::max());
    if (normalize)
    {
      m_min.push_back(0);
      m_max.push_back(255);
    }
    else
    {
      m_min.push_back(std::numeric_limits< double >::max());
      m_max.push_back(-std::numeric_limits< double >::max());
    }
    m_minerror.push_back(std::numeric_limits< double >::max());
    m_maxerror.push_back(-std::numeric_limits< double >::max());
  }

  if (m_transfMatrix.size1() == 0)
    generateTransformMatrix(inputRasterBands, inputSensorBands, components);
  else
    m_transposeA = boost::numeric::ublas::trans(m_transfMatrix);

  // calculate the product A' * A
  boost::numeric::ublas::matrix<double> productAtA = boost::numeric::ublas::matrix<double>(m_transposeA.size1(), m_transfMatrix.size2());
  productAtA = boost::numeric::ublas::prod(m_transposeA, m_transfMatrix);

  // calculate the inverse of A' * A
  boost::numeric::ublas::matrix<double> invertAtA = boost::numeric::ublas::matrix<double>(productAtA.size1(), productAtA.size2());
  te::common::GetInverseMatrix(productAtA, invertAtA);

  boost::numeric::ublas::matrix<double> matrixR = boost::numeric::ublas::matrix<double>(m_transfMatrix.size1(), 1);
  boost::numeric::ublas::matrix<double> productAtR = boost::numeric::ublas::matrix<double>(m_transposeA.size1(), matrixR.size2());
  boost::numeric::ublas::matrix<double> matrixX = boost::numeric::ublas::matrix<double>(invertAtA.size1(), productAtR.size2());
  boost::numeric::ublas::matrix<double> productAX = boost::numeric::ublas::matrix<double>(m_transfMatrix.size1(), matrixX.size2());
  boost::numeric::ublas::matrix<double> matrixE = boost::numeric::ublas::matrix<double>(matrixR.size1(), matrixR.size2());

  te::common::TaskProgress task(TE_TR("Linear Mixture Model"), te::common::TaskProgress::UNDEFINED, inputRaster.getNumberOfRows());
  double value;

  for (unsigned r = 0; r < inputRaster.getNumberOfRows(); r++)
  {
    for (unsigned c = 0; c < inputRaster.getNumberOfColumns(); c++)
    {
// read input values
      for (unsigned b = 0; b < nBands; b++)
      {
        inputRaster.getValue(c, r, value, inputRasterBands[b]);
        if (value == dummy)
          matrixR(b, 0) = 0.;
        else
          matrixR(b, 0) = value / Qmax[b];
      }
      matrixR(nBands, 0) = 0.0;
// calculate the product A' * R
      productAtR = boost::numeric::ublas::prod(m_transposeA, matrixR);
// compute matrix X
      matrixX = boost::numeric::ublas::prod(invertAtA, productAtR);
// write output fractions
      for (unsigned b = 0; b < nBands; b++)
      {
        value = matrixX(b + 1, 0);
        tmpRaster->setValue(c, r, value, b);
        if (value != dummy)
        {
          tmp_min[b] = value < tmp_min[b] ? value : tmp_min[b];
          tmp_max[b] = value > tmp_max[b] ? value : tmp_max[b];
        }
      }

// compute error matrix
      if (outputRaster.size() > 1)
      {
        productAX = boost::numeric::ublas::prod(m_transfMatrix, matrixX);
        matrixE = matrixR - productAX;

// write output errors
        for (unsigned b = 0, e = 0; b < nBands; b++, e++)
        {
          value = matrixE(e, 0);
          outputRaster[1]->setValue(c, r, value, b);
          m_minerror[b] = value < m_minerror[b] ? value : m_minerror[b];
          m_maxerror[b] = value > m_maxerror[b] ? value : m_maxerror[b];
        }
      }
    }
    task.pulse();
  }

  //Transform to output type
  for (unsigned r = 0; r < tmpRaster->getNumberOfRows(); r++)
  {
    for (unsigned c = 0; c < tmpRaster->getNumberOfColumns(); c++)
    {
      for (unsigned b = 0; b < nComponents; b++)
      {
        tmpRaster->getValue(c, r, value, b);
        if (value != dummy)
        {
          if (normalize) //to 8 bits
            value = 255 * (value - tmp_min[b]) / (tmp_max[b] - tmp_min[b]);
          else
            value = (Qmax[b] - Qmin[b]) * (value - tmp_min[b]) / (tmp_max[b] - tmp_min[b]);
        }
        outputRaster[0]->setValue(c, r, value, b);
        m_min[b] = value < m_min[b] ? value : m_min[b];
        m_max[b] = value > m_max[b] ? value : m_max[b];
      }
    }
  }

  return true;
}

bool te::rp::MixtureModelLinearStrategy::getTransformMatrix(
  boost::numeric::ublas::matrix<double>& matrix) const
{
  matrix = m_transfMatrix;
  return true;
}


bool te::rp::MixtureModelLinearStrategy::setTransformMatrix(
  boost::numeric::ublas::matrix<double>& matrix) {
  m_transfMatrix = matrix;
  return true;
}

bool te::rp::MixtureModelLinearStrategy::generateTransformMatrix(const std::vector<unsigned int>& inputRasterBands, const std::vector<std::string>& inputSensorBands,
  const std::map<std::string, std::vector<double> >& components) 
{
  const unsigned int nComponents = (unsigned int)components.size();
  const unsigned int nBands = (unsigned int)inputRasterBands.size();

  boost::numeric::ublas::matrix<double> transposeA = boost::numeric::ublas::matrix<double>(nComponents + 1, nBands + 1);
  boost::numeric::ublas::matrix<double> matrixA(transposeA.size2(), transposeA.size1());

  // by definition, A * X = R, we have A, so that X = inv(A' * A) * A' * R
  // A is the set of known reflectances for each component
  // R is the reflectances of a specific pixel
  // X is the proportion of each component

  std::vector<double> Qmax;
  for (unsigned i = 0; i < inputSensorBands.size(); i++)
  {
    Qmax.push_back(GetDigitalNumberBandMax(inputSensorBands[i]));
  }

  // retrieve the transposed A matrix, normalizing values using band/sensor info
  unsigned row = 0;
  for (unsigned j = 0; j < nBands; j++)
    transposeA(row, j) = 1.0;
  transposeA(row, nBands) = 0.0;

  std::map<std::string, std::vector<double> >::const_iterator it;
  for (row = 1, it = components.begin(); it != components.end(); row++, it++)
  {
    for (unsigned j = 0; j < nBands; j++)
      transposeA(row, j) = it->second[j] / Qmax[j];
    transposeA(row, nBands) = 1.0;
  }

  // calculate the matrixA from the transpose of A
  matrixA = boost::numeric::ublas::trans(transposeA);
  m_transfMatrix = matrixA;
  m_transposeA = transposeA;

  return true;
}

bool te::rp::MixtureModelLinearStrategy::getMinMax(std::vector<double>& min, std::vector<double>& max) const
{
  min = m_min;
  max = m_max;
  return true;
}

bool te::rp::MixtureModelLinearStrategy::getMinMaxError(std::vector<double>& min, std::vector<double>& max) const
{
  min = m_minerror;
  max = m_maxerror;
  return true;
}


te::rp::MixtureModelLinearStrategyFactory::MixtureModelLinearStrategyFactory()
  : te::rp::MixtureModelStrategyFactory("linear")
{
}

te::rp::MixtureModelLinearStrategyFactory::~MixtureModelLinearStrategyFactory()
{
}

te::rp::MixtureModelStrategy* te::rp::MixtureModelLinearStrategyFactory::build()
{
  return new te::rp::MixtureModelLinearStrategy();
}