29 #include "../raster/Interpolator.h"
30 #include "../raster/Enums.h"
31 #include "../raster/RasterFactory.h"
32 #include "../raster/Grid.h"
33 #include "../raster/Band.h"
34 #include "../raster/BandProperty.h"
35 #include "../raster/PositionIterator.h"
36 #include "../raster/Utils.h"
37 #include "../memory/CachedRaster.h"
38 #include "../geometry/Envelope.h"
39 #include "../geometry/GTFactory.h"
40 #include "../geometry/Polygon.h"
41 #include "../geometry/LinearRing.h"
42 #include "../geometry/MultiPolygon.h"
43 #include "../srs/Converter.h"
44 #include "../common/progress/TaskProgress.h"
46 #include <boost/shared_array.hpp>
47 #include <boost/filesystem.hpp>
78 m_feederRasterPtr = 0;
79 m_inputRastersBands.clear();
82 m_forceInputNoDataValue =
false;
84 m_autoEqualize =
true;
85 m_useRasterCache =
true;
86 m_enableProgress =
false;
87 m_enableMultiThread =
true;
134 m_outputRasterPtr.reset();
144 m_outputRasterPtr.reset();;
164 throw( te::rp::Exception )
176 std::auto_ptr< te::common::TaskProgress > progressPtr;
183 progressPtr->setMessage(
"Mosaic" );
188 double mosaicXResolution = 0.0;
189 double mosaicYResolution = 0.0;
190 double mosaicLLX = DBL_MAX;
191 double mosaicLLY = DBL_MAX;
192 double mosaicURX = -1.0 * DBL_MAX;
193 double mosaicURY = -1.0 * DBL_MAX;
195 std::vector< double > mosaicBandsRangeMin;
196 std::vector< double > mosaicBandsRangeMax;
198 std::vector< te::gm::Polygon > rastersBBoxes;
204 unsigned int inputRasterIdx = 0;
218 if( inputRasterIdx == 0 )
235 auxLinearRingPtr->
setPoint( 0, mosaicLLX, mosaicURY );
236 auxLinearRingPtr->
setPoint( 1, mosaicURX, mosaicURY );
237 auxLinearRingPtr->
setPoint( 2, mosaicURX, mosaicLLY );
238 auxLinearRingPtr->
setPoint( 3, mosaicLLX, mosaicLLY );
239 auxLinearRingPtr->
setPoint( 4, mosaicLLX, mosaicURY );
240 auxLinearRingPtr->
setSRID( mosaicSRID );
241 auxPolygon.
push_back( auxLinearRingPtr );
242 auxPolygon.
setSRID( mosaicSRID );
243 rastersBBoxes.push_back( auxPolygon );
248 *inputRasterPtr->
getGrid(), reprojectedExtent ),
249 "Detailed raster extent calcule error" );
253 reprojectedExtent.
transform( mosaicSRID );
268 auxPolygon.
push_back( auxLinearRingPtr );
269 auxPolygon.
setSRID( mosaicSRID );
270 rastersBBoxes.push_back( auxPolygon );
275 for( std::vector< unsigned int >::size_type inputRastersBandsIdx = 0 ;
276 inputRastersBandsIdx <
278 ++inputRastersBandsIdx )
280 const unsigned int& currBand =
294 progressPtr->pulse();
295 if( ! progressPtr->isActive() )
return false;
303 mosaicBandsRangeMin.resize(
305 mosaicBandsRangeMax.resize(
308 std::vector< te::rst::BandProperty* > bandsProperties;
309 for( std::vector< unsigned int >::size_type bandIdx = 0 ; bandIdx <
317 mosaicBandsRangeMin[ bandIdx ],
318 mosaicBandsRangeMax[ bandIdx ]);
323 mosaicURY ), mosaicSRID );
334 "Output raster creation error" );
339 std::auto_ptr< te::mem::CachedRaster > cachedRasterInstancePtr;
346 outputRasterPtr = cachedRasterInstancePtr.get();
355 unsigned int col = 0;
356 unsigned int row = 0;
357 unsigned int bandIdx = 0;
359 for( bandIdx = 0 ; bandIdx < nBands ; ++bandIdx )
363 for( row = 0 ; row < nRows ; ++row )
365 for( col = 0 ; col < nCols ; ++col )
375 progressPtr->pulse();
376 if( ! progressPtr->isActive() )
return false;
382 std::vector< double > mosaicTargetMeans( outputRasterPtr->
getNumberOfBands(), 0 );
383 std::vector< double > mosaicTargetVariances( outputRasterPtr->
getNumberOfBands(), 0 );
392 double inXStartGeo = 0;
393 double inYStartGeo = 0;
395 double outRowStartDouble = 0;
396 double outColStartDouble = 0;
398 outColStartDouble, outRowStartDouble );
400 const unsigned int outRowStart = (
unsigned int)std::max( 0.0, outRowStartDouble );
401 const unsigned int outColStart = (
unsigned int)std::max( 0.0, outColStartDouble );
402 const unsigned int outRowsBound = std::min( outRowStart +
405 const unsigned int outColsBound = std::min( outColStart +
409 const unsigned int nBands = (
unsigned int)
411 unsigned int outCol = 0;
412 unsigned int outRow = 0;
415 double bandNoDataValue = -1.0 * DBL_MAX;
416 std::complex< double > pixelCValue = 0;
419 unsigned int inputBandIdx = 0;
421 for(
unsigned int inputRastersBandsIdx = 0 ; inputRastersBandsIdx <
422 nBands ; ++inputRastersBandsIdx )
425 inputRastersBandsIdx ] ;
430 (*outputRasterPtr->
getBand( inputRastersBandsIdx ));
431 unsigned int validPixelsNumber = 0;
433 double& mean = mosaicTargetMeans[ inputRastersBandsIdx ];
436 for( outRow = outRowStart ; outRow < outRowsBound ; ++outRow )
438 inRow = ((double)outRow) - outRowStartDouble;
440 for( outCol = outColStart ; outCol < outColsBound ; ++outCol )
442 inCol = ((double)outCol) - outColStartDouble;
444 interpInstance.
getValue( inCol, inRow, pixelCValue, inputBandIdx );
446 if( pixelCValue.real() != bandNoDataValue )
448 outBand.
setValue( outCol, outRow, pixelCValue );
449 mean += pixelCValue.real();
455 mean /= ( (double)validPixelsNumber );
461 double& variance = mosaicTargetVariances[ inputRastersBandsIdx ];
464 double pixelValue = 0;
466 for( outRow = outRowStart ; outRow < outRowsBound ; ++outRow )
468 for( outCol = outColStart ; outCol < outColsBound ; ++outCol )
470 outBand.
getValue( outCol, outRow, pixelValue );
474 variance += ( ( pixelValue - mean ) * ( pixelValue -
475 mean ) ) / ( (
double)validPixelsNumber );
485 "Rasters bounding boxes number mismatch" );
489 progressPtr->pulse();
490 if( ! progressPtr->isActive() )
return false;
495 std::auto_ptr< te::gm::MultiPolygon > mosaicBBoxesUnionPtr(
497 outputRasterPtr->
getSRID(), 0 ) );
498 mosaicBBoxesUnionPtr->add( (
te::gm::Polygon*)rastersBBoxes[ 0 ].clone() );
509 std::vector< unsigned int > outputRasterBands;
510 std::vector< double > dummyRasterOffsets;
511 std::vector< double > dummyRasterScales;
512 for(
unsigned int outputRasterBandsIdx = 0 ; outputRasterBandsIdx <
515 outputRasterBands.push_back( outputRasterBandsIdx );
516 dummyRasterOffsets.push_back( 0.0 );
517 dummyRasterScales.push_back( 1.0 );
528 "Invalid boxes SRID" );
537 std::auto_ptr< te::rst::Raster > reprojectedInputRasterPtr;
539 std::string reprojectedRasterFileName;
541 if( outputRasterPtr->
getSRID() != inputRasterPtr->getSRID() )
543 reprojectedRasterFileName = boost::filesystem::unique_path(
544 boost::filesystem::temp_directory_path() /=
545 boost::filesystem::path(
"TerralibRGeoMosaic_%%%%-%%%%-%%%%-%%%%" ) ).string();
547 "Invalid temporary raster file name" );
548 reprojectedRasterFileName +=
".tif";
550 std::map< std::string, std::string > rInfo;
551 rInfo[
"URI" ] = reprojectedRasterFileName;
553 reprojectedInputRasterPtr.reset( inputRasterPtr->transform(
554 outputRasterPtr->
getSRID(), rInfo,
556 inputRasterPtr = reprojectedInputRasterPtr.get();
565 std::auto_ptr< te::mem::CachedRaster > cachedInputRasterPtr;
570 *inputRasterPtr, 20, 0 ) );
571 inputRasterPtr = cachedInputRasterPtr.get();
579 std::auto_ptr< te::gm::GeometricTransformation > transPtr;
589 outputRasterPtr->
getGrid()->
gridToGeo( auxTP.first.x, auxTP.first.y, auxX, auxY );
590 inputRasterPtr->getGrid()->geoToGrid( auxX, auxY, auxTP.second.x, auxTP.second.y );
595 outputRasterPtr->
getGrid()->
gridToGeo( auxTP.first.x, auxTP.first.y, auxX, auxY );
596 inputRasterPtr->getGrid()->geoToGrid( auxX, auxY, auxTP.second.x, auxTP.second.y );
607 outputRasterPtr->
getGrid()->
gridToGeo( auxTP.first.x, auxTP.first.y, auxX, auxY );
608 inputRasterPtr->getGrid()->geoToGrid( auxX, auxY, auxTP.second.x, auxTP.second.y );
615 "Could not initialize a geometric transformation" );
620 std::vector< double > currentRasterBandsOffsets;
621 std::vector< double > currentRasterBandsScales;
625 double currentRasterVariance = 0;
626 double currentRasterMean = 0;
628 for(
unsigned int inputRastersBandsIdx = 0 ; inputRastersBandsIdx <
630 ++inputRastersBandsIdx )
633 inputRastersBandsIdx ];
639 currentRasterVariance );
641 currentRasterBandsScales.push_back(
642 std::sqrt( mosaicTargetVariances[ inputRastersBandsIdx ] )
644 std::sqrt( currentRasterVariance ) );
645 currentRasterBandsOffsets.push_back( mosaicTargetMeans[ inputRastersBandsIdx ]
646 - ( currentRasterBandsScales[ inputRastersBandsIdx ] * currentRasterMean ) );
651 currentRasterBandsOffsets = dummyRasterOffsets;
652 currentRasterBandsScales = dummyRasterScales;
672 currentRasterBandsOffsets,
673 currentRasterBandsScales,
674 mosaicBBoxesUnionPtr.get(),
678 false ),
"Blender initiazing error" );
681 "Error blending images" );
689 std::auto_ptr< te::gm::Geometry > boxesUnionResultPtr;
690 boxesUnionResultPtr.reset( mosaicBBoxesUnionPtr->Union(
691 &( rastersBBoxes[ inputRasterIdx ] ) ) );
693 boxesUnionResultPtr->setSRID( outputRasterPtr->
getSRID() );
704 auxMultiPol->
add( boxesUnionResultPtr.release() );
706 mosaicBBoxesUnionPtr.reset( auxMultiPol );
715 cachedInputRasterPtr.reset();
719 if( reprojectedInputRasterPtr.get() )
721 reprojectedInputRasterPtr.reset();
722 remove( reprojectedRasterFileName.c_str() );
731 progressPtr->pulse();
732 if( ! progressPtr->isActive() )
return false;
738 if( cachedRasterInstancePtr.get() ) cachedRasterInstancePtr.reset();
750 throw( te::rp::Exception )
763 "Invalid m_feederRasterPtr" )
767 "Invalid number of rasters" )
776 for( std::vector< std::vector< unsigned int > >::size_type
777 inputRastersBandsIdx = 0 ; inputRastersBandsIdx <
781 inputRastersBandsIdx ].size() > 0,
"Invalid bands number" );
785 0 ].size(),
"Bands number mismatch" );
799 const bool& forceNoDataValue,
800 const double& noDataValue,
801 double& mean,
double& variance )
806 double internalNoDataValue = 0;
807 if( forceNoDataValue )
808 internalNoDataValue = noDataValue;
817 double pixelsNumber = 0;
819 unsigned int col = 0;
820 unsigned int line = 0;
822 for( line = 0 ; line < nLines ; ++line )
823 for( col = 0 ; col < nCols ; ++col )
827 if( value != internalNoDataValue )
834 if( pixelsNumber != 0.0 )
836 mean /= pixelsNumber;
838 for( line = 0 ; line < nLines ; ++line )
839 for( col = 0 ; col < nCols ; ++col )
843 if( value != internalNoDataValue )
845 variance += ( ( value - mean ) * ( value - mean ) ) / pixelsNumber;
void reset()
Clear all internal allocated resources and reset the parameters instance to its initial state...
virtual unsigned int getObjsCount() const =0
Return the total number of feeder objects.
void push_back(Curve *ring)
It adds the curve to the curve polygon.
MultiPolygon is a MultiSurface whose elements are Polygons.
bool initialize(const AlgorithmInputParameters &inputParams)
Initialize the algorithm instance making it ready for execution.
Near neighborhood interpolation method.
std::vector< TiePoint > m_tiePoints
Tie points.
void reset()
Clear all internal allocated objects and reset the algorithm to its initial state.
const OutputParameters & operator=(const OutputParameters ¶ms)
TERASTEREXPORT void GetDataTypeRanges(const int &dataType, double &min, double &max)
Return the values range of a given data type.
Index into a lookup table.
virtual void getValue(unsigned int c, unsigned int r, double &value) const =0
Returns the cell attribute value.
void setSRID(int srid)
It sets the Spatial Reference System ID of the geometry and all its parts if it is a GeometryCollecti...
AbstractParameters * clone() const
Create a clone copy of this instance.
bool isInitialized() const
Returns true if the algorithm instance is initialized and ready for execution.
A raster band description.
bool execute(AlgorithmOutputParameters &outputParams)
Executes the algorithm using the supplied parameters.
Blended pixel value calculation for two overlaped rasters.
int getSRID() const
Returns the grid spatial reference system identifier.
unsigned int getNumberOfColumns() const
Returns the raster number of columns.
int m_nblocksx
The number of blocks in x.
virtual const Band * getBand(std::size_t i) const =0
Returns the raster i-th band.
It interpolates one pixel based on a selected algorithm. Methods currently available are Nearest Neig...
int m_nblocksy
The number of blocks in y.
const double & getUpperRightX() const
It returns a constant refernce to the x coordinate of the upper right corner.
This class can be used to inform the progress of a task.
Raster Processing algorithm output parameters base interface.
double m_urx
Upper right corner x-coordinate.
const double & getLowerLeftY() const
It returns a constant refernce to the y coordinate of the lower left corner.
void getValue(const double &c, const double &r, std::complex< double > &v, const std::size_t &b)
Get the interpolated value at specific band.
double m_noDataValue
Value to indicate elements where there is no data, default is std::numeric_limits::max().
const double & getUpperRightY() const
It returns a constant refernce to the x coordinate of the upper right corner.
std::pair< Coord2D, Coord2D > TiePoint
Tie point type definition.
std::string m_rType
Output raster data source type (as described in te::raster::RasterFactory ).
void geoToGrid(const double &x, const double &y, double &col, double &row) const
Get the grid point associated to a spatial location.
bool blendIntoRaster1()
Execute blending of the given input rasters and write the result into raster1.
Raster Processing functions.
double getResolutionY() const
Returns the grid vertical (y-axis) resolution.
#define TERP_TRUE_OR_RETURN_FALSE(value, message)
Checks if value is true. For false values a warning message will be logged and a return of context wi...
const Algorithm & operator=(const Algorithm &)
te::common::AccessPolicy getAccessPolicy() const
Returns the raster access policy.
A LinearRing is a LineString that is both closed and simple.
double m_llx
Lower left corner x-coordinate.
void setPoint(std::size_t i, const double &x, const double &y)
It sets the value of the specified point.
virtual void reset()=0
Reset the feeder to the first position (subsequent accesses will start from the first sequence obejct...
An Envelope defines a 2D rectangular region.
An abstract class for raster data strucutures.
unsigned int getNumberOfRows() const
Returns the raster number of rows.
BandProperty * getProperty()
Returns the band property.
int m_blkw
Block width (pixels).
virtual std::size_t getNumberOfBands() const =0
Returns the number of bands (dimension of cells attribute values) in the raster.
Create a mosaic from a set of geo-referenced rasters.
double getResolutionX() const
Returns the grid horizontal (x-axis) resolution.
A RAM cache adaptor to an external existent raster that must always be avaliable. ...
A raster band description.
bool GetDetailedExtent(const te::rst::Grid &grid, te::gm::LinearRing &detailedExtent)
Create a datailed extent from the given grid.
Grid * getGrid()
It returns the raster grid.
static GeometricTransformation * make(const std::string &factoryKey)
It creates an object with the appropriated factory.
bool initialize(te::rst::Raster &raster1, const std::vector< unsigned int > &raster1Bands, const te::rst::Raster &raster2, const std::vector< unsigned int > &raster2Bands, const BlendMethod &blendMethod, const te::rst::Interpolator::Method &interpMethod1, const te::rst::Interpolator::Method &interpMethod2, const double &noDataValue, const bool forceRaster1NoDataValue, const bool forceRaster2NoDataValue, const std::vector< double > &pixelOffsets1, const std::vector< double > &pixelScales1, const std::vector< double > &pixelOffsets2, const std::vector< double > &pixelScales2, te::gm::MultiPolygon const *const r1ValidDataDelimiterPtr, te::gm::MultiPolygon const *const r2ValidDataDelimiterPtr, const te::gm::GeometricTransformation &geomTransformation, const unsigned int threadsNumber, const bool enableProgressInterface)
Inititate the blender instance.
double m_lly
Lower left corner y-coordinate.
virtual void setValue(unsigned int c, unsigned int r, const double value)=0
Sets the cell attribute value.
A Converter is responsible for the conversion of coordinates between different Coordinate Systems (CS...
Abstract parameters base interface.
std::auto_ptr< te::rst::Raster > m_outputRasterPtr
The generated output mosaic raster.
bool m_isInitialized
Tells if this instance is initialized.
int getSRID() const
Returns the raster spatial reference system identifier.
#define TERP_LOG_AND_RETURN_FALSE(message)
Logs a warning message will and return false.
Polygon is a subclass of CurvePolygon whose rings are defined by linear rings.
double m_ury
Upper right corner y-coordinate.
static void calcBandStatistics(const te::rst::Band &band, const bool &forceNoDataValue, const double &noDataValue, double &mean, double &variance)
Raster band statistics calcule.
GeoMosaic output parameters.
void add(Geometry *g)
It adds the geometry into the collection.
GeoMosaic::InputParameters m_inputParameters
Input execution parameters.
static Raster * make()
It creates and returns an empty raster with default raster driver.
const double & getLowerLeftX() const
It returns a constant reference to the x coordinate of the lower left corner.
#define TERP_DEBUG_TRUE_OR_RETURN_FALSE(value, message)
Checks if value is true. For false values a warning message will be logged and a return of context wi...
te::gm::Envelope * getExtent()
Returns the geographic extension of the grid.
void clear()
It deletes all the rings of the CurvePolygon and clear it.
void gridToGeo(const double &col, const double &row, double &x, double &y) const
Get the spatial location of a grid point.
void setSRID(int srid)
It sets the Spatial Reference System ID of the linestring.
2D Geometric transformation parameters.
virtual bool moveNext()=0
Advances to the next sequence obeject.
int m_blkh
Block height (pixels).
#define TERP_DEBUG_TRUE_OR_THROW(value, message)
Checks if value is true and throws an exception if not.
virtual te::rst::Raster const * getCurrentObj() const =0
Return the current sequence object.
A rectified grid is the spatial support for raster data.
const Envelope * getMBR() const
It returns the minimum bounding rectangle for the geometry in an internal representation.
virtual unsigned int getCurrentOffset() const =0
Return the index of the current object.
#define TERP_TRUE_OR_THROW(value, message)
Checks if value is true and throws an exception if not.
void transform(int srid)
It converts the coordinate values of the linestring to the new spatial reference system.
std::map< std::string, std::string > m_rInfo
The necessary information to create the output rasters (as described in te::raster::RasterFactory).