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 fa5dc5671d8d79b9750b81f5325384e5d72ed9f9
Author: Pieter Kempeneers <kempe...@gmail.com>
Date: Thu Mar 14 17:26:13 2013 +0100
pkfilter.cc supports fwhm and srf filtering
---
src/algorithms/Filter.cc | 86 +++++++++++-------
src/algorithms/Filter.h | 212 ++++++++++++++++++++++++++++---------------
src/algorithms/Makefile.am | 6 +-
src/algorithms/Makefile.in | 67 ++++++++------
src/algorithms/StatFactory.h | 60 ++++++------
src/apps/pkfilter.cc | 110 +++++++++++++++++++---
6 files changed, 364 insertions(+), 177 deletions(-)
diff --git a/src/algorithms/Filter.cc b/src/algorithms/Filter.cc
index b99a8df..21009de 100644
--- a/src/algorithms/Filter.cc
+++ b/src/algorithms/Filter.cc
@@ -104,38 +104,56 @@ void filter::Filter::doit(const ImgReaderGdal& input,
ImgWriterGdal& output, sho
}
}
-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;
- void* pProgressArg=NULL;
- GDALProgressFunc pfnProgress=GDALTermProgress;
- double progress=0;
- pfnProgress(progress,pszMessage,pProgressArg);
- for(int y=0;y<input.nrOfRow();++y){
- for(int iband=0;iband<input.nrOfBand();++iband)
- input.readData(lineInput[iband],GDT_Float64,y,iband);
- vector<double> pixelInput(input.nrOfBand());
- vector<double> pixelOutput;
- for(int x=0;x<input.nrOfCol();++x){
- 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];
- }
- for(int iband=0;iband<output.nrOfBand();++iband){
- try{
- output.writeData(lineOutput[iband],GDT_Float64,y,iband);
- }
- catch(string errorstring){
- cerr << errorstring << "in band " << iband << ", line " << y << endl;
- }
- }
- progress=(1.0+y)/output.nrOfRow();
- pfnProgress(progress,pszMessage,pProgressArg);
- }
-}
-
+// 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;
+// void* pProgressArg=NULL;
+// GDALProgressFunc pfnProgress=GDALTermProgress;
+// double progress=0;
+// pfnProgress(progress,pszMessage,pProgressArg);
+// for(int y=0;y<input.nrOfRow();++y){
+// for(int iband=0;iband<input.nrOfBand();++iband)
+// input.readData(lineInput[iband],GDT_Float64,y,iband);
+// applyFwhm<double>(wavelengthIn,lineInput,wavelengthOut,fwhm,
interpolationType, lineOutput, verbose);
+// for(int iband=0;iband<output.nrOfBand();++iband){
+// try{
+// output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+// }
+// catch(string errorstring){
+// cerr << errorstring << "in band " << iband << ", line " << y <<
endl;
+// }
+// }
+// progress=(1.0+y)/output.nrOfRow();
+// pfnProgress(progress,pszMessage,pProgressArg);
+// }
+// }
+// void filter::Filter::applySrf(const vector<double> &wavelengthIn, const
ImgReaderGdal& input, const vector< Vector2d<double> > &srf, const std::string&
interpolationType, ImgWriterGdal& output, bool verbose){
+// assert(output.nrOfBand()==srf.size());
+// double centreWavelength=0;
+// Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
+// const char* pszMessage;
+// void* pProgressArg=NULL;
+// GDALProgressFunc pfnProgress=GDALTermProgress;
+// double progress=0;
+// pfnProgress(progress,pszMessage,pProgressArg);
+// for(int y=0;y<input.nrOfRow();++y){
+// for(int iband=0;iband<input.nrOfBand();++iband)
+// input.readData(lineInput[iband],GDT_Float64,y,iband);
+// for(int isrf=0;isrf<srf.size();++isrf){
+// vector<double> lineOutput(input.nrOfCol());
+// centreWavelength=applySrf<double>(wavelengthIn,lineInput,srf[isrf],
interpolationType, lineOutput, verbose);
+// for(int iband=0;iband<output.nrOfBand();++iband){
+// try{
+// output.writeData(lineOutput,GDT_Float64,y,isrf);
+// }
+// catch(string errorstring){
+// cerr << errorstring << "in band " << iband << ", line " << y <<
endl;
+// }
+// }
+// }
+// progress=(1.0+y)/output.nrOfRow();
+// pfnProgress(progress,pszMessage,pProgressArg);
+// }
+// }
diff --git a/src/algorithms/Filter.h b/src/algorithms/Filter.h
index d20203e..75bce71 100644
--- a/src/algorithms/Filter.h
+++ b/src/algorithms/Filter.h
@@ -55,9 +55,11 @@ public:
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> double applySrf(const vector<double> &wavelengthIn, const
Vector2d<T>& input, const Vector2d<double>& srf, const std::string&
interpolationType, vector<T>& output, double delta=1.0, bool normalize=false,
bool verbose=false);
+ // void applySrf(const vector<double> &wavelengthIn, const ImgReaderGdal&
input, const vector< Vector2d<double> > &srf, const std::string&
interpolationType, ImgWriterGdal& output, 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);
+ 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, 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);
@@ -98,69 +100,88 @@ private:
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,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);
-// }
-// }
+//input[band][sample], output[sample]
+//returns wavelength for which srf is maximum
+ template<class T> double Filter::applySrf(const vector<double>
&wavelengthIn, const Vector2d<T>& input, const Vector2d<double>& srf, const
std::string& interpolationType, vector<T>& output, double delta, bool
normalize, bool verbose)
+{
+ assert(srf.size()==2);//[0]: wavelength, [1]: response function
+ int nband=srf[0].size();
+ unsigned int nsample=input[0].size();
+ output.resize(nsample);
+ double start=floor(wavelengthIn[0]);
+ double end=ceil(wavelengthIn.back());
+ if(verbose)
+ std::cout << "wavelengths in [" << start << "," << end << "]" << std::endl
<< std::flush;
+
+ statfactory::StatFactory stat;
+
+ gsl_interp_accel *acc;
+ stat.allocAcc(acc);
+ gsl_spline *spline;
+ stat.getSpline(interpolationType,nband,spline);
+ stat.initSpline(spline,&(srf[0][0]),&(srf[1][0]),nband);
+ // 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);
+ if(verbose)
+ std::cout << "calculating norm of srf" << std::endl << std::flush;
+ double norm=0;
+ norm=gsl_spline_eval_integ(spline,srf[0].front(),srf[0].back(),acc);
+ if(verbose)
+ std::cout << "norm of srf: " << norm << std::endl << std::flush;
+ gsl_spline_free(spline);
+ gsl_interp_accel_free(acc);
+ //interpolate input and srf to delta
+
+ vector<double> wavelength_fine;
+ for(double
win=floor(wavelengthIn[0]);win<=ceil(wavelengthIn.back());win+=delta)
+ wavelength_fine.push_back(win);
+
+ if(verbose)
+ std::cout << "interpolate wavelengths to " << wavelength_fine.size() << "
entries " << std::endl;
+ vector<double> srf_fine;//spectral response function, interpolated for
wavelength_fine
+
+
stat.interpolateUp(srf[0],srf[1],wavelength_fine,interpolationType,srf_fine,verbose);
+ assert(srf_fine.size()==wavelength_fine.size());
+
+ gsl_interp_accel *accOut;
+ stat.allocAcc(accOut);
+ gsl_spline *splineOut;
+ stat.getSpline(interpolationType,wavelength_fine.size(),splineOut);
+ assert(splineOut);
+
+ for(int isample=0;isample<nsample;++isample){
+ vector<T> inputValues;
+ input.selectCol(isample,inputValues);
+ assert(wavelengthIn.size()==inputValues.size());
+ vector<double> input_fine;
+ vector<double> product(wavelength_fine.size());
+
stat.interpolateUp(wavelengthIn,inputValues,wavelength_fine,interpolationType,input_fine,verbose);
+
+ for(int iband=0;iband<input_fine.size();++iband)
+ product[iband]=input_fine[iband]*srf_fine[iband];
+
+ assert(input_fine.size()==srf_fine.size());
+ assert(input_fine.size()==wavelength_fine.size());
+
stat.initSpline(splineOut,&(wavelength_fine[0]),&(product[0]),wavelength_fine.size());
+ if(normalize)
+ output[isample]=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
+ else
+ output[isample]=gsl_spline_eval_integ(splineOut,start,end,accOut);
+ }
+ gsl_spline_free(splineOut);
+ gsl_interp_accel_free(accOut);
+
+ 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;
+ }
+ }
+ return(srf[0][maxIndex]);
+}
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
@@ -200,18 +221,67 @@ template<class T> void Filter::applyFwhm(const
vector<double> &wavelengthIn, con
tf(indexIn,indexOut)/=stddev[indexOut];
norm+=tf(indexIn,indexOut);
}
- //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
- //todo: support more than one sample
output[indexOut]=0;
for(int indexIn=0;indexIn<nbandIn;++indexIn)
output[indexOut]+=input_fine[indexIn]*tf(indexIn,indexOut)/norm;
}
}
+
+ //input[inBand][sample], output[outBand][sample]
+template<class T> void Filter::applyFwhm(const vector<double> &wavelengthIn,
const Vector2d<T>& input, const vector<double> &wavelengthOut, const
vector<double> &fwhm, const std::string& interpolationType, Vector2d<T>&
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
+ statfactory::StatFactory stat;
+ vector<double> wavelength_fine;
+ for(double
win=floor(wavelengthIn[0]);win<=ceil(wavelengthIn.back());win+=delta)
+ wavelength_fine.push_back(win);
+ assert(wavelengthOut.size()==fwhm.size());
+ assert(wavelengthIn[0]<wavelengthOut[0]);
+ assert(wavelengthIn.back()>wavelengthOut.back());
+ 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;
+ }
+ int nbandIn=wavelength_fine.size();
+ int nbandOut=wavelengthOut.size();
+ output.resize(nbandOut,input[0].size());
+
+ gsl::matrix tf(nbandIn,nbandOut);
+ vector<double> norm(nbandOut);
+ for(int indexOut=0;indexOut<nbandOut;++indexOut){
+ norm[indexOut]=0;
+ for(int indexIn=0;indexIn<nbandIn;++indexIn){
+ tf(indexIn,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[indexOut]+=tf(indexIn,indexOut);
+ }
+ }
+
+ for(int isample=0;isample<input[0].size();++isample){
+ vector<T> inputValues;
+ input.selectCol(isample,inputValues);
+ assert(wavelengthIn.size()==inputValues.size());
+ for(int indexOut=0;indexOut<nbandOut;++indexOut){
+ vector<double> input_fine;
+
stat.interpolateUp(wavelengthIn,inputValues,wavelength_fine,interpolationType,input_fine,verbose);
+ output[indexOut][isample]=0;
+ for(int indexIn=0;indexIn<nbandIn;++indexIn){
+
output[indexOut][isample]+=input_fine[indexIn]*tf(indexIn,indexOut)/norm[indexOut];
+ }
+ }
+ }
+}
+
template<class T> void Filter::doit(const vector<T>& input, vector<T>& output,
int down, int offset)
{
output.resize((input.size()-offset+down-1)/down);
diff --git a/src/algorithms/Makefile.am b/src/algorithms/Makefile.am
index f44c11a..40d4bee 100644
--- a/src/algorithms/Makefile.am
+++ b/src/algorithms/Makefile.am
@@ -7,7 +7,7 @@ AM_CXXFLAGS = -I$(top_srcdir)/src $(GDAL_CFLAGS) @AM_CXXFLAGS@
# the program to build (the names of the final binaries)
#do not want to install pktestoption
-#noinst_PROGRAMS = pktestStat
+noinst_PROGRAMS = pktestStat
###############################################################################
# THE LIBRARIES TO BUILD
@@ -35,5 +35,5 @@ libalgorithms_a_SOURCES = $(libalgorithms_a_HEADERS) Egcs.cc
Filter2d.cc Filter.
###############################################################################
# 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
+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 f8258ad..a5e3600 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@am__append_1 = myfann_cpp.h
@USE_NLOPT_TRUE@am__append_2 = 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.@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 OptFactory.h
@@ -218,14 +228,6 @@ top_builddir = @top_builddir@
top_srcdir = @top_srcdir@
###############################################################################
-# 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
###############################################################################
@@ -242,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:
@@ -284,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)
@@ -294,6 +307,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/Egcs.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/Filter.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/Filter2d.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pktestStat.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/svm.Po@am__quote@
.cc.o:
@@ -428,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"; \
@@ -460,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)
@@ -530,23 +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
-###############################################################################
-
-# 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
# 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 350a49d..504c6c7 100644
--- a/src/algorithms/StatFactory.h
+++ b/src/algorithms/StatFactory.h
@@ -551,36 +551,36 @@ template<class T> void StatFactory::interpolateUp(const
vector<double>& waveleng
gsl_interp_accel_free(acc);
}
-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);
- 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);
-}
+// 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);
+// 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[isample].push_back(dout);
+// }
+// }
+// }
+// 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/pkfilter.cc b/src/apps/pkfilter.cc
index 0ee8774..cbf9ca0 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<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<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");
@@ -50,9 +50,10 @@ int main(int argc,char **argv) {
Optionpk<std::string> tap_opt("tap", "tap", "text file containing taps used
for spatial filtering (from ul to lr). Use dimX and dimY to specify tap
dimensions in x and y. Leave empty for not using taps");
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 ...)");
+ Optionpk<std::string> srf_opt("srf", "srf", "list of ASCII files containing
spectral response functions (two columns: wavelength response)");
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> interpolationType_opt("interp", "interp", "type of
interpolation for spectral filtering (see
http://www.gnu.org/software/gsl/manual/html_node/Interpolation-Types.html)","akima");
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)");
@@ -77,6 +78,7 @@ int main(int argc,char **argv) {
tap_opt.retrieveOption(argc,argv);
tapz_opt.retrieveOption(argc,argv);
fwhm_opt.retrieveOption(argc,argv);
+ srf_opt.retrieveOption(argc,argv);
wavelengthIn_opt.retrieveOption(argc,argv);
wavelengthOut_opt.retrieveOption(argc,argv);
down_opt.retrieveOption(argc,argv);
@@ -128,9 +130,11 @@ int main(int argc,char **argv) {
}
try{
assert(output_opt.size());
- if(wavelengthOut_opt.size())
+ if(fwhm_opt.size()||srf_opt.size()){
//todo: support down and offset
-
output.open(output_opt[0],input.nrOfCol(),input.nrOfRow(),wavelengthOut_opt.size(),theType,imageType,option_opt);
+ int nband=fwhm_opt.size()? fwhm_opt.size():srf_opt.size();
+
output.open(output_opt[0],input.nrOfCol(),input.nrOfRow(),nband,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);
}
@@ -153,10 +157,10 @@ int main(int argc,char **argv) {
else if(input.getColorTable()!=NULL)
output.setColorTable(input.getColorTable());
- // if(input.isGeoRef()){
- // output.setProjection(input.getProjection());
- // output.copyGeoTransform(input);
- // }
+ if(input.isGeoRef()){
+ output.setProjection(input.getProjection());
+ output.copyGeoTransform(input);
+ }
filter2d::Filter2d filter2d;
filter::Filter filter1d;
@@ -210,12 +214,96 @@ int main(int argc,char **argv) {
filter1d.setTaps(tapz_opt);
filter1d.doit(input,output,down_opt[0]);
}
- else if(wavelengthOut_opt.size()){
+ else if(fwhm_opt.size()){
if(verbose_opt[0])
- std::cout << "spectral filtering to " << wavelengthOut_opt.size() << "
bands with provided fwhm " << std::endl;
+ std::cout << "spectral filtering to " << fwhm_opt.size() << " bands with
provided fwhm " << std::endl;
assert(wavelengthOut_opt.size()==fwhm_opt.size());
assert(wavelengthIn_opt.size());
-
filter1d.applyFwhm(wavelengthIn_opt,input,wavelengthOut_opt,fwhm_opt,interpolationType_opt[0],output,verbose_opt[0]);
+
+ Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
+ Vector2d<double> lineOutput(wavelengthOut_opt.size(),input.nrOfCol());
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+ for(int y=0;y<input.nrOfRow();++y){
+ for(int iband=0;iband<input.nrOfBand();++iband)
+ input.readData(lineInput[iband],GDT_Float64,y,iband);
+
filter1d.applyFwhm<double>(wavelengthIn_opt,lineInput,wavelengthOut_opt,fwhm_opt,
interpolationType_opt[0], lineOutput, verbose_opt[0]);
+ for(int iband=0;iband<output.nrOfBand();++iband){
+ try{
+ output.writeData(lineOutput[iband],GDT_Float64,y,iband);
+ }
+ catch(string errorstring){
+ cerr << errorstring << "in band " << iband << ", line " << y << endl;
+ }
+ }
+ progress=(1.0+y)/output.nrOfRow();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+ //
filter1d.applyFwhm(wavelengthIn_opt,input,wavelengthOut_opt,fwhm_opt,interpolationType_opt[0],output,verbose_opt[0]);
+ }
+ else if(srf_opt.size()){
+ if(verbose_opt[0])
+ std::cout << "spectral filtering to " << srf_opt.size() << " bands with
provided SRF " << std::endl;
+ assert(wavelengthIn_opt.size());
+ vector< Vector2d<double> > srf(srf_opt.size());//[0] srf_nr, [1]:
wavelength, [2]: response
+ ifstream srfFile;
+ for(int isrf=0;isrf<srf_opt.size();++isrf){
+ srf[isrf].resize(2);
+ srfFile.open(srf_opt[isrf].c_str());
+ double v;
+ //add 0 to make sure srf is 0 at boundaries after interpolation step
+ srf[isrf][0].push_back(0);
+ srf[isrf][1].push_back(0);
+ srf[isrf][0].push_back(1);
+ srf[isrf][1].push_back(0);
+ while(srfFile >> v){
+ srf[isrf][0].push_back(v);
+ srfFile >> v;
+ srf[isrf][1].push_back(v);
+ }
+ srfFile.close();
+ //add 0 to make sure srf[isrf] is 0 at boundaries after interpolation
step
+ srf[isrf][0].push_back(srf[isrf][0].back()+1);
+ srf[isrf][1].push_back(0);
+ srf[isrf][0].push_back(srf[isrf][0].back()+1);
+ srf[isrf][1].push_back(0);
+ if(verbose_opt[0])
+ cout << "srf file details: " << srf[isrf][0].size() << " wavelengths
defined" << endl;
+ }
+ assert(output.nrOfBand()==srf.size());
+ double centreWavelength=0;
+ Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ double progress=0;
+ pfnProgress(progress,pszMessage,pProgressArg);
+ for(int y=0;y<input.nrOfRow();++y){
+ for(int iband=0;iband<input.nrOfBand();++iband)
+ input.readData(lineInput[iband],GDT_Float64,y,iband);
+ for(int isrf=0;isrf<srf.size();++isrf){
+ vector<double> lineOutput(input.nrOfCol());
+ double delta=1.0;
+ bool normalize=true;
+
centreWavelength=filter1d.applySrf<double>(wavelengthIn_opt,lineInput,srf[isrf],
interpolationType_opt[0], lineOutput, delta, normalize, verbose_opt[0]);
+ if(verbose_opt[0])
+ std::cout << "centre wavelength srf " << isrf << ": " <<
centreWavelength << std::endl;
+ try{
+ output.writeData(lineOutput,GDT_Float64,y,isrf);
+ }
+ catch(string errorstring){
+ cerr << errorstring << "in srf " << srf_opt[isrf] << ", line " << y
<< endl;
+ }
+
+ }
+ progress=(1.0+y)/output.nrOfRow();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+
+ //
filter1d.applySrf(wavelengthIn_opt,input,srf,interpolationType_opt[0],output,verbose_opt[0]);
}
else{
if(colorTable_opt.size())
--
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
