[pktools] 66/375: Filter1d and Filter2d with getFilterType returning enum FILTER_TYPE in function of string
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:53:59 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 9f33ca712fdfd38296dd5877e45b3510bf110260
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Fri Feb 22 17:52:09 2013 +0100
Filter1d and Filter2d with getFilterType returning enum FILTER_TYPE in function of string
---
src/algorithms/Filter.cc | 10 +-
src/algorithms/Filter.h | 241 ++++++++++++++++++++++----------------
src/algorithms/Filter2d.cc | 54 ++++-----
src/algorithms/Filter2d.h | 97 ++++++++++-----
src/algorithms/Makefile.am | 12 ++
src/algorithms/Makefile.in | 55 ++++++---
src/algorithms/StatFactory.h | 216 ++++++++++++++++++----------------
src/apps/pkcrop.cc | 84 +++++++++----
src/apps/pkdumpimg.cc | 180 ++++++++++++++--------------
src/apps/pkfilter.cc | 171 +++++++++++----------------
src/apps/pkinfo.cc | 8 +-
src/apps/pklas2img.cc | 16 +--
src/fileclasses/FileReaderAscii.h | 108 +++++++++++++++++
13 files changed, 748 insertions(+), 504 deletions(-)
diff --git a/src/algorithms/Filter.cc b/src/algorithms/Filter.cc
index f292314..b99a8df 100644
--- a/src/algorithms/Filter.cc
+++ b/src/algorithms/Filter.cc
@@ -38,7 +38,7 @@ void filter::Filter::setTaps(const vector<double> &taps)
assert(m_taps.size()%2);
}
-void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, int method, int dim, short down, int offset)
+void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down, int offset)
{
Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
@@ -104,7 +104,7 @@ void filter::Filter::doit(const ImgReaderGdal& input, ImgWriterGdal& output, sho
}
}
-void filter::Filter::applyFwhm(const ImgReaderGdal& input, ImgWriterGdal& output, const vector<double> &wavelengthIn, const vector<double> &wavelengthOut, const vector<double> &fwhm, bool verbose){
+void filter::Filter::applyFwhm(const vector<double> &wavelengthIn, const ImgReaderGdal& input, const vector<double> &wavelengthOut, const vector<double> &fwhm, const std::string& interpolationType, ImgWriterGdal& output, bool verbose){
Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
Vector2d<double> lineOutput(wavelengthOut.size(),input.nrOfCol());
const char* pszMessage;
@@ -118,8 +118,10 @@ void filter::Filter::applyFwhm(const ImgReaderGdal& input, ImgWriterGdal& output
vector<double> pixelInput(input.nrOfBand());
vector<double> pixelOutput;
for(int x=0;x<input.nrOfCol();++x){
- pixelInput=lineInput.selectCol(x);
- applyFwhm<double>(pixelInput,wavelengthIn,pixelOutput,wavelengthOut, fwhm, verbose);
+ for(int iband=0;iband<input.nrOfBand();++iband)
+ pixelInput[iband]=lineInput[iband][x];
+ applyFwhm<double>(wavelengthIn,pixelInput,wavelengthOut,fwhm, interpolationType, pixelOutput, verbose);
+ assert(pixelOutput.size()==wavelengthOut.size());
for(int iband=0;iband<pixelOutput.size();++iband)
lineOutput[iband][x]=pixelOutput[iband];
}
diff --git a/src/algorithms/Filter.h b/src/algorithms/Filter.h
index a4b6cb4..d20203e 100644
--- a/src/algorithms/Filter.h
+++ b/src/algorithms/Filter.h
@@ -33,141 +33,182 @@ 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};
+
class Filter
{
public:
- enum 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, SMOOTH=15, DENSITY=16, MAJORITY=17, MIXED=18};
Filter(void);
Filter(const vector<double> &taps);
virtual ~Filter(){};
+ FILTER_TYPE getFilterType(const std::string filterType){
+ std::map<std::string, FILTER_TYPE> m_filterMap;
+ initMap(m_filterMap);
+ return m_filterMap[filterType];
+ };
void setTaps(const vector<double> &taps);
void pushClass(short theClass=1){m_class.push_back(theClass);};
void pushMask(short theMask=0){m_mask.push_back(theMask);};
template<class T> void doit(const vector<T>& input, vector<T>& output, int down=1, int offset=0);
template<class T> void doit(T* input, int inputSize, vector<T>& output, int down=1, int offset=0);
- template<class T> void morphology(const vector<T>& input, vector<T>& output, int method, int dim, short down=1, int offset=0, bool verbose=0);
- void morphology(const ImgReaderGdal& input, ImgWriterGdal& output, int method, int dim, short down=1, int offset=0);
+ template<class T> void morphology(const vector<T>& input, vector<T>& output, const std::string& method, int dim, short down=1, int offset=0, bool verbose=0);
+ void morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down=1, int offset=0);
void doit(const ImgReaderGdal& input, ImgWriterGdal& output, short down=1, int offset=0);
template<class T> void applySrf(const Vector2d<T>& input, const Vector2d<double>& srf, Vector2d<T>& output, double delta=1, bool normalize=false, double centreWavelength=0, bool verbose=false);
- template<class T> void applyFwhm(const vector<double>& input, const vector<double> &wavelengthIn, vector<double>& output, const vector<double> &wavelengthOut, const vector<double> &fwhm, bool verbose=false);
- void applyFwhm(const ImgReaderGdal& input, ImgWriterGdal& output, const vector<double> &wavelengthIn, const vector<double> &wavelengthOut, const vector<double> &fwhm, bool verbose=false);
+ template<class T> void applyFwhm(const vector<double> &wavelengthIn, const vector<double>& input, const vector<double> &wavelengthOut, const vector<double> &fwhm, const std::string& interpolationType, vector<double>& output, bool verbose=false);
+ void applyFwhm(const vector<double> &wavelengthIn, const ImgReaderGdal& input, const vector<double> &wavelengthOut, const vector<double> &fwhm, const std::string& interpolationType, ImgWriterGdal& output, bool verbose=false);
// int fir(double* input, int nbandIn, vector<double>& output, int startBand, const string& wavelength, const string& fwhm, bool verbose);
// int fir(const vector<double>&input, vector<double>& output, int startBand, double fwhm, int ntaps, int down, int offset, bool verbose);
// int fir(double* input, int nbandIn, vector<double>& output, int startBand, double fwhm, int ntaps, int down, int offset, bool verbose);
private:
+ static void initMap(std::map<std::string, FILTER_TYPE>& m_filterMap){
+ //initialize selMap
+ m_filterMap["stdev"]=filter::stdev;
+ m_filterMap["var"]=filter::var;
+ m_filterMap["min"]=filter::min;
+ m_filterMap["max"]=filter::max;
+ m_filterMap["sum"]=filter::sum;
+ m_filterMap["mean"]=filter::mean;
+ m_filterMap["minmax"]=filter::minmax;
+ m_filterMap["dilate"]=filter::dilate;
+ m_filterMap["erode"]=filter::erode;
+ m_filterMap["close"]=filter::close;
+ m_filterMap["open"]=filter::open;
+ m_filterMap["homog"]=filter::homog;
+ m_filterMap["sobelx"]=filter::sobelx;
+ m_filterMap["sobely"]=filter::sobely;
+ m_filterMap["sobelxy"]=filter::sobelxy;
+ m_filterMap["sobelyx"]=filter::sobelyx;
+ m_filterMap["smooth"]=filter::smooth;
+ m_filterMap["density"]=filter::density;
+ m_filterMap["majority"]=filter::majority;
+ m_filterMap["mixed"]=filter::mixed;
+ m_filterMap["smoothnodata"]=filter::smoothnodata;
+ m_filterMap["threshold"]=filter::threshold;
+ m_filterMap["ismin"]=filter::ismin;
+ m_filterMap["ismax"]=filter::ismax;
+ m_filterMap["heterog"]=filter::heterog;
+ m_filterMap["order"]=filter::order;
+ m_filterMap["median"]=filter::median;
+ }
vector<double> m_taps;
vector<short> m_class;
vector<short> m_mask;
};
- template<class T> void Filter::applySrf(const Vector2d<T>& input, const Vector2d<double>& srf, Vector2d<T>& output, double delta, bool normalize, double centreWavelength, bool verbose)
-{
- output.resize(input.size());
- assert(srf.size()==2);//[0]: wavelength, [1]: response function
- assert(input.size()>1);//[0]: wavelength, [1],[2],...: value
- double start=floor(input[0][0]);
- double end=ceil(input[0].back());
- Vector2d<double> product(input.size());
- assert(input.size());
- assert(input[0].size()>1);
- int nband=srf[0].size();
- gsl_interp_accel *acc=gsl_interp_accel_alloc();
- gsl_spline *spline=gsl_spline_alloc(gsl_interp_linear,nband);
- gsl_spline_init(spline,&(srf[0][0]),&(srf[1][0]),nband);
-// double norm=gsl_spline_eval_integ(spline,start,end,acc);
- double norm=gsl_spline_eval_integ(spline,srf[0].front(),srf[0].back(),acc);
- gsl_spline_free(spline);
- gsl_interp_accel_free(acc);
- //interpolate input and srf to delta
- statfactory::StatFactory stat;
- Vector2d<double> input_d;
- Vector2d<double> srf_d;
- stat.interpolateUp(input,input_d,start,end,delta);
- stat.interpolateUp(srf,srf_d,start,end,delta);
- nband=input_d[0].size();
- if(verbose)
- cout << "number of interpolated bands: " << nband << endl;
- for(int isample=0;isample<input_d.size();++isample){
- product[isample].resize(nband);
- for(int iband=0;iband<nband;++iband){
- if(!isample)
- product[isample][iband]=input_d[isample][iband];
- else{
-// if(verbose&&isample==1)
-// cout << srf_d[0][iband] << " " << srf_d[1][iband] << endl;
- product[isample][iband]=input_d[isample][iband]*srf_d[1][iband];
- }
- }
- }
- output[0].resize(1);
- if(centreWavelength)
- output[0][0]=centreWavelength;
- else{
- double maxResponse=0;
- int maxIndex=0;
- for(int index=0;index<srf[1].size();++index){
- if(maxResponse<srf[1][index]){
- maxResponse=srf[1][index];
- maxIndex=index;
- }
- }
- output[0][0]=srf[0][maxIndex];
- }
- for(int isample=1;isample<input_d.size();++isample){
- output[isample].resize(1);
- gsl_interp_accel *acc=gsl_interp_accel_alloc();
- gsl_spline *spline=gsl_spline_alloc(gsl_interp_linear,nband);
- gsl_spline_init(spline,&(product[0][0]),&(product[isample][0]),nband);
- if(normalize)
- output[isample][0]=gsl_spline_eval_integ(spline,start,end,acc)/norm;
- else
- output[isample][0]=gsl_spline_eval_integ(spline,start,end,acc);
- gsl_spline_free(spline);
- gsl_interp_accel_free(acc);
- }
-}
+// template<class T> void Filter::applySrf(const Vector2d<T>& input, const Vector2d<double>& srf, Vector2d<T>& output, double delta, bool normalize, double centreWavelength, bool verbose)
+// {
+// output.resize(input.size());
+// assert(srf.size()==2);//[0]: wavelength, [1]: response function
+// assert(input.size()>1);//[0]: wavelength, [1],[2],...: value
+// double start=floor(input[0][0]);
+// double end=ceil(input[0].back());
+// Vector2d<double> product(input.size());
+// assert(input.size());
+// assert(input[0].size()>1);
+// int nband=srf[0].size();
+// gsl_interp_accel *acc=gsl_interp_accel_alloc();
+// gsl_spline *spline=gsl_spline_alloc(gsl_interp_linear,nband);
+// gsl_spline_init(spline,&(srf[0][0]),&(srf[1][0]),nband);
+// double norm=gsl_spline_eval_integ(spline,srf[0].front(),srf[0].back(),acc);
+// gsl_spline_free(spline);
+// gsl_interp_accel_free(acc);
+// //interpolate input and srf to delta
+// statfactory::StatFactory stat;
+// Vector2d<double> input_d;
+// Vector2d<double> srf_d;
+// stat.interpolateUp(input,input_d,start,end,delta,gsl_interp_polynomial);
+// stat.interpolateUp(srf,srf_d,start,end,delta,gsl_interp_polynomial);
+// nband=input_d[0].size();
+// if(verbose)
+// cout << "number of interpolated bands: " << nband << endl;
+// for(int isample=0;isample<input_d.size();++isample){
+// product[isample].resize(nband);
+// for(int iband=0;iband<nband;++iband){
+// if(!isample)
+// product[isample][iband]=input_d[isample][iband];
+// else{
+// product[isample][iband]=input_d[isample][iband]*srf_d[1][iband];
+// }
+// }
+// }
+// output[0].resize(1);
+// if(centreWavelength)
+// output[0][0]=centreWavelength;
+// else{
+// double maxResponse=0;
+// int maxIndex=0;
+// for(int index=0;index<srf[1].size();++index){
+// if(maxResponse<srf[1][index]){
+// maxResponse=srf[1][index];
+// maxIndex=index;
+// }
+// }
+// output[0][0]=srf[0][maxIndex];
+// }
+// for(int isample=1;isample<input_d.size();++isample){
+// output[isample].resize(1);
+// gsl_interp_accel *acc=gsl_interp_accel_alloc();
+// gsl_spline *spline=gsl_spline_alloc(gsl_interp_linear,nband);
+// gsl_spline_init(spline,&(product[0][0]),&(product[isample][0]),nband);
+// if(normalize)
+// output[isample][0]=gsl_spline_eval_integ(spline,start,end,acc)/norm;
+// else
+// output[isample][0]=gsl_spline_eval_integ(spline,start,end,acc);
+// gsl_spline_free(spline);
+// gsl_interp_accel_free(acc);
+// }
+// }
- template<class T> void Filter::applyFwhm(const vector<double>& input, const vector<double> &wavelengthIn, vector<double>& output, const vector<double> &wavelengthOut, const vector<double> &fwhm, bool verbose){
+template<class T> void Filter::applyFwhm(const vector<double> &wavelengthIn, const vector<double>& input, const vector<double> &wavelengthOut, const vector<double> &fwhm, const std::string& interpolationType, vector<double>& output, bool verbose){
double delta=1;//1 nm resolution
vector<double> stddev(fwhm.size());
for(int index=0;index<fwhm.size();++index)
stddev[index]=fwhm[index]/2.0/sqrt(2*log(2.0));//http://mathworld.wolfram.com/FullWidthatHalfMaximum.html
assert(wavelengthOut.size()==fwhm.size());
assert(wavelengthIn.size()==input.size());
- //todo densify input to 1 nm
assert(wavelengthIn[0]<wavelengthOut[0]);
- assert(wavelengthIn.back()<wavelengthOut.back());
+ assert(wavelengthIn.back()>wavelengthOut.back());
statfactory::StatFactory stat;
- Vector2d<double> input_course(1);
- input_course[0]=input;
- Vector2d<double> input_fine;
+ vector<double> input_fine;
vector<double> wavelength_fine;
- stat.interpolateUp(input_course,wavelengthIn,input_fine,wavelength_fine,wavelengthIn[0],wavelengthIn.back(),delta);
+ for(double win=floor(wavelengthIn[0]);win<=ceil(wavelengthIn.back());win+=delta)
+ wavelength_fine.push_back(win);
+ if(verbose){
+ for(int index=0;index<wavelength_fine.size();++index)
+ std::cout << " " << wavelength_fine[index];
+ std::cout << std::endl;
+ std::cout << "interpolate input wavelength to " << delta << " nm resolution (size=" << wavelength_fine.size() << ")" << std::endl;
+ }
+ stat.interpolateUp(wavelengthIn,input,wavelength_fine,interpolationType,input_fine,verbose);
int nbandIn=wavelength_fine.size();
+
int nbandOut=wavelengthOut.size();
+ output.resize(nbandOut);
gsl::matrix tf(nbandIn,nbandOut);
for(int indexOut=0;indexOut<nbandOut;++indexOut){
- for(int indexIn=0;indexIn<wavelengthIn.size();++indexIn){
+ double norm=0;
+ for(int indexIn=0;indexIn<nbandIn;++indexIn){
tf(indexIn,indexOut)=
- exp((wavelengthOut[indexOut]-wavelengthIn[indexIn])
- *(wavelengthIn[indexIn]-wavelengthOut[indexOut])
- /2.0/stddev[indexOut]
- /stddev[indexOut]);
+ exp((wavelengthOut[indexOut]-wavelength_fine[indexIn])
+ *(wavelength_fine[indexIn]-wavelengthOut[indexOut])
+ /2.0/stddev[indexOut]
+ /stddev[indexOut]);
tf(indexIn,indexOut)/=sqrt(2.0*M_PI);
tf(indexIn,indexOut)/=stddev[indexOut];
+ norm+=tf(indexIn,indexOut);
}
- double norm=exp(tf.LU_lndet());//(tf.Column(indexOut+1)).NormFrobenius();
- if(norm)
- tf.get_col(indexOut+1)/=norm;
- }
+ //todo: check how to normalize...
+ // double norm=exp((tf.transpose()*tf).LU_lndet());//(tf.Column(indexOut+1)).NormFrobenius();
+ // if(norm)
+ // tf.get_col(indexOut+1)/=norm;
//create filtered vector
- output.resize(nbandOut);
- for(int indexOut=0;indexOut<nbandOut;++indexOut){
+ //todo: support more than one sample
output[indexOut]=0;
for(int indexIn=0;indexIn<nbandIn;++indexIn)
- output[indexOut]+=input_fine[0][indexIn]*tf(indexIn,indexOut);
+ output[indexOut]+=input_fine[indexIn]*tf(indexIn,indexOut)/norm;
}
}
@@ -203,7 +244,7 @@ template<class T> void Filter::doit(const vector<T>& input, vector<T>& output, i
}
}
-template<class T> void Filter::morphology(const vector<T>& input, vector<T>& output, int method, int dim, short down, int offset, bool verbose)
+template<class T> void Filter::morphology(const vector<T>& input, vector<T>& output, const std::string& method, int dim, short down, int offset, bool verbose)
{
assert(dim);
output.resize((input.size()-offset+down-1)/down);
@@ -246,11 +287,11 @@ template<class T> void Filter::morphology(const vector<T>& input, vector<T>& out
statBuffer.clear();
continue;
}
- switch(method){
- case(DILATE):
+ switch(getFilterType(method)){
+ case(filter::dilate):
output[(i-offset+down-1)/down]=stat.max(statBuffer);
break;
- case(ERODE):
+ case(filter::erode):
output[(i-offset+down-1)/down]=stat.min(statBuffer);
break;
default:
@@ -286,11 +327,11 @@ template<class T> void Filter::morphology(const vector<T>& input, vector<T>& out
statBuffer.clear();
continue;
}
- switch(method){
- case(DILATE):
+ switch(getFilterType(method)){
+ case(filter::dilate):
output[(i-offset+down-1)/down]=stat.max(statBuffer);
break;
- case(ERODE):
+ case(filter::erode):
output[(i-offset+down-1)/down]=stat.min(statBuffer);
break;
default:
@@ -340,11 +381,11 @@ template<class T> void Filter::morphology(const vector<T>& input, vector<T>& out
statBuffer.clear();
continue;
}
- switch(method){
- case(DILATE):
+ switch(getFilterType(method)){
+ case(filter::dilate):
output[(i-offset+down-1)/down]=stat.max(statBuffer);
break;
- case(ERODE):
+ case(filter::erode):
output[(i-offset+down-1)/down]=stat.min(statBuffer);
break;
default:
diff --git a/src/algorithms/Filter2d.cc b/src/algorithms/Filter2d.cc
index 0ea019e..d4eb580 100644
--- a/src/algorithms/Filter2d.cc
+++ b/src/algorithms/Filter2d.cc
@@ -389,7 +389,7 @@ void filter2d::Filter2d::median(const string& inputFilename, const string& outpu
ImgWriterGdal output;
input.open(inputFilename);
output.open(outputFilename,input);
- doit(input,output,MEDIAN,dim,disc);
+ doit(input,output,"median",dim,disc);
}
void filter2d::Filter2d::var(const string& inputFilename, const string& outputFilename,int dim, bool disc)
@@ -398,15 +398,15 @@ void filter2d::Filter2d::var(const string& inputFilename, const string& outputFi
ImgWriterGdal output;
input.open(inputFilename);
output.open(outputFilename,input);
- doit(input,output,VAR,dim,disc);
+ doit(input,output,"var",dim,disc);
}
-void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output, int method, int dim, short down, bool disc)
+void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down, bool disc)
{
doit(input,output,method,dim,dim,down,disc);
}
-void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output, int method, int dimX, int dimY, short down, bool disc)
+void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, short down, bool disc)
{
assert(dimX);
assert(dimY);
@@ -501,49 +501,49 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
}
}
}
- switch(method){
- case(MEDIAN):
+ switch(getFilterType(method)){
+ case(filter2d::median):
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=stat.median(windowBuffer);
break;
- case(VAR):{
+ case(filter2d::var):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=stat.var(windowBuffer);
break;
}
- case(STDEV):{
+ case(filter2d::stdev):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=sqrt(stat.var(windowBuffer));
break;
}
- case(MEAN):{
+ case(filter2d::mean):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=stat.mean(windowBuffer);
break;
}
- case(MIN):{
+ case(filter2d::min):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=stat.min(windowBuffer);
break;
}
- case(ISMIN):{
+ case(filter2d::ismin):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=(stat.min(windowBuffer)==windowBuffer[dimX*dimY/2])? 1:0;
break;
}
- case(MINMAX):{
+ case(filter2d::minmax):{
double min=0;
double max=0;
if(windowBuffer.empty())
@@ -557,21 +557,21 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
}
break;
}
- case(MAX):{
+ case(filter2d::max):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=stat.max(windowBuffer);
break;
}
- case(ISMAX):{
+ case(filter2d::ismax):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=(stat.max(windowBuffer)==windowBuffer[dimX*dimY/2])? 1:0;
break;
}
- case(ORDER):{
+ case(filter2d::order):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else{
@@ -584,17 +584,17 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
}
break;
}
- case(SUM):{
+ case(filter2d::sum):{
outBuffer[x/down]=stat.sum(windowBuffer);
break;
}
- case(HOMOG):
+ case(filter2d::homog):
if(occurrence.size()==1)//all values in window must be the same
outBuffer[x/down]=inBuffer[dimY/2][x];
else//favorize original value in case of ties
outBuffer[x/down]=m_noValue;
break;
- case(HETEROG):{
+ case(filter2d::heterog):{
for(vector<double>::const_iterator wit=windowBuffer.begin();wit!=windowBuffer.end();++wit){
if(wit==windowBuffer.begin()+windowBuffer.size()/2)
continue;
@@ -607,7 +607,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
}
break;
}
- case(DENSITY):{
+ case(filter2d::density):{
if(windowBuffer.size()){
vector<short>::const_iterator vit=m_class.begin();
while(vit!=m_class.end())
@@ -617,7 +617,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
outBuffer[x/down]=m_noValue;
break;
}
- case(MAJORITY):{
+ case(filter2d::majority):{
if(occurrence.size()){
map<int,int>::const_iterator maxit=occurrence.begin();
for(map<int,int>::const_iterator mit=occurrence.begin();mit!=occurrence.end();++mit){
@@ -633,7 +633,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
outBuffer[x/down]=m_noValue;
break;
}
- case(THRESHOLD):{
+ case(filter2d::threshold):{
assert(m_class.size()==m_threshold.size());
if(windowBuffer.size()){
outBuffer[x/down]=inBuffer[dimY/2][x];//initialize with original value (in case thresholds not met)
@@ -646,7 +646,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
outBuffer[x/down]=m_noValue;
break;
}
- case(MIXED):{
+ case(filter2d::mixed):{
enum Type { BF=11, CF=12, MF=13, NF=20, W=30 };
double nBF=occurrence[BF];
double nCF=occurrence[CF];
@@ -799,7 +799,7 @@ void filter2d::Filter2d::doit(const ImgReaderGdal& input, ImgWriterGdal& output,
// output.close();
// }
-void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, int method, int dimX, int dimY, bool disc, double angle)
+void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, bool disc, double angle)
{
assert(dimX);
assert(dimY);
@@ -932,16 +932,16 @@ void filter2d::Filter2d::morphology(const ImgReaderGdal& input, ImgWriterGdal& o
}
}
if(statBuffer.size()){
- switch(method){
- case(DILATE):
+ switch(getFilterType(method)){
+ case(filter2d::dilate):
outBuffer[x]=stat.max(statBuffer);
break;
- case(ERODE):
+ case(filter2d::erode):
outBuffer[x]=stat.min(statBuffer);
break;
default:
ostringstream ess;
- ess << "Error: morphology method " << method << " not supported, choose " << DILATE << " (dilate) or " << ERODE << " (erode)" << endl;
+ ess << "Error: morphology method " << method << " not supported, choose " << filter2d::dilate << " (dilate) or " << filter2d::erode << " (erode)" << endl;
throw(ess.str());
break;
}
diff --git a/src/algorithms/Filter2d.h b/src/algorithms/Filter2d.h
index ed45d1f..fd1abc9 100644
--- a/src/algorithms/Filter2d.h
+++ b/src/algorithms/Filter2d.h
@@ -43,7 +43,7 @@ using namespace std;
// using namespace cimg_library;
namespace filter2d
{
- enum 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};
+ 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};
class Filter2d
{
@@ -51,6 +51,12 @@ public:
Filter2d(void);
Filter2d(const Vector2d<double> &taps);
virtual ~Filter2d(){};
+ static FILTER_TYPE getFilterType(const std::string filterType){
+ std::map<std::string, FILTER_TYPE> m_filterMap;
+ initMap(m_filterMap);
+ return m_filterMap[filterType];
+ };
+
void setTaps(const Vector2d<double> &taps);
void setNoValue(double noValue=0){m_noValue=noValue;};
void pushClass(short theClass=1){m_class.push_back(theClass);};
@@ -68,18 +74,49 @@ public:
template<class T1, class T2> void smooth(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector,int dimX, int dimY);
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, int method, int dim, short down=2, bool disc=false);
- void doit(const ImgReaderGdal& input, ImgWriterGdal& output, int method, int dimX, int dimY, short down=1, bool disc=false);
- template<class T1, class T2> void doit(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector, int method, int dimX, int dimY, short down=1, bool disc=false);
+ void doit(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down=2, bool disc=false);
+ void doit(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, short down=1, bool disc=false);
+ template<class T1, class T2> void doit(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector, const std::string& method, int dimX, int dimY, short down=1, bool disc=false);
void median(const string& inputFilename, const string& outputFilename, int dim, bool disc=false);
void var(const string& inputFilename, const string& outputFilename, int dim, bool disc=false);
- void morphology(const ImgReaderGdal& input, ImgWriterGdal& output, int method, int dimX, int dimY, bool disc=false, double angle=-190);
- template<class T> unsigned long int morphology(const Vector2d<T>& input, Vector2d<T>& output, int method, int dimX, int dimY, bool disc=false, double hThreshold=0);
+ void morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, bool disc=false, double angle=-190);
+ template<class T> unsigned long int morphology(const Vector2d<T>& input, Vector2d<T>& output, const std::string& method, int dimX, int dimY, bool disc=false, double hThreshold=0);
template<class T> void shadowDsm(const Vector2d<T>& input, Vector2d<T>& output, double sza, double saa, double pixelSize, short shadowFlag=1);
void shadowDsm(const ImgReaderGdal& input, ImgWriterGdal& output, double sza, double saa, double pixelSize, short shadowFlag=1);
void dwt_texture(const string& inputFilename, const string& outputFilename,int dim, int scale, int down=1, int iband=0, bool verbose=false);
private:
+ static void initMap(std::map<std::string, FILTER_TYPE>& m_filterMap){
+ //initialize selMap
+ m_filterMap["stdev"]=filter2d::stdev;
+ m_filterMap["var"]=filter2d::var;
+ m_filterMap["min"]=filter2d::min;
+ m_filterMap["max"]=filter2d::max;
+ m_filterMap["sum"]=filter2d::sum;
+ m_filterMap["mean"]=filter2d::mean;
+ m_filterMap["minmax"]=filter2d::minmax;
+ m_filterMap["dilate"]=filter2d::dilate;
+ m_filterMap["erode"]=filter2d::erode;
+ m_filterMap["close"]=filter2d::close;
+ m_filterMap["open"]=filter2d::open;
+ m_filterMap["homog"]=filter2d::homog;
+ m_filterMap["sobelx"]=filter2d::sobelx;
+ m_filterMap["sobely"]=filter2d::sobely;
+ m_filterMap["sobelxy"]=filter2d::sobelxy;
+ m_filterMap["sobelyx"]=filter2d::sobelyx;
+ m_filterMap["smooth"]=filter2d::smooth;
+ m_filterMap["density"]=filter2d::density;
+ m_filterMap["majority"]=filter2d::majority;
+ m_filterMap["mixed"]=filter2d::mixed;
+ m_filterMap["smoothnodata"]=filter2d::smoothnodata;
+ m_filterMap["threshold"]=filter2d::threshold;
+ m_filterMap["ismin"]=filter2d::ismin;
+ m_filterMap["ismax"]=filter2d::ismax;
+ m_filterMap["heterog"]=filter2d::heterog;
+ m_filterMap["order"]=filter2d::order;
+ m_filterMap["median"]=filter2d::median;
+ }
+
Vector2d<double> m_taps;
double m_noValue;
vector<short> m_class;
@@ -150,7 +187,7 @@ private:
}
}
-template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector, int method, int dimX, int dimY, short down, bool disc)
+template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector, const std::string& method, int dimX, int dimY, short down, bool disc)
{
statfactory::StatFactory stat;
outputVector.resize((inputVector.size()+down-1)/down);
@@ -227,49 +264,49 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
}
}
}
- switch(method){
- case(MEDIAN):
+ switch(getFilterType(method)){
+ case(filter2d::median):
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=stat.median(windowBuffer);
break;
- case(VAR):{
+ case(filter2d::var):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=stat.var(windowBuffer);
break;
}
- case(STDEV):{
+ case(filter2d::stdev):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=sqrt(stat.var(windowBuffer));
break;
}
- case(MEAN):{
+ case(filter2d::mean):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=stat.mean(windowBuffer);
break;
}
- case(MIN):{
+ case(filter2d::min):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=stat.min(windowBuffer);
break;
}
- case(ISMIN):{
+ case(filter2d::ismin):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=(stat.min(windowBuffer)==windowBuffer[dimX*dimY/2])? 1:0;
break;
}
- case(MINMAX):{
+ case(filter2d::minmax):{
double min=0;
double max=0;
if(windowBuffer.empty())
@@ -283,21 +320,21 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
}
break;
}
- case(MAX):{
+ case(filter2d::max):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=stat.max(windowBuffer);
break;
}
- case(ISMAX):{
+ case(filter2d::ismax):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else
outBuffer[x/down]=(stat.max(windowBuffer)==windowBuffer[dimX*dimY/2])? 1:0;
break;
}
- case(ORDER):{
+ case(filter2d::order):{
if(windowBuffer.empty())
outBuffer[x/down]=m_noValue;
else{
@@ -310,17 +347,17 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
}
break;
}
- case(SUM):{
+ case(filter2d::sum):{
outBuffer[x/down]=stat.sum(windowBuffer);
break;
}
- case(HOMOG):
+ case(filter2d::homog):
if(occurrence.size()==1)//all values in window must be the same
outBuffer[x/down]=inBuffer[dimY/2][x];
else//favorize original value in case of ties
outBuffer[x/down]=m_noValue;
break;
- case(HETEROG):{
+ case(filter2d::heterog):{
for(vector<double>::const_iterator wit=windowBuffer.begin();wit!=windowBuffer.end();++wit){
if(wit==windowBuffer.begin()+windowBuffer.size()/2)
continue;
@@ -333,7 +370,7 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
}
break;
}
- case(DENSITY):{
+ case(filter2d::density):{
if(windowBuffer.size()){
vector<short>::const_iterator vit=m_class.begin();
while(vit!=m_class.end())
@@ -343,7 +380,7 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
outBuffer[x/down]=m_noValue;
break;
}
- case(MAJORITY):{
+ case(filter2d::majority):{
if(occurrence.size()){
map<int,int>::const_iterator maxit=occurrence.begin();
for(map<int,int>::const_iterator mit=occurrence.begin();mit!=occurrence.end();++mit){
@@ -359,7 +396,7 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
outBuffer[x/down]=m_noValue;
break;
}
- case(THRESHOLD):{
+ case(filter2d::threshold):{
assert(m_class.size()==m_threshold.size());
if(windowBuffer.size()){
outBuffer[x/down]=inBuffer[dimY/2][x];//initialize with original value (in case thresholds not met)
@@ -372,7 +409,7 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
outBuffer[x/down]=m_noValue;
break;
}
- case(MIXED):{
+ case(filter2d::mixed):{
enum MixType { BF=11, CF=12, MF=13, NF=20, W=30 };
double nBF=occurrence[BF];
double nCF=occurrence[CF];
@@ -419,7 +456,7 @@ template<class T1, class T2> void Filter2d::doit(const Vector2d<T1>& inputVector
// }
// };
-template<class T> unsigned long int Filter2d::morphology(const Vector2d<T>& input, Vector2d<T>& output, int method, int dimX, int dimY, bool disc, double hThreshold)
+template<class T> unsigned long int Filter2d::morphology(const Vector2d<T>& input, Vector2d<T>& output, const std::string& method, int dimX, int dimY, bool disc, double hThreshold)
{
unsigned long int nchange=0;
assert(dimX);
@@ -510,14 +547,14 @@ template<class T> unsigned long int Filter2d::morphology(const Vector2d<T>& inpu
}
}
if(statBuffer.size()){
- switch(method){
- case(DILATE):
+ switch(getFilterType(method)){
+ case(filter2d::dilate):
if(output[y][x]<stat.max(statBuffer)-hThreshold){
output[y][x]=stat.max(statBuffer);
++nchange;
}
break;
- case(ERODE):
+ case(filter2d::erode):
if(output[y][x]>stat.min(statBuffer)+hThreshold){
output[y][x]=stat.min(statBuffer);
++nchange;
@@ -525,7 +562,7 @@ template<class T> unsigned long int Filter2d::morphology(const Vector2d<T>& inpu
break;
default:
ostringstream ess;
- ess << "Error: morphology method " << method << " not supported, choose " << DILATE << " (dilate) or " << ERODE << " (erode)" << endl;
+ ess << "Error: morphology method " << method << " not supported, choose " << filter2d::dilate << " (dilate) or " << filter2d::erode << " (erode)" << endl;
throw(ess.str());
break;
}
diff --git a/src/algorithms/Makefile.am b/src/algorithms/Makefile.am
index 3ef425a..be84cb4 100644
--- a/src/algorithms/Makefile.am
+++ b/src/algorithms/Makefile.am
@@ -2,6 +2,14 @@ AM_LDFLAGS = $(GSL_LIBS) $(GDAL_LDFLAGS) @AM_LDFLAGS@
AM_CXXFLAGS = -I$(top_srcdir)/src $(GDAL_CFLAGS) @AM_CXXFLAGS@
###############################################################################
+# THE PROGRAMS TO BUILD
+###############################################################################
+
+# the program to build (the names of the final binaries)
+#do not want to install pktestoption
+noinst_PROGRAMS = pktestStat
+
+###############################################################################
# THE LIBRARIES TO BUILD
###############################################################################
@@ -25,3 +33,7 @@ endif
# the sources to add to the library and to add to the source distribution
libalgorithms_a_SOURCES = $(libalgorithms_a_HEADERS) Egcs.cc Filter2d.cc Filter.cc ConfusionMatrix.cc svm.cpp
###############################################################################
+
+# list of sources for the binaries
+pktestStat_SOURCES = pktestStat.cc
+pktestStat_LDADD = $(GSL_LIBS) -lgslwrap $(top_builddir)/src/algorithms/libalgorithms.a $(top_builddir)/src/imageclasses/libimageClasses.a -lgslwrap
diff --git a/src/algorithms/Makefile.in b/src/algorithms/Makefile.in
index c6d7f0e..a0dda50 100644
--- a/src/algorithms/Makefile.in
+++ b/src/algorithms/Makefile.in
@@ -16,6 +16,7 @@
@SET_MAKE@
+
VPATH = @srcdir@
pkgdatadir = $(datadir)/@PACKAGE@
pkgincludedir = $(includedir)/@PACKAGE@
@@ -33,6 +34,7 @@ POST_INSTALL = :
NORMAL_UNINSTALL = :
PRE_UNINSTALL = :
POST_UNINSTALL = :
+noinst_PROGRAMS = pktestStat$(EXEEXT)
@USE_FANN_TRUE at am__append_1 = myfann_cpp.h
@USE_NLOPT_TRUE at am__append_2 = SensorModel.h OptFactory.h
subdir = src/algorithms
@@ -62,6 +64,13 @@ am_libalgorithms_a_OBJECTS = $(am__objects_2) Egcs.$(OBJEXT) \
Filter2d.$(OBJEXT) Filter.$(OBJEXT) ConfusionMatrix.$(OBJEXT) \
svm.$(OBJEXT)
libalgorithms_a_OBJECTS = $(am_libalgorithms_a_OBJECTS)
+PROGRAMS = $(noinst_PROGRAMS)
+am_pktestStat_OBJECTS = pktestStat.$(OBJEXT)
+pktestStat_OBJECTS = $(am_pktestStat_OBJECTS)
+am__DEPENDENCIES_1 =
+pktestStat_DEPENDENCIES = $(am__DEPENDENCIES_1) \
+ $(top_builddir)/src/algorithms/libalgorithms.a \
+ $(top_builddir)/src/imageclasses/libimageClasses.a
DEFAULT_INCLUDES = -I. at am__isrc@ -I$(top_builddir)
depcomp = $(SHELL) $(top_srcdir)/depcomp
am__depfiles_maybe = depfiles
@@ -75,8 +84,9 @@ COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \
$(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS)
CCLD = $(CC)
LINK = $(CCLD) $(AM_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) $(LDFLAGS) -o $@
-SOURCES = $(libalgorithms_a_SOURCES)
-DIST_SOURCES = $(am__libalgorithms_a_SOURCES_DIST)
+SOURCES = $(libalgorithms_a_SOURCES) $(pktestStat_SOURCES)
+DIST_SOURCES = $(am__libalgorithms_a_SOURCES_DIST) \
+ $(pktestStat_SOURCES)
am__libalgorithms_a_HEADERS_DIST = Egcs.h Filter2d.h Filter.h \
StatFactory.h ConfusionMatrix.h svm.h FeatureSelector.h \
myfann_cpp.h SensorModel.h OptFactory.h
@@ -234,6 +244,11 @@ libalgorithms_a_HEADERS = Egcs.h Filter2d.h Filter.h StatFactory.h \
# the sources to add to the library and to add to the source distribution
libalgorithms_a_SOURCES = $(libalgorithms_a_HEADERS) Egcs.cc Filter2d.cc Filter.cc ConfusionMatrix.cc svm.cpp
+###############################################################################
+
+# list of sources for the binaries
+pktestStat_SOURCES = pktestStat.cc
+pktestStat_LDADD = $(GSL_LIBS) -lgslwrap $(top_builddir)/src/algorithms/libalgorithms.a $(top_builddir)/src/imageclasses/libimageClasses.a -lgslwrap
all: all-am
.SUFFIXES:
@@ -276,6 +291,12 @@ libalgorithms.a: $(libalgorithms_a_OBJECTS) $(libalgorithms_a_DEPENDENCIES)
$(libalgorithms_a_AR) libalgorithms.a $(libalgorithms_a_OBJECTS) $(libalgorithms_a_LIBADD)
$(RANLIB) libalgorithms.a
+clean-noinstPROGRAMS:
+ -test -z "$(noinst_PROGRAMS)" || rm -f $(noinst_PROGRAMS)
+pktestStat$(EXEEXT): $(pktestStat_OBJECTS) $(pktestStat_DEPENDENCIES)
+ @rm -f pktestStat$(EXEEXT)
+ $(CXXLINK) $(pktestStat_OBJECTS) $(pktestStat_LDADD) $(LIBS)
+
mostlyclean-compile:
-rm -f *.$(OBJEXT)
@@ -286,6 +307,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/Egcs.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/Filter.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/Filter2d.Po at am__quote@
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pktestStat.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/svm.Po at am__quote@
.cc.o:
@@ -420,7 +442,7 @@ distdir: $(DISTFILES)
done
check-am: all-am
check: check-am
-all-am: Makefile $(LIBRARIES) $(HEADERS)
+all-am: Makefile $(LIBRARIES) $(PROGRAMS) $(HEADERS)
installdirs:
for dir in "$(DESTDIR)$(libalgorithms_adir)"; do \
test -z "$$dir" || $(MKDIR_P) "$$dir"; \
@@ -452,7 +474,8 @@ maintainer-clean-generic:
@echo "it deletes files that may require special tools to rebuild."
clean: clean-am
-clean-am: clean-generic clean-noinstLIBRARIES mostlyclean-am
+clean-am: clean-generic clean-noinstLIBRARIES clean-noinstPROGRAMS \
+ mostlyclean-am
distclean: distclean-am
-rm -rf ./$(DEPDIR)
@@ -522,19 +545,19 @@ uninstall-am: uninstall-libalgorithms_aHEADERS
.MAKE: install-am install-strip
.PHONY: CTAGS GTAGS all all-am check check-am clean clean-generic \
- clean-noinstLIBRARIES ctags distclean distclean-compile \
- distclean-generic distclean-tags distdir dvi dvi-am html \
- html-am info info-am install install-am install-data \
- install-data-am install-dvi install-dvi-am install-exec \
- install-exec-am install-html install-html-am install-info \
- install-info-am install-libalgorithms_aHEADERS install-man \
- install-pdf install-pdf-am install-ps install-ps-am \
- install-strip installcheck installcheck-am installdirs \
- maintainer-clean maintainer-clean-generic mostlyclean \
- mostlyclean-compile mostlyclean-generic pdf pdf-am ps ps-am \
- tags uninstall uninstall-am uninstall-libalgorithms_aHEADERS
+ clean-noinstLIBRARIES clean-noinstPROGRAMS ctags distclean \
+ distclean-compile distclean-generic distclean-tags distdir dvi \
+ dvi-am html html-am info info-am install install-am \
+ install-data install-data-am install-dvi install-dvi-am \
+ install-exec install-exec-am install-html install-html-am \
+ install-info install-info-am install-libalgorithms_aHEADERS \
+ install-man install-pdf install-pdf-am install-ps \
+ install-ps-am install-strip installcheck installcheck-am \
+ installdirs maintainer-clean maintainer-clean-generic \
+ mostlyclean mostlyclean-compile mostlyclean-generic pdf pdf-am \
+ ps ps-am tags uninstall uninstall-am \
+ uninstall-libalgorithms_aHEADERS
-###############################################################################
# Tell versions [3.59,3.63) of GNU make to not export all variables.
# Otherwise a system limit (for SysV at least) may be exceeded.
diff --git a/src/algorithms/StatFactory.h b/src/algorithms/StatFactory.h
index 2e38ef7..350a49d 100644
--- a/src/algorithms/StatFactory.h
+++ b/src/algorithms/StatFactory.h
@@ -22,6 +22,7 @@ along with pktools. If not, see <http://www.gnu.org/licenses/>.
#include <iostream>
#include <vector>
+#include <map>
#include <math.h>
#include <assert.h>
#include <string>
@@ -38,12 +39,57 @@ using namespace std;
namespace statfactory
{
- enum INTERPOLATION_TYPE {UNDEFINED=0,POLYNOMIAL=1,CSPLINE=2,PERIODIC=3,AKIMA=4,AKIMA_PERIODIC=5,LINEAR=6};
-
class StatFactory{
+
public:
+ enum INTERPOLATION_TYPE {undefined=0,linear=1,polynomial=2,cspline=3,cspline_periodic=4,akima=5,akima_periodic=6};
StatFactory(void){};
virtual ~StatFactory(void){};
+ INTERPOLATION_TYPE getInterpolationType(const std::string interpolationType){
+ std::map<std::string, INTERPOLATION_TYPE> m_interpMap;
+ initMap(m_interpMap);
+ // gsl_interp_accel *acc=gsl_interp_accel_alloc();
+ // gsl_spline *spline=gsl_spline_alloc(m_interpMap["interpolationType"],theSize);
+ return m_interpMap[interpolationType];
+ };
+ static void allocAcc(gsl_interp_accel *&acc){
+ acc = gsl_interp_accel_alloc ();
+ };
+
+ static void getSpline(const std::string type, int size, gsl_spline *& spline){
+ std::map<std::string, INTERPOLATION_TYPE> m_interpMap;
+ initMap(m_interpMap);
+ switch(m_interpMap[type]){
+ case(polynomial):
+ spline=gsl_spline_alloc(gsl_interp_polynomial,size);
+ break;
+ case(cspline):
+ spline=gsl_spline_alloc(gsl_interp_cspline,size);
+ break;
+ case(cspline_periodic):
+ spline=gsl_spline_alloc(gsl_interp_cspline_periodic,size);
+ break;
+ case(akima):
+ spline=gsl_spline_alloc(gsl_interp_akima,size);
+ break;
+ case(akima_periodic):
+ spline=gsl_spline_alloc(gsl_interp_akima_periodic,size);
+ break;
+ case(linear):
+ default:
+ spline=gsl_spline_alloc(gsl_interp_linear,size);
+ break;
+ }
+ assert(spline);
+ };
+ static void initSpline(gsl_spline *spline, const double *x, const double *y, int size){
+ gsl_spline_init (spline, x, y, size);
+ };
+ static double evalSpline(gsl_spline *spline, double x, gsl_interp_accel *acc){
+ return gsl_spline_eval (spline, x, acc);
+ };
+
+
template<class T> T max(const vector<T>& v) const;
template<class T> T min(const vector<T>& v) const;
// template<class T> typename vector<T>::const_iterator max(const vector<T>& v, typename vector<T>::const_iterator begin, typename vector<T>::const_iterator end) const;
@@ -75,16 +121,29 @@ public:
template<class T> double correlation(const vector<T>& x, const vector<T>& y, int delay=0) const;
template<class T> double cross_correlation(const vector<T>& x, const vector<T>& y, int maxdelay, vector<T>& z) const;
template<class T> double linear_regression(const vector<T>& x, const vector<T>& y, double &c0, double &c1) const;
- template<class T> void interpolateUp(const vector< vector<T> >& input, vector< vector<T> >& output, double start, double end, double step, int type=1);
- template<class T> void interpolateUp(const vector< vector<T> >& input, const vector<double>& wavelengthIn, vector< vector<T> >& output, vector<double>& wavelengthOut, double start, double end, double step, int type=1);
+ template<class T> void interpolateUp(const vector<double>& wavelengthIn, const vector<T>& input, const vector<double>& wavelengthOut, const std::string& type, vector<T>& output, bool verbose=false);
+ template<class T> void interpolateUp(const vector<double>& wavelengthIn, const vector< vector<T> >& input, const vector<double>& wavelengthOut, const std::string& type, vector< vector<T> >& output, bool verbose=false);
+ // template<class T> void interpolateUp(const vector< vector<T> >& input, vector< vector<T> >& output, double start, double end, double step, const gsl_interp_type* type);
+ // template<class T> void interpolateUp(const vector< vector<T> >& input, const vector<double>& wavelengthIn, vector< vector<T> >& output, vector<double>& wavelengthOut, double start, double end, double step, const gsl_interp_type* type);
template<class T> void interpolateUp(const vector<T>& input, vector<T>& output, int nbin);
template<class T> void interpolateUp(double* input, int dim, vector<T>& output, int nbin);
template<class T> void interpolateDown(const vector<T>& input, vector<T>& output, int nbin);
template<class T> void interpolateDown(double* input, int dim, vector<T>& output, int nbin);
private:
+ static void initMap(std::map<std::string, INTERPOLATION_TYPE>& m_interpMap){
+ //initialize selMap
+ m_interpMap["linear"]=linear;
+ m_interpMap["polynomial"]=polynomial;
+ m_interpMap["cspline"]=cspline;
+ m_interpMap["cspline_periodic"]=cspline_periodic;
+ m_interpMap["akima"]=akima;
+ m_interpMap["akima_periodic"]=akima_periodic;
+ }
+
};
+
template<class T> typename vector<T>::iterator StatFactory::max(const vector<T>& v, typename vector<T>::iterator begin, typename vector<T>::iterator end) const
{
typename vector<T>::iterator tmpIt=begin;
@@ -460,108 +519,67 @@ template<class T> double StatFactory::cross_correlation(const vector<T>& x, cons
//alternatively: use GNU scientific library:
// gsl_stats_correlation (const double data1[], const size_t stride1, const double data2[], const size_t stride2, const size_t n)
- template<class T> void StatFactory::interpolateUp(const vector< vector<T> >& input, const vector<double>& wavelengthIn, vector< vector<T> >& output, vector<double>& wavelengthOut, double start, double end, double step, int type){
- int nsample=input.size();//first sample contains wavelength
- int nband=wavelengthIn.size();
+template<class T> void StatFactory::interpolateUp(const vector<double>& wavelengthIn, const vector<T>& input, const vector<double>& wavelengthOut, const std::string& type, vector<T>& output, bool verbose){
+ assert(wavelengthIn.size());
+ assert(input.size()==wavelengthIn.size());
+ assert(wavelengthOut.size());
+ int nband=wavelengthIn.size();
output.clear();
- wavelengthOut.clear();
- output.resize(nsample);//first sample contains wavelength
- start=(start)?start:floor(wavelengthIn[0]);
- end=(end)?end:ceil(wavelengthIn.back());
- for(double xi=start;xi<=end;xi+=step)
- wavelengthOut.push_back(xi);
- for(int isample=0;isample<nsample;++isample){
- gsl_interp_accel *acc=gsl_interp_accel_alloc();
- gsl_spline *spline;
- switch(type){
- case(POLYNOMIAL):
- spline=gsl_spline_alloc(gsl_interp_polynomial,nband);
- break;
- case(CSPLINE):
- spline=gsl_spline_alloc(gsl_interp_cspline,nband);
- break;
- case(PERIODIC):
- spline=gsl_spline_alloc(gsl_interp_cspline_periodic,nband);
- break;
- case(AKIMA):
- spline=gsl_spline_alloc(gsl_interp_akima,nband);
- break;
- case(AKIMA_PERIODIC):
- spline=gsl_spline_alloc(gsl_interp_akima_periodic,nband);
- break;
- case(LINEAR):
- spline=gsl_spline_alloc(gsl_interp_linear,nband);
- break;
- case(UNDEFINED):
- default:{
- string errorString="Error: interpolation type not defined: ";
- errorString+=type;
- throw(errorString);
- break;
- }
- }
- gsl_spline_init(spline,&(wavelengthIn[0]),&(input[isample][0]),nband);
- for(double xi=start;xi<=end;xi+=step){
- if(!type&&xi>wavelengthIn.back())
- output[isample].push_back(output[isample].back());
- else
- output[isample].push_back(gsl_spline_eval(spline,xi,acc));
+ gsl_interp_accel *acc;
+ allocAcc(acc);
+ gsl_spline *spline;
+ getSpline(type,nband,spline);
+ assert(spline);
+ assert(&(wavelengthIn[0]));
+ assert(&(input[0]));
+ initSpline(spline,&(wavelengthIn[0]),&(input[0]),nband);
+ for(int index=0;index<wavelengthOut.size();++index){
+ if(type=="linear"){
+ if(wavelengthOut[index]<wavelengthIn.back()){
+ output.push_back(*(input.begin()));
+ continue;
+ }
+ else if(wavelengthOut[index]>wavelengthIn.back()){
+ output.push_back(input.back());
+ continue;
+ }
}
- gsl_spline_free(spline);
- gsl_interp_accel_free(acc);
- }
+ double dout=evalSpline(spline,wavelengthOut[index],acc);
+ output.push_back(dout);
}
+ gsl_spline_free(spline);
+ gsl_interp_accel_free(acc);
+}
- template<class T> void StatFactory::interpolateUp(const vector< vector<T> >& input, vector< vector<T> >& output, double start, double end, double step, int type)
-{
- int nsample=input.size()-1;//first sample contains wavelength
- int nband=input[0].size();
+template<class T> void StatFactory::interpolateUp(const vector<double>& wavelengthIn, const vector< vector<T> >& input, const vector<double>& wavelengthOut, const std::string& type, vector< vector<T> >& output, bool verbose){
+ assert(wavelengthIn.size());
+ assert(wavelengthOut.size());
+ int nsample=input.size();
+ int nband=wavelengthIn.size();
output.clear();
- output.resize(nsample+1);//first sample contains wavelength
- start=(start)?start:floor(input[0][0]);
- end=(end)?end:ceil(input[0].back());
- for(double xi=start;xi<=end;xi+=step)
- output[0].push_back(xi);
- for(int isample=1;isample<nsample+1;++isample){
- gsl_interp_accel *acc=gsl_interp_accel_alloc();
- gsl_spline *spline;
- switch(type){
- case(POLYNOMIAL):
- spline=gsl_spline_alloc(gsl_interp_polynomial,nband);
- break;
- case(CSPLINE):
- spline=gsl_spline_alloc(gsl_interp_cspline,nband);
- break;
- case(PERIODIC):
- spline=gsl_spline_alloc(gsl_interp_cspline_periodic,nband);
- break;
- case(AKIMA):
- spline=gsl_spline_alloc(gsl_interp_akima,nband);
- break;
- case(AKIMA_PERIODIC):
- spline=gsl_spline_alloc(gsl_interp_akima_periodic,nband);
- break;
- case(LINEAR):
- spline=gsl_spline_alloc(gsl_interp_linear,nband);
- break;
- case(UNDEFINED):
- default:{
- string errorString="Error: interpolation type not defined: ";
- errorString+=type;
- throw(errorString);
- break;
- }
- }
- gsl_spline_init(spline,&(input[0][0]),&(input[isample][0]),nband);
- for(double xi=start;xi<=end;xi+=step){
- if(!type&&xi>input[0].back())
- output[isample].push_back(output[isample].back());
- else
- output[isample].push_back(gsl_spline_eval(spline,xi,acc));
+ output.resize(nsample);
+ gsl_interp_accel *acc;
+ allocAcc(acc);
+ gsl_spline *spline;
+ getSpline(type,nband,spline);
+ for(int isample=0;isample<nsample;++isample){
+ assert(input[isample].size()==wavelengthIn.size());
+ initSpline(spline,&(wavelengthIn[0]),&(input[isample][0]),nband);
+ for(int index=0;index<wavelengthOut.size();++index){
+ if(type=="linear"){
+ if(wavelengthOut[index]<wavelengthIn.back())
+ output[isample].push_back(*(input.begin()));
+ else if(wavelengthOut[index]>wavelengthIn.back())
+ output[isample].push_back(input.back());
+ }
+ else{
+ double dout=evalSpline(spline,wavelengthOut[index],acc);
+ output.push_back(dout);
+ }
}
- gsl_spline_free(spline);
- gsl_interp_accel_free(acc);
}
+ gsl_spline_free(spline);
+ gsl_interp_accel_free(acc);
}
template<class T> void StatFactory::interpolateUp(const vector<T>& input, vector<T>& output, int nbin)
diff --git a/src/apps/pkcrop.cc b/src/apps/pkcrop.cc
index 158ed5f..54dbbec 100644
--- a/src/apps/pkcrop.cc
+++ b/src/apps/pkcrop.cc
@@ -39,8 +39,14 @@ int main(int argc, char *argv[])
Optionpk<double> uly_opt("uly", "uly", "Upper left y value bounding box (in geocoordinates if georef is true)", 0.0);
Optionpk<double> lrx_opt("lrx", "lrx", "Lower right x value bounding box (in geocoordinates if georef is true)", 0.0);
Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box (in geocoordinates if georef is true)", 0.0);
- Optionpk<double> dx_opt("dx", "dx", "Output resolution in x (in meter) (0.0: keep original resolution)", 0.0);
- Optionpk<double> dy_opt("dy", "dy", "Output resolution in y (in meter) (0.0: keep original resolution)", 0.0);
+ Optionpk<double> dx_opt("dx", "dx", "Output resolution in x (in meter) (0.0: keep original resolution)");
+ Optionpk<double> dy_opt("dy", "dy", "Output resolution in y (in meter) (0.0: keep original resolution)");
+ Optionpk<double> cx_opt("x", "x", "x-coordinate of image centre to crop (in meter)");
+ Optionpk<double> cy_opt("y", "y", "y-coordinate of image centre to crop (in meter)");
+ Optionpk<double> nx_opt("nx", "nx", "image size in x to crop (in meter)");
+ Optionpk<double> ny_opt("ny", "ny", "image size in y to crop (in meter)");
+ Optionpk<int> ns_opt("ns", "ns", "number of samples to crop (in pixels)");
+ Optionpk<int> nl_opt("nl", "nl", "number of lines to crop (in pixels)");
Optionpk<int> band_opt("b", "band", "band index to crop (-1: crop all bands)", -1);
Optionpk<double> scale_opt("s", "scale", "output=scale*input+offset", 1);
Optionpk<double> offset_opt("off", "offset", "output=scale*input+offset", 0);
@@ -72,6 +78,12 @@ int main(int argc, char *argv[])
colorTable_opt.retrieveOption(argc,argv);
dx_opt.retrieveOption(argc,argv);
dy_opt.retrieveOption(argc,argv);
+ cx_opt.retrieveOption(argc,argv);
+ cy_opt.retrieveOption(argc,argv);
+ nx_opt.retrieveOption(argc,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);
@@ -118,6 +130,27 @@ int main(int argc, char *argv[])
pfnProgress(progress,pszMessage,pProgressArg);
ImgReaderGdal imgReader;
ImgWriterGdal imgWriter;
+ //open input images to extract number of bands and spatial resolution
+ int ncropband=0;//total number of bands to write
+ double dx=(dx_opt.size())? dx_opt[0]:0;
+ double dy=(dy_opt.size())? dy_opt[0]:0;
+ for(int iimg=0;iimg<input_opt.size();++iimg){
+ imgReader.open(input_opt[iimg]);
+ if(dx_opt.empty()){
+ if(!iimg||imgReader.getDeltaX()<dx)
+ dx=imgReader.getDeltaX();
+ }
+ if(dy_opt.empty()){
+ if(!iimg||imgReader.getDeltaY()<dy)
+ dy=imgReader.getDeltaY();
+ }
+ if(band_opt[0]>=0)
+ ncropband+=band_opt.size();
+ else
+ ncropband+=imgReader.nrOfBand();
+ imgReader.close();
+ }
+
GDALDataType theType=GDT_Unknown;
if(verbose_opt[0])
cout << "possible output data types: ";
@@ -136,8 +169,6 @@ int main(int argc, char *argv[])
else
cout << "Output pixel type: " << GDALGetDataTypeName(theType) << endl;
}
- double dx=dx_opt[0];
- double dy=dy_opt[0];
//bounding box of cropped image
double cropulx=ulx_opt[0];
double cropuly=uly_opt[0];
@@ -152,23 +183,34 @@ int main(int argc, char *argv[])
cerr << "Error: could not get extent from " << extent_opt[0] << endl;
exit(1);
}
- if(ulx_opt[0]<cropulx)
- cropulx=ulx_opt[0];
- if(uly_opt[0]>cropuly)
- cropuly=uly_opt[0];
- if(lry_opt[0]<croplry)
- croplry=lry_opt[0];
- if(lrx_opt[0]>croplrx)
- croplrx=lrx_opt[0];
extentReader.close();
}
if(mask_opt[0])
extentReader.open(extent_opt[0]);
}
+ else if(cx_opt.size()&&cy_opt.size()&&nx_opt.size()&&ny_opt.size()){
+ ulx_opt[0]=cx_opt[0]-nx_opt[0]/2.0;
+ uly_opt[0]=cy_opt[0]+ny_opt[0]/2.0;
+ lrx_opt[0]=cx_opt[0]+nx_opt[0]/2.0;
+ lry_opt[0]=cy_opt[0]-ny_opt[0]/2.0;
+ }
+ else if(cx_opt.size()&&cy_opt.size()&&ns_opt.size()&&nl_opt.size()){
+ ulx_opt[0]=cx_opt[0]-ns_opt[0]*dx/2.0;
+ uly_opt[0]=cy_opt[0]+nl_opt[0]*dy/2.0;
+ lrx_opt[0]=cx_opt[0]+ns_opt[0]*dx/2.0;
+ lry_opt[0]=cy_opt[0]-nl_opt[0]*dy/2.0;
+ }
+ if(ulx_opt[0]<cropulx)
+ cropulx=ulx_opt[0];
+ if(uly_opt[0]>cropuly)
+ cropuly=uly_opt[0];
+ if(lry_opt[0]<croplry)
+ croplry=lry_opt[0];
+ if(lrx_opt[0]>croplrx)
+ croplrx=lrx_opt[0];
if(verbose_opt[0])
cout << "--ulx=" << ulx_opt[0] << " --uly=" << uly_opt[0] << " --lrx=" << lrx_opt[0] << " --lry=" << lry_opt[0] << endl;
//determine number of output bands
- int ncropband=0;//total number of bands to write
int writeBand=0;//write band
while(scale_opt.size()<band_opt.size())
@@ -177,14 +219,6 @@ int main(int argc, char *argv[])
offset_opt.push_back(offset_opt[0]);
for(int iimg=0;iimg<input_opt.size();++iimg){
- imgReader.open(input_opt[iimg]);
- if(band_opt[0]>=0)
- ncropband+=band_opt.size();
- else
- ncropband+=imgReader.nrOfBand();
- imgReader.close();
- }
- for(int iimg=0;iimg<input_opt.size();++iimg){
if(verbose_opt[0])
cout << "opening image " << input_opt[iimg] << endl;
imgReader.open(input_opt[iimg]);
@@ -203,10 +237,10 @@ int main(int argc, char *argv[])
int ncol=imgReader.nrOfCol();
int ncropcol=0;
int ncroprow=0;
- if(!dx||!dy){
- dx=imgReader.getDeltaX();
- dy=imgReader.getDeltaY();
- }
+ // if(!dx||!dy){
+ // dx=imgReader.getDeltaX();
+ // dy=imgReader.getDeltaY();
+ // }
if(verbose_opt[0])
cout << "size of " << input_opt[iimg] << ": " << ncol << " cols, "<< nrow << " rows" << endl;
double uli,ulj,lri,lrj;//image coordinates
diff --git a/src/apps/pkdumpimg.cc b/src/apps/pkdumpimg.cc
index ea7e52a..50a88e1 100644
--- a/src/apps/pkdumpimg.cc
+++ b/src/apps/pkdumpimg.cc
@@ -36,7 +36,7 @@ int main(int argc, char *argv[])
Optionpk<string> output_opt("o", "output", "Output ascii file (Default is empty: use stdout", "");
Optionpk<string> otype_opt("ot", "otype", "Data type for output ({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 format (matrix form or list (x,y,z) form. Default is matrix form", "matrix");
- Optionpk<int> band_opt("b", "band", "band index to crop (-1: crop all bands)", 0);
+ Optionpk<int> band_opt("b", "band", "band index to crop");
Optionpk<string> extent_opt("e", "extent", "get boundary from extent from polygons in vector file", "");
Optionpk<double> ulx_opt("ulx", "ulx", "Upper left x value bounding box (in geocoordinates if georef is true)",0.0);
Optionpk<double> uly_opt("uly", "uly", "Upper left y value bounding box (in geocoordinates if georef is true)",0.0);
@@ -213,105 +213,111 @@ int main(int argc, char *argv[])
//test
// vector<double> readBuffer(readncol);
vector<double> readBuffer(readncol+1);
- assert(band_opt[0]>=0);
- assert(band_opt[0]<imgReader.nrOfBand());
assert(imgWriter.isGeoRef());
- for(int irow=0;irow<ncroprow;++irow){
- if(verbose_opt[0])
- std::cout << irow << std::endl;
- double x=0;
- double y=0;
- //convert irow to geo
- imgWriter.image2geo(0,irow,x,y);
- //lookup corresponding row for irow in this file
- imgReader.geo2image(x,y,readCol,readRow);
- if(readRow<0||readRow>=imgReader.nrOfRow()){
- if(verbose_opt[0]>1)
- std::cout << "skipping row " << readRow << std::endl;
- continue;
- }
- try{
- //test
- // imgReader.readData(readBuffer,GDT_Float64,startCol,endCol,readRow,band_opt[0],theResample);
- if(endCol<imgReader.nrOfCol()-1)
- imgReader.readData(readBuffer,GDT_Float64,startCol,endCol+1,readRow,band_opt[0],theResample);
- else
- imgReader.readData(readBuffer,GDT_Float64,startCol,endCol,readRow,band_opt[0],theResample);
- for(int ib=0;ib<ncropcol;++ib){
- assert(imgWriter.image2geo(ib,irow,x,y));
- //lookup corresponding row for irow in this file
- imgReader.geo2image(x,y,readCol,readRow);
- if(readCol<0||readCol>=imgReader.nrOfCol()){
- if(oformat_opt[0]=="matrix"){
- if(output_opt[0].empty())
- std::cout << flag_opt[0] << " ";
- else
- outputStream << flag_opt[0] << " ";
- }
- else{
- if(output_opt[0].empty())
- std::cout << x << " " << y << " " << flag_opt[0] << endl;
- else
- outputStream << x << " " << y << " " << flag_opt[0] << endl;
- }
- }
- else{
- switch(theResample){
- case(BILINEAR):
- lowerCol=readCol-0.5;
- lowerCol=static_cast<int>(lowerCol);
- upperCol=readCol+0.5;
- upperCol=static_cast<int>(upperCol);
- if(lowerCol<0)
- lowerCol=0;
- if(upperCol>=imgReader.nrOfCol())
- upperCol=imgReader.nrOfCol()-1;
+ if(band_opt.empty()){
+ for(int iband=0;iband<imgReader.nrOfBand();++iband)
+ band_opt.push_back(iband);
+ }
+ for(int iband=0;iband<band_opt.size();++iband){
+ assert(band_opt[iband]>=0);
+ assert(band_opt[iband]<imgReader.nrOfBand());
+ for(int irow=0;irow<ncroprow;++irow){
+ if(verbose_opt[0])
+ std::cout << irow << std::endl;
+ double x=0;
+ double y=0;
+ //convert irow to geo
+ imgWriter.image2geo(0,irow,x,y);
+ //lookup corresponding row for irow in this file
+ imgReader.geo2image(x,y,readCol,readRow);
+ if(readRow<0||readRow>=imgReader.nrOfRow()){
+ if(verbose_opt[0]>1)
+ std::cout << "skipping row " << readRow << std::endl;
+ continue;
+ }
+ try{
+ //test
+ // imgReader.readData(readBuffer,GDT_Float64,startCol,endCol,readRow,band_opt[0],theResample);
+ if(endCol<imgReader.nrOfCol()-1)
+ imgReader.readData(readBuffer,GDT_Float64,startCol,endCol+1,readRow,band_opt[iband],theResample);
+ else
+ imgReader.readData(readBuffer,GDT_Float64,startCol,endCol,readRow,band_opt[iband],theResample);
+ for(int ib=0;ib<ncropcol;++ib){
+ assert(imgWriter.image2geo(ib,irow,x,y));
+ //lookup corresponding row for irow in this file
+ imgReader.geo2image(x,y,readCol,readRow);
+ if(readCol<0||readCol>=imgReader.nrOfCol()){
if(oformat_opt[0]=="matrix"){
if(output_opt[0].empty())
- std::cout << (readCol-0.5-lowerCol)*readBuffer[upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[lowerCol-startCol] << " ";
+ std::cout << flag_opt[0] << " ";
else
- outputStream << (readCol-0.5-lowerCol)*readBuffer[upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[lowerCol-startCol] << " ";
+ outputStream << flag_opt[0] << " ";
}
- else if(!imgReader.isNoData(readBuffer[upperCol-startCol])&&!imgReader.isNoData(readBuffer[lowerCol-startCol])){
+ else{
+ if(output_opt[0].empty())
+ std::cout << x << " " << y << " " << flag_opt[0] << endl;
+ else
+ outputStream << x << " " << y << " " << flag_opt[0] << endl;
+ }
+ }
+ else{
+ switch(theResample){
+ case(BILINEAR):
+ lowerCol=readCol-0.5;
+ lowerCol=static_cast<int>(lowerCol);
+ upperCol=readCol+0.5;
+ upperCol=static_cast<int>(upperCol);
+ if(lowerCol<0)
+ lowerCol=0;
+ if(upperCol>=imgReader.nrOfCol())
+ upperCol=imgReader.nrOfCol()-1;
+ if(oformat_opt[0]=="matrix"){
+ if(output_opt[0].empty())
+ std::cout << (readCol-0.5-lowerCol)*readBuffer[upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[lowerCol-startCol] << " ";
+ else
+ outputStream << (readCol-0.5-lowerCol)*readBuffer[upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[lowerCol-startCol] << " ";
+ }
+ else if(!imgReader.isNoData(readBuffer[upperCol-startCol])&&!imgReader.isNoData(readBuffer[lowerCol-startCol])){
if(output_opt[0].empty())
std::cout << x << " " << y << " " << (readCol-0.5-lowerCol)*readBuffer[upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[lowerCol-startCol] << " " << endl;
else
outputStream << x << " " << y << " " << (readCol-0.5-lowerCol)*readBuffer[upperCol-startCol]+(1-readCol+0.5+lowerCol)*readBuffer[lowerCol-startCol] << " " << endl;
+ }
+ break;
+ default:
+ readCol=static_cast<int>(readCol);
+ readCol-=startCol;//we only start reading from startCol
+ assert(readCol>=0);
+ if(readCol>=readBuffer.size())
+ std::cout << "Error: " << readCol << " >= " << readBuffer.size() << std::endl;
+ assert(readCol<readBuffer.size());
+ if(oformat_opt[0]=="matrix"){
+ if(output_opt[0].empty())
+ std::cout << readBuffer[readCol] << " ";
+ else
+ outputStream << readBuffer[readCol] << " ";
+ }
+ else if(!imgReader.isNoData(readBuffer[readCol])){
+ if(output_opt[0].empty())
+ std::cout << x << " " << y << " " << readBuffer[readCol] << std::endl;
+ else
+ outputStream << x << " " << y << " " << readBuffer[readCol] << std::endl;
+ }
+ break;
}
- break;
- default:
- readCol=static_cast<int>(readCol);
- readCol-=startCol;//we only start reading from startCol
- assert(readCol>=0);
- if(readCol>=readBuffer.size())
- std::cout << "Error: " << readCol << " >= " << readBuffer.size() << std::endl;
- assert(readCol<readBuffer.size());
- if(oformat_opt[0]=="matrix"){
- if(output_opt[0].empty())
- std::cout << readBuffer[readCol] << " ";
- else
- outputStream << readBuffer[readCol] << " ";
- }
- else if(!imgReader.isNoData(readBuffer[readCol])){
- if(output_opt[0].empty())
- std::cout << x << " " << y << " " << readBuffer[readCol] << std::endl;
- else
- outputStream << x << " " << y << " " << readBuffer[readCol] << std::endl;
- }
- break;
}
}
}
- }
- catch(string errorstring){
- cout << errorstring << endl;
- exit(1);
- }
- if(oformat_opt[0]=="matrix"){
- if(output_opt[0].empty())
- std::cout << std::endl;
- else
- outputStream << std::endl;
+ catch(string errorstring){
+ cout << errorstring << endl;
+ exit(1);
+ }
+ if(oformat_opt[0]=="matrix"){
+ if(output_opt[0].empty())
+ std::cout << std::endl;
+ else
+ outputStream << std::endl;
+ }
}
}
if(extent_opt[0]!=""){
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index bb20580..0ee8774 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -40,7 +40,7 @@ int main(int argc,char **argv) {
Optionpk<std::string> output_opt("o", "output", "Output image file");
Optionpk<bool> disc_opt("c", "circular", "circular disc kernel for dilation and erosion", false);
Optionpk<double> angle_opt("a", "angle", "angle used for directional filtering in dilation.");
- Optionpk<int> function_opt("f", "filter", "filter function (0: median, 1: variance, 2: min, 3: max, 4: sum, 5: mean, 6: minmax, 7: dilation, 8: erosion, 9: closing, 10: opening, 11: spatially homogeneous (central pixel must be identical to all other pixels within window), 12: SobelX edge detection in X, 13: SobelY edge detection in Y, 14: SobelXY, -14: SobelYX, 15: smooth, 16: density, 17: majority voting (only for classes), 18: forest aggregation (mixed), 19: smooth no data (mask) val [...]
+ Optionpk<std::string> method_opt("f", "filter", "filter function (median,variance,min,max,sum,mean,minmax,dilation,erosion,closing,opening,spatially homogeneous (central pixel must be identical to all other pixels within window),SobelX edge detection in X,SobelY edge detection in Y,SobelXY,SobelYX,smooth,density,majority voting (only for classes),forest aggregation (mixed),smooth no data (mask) values,threshold local filtering,ismin,ismax,heterogeneous (central pixel must be different [...]
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");
@@ -52,6 +52,7 @@ int main(int argc,char **argv) {
Optionpk<double> fwhm_opt("fwhm", "fwhm", "list of full width half to apply spectral filtering (-fwhm band1 -fwhm band2 ...)");
Optionpk<double> wavelengthIn_opt("win", "wavelengthIn", "list of wavelengths in input spectrum (-w band1 -w band2 ...)");
Optionpk<double> wavelengthOut_opt("wout", "wavelengthOut", "list of wavelengths in output spectrum (-w band1 -w band2 ...)");
+ Optionpk<std::string> interpolationType_opt("interp", "interp", "type of interpolation for spectral filtering (see http://www.gnu.org/software/gsl/manual/html_node/Interpolation-Types.html)","linear");
Optionpk<std::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)");
@@ -65,7 +66,7 @@ int main(int argc,char **argv) {
output_opt.retrieveOption(argc,argv);
disc_opt.retrieveOption(argc,argv);
angle_opt.retrieveOption(argc,argv);
- function_opt.retrieveOption(argc,argv);
+ method_opt.retrieveOption(argc,argv);
dimX_opt.retrieveOption(argc,argv);
dimY_opt.retrieveOption(argc,argv);
dimZ_opt.retrieveOption(argc,argv);
@@ -79,6 +80,7 @@ int main(int argc,char **argv) {
wavelengthIn_opt.retrieveOption(argc,argv);
wavelengthOut_opt.retrieveOption(argc,argv);
down_opt.retrieveOption(argc,argv);
+ interpolationType_opt.retrieveOption(argc,argv);
otype_opt.retrieveOption(argc,argv);
oformat_opt.retrieveOption(argc,argv);
colorTable_opt.retrieveOption(argc,argv);
@@ -126,7 +128,11 @@ int main(int argc,char **argv) {
}
try{
assert(output_opt.size());
- output.open(output_opt[0],(input.nrOfCol()+down_opt[0]-1)/down_opt[0],(input.nrOfRow()+down_opt[0]-1)/down_opt[0],input.nrOfBand(),theType,imageType,option_opt);
+ if(wavelengthOut_opt.size())
+ //todo: support down and offset
+ output.open(output_opt[0],input.nrOfCol(),input.nrOfRow(),wavelengthOut_opt.size(),theType,imageType,option_opt);
+ else
+ output.open(output_opt[0],(input.nrOfCol()+down_opt[0]-1)/down_opt[0],(input.nrOfRow()+down_opt[0]-1)/down_opt[0],input.nrOfBand(),theType,imageType,option_opt);
}
catch(string errorstring){
cout << errorstring << endl;
@@ -154,27 +160,29 @@ int main(int argc,char **argv) {
filter2d::Filter2d filter2d;
filter::Filter filter1d;
- if(class_opt.size()&&verbose_opt[0])
+ if(class_opt.size()&&verbose_opt[0]){
std::cout<< "class values: ";
- for(int iclass=0;iclass<class_opt.size();++iclass){
- if(!dimZ_opt.size())
- filter2d.pushClass(class_opt[iclass]);
- else
- filter1d.pushClass(class_opt[iclass]);
+ for(int iclass=0;iclass<class_opt.size();++iclass){
+ if(!dimZ_opt.size())
+ filter2d.pushClass(class_opt[iclass]);
+ else
+ filter1d.pushClass(class_opt[iclass]);
+ if(verbose_opt[0])
+ std::cout<< class_opt[iclass] << " ";
+ }
if(verbose_opt[0])
- std::cout<< class_opt[iclass] << " ";
+ std::cout<< std::endl;
}
- if(verbose_opt[0])
- std::cout<< std::endl;
- if(verbose_opt[0])
+ if(mask_opt.size()&&verbose_opt[0]){
std::cout<< "mask values: ";
- for(int imask=0;imask<mask_opt.size();++imask){
+ for(int imask=0;imask<mask_opt.size();++imask){
+ if(verbose_opt[0])
+ std::cout<< mask_opt[imask] << " ";
+ filter2d.pushMask(mask_opt[imask]);
+ }
if(verbose_opt[0])
- std::cout<< mask_opt[imask] << " ";
- filter2d.pushMask(mask_opt[imask]);
+ std::cout<< std::endl;
}
- if(verbose_opt[0])
- std::cout<< std::endl;
if(tap_opt.size()){
ifstream tapfile(tap_opt[0].c_str());
assert(tapfile);
@@ -203,95 +211,56 @@ int main(int argc,char **argv) {
filter1d.doit(input,output,down_opt[0]);
}
else if(wavelengthOut_opt.size()){
+ if(verbose_opt[0])
+ std::cout << "spectral filtering to " << wavelengthOut_opt.size() << " bands with provided fwhm " << std::endl;
assert(wavelengthOut_opt.size()==fwhm_opt.size());
assert(wavelengthIn_opt.size());
- filter1d.applyFwhm(input,output,wavelengthIn_opt,wavelengthOut_opt,fwhm_opt,verbose_opt[0]);
+ filter1d.applyFwhm(wavelengthIn_opt,input,wavelengthOut_opt,fwhm_opt,interpolationType_opt[0],output,verbose_opt[0]);
}
else{
if(colorTable_opt.size())
output.setColorTable(colorTable_opt[0]);
- switch(function_opt[0]){
- case(filter2d::MEDIAN):
- filter2d.doit(input,output,filter2d::MEDIAN,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- case(filter2d::VAR):
- filter2d.doit(input,output,filter2d::VAR,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- case(filter2d::STDEV):
- filter2d.doit(input,output,filter2d::STDEV,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- case(filter2d::MIN):
- filter2d.doit(input,output,filter2d::MIN,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- case(filter2d::ISMIN):
- filter2d.doit(input,output,filter2d::ISMIN,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- case(filter2d::MAX):
- filter2d.doit(input,output,filter2d::MAX,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- case(filter2d::ISMAX):
- filter2d.doit(input,output,filter2d::ISMAX,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- case(filter2d::MINMAX):
- filter2d.doit(input,output,filter2d::MINMAX,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- case(filter2d::SUM):
- filter2d.doit(input,output,filter2d::SUM,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- case(filter2d::MEAN):
- filter2d.doit(input,output,filter2d::MEAN,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- case(filter2d::MAJORITY):
- filter2d.doit(input,output,filter2d::MAJORITY,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- case(filter2d::THRESHOLD):
- filter2d.setThresholds(threshold_opt);
- filter2d.setClasses(class_opt);
- filter2d.doit(input,output,filter2d::THRESHOLD,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- case(filter2d::MIXED):
- filter2d.doit(input,output,filter2d::MIXED,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- case(filter2d::DILATE):
+ switch(filter2d::Filter2d::getFilterType(method_opt[0])){
+ case(filter2d::dilate):
if(dimZ_opt.size()){
if(verbose_opt[0])
std::cout<< "1-D filtering: dilate" << std::endl;
- filter1d.morphology(input,output,filter::Filter::DILATE,dimZ_opt[0]);
+ filter1d.morphology(input,output,"dilate",dimZ_opt[0]);
}
else{
if(angle_opt.size())
- filter2d.morphology(input,output,filter2d::DILATE,dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
+ filter2d.morphology(input,output,"dilate",dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
else
- filter2d.morphology(input,output,filter2d::DILATE,dimX_opt[0],dimY_opt[0],disc_opt[0]);
+ filter2d.morphology(input,output,"dilate",dimX_opt[0],dimY_opt[0],disc_opt[0]);
}
break;
- case(filter2d::ERODE):
+ case(filter2d::erode):
if(dimZ_opt[0]>0){
if(verbose_opt[0])
std::cout<< "1-D filtering: dilate" << std::endl;
- filter1d.morphology(input,output,filter::Filter::ERODE,dimZ_opt[0]);
+ filter1d.morphology(input,output,"erode",dimZ_opt[0]);
}
else{
if(angle_opt.size())
- filter2d.morphology(input,output,filter2d::ERODE,dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
+ filter2d.morphology(input,output,"erode",dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
else
- filter2d.morphology(input,output,filter2d::ERODE,dimX_opt[0],dimY_opt[0],disc_opt[0]);
+ filter2d.morphology(input,output,"erode",dimX_opt[0],dimY_opt[0],disc_opt[0]);
}
break;
- case(filter2d::CLOSE):{//closing
+ case(filter2d::close):{//closing
ostringstream tmps;
tmps << "/tmp/dilation_" << getpid() << ".tif";
ImgWriterGdal tmpout;
tmpout.open(tmps.str(),input);
try{
if(dimZ_opt.size()){
- filter1d.morphology(input,tmpout,filter::Filter::DILATE,dimZ_opt[0]);
+ filter1d.morphology(input,tmpout,"dilate",dimZ_opt[0]);
}
else{
if(angle_opt.size())
- filter2d.morphology(input,tmpout,filter2d::DILATE,dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
+ filter2d.morphology(input,tmpout,"dilate",dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
else
- filter2d.morphology(input,tmpout,filter2d::DILATE,dimX_opt[0],dimY_opt[0],disc_opt[0]);
+ filter2d.morphology(input,tmpout,"dilate",dimX_opt[0],dimY_opt[0],disc_opt[0]);
}
}
catch(std::string errorString){
@@ -302,13 +271,13 @@ int main(int argc,char **argv) {
ImgReaderGdal tmpin;
tmpin.open(tmps.str());
if(dimZ_opt.size()){
- filter1d.morphology(tmpin,output,filter::Filter::ERODE,dimZ_opt[0]);
+ filter1d.morphology(tmpin,output,"erode",dimZ_opt[0]);
}
else{
if(angle_opt.size())
- filter2d.morphology(tmpin,output,filter2d::ERODE,dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
+ filter2d.morphology(tmpin,output,"erode",dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
else
- filter2d.morphology(tmpin,output,filter2d::ERODE,dimX_opt[0],dimY_opt[0],disc_opt[0]);
+ filter2d.morphology(tmpin,output,"erode",dimX_opt[0],dimY_opt[0],disc_opt[0]);
}
tmpin.close();
if(remove(tmps.str().c_str( )) !=0){
@@ -316,31 +285,31 @@ int main(int argc,char **argv) {
}
break;
}
- case(filter2d::OPEN):{//opening
+ case(filter2d::open):{//opening
ostringstream tmps;
tmps << "/tmp/erosion_" << getpid() << ".tif";
ImgWriterGdal tmpout;
tmpout.open(tmps.str(),input);
if(dimZ_opt.size()){
- filter1d.morphology(input,tmpout,filter::Filter::ERODE,dimZ_opt[0]);
+ filter1d.morphology(input,tmpout,"erode",dimZ_opt[0]);
}
else{
if(angle_opt.size())
- filter2d.morphology(input,tmpout,filter2d::ERODE,dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
+ filter2d.morphology(input,tmpout,"erode",dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
else
- filter2d.morphology(input,tmpout,filter2d::ERODE,dimX_opt[0],dimY_opt[0],disc_opt[0]);
+ filter2d.morphology(input,tmpout,"erode",dimX_opt[0],dimY_opt[0],disc_opt[0]);
}
tmpout.close();
ImgReaderGdal tmpin;
tmpin.open(tmps.str());
if(dimZ_opt.size()){
- filter1d.morphology(tmpin,output,filter::Filter::DILATE,dimZ_opt[0]);
+ filter1d.morphology(tmpin,output,"dilate",dimZ_opt[0]);
}
else{
if(angle_opt.size())
- filter2d.morphology(tmpin,output,filter2d::DILATE,dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
+ filter2d.morphology(tmpin,output,"dilate",dimX_opt[0],dimY_opt[0],disc_opt[0],angle_opt[0]);
else
- filter2d.morphology(tmpin,output,filter2d::DILATE,dimX_opt[0],dimY_opt[0],disc_opt[0]);
+ filter2d.morphology(tmpin,output,"dilate",dimX_opt[0],dimY_opt[0],disc_opt[0]);
}
tmpin.close();
if(remove(tmps.str().c_str( )) !=0){
@@ -348,16 +317,16 @@ int main(int argc,char **argv) {
}
break;
}
- case(filter2d::HOMOG):{//spatially homogeneous
- filter2d.doit(input,output,filter2d::HOMOG,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
+ case(filter2d::homog):{//spatially homogeneous
+ filter2d.doit(input,output,"homog",dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
// filter2d.homogeneousSpatial(input,output,dimX_opt[0],disc_opt[0]);
break;
}
- case(filter2d::HETEROG):{//spatially heterogeneous
- filter2d.doit(input,output,filter2d::HETEROG,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
+ case(filter2d::heterog):{//spatially heterogeneous
+ filter2d.doit(input,output,"heterog",dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
break;
}
- case(filter2d::SOBELX):{//Sobel edge detection in X
+ case(filter2d::sobelx):{//Sobel edge detection in X
Vector2d<double> theTaps(3,3);
theTaps[0][0]=-1.0;
theTaps[0][1]=0.0;
@@ -372,7 +341,7 @@ int main(int argc,char **argv) {
filter2d.filter(input,output,true);
break;
}
- case(filter2d::SOBELY):{//Sobel edge detection in Y
+ case(filter2d::sobely):{//Sobel edge detection in Y
Vector2d<double> theTaps(3,3);
theTaps[0][0]=1.0;
theTaps[0][1]=2.0;
@@ -387,7 +356,7 @@ int main(int argc,char **argv) {
filter2d.filter(input,output,true);
break;
}
- case(filter2d::SOBELXY):{//Sobel edge detection in XY
+ case(filter2d::sobelxy):{//Sobel edge detection in XY
Vector2d<double> theTaps(3,3);
theTaps[0][0]=0.0;
theTaps[0][1]=1.0;
@@ -402,7 +371,7 @@ int main(int argc,char **argv) {
filter2d.filter(input,output,true);
break;
}
- case(filter2d::SOBELYX):{//Sobel edge detection in XY
+ case(filter2d::sobelyx):{//Sobel edge detection in XY
Vector2d<double> theTaps(3,3);
theTaps[0][0]=2.0;
theTaps[0][1]=1.0;
@@ -417,26 +386,20 @@ int main(int argc,char **argv) {
filter2d.filter(input,output,true);
break;
}
- case(filter2d::SMOOTH):{//Smoothing filter
+ case(filter2d::smooth):{//Smoothing filter
filter2d.smooth(input,output,dimX_opt[0],dimY_opt[0]);
break;
}
- case(filter2d::SMOOTHNODATA):{//Smoothing filter
+ case(filter2d::smoothnodata):{//Smoothing filter
filter2d.smoothNoData(input,output,dimX_opt[0],dimY_opt[0]);
break;
}
- case(filter2d::DENSITY):{//estimation of forest density
- filter2d.doit(input,output,filter2d::DENSITY,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- }
- case(filter2d::ORDER):{//order centre pixel with respect to values in window
- assert(dimX_opt[0]);
- filter2d.doit(input,output,filter2d::ORDER,dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
- break;
- }
+ case(filter2d::threshold):
+ filter2d.setThresholds(threshold_opt);
+ filter2d.setClasses(class_opt);//deliberate fall through
default:
- cerr << "filter function " << function_opt[0] << " not supported" << std::endl;
- return 1;
+ filter2d.doit(input,output,method_opt[0],dimX_opt[0],dimY_opt[0],down_opt[0],disc_opt[0]);
+ break;
}
}
input.close();
diff --git a/src/apps/pkinfo.cc b/src/apps/pkinfo.cc
index 7776a10..2a89043 100644
--- a/src/apps/pkinfo.cc
+++ b/src/apps/pkinfo.cc
@@ -31,8 +31,8 @@ int main(int argc, char *argv[])
Optionpk<bool> bbox_te_opt("te", "te", "Shows bounding box in GDAL format: xmin ymin xmax ymax ", false,0);
Optionpk<bool> centre_opt("c", "centre", "Image centre in projected X,Y coordinates ", false,0);
Optionpk<bool> colorTable_opt("ct", "colourtable", "Shows colour table ", false,0);
- Optionpk<bool> samples_opt("s", "samples", "Number of samples in image ", false,0);
- Optionpk<bool> lines_opt("l", "lines", "Number of lines in image ", false,0);
+ Optionpk<bool> samples_opt("ns", "ns", "Number of samples in image ", false,0);
+ Optionpk<bool> lines_opt("nl", "nl", "Number of lines in image ", false,0);
Optionpk<bool> nband_opt("nb", "nband", "Show number of bands in image", false,0);
Optionpk<short> band_opt("b", "band", "Band specific information", 0,0);
Optionpk<bool> dx_opt("dx", "dx", "Gets resolution in x (in m)", false,0);
@@ -228,9 +228,9 @@ int main(int argc, char *argv[])
std::cout << "--ct none ";
}
if(samples_opt[0])
- std::cout << "--samples " << imgReader.nrOfCol() << " ";
+ std::cout << "--ns " << imgReader.nrOfCol() << " ";
if(lines_opt[0])
- std::cout << "--rows " << imgReader.nrOfRow() << " ";
+ std::cout << "--nl " << imgReader.nrOfRow() << " ";
if(nband_opt[0])
std::cout << "--nband " << imgReader.nrOfBand() << " ";
double minValue=0;
diff --git a/src/apps/pklas2img.cc b/src/apps/pklas2img.cc
index a19cd97..fecb07a 100644
--- a/src/apps/pklas2img.cc
+++ b/src/apps/pklas2img.cc
@@ -412,7 +412,7 @@ int main(int argc,char **argv) {
while(nchange&&iteration<maxIter_opt[0]){
double hThreshold=maxSlope_opt[0]*dimx;
Vector2d<float> newOutput;
- nchange=morphFilter.morphology(currentOutput,newOutput,filter2d::ERODE,dimx,dimy,disc_opt[0],hThreshold);
+ nchange=morphFilter.morphology(currentOutput,newOutput,"erode",dimx,dimy,disc_opt[0],hThreshold);
currentOutput=newOutput;
dimx+=2;//change from theory: originally double cellCize
dimy+=2;//change from theory: originally double cellCize
@@ -439,10 +439,10 @@ int main(int argc,char **argv) {
std::cout << "iteration " << iteration << " with window size " << dimx << " and dh_max: " << hThreshold << std::endl;
Vector2d<float> newOutput;
try{
- theFilter.morphology(outputData,currentOutput,filter2d::ERODE,dimx,dimy,disc_opt[0],maxSlope_opt[0]);
- theFilter.morphology(currentOutput,outputData,filter2d::DILATE,dimx,dimy,disc_opt[0],maxSlope_opt[0]);
+ theFilter.morphology(outputData,currentOutput,"erode",dimx,dimy,disc_opt[0],maxSlope_opt[0]);
+ theFilter.morphology(currentOutput,outputData,"dilate",dimx,dimy,disc_opt[0],maxSlope_opt[0]);
if(postFilter_opt[0]=="bunting"){//todo: implement doit in Filter2d on Vector2d
- theFilter.doit(outputData,currentOutput,filter2d::MEDIAN,dimx,dimy,1,disc_opt[0]);
+ theFilter.doit(outputData,currentOutput,"median",dimx,dimy,1,disc_opt[0]);
outputData=currentOutput;
}
}
@@ -468,8 +468,8 @@ int main(int argc,char **argv) {
morphFilter.setNoValue(0);
Vector2d<float> filterInput=outputData;
try{
- morphFilter.morphology(outputData,filterInput,filter2d::ERODE,dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
- morphFilter.morphology(filterInput,outputData,filter2d::DILATE,dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
+ morphFilter.morphology(outputData,filterInput,"erode",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
+ morphFilter.morphology(filterInput,outputData,"dilate",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
}
catch(std::string errorString){
cout << errorString << endl;
@@ -483,8 +483,8 @@ int main(int argc,char **argv) {
morphFilter.setNoValue(0);
Vector2d<float> filterInput=outputData;
try{
- morphFilter.morphology(outputData,filterInput,filter2d::DILATE,dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
- morphFilter.morphology(filterInput,outputData,filter2d::ERODE,dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
+ morphFilter.morphology(outputData,filterInput,"dilate",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
+ morphFilter.morphology(filterInput,outputData,"erode",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
}
catch(std::string errorString){
cout << errorString << endl;
diff --git a/src/fileclasses/FileReaderAscii.h b/src/fileclasses/FileReaderAscii.h
index 384f558..c106aca 100644
--- a/src/fileclasses/FileReaderAscii.h
+++ b/src/fileclasses/FileReaderAscii.h
@@ -32,6 +32,7 @@ public:
FileReaderAscii(const std::string& filename);
FileReaderAscii(const std::string& filename, const char& fieldseparator);
~FileReaderAscii(void);
+ void reset(){m_ifstream.clear();m_ifstream.seekg(0,ios::beg);};
void open(const std::string& filename);
void close(void);
void setFieldSeparator(const char& fieldseparator){m_fs=fieldseparator;};
@@ -39,6 +40,7 @@ public:
void setMaxRow(int maxRow){m_maxRow=maxRow;};
void setComment(char comment){m_comment=comment;};
template<class T> unsigned int readData(vector<vector<T> > &dataVector, const vector<int> &cols);
+ template<class T> unsigned int readData(vector<T> &dataVector, int col);
protected:
std::string m_filename;
std::ifstream m_ifstream;
@@ -79,6 +81,112 @@ void FileReaderAscii::close(){
// m_ifstream.clear();
}
+template<class T> unsigned int FileReaderAscii::readData(vector<T> &dataVector, int col){
+ bool verbose=false;
+ dataVector.clear();
+ int nrow=0;
+ bool withinRange=true;
+ if(m_fs>' '&&m_fs<='~'){//field separator is a regular character (minimum ASCII code is space, maximum ASCII code is tilde)
+ if(verbose)
+ cout << "reading csv file " << m_filename << endl;
+ string csvRecord;
+ while(getline(m_ifstream,csvRecord)){//read a line
+ withinRange=true;
+ if(nrow<m_minRow)
+ withinRange=false;
+ if(m_maxRow>m_minRow)
+ if(nrow>m_maxRow)
+ withinRange=false;
+ if(withinRange){
+ istringstream csvstream(csvRecord);
+ string item;
+ int ncol=0;
+ bool isComment=false;
+ while(getline(csvstream,item,m_fs)){//read a column
+ if(verbose)
+ cout << item << " ";
+ unsigned pos=item.find(m_comment);
+ if(pos!=std::string::npos){
+ if(pos>0)
+ item=item.substr(0,pos-1);
+ else
+ break;
+ if(verbose)
+ std::cout << "comment found, string is " << item << std::endl;
+ isComment=true;
+ }
+ if(ncol==col){
+ T value=string2type<T>(item);
+ if((value>=m_min&&value<=m_max)||m_max<=m_min)
+ dataVector.push_back(value);
+ }
+ ++ncol;
+ if(isComment)
+ break;
+ }
+ if(verbose)
+ cout << endl;
+ if(dataVector.size())
+ assert(ncol>=col);
+ }
+ ++nrow;
+ }
+ assert(dataVector.size());
+ }
+ else{//space or tab delimited fields
+ if(verbose)
+ std::cout << "space or tab delimited fields" << std::endl;
+ string spaceRecord;
+ while(!getline(m_ifstream, spaceRecord).eof()){
+ withinRange=true;
+ if(nrow<m_minRow)
+ withinRange=false;
+ if(m_maxRow>m_minRow)
+ if(nrow>m_maxRow)
+ withinRange=false;
+ if(withinRange){
+ if(verbose>1)
+ cout << spaceRecord << endl;
+ istringstream lineStream(spaceRecord);
+ string item;
+ int ncol=0;
+ bool isComment=false;
+ while(lineStream >> item){
+ if(verbose)
+ cout << item << " ";
+ // istringstream itemStream(item);
+ unsigned pos=item.find(m_comment);
+ if(pos!=std::string::npos){
+ if(pos>0)
+ item=item.substr(0,pos-1);
+ else
+ break;
+ if(verbose)
+ std::cout << "comment found, string is " << item << std::endl;
+ isComment=true;
+ }
+ T value=string2type<T>(item);
+ if(ncol==col){
+ if((value>=m_min&&value<=m_max)||m_max<=m_min)
+ dataVector.push_back(value);
+ }
+ ++ncol;
+ if(isComment)
+ break;
+ }
+ if(verbose>1)
+ cout << endl;
+ if(verbose)
+ cout << "number of columns: " << ncol << endl;
+ if(dataVector.size())
+ assert(ncol>=col);
+ }
+ ++nrow;
+ }
+ }
+ return dataVector.size();
+}
+
template<class T> unsigned int FileReaderAscii::readData(vector<vector<T> > &dataVector, const vector<int> &cols){
bool verbose=false;
dataVector.clear();
--
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