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_ptr.hpp>
47 #include <boost/scoped_array.hpp>
77 m_feederRasterPtr = 0;
78 m_inputRastersBands.clear();
80 m_geomTransfName =
"Affine";
83 m_forceInputNoDataValue =
false;
85 m_autoEqualize =
true;
86 m_useRasterCache =
true;
87 m_enableProgress =
false;
88 m_enableMultiThread =
true;
137 m_outputRasterPtr.reset();
166 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;
201 std::vector< boost::shared_ptr< te::gm::GeometricTransformation > >
202 eachRasterPixelToFirstRasterPixelGeomTransfms;
207 std::vector< te::rst::Grid > rastersGrids;
210 unsigned int inputRasterIdx = 0;
221 rastersGrids.push_back( (*inputRasterPtr->
getGrid()) );
225 if( inputRasterIdx == 0 )
242 auxLinearRingPtr->
setPoint( 0, mosaicLLX, mosaicURY );
243 auxLinearRingPtr->
setPoint( 1, mosaicURX, mosaicURY );
244 auxLinearRingPtr->
setPoint( 2, mosaicURX, mosaicLLY );
245 auxLinearRingPtr->
setPoint( 3, mosaicLLX, mosaicLLY );
246 auxLinearRingPtr->
setPoint( 4, mosaicLLX, mosaicURY );
247 auxLinearRingPtr->
setSRID( mosaicSRID );
250 auxPolygon.
push_back( auxLinearRingPtr );
251 auxPolygon.
setSRID( mosaicSRID );
252 rastersBBoxes.push_back( auxPolygon );
258 if( ( inputRasterIdx == 1 ) ||
269 const std::vector< te::gm::GTParameters::TiePoint >& inputTPs =
271 const unsigned int inputTPsSize = inputTPs.size();
273 (*eachRasterPixelToFirstRasterPixelGeomTransfms[ inputRasterIdx - 2 ].get());
275 for(
unsigned int inputTPsIdx = 0 ; inputTPsIdx < inputTPsSize ;
278 auxTP.second = inputTPs[ inputTPsIdx ].second;
279 lastTransf.
inverseMap( inputTPs[ inputTPsIdx ].first, auxTP.first );
287 boost::shared_ptr< te::gm::GeometricTransformation > auxTransPtr(
290 "Geometric transformation instatiation error" );
292 "Geometric transformation parameters calcule error" );
293 eachRasterPixelToFirstRasterPixelGeomTransfms.push_back( auxTransPtr );
299 *inputRasterPtr->
getGrid(), inRasterIndexedDetailedExtent ),
300 "Error creating the raster detailed extent" );
314 for(
unsigned int inRasterDetExtentIdx = 0 ; inRasterDetExtentIdx <
315 inRasterIndexedDetailedExtent.
size() ; ++inRasterDetExtentIdx )
317 auxTransPtr->inverseMap(
318 inRasterIndexedDetailedExtent.
getX( inRasterDetExtentIdx ),
319 inRasterIndexedDetailedExtent.
getY( inRasterDetExtentIdx ),
323 rastersGrids[ 0 ].gridToGeo( mappedX, mappedY, geoX,
326 inRasterDetailedExtent.
setPoint( inRasterDetExtentIdx, geoX, geoY );
342 auxPolygon.setSRID( mosaicSRID );
343 rastersBBoxes.push_back( auxPolygon );
348 for( std::vector< unsigned int >::size_type inputRastersBandsIdx = 0 ;
349 inputRastersBandsIdx <
351 ++inputRastersBandsIdx )
353 const unsigned int& currBand =
369 progressPtr->pulse();
370 if( ! progressPtr->isActive() )
return false;
378 mosaicBandsRangeMin.resize(
380 mosaicBandsRangeMax.resize(
383 std::vector< te::rst::BandProperty* > bandsProperties;
384 for( std::vector< unsigned int >::size_type bandIdx = 0 ; bandIdx <
392 mosaicBandsRangeMin[ bandIdx ],
393 mosaicBandsRangeMax[ bandIdx ]);
398 mosaicURY ), mosaicSRID );
409 "Output raster creation error" );
414 std::auto_ptr< te::mem::CachedRaster > cachedOutputRasterInstancePtr;
421 outputRasterPtr = cachedOutputRasterInstancePtr.get();
426 progressPtr->pulse();
427 if( ! progressPtr->isActive() )
return false;
434 std::vector< boost::shared_ptr< te::gm::GeometricTransformation > >
435 eachRasterPixelToMosaicRasterPixelGeomTransfms;
438 const double firstRasterColOffset = std::abs( rastersBBoxes[ 0 ].getMBR()->m_llx -
441 const double firstRasterLinOffset = std::abs( rastersBBoxes[ 0 ].getMBR()->m_ury -
451 const double prevRasterColOffset = std::abs( rastersBBoxes[ tiePointsIdx ].getMBR()->m_llx -
454 const double prevRasterLinOffset = std::abs( rastersBBoxes[ tiePointsIdx ].getMBR()->m_ury -
458 for(
unsigned int tpIdx = 0 ; tpIdx < transfParams.
m_tiePoints.size() ;
463 transfParams.
m_tiePoints[ tpIdx ].first.x += firstRasterColOffset;
464 transfParams.
m_tiePoints[ tpIdx ].first.y += firstRasterLinOffset;
468 transfParams.
m_tiePoints[ tpIdx ].first.x += prevRasterColOffset;
469 transfParams.
m_tiePoints[ tpIdx ].first.y += prevRasterLinOffset;
473 boost::shared_ptr< te::gm::GeometricTransformation > auxTransPtr(
476 "Geometric transformation instatiation error" );
478 "Geometric transformation parameters calcule error" );
479 eachRasterPixelToMosaicRasterPixelGeomTransfms.push_back( auxTransPtr );
485 progressPtr->pulse();
486 if( ! progressPtr->isActive() )
return false;
495 unsigned int col = 0;
496 unsigned int row = 0;
497 unsigned int bandIdx = 0;
499 for( bandIdx = 0 ; bandIdx < nBands ; ++bandIdx )
503 for( row = 0 ; row < nRows ; ++row )
505 for( col = 0 ; col < nCols ; ++col )
515 progressPtr->pulse();
516 if( ! progressPtr->isActive() )
return false;
522 std::vector< double > mosaicTargetMeans( outputRasterPtr->
getNumberOfBands(), 0 );
523 std::vector< double > mosaicTargetVariances( outputRasterPtr->
getNumberOfBands(), 0 );
532 double inXStartGeo = 0;
533 double inYStartGeo = 0;
535 double outRowStartDouble = 0;
536 double outColStartDouble = 0;
538 outColStartDouble, outRowStartDouble );
540 const unsigned int outRowStart = (
unsigned int)std::max( 0.0, outRowStartDouble );
541 const unsigned int outColStart = (
unsigned int)std::max( 0.0, outColStartDouble );
542 const unsigned int outRowsBound = std::min( outRowStart +
545 const unsigned int outColsBound = std::min( outColStart +
549 const unsigned int nBands = (
unsigned int)
551 unsigned int outCol = 0;
552 unsigned int outRow = 0;
555 double bandNoDataValue = -1.0 * DBL_MAX;
556 std::complex< double > pixelCValue = 0;
559 unsigned int inputBandIdx = 0;
561 for(
unsigned int inputRastersBandsIdx = 0 ; inputRastersBandsIdx <
562 nBands ; ++inputRastersBandsIdx )
565 inputRastersBandsIdx ] ;
570 (*outputRasterPtr->
getBand( inputRastersBandsIdx ));
571 unsigned int validPixelsNumber = 0;
573 double& mean = mosaicTargetMeans[ inputRastersBandsIdx ];
576 for( outRow = outRowStart ; outRow < outRowsBound ; ++outRow )
578 inRow = ((double)outRow) - outRowStartDouble;
580 for( outCol = outColStart ; outCol < outColsBound ; ++outCol )
582 inCol = ((double)outCol) - outColStartDouble;
584 interpInstance.
getValue( inCol, inRow, pixelCValue, inputBandIdx );
586 if( pixelCValue.real() != bandNoDataValue )
588 outBand.
setValue( outCol, outRow, pixelCValue );
589 mean += pixelCValue.real();
595 mean /= ( (double)validPixelsNumber );
601 double& variance = mosaicTargetVariances[ inputRastersBandsIdx ];
604 double pixelValue = 0;
606 for( outRow = outRowStart ; outRow < outRowsBound ; ++outRow )
608 for( outCol = outColStart ; outCol < outColsBound ; ++outCol )
610 outBand.
getValue( outCol, outRow, pixelValue );
614 variance += ( ( pixelValue - mean ) * ( pixelValue -
615 mean ) ) / ( (
double)validPixelsNumber );
625 progressPtr->pulse();
626 if( ! progressPtr->isActive() )
return false;
631 std::auto_ptr< te::gm::MultiPolygon > mosaicBBoxesUnionPtr(
633 rastersBBoxes[ 0 ].getSRID(), 0 ) );
635 mosaicBBoxesUnionPtr->add( (
te::gm::Polygon*)( rastersBBoxes[ 0 ].clone() ) );
639 std::vector< unsigned int > outputRasterBands;
640 std::vector< double > dummyRasterOffsets;
641 std::vector< double > dummyRasterScales;
642 for(
unsigned int outputRasterBandsIdx = 0 ; outputRasterBandsIdx <
645 outputRasterBands.push_back( outputRasterBandsIdx );
646 dummyRasterOffsets.push_back( 0.0 );
647 dummyRasterScales.push_back( 1.0 );
663 std::auto_ptr< te::mem::CachedRaster > cachedInputRasterPtr;
667 *nonCachedInputRasterPtr, 25, 0 ) );
668 inputRasterPtr = cachedInputRasterPtr.get();
673 std::vector< double > currentRasterBandsOffsets;
674 std::vector< double > currentRasterBandsScales;
678 double currentRasterVariance = 0;
679 double currentRasterMean = 0;
681 for(
unsigned int inputRastersBandsIdx = 0 ; inputRastersBandsIdx <
683 ++inputRastersBandsIdx )
686 inputRastersBandsIdx ];
692 currentRasterVariance );
694 currentRasterBandsScales.push_back(
695 std::sqrt( mosaicTargetVariances[ inputRastersBandsIdx ] )
697 std::sqrt( currentRasterVariance ) );
698 currentRasterBandsOffsets.push_back( mosaicTargetMeans[ inputRastersBandsIdx ]
699 - ( currentRasterBandsScales[ inputRastersBandsIdx ] * currentRasterMean ) );
704 currentRasterBandsOffsets = dummyRasterOffsets;
705 currentRasterBandsScales = dummyRasterScales;
725 currentRasterBandsOffsets,
726 currentRasterBandsScales,
727 mosaicBBoxesUnionPtr.get(),
729 *( eachRasterPixelToMosaicRasterPixelGeomTransfms[ inputRasterIdx - 1 ] ),
732 "Blender initiazing error" );
735 "Error blending images" );
743 std::auto_ptr< te::gm::Geometry > unionMultiPolPtr;
744 unionMultiPolPtr.reset( mosaicBBoxesUnionPtr->Union(
745 &( rastersBBoxes[ inputRasterIdx ] ) ) );
747 unionMultiPolPtr->setSRID( mosaicBBoxesUnionPtr->getSRID() );
756 unionMultiPolPtr->getSRID(), 0 );
758 mosaicBBoxesUnionPtr.reset( mPolPtr );
768 cachedInputRasterPtr.reset();
774 progressPtr->pulse();
775 if( ! progressPtr->isActive() )
return false;
789 throw( te::rp::Exception )
802 "Invalid m_feederRasterPtr" )
806 "Invalid number of rasters" )
815 for( std::vector< std::vector< unsigned int > >::size_type
816 inputRastersBandsIdx = 0 ; inputRastersBandsIdx <
820 inputRastersBandsIdx ].size() > 0,
"Invalid bands number" );
824 0 ].size(),
"Bands number mismatch" );
845 const bool& forceNoDataValue,
846 const double& noDataValue,
847 double& mean,
double& variance )
852 double internalNoDataValue = 0;
853 if( forceNoDataValue )
854 internalNoDataValue = noDataValue;
863 double pixelsNumber = 0;
865 unsigned int col = 0;
866 unsigned int line = 0;
868 for( line = 0 ; line < nLines ; ++line )
869 for( col = 0 ; col < nCols ; ++col )
873 if( value != internalNoDataValue )
880 if( pixelsNumber != 0.0 )
882 mean /= pixelsNumber;
884 for( line = 0 ; line < nLines ; ++line )
885 for( col = 0 ; col < nCols ; ++col )
889 if( value != internalNoDataValue )
891 variance += ( ( value - mean ) * ( value - mean ) ) / pixelsNumber;
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.
std::map< std::string, std::string > m_rInfo
The necessary information to create the output rasters (as described in te::raster::RasterFactory).
Near neighborhood interpolation method.
const OutputParameters & operator=(const OutputParameters ¶ms)
std::vector< TiePoint > m_tiePoints
Tie points.
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 reset()
Clear all internal allocated objects and reset the algorithm to its initial state.
void setSRID(int srid)
It sets the Spatial Reference System ID of the geometry and all its parts if it is a GeometryCollecti...
bool initialize(const AlgorithmInputParameters &inputParams)
Initialize the algorithm instance making it ready for execution.
A raster band description.
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.
AbstractParameters * clone() const
Create a clone copy of this instance.
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.
static void calcBandStatistics(const te::rst::Band &band, const bool &forceNoDataValue, const double &noDataValue, double &mean, double &variance)
Raster band statistics calcule.
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.
bool m_isInitialized
Tells if this instance is initialized.
std::pair< Coord2D, Coord2D > TiePoint
Tie point type definition.
void geoToGrid(const double &x, const double &y, double &col, double &row) const
Get the grid point associated to a spatial location.
const double & getY(std::size_t i) const
It returns the n-th y coordinate value.
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.
bool execute(AlgorithmOutputParameters &outputParams)
Executes the algorithm using the supplied parameters.
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.
const double & getX(std::size_t i) const
It returns the n-th x coordinate value.
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.
double getResolutionX() const
Returns the grid horizontal (x-axis) resolution.
A RAM cache adaptor to an external existent raster that must always be avaliable. ...
std::string m_rType
Output raster data source type (as described in te::raster::RasterFactory ).
A raster band description.
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.
#define TERP_LOG_AND_RETURN_FALSE(message)
Logs a warning message will and return false.
bool isInitialized() const
Returns true if the algorithm instance is initialized and ready for execution.
Polygon is a subclass of CurvePolygon whose rings are defined by linear rings.
double m_ury
Upper right corner y-coordinate.
void add(Geometry *g)
It adds the geometry into the collection.
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.
Mosaic output parameters.
bool GetIndexedDetailedExtent(const te::rst::Grid &grid, te::gm::LinearRing &indexedDetailedExtent)
Create a indexed (lines,columns) datailed extent from the given grid.
#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 gridToGeo(const double &col, const double &row, double &x, double &y) const
Get the spatial location of a grid point.
std::auto_ptr< te::rst::Raster > m_outputRasterPtr
The generated output mosaic raster.
void setSRID(int srid)
It sets the Spatial Reference System ID of the linestring.
void reset()
Clear all internal allocated resources and reset the parameters instance to its initial state...
2D Geometric transformation parameters.
Create a mosaic from a set of rasters using tie-points.
TiePointsMosaic::InputParameters m_inputParameters
Input execution parameters.
virtual bool moveNext()=0
Advances to the next sequence obeject.
int m_blkh
Block height (pixels).
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.
std::size_t size() const
It returns the number of points (vertexes) in the geometry.