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 <kempe...@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 _______________________________________________ Pkg-grass-devel mailing list Pkg-grass-devel@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel