pktools  2.6.7
Processing Kernel for geospatial data
Filter.cc
1 /**********************************************************************
2 Filter.cc: class for filtering
3 Copyright (C) 2008-2012 Pieter Kempeneers
4 
5 This file is part of pktools
6 
7 pktools is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 pktools is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with pktools. If not, see <http://www.gnu.org/licenses/>.
19 ***********************************************************************/
20 #include "Filter.h"
21 #include <assert.h>
22 #include <math.h>
23 #include <iostream>
24 
25 using namespace std;
26 
27 filter::Filter::Filter(void)
28  : m_padding("symmetric")
29 {
30 }
31 
32 
33 filter::Filter::Filter(const vector<double> &taps)
34  : m_padding("symmetric")
35 {
36  setTaps(taps);
37 }
38 
39 void filter::Filter::setTaps(const vector<double> &taps, bool normalize)
40 {
41  m_taps.resize(taps.size());
42  double norm=0;
43  for(int itap=0;itap<taps.size();++itap)
44  norm+=taps[itap];
45  if(norm){
46  for(int itap=0;itap<taps.size();++itap)
47  m_taps[itap]=taps[itap]/norm;
48  }
49  else
50  m_taps=taps;
51  assert(m_taps.size()%2);
52 }
53 
54 unsigned int filter::Filter::pushNoDataValue(double noDataValue){
55  if(find(m_noDataValues.begin(),m_noDataValues.end(),noDataValue)==m_noDataValues.end())
56  m_noDataValues.push_back(noDataValue);
57  return m_noDataValues.size();
58 };
59 
60 unsigned int filter::Filter::setNoDataValues(std::vector<double> vnodata){
61  m_noDataValues=vnodata;
62  return m_noDataValues.size();
63 };
64 
65 void filter::Filter::dwtForward(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family){
66  const char* pszMessage;
67  void* pProgressArg=NULL;
68  GDALProgressFunc pfnProgress=GDALTermProgress;
69  double progress=0;
70  pfnProgress(progress,pszMessage,pProgressArg);
71  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
72  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
73  for(int y=0;y<input.nrOfRow();++y){
74  for(int iband=0;iband<input.nrOfBand();++iband)
75  input.readData(lineInput[iband],y,iband);
76  vector<double> pixelInput(input.nrOfBand());
77  for(int x=0;x<input.nrOfCol();++x){
78  pixelInput=lineInput.selectCol(x);
79  dwtForward(pixelInput,wavelet_type,family);
80  for(int iband=0;iband<input.nrOfBand();++iband)
81  lineOutput[iband][x]=pixelInput[iband];
82  }
83  for(int iband=0;iband<input.nrOfBand();++iband){
84  try{
85  output.writeData(lineOutput[iband],y,iband);
86  }
87  catch(string errorstring){
88  cerr << errorstring << "in band " << iband << ", line " << y << endl;
89  }
90  }
91  progress=(1.0+y)/output.nrOfRow();
92  pfnProgress(progress,pszMessage,pProgressArg);
93  }
94 }
95 
96 void filter::Filter::dwtInverse(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family){
97  const char* pszMessage;
98  void* pProgressArg=NULL;
99  GDALProgressFunc pfnProgress=GDALTermProgress;
100  double progress=0;
101  pfnProgress(progress,pszMessage,pProgressArg);
102  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
103  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
104  for(int y=0;y<input.nrOfRow();++y){
105  for(int iband=0;iband<input.nrOfBand();++iband)
106  input.readData(lineInput[iband],y,iband);
107  vector<double> pixelInput(input.nrOfBand());
108  for(int x=0;x<input.nrOfCol();++x){
109  pixelInput=lineInput.selectCol(x);
110  dwtInverse(pixelInput,wavelet_type,family);
111  for(int iband=0;iband<input.nrOfBand();++iband)
112  lineOutput[iband][x]=pixelInput[iband];
113  }
114  for(int iband=0;iband<input.nrOfBand();++iband){
115  try{
116  output.writeData(lineOutput[iband],y,iband);
117  }
118  catch(string errorstring){
119  cerr << errorstring << "in band " << iband << ", line " << y << endl;
120  }
121  }
122  progress=(1.0+y)/output.nrOfRow();
123  pfnProgress(progress,pszMessage,pProgressArg);
124  }
125 }
126 
127 void filter::Filter::dwtCut(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, double cut){
128  const char* pszMessage;
129  void* pProgressArg=NULL;
130  GDALProgressFunc pfnProgress=GDALTermProgress;
131  double progress=0;
132  pfnProgress(progress,pszMessage,pProgressArg);
133  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
134  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
135  for(int y=0;y<input.nrOfRow();++y){
136  for(int iband=0;iband<input.nrOfBand();++iband)
137  input.readData(lineInput[iband],y,iband);
138  vector<double> pixelInput(input.nrOfBand());
139  for(int x=0;x<input.nrOfCol();++x){
140  pixelInput=lineInput.selectCol(x);
141  dwtCut(pixelInput,wavelet_type,family,cut);
142  for(int iband=0;iband<input.nrOfBand();++iband)
143  lineOutput[iband][x]=pixelInput[iband];
144  }
145  for(int iband=0;iband<input.nrOfBand();++iband){
146  try{
147  output.writeData(lineOutput[iband],y,iband);
148  }
149  catch(string errorstring){
150  cerr << errorstring << "in band " << iband << ", line " << y << endl;
151  }
152  }
153  progress=(1.0+y)/output.nrOfRow();
154  pfnProgress(progress,pszMessage,pProgressArg);
155  }
156 }
157 
158 void filter::Filter::dwtCutFrom(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, int band){
159  const char* pszMessage;
160  void* pProgressArg=NULL;
161  GDALProgressFunc pfnProgress=GDALTermProgress;
162  double progress=0;
163  pfnProgress(progress,pszMessage,pProgressArg);
164  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
165  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
166  for(int y=0;y<input.nrOfRow();++y){
167  for(int iband=0;iband<input.nrOfBand();++iband)
168  input.readData(lineInput[iband],y,iband);
169  vector<double> pixelInput(input.nrOfBand());
170  for(int x=0;x<input.nrOfCol();++x){
171  pixelInput=lineInput.selectCol(x);
172  dwtForward(pixelInput,wavelet_type,family);
173  for(int iband=0;iband<input.nrOfBand();++iband){
174  if(iband>=band)
175  pixelInput[iband]=0;
176  }
177  dwtInverse(pixelInput,wavelet_type,family);
178  for(int iband=0;iband<input.nrOfBand();++iband)
179  lineOutput[iband][x]=pixelInput[iband];
180  }
181  for(int iband=0;iband<input.nrOfBand();++iband){
182  try{
183  output.writeData(lineOutput[iband],y,iband);
184  }
185  catch(string errorstring){
186  cerr << errorstring << "in band " << iband << ", line " << y << endl;
187  }
188  }
189  progress=(1.0+y)/output.nrOfRow();
190  pfnProgress(progress,pszMessage,pProgressArg);
191  }
192 }
193 
194 //todo: support different padding strategies
195 void filter::Filter::dwtForward(std::vector<double>& data, const std::string& wavelet_type, int family){
196  int origsize=data.size();
197  //make sure data size if power of 2
198  while(data.size()&(data.size()-1))
199  data.push_back(data.back());
200 
201  int nsize=data.size();
202  gsl_wavelet *w;
203  gsl_wavelet_workspace *work;
204  assert(nsize);
205  w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
206  work=gsl_wavelet_workspace_alloc(nsize);
207  gsl_wavelet_transform_forward(w,&(data[0]),1,nsize,work);
208  data.erase(data.begin()+origsize,data.end());
209  gsl_wavelet_free (w);
210  gsl_wavelet_workspace_free (work);
211 }
212 
213 //todo: support different padding strategies
214 void filter::Filter::dwtInverse(std::vector<double>& data, const std::string& wavelet_type, int family){
215  int origsize=data.size();
216  //make sure data size if power of 2
217  while(data.size()&(data.size()-1))
218  data.push_back(data.back());
219  int nsize=data.size();
220  gsl_wavelet *w;
221  gsl_wavelet_workspace *work;
222  assert(nsize);
223  w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
224  work=gsl_wavelet_workspace_alloc(nsize);
225  gsl_wavelet_transform_inverse(w,&(data[0]),1,nsize,work);
226  data.erase(data.begin()+origsize,data.end());
227  gsl_wavelet_free (w);
228  gsl_wavelet_workspace_free (work);
229 }
230 
231 //todo: support different padding strategies
232 void filter::Filter::dwtCut(std::vector<double>& data, const std::string& wavelet_type, int family, double cut){
233  int origsize=data.size();
234  //make sure data size if power of 2
235  while(data.size()&(data.size()-1))
236  data.push_back(data.back());
237  int nsize=data.size();
238  gsl_wavelet *w;
239  gsl_wavelet_workspace *work;
240  assert(nsize);
241  w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
242  work=gsl_wavelet_workspace_alloc(nsize);
243  gsl_wavelet_transform_forward(w,&(data[0]),1,nsize,work);
244  std::vector<double> abscoeff(data.size());
245  size_t* p=new size_t[data.size()];
246  for(int index=0;index<data.size();++index){
247  abscoeff[index]=fabs(data[index]);
248  }
249  int nc=(100-cut)/100.0*nsize;
250  gsl_sort_index(p,&(abscoeff[0]),1,nsize);
251  for(int i=0;(i+nc)<nsize;i++)
252  data[p[i]]=0;
253  gsl_wavelet_transform_inverse(w,&(data[0]),1,nsize,work);
254  data.erase(data.begin()+origsize,data.end());
255  delete[] p;
256  gsl_wavelet_free (w);
257  gsl_wavelet_workspace_free (work);
258 }
259 
260 void filter::Filter::morphology(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short verbose)
261 {
262  // bool bverbose=(verbose>1)? true:false;
263  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
264  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
265  const char* pszMessage;
266  void* pProgressArg=NULL;
267  GDALProgressFunc pfnProgress=GDALTermProgress;
268  double progress=0;
269  pfnProgress(progress,pszMessage,pProgressArg);
270  for(int y=0;y<input.nrOfRow();++y){
271  for(int iband=0;iband<input.nrOfBand();++iband)
272  input.readData(lineInput[iband],y,iband);
273  vector<double> pixelInput(input.nrOfBand());
274  vector<double> pixelOutput(input.nrOfBand());
275  for(int x=0;x<input.nrOfCol();++x){
276  pixelInput=lineInput.selectCol(x);
277  filter(pixelInput,pixelOutput,method,dim);
278  // morphology(pixelInput,pixelOutput,method,dim,bverbose);
279  for(int iband=0;iband<input.nrOfBand();++iband)
280  lineOutput[iband][x]=pixelOutput[iband];
281  }
282  for(int iband=0;iband<input.nrOfBand();++iband){
283  try{
284  output.writeData(lineOutput[iband],y,iband);
285  }
286  catch(string errorstring){
287  cerr << errorstring << "in band " << iband << ", line " << y << endl;
288  }
289  }
290  progress=(1.0+y)/output.nrOfRow();
291  pfnProgress(progress,pszMessage,pProgressArg);
292  }
293 }
294 
295 void filter::Filter::smoothNoData(ImgReaderGdal& input, const std::string& interpolationType, ImgWriterGdal& output)
296 {
297  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
298  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
299  const char* pszMessage;
300  void* pProgressArg=NULL;
301  GDALProgressFunc pfnProgress=GDALTermProgress;
302  double progress=0;
303  pfnProgress(progress,pszMessage,pProgressArg);
304  for(int y=0;y<input.nrOfRow();++y){
305  for(int iband=0;iband<input.nrOfBand();++iband)
306  input.readData(lineInput[iband],y,iband);
307  vector<double> pixelInput(input.nrOfBand());
308  vector<double> pixelOutput(input.nrOfBand());
309  for(int x=0;x<input.nrOfCol();++x){
310  pixelInput=lineInput.selectCol(x);
311  smoothNoData(pixelInput,interpolationType,pixelOutput);
312  for(int iband=0;iband<input.nrOfBand();++iband)
313  lineOutput[iband][x]=pixelOutput[iband];
314  }
315  for(int iband=0;iband<input.nrOfBand();++iband){
316  try{
317  output.writeData(lineOutput[iband],y,iband);
318  }
319  catch(string errorstring){
320  cerr << errorstring << "in band " << iband << ", line " << y << endl;
321  }
322  }
323  progress=(1.0+y)/output.nrOfRow();
324  pfnProgress(progress,pszMessage,pProgressArg);
325  }
326 }
327 
328 void filter::Filter::smooth(ImgReaderGdal& input, ImgWriterGdal& output, short dim)
329 {
330  assert(dim>0);
331  m_taps.resize(dim);
332  for(int itap=0;itap<dim;++itap)
333  m_taps[itap]=1.0/dim;
334  filter(input,output);
335 }
336 
337 void filter::Filter::filter(ImgReaderGdal& input, ImgWriterGdal& output)
338 {
339  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
340  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
341  const char* pszMessage;
342  void* pProgressArg=NULL;
343  GDALProgressFunc pfnProgress=GDALTermProgress;
344  double progress=0;
345  pfnProgress(progress,pszMessage,pProgressArg);
346  for(int y=0;y<input.nrOfRow();++y){
347  for(int iband=0;iband<input.nrOfBand();++iband)
348  input.readData(lineInput[iband],y,iband);
349  vector<double> pixelInput(input.nrOfBand());
350  vector<double> pixelOutput(input.nrOfBand());
351  for(int x=0;x<input.nrOfCol();++x){
352  pixelInput=lineInput.selectCol(x);
353  filter(pixelInput,pixelOutput);
354  for(int iband=0;iband<input.nrOfBand();++iband)
355  lineOutput[iband][x]=pixelOutput[iband];
356  }
357  for(int iband=0;iband<input.nrOfBand();++iband){
358  try{
359  output.writeData(lineOutput[iband],y,iband);
360  }
361  catch(string errorstring){
362  cerr << errorstring << "in band " << iband << ", line " << y << endl;
363  }
364  }
365  progress=(1.0+y)/output.nrOfRow();
366  pfnProgress(progress,pszMessage,pProgressArg);
367  }
368 }
369 
370 void filter::Filter::stat(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method)
371 {
372  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
373  assert(output.nrOfCol()==input.nrOfCol());
374  vector<double> lineOutput(output.nrOfCol());
376  stat.setNoDataValues(m_noDataValues);
377  const char* pszMessage;
378  void* pProgressArg=NULL;
379  GDALProgressFunc pfnProgress=GDALTermProgress;
380  double progress=0;
381  pfnProgress(progress,pszMessage,pProgressArg);
382  for(int y=0;y<input.nrOfRow();++y){
383  for(int iband=0;iband<input.nrOfBand();++iband)
384  input.readData(lineInput[iband],y,iband);
385  vector<double> pixelInput(input.nrOfBand());
386  for(int x=0;x<input.nrOfCol();++x){
387  pixelInput=lineInput.selectCol(x);
388  switch(getFilterType(method)){
389  case(filter::median):
390  lineOutput[x]=stat.median(pixelInput);
391  break;
392  case(filter::min):
393  lineOutput[x]=stat.mymin(pixelInput);
394  break;
395  case(filter::max):
396  lineOutput[x]=stat.mymax(pixelInput);
397  break;
398  case(filter::sum):
399  lineOutput[x]=stat.sum(pixelInput);
400  break;
401  case(filter::var):
402  lineOutput[x]=stat.var(pixelInput);
403  break;
404  case(filter::stdev):
405  lineOutput[x]=sqrt(stat.var(pixelInput));
406  break;
407  case(filter::mean):
408  lineOutput[x]=stat.mean(pixelInput);
409  break;
410  case(filter::percentile):
411  assert(m_threshold.size());
412  lineOutput[x]=stat.percentile(pixelInput,pixelInput.begin(),pixelInput.end(),m_threshold[0]);
413  break;
414  default:
415  std::string errorString="method not supported";
416  throw(errorString);
417  break;
418  }
419  }
420  try{
421  output.writeData(lineOutput,y);
422  }
423  catch(string errorstring){
424  cerr << errorstring << "in line " << y << endl;
425  }
426  progress=(1.0+y)/output.nrOfRow();
427  pfnProgress(progress,pszMessage,pProgressArg);
428  }
429 }
430 
431 void filter::Filter::stats(ImgReaderGdal& input, ImgWriterGdal& output, const vector<std::string>& methods)
432 {
433  assert(output.nrOfBand()==methods.size());
434  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
435  assert(output.nrOfCol()==input.nrOfCol());
436  Vector2d<double> lineOutput(methods.size(),output.nrOfCol());
438  stat.setNoDataValues(m_noDataValues);
439  const char* pszMessage;
440  void* pProgressArg=NULL;
441  GDALProgressFunc pfnProgress=GDALTermProgress;
442  double progress=0;
443  pfnProgress(progress,pszMessage,pProgressArg);
444  for(int y=0;y<input.nrOfRow();++y){
445  for(int iband=0;iband<input.nrOfBand();++iband)
446  input.readData(lineInput[iband],y,iband);
447  vector<double> pixelInput(input.nrOfBand());
448  for(int x=0;x<input.nrOfCol();++x){
449  try{
450  pixelInput=lineInput.selectCol(x);
451  }
452  catch(...){
453  std::cout << "error is caught" << std::endl;
454  exit(1);
455  }
456  int ithreshold=0;//threshold to use for percentiles
457  try{
458  for(int imethod=0;imethod<methods.size();++imethod){
459  switch(getFilterType(methods[imethod])){
460  case(filter::nvalid):
461  lineOutput[imethod][x]=stat.nvalid(pixelInput);
462  break;
463  case(filter::median):
464  lineOutput[imethod][x]=stat.median(pixelInput);
465  break;
466  case(filter::min):
467  lineOutput[imethod][x]=stat.mymin(pixelInput);
468  break;
469  case(filter::max):
470  lineOutput[imethod][x]=stat.mymax(pixelInput);
471  break;
472  case(filter::sum):
473  lineOutput[imethod][x]=stat.sum(pixelInput);
474  break;
475  case(filter::var):
476  lineOutput[imethod][x]=stat.var(pixelInput);
477  break;
478  case(filter::stdev):
479  lineOutput[imethod][x]=sqrt(stat.var(pixelInput));
480  break;
481  case(filter::mean):
482  lineOutput[imethod][x]=stat.mean(pixelInput);
483  break;
484  case(filter::percentile):{
485  assert(m_threshold.size());
486  double threshold=(ithreshold<m_threshold.size())? m_threshold[ithreshold] : m_threshold[0];
487  lineOutput[imethod][x]=stat.percentile(pixelInput,pixelInput.begin(),pixelInput.end(),threshold);
488  ++ithreshold;
489  break;
490  }
491  default:
492  std::string errorString="method not supported";
493  throw(errorString);
494  break;
495  }
496  }
497  }
498  catch(...){
499  cerr << "An Error in line " << y << endl;
500  }
501  }
502  for(int imethod=0;imethod<methods.size();++imethod)
503  output.writeData(lineOutput[imethod],y,imethod);
504  progress=(1.0+y)/output.nrOfRow();
505  pfnProgress(progress,pszMessage,pProgressArg);
506  }
507 }
508 
509 void filter::Filter::filter(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim)
510 {
511  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
512  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());;
513  const char* pszMessage;
514  void* pProgressArg=NULL;
515  GDALProgressFunc pfnProgress=GDALTermProgress;
516  double progress=0;
517  pfnProgress(progress,pszMessage,pProgressArg);
518  for(int y=0;y<input.nrOfRow();++y){
519  for(int iband=0;iband<input.nrOfBand();++iband)
520  input.readData(lineInput[iband],y,iband);
521  vector<double> pixelInput(input.nrOfBand());
522  vector<double> pixelOutput;
523  for(int x=0;x<input.nrOfCol();++x){
524  pixelInput=lineInput.selectCol(x);
525  filter(pixelInput,pixelOutput,method,dim);
526  for(int iband=0;iband<pixelOutput.size();++iband){
527  lineOutput[iband][x]=pixelOutput[iband];
528  // if(pixelInput[iband]!=0)
529  // assert(pixelOutput[iband]!=0);
530  }
531  }
532  for(int iband=0;iband<input.nrOfBand();++iband){
533  try{
534  output.writeData(lineOutput[iband],y,iband);
535  }
536  catch(string errorstring){
537  cerr << errorstring << "in band " << iband << ", line " << y << endl;
538  }
539  }
540  progress=(1.0+y)/output.nrOfRow();
541  pfnProgress(progress,pszMessage,pProgressArg);
542  }
543 }
544 
545 void filter::Filter::getSavGolayCoefficients(vector<double> &tapz, int np, int nl, int nr, int ld, int m) {
546  int j, k, imj, ipj, kk, mm;
547  double d, fac, sum;
548 
549  // c.resize(nl+1+nr);
550  vector<double> tmpc(np);
551  if(np < nl + nr + 1 || nl < 0 || nr < 0 || ld > m || nl + nr < m) {
552  cerr << "bad args in savgol" << endl;
553  return;
554  }
555  vector<int> indx(m + 1, 0);
556  vector<double> a((m + 1) * (m + 1), 0.0);
557  vector<double> b(m + 1, 0.0);
558 
559  for(ipj = 0; ipj <= (m << 1); ++ipj) {
560  sum = (ipj ? 0.0 : 1.0);
561  for(k = 1; k <= nr; ++k)
562  sum += (int)pow((double)k, (double)ipj);
563  for(k = 1; k <= nl; ++k)
564  sum += (int)pow((double) - k, (double)ipj);
565  mm = (ipj < 2 * m - ipj ? ipj : 2 * m - ipj);
566  for(imj = -mm; imj <= mm; imj += 2)
567  a[(ipj + imj) / 2 * (m + 1) + (ipj - imj) / 2] = sum;
568  }
569  ludcmp(a, indx, d);
570 
571  for(j = 0; j < m + 1; ++j)
572  b[j] = 0.0;
573  b[ld] = 1.0;
574 
575  lubksb(a, indx, b);
576  // for(kk = 0; kk < np; ++kk)
577  // c[kk] = 0.0;
578  for(k = -nl; k <= nr; ++k) {
579  // for(k = -nl; k < nr; ++k) {
580  sum = b[0];
581  fac = 1.0;
582  for(mm = 1; mm <= m; ++mm)
583  sum += b[mm] * (fac *= k);
584  // store in wrap=around order
585  kk = (np - k) % np;
586  //re-order c as I need for taps
587  // kk=k+nl;
588  tmpc[kk] = sum;
589  }
590  tapz.resize(nl+1+nr);
591  // for(k=0;k<nl+1+nr)
592  tapz[tapz.size()/2]=tmpc[0];
593  //past data points
594  for(k=1;k<=tapz.size()/2;++k)
595  tapz[tapz.size()/2-k]=tmpc[k];
596  //future data points
597  for(k=1;k<=tapz.size()/2;++k)
598  tapz[tapz.size()/2+k]=tmpc[np-k];
599 }
600 
601 void filter::Filter::ludcmp(vector<double> &a, vector<int> &indx, double &d) {
602  const double TINY = 1.0e-20;
603  int i, imax = -1, j, k;
604  double big, dum, sum, temp;
605 
606  int n = indx.size();
607  vector<double> vv(n, 0.0);
608 
609  d = 1.0;
610  for(i = 0; i < n; ++i) {
611  big = 0.0;
612  for(j = 0; j < n; ++j)
613  if((temp = fabs(a[i * n + j])) > big) big = temp;
614 
615  if(big == 0.0) {
616  cerr << "Singular matrix in routine ludcmp" << endl;
617  return;
618  }
619  vv[i] = 1. / big;
620  }
621 
622  for(j = 0; j < n; ++j) {
623  for(i = 0; i < j; ++i) {
624  sum = a[i * n + j];
625  for(k = 0; k < i; ++k)
626  sum -= a[i * n + k] * a[k * n + j];
627  a[i * n + j] = sum;
628  }
629  big = 0.0;
630  for(i = j; i < n; ++i) {
631  sum = a[i * n + j];
632  for(k = 0; k < j; ++k)
633  sum -= a[i * n + k] * a[k * n + j];
634  a[i * n + j] = sum;
635  if((dum = vv[i] * fabs(sum)) >= big) {
636  big = dum;
637  imax = i;
638  }
639  }
640 
641  if(j != imax) {
642  for(k = 0; k < n; ++k) {
643  dum = a[imax * n + k];
644  a[imax * n + k] = a[j * n + k];
645  a[j * n + k] = dum;
646  }
647  d = -d;
648  vv[imax] = vv[j];
649  }
650  indx[j] = imax;
651  if(a[j * n + j] == 0.0) a[j * n + j] = TINY;
652  if(j != n - 1) {
653  dum = 1. / a[j * n + j];
654  for(i = j + 1; i < n; ++i)
655  a[i * n + j] *= dum;
656  }
657  }
658 }
659 
660 void filter::Filter::lubksb(vector<double> &a, vector<int> &indx, vector<double> &b) {
661  int i, ii = 0, ip, j;
662  double sum;
663  int n = indx.size();
664 
665  for(i = 0; i < n; ++i) {
666  ip = indx[i];
667  sum = b[ip];
668  b[ip] = b[i];
669  if(ii != 0)
670  for(j = ii - 1; j < i; ++j)
671  sum -= a[i * n + j] * b[j];
672  else if(sum != 0.0)
673  ii = i + 1;
674  b[i] = sum;
675  }
676  for(i = n - 1; i >= 0; --i) {
677  sum = b[i];
678  for(j = i + 1; j < n; ++j)
679  sum -= a[i * n + j] * b[j];
680  b[i] = sum / a[i * n + i];
681  }
682 }
683 
684 double filter::Filter::getCentreWavelength(const std::vector<double> &wavelengthIn, const Vector2d<double>& srf, const std::string& interpolationType, double delta, bool verbose)
685 {
686  assert(srf.size()==2);//[0]: wavelength, [1]: response function
687  int nband=srf[0].size();
688  double start=floor(wavelengthIn[0]);
689  double end=ceil(wavelengthIn.back());
690  if(verbose)
691  std::cout << "wavelengths in [" << start << "," << end << "]" << std::endl << std::flush;
692 
694 
695  gsl_interp_accel *acc;
696  stat.allocAcc(acc);
697  gsl_spline *spline;
698  stat.getSpline(interpolationType,nband,spline);
699  stat.initSpline(spline,&(srf[0][0]),&(srf[1][0]),nband);
700  if(verbose)
701  std::cout << "calculating norm of srf" << std::endl << std::flush;
702  double norm=0;
703  norm=gsl_spline_eval_integ(spline,srf[0].front(),srf[0].back(),acc);
704  if(verbose)
705  std::cout << "norm of srf: " << norm << std::endl << std::flush;
706  gsl_spline_free(spline);
707  gsl_interp_accel_free(acc);
708  std::vector<double> wavelength_fine;
709  for(double win=floor(wavelengthIn[0]);win<=ceil(wavelengthIn.back());win+=delta)
710  wavelength_fine.push_back(win);
711 
712  if(verbose)
713  std::cout << "interpolate wavelengths to " << wavelength_fine.size() << " entries " << std::endl;
714  std::vector<double> srf_fine;//spectral response function, interpolated for wavelength_fine
715 
716  stat.interpolateUp(srf[0],srf[1],wavelength_fine,interpolationType,srf_fine,verbose);
717  assert(srf_fine.size()==wavelength_fine.size());
718 
719  gsl_interp_accel *accOut;
720  stat.allocAcc(accOut);
721  gsl_spline *splineOut;
722  stat.getSpline(interpolationType,wavelength_fine.size(),splineOut);
723  assert(splineOut);
724 
725  std::vector<double> wavelengthOut(wavelength_fine.size());
726 
727  for(int iband=0;iband<wavelengthOut.size();++iband)
728  wavelengthOut[iband]=wavelength_fine[iband]*srf_fine[iband];
729 
730  stat.initSpline(splineOut,&(wavelength_fine[0]),&(wavelengthOut[0]),wavelength_fine.size());
731  double centreWavelength=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
732 
733  gsl_spline_free(splineOut);
734  gsl_interp_accel_free(accOut);
735 
736  return(centreWavelength);
737 }
738 
739 // 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){
740 // Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
741 // Vector2d<double> lineOutput(wavelengthOut.size(),input.nrOfCol());
742 // const char* pszMessage;
743 // void* pProgressArg=NULL;
744 // GDALProgressFunc pfnProgress=GDALTermProgress;
745 // double progress=0;
746 // pfnProgress(progress,pszMessage,pProgressArg);
747 // for(int y=0;y<input.nrOfRow();++y){
748 // for(int iband=0;iband<input.nrOfBand();++iband)
749 // input.readData(lineInput[iband],GDT_Float64,y,iband);
750 // applyFwhm<double>(wavelengthIn,lineInput,wavelengthOut,fwhm, interpolationType, lineOutput, verbose);
751 // for(int iband=0;iband<output.nrOfBand();++iband){
752 // try{
753 // output.writeData(lineOutput[iband],GDT_Float64,y,iband);
754 // }
755 // catch(string errorstring){
756 // cerr << errorstring << "in band " << iband << ", line " << y << endl;
757 // }
758 // }
759 // progress=(1.0+y)/output.nrOfRow();
760 // pfnProgress(progress,pszMessage,pProgressArg);
761 // }
762 // }
763 
764 // void filter::Filter::applySrf(const vector<double> &wavelengthIn, const ImgReaderGdal& input, const vector< Vector2d<double> > &srf, const std::string& interpolationType, ImgWriterGdal& output, bool verbose){
765 // assert(output.nrOfBand()==srf.size());
766 // double centreWavelength=0;
767 // Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
768 // const char* pszMessage;
769 // void* pProgressArg=NULL;
770 // GDALProgressFunc pfnProgress=GDALTermProgress;
771 // double progress=0;
772 // pfnProgress(progress,pszMessage,pProgressArg);
773 // for(int y=0;y<input.nrOfRow();++y){
774 // for(int iband=0;iband<input.nrOfBand();++iband)
775 // input.readData(lineInput[iband],GDT_Float64,y,iband);
776 // for(int isrf=0;isrf<srf.size();++isrf){
777 // vector<double> lineOutput(input.nrOfCol());
778 // centreWavelength=applySrf<double>(wavelengthIn,lineInput,srf[isrf], interpolationType, lineOutput, verbose);
779 // for(int iband=0;iband<output.nrOfBand();++iband){
780 // try{
781 // output.writeData(lineOutput,GDT_Float64,y,isrf);
782 // }
783 // catch(string errorstring){
784 // cerr << errorstring << "in band " << iband << ", line " << y << endl;
785 // }
786 // }
787 // }
788 // progress=(1.0+y)/output.nrOfRow();
789 // pfnProgress(progress,pszMessage,pProgressArg);
790 // }
791 // }
STL namespace.
Definition: Filter.h:33
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
void readData(T &value, int col, int row, int band=0)
Read a single pixel cell value at a specific column and row for a specific band (all indices start co...
Definition: ImgReaderGdal.h:95
int nrOfRow(void) const
Get the number of rows of this dataset.
bool writeData(T &value, int col, int row, int band=0)
Write a single pixel cell value at a specific column and row for a specific band (all indices start c...
Definition: ImgWriterGdal.h:96
int nrOfBand(void) const
Get the number of bands of this dataset.