GeoMosaic.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/GeoMosaic.cpp
  \brief Create a mosaic from a set of geo-referenced rasters.
*/

#include "GeoMosaic.h"

#include "Macros.h"
#include "Functions.h"
#include "../raster/Interpolator.h"
#include "../raster/Enums.h"
#include "../raster/RasterFactory.h"
#include "../raster/Grid.h"
#include "../raster/Band.h"
#include "../raster/BandProperty.h"
#include "../raster/PositionIterator.h"
#include "../raster/Utils.h"
#include "../memory/CachedRaster.h"
#include "../geometry/Envelope.h"
#include "../geometry/GTFactory.h"
#include "../geometry/Polygon.h"
#include "../geometry/LinearRing.h"
#include "../geometry/MultiPolygon.h"
#include "../srs/Converter.h"
#include "../common/progress/TaskProgress.h"
#include "../common/MathUtils.h"

#include <boost/shared_array.hpp>
#include <boost/filesystem.hpp>
#include <boost/lexical_cast.hpp>

#include <climits>
#include <cfloat>
#include <cmath>
#include <memory>
#include <cstdio>

namespace te
{
  namespace rp
  {

    GeoMosaic::InputParameters::InputParameters()
    {
      reset();
    }

    GeoMosaic::InputParameters::InputParameters( const InputParameters& other )
    {
      reset();
      operator=( other );
    }

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

    void GeoMosaic::InputParameters::reset() _NOEXCEPT_OP(false)
    {
      m_feederRasterPtr = nullptr;
      m_inputRastersBands.clear();
      m_interpMethod = te::rst::NearestNeighbor;
      m_noDataValue = 0.0;
      m_forceInputNoDataValue = false;
      m_blendMethod = te::rp::Blender::NoBlendMethod;
      m_autoEqualize = true;
      m_useRasterCache = true;
      m_enableProgress = false;
      m_enableMultiThread = true;
    }

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

      m_feederRasterPtr = params.m_feederRasterPtr;
      m_inputRastersBands = params.m_inputRastersBands;
      m_interpMethod = params.m_interpMethod;
      m_noDataValue = params.m_noDataValue;
      m_forceInputNoDataValue = params.m_forceInputNoDataValue;
      m_blendMethod = params.m_blendMethod;
      m_autoEqualize = params.m_autoEqualize;
      m_useRasterCache = params.m_useRasterCache;
      m_enableProgress = params.m_enableProgress;
      m_enableMultiThread = params.m_enableMultiThread;

      return *this;
    }

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

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

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

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

    void GeoMosaic::OutputParameters::reset() _NOEXCEPT_OP(false)
    {
      m_rType.clear();
      m_rInfo.clear();
      m_outputRasterPtr.reset();
    }

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

      m_rType = params.m_rType;
      m_rInfo = params.m_rInfo;
      m_outputRasterPtr.reset();;

