ClassifierSAMStrategy.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/ClassifierSAMStrategy.cpp

  \brief Maximum a posteriori probability strategy.
*/

// TerraLib

#include "ClassifierSAMStrategy.h"
#include "Macros.h"
#include "../common/progress/TaskProgress.h"
#include "../raster/Band.h"

#include <cfloat>
#include <cmath>

namespace
{
  static te::rp::ClassifierSAMStrategyFactory classifierMAPStrategyFactoryInstance;
}

//-----------------------------------------------------------------------------

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

te::rp::ClassifierSAMStrategy::Parameters::~Parameters() = default;

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

  m_trainSamplesPtr = rhs.m_trainSamplesPtr;
  m_maxAngularDistances = rhs.m_maxAngularDistances;

  return *this;
}

void te::rp::ClassifierSAMStrategy::Parameters::reset() throw(te::rp::Exception)
{
  m_trainSamplesPtr.reset();
  m_maxAngularDistances.clear();
}

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

//-----------------------------------------------------------------------------

te::rp::ClassifierSAMStrategy::ClassifierSAMStrategy()
{
  m_isInitialized = false;
  m_highestTrainClassID = 0;
}

te::rp::ClassifierSAMStrategy::~ClassifierSAMStrategy() = default;

bool te::rp::ClassifierSAMStrategy::initialize(
  te::rp::ClassifierStrategyParameters const* const strategyParams) throw(te::rp::Exception)
{
  m_isInitialized = false;
  m_initParams.reset();
  m_classesMeans.clear();
  m_classesIndex2ID.clear();
  
  ClassifierSAMStrategy::Parameters const * const castParamsPtr = 
    dynamic_cast< ClassifierSAMStrategy::Parameters const * >( strategyParams );
  TERP_INSTANCE_TRUE_OR_RETURN_FALSE( castParamsPtr, "Invalid parameters" );
  
  // Checking the angluar distances
  
  TERP_INSTANCE_TRUE_OR_RETURN_FALSE( castParamsPtr->m_maxAngularDistances.size() ==
    castParamsPtr->m_trainSamplesPtr->size(),
    "Mismatch between the classes number and the max angular distances" );
    
  for( unsigned int maxAngularDistancesIdx = 0 ; maxAngularDistancesIdx <
    castParamsPtr->m_maxAngularDistances.size() ; ++ maxAngularDistancesIdx )
  {
    TERP_INSTANCE_TRUE_OR_RETURN_FALSE( ( castParamsPtr->m_maxAngularDistances[
      maxAngularDistancesIdx ] >= 0.0 ),
      "Invalid max angular distance" );
  }  
  
  // Calculating m_classesMeans and m_classesIndex2ID
  
  TERP_INSTANCE_TRUE_OR_RETURN_FALSE( castParamsPtr->m_trainSamplesPtr,
    "Invalid classes samples pointer" )
  TERP_INSTANCE_TRUE_OR_RETURN_FALSE( castParamsPtr->m_trainSamplesPtr->size() > 0,
    "Invalid classes samples number" )
  TERP_INSTANCE_TRUE_OR_RETURN_FALSE( castParamsPtr->m_trainSamplesPtr->begin()->second.size() > 0,
    "Invalid classe samples number" )    
  const unsigned int dimsNumber = (unsigned int)
    castParamsPtr->m_trainSamplesPtr->begin()->second.operator[]( 0 ).size();
  TERP_INSTANCE_TRUE_OR_RETURN_FALSE( dimsNumber > 0,
    "Invalid dimensions number" )       
  
  {
    ClassesSamplesT::const_iterator classesIt = 
      castParamsPtr->m_trainSamplesPtr->begin();
    const ClassesSamplesT::const_iterator classesItE = 
      castParamsPtr->m_trainSamplesPtr->end();
    unsigned int dimIdx = 0;
    SamplesT::const_iterator samplesIt;
    SamplesT::const_iterator samplesItE;

    m_highestTrainClassID = 0;
      
    while( classesIt != classesItE )
    {
      const ClassIDT& classID = classesIt->first;
      TERP_INSTANCE_TRUE_OR_RETURN_FALSE( classID > 0, "Invalid class ID" );
      
      const SamplesT& classSamples = classesIt->second;
      TERP_INSTANCE_TRUE_OR_RETURN_FALSE( classSamples.size() > 0,
        "Invalid class samples number" );

      m_highestTrainClassID = std::max( m_highestTrainClassID, classID );
        
      SampleT classMeans( dimsNumber );      
      
      for( dimIdx = 0 ; dimIdx < dimsNumber ; ++dimIdx )
      {
        double& dimMean = classMeans[ dimIdx ];
        dimMean = 0.0;
        
        samplesIt = classSamples.begin();
        samplesItE = classSamples.end();
          
        while( samplesIt != samplesItE )
        {
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( samplesIt->size() == dimsNumber,
            "Sample size mismatch" )
          dimMean += samplesIt->operator[]( dimIdx );
          
          ++samplesIt;
        }
        
        dimMean /= (double)( classSamples.size() );
      }
      
      m_classesMeans.push_back( classMeans );
      m_classesIndex2ID.push_back( classID );
      
      ++classesIt;
    }
    
  }
  
  // Finalizing
  
  m_initParams = (*castParamsPtr);
  
  m_isInitialized = true;
  
  return true;
}

