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

  \brief Raster ISOSeg strategy for segmentation-based classification.
*/

// TerraLib
#include "../common/MatrixUtils.h"
#include "../common/progress/TaskProgress.h"
#include "../geometry/Coord2D.h"
#include "../geometry/Envelope.h"
#include "../geometry/Point.h"
#include "../geometry/Polygon.h"
#include "../raster/Grid.h"
#include "../raster/PositionIterator.h"
#include "../raster/Utils.h"
#include "../raster/SynchronizedRaster.h"
#include "../memory/CachedRaster.h"
#include "ClassifierISOSegStrategy.h"
#include "Functions.h"
#include "Macros.h"

// STL
#include <iostream>
#include <set>
#include <stdlib.h>
#include <thread>

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

// STL
#include <cfloat>

namespace
{
  static te::rp::ClassifierISOSegStrategyFactory classifierISOSegStrategyFactoryInstance;
}

const static double classifierISOSegStrategyChiTable[10][6] =
{
//   75.0%    90.0%    95.0%    99.0%    99.9%    100.%
  { 1.32,    2.71,    3.84,    6.64,   10.83,   std::numeric_limits<double>::max()},
  { 2.77,    4.61,    5.99,    9.21,   13.82,   std::numeric_limits<double>::max()},
  { 4.11,    6.25,    7.82,   11.35,   16.27,   std::numeric_limits<double>::max()},
  { 5.39,    7.78,    9.49,   13.28,   18.47,   std::numeric_limits<double>::max()},
  { 6.63,    9.24,   11.07,   15.09,   20.52,   std::numeric_limits<double>::max()},
  { 7.84,   10.65,   12.59,   16.81,   22.46,   std::numeric_limits<double>::max()},
  { 9.04,   12.02,   14.07,   18.48,   24.32,   std::numeric_limits<double>::max()},
  {10.22,   13.36,   15.51,   20.09,   26.13,   std::numeric_limits<double>::max()},
  {11.39,   14.68,   16.92,   21.67,   27.88,   std::numeric_limits<double>::max()},
  {12.55,   15.99,   18.31,   23.21,   29.59,   std::numeric_limits<double>::max()}
};

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

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

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

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

  m_distanceType = rhs.m_distanceType;
  m_acceptanceThreshold = rhs.m_acceptanceThreshold;

  return *this;
}

void te::rp::ClassifierISOSegStrategy::Parameters::reset() throw(te::rp::Exception)
{
  m_distanceType = InvalidDistanceType;
  m_acceptanceThreshold = 0;
}

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

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

te::rp::ClassifierISOSegStrategy::Pattern::Pattern()
{
  reset();
}

te::rp::ClassifierISOSegStrategy::Pattern::Pattern(int i, double a, std::vector<double> mv, boost::numeric::ublas::matrix<double> cm)
  : m_id(i),
    m_myCluster(this),
    m_area(a),
    m_meanVector(mv),
    m_covarianceMatrix(cm)
{
  m_covarianceInversion = boost::numeric::ublas::matrix<double>(m_covarianceMatrix.size1(), m_covarianceMatrix.size2());
  
  if( ! te::common::GetDeterminant( m_covarianceMatrix, m_covarianceMatrixDet ) )
  {
    m_covarianceMatrixDet = 0;
  }
  
  te::common::GetInverseMatrix(m_covarianceMatrix, m_covarianceInversion);
}

te::rp::ClassifierISOSegStrategy::Pattern::Pattern( const Pattern& rhs)
{
  operator=( rhs );
}

te::rp::ClassifierISOSegStrategy::Pattern::~Pattern() = default;

double te::rp::ClassifierISOSegStrategy::Pattern::getMahalanobisDistance(
  Pattern const * const p ) const 
{
  assert(m_meanVector.size() == p->m_meanVector.size());

  unsigned int nBands = static_cast<unsigned int>(m_meanVector.size());

  boost::numeric::ublas::matrix<double> term1(1, nBands);
  boost::numeric::ublas::matrix<double> term2(nBands, 1);

  for (unsigned int i = 0; i < nBands; i++)
  {
    term1(0, i) = m_meanVector[i] - p->m_meanVector[i];

    term2(i, 0) = term1(0, i);
  }

  term1 = prod(term1, m_covarianceInversion);
  term1 = prod(term1, term2);

  if (term1(0, 0) < 0)
    return std::numeric_limits<double>::max();

  return term1(0, 0);
}

