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 f659640027f027ed00dec14753180ed478beb0e5
Author: Pieter Kempeneers <kempe...@gmail.com>
Date: Tue Jul 9 10:43:46 2013 +0200
added support for discrete wavelet transform in filters
---
src/algorithms/Filter.cc | 43 +++++++++-------
src/algorithms/Filter.h | 25 ++++-----
src/algorithms/Filter2d.cc | 19 +++++++
src/algorithms/Filter2d.h | 106 +++++++++++++++++++++++++++++++++++++-
src/apps/Makefile.am | 4 +-
src/apps/pkcrop.cc | 6 +--
src/apps/pkextract.cc | 50 +++++++++---------
src/apps/pkfilter.cc | 14 +++++
src/apps/pkfilterascii.cc | 61 +++++++++++++++-------
src/apps/pkinfo.cc | 70 +++++++++++++------------
src/imageclasses/ImgReaderGdal.cc | 39 ++++++++++++++
src/imageclasses/ImgReaderGdal.h | 1 +
src/imageclasses/Makefile.am | 4 +-
13 files changed, 325 insertions(+), 117 deletions(-)
diff --git a/src/algorithms/Filter.cc b/src/algorithms/Filter.cc
index fc7cb2d..56d3bf0 100644
--- a/src/algorithms/Filter.cc
+++ b/src/algorithms/Filter.cc
@@ -19,6 +19,7 @@ along with pktools. If not, see
<http://www.gnu.org/licenses/>.
***********************************************************************/
#include "Filter.h"
#include <assert.h>
+#include <math.h>
#include <iostream>
filter::Filter::Filter(void)
@@ -38,25 +39,31 @@ void filter::Filter::setTaps(const vector<double> &taps)
assert(m_taps.size()%2);
}
-// void filter::Filter::dwtForward(std::vector<double>& data, const
std::string& wavelet_type, int family){
-// int nsize=data.size();
-// assert(nsize);
-// gsl_wavelet *w;
-// gsl_wavelet_workspace *work;
-// w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
-// work=gsl_wavelet_workspace_alloc(nsize);
-// gsl_wavelet_transform_forward(w,&(data[0]),1,nsize,work);
-// }
+void filter::Filter::dwtForward(std::vector<double>& data, const std::string&
wavelet_type, int family){
+ //make sure data size if power of 2
+ while(data.size()&(data.size()-1))
+ data.push_back(data.back());
+ int nsize=data.size();
+ gsl_wavelet *w;
+ gsl_wavelet_workspace *work;
+ assert(nsize);
+ w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
+ work=gsl_wavelet_workspace_alloc(nsize);
+ gsl_wavelet_transform_forward(w,&(data[0]),1,nsize,work);
+}
-// void filter::Filter::dwtInverse(std::vector<double>& data, const
std::string& wavelet_type, int family){
-// int nsize=data.size();
-// assert(nsize);
-// gsl_wavelet *w;
-// gsl_wavelet_workspace *work;
-// w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
-// work=gsl_wavelet_workspace_alloc(nsize);
-// gsl_wavelet_transform_inverse(w,&(data[0]),1,nsize,work);
-// }
+void filter::Filter::dwtInverse(std::vector<double>& data, const std::string&
wavelet_type, int family){
+ //make sure data size if power of 2
+ while(data.size()&(data.size()-1))
+ data.push_back(data.back());
+ int nsize=data.size();
+ assert(nsize);
+ gsl_wavelet *w;
+ gsl_wavelet_workspace *work;
+ w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
+ work=gsl_wavelet_workspace_alloc(nsize);
+ gsl_wavelet_transform_inverse(w,&(data[0]),1,nsize,work);
+}
void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal&
output, const std::string& method, int dim, short down, int offset)
{
diff --git a/src/algorithms/Filter.h b/src/algorithms/Filter.h
index ef89a01..1b6d30e 100644
--- a/src/algorithms/Filter.h
+++ b/src/algorithms/Filter.h
@@ -31,7 +31,7 @@ using namespace std;
namespace filter
{
- enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6,
dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13,
sobelxy=14, sobelyx=-14, smooth=15, density=16, majority=17, mixed=18,
smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24,
stdev=25, dwtForward=26, dwtInverse=27};
+ enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6,
dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13,
sobelxy=14, sobelyx=-14, smooth=15, density=16, majority=17, mixed=18,
smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24,
stdev=25, dwtForward=26, dwtInverse=27, dwtQuantize=28};
class Filter
{
@@ -40,9 +40,12 @@ public:
Filter(const vector<double> &taps);
virtual ~Filter(){};
static const gsl_wavelet_type* getWaveletType(const std::string waveletType){
- std::map<std::string, const gsl_wavelet_type* > m_waveletMap;
- initWaveletMap(m_waveletMap);
- return m_waveletMap[waveletType];
+ if(waveletType=="daubechies") return(gsl_wavelet_daubechies);
+ if(waveletType=="daubechies_centered")
return(gsl_wavelet_daubechies_centered);
+ if(waveletType=="haar") return(gsl_wavelet_haar);
+ if(waveletType=="haar_centered") return(gsl_wavelet_haar_centered);
+ if(waveletType=="bspline") return(gsl_wavelet_bspline);
+ if(waveletType=="bspline_centered") return(gsl_wavelet_bspline_centered);
}
static FILTER_TYPE getFilterType(const std::string filterType){
std::map<std::string, FILTER_TYPE> m_filterMap;
@@ -63,24 +66,16 @@ public:
template<class T> void applyFwhm(const vector<double> &wavelengthIn, const
vector<T>& input, const vector<double> &wavelengthOut, const vector<double>
&fwhm, const std::string& interpolationType, vector<T>& output, bool
verbose=false);
template<class T> void applyFwhm(const vector<double> &wavelengthIn, const
Vector2d<T>& input, const vector<double> &wavelengthOut, const vector<double>
&fwhm, const std::string& interpolationType, Vector2d<T>& output, int down=1,
bool verbose=false);
- // void dwtForward(std::vector<double>& data, const std::string&
wavelet_type, int family);
- // void dwtInverse(std::vector<double>& data, const std::string&
wavelet_type, int family);
+ void dwtForward(std::vector<double>& data, const std::string& wavelet_type,
int family);
+ void dwtInverse(std::vector<double>& data, const std::string& wavelet_type,
int family);
private:
- static void initWaveletMap(std::map<std::string, const gsl_wavelet_type*>
m_waveletMap){
- //initialize Map
- m_waveletMap["daubechies"]=gsl_wavelet_daubechies;
- m_waveletMap["daubechies_centered"]=gsl_wavelet_daubechies_centered;
- m_waveletMap["haar"]=gsl_wavelet_haar;
- m_waveletMap["haar_centered"]=gsl_wavelet_haar_centered;
- m_waveletMap["bspline"]=gsl_wavelet_bspline;
- m_waveletMap["bspline_centered"]=gsl_wavelet_bspline_centered;
- }
static void initFilterMap(std::map<std::string, FILTER_TYPE>& m_filterMap){
//initialize Map
m_filterMap["dwtForward"]=filter::dwtForward;
m_filterMap["dwtInverse"]=filter::dwtInverse;
+ m_filterMap["dwtQuantize"]=filter::dwtQuantize;
m_filterMap["stdev"]=filter::stdev;
m_filterMap["var"]=filter::var;
m_filterMap["min"]=filter::min;
diff --git a/src/algorithms/Filter2d.cc b/src/algorithms/Filter2d.cc
index 11a2d88..3cb5e73 100644
--- a/src/algorithms/Filter2d.cc
+++ b/src/algorithms/Filter2d.cc
@@ -1106,3 +1106,22 @@ void filter2d::Filter2d::shadowDsm(const ImgReaderGdal&
input, ImgWriterGdal& ou
shadowDsm(inputBuffer, outputBuffer, sza, saa, pixelSize, shadowFlag);
output.writeDataBlock(outputBuffer,GDT_Float32,0,output.nrOfCol()-1,0,output.nrOfRow()-1,0);
}
+
+void filter2d::Filter2d::dwtForward(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);
+ dwtForward(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){
+ input.readDataBlock(theBuffer, GDT_Float32, 0, input.nrOfCol()-1, 0,
input.nrOfRow()-1, iband);
+ std::cout << "filtering band " << iband << std::endl << std::flush;
+ dwtQuantize(theBuffer, wavelet_type, family, quantize);
+
output.writeDataBlock(theBuffer,GDT_Float32,0,output.nrOfCol()-1,0,output.nrOfRow()-1,iband);
+ }
+}
diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h
index fb2d3ea..b136eac 100644
--- a/src/algorithms/Filter2d.h
+++ b/src/algorithms/Filter2d.h
@@ -34,7 +34,11 @@ along with pktools. If not, see
<http://www.gnu.org/licenses/>.
#include <vector>
#include <string>
#include <map>
+#include <gsl/gsl_sort.h>
+#include <gsl/gsl_wavelet.h>
+#include <gsl/gsl_wavelet2d.h>
#include "base/Vector2d.h"
+#include "Filter.h"
#include "imageclasses/ImgReaderGdal.h"
#include "imageclasses/ImgWriterGdal.h"
#include "algorithms/StatFactory.h"
@@ -43,7 +47,7 @@ using namespace std;
// using namespace cimg_library;
namespace filter2d
{
- enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6,
dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13,
sobelxy=14, sobelyx=-14, smooth=15, density=16, majority=17, mixed=18,
smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24,
stdev=25, mrf=26};
+ enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6,
dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13,
sobelxy=14, sobelyx=-14, smooth=15, density=16, majority=17, mixed=18,
smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24,
stdev=25, mrf=26, dwtForward=27, dwtInverse=28, dwtQuantize=29};
class Filter2d
{
@@ -72,6 +76,11 @@ public:
template<class T1, class T2> void filter(const Vector2d<T1>& inputVector,
Vector2d<T2>& outputVector);
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 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);
void majorVoting(const string& inputFilename, const string&
outputFilename,int dim=0,const vector<int> &prior=vector<int>());
/* void homogeneousSpatial(const string& inputFilename, const 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);
@@ -118,6 +127,7 @@ private:
m_filterMap["heterog"]=filter2d::heterog;
m_filterMap["order"]=filter2d::order;
m_filterMap["median"]=filter2d::median;
+ m_filterMap["dwtQuantize"]=filter2d::dwtQuantize;
}
Vector2d<double> m_taps;
@@ -628,6 +638,100 @@ 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){
+ //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_forward (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];
+ }
+ }
+}
+
+template<class T> void Filter2d::dwtQuantize(Vector2d<T>& theBuffer, const
std::string& wavelet_type, int family, double quantize){
+ //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=new double[theBuffer.size()*theBuffer[0].size()];
+ double* abscoeff=new double[theBuffer.size()*theBuffer[0].size()];
+ size_t* p=new size_t[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;
+ assert(index<theBuffer.size()*theBuffer[0].size());
+ 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_forward (w, data, theBuffer.size(),
theBuffer[0].size(),theBuffer[0].size(), work);
+ for(int irow=0;irow<theBuffer.size();++irow){
+ for(int icol=0;icol<theBuffer[0].size();++icol){
+ int index=irow*theBuffer[0].size()+icol;
+ abscoeff[index]=fabs(data[index]);
+ if(quantize<0){//absolute threshold
+ if(abscoeff[index]<-quantize)
+ data[index]=0;
+ }
+ }
+ }
+ if(quantize>0){//percentual threshold
+ int nc=quantize/100.0*nsize;
+ gsl_sort_index(p,abscoeff,1,nsize);
+ for(int i=0;(i+nc)<nsize;i++)
+ data[p[i]]=0;
+ }
+ gsl_wavelet2d_nstransform_inverse (w, data, theBuffer.size(),
theBuffer[0].size(),theBuffer[0].size(), work);
+ for(int irow=0;irow<theBuffer.size();++irow){
+ for(int icol=0;icol<theBuffer[irow].size();++icol){
+ int index=irow*theBuffer[irow].size()+icol;
+ theBuffer[irow][icol]=data[index];
+ }
+ }
+ 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());
+ delete[] data;
+ delete[] abscoeff;
+ delete[] p;
+}
+
}
#endif /* _MYFILTER_H_ */
diff --git a/src/apps/Makefile.am b/src/apps/Makefile.am
index 6eeebb2..2de9f1a 100644
--- a/src/apps/Makefile.am
+++ b/src/apps/Makefile.am
@@ -6,7 +6,7 @@ LDADD = $(GDAL_LDFLAGS)
$(top_builddir)/src/algorithms/libalgorithms.la $(top_bu
###############################################################################
# the program to build and install (the names of the final binaries)
-bin_PROGRAMS = pkinfo pkcrop pkreclass pkgetmask pksetmask pkcreatect
pkdumpimg pkdumpogr pksieve pkstat pkstatogr pkegcs pkextract pkfillnodata
pkfilter pkfilterascii pkdsm2shadow pkmosaic pkndvi pkpolygonize pkascii2img
pkdiff pkclassify_svm pkfs_svm pkascii2ogr pkeditogr
+bin_PROGRAMS = pkinfo pkcrop pkreclass pkgetmask pksetmask pkcreatect
pkdumpimg pkdumpogr pksieve pkstat pkstatogr pkegcs pkextract pkfillnodata
pkfilter pkenhance pkfilterascii pkdsm2shadow pkmosaic pkndvi pkpolygonize
pkascii2img pkdiff pkclassify_svm pkfs_svm pkascii2ogr pkeditogr
# the program to build but not install (the names of the final binaries)
#noinst_PROGRAMS = pkxcorimg pkgeom
@@ -53,6 +53,8 @@ pkextract_SOURCES = pkextract.cc
pkfillnodata_SOURCES = pkfillnodata.cc
pkfilter_SOURCES = pkfilter.cc
pkfilter_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lgsl -lgdal
+pkenhance_SOURCES = pkenhance.cc
+pkenhance_LDADD = $(AM_LDFLAGS) -lgdal
pkfilterascii_SOURCES = pkfilterascii.cc
pkfilterascii_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lgsl -lgdal
pkdsm2shadow_SOURCES = pkdsm2shadow.cc
diff --git a/src/apps/pkcrop.cc b/src/apps/pkcrop.cc
index ea65ca2..6f86cb2 100644
--- a/src/apps/pkcrop.cc
+++ b/src/apps/pkcrop.cc
@@ -53,8 +53,8 @@ int main(int argc, char *argv[])
Optionpk<double> offset_opt("off", "offset", "output=scale*input+offset", 0);
Optionpk<string> otype_opt("ot", "otype", "Data type for output image
({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}).
Empty string: inherit type from input image","");
Optionpk<string> oformat_opt("of", "oformat", "Output image format (see
also gdal_translate). Empty string: inherit from input image");
- Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5
columns: id R G B ALFA (0: transparent, 255: solid)");
Optionpk<string> option_opt("co", "co", "options: NAME=VALUE [-co
COMPRESS=LZW] [-co INTERLEAVE=BAND]");
+ Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5
columns: id R G B ALFA (0: transparent, 255: solid)");
Optionpk<short> flag_opt("f", "flag", "Flag value to put in image if out of
bounds.", 0);
Optionpk<string> resample_opt("r", "resampling-method", "Resampling method
(near: nearest neighbour, bilinear: bi-linear interpolation).", "near");
Optionpk<string> description_opt("d", "description", "Set image
description", "");
@@ -77,6 +77,7 @@ int main(int argc, char *argv[])
offset_opt.retrieveOption(argc,argv);
otype_opt.retrieveOption(argc,argv);
oformat_opt.retrieveOption(argc,argv);
+ option_opt.retrieveOption(argc,argv);
colorTable_opt.retrieveOption(argc,argv);
dx_opt.retrieveOption(argc,argv);
dy_opt.retrieveOption(argc,argv);
@@ -86,7 +87,6 @@ int main(int argc, char *argv[])
ny_opt.retrieveOption(argc,argv);
ns_opt.retrieveOption(argc,argv);
nl_opt.retrieveOption(argc,argv);
- option_opt.retrieveOption(argc,argv);
flag_opt.retrieveOption(argc,argv);
resample_opt.retrieveOption(argc,argv);
description_opt.retrieveOption(argc,argv);
@@ -105,7 +105,7 @@ int main(int argc, char *argv[])
exit(0);//help was invoked, stop processing
}
if(output_opt.empty()){
- std::cerr << "No output file provided (use option -i). Use pkinfo --help
for help information" << std::endl;
+ std::cerr << "No output file provided (use option -o). Use pkinfo --help
for help information" << std::endl;
exit(0);//help was invoked, stop processing
}
diff --git a/src/apps/pkextract.cc b/src/apps/pkextract.cc
index 007c504..90b3028 100644
--- a/src/apps/pkextract.cc
+++ b/src/apps/pkextract.cc
@@ -125,7 +125,7 @@ int main(int argc, char *argv[])
}
short theDim=boundary_opt[0];
if(verbose_opt[0]>1)
- std::cout << boundary_opt[0] << std::endl;
+ std::cout << "boundary: " << boundary_opt[0] << std::endl;
ImgReaderGdal imgReader;
try{
imgReader.open(image_opt[0]);
@@ -398,18 +398,20 @@ int main(int argc, char *argv[])
if(p>theThreshold)
continue;//do not select for now, go to next column
}
- else{//absolute value
+ else if(nvalid.size()>processClass){//absolute value
if(nvalid[processClass]>-theThreshold)
- continue;//do not select any more pixels for this class, go
to next column to search for other classes
+ continue;//do not select any more pixels for this class, go
to next column to search for other classes
}
writeBuffer.push_back(sample);
writeBufferClass.push_back(theClass);
++ntotalvalid;
- ++(nvalid[processClass]);
+ if(nvalid.size()>processClass)
+ ++(nvalid[processClass]);
}
else{
++ntotalinvalid;
- ++(ninvalid[processClass]);
+ if(ninvalid.size()>processClass)
+ ++(ninvalid[processClass]);
}
}//processClass
}//icol
@@ -468,8 +470,10 @@ int main(int argc, char *argv[])
nsample=writeBuffer.size();
if(verbose_opt[0]){
std::cout << "total number of samples written: " << nsample <<
std::endl;
- for(int iclass=0;iclass<class_opt.size();++iclass)
- std::cout << "class " << class_opt[iclass] << " has " <<
nvalid[iclass] << " samples" << std::endl;
+ if(nvalid.size()==class_opt.size()){
+ for(int iclass=0;iclass<class_opt.size();++iclass)
+ std::cout << "class " << class_opt[iclass] << " has " <<
nvalid[iclass] << " samples" << std::endl;
+ }
}
}
else{//class_opt[0]!=0
@@ -518,7 +522,7 @@ int main(int argc, char *argv[])
int theClass=0;
// double theClass=0;
int processClass=-1;
- if(class_opt[0]<0){//process every class except 0
+ if(class_opt.empty()){//process every class
if(classBuffer[icol]){
processClass=0;
theClass=classBuffer[icol];
@@ -643,7 +647,7 @@ int main(int argc, char *argv[])
if(p>theThreshold)
continue;//do not select for now, go to next column
}
- else{//absolute value
+ else if(nvalid.size()>processClass){//absolute value
if(nvalid[processClass]>-theThreshold)
continue;//do not select any more pixels for this class, go
to next column to search for other classes
}
@@ -651,11 +655,13 @@ int main(int argc, char *argv[])
//
writeBufferClass.push_back(class_opt[processClass]);
writeBufferClass.push_back(theClass);
++ntotalvalid;
- ++(nvalid[processClass]);
+ if(nvalid.size()>processClass)
+ ++(nvalid[processClass]);
}
else{
++ntotalinvalid;
- ++(ninvalid[processClass]);
+ if(ninvalid.size()>processClass)
+ ++(ninvalid[processClass]);
}
}//processClass
}//icol
@@ -716,8 +722,10 @@ int main(int argc, char *argv[])
nsample=writeBuffer.size();
if(verbose_opt[0]){
std::cout << "total number of samples written: " << nsample <<
std::endl;
- for(int iclass=0;iclass<class_opt.size();++iclass)
- std::cout << "class " << class_opt[iclass] << " has " <<
nvalid[iclass] << " samples" << std::endl;
+ if(nvalid.size()==class_opt.size()){
+ for(int iclass=0;iclass<class_opt.size();++iclass)
+ std::cout << "class " << class_opt[iclass] << " has " <<
nvalid[iclass] << " samples" << std::endl;
+ }
}
}
}
@@ -768,6 +776,7 @@ int main(int argc, char *argv[])
switch(ruleMap[rule_opt[0]]){
// switch(rule_opt[0]){
case(rule::proportion):{//proportion for each class
+ assert(class_opt.size());
for(int iclass=0;iclass<class_opt.size();++iclass){
ostringstream cs;
cs << class_opt[iclass];
@@ -781,9 +790,10 @@ int main(int argc, char *argv[])
case(rule::minimum):
case(rule::maximum):
case(rule::maxvote):
+ assert(class_opt.size());
ogrWriter.createField(label_opt[0],fieldType);
- if(test_opt.size())
- ogrTestWriter.createField(label_opt[0],fieldType);
+ if(test_opt.size())
+ ogrTestWriter.createField(label_opt[0],fieldType);
break;
case(rule::point):
case(rule::mean):
@@ -1218,9 +1228,6 @@ int main(int argc, char *argv[])
imgReader.image2geo(i,j,x,y);
thePoint.setX(x);
thePoint.setY(y);
- // //test
- // writeRing.addPoint(&thePoint);
- // if(thePoint.Within(&readPolygon)){
if(readPolygon.Contains(&thePoint)){
bool valid=true;
for(int imask=0;imask<mask_opt.size();++imask){
@@ -1432,8 +1439,6 @@ int main(int argc, char *argv[])
}
}
}
- // //test
- // std::cout << "before write polygon" << std::endl;
if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean){
//do not create if no points found within polygon
if(!nPointPolygon)
@@ -1832,9 +1837,6 @@ int main(int argc, char *argv[])
imgReader.image2geo(i,j,x,y);
thePoint.setX(x);
thePoint.setY(y);
- // //test
- // writeRing.addPoint(&thePoint);
- // if(thePoint.Within(&readPolygon)){
if(readPolygon.Contains(&thePoint)){
bool valid=true;
for(int imask=0;imask<mask_opt.size();++imask){
@@ -2044,8 +2046,6 @@ int main(int argc, char *argv[])
}
}
}
- // //test
- // std::cout << "before write polygon" << std::endl;
if(polygon_opt[0]||ruleMap[rule_opt[0]]==rule::mean){
//add ring to polygon
if(polygon_opt[0]){
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 689cb42..17c46ac 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -45,6 +45,9 @@ int main(int argc,char **argv) {
Optionpk<int> dimX_opt("dx", "dx", "filter kernel size in x, must be odd",
3);
Optionpk<int> dimY_opt("dy", "dy", "filter kernel size in y, must be odd",
3);
Optionpk<int> dimZ_opt("dz", "dz", "filter kernel size in z (band or
spectral dimension), must be odd (example: 3).. Set dz>0 if 1-D filter must be
used in band domain");
+ Optionpk<std::string> wavelet_type_opt("wt", "wavelet", "wavelet type:
daubechies,daubechies_centered, haar, haar_centered, bspline,
bspline_centered", "daubechies");
+ Optionpk<int> family_opt("wf", "family", "wavelet family (vanishing moment,
see also
http://www.gnu.org/software/gsl/manual/html_node/DWT-Initialization.html)", 4);
+ Optionpk<double> quantize_opt("q", "quantize", "Quantize threshold (positive
for percentage value, negative for absolute value)",0);
Optionpk<short> class_opt("class", "class", "class value(s) to use for
density, erosion, dilation, openening and closing, thresholding");
Optionpk<double> threshold_opt("t", "threshold", "threshold value(s) to use
for threshold filter (one for each class)", 0);
Optionpk<short> mask_opt("m", "mask", "mask value(s) ");
@@ -74,6 +77,9 @@ int main(int argc,char **argv) {
dimY_opt.retrieveOption(argc,argv);
dimZ_opt.retrieveOption(argc,argv);
option_opt.retrieveOption(argc,argv);
+ wavelet_type_opt.retrieveOption(argc,argv);
+ family_opt.retrieveOption(argc,argv);
+ quantize_opt.retrieveOption(argc,argv);
class_opt.retrieveOption(argc,argv);
threshold_opt.retrieveOption(argc,argv);
mask_opt.retrieveOption(argc,argv);
@@ -506,6 +512,14 @@ int main(int argc,char **argv) {
filter2d.smoothNoData(input,output,dimX_opt[0],dimY_opt[0]);
break;
}
+ case(filter2d::dwtForward):
+ filter2d.dwtForward(input, output, wavelet_type_opt[0], family_opt[0]);
+ break;
+ case(filter2d::dwtQuantize):
+ if(verbose_opt[0])
+ std::cout << "Quantization filtering" << std::endl;
+ filter2d.dwtQuantize(input, output, wavelet_type_opt[0], family_opt[0],
quantize_opt[0]);
+ break;
case(filter2d::threshold):
filter2d.setThresholds(threshold_opt);
filter2d.setClasses(class_opt);//deliberate fall through
diff --git a/src/apps/pkfilterascii.cc b/src/apps/pkfilterascii.cc
index 434a5a0..af3f787 100644
--- a/src/apps/pkfilterascii.cc
+++ b/src/apps/pkfilterascii.cc
@@ -36,9 +36,10 @@ int main(int argc,char **argv) {
Optionpk<std::string> input_opt("i","input","input ASCII file");
Optionpk<std::string> output_opt("o", "output", "Output ASCII file");
Optionpk<int> inputCols_opt("ic", "inputCols", "input columns (e.g., for
three dimensional input data in first three columns use: -ic 0 -ic 1 -ic 2");
- Optionpk<std::string> method_opt("f", "filter", "filter function (to be
implemented: dwtForward, dwtInverse)");
- Optionpk<std::string> wavelet_type_opt("wt", "wavelet", "wavelet type:
daubechies,daubechies_centered, haar, haar_centered, bspline,
bspline_centered", "daubechies");
- Optionpk<int> family_opt("wf", "family", "wavelet family (vanishing
moment, see also
http://www.gnu.org/software/gsl/manual/html_node/DWT-Initialization.html)", 4);
+ Optionpk<std::string> method_opt("f", "filter", "filter function (to be
implemented: dwtForward, dwtInverse,dwtQuantize)");
+ Optionpk<std::string> wavelet_type_opt("wt", "wavelet", "wavelet type:
daubechies,daubechies_centered, haar, haar_centered, bspline,
bspline_centered", "daubechies");
+ Optionpk<int> family_opt("wf", "family", "wavelet family (vanishing moment,
see also
http://www.gnu.org/software/gsl/manual/html_node/DWT-Initialization.html)", 4);
+ Optionpk<double> quantize_opt("q", "quantize", "Quantize threshold",0);
Optionpk<int> dimZ_opt("dz", "dz", "filter kernel size in z (band or
spectral dimension), must be odd (example: 3).. Set dz>0 if 1-D filter must be
used in band domain");
Optionpk<double> tapz_opt("tapz", "tapz", "taps used for spectral
filtering");
Optionpk<double> fwhm_opt("fwhm", "fwhm", "list of full width half to apply
spectral filtering (-fwhm band1 -fwhm band2 ...)");
@@ -55,6 +56,9 @@ int main(int argc,char **argv) {
output_opt.retrieveOption(argc,argv);
inputCols_opt.retrieveOption(argc,argv);
method_opt.retrieveOption(argc,argv);
+ wavelet_type_opt.retrieveOption(argc,argv);
+ family_opt.retrieveOption(argc,argv);
+ quantize_opt.retrieveOption(argc,argv);
dimZ_opt.retrieveOption(argc,argv);
tapz_opt.retrieveOption(argc,argv);
fwhm_opt.retrieveOption(argc,argv);
@@ -158,7 +162,7 @@ int main(int argc,char **argv) {
else{//no filtering
if(verbose_opt[0])
std::cout << "no filtering selected" << std::endl;
- wavelengthOut=wavelengthIn;
+ // wavelengthOut=wavelengthIn;
for(int icol=0;icol<inputCols_opt.size();++icol)
filteredData[icol]=inputData[icol];
}
@@ -166,12 +170,24 @@ int main(int argc,char **argv) {
if(method_opt.size()){
for(int icol=0;icol<inputCols_opt.size();++icol){
switch(filter::Filter::getFilterType(method_opt[0])){
- // case(filter::dwtForward):
- //
filter1d.dwtForward(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
- // break;
- // case(filter::dwtInverse):
- //
filter1d.dwtInverse(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
- // break;
+ case(filter::dwtForward):
+
filter1d.dwtForward(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
+ break;
+ case(filter::dwtInverse):
+
filter1d.dwtInverse(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
+ break;
+ case(filter::dwtQuantize):{
+ int origSize=filteredData[icol].size();
+
filter1d.dwtForward(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
+ for(int iband=0;iband<filteredData[icol].size();++iband){
+ if(fabs(filteredData[icol][iband])<quantize_opt[0])
+ filteredData[icol][iband]=0;
+ }
+
filter1d.dwtInverse(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
+ //remove extended samples
+
filteredData[icol].erase(filteredData[icol].begin()+origSize,filteredData[icol].end());
+ break;
+ }
default:
if(verbose_opt[0])
std::cout << "method to be implemented" << std::endl;
@@ -185,17 +201,17 @@ int main(int argc,char **argv) {
if(transpose_opt[0]){
for(int icol=0;icol<inputCols_opt.size();++icol){
- for(int iband=0;iband<wavelengthOut.size();++iband){
+ for(int iband=0;iband<filteredData[icol].size();++iband){
if(!output_opt.empty()){
outputStream << filteredData[icol][iband];
- if(iband<wavelengthOut.size()-1)
+ if(iband<filteredData[icol].size()-1)
outputStream << " ";
else
outputStream << std::endl;
}
else{
std::cout << filteredData[icol][iband];
- if(iband<wavelengthOut.size()-1)
+ if(iband<filteredData[icol].size()-1)
std::cout << " ";
else
std::cout << std::endl;
@@ -204,11 +220,20 @@ int main(int argc,char **argv) {
}
}
else{
- for(int iband=0;iband<wavelengthOut.size();++iband){
- if(!output_opt.empty())
- outputStream << wavelengthOut[iband] << " ";
- else
- std::cout << wavelengthOut[iband] << " ";
+ int nband=wavelengthOut.size()? wavelengthOut.size() :
filteredData[0].size();
+ for(int iband=0;iband<nband;++iband){
+ if(!output_opt.empty()){
+ if(wavelengthOut.size())
+ outputStream << wavelengthOut[iband] << " ";
+ else
+ outputStream << iband << " ";
+ }
+ else{
+ if(wavelengthOut.size())
+ std::cout << wavelengthOut[iband] << " ";
+ else
+ std::cout << iband << " ";
+ }
for(int icol=0;icol<inputCols_opt.size();++icol){
if(!output_opt.empty()){
outputStream << filteredData[icol][iband];
diff --git a/src/apps/pkinfo.cc b/src/apps/pkinfo.cc
index a8f15ee..cfba529 100644
--- a/src/apps/pkinfo.cc
+++ b/src/apps/pkinfo.cc
@@ -261,44 +261,46 @@ int main(int argc, char *argv[])
}
if(hist_opt[0]){
assert(band_opt[0]<imgReader.nrOfBand());
- imgReader.getMinMax(minValue,maxValue,band_opt[0]);
+ int nbin=nbin_opt[0];
+ // imgReader.getMinMax(minValue,maxValue,band_opt[0]);
+ // if(min_opt.size())
+ // minValue=min_opt[0];
+ // if(max_opt.size())
+ // maxValue=max_opt[0];
+ // if(nbin_opt[0]==0)
+ // nbin=maxValue-minValue+1;
+ // assert(nbin>0);
+ // std::vector<unsigned long int> output(nbin);
+ // unsigned long int nsample=0;
+ // unsigned long int ninvalid=0;
+ // std::vector<double> lineBuffer(imgReader.nrOfCol());
+ // for(int i=0;i<nbin;output[i++]=0);
+ // for(int irow=0;irow<imgReader.nrOfRow();++irow){
+ // imgReader.readData(lineBuffer,GDT_Float64,irow,band_opt[0]);
+ // for(int icol=0;icol<imgReader.nrOfCol();++icol){
+ // if(imgReader.isNoData(lineBuffer[icol]))
+ // ++ninvalid;
+ // else if(lineBuffer[icol]>maxValue)
+ // ++ninvalid;
+ // else if(lineBuffer[icol]<minValue)
+ // ++ninvalid;
+ // else if(lineBuffer[icol]==maxValue)
+ // ++output[nbin-1];
+ // else
+ //
++output[static_cast<int>(static_cast<double>(lineBuffer[icol]-minValue)/(maxValue-minValue)*nbin)];
+ // }
+ // }
+ std::vector<unsigned long int> output(nbin_opt[0]);
+ minValue=0;
+ maxValue=0;
if(min_opt.size())
- minValue=min_opt[0];
+ minValue=min_opt.size();
if(max_opt.size())
- maxValue=max_opt[0];
- int nbin=nbin_opt[0];
- if(nbin_opt[0]==0)
- nbin=maxValue-minValue+1;
- assert(nbin>0);
- std::vector<unsigned long int> output(nbin);
- unsigned long int nsample=0;
- unsigned long int ninvalid=0;
- std::vector<double> lineBuffer(imgReader.nrOfCol());
- for(int i=0;i<nbin;output[i++]=0);
- for(int irow=0;irow<imgReader.nrOfRow();++irow){
- imgReader.readData(lineBuffer,GDT_Float64,irow,band_opt[0]);
- for(int icol=0;icol<imgReader.nrOfCol();++icol){
- if(imgReader.isNoData(lineBuffer[icol]))
- ++ninvalid;
- else if(lineBuffer[icol]>maxValue)
- ++ninvalid;
- else if(lineBuffer[icol]<minValue)
- ++ninvalid;
- else if(lineBuffer[icol]==maxValue)
- ++output[nbin-1];
- // else
if(static_cast<double>(lineBuffer[icol]-minValue)/(maxValue-minValue)*nbin>=nbin){
- // //test
- // std::cout << "..." << lineBuffer[icol] << std::endl;
- // ++output[nbin-1];
- // }
- else
-
++output[static_cast<int>(static_cast<double>(lineBuffer[icol]-minValue)/(maxValue-minValue)*nbin)];
- }
- }
- nsample=imgReader.nrOfCol()*imgReader.nrOfRow()-ninvalid;
+ maxValue=max_opt.size();
+ unsigned long int
nsample=imgReader.getHistogram(output,minValue,maxValue,nbin,band_opt[0]);
std::cout.precision(10);
for(int bin=0;bin<nbin;++bin){
- nsample+=output[bin];
+ // nsample+=output[bin];
if(output[bin]>0){
std::cout << (maxValue-minValue)*bin/(nbin-1)+minValue << " ";
if(relative_opt[0])
diff --git a/src/imageclasses/ImgReaderGdal.cc
b/src/imageclasses/ImgReaderGdal.cc
index f28c076..7c1991c 100644
--- a/src/imageclasses/ImgReaderGdal.cc
+++ b/src/imageclasses/ImgReaderGdal.cc
@@ -445,6 +445,45 @@ void ImgReaderGdal::getMinMax(double& minValue, double&
maxValue, int band, bool
}
}
+unsigned long int ImgReaderGdal::getHistogram(vector<unsigned long int>&
histvector, double& min, double& max, int& nbin, int theBand) const{
+ double minValue=0;
+ double maxValue=0;
+ getMinMax(minValue,maxValue,theBand);
+ if(min<max&&min>minValue)
+ minValue=min;
+ else
+ min=minValue;
+ if(min<max&&max<maxValue)
+ maxValue=max;
+ else
+ max=maxValue;
+ if(nbin==0)
+ nbin=maxValue-minValue+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);
+ for(int irow=0;irow<nrOfRow();++irow){
+ readData(lineBuffer,GDT_Float64,irow,theBand);
+ for(int icol=0;icol<nrOfCol();++icol){
+ if(isNoData(lineBuffer[icol]))
+ ++ninvalid;
+ else if(lineBuffer[icol]>maxValue)
+ ++ninvalid;
+ else if(lineBuffer[icol]<minValue)
+ ++ninvalid;
+ else if(lineBuffer[icol]==maxValue)
+ ++histvector[nbin-1];
+ else
+
++histvector[static_cast<int>(static_cast<double>(lineBuffer[icol]-minValue)/(maxValue-minValue)*(nbin-1))];
+ }
+ }
+ unsigned long int nvalid=nrOfCol()*nrOfRow()-ninvalid;
+ return nvalid;
+}
+
void ImgReaderGdal::getRange(vector<short>& range, int band) const
{
vector<short> lineBuffer(nrOfCol());
diff --git a/src/imageclasses/ImgReaderGdal.h b/src/imageclasses/ImgReaderGdal.h
index f11c7e2..cb053d6 100644
--- a/src/imageclasses/ImgReaderGdal.h
+++ b/src/imageclasses/ImgReaderGdal.h
@@ -78,6 +78,7 @@ public:
void getMinMax(int startCol, int endCol, int startRow, int endRow, int band,
double& minValue, double& maxValue) const;
void getMinMax(double& minValue, double& maxValue, int band=0, bool
exhaustiveSearch=false) const;
double getMin(int& col, int& row, int band=0) const;
+ unsigned long int getHistogram(vector<unsigned long int>& histvector,
double& min, double& max,int& nbin, int theBand=0) const;
double getMax(int& col, int& row, int band=0) const;
void getRefPix(double& refX, double &refY, int band=0) const;
void getRange(vector<short>& range, int Band=0) const;
diff --git a/src/imageclasses/Makefile.am b/src/imageclasses/Makefile.am
index fceea0e..4af4a99 100644
--- a/src/imageclasses/Makefile.am
+++ b/src/imageclasses/Makefile.am
@@ -1,5 +1,5 @@
-AM_LDFLAGS = $(GDAL_LDFLAGS) @AM_LDFLAGS@
-AM_CXXFLAGS = -I$(top_srcdir)/src $(GDAL_CFLAGS) @AM_CXXFLAGS@
+AM_LDFLAGS = $(GSL_LDFLAGS) $(GDAL_LDFLAGS) @AM_LDFLAGS@
+AM_CXXFLAGS = -I$(top_srcdir)/src $(GSL_CFLAGS) $(GDAL_CFLAGS) @AM_CXXFLAGS@
###############################################################################
# THE LIBRARIES TO BUILD
--
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
