[pktools] 137/375: central wavelength SRF in Filter.h

Bas Couwenberg sebastic at xs4all.nl
Wed Dec 3 21:54:07 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 8c5dd47b2cdd77ea0063f93e6486450f62d0ca0c
Author: Pieter Kempeneers <kempenep at gmail.com>
Date:   Tue Nov 5 16:30:26 2013 +0100

    central wavelength SRF in Filter.h
---
 src/algorithms/Filter.h   | 59 +++++++++++++++++++++++++++++++----------------
 src/apps/pkfilterascii.cc | 38 +++++++++++++++++++-----------
 2 files changed, 63 insertions(+), 34 deletions(-)

diff --git a/src/algorithms/Filter.h b/src/algorithms/Filter.h
index 65d0028..26e530f 100644
--- a/src/algorithms/Filter.h
+++ b/src/algorithms/Filter.h
@@ -158,12 +158,15 @@ private:
   assert(wavelengthIn.size()==input.size());
   std::vector<double> input_fine;
   std::vector<double> product(wavelength_fine.size());
+  std::vector<double> wavelengthOut(wavelength_fine.size());
   stat.interpolateUp(wavelengthIn,input,wavelength_fine,interpolationType,input_fine,verbose);
 
   if(verbose)
     std::cout << "input_fine.size(): " << input_fine.size() << std::endl;
-  for(int iband=0;iband<input_fine.size();++iband)
+  for(int iband=0;iband<input_fine.size();++iband){
     product[iband]=input_fine[iband]*srf_fine[iband];
+    wavelengthOut[iband]=wavelength_fine[iband]*srf_fine[iband];
+  }
 
   assert(input_fine.size()==srf_fine.size());
   assert(input_fine.size()==wavelength_fine.size());
@@ -172,18 +175,23 @@ private:
     output=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
   else
     output=gsl_spline_eval_integ(splineOut,start,end,accOut);
+
+  stat.initSpline(splineOut,&(wavelength_fine[0]),&(wavelengthOut[0]),wavelength_fine.size());
+  double centreWavelength=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
+  
   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]);
+  // 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]);
+  return(centreWavelength);
 }
 
 //input[band][sample], output[sample] (if !transposeInput)
@@ -233,6 +241,8 @@ private:
   stat.getSpline(interpolationType,wavelength_fine.size(),splineOut);
   assert(splineOut);
 
+  std::vector<double> wavelengthOut;
+  double centreWavelength=0;
   for(int isample=0;isample<nsample;++isample){
     if((isample+1+down/2)%down)
       continue;
@@ -246,8 +256,11 @@ private:
     std::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)
+    for(int iband=0;iband<input_fine.size();++iband){
       product[iband]=input_fine[iband]*srf_fine[iband];
+      if(wavelengthOut.size()<input_fine.size())
+        wavelengthOut.push_back(wavelength_fine[iband]*srf_fine[iband]);
+    }
 
     assert(input_fine.size()==srf_fine.size());
     assert(input_fine.size()==wavelength_fine.size());
@@ -256,19 +269,25 @@ private:
       output[isample/down]=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
     else
       output[isample/down]=gsl_spline_eval_integ(splineOut,start,end,accOut);
+
+    stat.initSpline(splineOut,&(wavelength_fine[0]),&(wavelengthOut[0]),wavelength_fine.size());
+    if(centreWavelength>0);
+    else
+      centreWavelength=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
   }
   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]);
+  // 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]);
+  return(centreWavelength);
 }
 
 template<class T> void Filter::applyFwhm(const std::vector<double> &wavelengthIn, const std::vector<T>& input, const std::vector<double> &wavelengthOut, const std::vector<double> &fwhm, const std::string& interpolationType, std::vector<T>& output, bool verbose){
diff --git a/src/apps/pkfilterascii.cc b/src/apps/pkfilterascii.cc
index 3c393cf..87ca91b 100644
--- a/src/apps/pkfilterascii.cc
+++ b/src/apps/pkfilterascii.cc
@@ -229,6 +229,8 @@ int main(int argc,char **argv) {
   if(!output_opt.empty())
     outputStream.open(output_opt[0].c_str(),ios::out);
 
+  if(verbose_opt[0])
+    std::cout << "stream to output" << std::endl;
   if(transpose_opt[0]){
     for(int icol=0;icol<inputCols_opt.size();++icol){
       for(int iband=0;iband<filteredData[icol].size();++iband){
@@ -252,21 +254,29 @@ int main(int argc,char **argv) {
   else{
     // int nband=wavelengthOut.size()? wavelengthOut.size() : filteredData[0].size();
     int nband=0;
-    switch(filter::Filter::getFilterType(method_opt[0])){
-    case(filter::dwtForward):
-      nband=filteredData[0].size();
-      break;
-    case(filter::dwtInverse):
-      nband=filteredData[0].size();
-      break;
-    case(filter::dwtQuantize):
-      nband=wavelengthOut.size();
-      assert(wavelengthOut.size()==nband);
-      assert(filteredData[0].size()==nband);
-      break;
-    default:
+    if(method_opt.size()){
+      switch(filter::Filter::getFilterType(method_opt[0])){
+      case(filter::dwtForward):
+        nband=filteredData[0].size();
+        break;
+      case(filter::dwtInverse):
+        nband=filteredData[0].size();
+        break;
+      case(filter::dwtQuantize):
+        nband=wavelengthOut.size();
+        assert(wavelengthOut.size()==nband);
+        assert(filteredData[0].size()==nband);
+        break;
+      default:
+        nband=wavelengthOut.size();
+        break;
+      }
+    }
+    else
       nband=wavelengthOut.size();
-      break;
+    if(verbose_opt[0]){
+      std::cout << "number of bands: " << nband << std::endl;
+      std::cout << "wavelengthOut.size(): " << wavelengthOut.size() << std::endl;
     }
     for(int iband=0;iband<nband;++iband){
       if(!output_opt.empty()){

-- 
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