src/terralib/rp/Contrast.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/Contrast.cpp
 
  \brief Contrast enhancement.
*/

// TerraLib
#include "../raster/Raster.h"
#include "../raster/Band.h"
#include "../raster/BandIterator.h"
#include "../raster/BandProperty.h"
#include "../raster/Grid.h"
#include "../raster/RasterProperty.h"
#include "../raster/RasterSummaryManager.h"
#include "../raster/RasterFactory.h"
#include "../memory/ExpansibleRaster.h"
#include "../common/progress/TaskProgress.h"
#include "Contrast.h"
#include "Functions.h"
#include "Macros.h"

// STL
#include <cfloat>

// Boost
#include <boost/concept_check.hpp>

namespace te
{
  namespace rp
  {
    Contrast::InputParameters::InputParameters()
    {
      reset();
    }

    Contrast::InputParameters::~InputParameters()
    {
      reset();
    }

    void Contrast::InputParameters::reset() throw( te::rp::Exception )
    {
      m_type = InputParameters::InvalidContrastT;
      
      m_lCMinInput.clear();
      m_lCMaxInput.clear();
      
      m_hECMaxInput.clear();
      
      m_squareCMinInput.clear();
      m_squareCMaxInput.clear();
      
      m_squareRootCMinInput.clear();
      m_squareRootCMaxInput.clear();
      
      m_logCMinInput.clear();
      m_logCMaxInput.clear();
      
      m_sMASCMeanInput.clear();
      m_sMASCStdInput.clear();

      m_inRasterPtr = nullptr;
      m_inRasterBands.clear();
      m_enableProgress = false;

      m_outRangeMin = -1.0 * std::numeric_limits<double>::max();
      m_outRangeMax = std::numeric_limits<double>::max();
    }

    const Contrast::InputParameters& Contrast::InputParameters::operator=(
      const Contrast::InputParameters& params )
    {
      reset();

      m_type = params.m_type;
      
      m_lCMinInput = params.m_lCMinInput;
      m_lCMaxInput = params.m_lCMaxInput;
      
      m_hECMaxInput = params.m_hECMaxInput;
      
      m_squareCMinInput = params.m_squareCMinInput;
      m_squareCMaxInput = params.m_squareCMaxInput;
      
      m_squareRootCMinInput = params.m_squareRootCMinInput;
      m_squareRootCMaxInput = params.m_squareRootCMaxInput;
      
      m_logCMinInput = params.m_logCMinInput;
      m_logCMaxInput = params.m_logCMaxInput;
      
      m_sMASCMeanInput = params.m_sMASCMeanInput;
      m_sMASCStdInput = params.m_sMASCStdInput;

      m_inRasterPtr = params.m_inRasterPtr;
      m_inRasterBands = params.m_inRasterBands;
      m_enableProgress = params.m_enableProgress;

      m_outRangeMin = params.m_outRangeMin;
      m_outRangeMax = params.m_outRangeMax;

      return *this;
    }

    te::common::AbstractParameters* Contrast::InputParameters::clone() const
    {
      return new InputParameters( *this );
    }

    Contrast::OutputParameters::OutputParameters()
    {
      reset();
    }

    Contrast::OutputParameters::OutputParameters( const OutputParameters& other )
    {
      reset();
      operator=( other );
    }

    Contrast::OutputParameters::~OutputParameters()
    {
      reset();
    }

    void Contrast::OutputParameters::reset() throw( te::rp::Exception )
    {
      m_outRasterPtr = nullptr;
      m_createdOutRasterPtr.reset();
      m_outRasterBands.clear();
      m_createdOutRasterDSType.clear();
      m_createdOutRasterInfo.clear();
    }

    const Contrast::OutputParameters& Contrast::OutputParameters::operator=(
      const Contrast::OutputParameters& params )
    {
      reset();

      m_outRasterPtr = params.m_outRasterPtr;
      m_outRasterBands = params.m_outRasterBands;
      m_createdOutRasterDSType = params.m_createdOutRasterDSType;
      m_createdOutRasterInfo = params.m_createdOutRasterInfo;

      return *this;
    }

    te::common::AbstractParameters* Contrast::OutputParameters::clone() const
    {
      return new OutputParameters( *this );
    }

    Contrast::Contrast()
    {
      reset();
    }

    Contrast::~Contrast() = default;