double te::rp::ClassifierISOSegStrategy::Pattern::getBhattacharyyaDistance(
  Pattern const * const p) const
{
  assert(m_meanVector.size() == p->m_meanVector.size());
  
  if( ( m_covarianceMatrixDet == 0 ) || ( p->m_covarianceMatrixDet == 0 ) )
  {
    return std::numeric_limits< double >::max();
  }

  m_getBhattacharyyaDistance_nBands = static_cast<unsigned int>(m_meanVector.size());

  m_getBhattacharyyaDistance_meansDif.resize(m_getBhattacharyyaDistance_nBands,
    1 );
  
  for ( m_getBhattacharyyaDistance_bandIdx = 0 ; m_getBhattacharyyaDistance_bandIdx
    < m_getBhattacharyyaDistance_nBands; ++m_getBhattacharyyaDistance_bandIdx)
  {
    m_getBhattacharyyaDistance_meansDif(m_getBhattacharyyaDistance_bandIdx, 0 ) = 
      m_meanVector[ m_getBhattacharyyaDistance_bandIdx] - 
      p->m_meanVector[m_getBhattacharyyaDistance_bandIdx];
  }
  
  m_getBhattacharyyaDistance_covsMege = ( m_covarianceMatrix
    + p->m_covarianceMatrix ) / 2.0;
    
  if( ! te::common::GetDeterminant( m_getBhattacharyyaDistance_covsMege,
    m_getBhattacharyyaDistance_covsMegeDet ) )
  {
    return std::numeric_limits< double >::max();
  }
    
  if( ! te::common::GetInverseMatrix( m_getBhattacharyyaDistance_covsMege, 
    m_getBhattacharyyaDistance_covsMegeInv ) )
  {
    return std::numeric_limits< double >::max();
  }
  
  m_getBhattacharyyaDistance_distanceTerm1 = 
    boost::numeric::ublas::prod( 
      boost::numeric::ublas::trans( m_getBhattacharyyaDistance_meansDif )
      ,
      m_getBhattacharyyaDistance_covsMegeInv
    );
  m_getBhattacharyyaDistance_distanceTerm1 = 
    boost::numeric::ublas::prod( 
      m_getBhattacharyyaDistance_distanceTerm1
      ,
      m_getBhattacharyyaDistance_meansDif
    );     
    
  m_getBhattacharyyaDistance_distanceValue = 
    m_getBhattacharyyaDistance_distanceTerm1( 0, 0 )
    /
    8.0;
  m_getBhattacharyyaDistance_distanceValue +=
    (
      0.5
      *
      std::log(
        m_getBhattacharyyaDistance_covsMegeDet
        /
        (
          std::sqrt( m_covarianceMatrixDet * p->m_covarianceMatrixDet )
        )
      )  
    );
    
  return m_getBhattacharyyaDistance_distanceValue;
}

bool te::rp::ClassifierISOSegStrategy::Pattern::operator==(const Pattern& rhs) const
{
  assert(m_myCluster == nullptr);
  assert(rhs.m_myCluster == nullptr);

  return m_id == rhs.m_id;
}

te::rp::ClassifierISOSegStrategy::Pattern& te::rp::ClassifierISOSegStrategy::Pattern::operator=( const Pattern& rhs )
{
  m_id = rhs.m_id;
  m_myCluster = rhs.m_myCluster;
  m_area = rhs.m_area;
  m_covarianceMatrixDet = rhs.m_covarianceMatrixDet;
  m_meanVector = rhs.m_meanVector;
  m_covarianceMatrix = rhs.m_covarianceMatrix;
  m_covarianceInversion = rhs.m_covarianceInversion;

  return *this;
}

bool te::rp::ClassifierISOSegStrategy::Pattern::operator<(const Pattern& rhs) const
{
  return m_area < rhs.m_area;
}