      return *this;
    }

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

    GeoMosaic::GeoMosaic()
    {
      reset();
    }

    GeoMosaic::~GeoMosaic() = default;

    bool GeoMosaic::execute( AlgorithmOutputParameters& outputParams )
      _NOEXCEPT_OP(false)
    {
      if( ! m_isInitialized ) return false;
      
      GeoMosaic::OutputParameters* outParamsPtr = dynamic_cast<
        GeoMosaic::OutputParameters* >( &outputParams );
      TERP_TRUE_OR_THROW( outParamsPtr, "Invalid paramters" );
      
      assert( m_inputParameters.m_feederRasterPtr->getObjsCount() > 1 );
      
      // progress
      
      std::unique_ptr< te::common::TaskProgress > progressPtr;
      if( m_inputParameters.m_enableProgress )
      {
        progressPtr.reset( new te::common::TaskProgress );
        
        progressPtr->setTotalSteps( static_cast<int>(2 + m_inputParameters.m_feederRasterPtr->getObjsCount()) );
        
        progressPtr->setMessage( "Mosaic" );
      }       
      
      // First pass: getting global mosaic info

      double mosaicXResolution = 0.0;
      double mosaicYResolution = 0.0;
      double mosaicLLX = DBL_MAX; // world coords
      double mosaicLLY = DBL_MAX; // world coords
      double mosaicURX = -1.0 * DBL_MAX; // world coords
      double mosaicURY = -1.0 * DBL_MAX; // world coords
      int mosaicSRID = 0;     
      te::rst::BandProperty mosaicBaseBandProperties( 0, 0, "" );
      std::vector< te::gm::Polygon > rastersBBoxes; // all rasters bounding boxes (under the first raster world coords.

      {
        te::rst::Raster const* inputRasterPtr = nullptr;
        unsigned int inputRasterIdx = 0;
        std::unique_ptr< te::gm::Polygon > rasterExtentPolPtr;

        m_inputParameters.m_feederRasterPtr->reset();
        while( ( inputRasterPtr = m_inputParameters.m_feederRasterPtr->getCurrentObj() ) )
        {
          inputRasterIdx = m_inputParameters.m_feederRasterPtr->getCurrentOffset();
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE(
            inputRasterPtr->getAccessPolicy() & te::common::RAccess,
            "Invalid raster" );
          
          // checking the input bands

          for( std::vector< unsigned int >::size_type inputRastersBandsIdx = 0 ;
            inputRastersBandsIdx <
            m_inputParameters.m_inputRastersBands[ inputRasterIdx ].size() ;
            ++inputRastersBandsIdx )
          {
            const unsigned int& currBand =
              m_inputParameters.m_inputRastersBands[ inputRasterIdx ][ inputRastersBandsIdx ];

            TERP_INSTANCE_TRUE_OR_RETURN_FALSE( currBand < inputRasterPtr->getNumberOfBands(),
              "Invalid band" )
          }          
            
          // Defining the base mosaic info

          if( inputRasterIdx == 0 )
          {
            mosaicXResolution = inputRasterPtr->getGrid()->getResolutionX();
            mosaicYResolution = inputRasterPtr->getGrid()->getResolutionY();
            mosaicSRID = inputRasterPtr->getGrid()->getSRID();
            mosaicBaseBandProperties = *inputRasterPtr->getBand( 
              m_inputParameters.m_inputRastersBands[ inputRasterIdx ][ 0 ] )->getProperty();
          }

          rasterExtentPolPtr.reset( static_cast<te::gm::Polygon*>(te::gm::GetGeomFromEnvelope(
            inputRasterPtr->getGrid()->getExtent(), inputRasterPtr->getGrid()->getSRID() ) ) );
          
          if( mosaicSRID != inputRasterPtr->getGrid()->getSRID() )
          {
            rasterExtentPolPtr->transform( mosaicSRID );
          }

          rastersBBoxes.push_back( *rasterExtentPolPtr );

          // expanding mosaic area

          mosaicLLX = std::min( mosaicLLX, rasterExtentPolPtr->getMBR()->getLowerLeftX() );
          mosaicLLY = std::min( mosaicLLY, rasterExtentPolPtr->getMBR()->getLowerLeftY() );

          mosaicURX = std::max( mosaicURX, rasterExtentPolPtr->getMBR()->getUpperRightX() );
          mosaicURY = std::max( mosaicURY, rasterExtentPolPtr->getMBR()->getUpperRightY() );

          // move to the next raster

          m_inputParameters.m_feederRasterPtr->moveNext();
        }
      }
      
      if( m_inputParameters.m_enableProgress )
      {
        progressPtr->pulse();
        if( ! progressPtr->isActive() ) return false;
      }        

      // creating the output raster
      
      te::rst::Raster* outputRasterPtr = nullptr;
      std::vector< double > mosaicBandsRangeMin;
      std::vector< double > mosaicBandsRangeMax;        

      {
        mosaicBandsRangeMin.resize( 
          m_inputParameters.m_inputRastersBands[ 0 ].size(), 0 );
        mosaicBandsRangeMax.resize( 
          m_inputParameters.m_inputRastersBands[ 0 ].size(), 0 );           
        
        std::vector< te::rst::BandProperty* > bandsProperties;
        for( std::vector< unsigned int >::size_type bandIdx = 0 ;  bandIdx <
          m_inputParameters.m_inputRastersBands[ 0 ].size() ; ++bandIdx )
        {
          bandsProperties.push_back( new te::rst::BandProperty( mosaicBaseBandProperties ) );
          bandsProperties[ bandIdx ]->m_colorInterp = te::rst::GrayIdxCInt;
          bandsProperties[ bandIdx ]->m_noDataValue = m_inputParameters.m_noDataValue;
          
          te::rst::GetDataTypeRanges( bandsProperties[ bandIdx ]->m_type,
            mosaicBandsRangeMin[ bandIdx ],
            mosaicBandsRangeMax[ bandIdx ]);           
        }

        te::rst::Grid* outputGrid = new te::rst::Grid( mosaicXResolution,
          mosaicYResolution,  new te::gm::Envelope( mosaicLLX, mosaicLLY, mosaicURX,
          mosaicURY ), mosaicSRID );

        outParamsPtr->m_outputRasterPtr.reset(
          te::rst::RasterFactory::make(
            outParamsPtr->m_rType,
            outputGrid,
            bandsProperties,
            outParamsPtr->m_rInfo,
            nullptr,
            nullptr ) );
        TERP_INSTANCE_TRUE_OR_RETURN_FALSE( outParamsPtr->m_outputRasterPtr.get(),
          "Output raster creation error" );
          
        outputRasterPtr = outParamsPtr->m_outputRasterPtr.get();
      }

      std::unique_ptr< te::mem::CachedRaster > outputCachedRasterInstancePtr;
      
      if( m_inputParameters.m_useRasterCache )
      {
        outputCachedRasterInstancePtr.reset( new te::mem::CachedRaster(
          *(outParamsPtr->m_outputRasterPtr.get()), 20, 0 ) );   
          
        outputRasterPtr = outputCachedRasterInstancePtr.get();
      }

      // fill output with no data values

      {
        const unsigned int nBands = static_cast<unsigned int>(outputRasterPtr->getNumberOfBands());
        const unsigned int nRows = outputRasterPtr->getNumberOfRows();
        const unsigned int nCols = outputRasterPtr->getNumberOfColumns();
        unsigned int col = 0;
        unsigned int row = 0;
        unsigned int bandIdx = 0;

        for( bandIdx = 0 ; bandIdx < nBands ; ++bandIdx )
        {
          te::rst::Band& outBand = ( *( outputRasterPtr->getBand( bandIdx ) ) );

          for( row = 0 ; row < nRows ; ++row )
          {
            for( col = 0 ; col < nCols ; ++col )
            {
              outBand.setValue( col, row, m_inputParameters.m_noDataValue );
            }
          }
        }
      }
      
      if( m_inputParameters.m_enableProgress )
      {
        progressPtr->pulse();
        if( ! progressPtr->isActive() ) return false;
      }        
      
      // Copying the first image data to the output mosaic
      // and find the base mosaic mean and offset values

      std::vector< double > mosaicTargetMeans( outputRasterPtr->getNumberOfBands(), 0 );
      std::vector< double > mosaicTargetVariances( outputRasterPtr->getNumberOfBands(), 0 );

      {
        m_inputParameters.m_feederRasterPtr->reset();

        te::rst::Raster const* inputRasterPtr =
          m_inputParameters.m_feederRasterPtr->getCurrentObj();
        TERP_DEBUG_TRUE_OR_RETURN_FALSE( inputRasterPtr, "Invalid raster pointer" );

        double inXStartGeo = 0;
        double inYStartGeo = 0;
        inputRasterPtr->getGrid()->gridToGeo( 0.0, 0.0, inXStartGeo, inYStartGeo );
        
        double outFirstRowDouble = 0;
        double outFirstColDouble = 0;
        outputRasterPtr->getGrid()->geoToGrid( inXStartGeo, inYStartGeo,
          outFirstColDouble, outFirstRowDouble );
        
        const double outRowsBoundDouble = outFirstRowDouble +
          static_cast<double>( inputRasterPtr->getNumberOfRows() );
        const double outColsBoundDouble = outFirstColDouble +
          static_cast<double>( inputRasterPtr->getNumberOfColumns() );
          
        const unsigned int outFirstRow = static_cast<unsigned int>(std::max( 0u,
          te::common::Round< double, unsigned int >( outFirstRowDouble ) ) );
        const unsigned int outFirstCol = static_cast<unsigned int>(std::max( 0u,
          te::common::Round< double, unsigned int >( outFirstColDouble ) ) );
        
        const unsigned int outRowsBound = static_cast<unsigned int>(std::min(
          te::common::Round< double, unsigned int >( outRowsBoundDouble ),
          outputRasterPtr->getNumberOfRows() ) );
        const unsigned int outColsBound = static_cast<unsigned int>(std::min(
          te::common::Round< double, unsigned int >( outColsBoundDouble ),
          outputRasterPtr->getNumberOfColumns() ) );

        const unsigned int nBands = static_cast<unsigned int>(m_inputParameters.m_inputRastersBands[ 0 ].size());
        unsigned int outCol = 0;
        unsigned int outRow = 0;
        double inCol = 0;
        double inRow = 0;
        double bandNoDataValue = -1.0 * DBL_MAX;
        std::complex< double > pixelCValue = 0;
        te::rst::Interpolator interpInstance( inputRasterPtr,
          m_inputParameters.m_interpMethod );
        unsigned int inputBandIdx = 0;

        for( unsigned int inputRastersBandsIdx = 0 ; inputRastersBandsIdx <
          nBands ; ++inputRastersBandsIdx )
        {
          inputBandIdx =  m_inputParameters.m_inputRastersBands[ 0 ][
            inputRastersBandsIdx ] ;
          bandNoDataValue = m_inputParameters.m_forceInputNoDataValue ?
            m_inputParameters.m_noDataValue : inputRasterPtr->getBand( inputBandIdx
            )->getProperty()->m_noDataValue;
          te::rst::Band& outBand =
            (*outputRasterPtr->getBand( inputRastersBandsIdx ));
          unsigned long int validPixelsNumber = 0;

          double mean = 0;

          for( outRow = outFirstRow ; outRow < outRowsBound ; ++outRow )
          {
            inRow = (static_cast<double>(outRow)) - outFirstRowDouble;

            for( outCol = outFirstCol ; outCol < outColsBound ; ++outCol )
            {
              inCol = (static_cast<double>(outCol)) - outFirstColDouble;

              interpInstance.getValue( inCol, inRow, pixelCValue, inputBandIdx );

              if( pixelCValue.real() != bandNoDataValue )
              {
                outBand.setValue( outCol, outRow, pixelCValue );
                mean += pixelCValue.real();
                ++validPixelsNumber;
              }
            }
          }

          // variance calcule

          if( m_inputParameters.m_autoEqualize && ( validPixelsNumber > 0 ) )
          {
            mean /= ( static_cast<double>(validPixelsNumber) );
            mosaicTargetMeans[ inputRastersBandsIdx ] = mean;
            
            double& variance = mosaicTargetVariances[ inputRastersBandsIdx ];
            variance = 0;

            double pixelValue = 0;

            for( outRow = outFirstRow ; outRow < outRowsBound ; ++outRow )
            {
              for( outCol = outFirstCol ; outCol < outColsBound ; ++outCol )
              {
                outBand.getValue( outCol, outRow, pixelValue );

                if( pixelValue != m_inputParameters.m_noDataValue )
                {
                  variance += ( ( pixelValue - mean ) * ( pixelValue -
                    mean ) ) / ( static_cast<double>(validPixelsNumber) );
                }
              }
            }
          }
        }
      }
      
      TERP_DEBUG_TRUE_OR_THROW( rastersBBoxes.size() ==
        m_inputParameters.m_feederRasterPtr->getObjsCount(),
        "Rasters bounding boxes number mismatch" );
        
      if( m_inputParameters.m_enableProgress )
      {
        progressPtr->pulse();
        if( ! progressPtr->isActive() ) return false;
      }          

      // Initiating the mosaic bounding boxes union

      std::unique_ptr< te::gm::MultiPolygon > mosaicBBoxesUnionPtr(
        new te::gm::MultiPolygon( 0, te::gm::MultiPolygonType,
        outputRasterPtr->getSRID(), nullptr ) );
      mosaicBBoxesUnionPtr->add( static_cast<te::gm::Polygon*>(rastersBBoxes[ 0 ].clone()) );

      // skipping the first raster

      m_inputParameters.m_feederRasterPtr->reset();
      m_inputParameters.m_feederRasterPtr->moveNext();

      // iterating over the other rasters
      
      te::rst::Raster const* originalInputRasterPtr = nullptr;

      std::vector< unsigned int > outputRasterBands;
      std::vector< double > dummyRasterOffsets;
      std::vector< double > dummyRasterScales;
      for( unsigned int outputRasterBandsIdx = 0 ; outputRasterBandsIdx <
        outputRasterPtr->getNumberOfBands() ; ++outputRasterBandsIdx )
      {
        outputRasterBands.push_back( outputRasterBandsIdx );
        dummyRasterOffsets.push_back( 0.0 );
        dummyRasterScales.push_back( 1.0 );
      }      

      while( ( originalInputRasterPtr = m_inputParameters.m_feederRasterPtr->getCurrentObj() ) )
      {
        const unsigned int inputRasterIdx = m_inputParameters.m_feederRasterPtr->getCurrentOffset();
        
//         Copy2DiskRaster( outputRaster, boost::lexical_cast< std::string >( inputRasterIdx ) +
//           "_output_raster_begininng.tif" );

        TERP_DEBUG_TRUE_OR_THROW( rastersBBoxes[ inputRasterIdx ].getSRID() == outputRasterPtr->getSRID(),
          "Invalid boxes SRID" );
          
        te::rst::Raster const* inputRasterPtr = originalInputRasterPtr;
        
//         Copy2DiskRaster( *inputRasterPtr, boost::lexical_cast< std::string >( inputRasterIdx ) +
//           "_input_raster.tif" );

        // reprojection / caching issues
        
        std::unique_ptr< te::rst::Raster > reprojectedInputRasterPtr;
        
        if( outputRasterPtr->getSRID() != inputRasterPtr->getSRID() )
        {
          std::map< std::string, std::string > rInfo;
          rInfo[ "RTYPE" ] = "EXPANSIBLE";
          
          reprojectedInputRasterPtr.reset( inputRasterPtr->transform( 
            outputRasterPtr->getSRID(), rInfo,
            m_inputParameters.m_interpMethod) );
          inputRasterPtr = reprojectedInputRasterPtr.get();
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( inputRasterPtr, "Reprojection error" );
          
//           Copy2DiskRaster( *inputRasterPtr, boost::lexical_cast< std::string >( inputRasterIdx ) +
//             "reprojected_input_raster.tif" );
        }
        
        // Caching issues 
        
        std::unique_ptr< te::mem::CachedRaster > cachedInputRasterPtr;
        
        if(
            m_inputParameters.m_useRasterCache
           &&
           ( reprojectedInputRasterPtr.get() == nullptr )
           // Reprojected rasters already have internal cache
          )
        {
          cachedInputRasterPtr.reset( new te::mem::CachedRaster(
            *inputRasterPtr, 20, 0 ) );
          inputRasterPtr = cachedInputRasterPtr.get();
        }         
        
        // The transformation mapping outputRaster pixels 
        // ( te::gm::GTParameters::TiePoint::first ) to input pixels 
        // ( te::gm::GTParameters::TiePoint::second ) 
        // (Note: all coords are indexed by lines/columns).
        
        std::unique_ptr< te::gm::GeometricTransformation > transPtr;
        
        {
          te::gm::GTParameters transParams;
          double auxX;
          double auxY;
          te::gm::GTParameters::TiePoint auxTP;
          
          auxTP.first.x = 0;
          auxTP.first.y = 0;
          outputRasterPtr->getGrid()->gridToGeo( auxTP.first.x, auxTP.first.y, auxX, auxY );
          inputRasterPtr->getGrid()->geoToGrid( auxX, auxY, auxTP.second.x, auxTP.second.y );
          transParams.m_tiePoints.push_back( auxTP );
          
          auxTP.first.x = 0;
          auxTP.first.y = static_cast<double>( outputRasterPtr->getNumberOfRows() - 1 );
          outputRasterPtr->getGrid()->gridToGeo( auxTP.first.x, auxTP.first.y, auxX, auxY );
          inputRasterPtr->getGrid()->geoToGrid( auxX, auxY, auxTP.second.x, auxTP.second.y );          
          transParams.m_tiePoints.push_back( auxTP );

          auxTP.first.x = static_cast<double>( outputRasterPtr->getNumberOfColumns() - 1);
          auxTP.first.y = 0;
          outputRasterPtr->getGrid()->gridToGeo( auxTP.first.x, auxTP.first.y, auxX, auxY );
          inputRasterPtr->getGrid()->geoToGrid( auxX, auxY, auxTP.second.x, auxTP.second.y );
          transParams.m_tiePoints.push_back( auxTP );
          
          auxTP.first.x = static_cast<double>( outputRasterPtr->getNumberOfColumns() - 1 );
          auxTP.first.y = static_cast<double>( outputRasterPtr->getNumberOfRows() - 1 );
          outputRasterPtr->getGrid()->gridToGeo( auxTP.first.x, auxTP.first.y, auxX, auxY );
          inputRasterPtr->getGrid()->geoToGrid( auxX, auxY, auxTP.second.x, auxTP.second.y );                    
          transParams.m_tiePoints.push_back( auxTP );
          
          transPtr.reset( te::gm::GTFactory::make( "Affine" ) );
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( transPtr.get(), "Could not instantiate a geometric transformation" );
          
          TERP_INSTANCE_TRUE_OR_RETURN_FALSE( transPtr->initialize( transParams ),
            "Could not initialize a geometric transformation" );
        }

        // Generating the offset and gain info for eath band from the current raster

        std::vector< double > currentRasterBandsOffsets = dummyRasterOffsets;
        std::vector< double > currentRasterBandsScales = dummyRasterScales;

        if( m_inputParameters.m_autoEqualize )
        {
          double currentRasterVariance = 0;
          double currentRasterMean = 0;

          for( unsigned int inputRastersBandsIdx = 0 ; inputRastersBandsIdx <
            m_inputParameters.m_inputRastersBands[ inputRasterIdx ].size() ;
            ++inputRastersBandsIdx )
          {
            const unsigned int inputBandIdx = m_inputParameters.m_inputRastersBands[ inputRasterIdx ][
              inputRastersBandsIdx ];

            if( 
                ( mosaicTargetMeans[ inputRastersBandsIdx ] != 0.0 )
                &&
                ( mosaicTargetVariances[ inputRastersBandsIdx ] != 0.0 )
              )
            {
              calcBandStatistics( (*inputRasterPtr->getBand( inputBandIdx ) ),
                m_inputParameters.m_forceInputNoDataValue,
                m_inputParameters.m_noDataValue,
                currentRasterMean,
                currentRasterVariance );

              currentRasterBandsScales[ inputRastersBandsIdx ] = 
                ( 
                  std::sqrt( mosaicTargetVariances[ inputRastersBandsIdx ] )
                   /
                  std::sqrt( currentRasterVariance ) 
                );
              currentRasterBandsOffsets[ inputRastersBandsIdx ] =
                ( 
                  mosaicTargetMeans[ inputRastersBandsIdx ]
                   - 
                  ( 
                    currentRasterBandsScales[ inputRastersBandsIdx ] 
                    * 
                    currentRasterMean 
                  ) 
                );
            }
          }
        }

         // blending

        te::rp::Blender blenderInstance;

        TERP_INSTANCE_TRUE_OR_RETURN_FALSE( blenderInstance.initialize(
          *outputRasterPtr,
          outputRasterBands,
          *inputRasterPtr,
          m_inputParameters.m_inputRastersBands[ inputRasterIdx ],
          m_inputParameters.m_blendMethod,
          te::rst::NearestNeighbor,
          m_inputParameters.m_interpMethod,
          m_inputParameters.m_noDataValue,
          false,
          m_inputParameters.m_forceInputNoDataValue,
          dummyRasterOffsets,
          dummyRasterScales,
          currentRasterBandsOffsets,
          currentRasterBandsScales,
          mosaicBBoxesUnionPtr.get(),
          nullptr,
          *transPtr,
          m_inputParameters.m_enableMultiThread ? 0 : 1,
          false ), "Blender initiazing error" );
        
        TERP_INSTANCE_TRUE_OR_RETURN_FALSE( blenderInstance.blendIntoRaster1(),
          "Error blending images" );

//           Copy2DiskRaster( outputRaster, boost::lexical_cast< std::string >( inputRasterIdx ) +
//              "output_raster_after_blending_" + 
//             boost::lexical_cast< std::string >( mosaicBBoxesUnionIdx ) + ".tif" );

        // updating the  gloabal mosaic boxes

        std::unique_ptr< te::gm::Geometry > boxesUnionResultPtr; // under the mosaic SRID
        TERP_INSTANCE_TRUE_OR_RETURN_FALSE( mosaicBBoxesUnionPtr->isValid(),
          "Invalid mosaic bounding boxes union geometry" );
        TERP_INSTANCE_TRUE_OR_RETURN_FALSE( rastersBBoxes[ inputRasterIdx ].isValid(),
          "Invalid raster bounding boxes union geometry (raster index:"
          + boost::lexical_cast< std::string >( inputRasterIdx ) + ")" );  
        
        try
        {
          boxesUnionResultPtr.reset( mosaicBBoxesUnionPtr->Union(
            &( rastersBBoxes[ inputRasterIdx ] ) ) );
        }
        catch( const std::exception&/* e */)
        {
          TERP_INSTANCE_LOG_AND_RETURN_FALSE( "Mosaic bounding boxes union error" );
        }
        
        TERP_TRUE_OR_THROW( boxesUnionResultPtr.get(), "Invalid pointer" );
        
        boxesUnionResultPtr->setSRID( outputRasterPtr->getSRID() );

        if( 
            ( boxesUnionResultPtr->getGeomTypeId() == te::gm::MultiPolygonType )
            ||
            ( boxesUnionResultPtr->getGeomTypeId() == te::gm::MultiPolygonZType )
            ||
            ( boxesUnionResultPtr->getGeomTypeId() == te::gm::MultiPolygonMType )
            ||
            ( boxesUnionResultPtr->getGeomTypeId() == te::gm::MultiPolygonZMType )
          )
        {
          mosaicBBoxesUnionPtr.reset( static_cast<te::gm::MultiPolygon*>(boxesUnionResultPtr.release()) );
        }
        else if( 
                 ( boxesUnionResultPtr->getGeomTypeId() == te::gm::PolygonType )
                 ||
                 ( boxesUnionResultPtr->getGeomTypeId() == te::gm::PolygonZType )
                 ||
                 ( boxesUnionResultPtr->getGeomTypeId() == te::gm::PolygonMType )
                 ||
                 ( boxesUnionResultPtr->getGeomTypeId() == te::gm::PolygonZMType )
               )
        {
          // transforming it into a te::gm::MultiPolygon
          te::gm::MultiPolygon* auxMultiPol = new te::gm::MultiPolygon( 0,
            te::gm::MultiPolygonType, boxesUnionResultPtr->getSRID(), nullptr );
          auxMultiPol->add( boxesUnionResultPtr.release() );

          mosaicBBoxesUnionPtr.reset( auxMultiPol );
        }
        else
        {
          TERP_LOGWARN( "Invalid union geometry type" );
        }
        
        // Reseting the input cache
        
        cachedInputRasterPtr.reset();
        
        // moving to the next raster

        m_inputParameters.m_feederRasterPtr->moveNext();
        
        if( m_inputParameters.m_enableProgress )
        {
          progressPtr->pulse();
          if( ! progressPtr->isActive() ) return false;
        }          
      }
      
      // reseting the output cache
      
      if( outputCachedRasterInstancePtr.get() ) outputCachedRasterInstancePtr.reset();
      
      return true;
    }

    void GeoMosaic::reset() _NOEXCEPT_OP(false)
    {
      Algorithm::reset();
      
      m_inputParameters.reset();
      m_isInitialized = false;
    }

    bool GeoMosaic::initialize( const AlgorithmInputParameters& inputParams )
      _NOEXCEPT_OP(false)
    {
      reset();

      GeoMosaic::InputParameters const* inputParamsPtr = dynamic_cast<
        GeoMosaic::InputParameters const* >( &inputParams );
      TERP_TRUE_OR_THROW( inputParamsPtr, "Invalid paramters pointer" );

      m_inputParameters = *inputParamsPtr;

      // Checking the feeder

      TERP_INSTANCE_TRUE_OR_RETURN_FALSE( m_inputParameters.m_feederRasterPtr,
        "Invalid m_feederRasterPtr" )

      TERP_INSTANCE_TRUE_OR_RETURN_FALSE(
        m_inputParameters.m_feederRasterPtr->getObjsCount() > 0,
        "Invalid number of rasters" )

      // checking m_inputRastersBands

      TERP_INSTANCE_TRUE_OR_RETURN_FALSE(
        (static_cast<unsigned int>(m_inputParameters.m_inputRastersBands.size())) ==
        m_inputParameters.m_feederRasterPtr->getObjsCount(),
        "Bands mismatch" );

      for( std::vector< std::vector< unsigned int > >::size_type
        inputRastersBandsIdx = 0 ;  inputRastersBandsIdx <
        m_inputParameters.m_inputRastersBands.size() ; ++inputRastersBandsIdx )
      {
        TERP_INSTANCE_TRUE_OR_RETURN_FALSE( m_inputParameters.m_inputRastersBands[
          inputRastersBandsIdx ].size() > 0, "Invalid bands number" );

        TERP_INSTANCE_TRUE_OR_RETURN_FALSE( m_inputParameters.m_inputRastersBands[
          inputRastersBandsIdx ].size() ==  m_inputParameters.m_inputRastersBands[
          0 ].size(), "Bands number mismatch" );
      }

      m_isInitialized = true;

      return true;
    }

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

    void GeoMosaic::calcBandStatistics( const te::rst::Band& band,
      const bool& forceNoDataValue,
      const double& noDataValue,
      double& mean, double& variance )
    {
      mean = 0;
      variance = 0;

      double internalNoDataValue = 0;
      if( forceNoDataValue )
        internalNoDataValue = noDataValue;
      else
        internalNoDataValue = band.getProperty()->m_noDataValue;

      const unsigned int nCols = static_cast<unsigned int>(band.getProperty()->m_blkw * band.getProperty()->m_nblocksx);
      const unsigned int nLines = static_cast<unsigned int>(band.getProperty()->m_blkh * band.getProperty()->m_nblocksy);

      double pixelsNumber = 0;
      double value = 0;
      unsigned int col = 0;
      unsigned int line = 0;

      for( line = 0 ; line < nLines ; ++line )
        for( col = 0 ; col < nCols ; ++col )
        {
          band.getValue( col, line, value );

          if( value != internalNoDataValue )
          {
            mean += value;
            ++pixelsNumber;
          }
        }

      if( pixelsNumber != 0.0 )
      {
        mean /= pixelsNumber;

        for( line = 0 ; line < nLines ; ++line )
          for( col = 0 ; col < nCols ; ++col )
          {
            band.getValue( col, line, value );

            if( value != internalNoDataValue )
            {
              variance += ( ( value - mean ) * ( value - mean ) ) / pixelsNumber;
            }
          }

      }
    }

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