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;
221 if( inputRasterIdx == 0 )
238 auxLinearRingPtr->
setPoint( 0, mosaicLLX, mosaicURY );
239 auxLinearRingPtr->
setPoint( 1, mosaicURX, mosaicURY );
240 auxLinearRingPtr->
setPoint( 2, mosaicURX, mosaicLLY );
241 auxLinearRingPtr->
setPoint( 3, mosaicLLX, mosaicLLY );
242 auxLinearRingPtr->
setPoint( 4, mosaicLLX, mosaicURY );
243 auxPolygon.
push_back( auxLinearRingPtr );
244 auxPolygon.
setSRID( mosaicSRID );
245 rastersBBoxes.push_back( auxPolygon );
275 mosaicLLX = std::min( mosaicLLX, urCoord1.
x );
276 mosaicLLX = std::min( mosaicLLX, llCoord1.
x );
278 mosaicLLY = std::min( mosaicLLY, urCoord1.
y );
279 mosaicLLY = std::min( mosaicLLY, llCoord1.
y );
281 mosaicURX = std::max( mosaicURX, urCoord1.
x );
282 mosaicURX = std::max( mosaicURX, llCoord1.
x );
284 mosaicURY = std::max( mosaicURY, urCoord1.
y );
285 mosaicURY = std::max( mosaicURY, llCoord1.
y );
291 auxLinearRingPtr->
setPoint( 0, llCoord1.
x, urCoord1.
y );
292 auxLinearRingPtr->
setPoint( 1, urCoord1.
x, urCoord1.
y );
293 auxLinearRingPtr->
setPoint( 2, urCoord1.
x, llCoord1.
y );
294 auxLinearRingPtr->
setPoint( 3, llCoord1.
x, llCoord1.
y );
295 auxLinearRingPtr->
setPoint( 4, llCoord1.
x, urCoord1.
y );
296 auxPolygon.
push_back( auxLinearRingPtr );
297 auxPolygon.
setSRID( mosaicSRID );
298 rastersBBoxes.push_back( auxPolygon );
303 for( std::vector< unsigned int >::size_type inputRastersBandsIdx = 0 ;
304 inputRastersBandsIdx <
306 ++inputRastersBandsIdx )
308 const unsigned int& currBand =
322 progressPtr->pulse();
323 if( ! progressPtr->isActive() )
return false;
331 mosaicBandsRangeMin.resize(
333 mosaicBandsRangeMax.resize(
336 std::vector< te::rst::BandProperty* > bandsProperties;
337 for( std::vector< unsigned int >::size_type bandIdx = 0 ; bandIdx <
345 mosaicBandsRangeMin[ bandIdx ],
346 mosaicBandsRangeMax[ bandIdx ]);
351 mosaicURY ), mosaicSRID );
362 "Output raster creation error" );
367 std::auto_ptr< te::mem::CachedRaster > cachedRasterInstancePtr;
374 outputRasterPtr = cachedRasterInstancePtr.get();
383 unsigned int col = 0;
384 unsigned int row = 0;
385 unsigned int bandIdx = 0;
387 for( bandIdx = 0 ; bandIdx < nBands ; ++bandIdx )
391 for( row = 0 ; row < nRows ; ++row )
393 for( col = 0 ; col < nCols ; ++col )
403 progressPtr->pulse();
404 if( ! progressPtr->isActive() )
return false;
410 std::vector< double > mosaicTargetMeans( outputRasterPtr->
getNumberOfBands(), 0 );
411 std::vector< double > mosaicTargetVariances( outputRasterPtr->
getNumberOfBands(), 0 );
420 double inXStartGeo = 0;
421 double inYStartGeo = 0;
423 double outRowStartDouble = 0;
424 double outColStartDouble = 0;
426 outColStartDouble, outRowStartDouble );
428 const unsigned int outRowStart = (
unsigned int)std::max( 0.0, outRowStartDouble );
429 const unsigned int outColStart = (
unsigned int)std::max( 0.0, outColStartDouble );
430 const unsigned int outRowsBound = std::min( outRowStart +
433 const unsigned int outColsBound = std::min( outColStart +
437 const unsigned int nBands = (
unsigned int)
439 unsigned int outCol = 0;
440 unsigned int outRow = 0;
443 double bandNoDataValue = -1.0 * DBL_MAX;
444 std::complex< double > pixelCValue = 0;
447 unsigned int inputBandIdx = 0;
449 for(
unsigned int inputRastersBandsIdx = 0 ; inputRastersBandsIdx <
450 nBands ; ++inputRastersBandsIdx )
453 inputRastersBandsIdx ] ;
458 (*outputRasterPtr->
getBand( inputRastersBandsIdx ));
459 unsigned int validPixelsNumber = 0;
461 double& mean = mosaicTargetMeans[ inputRastersBandsIdx ];
464 for( outRow = outRowStart ; outRow < outRowsBound ; ++outRow )
466 inRow = ((double)outRow) - outRowStartDouble;
468 for( outCol = outColStart ; outCol < outColsBound ; ++outCol )
470 inCol = ((double)outCol) - outColStartDouble;
472 interpInstance.
getValue( inCol, inRow, pixelCValue, inputBandIdx );
474 if( pixelCValue.real() != bandNoDataValue )
476 outBand.
setValue( outCol, outRow, pixelCValue );
477 mean += pixelCValue.real();
483 mean /= ( (double)validPixelsNumber );
489 double& variance = mosaicTargetVariances[ inputRastersBandsIdx ];
492 double pixelValue = 0;
494 for( outRow = outRowStart ; outRow < outRowsBound ; ++outRow )
496 for( outCol = outColStart ; outCol < outColsBound ; ++outCol )
498 outBand.
getValue( outCol, outRow, pixelValue );
502 variance += ( ( pixelValue - mean ) * ( pixelValue -
503 mean ) ) / ( (
double)validPixelsNumber );
513 "Rasters bounding boxes number mismatch" );
517 progressPtr->pulse();
518 if( ! progressPtr->isActive() )
return false;
523 std::auto_ptr< te::gm::MultiPolygon > mosaicBBoxesUnionPtr(
525 outputRasterPtr->
getSRID(), 0 ) );
526 mosaicBBoxesUnionPtr->add( (
te::gm::Polygon*)rastersBBoxes[ 0 ].clone() );
537 std::vector< unsigned int > outputRasterBands;
538 std::vector< double > dummyRasterOffsets;
539 std::vector< double > dummyRasterScales;
540 for(
unsigned int outputRasterBandsIdx = 0 ; outputRasterBandsIdx <
543 outputRasterBands.push_back( outputRasterBandsIdx );
544 dummyRasterOffsets.push_back( 0.0 );
545 dummyRasterScales.push_back( 1.0 );
556 "Invalid boxes SRID" );
565 std::auto_ptr< te::rst::Raster > reprojectedInputRasterPtr;
567 std::string reprojectedRasterFileName;
569 if( outputRasterPtr->
getSRID() != inputRasterPtr->getSRID() )
571 reprojectedRasterFileName = boost::filesystem::unique_path(
572 boost::filesystem::temp_directory_path() /=
573 boost::filesystem::path(
"TerralibRGeoMosaic_%%%%-%%%%-%%%%-%%%%" ) ).string();
575 "Invalid temporary raster file name" );
576 reprojectedRasterFileName +=
".tif";
578 std::map< std::string, std::string > rInfo;
579 rInfo[
"URI" ] = reprojectedRasterFileName;
581 reprojectedInputRasterPtr.reset( inputRasterPtr->transform(
582 outputRasterPtr->
getSRID(), rInfo,
584 inputRasterPtr = reprojectedInputRasterPtr.get();
593 std::auto_ptr< te::mem::CachedRaster > cachedInputRasterPtr;
598 *inputRasterPtr, 20, 0 ) );
599 inputRasterPtr = cachedInputRasterPtr.get();
607 std::auto_ptr< te::gm::GeometricTransformation > transPtr;
617 outputRasterPtr->
getGrid()->
gridToGeo( auxTP.first.x, auxTP.first.y, auxX, auxY );
618 inputRasterPtr->getGrid()->geoToGrid( auxX, auxY, auxTP.second.x, auxTP.second.y );
623 outputRasterPtr->
getGrid()->
gridToGeo( auxTP.first.x, auxTP.first.y, auxX, auxY );
624 inputRasterPtr->getGrid()->geoToGrid( auxX, auxY, auxTP.second.x, auxTP.second.y );
629 outputRasterPtr->
getGrid()->
gridToGeo( auxTP.first.x, auxTP.first.y, auxX, auxY );
630 inputRasterPtr->getGrid()->geoToGrid( auxX, auxY, auxTP.second.x, auxTP.second.y );
635 outputRasterPtr->
getGrid()->
gridToGeo( auxTP.first.x, auxTP.first.y, auxX, auxY );
636 inputRasterPtr->getGrid()->geoToGrid( auxX, auxY, auxTP.second.x, auxTP.second.y );
643 "Could not initialize a geometric transformation" );
648 std::vector< double > currentRasterBandsOffsets;
649 std::vector< double > currentRasterBandsScales;
653 double currentRasterVariance = 0;
654 double currentRasterMean = 0;
656 for(
unsigned int inputRastersBandsIdx = 0 ; inputRastersBandsIdx <
658 ++inputRastersBandsIdx )
661 inputRastersBandsIdx ];
667 currentRasterVariance );
669 currentRasterBandsScales.push_back( std::sqrt( mosaicTargetVariances[ inputRastersBandsIdx ] /
670 currentRasterVariance ) );
671 currentRasterBandsOffsets.push_back( mosaicTargetMeans[ inputRastersBandsIdx ] -
672 ( currentRasterBandsScales[ inputRastersBandsIdx ] * currentRasterMean ) );
677 currentRasterBandsOffsets = dummyRasterOffsets;
678 currentRasterBandsScales = dummyRasterScales;
697 currentRasterBandsOffsets,
698 currentRasterBandsScales,
699 mosaicBBoxesUnionPtr.get(),
703 false ),
"Blender initiazing error" );
706 "Error blending images" );
714 std::auto_ptr< te::gm::Geometry > boxesUnionResultPtr;
715 boxesUnionResultPtr.reset( mosaicBBoxesUnionPtr->Union(
716 &( rastersBBoxes[ inputRasterIdx ] ) ) );
718 boxesUnionResultPtr->setSRID( outputRasterPtr->
getSRID() );
729 auxMultiPol->
add( boxesUnionResultPtr.release() );
731 mosaicBBoxesUnionPtr.reset( auxMultiPol );
740 cachedInputRasterPtr.reset();
744 if( reprojectedInputRasterPtr.get() )
746 reprojectedInputRasterPtr.reset();
747 remove( reprojectedRasterFileName.c_str() );
756 progressPtr->pulse();
757 if( ! progressPtr->isActive() )
return false;
763 if( cachedRasterInstancePtr.get() ) cachedRasterInstancePtr.reset();
775 throw( te::rp::Exception )
788 "Invalid m_feederRasterPtr" )
792 "Invalid number of rasters" )
801 for( std::vector< std::vector< unsigned int > >::size_type
802 inputRastersBandsIdx = 0 ; inputRastersBandsIdx <
806 inputRastersBandsIdx ].size() > 0,
"Invalid bands number" );
810 0 ].size(),
"Bands number mismatch" );
824 const bool& forceNoDataValue,
825 const double& noDataValue,
826 double& mean,
double& variance )
831 double internalNoDataValue = 0;
832 if( forceNoDataValue )
833 internalNoDataValue = noDataValue;
842 double pixelsNumber = 0;
844 unsigned int col = 0;
845 unsigned int line = 0;
847 for( line = 0 ; line < nLines ; ++line )
848 for( col = 0 ; col < nCols ; ++col )
852 if( value != internalNoDataValue )
859 if( pixelsNumber != 0.0 )
861 mean /= pixelsNumber;
863 for( line = 0 ; line < nLines ; ++line )
864 for( col = 0 ; col < nCols ; ++col )
868 if( value != internalNoDataValue )
870 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.
bool convert(double *xIn, double *yIn, double *xOut, double *yOut, long numCoord, int coordOffset=1) const
Converts a vector of coordinates from source SRS to target SRS.
MultiPolygon is a MultiSurface whose elements are Polygons.
bool initialize(const AlgorithmInputParameters &inputParams)
Initialize the algorithm instance making it ready for execution.
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.
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.
An utility struct for representing 2D coordinates.
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().
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.
void setSourceSRID(int sourceSRID)
Sets the source SRS identifier.
#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.
void setTargetSRID(int targetSRID)
Sets the target SRS identifier.
A RAM cache adaptor to an external existent raster that must always be avaliable. ...
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.
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.
#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.
2D Geometric transformation parameters.
virtual bool moveNext()=0
Advances to the next sequence obeject.
int m_blkh
Block height (pixels).
Near neighborhood interpolation method.
#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.
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.
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 forceInputNoDataValue, 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.
std::map< std::string, std::string > m_rInfo
The necessary information to create the output rasters (as described in te::raster::RasterFactory).