void te::rp::ClassifierISOSegStrategy::Pattern::reset()
{
  m_id = -1;
  m_myCluster = this;
  m_area = 0.0;
  m_covarianceMatrixDet = 0;
  m_meanVector.clear();
  m_covarianceMatrix.clear();
  m_covarianceInversion.clear();
}

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

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

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

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

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

  if(!paramsPtr)
    return false;

  m_parameters = *(paramsPtr);
  
  TERP_INSTANCE_TRUE_OR_RETURN_FALSE( ( m_parameters.m_distanceType ==
    Parameters::BhattacharyyaDistanceType ) || ( m_parameters.m_distanceType == 
    Parameters::MahalanobisDistanceType ),
    "Invalid distance type" );

  TERP_INSTANCE_TRUE_OR_RETURN_FALSE(m_parameters.m_acceptanceThreshold > 0 && m_parameters.m_acceptanceThreshold <= 100, "Invalid acceptance threshold [0, 100].")

  m_isInitialized = true;

  return true;
}

double te::rp::ClassifierISOSegStrategy::getThreshold(const double acceptanceThreshold, const unsigned nBands) const
{
  if(acceptanceThreshold < 90.0)
    return classifierISOSegStrategyChiTable[ std::min( nBands, 6u ) -1 ][0];
  else if(acceptanceThreshold < 95.0)
    return classifierISOSegStrategyChiTable[ std::min( nBands, 6u ) -1 ][1];
  else if(acceptanceThreshold < 99.0)
    return classifierISOSegStrategyChiTable[ std::min( nBands, 6u ) -1 ][2];
  else if(acceptanceThreshold < 99.9)
    return classifierISOSegStrategyChiTable[ std::min( nBands, 6u ) -1 ][3];
  else if(acceptanceThreshold < 100.0)
    return classifierISOSegStrategyChiTable[ std::min( nBands, 6u ) -1 ][4];
  else
    return classifierISOSegStrategyChiTable[ std::min( nBands, 6u ) -1 ][5];
}

