[pktools] 72/375: pkfilter.cc supports fwhm and srf filtering

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 fa5dc5671d8d79b9750b81f5325384e5d72ed9f9
Author: Pieter Kempeneers <kempenep at 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 at am__append_1 = myfann_cpp.h
 @USE_NLOPT_TRUE at 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. 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 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 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:
@@ -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



More information about the Pkg-grass-devel mailing list