[pktools] 161/375: pkfilter 2d wavelet forward and inverse, pkinfo hist when min=max
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:54:10 UTC 2014
This is an automated email from the git hooks/post-receive script.
sebastic-guest pushed a commit to branch upstream-master
in repository pktools.
commit 0ae2e5812c4d3d775f0857366419466b5f2ef05d
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Fri Dec 20 23:36:19 2013 +0100
pkfilter 2d wavelet forward and inverse, pkinfo hist when min=max
---
ChangeLog | 5 ++
src/algorithms/Filter2d.cc | 11 ++++
src/algorithms/Filter2d.h | 80 +++++++++++++++++++++++--
src/apps/pkcrop.cc | 4 +-
src/apps/pkdiff.cc | 72 +++++++++++++----------
src/apps/pkfilter.cc | 3 +
src/apps/pkfilterascii.cc | 5 +-
src/apps/pkinfo.cc | 12 +++-
src/apps/pkmosaic.cc | 4 ++
src/imageclasses/ImgReaderGdal.cc | 119 +++++++++++++++++++++++++++-----------
src/imageclasses/ImgWriterGdal.cc | 103 ++++++++++++++++++++++++---------
11 files changed, 316 insertions(+), 102 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index 3f3e67d..ae67916 100755
--- a/ChangeLog
+++ b/ChangeLog
@@ -208,3 +208,8 @@ version 2.4.3
- pksetmask
option msknodata to set nodata value in mask
+version 2.4.4
+ - pkinfo
+ hist: corrected nan when min=max
+ - pkfilter
+ corrected bug in 2d wavelet transform forward and inverse
diff --git a/src/algorithms/Filter2d.cc b/src/algorithms/Filter2d.cc
index e11c021..91d2695 100644
--- a/src/algorithms/Filter2d.cc
+++ b/src/algorithms/Filter2d.cc
@@ -1097,11 +1097,22 @@ void filter2d::Filter2d::dwtForward(const ImgReaderGdal& input, ImgWriterGdal& o
Vector2d<float> theBuffer;
for(int iband=0;iband<input.nrOfBand();++iband){
input.readDataBlock(theBuffer, GDT_Float32, 0, input.nrOfCol()-1, 0, input.nrOfRow()-1, iband);
+ std::cout << "filtering band " << iband << std::endl << std::flush;
dwtForward(theBuffer, wavelet_type, family);
output.writeDataBlock(theBuffer,GDT_Float32,0,output.nrOfCol()-1,0,output.nrOfRow()-1,iband);
}
}
+void filter2d::Filter2d::dwtInverse(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family){
+ Vector2d<float> theBuffer;
+ for(int iband=0;iband<input.nrOfBand();++iband){
+ input.readDataBlock(theBuffer, GDT_Float32, 0, input.nrOfCol()-1, 0, input.nrOfRow()-1, iband);
+ std::cout << "filtering band " << iband << std::endl << std::flush;
+ dwtInverse(theBuffer, wavelet_type, family);
+ output.writeDataBlock(theBuffer,GDT_Float32,0,output.nrOfCol()-1,0,output.nrOfRow()-1,iband);
+ }
+}
+
void filter2d::Filter2d::dwtQuantize(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, double quantize, bool verbose){
Vector2d<float> theBuffer;
for(int iband=0;iband<input.nrOfBand();++iband){
diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h
index 905f596..1718fda 100644
--- a/src/algorithms/Filter2d.h
+++ b/src/algorithms/Filter2d.h
@@ -95,10 +95,11 @@ public:
template<class T1, class T2> void smooth(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector,int dim);
template<class T1, class T2> void smooth(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector,int dimX, int dimY);
void dwtForward(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family);
+ void dwtInverse(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family);
void dwtQuantize(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, double quantize, bool verbose=false);
template<class T> void dwtForward(Vector2d<T>& data, const std::string& wavelet_type, int family);
- template<class T> void dwtQuantize(Vector2d<T>& data, const std::string& wavelet_type, int family, double quantize);
template<class T> void dwtInverse(Vector2d<T>& data, const std::string& wavelet_type, int family);
+ template<class T> void dwtQuantize(Vector2d<T>& data, const std::string& wavelet_type, int family, double quantize);
void majorVoting(const std::string& inputFilename, const std::string& outputFilename,int dim=0,const std::vector<int> &prior=std::vector<int>());
/* void homogeneousSpatial(const std::string& inputFilename, const std::string& outputFilename, int dim, bool disc=false, int noValue=0); */
void doit(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down=2, bool disc=false);
@@ -121,8 +122,7 @@ public:
private:
static void initMap(std::map<std::string, FILTER_TYPE>& m_filterMap){
//initialize selMap
- m_filterMap["mrf"]=filter2d::mrf;
- m_filterMap["stdev"]=filter2d::stdev;
+ m_filterMap["median"]=filter2d::median;
m_filterMap["var"]=filter2d::var;
m_filterMap["min"]=filter2d::min;
m_filterMap["max"]=filter2d::max;
@@ -148,7 +148,10 @@ private:
m_filterMap["ismax"]=filter2d::ismax;
m_filterMap["heterog"]=filter2d::heterog;
m_filterMap["order"]=filter2d::order;
- m_filterMap["median"]=filter2d::median;
+ m_filterMap["stdev"]=filter2d::stdev;
+ m_filterMap["mrf"]=filter2d::mrf;
+ m_filterMap["dwtForward"]=filter2d::dwtForward;
+ m_filterMap["dwtInverse"]=filter2d::dwtInverse;
m_filterMap["dwtQuantize"]=filter2d::dwtQuantize;
m_filterMap["scramble"]=filter2d::scramble;
m_filterMap["shift"]=filter2d::shift;
@@ -772,6 +775,12 @@ template<class T> unsigned long int Filter2d::morphology(const Vector2d<T>& inpu
}
template<class T> void Filter2d::dwtForward(Vector2d<T>& theBuffer, const std::string& wavelet_type, int family){
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+
//make sure data size if power of 2
int nRow=theBuffer.size();
assert(nRow);
@@ -803,10 +812,67 @@ template<class T> void Filter2d::dwtForward(Vector2d<T>& theBuffer, const std::s
int index=irow*theBuffer[irow].size()+icol;
theBuffer[irow][icol]=data[index];
}
+ progress=(1.0+irow);
+ progress/=theBuffer.nRows();
+ pfnProgress(progress,pszMessage,pProgressArg);
}
+ gsl_wavelet_free (w);
+ gsl_wavelet_workspace_free (work);
+}
+
+template<class T> void Filter2d::dwtInverse(Vector2d<T>& theBuffer, const std::string& wavelet_type, int family){
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+
+ //make sure data size if power of 2
+ int nRow=theBuffer.size();
+ assert(nRow);
+ int nCol=theBuffer[0].size();
+ assert(nCol);
+ while(theBuffer.size()&(theBuffer.size()-1))
+ theBuffer.push_back(theBuffer.back());
+ for(int irow=0;irow<theBuffer.size();++irow)
+ while(theBuffer[irow].size()&(theBuffer[irow].size()-1))
+ theBuffer[irow].push_back(theBuffer[irow].back());
+ double data[theBuffer.size()*theBuffer[0].size()];
+ for(int irow=0;irow<theBuffer.size();++irow){
+ for(int icol=0;icol<theBuffer[0].size();++icol){
+ int index=irow*theBuffer[0].size()+icol;
+ data[index]=theBuffer[irow][icol];
+ }
+ }
+ int nsize=theBuffer.size()*theBuffer[0].size();
+ gsl_wavelet *w;
+ gsl_wavelet_workspace *work;
+ assert(nsize);
+ w=gsl_wavelet_alloc(filter::Filter::getWaveletType(wavelet_type),family);
+ work=gsl_wavelet_workspace_alloc(nsize);
+ gsl_wavelet2d_nstransform_inverse (w, data, theBuffer.size(), theBuffer.size(),theBuffer[0].size(), work);
+ theBuffer.erase(theBuffer.begin()+nRow,theBuffer.end());
+ for(int irow=0;irow<theBuffer.size();++irow){
+ theBuffer[irow].erase(theBuffer[irow].begin()+nCol,theBuffer[irow].end());
+ for(int icol=0;icol<theBuffer[irow].size();++icol){
+ int index=irow*theBuffer[irow].size()+icol;
+ theBuffer[irow][icol]=data[index];
+ }
+ progress=(1.0+irow);
+ progress/=theBuffer.nRows();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+ gsl_wavelet_free (w);
+ gsl_wavelet_workspace_free (work);
}
template<class T> void Filter2d::dwtQuantize(Vector2d<T>& theBuffer, const std::string& wavelet_type, int family, double quantize){
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+
//make sure data size if power of 2
int nRow=theBuffer.size();
assert(nRow);
@@ -856,6 +922,9 @@ template<class T> void Filter2d::dwtQuantize(Vector2d<T>& theBuffer, const std::
int index=irow*theBuffer[irow].size()+icol;
theBuffer[irow][icol]=data[index];
}
+ progress=(1.0+irow);
+ progress/=theBuffer.nRows();
+ pfnProgress(progress,pszMessage,pProgressArg);
}
theBuffer.erase(theBuffer.begin()+nRow,theBuffer.end());
for(int irow=0;irow<theBuffer.size();++irow)
@@ -863,6 +932,9 @@ template<class T> void Filter2d::dwtQuantize(Vector2d<T>& theBuffer, const std::
delete[] data;
delete[] abscoeff;
delete[] p;
+ gsl_wavelet_free (w);
+ gsl_wavelet_workspace_free (work);
+
}
}
diff --git a/src/apps/pkcrop.cc b/src/apps/pkcrop.cc
index e95a3ae..e1cd306 100644
--- a/src/apps/pkcrop.cc
+++ b/src/apps/pkcrop.cc
@@ -277,8 +277,8 @@ int main(int argc, char *argv[])
}
imgReader.geo2image(cropulx+(magicX-1.0)*imgReader.getDeltaX(),cropuly-(magicY-1.0)*imgReader.getDeltaY(),uli,ulj);
imgReader.geo2image(croplrx+(magicX-2.0)*imgReader.getDeltaX(),croplry-(magicY-2.0)*imgReader.getDeltaY(),lri,lrj);
- ncropcol=ceil((croplrx-cropulx)/dx);
- ncroprow=ceil((cropuly-croplry)/dy);
+ // ncropcol=ceil((croplrx-cropulx)/dx);
+ // ncroprow=ceil((cropuly-croplry)/dy);
}
else{
double magicX=1,magicY=1;
diff --git a/src/apps/pkdiff.cc b/src/apps/pkdiff.cc
index aa46af2..4c6d011 100644
--- a/src/apps/pkdiff.cc
+++ b/src/apps/pkdiff.cc
@@ -36,7 +36,7 @@ int main(int argc, char *argv[])
Optionpk<short> valueE_opt("\0", "correct", "Value for correct pixels (0)", 0,1);
Optionpk<short> valueO_opt("\0", "omission", "Value for omission errors: input label > reference label (default value is 1)", 1,1);
Optionpk<short> valueC_opt("\0", "commission", "Value for commission errors: input label < reference label (default value is 2)", 2,1);
- Optionpk<short> nodata_opt("nodata", "nodata", "No value flag(s)", 0);
+ Optionpk<short> nodata_opt("nodata", "nodata", "No value flag(s)");
Optionpk<short> band_opt("b", "band", "Band to extract (0)", 0);
Optionpk<bool> confusion_opt("cm", "confusion", "create confusion matrix (to std out) (default value is 0)", false);
Optionpk<string> labelref_opt("lr", "lref", "name of the reference label in case reference is shape file(default is label)", "label");
@@ -416,10 +416,11 @@ int main(int argc, char *argv[])
}
}
if(inputValue==referenceValue){//correct
- if(valueE_opt[0]!=nodata_opt[0])
- outputValue=valueE_opt[0];
- else
- outputValue=inputValue;
+ outputValue=valueE_opt[0];
+ if(nodata_opt.size()){
+ if(valueE_opt[0]==nodata_opt[0])
+ outputValue=inputValue;
+ }
}
else if(inputValue>referenceValue)//1=forest,2=non-forest
outputValue=valueO_opt[0];//omission error
@@ -473,10 +474,11 @@ int main(int argc, char *argv[])
}
}
if(inputValue==referenceValue){//correct
- if(valueE_opt[0]!=nodata_opt[0])
- outputValue=valueE_opt[0];
- else
- outputValue=inputValue;
+ outputValue=valueE_opt[0];
+ if(nodata_opt.size()){
+ if(valueE_opt[0]==nodata_opt[0])
+ outputValue=inputValue;
+ }
}
else if(inputValue>referenceValue)//1=forest,2=non-forest
outputValue=valueO_opt[0];//omission error
@@ -523,7 +525,8 @@ int main(int argc, char *argv[])
option_opt.push_back(theInterleave);
}
imgWriter.open(output_opt[0],inputReader.nrOfCol(),inputReader.nrOfRow(),1,inputReader.getDataType(),inputReader.getImageType(),option_opt);
- imgWriter.GDALSetNoDataValue(nodata_opt[0]);
+ if(nodata_opt.size())
+ imgWriter.GDALSetNoDataValue(nodata_opt[0]);
if(inputReader.isGeoRef()){
imgWriter.copyGeoTransform(inputReader);
}
@@ -653,10 +656,11 @@ int main(int argc, char *argv[])
}
if(lineInput[icol]==lineReference[ireference]){//correct
if(output_opt.size()){
- if(valueE_opt[0]!=nodata_opt[0])
- lineOutput[icol]=valueE_opt[0];
- else
- lineOutput[icol]=lineInput[icol];
+ lineOutput[icol]=valueE_opt[0];
+ if(nodata_opt.size()){
+ if(valueE_opt[0]==nodata_opt[0])
+ lineOutput[icol]=lineInput[icol];
+ }
}
}
else{//error
@@ -665,28 +669,32 @@ int main(int argc, char *argv[])
break;
}
if(output_opt.size()){
- if(lineInput[icol]<20){//forest
- if(lineReference[icol]>=20)//gain
- lineOutput[icol]=lineInput[icol]*10+1;//GAIN is 111,121,131
- else//forest type changed: mixed
- lineOutput[icol]=130;//MIXED FOREST
- }
- else if(lineReference[icol]<20){//loss
- lineOutput[icol]=20*10+lineReference[icol];//LOSS is 211 212 213
- }
- else//no forest
- lineOutput[icol]=20*10;//NON FOREST is 200
- // if(lineInput[icol]>lineReference[ireference])//1=forest,2=non-forest
- // lineOutput[icol]=valueO_opt[0];//omission error
- // else
- // lineOutput[icol]=valueC_opt[0];//commission error
+ // if(lineInput[icol]<20){//forest
+ // if(lineReference[icol]>=20)//gain
+ // lineOutput[icol]=lineInput[icol]*10+1;//GAIN is 111,121,131
+ // else//forest type changed: mixed
+ // lineOutput[icol]=130;//MIXED FOREST
+ // }
+ // else if(lineReference[icol]<20){//loss
+ // lineOutput[icol]=20*10+lineReference[icol];//LOSS is 211 212 213
+ // }
+ // else//no forest
+ // lineOutput[icol]=20*10;//NON FOREST is 200
+ if(lineInput[icol]>lineReference[ireference])
+ lineOutput[icol]=valueO_opt[0];//omission error
+ else
+ lineOutput[icol]=valueC_opt[0];//commission error
}
}
- }
+ }
else{
++nflagged;
- if(output_opt.size())
- lineOutput[icol]=nodata_opt[0];
+ if(output_opt.size()){
+ if(nodata_opt.size())
+ lineOutput[icol]=nodata_opt[0];
+ else //should never occur?
+ lineOutput[icol]=0;
+ }
}
}
if(output_opt.size()){
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 13995e3..fa1a94c 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -566,6 +566,9 @@ int main(int argc,char **argv) {
case(filter2d::dwtForward):
filter2d.dwtForward(input, output, wavelet_type_opt[0], family_opt[0]);
break;
+ case(filter2d::dwtInverse):
+ filter2d.dwtInverse(input, output, wavelet_type_opt[0], family_opt[0]);
+ break;
case(filter2d::dwtQuantize):
if(verbose_opt[0])
std::cout << "Quantization filtering" << std::endl;
diff --git a/src/apps/pkfilterascii.cc b/src/apps/pkfilterascii.cc
index 87ca91b..2052da2 100644
--- a/src/apps/pkfilterascii.cc
+++ b/src/apps/pkfilterascii.cc
@@ -24,12 +24,15 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include <math.h>
#include <sys/types.h>
#include <stdio.h>
-#include <gsl/gsl_sort.h>
#include "base/Optionpk.h"
#include "base/Vector2d.h"
#include "algorithms/Filter.h"
#include "fileclasses/FileReaderAscii.h"
+extern "C" {
+#include <gsl/gsl_sort.h>
+}
+
/*------------------
Main procedure
----------------*/
diff --git a/src/apps/pkinfo.cc b/src/apps/pkinfo.cc
index 87819f0..00823fc 100644
--- a/src/apps/pkinfo.cc
+++ b/src/apps/pkinfo.cc
@@ -279,7 +279,7 @@ int main(int argc, char *argv[])
if(hist_opt[0]){
assert(band_opt[0]<imgReader.nrOfBand());
int nbin=nbin_opt[0];
- std::vector<unsigned long int> output(nbin_opt[0]);
+ std::vector<unsigned long int> output;
minValue=0;
maxValue=0;
//todo: optimize such that getMinMax is only called once...
@@ -293,13 +293,21 @@ int main(int argc, char *argv[])
std::cout.precision(10);
for(int bin=0;bin<nbin;++bin){
// nsample+=output[bin];
- if(output[bin]>0){
+ if(nbin>1){
std::cout << (maxValue-minValue)*bin/(nbin-1)+minValue << " ";
if(relative_opt[0])
std::cout << 100.0*static_cast<double>(output[bin])/static_cast<double>(nsample) << std::endl;
else
std::cout << static_cast<double>(output[bin]) << std::endl;
}
+ else{
+ std::cout << (maxValue-minValue)*bin/(2-1)+minValue << " ";
+ if(relative_opt[0])
+ std::cout << 100.0*static_cast<double>(output[bin])/static_cast<double>(nsample) << std::endl;
+ else
+ std::cout << static_cast<double>(output[bin]) << std::endl;
+ }
+
}
}
else{
diff --git a/src/apps/pkmosaic.cc b/src/apps/pkmosaic.cc
index bb6e6c5..a5a904b 100644
--- a/src/apps/pkmosaic.cc
+++ b/src/apps/pkmosaic.cc
@@ -221,6 +221,10 @@ int main(int argc, char *argv[])
}
double theULX, theULY, theLRX, theLRY;
imgReader.getBoundingBox(theULX,theULY,theLRX,theLRY);
+ if(theLRY>theULY){
+ cerr << "Error: " << input_opt[ifile] << " is not georeferenced, only referenced images are supported for pkmosaic " << endl;
+ exit(1);
+ }
if(verbose_opt[0])
cout << "Bounding Box (ULX ULY LRX LRY): " << fixed << setprecision(6) << theULX << " " << theULY << " " << theLRX << " " << theLRY << endl;
if(!init){
diff --git a/src/imageclasses/ImgReaderGdal.cc b/src/imageclasses/ImgReaderGdal.cc
index e0579d3..963012d 100644
--- a/src/imageclasses/ImgReaderGdal.cc
+++ b/src/imageclasses/ImgReaderGdal.cc
@@ -217,36 +217,58 @@ std::string ImgReaderGdal::getCompression() const
bool ImgReaderGdal::getBoundingBox(double& ulx, double& uly, double& lrx, double& lry) const
{
+ double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
+ m_gds->GetGeoTransform(gt);
+
+ //assuming
+ //adfGeotransform[0]: ULX (upper left X coordinate)
+ //adfGeotransform[1]: $cos(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[2]: $-sin(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[3]: ULY (upper left Y coordinate)
+ //adfGeotransform[4]: $-sin(\alpha)\cdot\textrm{Yres}$
+ //adfGeotransform[5]: $-cos(\alpha)\cdot\textrm{Yres}$
+ ulx=gt[0];
+ uly=gt[3];
+ lrx=gt[0]+nrOfCol()*gt[1]+nrOfRow()*gt[2];
+ lry=gt[3]+nrOfCol()*gt[4]+nrOfRow()*gt[5];
if(m_isGeoRef){
- // ulx=m_ulx-(m_magic_x-1.0)*m_delta_x;
- // uly=m_uly+(m_magic_y-1.0)*m_delta_y;
- ulx=m_ulx;
- uly=m_uly;
- lrx=ulx+nrOfCol()*m_delta_x;
- lry=uly-nrOfRow()*m_delta_y;
+ // ulx=m_ulx;
+ // uly=m_uly;
+ // lrx=ulx+nrOfCol()*m_delta_x;
+ // lry=uly-nrOfRow()*m_delta_y;
return true;
}
else{
- ulx=0;
- uly=nrOfRow()-1;
- lrx=nrOfCol()-1;
- lry=0;
+ // ulx=0;
+ // uly=nrOfRow()-1;
+ // lrx=nrOfCol()-1;
+ // lry=0;
return false;
}
}
bool ImgReaderGdal::getCentrePos(double& x, double& y) const
{
+ double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
+ m_gds->GetGeoTransform(gt);
+
+ //assuming
+ //adfGeotransform[0]: ULX (upper left X coordinate)
+ //adfGeotransform[1]: $cos(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[2]: $-sin(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[3]: ULY (upper left Y coordinate)
+ //adfGeotransform[4]: $-sin(\alpha)\cdot\textrm{Yres}$
+ //adfGeotransform[5]: $-cos(\alpha)\cdot\textrm{Yres}$
+ x=gt[0]+(nrOfCol()/2.0)*gt[1]+(nrOfRow()/2.0)*gt[2];
+ y=gt[3]+(nrOfCol()/2.0)*gt[4]+(nrOfRow()/2.0)*gt[5];
if(m_isGeoRef){
-// x=m_ulx+(nrOfCol()/2.0-(m_magic_x-1.0))*m_delta_x;
-// y=m_uly-(nrOfRow()/2.0-(m_magic_y-1.0))*m_delta_y;
- x=m_ulx+(nrOfCol()/2.0)*m_delta_x;
- y=m_uly-(nrOfRow()/2.0)*m_delta_y;
+ // x=m_ulx+(nrOfCol()/2.0)*m_delta_x;
+ // y=m_uly-(nrOfRow()/2.0)*m_delta_y;
return true;
}
else{
- x=nrOfCol()/2.0;
- y=nrOfRow()/2.0;
+ // x=nrOfCol()/2.0;
+ // y=nrOfRow()/2.0;
return false;
}
}
@@ -255,18 +277,31 @@ bool ImgReaderGdal::getCentrePos(double& x, double& y) const
bool ImgReaderGdal::geo2image(double x, double y, double& i, double& j) const
{
//double values are returned, caller is responsible for interpolation step
+ double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
+ m_gds->GetGeoTransform(gt);
+ //assuming
+ //adfGeotransform[0]: ULX (upper left X coordinate)
+ //adfGeotransform[1]: $cos(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[2]: $-sin(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[3]: ULY (upper left Y coordinate)
+ //adfGeotransform[4]: $-sin(\alpha)\cdot\textrm{Yres}$
+ //adfGeotransform[5]: $-cos(\alpha)\cdot\textrm{Yres}$
+ double denom=(gt[1]-gt[2]*gt[4]/gt[5]);
+ double eps=0.00001;
+ if(denom>eps){
+ i=(x-gt[0]-gt[2]/gt[5]*(y-gt[3]))/denom;
+ j=(y-gt[3]-gt[4]*(x-gt[0]-gt[2]/gt[5]*(y-gt[3]))/denom)/gt[5];
+ }
if(m_isGeoRef){
-// double ulx=m_ulx-(m_magic_x-1.0)*m_delta_x;
-// double uly=m_uly+(m_magic_y-1.0)*m_delta_y;
- double ulx=m_ulx;
- double uly=m_uly;
- i=(x-ulx)/m_delta_x;
- j=(uly-y)/m_delta_y;
+ // double ulx=m_ulx;
+ // double uly=m_uly;
+ // i=(x-ulx)/m_delta_x;
+ // j=(uly-y)/m_delta_y;
return true;
}
else{
- i=x;
- j=nrOfRow()-y;
+ // i=x;
+ // j=nrOfRow()-y;
return false;
}
}
@@ -274,16 +309,27 @@ bool ImgReaderGdal::geo2image(double x, double y, double& i, double& j) const
//x and y represent centre of pixel, return true if image is georeferenced
bool ImgReaderGdal::image2geo(double i, double j, double& x, double& y) const
{
+ double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
+ m_gds->GetGeoTransform(gt);
+
+ //assuming
+ //adfGeotransform[0]: ULX (upper left X coordinate)
+ //adfGeotransform[1]: $cos(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[2]: $-sin(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[3]: ULY (upper left Y coordinate)
+ //adfGeotransform[4]: $-sin(\alpha)\cdot\textrm{Yres}$
+ //adfGeotransform[5]: $-cos(\alpha)\cdot\textrm{Yres}$
+
+ x=gt[0]+(0.5+i)*gt[1]+(0.5+j)*gt[2];
+ y=gt[3]+(0.5+i)*gt[4]+(0.5+j)*gt[5];
if(m_isGeoRef){
-// x=m_ulx+(1.5-m_magic_x+i)*m_delta_x;
-// y=m_uly-(1.5-m_magic_y+j)*m_delta_y;
- x=m_ulx+(0.5+i)*m_delta_x;
- y=m_uly-(0.5+j)*m_delta_y;
+ // x=m_ulx+(0.5+i)*m_delta_x;
+ // y=m_uly-(0.5+j)*m_delta_y;
return true;
}
else{
- x=0.5+i;
- y=nrOfRow()-(0.5+j);
+ // x=0.5+i;
+ // y=nrOfRow()-(0.5+j);
return false;
}
}
@@ -455,15 +501,20 @@ unsigned long int ImgReaderGdal::getHistogram(std::vector<unsigned long int>& hi
maxValue=max;
min=minValue;
max=maxValue;
- if(nbin==0)
- nbin=maxValue-minValue+1;
+ double scale=0;
+ if(maxValue>minValue){
+ if(nbin==0)
+ nbin=maxValue-minValue+1;
+ scale=static_cast<double>(nbin-1)/(maxValue-minValue);
+ }
+ else
+ nbin=1;
assert(nbin>0);
histvector.resize(nbin);
unsigned long int nsample=0;
unsigned long int ninvalid=0;
std::vector<double> lineBuffer(nrOfCol());
for(int i=0;i<nbin;histvector[i++]=0);
- double scale=static_cast<double>(nbin-1)/(maxValue-minValue);
for(int irow=0;irow<nrOfRow();++irow){
readData(lineBuffer,GDT_Float64,irow,theBand);
for(int icol=0;icol<nrOfCol();++icol){
@@ -473,6 +524,8 @@ unsigned long int ImgReaderGdal::getHistogram(std::vector<unsigned long int>& hi
++ninvalid;
else if(lineBuffer[icol]<minValue)
++ninvalid;
+ else if(nbin==1)
+ ++histvector[0];
else{//scale to [0:nbin]
int theBin=static_cast<unsigned long int>(scale*(lineBuffer[icol]-minValue));
assert(theBin>=0);
diff --git a/src/imageclasses/ImgWriterGdal.cc b/src/imageclasses/ImgWriterGdal.cc
index a640268..0b8da57 100644
--- a/src/imageclasses/ImgWriterGdal.cc
+++ b/src/imageclasses/ImgWriterGdal.cc
@@ -355,33 +355,54 @@ string ImgWriterGdal::setProjection(void)
bool ImgWriterGdal::getBoundingBox(double& ulx, double& uly, double& lrx, double& lry) const
{
+ double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
+ m_gds->GetGeoTransform(gt);
+
+ //assuming
+ //adfGeotransform[0]: ULX (upper left X coordinate)
+ //adfGeotransform[1]: $cos(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[2]: $-sin(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[3]: ULY (upper left Y coordinate)
+ //adfGeotransform[4]: $-sin(\alpha)\cdot\textrm{Yres}$
+ //adfGeotransform[5]: $-cos(\alpha)\cdot\textrm{Yres}$
+
+ ulx=gt[0];
+ uly=gt[3];
+ lrx=gt[0]+nrOfCol()*gt[1]+nrOfRow()*gt[2];
+ lry=gt[3]+nrOfCol()*gt[4]+nrOfRow()*gt[5];
if(m_isGeoRef){
-// ulx=m_ulx-(m_magic_x-1.0)*m_delta_x;
-// uly=m_uly+(m_magic_y-1.0)*m_delta_y;
-// lrx=ulx+(nrOfCol()+1.0-m_magic_x)*m_delta_x;
-// lry=uly-(nrOfRow()+1.0-m_magic_y)*m_delta_y;
- ulx=m_ulx;
- uly=m_uly;
- lrx=ulx+nrOfCol()*m_delta_x;
- lry=uly-nrOfRow()*m_delta_y;
+ // ulx=m_ulx;
+ // uly=m_uly;
+ // lrx=ulx+nrOfCol()*m_delta_x;
+ // lry=uly-nrOfRow()*m_delta_y;
return true;
}
else{
- ulx=0;
- uly=nrOfRow()-1;
- lrx=nrOfCol()-1;
- lry=0;
+ // ulx=0;
+ // uly=nrOfRow()-1;
+ // lrx=nrOfCol()-1;
+ // lry=0;
return false;
}
}
bool ImgWriterGdal::getCentrePos(double& x, double& y) const
{
+ double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
+ m_gds->GetGeoTransform(gt);
+
+ //assuming
+ //adfGeotransform[0]: ULX (upper left X coordinate)
+ //adfGeotransform[1]: $cos(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[2]: $-sin(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[3]: ULY (upper left Y coordinate)
+ //adfGeotransform[4]: $-sin(\alpha)\cdot\textrm{Yres}$
+ //adfGeotransform[5]: $-cos(\alpha)\cdot\textrm{Yres}$
+ x=gt[0]+(nrOfCol()/2.0)*gt[1]+(nrOfRow()/2.0)*gt[2];
+ y=gt[3]+(nrOfCol()/2.0)*gt[4]+(nrOfRow()/2.0)*gt[5];
if(m_isGeoRef){
-// x=m_ulx+(nrOfCol()/2.0-(m_magic_x-1.0))*m_delta_x;
-// y=m_uly-(nrOfRow()/2.0+(m_magic_y-1.0))*m_delta_y;
- x=m_ulx+nrOfCol()/2.0*m_delta_x;
- y=m_uly-nrOfRow()/2.0*m_delta_y;
+ // x=m_ulx+nrOfCol()/2.0*m_delta_x;
+ // y=m_uly-nrOfRow()/2.0*m_delta_y;
return true;
}
else
@@ -391,18 +412,32 @@ bool ImgWriterGdal::getCentrePos(double& x, double& y) const
bool ImgWriterGdal::geo2image(double x, double y, double& i, double& j) const
{
//double values are returned, caller is responsible for interpolation step
+ //double values are returned, caller is responsible for interpolation step
+ double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
+ m_gds->GetGeoTransform(gt);
+ //assuming
+ //adfGeotransform[0]: ULX (upper left X coordinate)
+ //adfGeotransform[1]: $cos(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[2]: $-sin(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[3]: ULY (upper left Y coordinate)
+ //adfGeotransform[4]: $-sin(\alpha)\cdot\textrm{Yres}$
+ //adfGeotransform[5]: $-cos(\alpha)\cdot\textrm{Yres}$
+ double denom=(gt[1]-gt[2]*gt[4]/gt[5]);
+ double eps=0.00001;
+ if(denom>eps){
+ i=(x-gt[0]-gt[2]/gt[5]*(y-gt[3]))/denom;
+ j=(y-gt[3]-gt[4]*(x-gt[0]-gt[2]/gt[5]*(y-gt[3]))/denom)/gt[5];
+ }
if(m_isGeoRef){
-// double ulx=m_ulx-(m_magic_x-1.0)*m_delta_x;
-// double uly=m_uly+(m_magic_y-1.0)*m_delta_y;
- double ulx=m_ulx;
- double uly=m_uly;
- i=(x-ulx)/m_delta_x;
- j=(uly-y)/m_delta_y;
+ // double ulx=m_ulx;
+ // double uly=m_uly;
+ // i=(x-ulx)/m_delta_x;
+ // j=(uly-y)/m_delta_y;
return true;
}
else{
- i=x;
- j=nrOfRow()-y;
+ // i=x;
+ // j=nrOfRow()-y;
return false;
}
}
@@ -410,11 +445,23 @@ bool ImgWriterGdal::geo2image(double x, double y, double& i, double& j) const
//centre of pixel is always returned (regardless of magic pixel reference)!
bool ImgWriterGdal::image2geo(double i, double j, double& x, double& y) const
{
+ double gt[6];// { 444720, 30, 0, 3751320, 0, -30 };
+ m_gds->GetGeoTransform(gt);
+
+ //assuming
+ //adfGeotransform[0]: ULX (upper left X coordinate)
+ //adfGeotransform[1]: $cos(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[2]: $-sin(\alpha)\cdot\textrm{Xres}$
+ //adfGeotransform[3]: ULY (upper left Y coordinate)
+ //adfGeotransform[4]: $-sin(\alpha)\cdot\textrm{Yres}$
+ //adfGeotransform[5]: $-cos(\alpha)\cdot\textrm{Yres}$
+
+ x=gt[0]+(0.5+i)*gt[1]+(0.5+j)*gt[2];
+ y=gt[3]+(0.5+i)*gt[4]+(0.5+j)*gt[5];
+
if(m_isGeoRef){
-// x=m_ulx+(1.5-m_magic_x+i)*m_delta_x;
-// y=m_uly-(1.5-m_magic_y+j)*m_delta_y;
- x=m_ulx+(0.5+i)*m_delta_x;
- y=m_uly-(0.5+j)*m_delta_y;
+ // x=m_ulx+(0.5+i)*m_delta_x;
+ // y=m_uly-(0.5+j)*m_delta_y;
return true;
}
else
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pktools.git
More information about the Pkg-grass-devel
mailing list