bool te::rp::ClassifierISOSegStrategy::execute() throw(te::rp::Exception)
{
  TERP_INSTANCE_TRUE_OR_RETURN_FALSE(m_isInitialized, "Instance not initialized")
  TERP_INSTANCE_TRUE_OR_RETURN_FALSE(m_inputPolygonsPtr != nullptr, "ISOSeg algorithm needs polygons")
  TERP_INSTANCE_TRUE_OR_RETURN_FALSE(m_inputPolygonsPtr->size() > 0, "ISOSeg algorithm needs polygons")

  // globals

  // raster cache

  te::rst::Raster const* inputRasterPtr = m_inputRasterPtr;
  std::unique_ptr< te::mem::CachedRaster > cachedRasterPtr;

  if( m_enableRasterCache && ( ! m_enableMultiThread ) )
  {
    cachedRasterPtr.reset( new te::mem::CachedRaster( *m_inputRasterPtr, 25, 0 ) );
    inputRasterPtr = cachedRasterPtr.get();
  }

  // pre-processing polygons

  std::vector<te::gm::Polygon*> const * inputPolygonsPtr = m_inputPolygonsPtr;
  boost::ptr_vector< te::gm::Polygon > reprojectedPolsHandler;
  std::vector< te::gm::Polygon* > reprojectedPolsNPtrs;

  if(  inputRasterPtr->getSRID() != m_inputPolygonsPtr->operator[]( 0 )->getSRID() )
  {
    const std::size_t inputPolygonsPtrSize = m_inputPolygonsPtr->size();
    const int rasterSRID = inputRasterPtr->getSRID();
    std::unique_ptr< te::gm::Geometry > transformedGeomPtr;

    std::unique_ptr< te::common::TaskProgress > taskPtr;
    if( m_progressInterfaceEnabled )
    {
      taskPtr.reset( new te::common::TaskProgress(
        TE_TR("ISOSeg algorithm - processing polygons"), te::common::TaskProgress::UNDEFINED,
        static_cast<int>(inputPolygonsPtrSize)));
    };

    for( std::size_t inputPolygonsPtrIdx = 0 ; inputPolygonsPtrIdx <
      inputPolygonsPtrSize ; ++inputPolygonsPtrIdx )
    {
      transformedGeomPtr.reset( (te::gm::Geometry*)
        m_inputPolygonsPtr->operator[]( inputPolygonsPtrIdx )->clone() );
      transformedGeomPtr->transform( rasterSRID );
      reprojectedPolsNPtrs.push_back( (te::gm::Polygon*)transformedGeomPtr.get() );
      reprojectedPolsHandler.push_back( (te::gm::Polygon*)transformedGeomPtr.release() );

      if( m_progressInterfaceEnabled )
      {
        taskPtr->pulse();
        if( !taskPtr->isActive() ) return false;
      }
    }

    inputPolygonsPtr = &reprojectedPolsNPtrs;
  }

  // fill regions, in the beginning, each region is a cluster

  std::multimap< double, Pattern > regions; //!< A descriptive set of regions (area, features).

  {
    std::mutex globalMutex;
    std::mutex condVarMutex;
    std::condition_variable condVar;
    te::rst::RasterSynchronizer inputRasterSync( *((te::rst::Raster*)inputRasterPtr),
      te::common::RAccess );
    std::size_t nextInputPolIdx2Process = 0;
    unsigned int runningThreadsCounter = 0;
    bool continueProcessingFlag = true;

    FillRegionsThreadEntryParams threadParams;
    threadParams.m_globalMutexPtr = &globalMutex;
    threadParams.m_condVarMutexPtr = &condVarMutex;
    threadParams.m_condVarPtr = &condVar;
    threadParams.m_inputRasterSyncPtr = &inputRasterSync;
    threadParams.m_inputPolygonsPtr = inputPolygonsPtr;
    threadParams.m_inputRasterBandsPtr = &m_inputRasterBands;
    threadParams.m_outputRegionsPtr = &regions;
    threadParams.m_nextInputPolIdx2ProcessPtr = &nextInputPolIdx2Process;
    threadParams.m_runningThreadsCounterPtr = &runningThreadsCounter;
    threadParams.m_enableProgress = false;
    threadParams.m_continueProcessingFlagPtr = &continueProcessingFlag;

    if( m_enableMultiThread )
    {
      // creating threads

      const unsigned int maxSegThreads = runningThreadsCounter =
        te::common::GetPhysProcNumber();
      const std::size_t inputPolygonsSize = inputPolygonsPtr->size();
      std::vector< std::shared_ptr< std::thread > > threadsPtrs;
      std::shared_ptr< std::thread > threadPtr;

      std::unique_ptr< te::common::TaskProgress > taskPtr;
      if( m_progressInterfaceEnabled )
      {
        taskPtr.reset( new te::common::TaskProgress(
          TE_TR("ISOSeg algorithm - Extracting features"), te::common::TaskProgress::UNDEFINED,
          (int)inputPolygonsSize) );
      };

      threadParams.m_enableProgress = false;

      for( unsigned int threadIdx = 0 ; threadIdx < maxSegThreads ; ++threadIdx )
      {
        threadPtr.reset( new std::thread( fillRegionsThreadEntry, &threadParams ) );
        threadsPtrs.push_back( threadPtr );
      }

      // waiting all threads to finish

      while( true )
      {
        std::unique_lock<std::mutex> condVarLocker( condVarMutex );
        condVar.wait_for( condVarLocker, std::chrono::seconds( 1 ) );
        condVarLocker.unlock();

//        std::cout << std::endl << "Woke up" << std::endl;

        std::unique_lock<std::mutex>  globalLocker( globalMutex );

        if( m_progressInterfaceEnabled )
        {
          taskPtr->setCurrentStep( nextInputPolIdx2Process );
          if( !taskPtr->isActive() )
          {
            continueProcessingFlag = false;
            break;
          }
        }

        if( (!continueProcessingFlag) || ( runningThreadsCounter == 0 ) )
        {
          break;
        }
      }

      // joining all threads

      for( unsigned int threadIdx = 0 ; threadIdx < maxSegThreads ; ++threadIdx )
      {
        threadsPtrs[ threadIdx ]->join();
      }
    }
    else
    {
      threadParams.m_enableProgress = m_progressInterfaceEnabled;
      runningThreadsCounter = 1;
      fillRegionsThreadEntry( &threadParams );
    }
  }

    // each region is a cluster at the beginning
    // Zero determinant regions will be grouped inside the same cluster

  {
    Pattern* firstZeroDeterminantClusterPtr = 0;
    std::multimap< double, Pattern >::iterator regionsIt = regions.begin();
    const std::multimap< double, Pattern >::iterator regionsItEnd = regions.end();
    std::unique_ptr< te::common::TaskProgress > taskPtr;
    if( m_progressInterfaceEnabled )
    {
      taskPtr.reset( new te::common::TaskProgress(
        TE_TR("ISOSeg algorithm - initializing clusters"), te::common::TaskProgress::UNDEFINED,
        regions.size() ) );
    };

    while( regionsIt != regionsItEnd )
    {
      if( regionsIt->second.m_covarianceMatrixDet == 0 )
      {
        if( firstZeroDeterminantClusterPtr )
        {
          regionsIt->second.m_myCluster = firstZeroDeterminantClusterPtr;
        }
        else
        {
          firstZeroDeterminantClusterPtr = &( regionsIt->second );
          regionsIt->second.m_myCluster = &( regionsIt->second );
        }
      }
      else
      {
//         printMatrix( regionsIt->second.m_covarianceMatrix );
//         printMatrix( regionsIt->second.m_covarianceInversion );
        regionsIt->second.m_myCluster = &( regionsIt->second );
      }

      if( m_progressInterfaceEnabled )
      {
        taskPtr->pulse();
        if( ! taskPtr->isActive() ) return false;
      }

      ++regionsIt;
    }
  }

  // merge similar clusters
  // loops until all the regions are classified
  // Initiating with the greater area region

  {
    double distance = 0;
    double distance2 = 0;
    std::multimap<double, Pattern >::reverse_iterator regionsRIt1;
    std::multimap<double, Pattern >::reverse_iterator regionsRIt2;
    std::multimap<double, Pattern >::reverse_iterator regionsRIt3;
    const std::multimap<double, Pattern >::reverse_iterator regionsItRBegin = regions.rbegin();
    const std::multimap<double, Pattern >::reverse_iterator regionsItREnd = regions.rend();
    bool stable = false;
    int oldid = 0;
    std::set<std::pair<unsigned int, unsigned int> > compared;

    std::unique_ptr< te::common::TaskProgress > taskPtr;
    if( m_progressInterfaceEnabled )
    {
      taskPtr.reset( new te::common::TaskProgress(
        TE_TR("ISOSeg algorithm - detecting clusters") ) );
    };

    double maxDistance = getThreshold(m_parameters.m_acceptanceThreshold, static_cast<unsigned int>(m_inputRasterBands.size()));
    
    while (!stable)
    {
      stable = true;

      compared.clear();

      // analyze all regions

      for (regionsRIt1 = regionsItRBegin; regionsRIt1 != regionsItREnd; ++regionsRIt1)
      {
        // compare all clusters to detect if they are still similar

        for (regionsRIt2 = regionsRIt1; regionsRIt2 != regionsItREnd; ++regionsRIt2)
        {
          if (regionsRIt1->second.m_myCluster != regionsRIt2->second.m_myCluster)
          {
            if (compared.find(std::pair<unsigned int, unsigned int>(
              regionsRIt1->second.m_myCluster->m_id,
              regionsRIt2->second.m_myCluster->m_id)) == compared.end())
            {
              compared.insert(std::pair<unsigned int, unsigned int>(
                regionsRIt1->second.m_myCluster->m_id,
                regionsRIt2->second.m_myCluster->m_id));
              compared.insert(std::pair<unsigned int, unsigned int>(
                regionsRIt2->second.m_myCluster->m_id,
                regionsRIt1->second.m_myCluster->m_id));

              // use the smallest distance between region1 <-> cluster2, and region2 <-> cluster1

              switch( m_parameters.m_distanceType )
              {
                case Parameters::MahalanobisDistanceType :
                  distance = regionsRIt1->second.m_myCluster->getMahalanobisDistance(
                    regionsRIt2->second.m_myCluster );
                  distance2 = regionsRIt2->second.m_myCluster->getMahalanobisDistance(
                    regionsRIt1->second.m_myCluster );
                  distance = (distance2 < distance? distance2: distance);
                  break;
                case Parameters::BhattacharyyaDistanceType :
                  distance = regionsRIt1->second.m_myCluster->getBhattacharyyaDistance(
                    regionsRIt2->second.m_myCluster );
                  break;
                default :
                  TERP_LOG_AND_THROW( "Invalid distance" );
              }

              if (distance <= maxDistance)
              {
                // two clusters are similar and must be merged
                stable = false;

                oldid = regionsRIt2->second.m_myCluster->m_id;

                for (regionsRIt3 = regionsRIt2; regionsRIt3 != regionsItREnd; ++regionsRIt3)
                {
                  // find all clusters with riti class and merge with rit cluster
                  if (regionsRIt3->second.m_myCluster->m_id == oldid)
                  {
                    regionsRIt3->second.m_myCluster = regionsRIt1->second.m_myCluster;
                  }
                }
              }
            }
          }
        }
      }

      if( m_progressInterfaceEnabled )
      {
        taskPtr->pulse();
        if( ! taskPtr->isActive() ) return false;
      }
    }
  }

  // remap cluster values to 1 -> N

  unsigned int maxColorMapOutValue = 0;
  std::map< int, unsigned int > clusterId2OutputValueMap;

  {
    std::multimap<double, Pattern >::iterator regionsIt;
    const std::multimap<double, Pattern >::iterator regionsItEnd = regions.end();
    unsigned int outputValue = 1;

    for (regionsIt = regions.begin(); regionsIt != regionsItEnd; ++regionsIt)
    {
      if( clusterId2OutputValueMap.find( regionsIt->second.m_myCluster->m_id ) ==
        clusterId2OutputValueMap.end() )
      {
        clusterId2OutputValueMap[ regionsIt->second.m_myCluster->m_id ] = outputValue;
        ++outputValue;
      }
    }

    maxColorMapOutValue = outputValue - 1;
  }

  // Creating the output raster

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

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

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

    TERP_INSTANCE_TRUE_OR_RETURN_FALSE( setOutputRasterPalette( maxColorMapOutValue + 1 ),
      "Output raster palette creation error" );

    // no-data fill

    const double outNoDataValue = m_outputRasterPtr->getBand( 0 )->getProperty()->m_noDataValue;
    const unsigned int nCols = m_outputRasterPtr->getNumberOfColumns();
    const unsigned int nRows = m_outputRasterPtr->getNumberOfRows();
    unsigned int col = 0;
    te::rst::Band& outBand = *m_outputRasterPtr->getBand( 0 );

    for( unsigned int row = 0 ; row < nRows ; ++row )
    {
      for( col = 0 ; col < nCols ; ++col )
      {
        outBand.setValue( col, row, outNoDataValue );
      }
    }
  }

  // classify output image

  {
    const std::multimap<double, Pattern >::iterator regionsItEnd = regions.end();
    double outputValue = 0;

    std::unique_ptr< te::common::TaskProgress > taskPtr;
    if( m_progressInterfaceEnabled )
    {
      taskPtr.reset( new te::common::TaskProgress(
        TE_TR("ISOSeg algorithm - classifying image"), te::common::TaskProgress::UNDEFINED,
        static_cast<int>(regions.size())));
    };

    for (std::multimap<double, Pattern >::iterator regionsIt = regions.begin();
      regionsIt != regionsItEnd; ++regionsIt)
    {
      te::gm::Polygon* polygon = inputPolygonsPtr->operator[]( regionsIt->second.m_id );
      outputValue = clusterId2OutputValueMap[regionsIt->second.m_myCluster->m_id ];

  // iterate over polygon to classify output image
      te::rst::PolygonIterator<unsigned> it = te::rst::PolygonIterator<unsigned>::begin(
        (te::rst::Raster*)inputRasterPtr, polygon);
      te::rst::PolygonIterator<unsigned> itend = te::rst::PolygonIterator<unsigned>::end(
        (te::rst::Raster*)inputRasterPtr, polygon);

      while (it != itend)
      {
        m_outputRasterPtr->setValue(it.getColumn(), it.getRow(), outputValue, 0);

        ++it;
      }

      if( m_progressInterfaceEnabled )
      {
        taskPtr->pulse();
        if( ! taskPtr->isActive() ) return false;
      }
    }
  }

  return true;
}