    bool Contrast::execute( AlgorithmOutputParameters& outputParams )
      throw( te::rp::Exception )
    {
      TERP_INSTANCE_TRUE_OR_RETURN_FALSE( m_isInitialized, "Algoritm not initialized" );

      // Initializing the output raster

      m_outputParametersPtr = dynamic_cast<
        Contrast::OutputParameters* >( &outputParams );
      TERP_INSTANCE_TRUE_OR_RETURN_FALSE( m_outputParametersPtr, "Invalid parameters" );

      if( m_outputParametersPtr->m_outRasterPtr == nullptr )
      {
        m_outputParametersPtr->m_outRasterBands.clear();

        std::vector< te::rst::BandProperty* > bandsProperties;
        for( unsigned int inRasterBandsIdx = 0 ; inRasterBandsIdx <
          m_inputParameters.m_inRasterBands.size() ; ++inRasterBandsIdx )
        {
          assert( m_inputParameters.m_inRasterBands[ inRasterBandsIdx ] <
            m_inputParameters.m_inRasterPtr->getNumberOfBands() );

          bandsProperties.push_back( new te::rst::BandProperty(
            *( m_inputParameters.m_inRasterPtr->getBand(
            m_inputParameters.m_inRasterBands[
            inRasterBandsIdx ] )->getProperty() ) ) );

          m_outputParametersPtr->m_outRasterBands.push_back( inRasterBandsIdx );
        }

        m_outputParametersPtr->m_createdOutRasterPtr.reset(
          te::rst::RasterFactory::make(
            m_outputParametersPtr->m_createdOutRasterDSType,
            new te::rst::Grid( *( m_inputParameters.m_inRasterPtr->getGrid() ) ),
            bandsProperties,
            m_outputParametersPtr->m_createdOutRasterInfo,
            nullptr,
            nullptr ) );
        TERP_INSTANCE_TRUE_OR_RETURN_FALSE( m_outputParametersPtr->m_createdOutRasterPtr.get(),
          "Output raster creation error" );

        m_outputParametersPtr->m_outRasterPtr =
          m_outputParametersPtr->m_createdOutRasterPtr.get();
      }
      else
      {
        if( ( m_outputParametersPtr->m_outRasterPtr->getAccessPolicy() & te::common::WAccess )
          && ( m_outputParametersPtr->m_outRasterPtr->getNumberOfColumns() ==
          m_inputParameters.m_inRasterPtr->getNumberOfColumns() )
          && ( m_outputParametersPtr->m_outRasterPtr->getNumberOfRows() ==
          m_inputParameters.m_inRasterPtr->getNumberOfRows() ) )
        {
          for( unsigned int inRasterBandsIdx = 0 ; inRasterBandsIdx <
            m_inputParameters.m_inRasterBands.size() ; ++inRasterBandsIdx )
          {
            TERP_INSTANCE_TRUE_OR_RETURN_FALSE(
              m_outputParametersPtr->m_outRasterBands[ inRasterBandsIdx ] <
              m_outputParametersPtr->m_outRasterPtr->getNumberOfBands(),
              "Invalid output raster band" )
          }
        }
        else
        {
          TERP_INSTANCE_LOG_AND_RETURN_FALSE( "Invalid output raster" );
        }
      }

      // Executing the contrast on the selected bands

      switch( m_inputParameters.m_type )
      {
        case InputParameters::LinearContrastT :
        {
          return execLinearContrast();
          break;
        }
        case InputParameters::HistogramEqualizationContrastT :
        {
          return execHistogramEqualizationContrast();
          break;
        }
        case InputParameters::SquareContrastT :
        {
          return execSquareContrast();
          break;
        }     
        case InputParameters::SquareRootContrastT :
        {
          return execSquareRootContrast();
          break;
        }            
        case InputParameters::LogContrastT :
        {
          return execLogContrast();
          break;
        }         
        case InputParameters::MeanAndStdContrastT :
        {
          return execSetMeanAndStdContrast();
          break;
        }
        case InputParameters::DecorrelationEnhancementT :
        {
          return execDecorrelationEnhancement();
          break;
        }        
        default :
        {
          TERP_LOG_AND_THROW( "Invalid contrast type" );
          break;
        }
      }
    }

    void Contrast::reset() throw( te::rp::Exception )
    {
      Algorithm::reset();

      m_inputParameters.reset();
      m_outputParametersPtr = nullptr;
      m_isInitialized = false;
    }

    bool Contrast::initialize( const AlgorithmInputParameters& inputParams )
      throw( te::rp::Exception )
    {
      reset();

      Contrast::InputParameters const* inputParamsPtr = dynamic_cast<
        Contrast::InputParameters const* >( &inputParams );
      TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputParamsPtr, "Invalid parameters" );

      TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputParamsPtr->m_inRasterPtr,
        "Invalid raster pointer" );
      TERP_INSTANCE_TRUE_OR_RETURN_FALSE(
        inputParamsPtr->m_inRasterPtr->getAccessPolicy() & te::common::RAccess,
        "Invalid raster" );
        
      TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputParamsPtr->m_inRasterBands.size() > 0,
        "Invalid bands number" );

      if (inputParamsPtr->m_outRangeMin != ( -1.0 * std::numeric_limits<double>::max() ) &&
        inputParamsPtr->m_outRangeMax != std::numeric_limits<double>::max())
      {
        m_inputParameters.m_outRangeMin = inputParamsPtr->m_outRangeMin;
        m_inputParameters.m_outRangeMax = inputParamsPtr->m_outRangeMax;
      }

      for( unsigned int inRasterBandsIdx = 0 ; inRasterBandsIdx <
        inputParamsPtr->m_inRasterBands.size() ; ++inRasterBandsIdx )
      {
        TERP_INSTANCE_TRUE_OR_RETURN_FALSE(
          inputParamsPtr->m_inRasterBands[ inRasterBandsIdx ] <
          inputParamsPtr->m_inRasterPtr->getNumberOfBands(),
          "Invalid input raster band" );
      }
      
      // Checking contrast specifi parameters
      
      switch( inputParamsPtr->m_type )
      {
        case InputParameters::LinearContrastT :
        {
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputParamsPtr->m_lCMinInput.size() ==
            inputParamsPtr->m_inRasterBands.size(), 
            "Invalid parameter m_lCMinInput" );
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputParamsPtr->m_lCMaxInput.size() ==
            inputParamsPtr->m_inRasterBands.size(), 
            "Invalid parameter m_lCMaxInput" );
            
          break;
        }
        case InputParameters::HistogramEqualizationContrastT :
        {
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputParamsPtr->m_hECMaxInput.size() ==
            inputParamsPtr->m_inRasterBands.size(), 
            "Invalid parameter m_hECMaxInput" );
            
          break;
        }
        case InputParameters::SquareContrastT :
        {
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputParamsPtr->m_squareCMinInput.size() ==
            inputParamsPtr->m_inRasterBands.size(), 
            "Invalid parameter m_squareCMinInput" );
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputParamsPtr->m_squareCMaxInput.size() ==
            inputParamsPtr->m_inRasterBands.size(), 
            "Invalid parameter m_squareCMaxInput" );
            
          break;
        }  
        case InputParameters::SquareRootContrastT :
        {
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputParamsPtr->m_squareRootCMinInput.size() ==
            inputParamsPtr->m_inRasterBands.size(), 
            "Invalid parameter m_squareRootCMinInput" );
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputParamsPtr->m_squareRootCMaxInput.size() ==
            inputParamsPtr->m_inRasterBands.size(), 
            "Invalid parameter m_squareRootRCMaxInput" );
            
          break;
        }         
        case InputParameters::LogContrastT :
        {
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputParamsPtr->m_logCMinInput.size() ==
            inputParamsPtr->m_inRasterBands.size(), 
            "Invalid parameter m_logCMinInput" );
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputParamsPtr->m_logCMaxInput.size() ==
            inputParamsPtr->m_inRasterBands.size(), 
            "Invalid parameter m_logCMaxInput" );
            
          break;
        }                
        case InputParameters::MeanAndStdContrastT :
        {
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputParamsPtr->m_sMASCMeanInput.size() ==
            inputParamsPtr->m_inRasterBands.size(), 
            "Invalid parameter m_sMASCMeanInput" );
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputParamsPtr->m_sMASCStdInput.size() ==
            inputParamsPtr->m_inRasterBands.size(), 
            "Invalid parameter m_sMASCStdInput" );            
          break;
        }
        case InputParameters::DecorrelationEnhancementT :
        {
          break;
        }        
        default :
        {
          TERP_LOG_AND_THROW( "Invalid contrast type" );
          break;
        }
      }

      m_inputParameters = *inputParamsPtr;
      m_isInitialized = true;

      return true;
    }

    bool Contrast::isInitialized() const
    {
      return m_isInitialized;
    }

    bool Contrast::getGainAndOffset( const InputParameters::ContrastType& type,
      const double& inRangeMin, const double& inRangeMax,
      const double& outRangeMin, const double& outRangeMax,
      double& gain, double& offset1, double& offset2 )
    {
      const double outRange = outRangeMax - outRangeMin;
      const double inRange = inRangeMax - inRangeMin;
      TERP_TRUE_OR_RETURN_FALSE( outRange >= 0.0, "Invalid values range" );
      TERP_TRUE_OR_RETURN_FALSE( inRange >= 0.0, "Invalid values range" );

      switch( type )
      {
        case InputParameters::LinearContrastT :
        {
          if( inRange == 0.0 )
          {
            gain = 1.0;
          }
          else
          {
            gain = outRange / inRange;
          }

          offset1 = (-1.0) * inRangeMin;
          offset2 = outRangeMin;

          break;
        }
        case InputParameters::SquareContrastT :
        {
          if( inRange == 0.0 )
          {
            gain = 0.0;
          }
          else
          {
            gain = outRange / std::pow( inRange, 2.0 );
          }

          offset1 = 0.0;
          offset2 = 0.0;

          break;
        }
        case InputParameters::SquareRootContrastT :
        {
          if( inRange == 0.0 )
          {
            gain = 0.0;
          }
          else
          {
            gain = outRange / std::sqrt( inRange );
          }

          offset1 = 0.0;
          offset2 = 0.0;

          break;
        }
        case InputParameters::LogContrastT :
        {
          if( inRange == 0.0 )
          {
            gain = 0.0;
          }
          else
          {
            gain = outRange / std::log10( inRange + 1.0 );
          }

          offset1 = -1.0 * inRangeMin;
          offset2 = 0.0;

          break;
        }
        default :
        {
          TERP_LOG_AND_RETURN_FALSE( "Invalid contrast type" );
          break;
        }
      }

      return true;
    }

    bool Contrast::execLinearContrast()
    {
      assert( m_inputParameters.m_inRasterPtr );
      assert( m_outputParametersPtr->m_outRasterPtr );
      
      const bool enableGlobalProgress = m_inputParameters.m_enableProgress &&
         ( m_inputParameters.m_inRasterBands.size() > 1 );
      std::unique_ptr< te::common::TaskProgress > progressPtr;
      if( enableGlobalProgress )
      {
        progressPtr.reset( new te::common::TaskProgress );
        progressPtr->setMessage( "Contrast" );
        progressPtr->setTotalSteps( static_cast<int>(m_inputParameters.m_inRasterBands.size()) );
      }

      for( unsigned int inRasterBandsIdx = 0 ; inRasterBandsIdx <
        m_inputParameters.m_inRasterBands.size() ; ++inRasterBandsIdx )
      {
        te::rst::Band const* inRasterBandPtr = m_inputParameters.m_inRasterPtr->getBand(
          m_inputParameters.m_inRasterBands[ inRasterBandsIdx ] );
        te::rst::Band* outRasterBandPtr =m_outputParametersPtr->m_outRasterPtr->getBand(
         m_outputParametersPtr->m_outRasterBands[ inRasterBandsIdx ] );


        if (m_inputParameters.m_outRangeMin == ( -1.0 * std::numeric_limits<double>::max() ) &&
          m_inputParameters.m_outRangeMax == std::numeric_limits<double>::max())
        {
          GetDataTypeRange(outRasterBandPtr->getProperty()->getType(),
            m_inputParameters.m_outRangeMin, m_inputParameters.m_outRangeMax);
        }
        
        if( ! getGainAndOffset( InputParameters::LinearContrastT,
            m_inputParameters.m_lCMinInput[ inRasterBandsIdx ],
            m_inputParameters.m_lCMaxInput[ inRasterBandsIdx ],
            m_inputParameters.m_outRangeMin,
            m_inputParameters.m_outRangeMax,
            m_offSetGainRemap_gain,
            m_offSetGainRemap_offset1,
            m_offSetGainRemap_offset2 )
          )
        {
          return false;
        }

        if( ! remapBandLevels( *inRasterBandPtr, *outRasterBandPtr,
          &Contrast::offSetGainRemap, m_inputParameters.m_enableProgress &&
          (!enableGlobalProgress) ) )
        {
          return false;
        }
        
        if( enableGlobalProgress )
        {
          progressPtr->pulse();
          if( ! progressPtr->isActive() ) return false;
        }        
      }

      return true;
    }

    bool Contrast::execHistogramEqualizationContrast()
    {
      assert( m_inputParameters.m_inRasterPtr );
      assert(m_outputParametersPtr->m_outRasterPtr );

      te::rst::RasterSummary const * const rasterSummaryPtr =
        te::rst::RasterSummaryManager::getInstance().get(
        m_inputParameters.m_inRasterPtr, te::rst::SUMMARY_R_HISTOGRAM, true);

      double newvalue = 0;
      unsigned int totalPixelsNumber = m_inputParameters.m_inRasterPtr->getNumberOfRows() *
        m_inputParameters.m_inRasterPtr->getNumberOfColumns();
        
      const bool enableGlobalProgress = m_inputParameters.m_enableProgress &&
         ( m_inputParameters.m_inRasterBands.size() > 1 );
      std::unique_ptr< te::common::TaskProgress > progressPtr;
      if( enableGlobalProgress )
      {
        progressPtr.reset( new te::common::TaskProgress );
        progressPtr->setMessage( "Contrast" );
        progressPtr->setTotalSteps( static_cast<int>(m_inputParameters.m_inRasterBands.size()) );
      }

      for( unsigned int b = 0 ; b < m_inputParameters.m_inRasterBands.size() ; b++ )
      {
        const unsigned int inputBandIdx = m_inputParameters.m_inRasterBands[b];
        const unsigned int outputBandIdx = m_outputParametersPtr->m_outRasterBands[b];
        const te::rst::Band& inputBand = *m_inputParameters.m_inRasterPtr->getBand(inputBandIdx);
        te::rst::Band& outputBand = *m_outputParametersPtr->m_outRasterPtr->getBand(outputBandIdx);

        // lut

//         const double lutfactor =
//           m_inputParameters.m_hECMaxInput[ b ] / (double) totalPixelsNumber;
        const int inputDataType = inputBand.getProperty()->getType();
        const int outputDataType = outputBand.getProperty()->getType();

        if (m_inputParameters.m_outRangeMin == ( -1.0 * std::numeric_limits<double>::max() ) &&
          m_inputParameters.m_outRangeMax == std::numeric_limits<double>::max())
        {
          GetDataTypeRange(outputDataType, m_inputParameters.m_outRangeMin,
            m_inputParameters.m_outRangeMax);
        }

        const double minp = rasterSummaryPtr->at(inputBandIdx).m_minVal->real();
        const double maxp = rasterSummaryPtr->at(inputBandIdx).m_maxVal->real();

        const bool isRealValuesBand = ( inputDataType == te::dt::FLOAT_TYPE ) ||
            ( inputDataType == te::dt::DOUBLE_TYPE );

        const double pixelValueStep = isRealValuesBand ? ( ( maxp - minp ) /
          ((double)te::rst::RasterSummaryManager::getInstance().getRealValuesRasterHistSize() ) )
          : 1.0;

        // cdf

        m_histogramEqualizationRemap_cdf.clear();

        if( isRealValuesBand )
        {
          for (double pixel = minp; pixel <= maxp; pixel += pixelValueStep )
          {
            newvalue = 0.0;

            for (double nj = minp; nj <= pixel; nj++)
            {
              newvalue += rasterSummaryPtr->at(inputBandIdx).m_histogramR->lower_bound( nj )->second;
            }

            m_histogramEqualizationRemap_cdf[pixel] = newvalue;
          }
        }
        else
        {
          for (double pixel = minp; pixel <= maxp; pixel += pixelValueStep )
          {
            newvalue = 0.0;

            for (double nj = minp; nj <= pixel; nj++)
            {
              newvalue += rasterSummaryPtr->at(inputBandIdx).m_histogramR->operator[](nj);
            }

            m_histogramEqualizationRemap_cdf[pixel] = newvalue;
          }
        }

        m_histogramEqualizationRemap_cDFMin = m_histogramEqualizationRemap_cdf.begin()->second;

        m_histogramEqualizationRemap_gain = m_inputParameters.m_hECMaxInput[ b ]
          / ((double)totalPixelsNumber );

        if( ! remapBandLevels( inputBand, outputBand,
          &Contrast::histogramEqualizationRemap, m_inputParameters.m_enableProgress &&
          (!enableGlobalProgress) ) )
        {
          return false;
        }

        if( enableGlobalProgress )
        {
          progressPtr->pulse();
          if( ! progressPtr->isActive() ) return false;
        }
      }

      return true;
    }
    
    bool Contrast::execSquareContrast()
    {
      assert( m_inputParameters.m_inRasterPtr );
      assert( m_outputParametersPtr->m_outRasterPtr );
      
      const bool enableGlobalProgress = m_inputParameters.m_enableProgress &&
         ( m_inputParameters.m_inRasterBands.size() > 1 );
      std::unique_ptr< te::common::TaskProgress > progressPtr;
      if( enableGlobalProgress )
      {
        progressPtr.reset( new te::common::TaskProgress );
        progressPtr->setMessage( "Contrast" );
        progressPtr->setTotalSteps( static_cast<int>(m_inputParameters.m_inRasterBands.size()) );
      }

      for( unsigned int inRasterBandsIdx = 0 ; inRasterBandsIdx <
        m_inputParameters.m_inRasterBands.size() ; ++inRasterBandsIdx )
      {
        te::rst::Band const* inRasterBandPtr = m_inputParameters.m_inRasterPtr->getBand(
          m_inputParameters.m_inRasterBands[ inRasterBandsIdx ] );
        te::rst::Band* outRasterBandPtr =m_outputParametersPtr->m_outRasterPtr->getBand(
         m_outputParametersPtr->m_outRasterBands[ inRasterBandsIdx ] );

        if (m_inputParameters.m_outRangeMin == ( -1.0 * std::numeric_limits<double>::max() ) &&
          m_inputParameters.m_outRangeMax == std::numeric_limits<double>::max())
        {
          GetDataTypeRange(outRasterBandPtr->getProperty()->getType(),
            m_inputParameters.m_outRangeMin, m_inputParameters.m_outRangeMax);
        }

        if( ! getGainAndOffset( InputParameters::SquareContrastT,
            m_inputParameters.m_squareCMinInput[ inRasterBandsIdx ],
            m_inputParameters.m_squareCMaxInput[ inRasterBandsIdx ],
            m_inputParameters.m_outRangeMin,
            m_inputParameters.m_outRangeMax,
            m_squareRemap_gain,
            m_offSetGainRemap_offset1,
            m_offSetGainRemap_offset2 )
          )
        {
          return false;
        }

        if( ! remapBandLevels( *inRasterBandPtr, *outRasterBandPtr,
          &Contrast::squareRemap, m_inputParameters.m_enableProgress &&
          (!enableGlobalProgress) ) )
        {
          return false;
        }
        
        if( enableGlobalProgress )
        {
          progressPtr->pulse();
          if( ! progressPtr->isActive() ) return false;
        }        
      }

      return true;
    }    

    bool Contrast::execSquareRootContrast()
    {
      assert( m_inputParameters.m_inRasterPtr );
      assert( m_outputParametersPtr->m_outRasterPtr );
      
      const bool enableGlobalProgress = m_inputParameters.m_enableProgress &&
         ( m_inputParameters.m_inRasterBands.size() > 1 );
      std::unique_ptr< te::common::TaskProgress > progressPtr;
      if( enableGlobalProgress )
      {
        progressPtr.reset( new te::common::TaskProgress );
        progressPtr->setMessage( "Contrast" );
        progressPtr->setTotalSteps( static_cast<int>(m_inputParameters.m_inRasterBands.size()) );
      }

      for( unsigned int inRasterBandsIdx = 0 ; inRasterBandsIdx <
        m_inputParameters.m_inRasterBands.size() ; ++inRasterBandsIdx )
      {
        te::rst::Band const* inRasterBandPtr = m_inputParameters.m_inRasterPtr->getBand(
          m_inputParameters.m_inRasterBands[ inRasterBandsIdx ] );
        te::rst::Band* outRasterBandPtr =m_outputParametersPtr->m_outRasterPtr->getBand(
         m_outputParametersPtr->m_outRasterBands[ inRasterBandsIdx ] );

        if (m_inputParameters.m_outRangeMin == ( -1.0 * std::numeric_limits<double>::max() ) &&
          m_inputParameters.m_outRangeMax == std::numeric_limits<double>::max())
        {
          GetDataTypeRange(outRasterBandPtr->getProperty()->getType(),
            m_inputParameters.m_outRangeMin, m_inputParameters.m_outRangeMax);
        }

        if( ! getGainAndOffset( InputParameters::SquareRootContrastT,
            m_inputParameters.m_squareRootCMinInput[ inRasterBandsIdx ],
            m_inputParameters.m_squareRootCMaxInput[ inRasterBandsIdx ],
            m_inputParameters.m_outRangeMin,
            m_inputParameters.m_outRangeMax,
            m_squareRootRemap_gain,
            m_offSetGainRemap_offset1,
            m_offSetGainRemap_offset2 )
          )
        {
          return false;
        }

        if( ! remapBandLevels( *inRasterBandPtr, *outRasterBandPtr,
          &Contrast::squareRootRemap, m_inputParameters.m_enableProgress &&
          (!enableGlobalProgress) ) )
        {
          return false;
        }
        
        if( enableGlobalProgress )
        {
          progressPtr->pulse();
          if( ! progressPtr->isActive() ) return false;
        }        
      }

      return true;
    }
    
    bool Contrast::execLogContrast()
    {
      assert( m_inputParameters.m_inRasterPtr );
      assert( m_outputParametersPtr->m_outRasterPtr );
      
      const bool enableGlobalProgress = m_inputParameters.m_enableProgress &&
         ( m_inputParameters.m_inRasterBands.size() > 1 );
      std::unique_ptr< te::common::TaskProgress > progressPtr;
      if( enableGlobalProgress )
      {
        progressPtr.reset( new te::common::TaskProgress );
        progressPtr->setMessage( "Contrast" );
        progressPtr->setTotalSteps( static_cast<int>(m_inputParameters.m_inRasterBands.size()) );
      }

      for( unsigned int inRasterBandsIdx = 0 ; inRasterBandsIdx <
        m_inputParameters.m_inRasterBands.size() ; ++inRasterBandsIdx )
      {
        te::rst::Band const* inRasterBandPtr = m_inputParameters.m_inRasterPtr->getBand(
          m_inputParameters.m_inRasterBands[ inRasterBandsIdx ] );
        te::rst::Band* outRasterBandPtr =m_outputParametersPtr->m_outRasterPtr->getBand(
         m_outputParametersPtr->m_outRasterBands[ inRasterBandsIdx ] );

        if (m_inputParameters.m_outRangeMin == ( -1.0 * std::numeric_limits<double>::max() ) &&
          m_inputParameters.m_outRangeMax == std::numeric_limits<double>::max())
        {
          GetDataTypeRange(outRasterBandPtr->getProperty()->getType(),
            m_inputParameters.m_outRangeMin, m_inputParameters.m_outRangeMax);
        }

        if( ! getGainAndOffset( InputParameters::LogContrastT,
            m_inputParameters.m_logCMinInput[ inRasterBandsIdx ],
            m_inputParameters.m_logCMaxInput[ inRasterBandsIdx ],
            m_inputParameters.m_outRangeMin,
            m_inputParameters.m_outRangeMax,
            m_logRemap_gain,
            m_logRemap_offset,
            m_offSetGainRemap_offset2 )
          )
        {
          return false;
        }

        if( ! remapBandLevels( *inRasterBandPtr, *outRasterBandPtr,
          &Contrast::logRemap, m_inputParameters.m_enableProgress &&
          (!enableGlobalProgress) ) )
        {
          return false;
        }
        
        if( enableGlobalProgress )
        {
          progressPtr->pulse();
          if( ! progressPtr->isActive() ) return false;
        }        
      }

      return true;
    }    
    
    bool Contrast::execSetMeanAndStdContrast()
    {
      assert( m_inputParameters.m_inRasterPtr );
      assert(m_outputParametersPtr->m_outRasterPtr );

      const te::rst::RasterSummary* rsummary =
        te::rst::RasterSummaryManager::getInstance().get(
        m_inputParameters.m_inRasterPtr,
        (te::rst::SummaryTypes) (te::rst::SUMMARY_MEAN | te::rst::SUMMARY_STD));
        
      std::unique_ptr< te::common::TaskProgress > progressPtr;
      if( m_inputParameters.m_enableProgress )
      {
        progressPtr.reset( new te::common::TaskProgress );
        progressPtr->setMessage( "Contrast" );
        progressPtr->setTotalSteps( static_cast<int>(m_inputParameters.m_inRasterBands.size() *
           m_inputParameters.m_inRasterPtr->getNumberOfRows()) );
      }          

      for( unsigned int b = 0 ; b < m_inputParameters.m_inRasterBands.size() ; b++ )
      {
        unsigned int niband = m_inputParameters.m_inRasterBands[b];
        unsigned int noband = m_outputParametersPtr->m_outRasterBands[b];

        const te::rst::Band* iband = m_inputParameters.m_inRasterPtr->getBand(niband);
        te::rst::Band* oband = m_outputParametersPtr->m_outRasterPtr->getBand(noband);
        
        if (m_inputParameters.m_outRangeMin == ( -1.0 * std::numeric_limits<double>::max() ) &&
          m_inputParameters.m_outRangeMax == std::numeric_limits<double>::max())
        {
          GetDataTypeRange(oband->getProperty()->getType(),
            m_inputParameters.m_outRangeMin, m_inputParameters.m_outRangeMax);
        }

        te::rst::ConstBandIterator<double> ibandit =
          te::rst::ConstBandIterator<double>::begin(iband);
        te::rst::ConstBandIterator<double> ibanditend =
          te::rst::ConstBandIterator<double>::end(iband);

        const double meanp = rsummary->at(niband).m_meanVal->real();
        const double stdp = rsummary->at(niband).m_stdVal->real();
        const double offset = m_inputParameters.m_sMASCMeanInput[ b ] /
          m_inputParameters.m_sMASCStdInput[ b ];
        const double c1 = m_inputParameters.m_sMASCStdInput[ b ] / stdp;
        const double c2 = offset * m_inputParameters.m_sMASCStdInput[ b ];

        double value;
        double newvalue;

        for (unsigned r = 0; r < iband->getRaster()->getNumberOfRows(); r++)
        {
          for (unsigned c = 0; c < iband->getRaster()->getNumberOfColumns(); c++)
          {
            iband->getValue(c, r, value);
            newvalue = (value - meanp) * c1 + c2;
            newvalue = std::max(m_inputParameters.m_outRangeMin, newvalue);
            newvalue = std::min(m_inputParameters.m_outRangeMax, newvalue);
            oband->setValue(c, r, newvalue);
          }
          
          if( m_inputParameters.m_enableProgress )
          {
            progressPtr->pulse();
            if( ! progressPtr->isActive() ) return false;
          }            
        }
      }

      return true;
    }
    
    bool Contrast::execDecorrelationEnhancement()
    {
      assert( m_inputParameters.m_inRasterPtr );
      assert( m_outputParametersPtr->m_outRasterPtr );

      // Progress interface
      
      std::unique_ptr< te::common::TaskProgress > progressPtr;
      if( m_inputParameters.m_enableProgress )
      {
        progressPtr.reset( new te::common::TaskProgress );
        progressPtr->setMessage( "Contrast" );
        progressPtr->setTotalSteps( static_cast<int>(m_inputParameters.m_inRasterBands.size() +
          2) );
      }
      
      // Compute PCA space raster
      
      std::unique_ptr< te::rst::Raster > pcaRasterPtr;
      boost::numeric::ublas::matrix< double > pcaMatrix;
      
      {
        te::rst::Grid* gridPtr = new te::rst::Grid( 
          *m_inputParameters.m_inRasterPtr->getGrid() );
        std::vector< te::rst::BandProperty * > bandProperties;
        std::vector< unsigned int > pcaRasterBands;
          
        for( unsigned int bandIdx = 0 ; bandIdx <
          m_inputParameters.m_inRasterBands.size() ; ++bandIdx )
        {
          bandProperties.push_back( new te::rst::BandProperty( 
            *m_inputParameters.m_inRasterPtr->getBand( bandIdx )->getProperty() ) );
          bandProperties[ bandIdx ]->m_type = te::dt::DOUBLE_TYPE;
          
          pcaRasterBands.push_back( bandIdx );
        } 
        
        try
        {
          pcaRasterPtr.reset( new te::mem::ExpansibleRaster( 25, gridPtr,
            bandProperties ) );    
        }
        catch( te::common::Exception& e )
        {
          TERP_LOG_AND_RETURN_FALSE( e.what() );
        }    
        
        TERP_TRUE_OR_RETURN_FALSE( te::rp::DirectPrincipalComponents( 
          *m_inputParameters.m_inRasterPtr, m_inputParameters.m_inRasterBands,
          pcaMatrix, *pcaRasterPtr, pcaRasterBands, 0 ),
          "Principal components generation error" );
          
//        te::rp::Copy2DiskRaster( *pcaRasterPtr, "pcaRaster.tif" );
      }      
      
      if( m_inputParameters.m_enableProgress )
      {
        progressPtr->pulse();
        if( ! progressPtr->isActive() ) return false;
      }          
      
      // The reference for enhancement: mean and std dev from the first PCA band
      
      double pca0MeanValue = 0;
      TERP_TRUE_OR_RETURN_FALSE( te::rp::GetMeanValue( 
        *pcaRasterPtr->getBand( 0 ), 0, pca0MeanValue ),
        "Mean value calcule error" );
      
      double pca0StdDevValue = 0;
      TERP_TRUE_OR_RETURN_FALSE( te::rp::GetStdDevValue( 
        *pcaRasterPtr->getBand( 0 ), 0, &pca0MeanValue,
        pca0StdDevValue ), "StdDev value calcule error" );      
      
      // Enhancing PCA

      for( unsigned int pcaBandIdx = 1 ; pcaBandIdx <
        m_inputParameters.m_inRasterBands.size() ; ++pcaBandIdx )
      {
        te::rst::Band& pcaBand = *pcaRasterPtr->getBand( pcaBandIdx );

        double meanValue = 0;
        TERP_TRUE_OR_RETURN_FALSE( te::rp::GetMeanValue( pcaBand, 0, meanValue ),
          "Mean value calcule error" );
        
        double stdDevValue = 0;
        TERP_TRUE_OR_RETURN_FALSE( te::rp::GetStdDevValue( pcaBand, 0, &meanValue,
          stdDevValue ), "StdDev value calcule error" );           

        m_decorrelationEnhancementRemap_offset1 = -1.0 * meanValue;
        m_decorrelationEnhancementRemap_offset2 = pca0MeanValue;
        m_decorrelationEnhancementRemap_gain = ( stdDevValue == 0.0 ) ? 0.0 : ( pca0StdDevValue /
          stdDevValue );

        if( ! remapBandLevels( pcaBand, pcaBand,
          &Contrast::DecorrelationEnhancementRemap, 0 ) )
        {
          return false;
        }
        
        if( m_inputParameters.m_enableProgress )
        {
          progressPtr->pulse();
          if( ! progressPtr->isActive() ) return false;
        }        
      }
      
      // Generating the output enhanced raster
      
      TERP_TRUE_OR_RETURN_FALSE( InversePrincipalComponents( *pcaRasterPtr,
        pcaMatrix, *m_outputParametersPtr->m_outRasterPtr, 
        m_outputParametersPtr->m_outRasterBands, 0 ),
        "Inverse PCA error" );
      
      if( m_inputParameters.m_enableProgress )
      {
        progressPtr->pulse();
        if( ! progressPtr->isActive() ) return false;
      }       

      return true;
    }    

    bool Contrast::remapBandLevels( const te::rst::Band& inRasterBand,
      te::rst::Band& outRasterBand, RemapFuncPtrT remapFuncPtr,
      const bool enableProgress )
    {
      const unsigned int inBlkWidthPixels = inRasterBand.getProperty()->m_blkw;
      const unsigned int inBlkHeightPixels = inRasterBand.getProperty()->m_blkh;
      const unsigned int outBlkWidthPixels = outRasterBand.getProperty()->m_blkw;
      const unsigned int outBlkHeightPixels = outRasterBand.getProperty()->m_blkh;
      
      if (m_inputParameters.m_outRangeMin == ( -1.0 * std::numeric_limits<double>::max() ) &&
        m_inputParameters.m_outRangeMax == std::numeric_limits<double>::max())
      {
        GetDataTypeRange(outRasterBand.getProperty()->getType(),
          m_inputParameters.m_outRangeMin, m_inputParameters.m_outRangeMax);
      }

      std::unique_ptr< te::common::TaskProgress > progressPtr;
      if( enableProgress )
      {
        progressPtr.reset( new te::common::TaskProgress );
        progressPtr->setMessage( "Contrast" );
      }          

      if( ( inBlkWidthPixels == outBlkWidthPixels ) &&
        ( inBlkHeightPixels == outBlkHeightPixels ) )
      { // block version
        const unsigned int blockSizePixels = inBlkWidthPixels *
          inBlkHeightPixels;
        const unsigned int xBlocksNmb = inRasterBand.getProperty()->m_nblocksx;
        const unsigned int yBlocksNmb = inRasterBand.getProperty()->m_nblocksy;
        unsigned char* inputBuffer = new unsigned char[ inRasterBand.getBlockSize() ];
        unsigned char* outputBuffer = new unsigned char[ outRasterBand.getBlockSize() ];
        double* inDoublesBuffer = new double[ blockSizePixels ];
        double* outDoublesBuffer = new double[ blockSizePixels ];
        const int inputVectorDataType = inRasterBand.getProperty()->getType();
        const int outputVectorDataType = outRasterBand.getProperty()->getType();
        unsigned int blockXIndex = 0;
        unsigned int blockOffset = 0;
        double outValue = 0;
        
        if( enableProgress ) progressPtr->setTotalSteps( yBlocksNmb * xBlocksNmb );

        for( unsigned int blockYIndex = 0 ; blockYIndex < yBlocksNmb ;
          ++blockYIndex )
        {
          for( blockXIndex = 0 ; blockXIndex < xBlocksNmb ;
            ++blockXIndex )
          {
            inRasterBand.read( blockXIndex, blockYIndex, inputBuffer );

            Convert2DoublesVector( inputBuffer, inputVectorDataType,
              blockSizePixels, inDoublesBuffer );

            for( blockOffset = 0 ; blockOffset < blockSizePixels ; ++blockOffset )
            {
              (this->*remapFuncPtr)( inDoublesBuffer[blockOffset],
                outValue);
              outDoublesBuffer[blockOffset] = std::max(m_inputParameters.m_outRangeMin, std::min(
                m_inputParameters.m_outRangeMax, outValue));
            }

            ConvertDoublesVector( outDoublesBuffer, blockSizePixels,
              outputVectorDataType, outputBuffer );

            outRasterBand.write( blockXIndex, blockYIndex, outputBuffer );
            
            if( enableProgress )
            {
              progressPtr->pulse();
              if( ! progressPtr->isActive() ) return false;
            }
          }
        }

        delete[] inputBuffer;
        delete[] outputBuffer;
        delete[] inDoublesBuffer;
        delete[] outDoublesBuffer;
      }
      else
      { // pixel by pixel version
        const unsigned int linesNumber = m_inputParameters.m_inRasterPtr->getNumberOfRows();
        const unsigned int columnsNumber =  m_inputParameters.m_inRasterPtr->getNumberOfColumns();
        
        if( enableProgress ) progressPtr->setTotalSteps( linesNumber );

        unsigned int col = 0;
        double inputValue = 0;
        double outputValue = 0;

        for( unsigned int line = 0 ; line < linesNumber ; ++line )
        {
          for( col = 0 ; col < columnsNumber ; ++col )
          {
            inRasterBand.getValue( col, line, inputValue );
            (this->*remapFuncPtr)( inputValue, outputValue );
            outRasterBand.setValue(col, line, std::max(m_inputParameters.m_outRangeMin, std::min(
              m_inputParameters.m_outRangeMax, outputValue)));
          }
          
          if( enableProgress )
          {
            progressPtr->pulse();
            if( ! progressPtr->isActive() ) return false;
          }          
        }
      }

      return true;
    }

  } // end namespace rp
}   // end namespace te