23#include <gsl/gsl_cdf.h>
25#include "ImgReaderGdal.h"
54#if GDAL_VERSION_MAJOR < 2
59 if(readMode==GA_ReadOnly)
60 m_gds = (GDALDataset*) GDALOpenEx(
m_filename.c_str(), GDAL_OF_READONLY|GDAL_OF_RASTER, NULL, NULL, NULL);
61 else if(readMode==GA_Update)
62 m_gds = (GDALDataset*) GDALOpenEx(
m_filename.c_str(), GDAL_OF_UPDATE|GDAL_OF_RASTER, NULL, NULL, NULL);
66 std::string errorString=
"FileOpenError";
72 double adfGeoTransform[6];
73 m_gds->GetGeoTransform( adfGeoTransform );
74 m_gt[0]=adfGeoTransform[0];
75 m_gt[1]=adfGeoTransform[1];
76 m_gt[2]=adfGeoTransform[2];
77 m_gt[3]=adfGeoTransform[3];
78 m_gt[4]=adfGeoTransform[4];
79 m_gt[5]=adfGeoTransform[5];
80 m_projection=
m_gds->GetProjectionRef();
91 std::vector<double> lineBuffer(
nrOfCol());
93 for(
int irow=0;irow<
nrOfRow();++irow){
95 for(
int icol=0;icol<
nrOfCol();++icol){
99 if(lineBuffer[icol]<minValue){
102 minValue=lineBuffer[icol];
108 minValue=lineBuffer[icol];
116 throw(
static_cast<std::string
>(
"Warning: not initialized"));
127 std::vector<double> lineBuffer(
nrOfCol());
129 for(
int irow=0;irow<
nrOfRow();++irow){
131 for(
int icol=0;icol<
nrOfCol();++icol){
135 if(lineBuffer[icol]>maxValue){
138 maxValue=lineBuffer[icol];
144 maxValue=lineBuffer[icol];
152 throw(
static_cast<std::string
>(
"Warning: not initialized"));
163 bool isConstraint=(maxValue>minValue);
164 double minConstraint=minValue;
165 double maxConstraint=maxValue;
166 std::vector<double> lineBuffer(endCol-startCol+1);
169 for(
int irow=startCol;irow<endRow+1;++irow){
170 readData(lineBuffer,startCol,endCol,irow,band);
171 for(
int icol=0;icol<lineBuffer.size();++icol){
176 if(lineBuffer[icol]<minConstraint)
178 if(lineBuffer[icol]>maxConstraint)
181 if(lineBuffer[icol]<minValue)
182 minValue=lineBuffer[icol];
183 if(lineBuffer[icol]>maxValue)
184 maxValue=lineBuffer[icol];
188 if(lineBuffer[icol]<minConstraint)
190 if(lineBuffer[icol]>maxConstraint)
193 minValue=lineBuffer[icol];
194 maxValue=lineBuffer[icol];
200 throw(
static_cast<std::string
>(
"Warning: not initialized"));
210 bool isConstraint=(maxValue>minValue);
211 double minConstraint=minValue;
212 double maxConstraint=maxValue;
213 std::vector<double> lineBuffer(
nrOfCol());
215 for(
int irow=0;irow<
nrOfRow();++irow){
217 for(
int icol=0;icol<
nrOfCol();++icol){
222 if(lineBuffer[icol]<minConstraint)
224 if(lineBuffer[icol]>maxConstraint)
227 if(lineBuffer[icol]<minValue)
228 minValue=lineBuffer[icol];
229 if(lineBuffer[icol]>maxValue)
230 maxValue=lineBuffer[icol];
234 if(lineBuffer[icol]<minConstraint)
236 if(lineBuffer[icol]>maxConstraint)
239 minValue=lineBuffer[icol];
240 maxValue=lineBuffer[icol];
246 throw(
static_cast<std::string
>(
"Warning: not initialized"));
267 if(min<max&&min>minValue)
269 if(min<max&&max<maxValue)
278 GDALProgressFunc pfnProgress;
280 GDALRasterBand* rasterBand;
282 rasterBand->ComputeStatistics(0,&minValue,&maxValue,&meanValue,&stdDev,pfnProgress,pProgressData);
290 sigma=1.06*stdDev*pow(
getNvalid(theBand),-0.2);
294 if(maxValue>minValue){
296 nbin=maxValue-minValue+1;
297 scale=
static_cast<double>(nbin-1)/(maxValue-minValue);
302 if(histvector.size()!=nbin){
303 histvector.resize(nbin);
304 for(
int i=0;i<nbin;histvector[i++]=0);
307 unsigned long int nsample=0;
308 unsigned long int ninvalid=0;
309 std::vector<double> lineBuffer(
nrOfCol());
310 for(
int irow=0;irow<
nrOfRow();++irow){
312 for(
int icol=0;icol<
nrOfCol();++icol){
315 else if(lineBuffer[icol]>maxValue)
317 else if(lineBuffer[icol]<minValue)
326 for(
int ibin=0;ibin<nbin;++ibin){
327 double icenter=minValue+
static_cast<double>(maxValue-minValue)*(ibin+0.5)/nbin;
328 double thePdf=gsl_ran_gaussian_pdf(lineBuffer[icol]-icenter, sigma);
329 histvector[ibin]+=thePdf;
334 int theBin=
static_cast<unsigned long int>(scale*(lineBuffer[icol]-minValue));
337 ++histvector[theBin];
357 std::vector<short> lineBuffer(
nrOfCol());
359 for(
int irow=0;irow<
nrOfRow();++irow){
361 for(
int icol=0;icol<
nrOfCol();++icol){
362 if(find(range.begin(),range.end(),lineBuffer[icol])==range.end())
363 range.push_back(lineBuffer[icol]);
366 sort(range.begin(),range.end());
375 unsigned long int nvalid=0;
377 std::vector<double> lineBuffer(
nrOfCol());
378 for(
int irow=0;irow<
nrOfRow();++irow){
380 for(
int icol=0;icol<
nrOfCol();++icol){
400 std::vector<double> lineBuffer(
nrOfCol());
405 for(
int irow=0;irow<
nrOfRow();++irow){
407 for(
int icol=0;icol<
nrOfCol();++icol){
421 double cgravi=validCol/nvalidCol-0.5;
422 double cgravj=validRow/nvalidRow-0.5;
423 double refpixeli=floor(cgravi);
424 double refpixelj=ceil(cgravj-1);
426 image2geo(refpixeli,refpixelj,refX,refY);
432 refX=floor(validCol/nvalidCol-0.5);
433 refY=floor(validRow/nvalidRow-0.5);
double getDeltaY(void) const
Get the pixel cell spacing in y.
virtual void close(void)
Close the image.
int nrOfRow(void) const
Get the number of rows of this dataset.
bool isGeoRef() const
Is this dataset georeferenced (pixel size in y must be negative) ?
int m_nband
number of bands in this dataset
int m_ncol
number of columns in this dataset
bool isNoData(double value) const
Check if value is nodata in this dataset.
std::string m_filename
filename of this dataset
GDALDataset * m_gds
instance of the GDAL dataset of this dataset
int m_nrow
number of rows in this dataset
int nrOfCol(void) const
Get the number of columns of this dataset.
double getDeltaX(void) const
Get the pixel cell spacing in x.
std::vector< double > m_noDataValues
no data values for this dataset
GDALRasterBand * getRasterBand(int band=0) const
Get the GDAL rasterband for this dataset.
std::vector< double > m_scale
Vector containing the scale factor to be applied (one scale value for each band)
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y)
double m_gt[6]
geotransform information of this dataset
void setCodec(const GDALAccess &readMode=GA_ReadOnly)
Set GDAL dataset number of columns, rows, bands and geotransform.
unsigned long int getNvalid(int band)
Calculate the number of valid pixels (with a value not defined as no data).
~ImgReaderGdal(void)
destructor
void readData(T &value, int col, int row, int band=0)
Read a single pixel cell value at a specific column and row for a specific band (all indices start co...
double getMax(int &col, int &row, int band=0)
Get the maximum cell values for a specific band and report the column and row in which the maximum va...
void getMinMax(int startCol, int endCol, int startRow, int endRow, int band, double &minValue, double &maxValue)
Get the minimum and maximum cell values for a specific band in a region of interest defined by startC...
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
void getRefPix(double &refX, double &refY, int band=0)
Calculate the reference pixel as the centre of gravity pixel (weighted average of all values not taki...
ImgReaderGdal(void)
default constructor. Image needs to be opened later with one of the open methods.
void getRange(std::vector< short > &range, int Band=0)
Calculate the range of cell values in the image for a specific band (start counting from 0).
double getHistogram(std::vector< double > &histvector, double &min, double &max, unsigned int &nbin, int theBand=0, bool kde=false)
Calculate the image histogram for a specific band using a defined number of bins and constrained by a...
double getMin(int &col, int &row, int band=0)
Get the minimum cell values for a specific band and report the column and row in which the minimum va...