void te::rp::ClassifierISOSegStrategy::fillRegionsThreadEntry(FillRegionsThreadEntryParams* paramsPtr)
{
  std::unique_lock<std::mutex>  locker( *paramsPtr->m_globalMutexPtr );

  const std::size_t inputPolygonsSize = paramsPtr->m_inputPolygonsPtr->size();
  const std::vector< unsigned int > inputRasterBands = *paramsPtr->m_inputRasterBandsPtr;

  locker.unlock();

  const std::size_t inputRasterBandsSize = inputRasterBands.size();
  const te::rst::SynchronizedRaster inputRaster(  *(paramsPtr->m_inputRasterSyncPtr),  25 );
  std::size_t inputRasterBandsIdx1 = 0;
  std::size_t inputRasterBandsIdx2 = 0;
  std::size_t currentPolIdx = 0;
  te::gm::Polygon* currentPolPtr = nullptr;
  unsigned int validPixelsCount = 0;
  std::vector< double > pixelValues( inputRasterBandsSize );
  bool isValidPixel = false;
  double covariance = 0;
  double pixelValue1 = 0;
  double pixelValue2 = 0;

  std::pair<double, Pattern> auxPair;
  auxPair.second.m_id = -1;
  auxPair.second.m_myCluster = 0;
  auxPair.second.m_area = 0;
  auxPair.second.m_covarianceMatrixDet = 0;
  auxPair.second.m_meanVector.resize( inputRasterBandsSize );
  auxPair.second.m_covarianceMatrix.resize( inputRasterBandsSize, inputRasterBandsSize );
  auxPair.second.m_covarianceInversion.resize( inputRasterBandsSize, inputRasterBandsSize );

  std::vector< double > bandsNoDataValues( inputRasterBandsSize );
  for( inputRasterBandsIdx1 = 0 ; inputRasterBandsIdx1 <
    inputRasterBandsSize ; ++inputRasterBandsIdx1 )
  {
    bandsNoDataValues[ inputRasterBandsIdx1 ] =
      inputRaster.getBand( inputRasterBands[ inputRasterBandsIdx1 ] )->getProperty()->m_noDataValue;
  }

  std::unique_ptr< te::common::TaskProgress > taskPtr;
  if( paramsPtr->m_enableProgress )
  {
    taskPtr.reset( new te::common::TaskProgress(
      TE_TR("ISOSeg algorithm - Extracting features"), te::common::TaskProgress::UNDEFINED,
      (int)inputPolygonsSize) );
  };

  while( true )
  {
    locker.lock();

    if( ! paramsPtr->m_continueProcessingFlagPtr )
    {
      --( *( paramsPtr->m_runningThreadsCounterPtr ) );
      return;
    }

    currentPolIdx = (*paramsPtr->m_nextInputPolIdx2ProcessPtr);

    if( currentPolIdx >= inputPolygonsSize )
    {
      --( *( paramsPtr->m_runningThreadsCounterPtr ) );
      return;
    }

    ++( *paramsPtr->m_nextInputPolIdx2ProcessPtr );

    currentPolPtr = paramsPtr->m_inputPolygonsPtr->operator[]( currentPolIdx );

    locker.unlock();

    // Is this a valid polygon ?

    te::rst::PolygonIterator<double> polIt = te::rst::PolygonIterator<double>::begin(
      &inputRaster, currentPolPtr );
    const te::rst::PolygonIterator<double> polItEnd = te::rst::PolygonIterator<double>::end(
      &inputRaster, currentPolPtr );

    if ( polIt != polItEnd )
    {
      auxPair.second.m_id = currentPolIdx;

      // means calcule

      for( inputRasterBandsIdx1 = 0 ; inputRasterBandsIdx1 <
        inputRasterBandsSize ; ++inputRasterBandsIdx1 )
      {
        auxPair.second.m_meanVector[ inputRasterBandsIdx1 ] = 0.0;
      }

      validPixelsCount = 0;

      while ( polIt != polItEnd )
      {
        isValidPixel = true;

        for( inputRasterBandsIdx1 = 0 ; inputRasterBandsIdx1 <
          inputRasterBandsSize ; ++inputRasterBandsIdx1 )
        {
          inputRaster.getValue( polIt.getColumn(), polIt.getRow(),
            pixelValues[ inputRasterBandsIdx1 ],  inputRasterBands[
            inputRasterBandsIdx1 ] );

          if(
              ( inputRasterBandsIdx1 == 0 )
              &&
              ( pixelValues[ 0 ] == bandsNoDataValues[ 0 ] )
            )
          {
            isValidPixel = false;
            break;
          }
        }

        if( isValidPixel )
        {
          for( inputRasterBandsIdx1 = 0 ; inputRasterBandsIdx1 <
            inputRasterBandsSize ; ++inputRasterBandsIdx1 )
          {
            auxPair.second.m_meanVector[ inputRasterBandsIdx1 ] +=
              pixelValues[ inputRasterBandsIdx1 ];
          }

          ++validPixelsCount;
        }

        ++polIt;
      }

      if( validPixelsCount )
      {
        for( inputRasterBandsIdx1 = 0 ; inputRasterBandsIdx1 <
          inputRasterBandsSize ; ++inputRasterBandsIdx1 )
        {
          auxPair.second.m_meanVector[ inputRasterBandsIdx1 ] /= (double)validPixelsCount;
        }
      }

      // Region real area

      auxPair.first = (double)validPixelsCount;
      auxPair.second.m_area = auxPair.first;

      // covariance matrix

      if( validPixelsCount > 2 )
      {
        for( inputRasterBandsIdx1 = 0 ; inputRasterBandsIdx1 <
          inputRasterBandsSize ; ++inputRasterBandsIdx1 )
        {
          for( inputRasterBandsIdx2 = 0 ; inputRasterBandsIdx2 <
            inputRasterBandsSize ; ++inputRasterBandsIdx2 )
          {
            if( inputRasterBandsIdx1 <= inputRasterBandsIdx2 )
            {
              polIt = te::rst::PolygonIterator<double>::begin(
                 &inputRaster, currentPolPtr );
              covariance = 0;

              while ( polIt != polItEnd )
              {
                inputRaster.getValue( polIt.getColumn(), polIt.getRow(),
                  pixelValue1,  inputRasterBands[ inputRasterBandsIdx1 ] );
                inputRaster.getValue( polIt.getColumn(), polIt.getRow(),
                  pixelValue2,  inputRasterBands[ inputRasterBandsIdx2 ] );

                if(
                    ( pixelValue1 != bandsNoDataValues[ inputRasterBandsIdx1 ]  )
                  )
                {
                  covariance +=
                    (
                      (
                        pixelValue1
                        -
                        auxPair.second.m_meanVector[ inputRasterBandsIdx1 ]
                      )
                      *
                      (
                        pixelValue2
                        -
                        auxPair.second.m_meanVector[ inputRasterBandsIdx2 ]
                      )
                    );
                }

                ++polIt;
              }

              covariance /= (double)( validPixelsCount - 1 );

              auxPair.second.m_covarianceMatrix( inputRasterBandsIdx1, inputRasterBandsIdx2 ) = covariance;
              auxPair.second.m_covarianceMatrix( inputRasterBandsIdx2, inputRasterBandsIdx1 ) = covariance;
            }
          }
        }

//        printMatrix( auxPair.second.m_covarianceMatrix );

        // inverse covariance matrix

        if( te::common::GetDeterminant( auxPair.second.m_covarianceMatrix, auxPair.second.m_covarianceMatrixDet ) )
        {
          if( ! te::common::GetInverseMatrix(auxPair.second.m_covarianceMatrix, auxPair.second.m_covarianceInversion) )
          {
            auxPair.second.m_covarianceMatrixDet = 0;
          }
        }
        else
        {
          auxPair.second.m_covarianceMatrixDet = 0;
        }

//        printMatrix( auxPair.second.m_covarianceInversion );
      }
      else
      {
        auxPair.second.m_covarianceMatrixDet = 0;
      }

      locker.lock();
      paramsPtr->m_outputRegionsPtr->insert( auxPair );
      locker.unlock();
    }

    if( paramsPtr->m_enableProgress )
    {
      taskPtr->setCurrentStep( currentPolIdx );
      
      if( !taskPtr->isActive() )
      {
        locker.lock();
        (*paramsPtr->m_continueProcessingFlagPtr) = false;
        --( *( paramsPtr->m_runningThreadsCounterPtr ) );
        return;
      }
    }
  }

  locker.lock();
  --( *( paramsPtr->m_runningThreadsCounterPtr ) );
}

void te::rp::ClassifierISOSegStrategy::printMatrix(const boost::numeric::ublas::matrix<double>& matrix)
{
  unsigned int size1 = matrix.size1();
  unsigned int size2 = matrix.size2();

  std::cout << std::endl << "size1=" << size1;
  std::cout << std::endl << "size2=" << size2;

  for (unsigned int i = 0; i < size1; ++ i)
  {
    std::cout << std::endl << "|";

    for (unsigned int j = 0; j < size2; ++ j)
    {
      std::cout << " " << matrix (i, j);
    }

    std::cout << " |";
  }

  std::cout << std::endl;
}

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

te::rp::ClassifierISOSegStrategyFactory::ClassifierISOSegStrategyFactory()
  : te::rp::ClassifierStrategyFactory("isoseg")
{
}

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

te::rp::ClassifierStrategy* te::rp::ClassifierISOSegStrategyFactory::build()
{
  return new te::rp::ClassifierISOSegStrategy();
}