bool te::rp::ClassifierSAMStrategy::execute() throw(te::rp::Exception)
{
  TERP_INSTANCE_TRUE_OR_RETURN_FALSE( m_isInitialized,
    "Classification strategy not initialized" );
  TERP_INSTANCE_TRUE_OR_RETURN_FALSE( m_inputRasterBands.size() > 0, "Invalid input bands" );
  TERP_INSTANCE_TRUE_OR_RETURN_FALSE( m_inputRasterBands.size() == m_classesMeans[ 0 ].size(),
    "Invalid input bands" );

  // Globals

  const unsigned int classesNumber = (unsigned int)m_classesMeans.size();
  const unsigned int nRows = (unsigned int)m_inputRasterPtr->getNumberOfRows();
  const unsigned int nCols = (unsigned int)m_inputRasterPtr->getNumberOfColumns();
  const unsigned int nDims = (unsigned int)m_classesMeans[ 0 ].size();

  // Creating the output raster

  {
    std::vector< int > dt;
    dt.push_back( te::dt::UINT32_TYPE );

    std::vector< double > noDataValues;
    noDataValues.push_back( (double)( m_highestTrainClassID + 2 ) );

    TERP_INSTANCE_TRUE_OR_RETURN_FALSE( createOutputRaster( dt, noDataValues ),
      "Output raster creation error" );

    TERP_INSTANCE_TRUE_OR_RETURN_FALSE( setOutputRasterPalette( m_highestTrainClassID + 2 ),
      "Output raster palette creation error" );
  }
    
  // progress
  
  std::unique_ptr< te::common::TaskProgress > progressPtr;
  if( m_progressInterfaceEnabled )
  {
    progressPtr.reset( new te::common::TaskProgress );
    progressPtr->setTotalSteps( m_inputRasterPtr->getNumberOfRows() );
    progressPtr->setMessage( "Classifying" );
  }    
    
  // Classifying

  unsigned int col = 0;
  unsigned int dim = 0;
  unsigned int classIdx = 0;
  double angularTR = 0;
  double angularTT = 0;
  double angularRR = 0;  
  double angularDist = 0;
  double minAngularDist = 0;
  unsigned int minAngularDistClassIdx = 0;
  const double outNoDataValue = m_outputRasterPtr->getBand( 0 )->getProperty()->m_noDataValue;
  const double inNoDataValue = m_inputRasterPtr->getBand( m_inputRasterBands[ 0 ] )->getProperty()->m_noDataValue;
  std::vector< double > currentPixelValues( nDims );
  
  for( unsigned int row = 0 ; row < nRows ; ++row )
  {
    for( col = 0 ; col < nCols ; ++col )
    {
      // Is this a valid pixel ?

      for( dim = 0 ; dim < nDims ; ++dim )
      {
        m_inputRasterPtr->getValue( col, row, currentPixelValues[ dim ],
          m_inputRasterBands[ dim ] );
      }

      if( currentPixelValues[ 0 ] == inNoDataValue )
      {
        m_outputRasterPtr->setValue( col, row, outNoDataValue, 0 );
      }
      else
      {
        // Looking the the closest class

        minAngularDist = DBL_MAX;
        minAngularDistClassIdx = classesNumber;

        for( classIdx = 0 ; classIdx < classesNumber ; ++classIdx )
        {
          const SampleT& classMeans = m_classesMeans[ classIdx ];

          angularTR = 0;
          angularTT = 0;
          angularRR = 0;

          for( dim = 0 ; dim < nDims ; ++dim )
          {
            const double& meanValue = classMeans[ dim ];

            const double& readedValue = currentPixelValues[ dim ];;

            angularTR += readedValue * meanValue;
            angularRR += meanValue * meanValue;
            angularTT += readedValue * readedValue;
          }

          angularDist = angularTR / ( sqrt( angularTT ) *
            sqrt( angularRR ) );

          if( std::abs( angularDist ) > 1.0 )
          {
            angularDist = DBL_MAX;
          }
          else
          {
            angularDist = acos( angularDist );
          }

          if( ( angularDist < minAngularDist ) && ( angularDist <
            m_initParams.m_maxAngularDistances[ classIdx ] ) )
          {
            minAngularDist = angularDist;
            minAngularDistClassIdx = classIdx;
          }
        }

        if( minAngularDistClassIdx == classesNumber )
          m_outputRasterPtr->setValue( col, row, 0, 0 );
        else
          m_outputRasterPtr->setValue( col, row, m_classesIndex2ID[ minAngularDistClassIdx ],
            0 );
      }
    }    
    
    if( m_progressInterfaceEnabled )
    {
      progressPtr->pulse();
      if( ! progressPtr->isActive() ) return false;
    }         
  }    

  return true;
}

//-----------------------------------------------------------------------------

te::rp::ClassifierSAMStrategyFactory::ClassifierSAMStrategyFactory()
  : te::rp::ClassifierStrategyFactory("sam")
{
}

te::rp::ClassifierSAMStrategyFactory::~ClassifierSAMStrategyFactory() = default;

te::rp::ClassifierStrategy* te::rp::ClassifierSAMStrategyFactory::build()
{
  return new te::rp::ClassifierSAMStrategy();
}