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],GDT_Float64,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],GDT_Float64,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],GDT_Float64,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],GDT_Float64,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],GDT_Float64,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],GDT_Float64,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],GDT_Float64,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],GDT_Float64,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],GDT_Float64,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],GDT_Float64,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],GDT_Float64,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],GDT_Float64,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],GDT_Float64,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],GDT_Float64,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],GDT_Float64,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,GDT_Float64,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],GDT_Float64,y,iband);
447  vector<double> pixelInput(input.nrOfBand());
448  for(int x=0;x<input.nrOfCol();++x){
449  pixelInput=lineInput.selectCol(x);
450  int ithreshold=0;//threshold to use for percentiles
451  for(int imethod=0;imethod<methods.size();++imethod){
452  switch(getFilterType(methods[imethod])){
453  case(filter::nvalid):
454  lineOutput[imethod][x]=stat.nvalid(pixelInput);
455  break;
456  case(filter::median):
457  lineOutput[imethod][x]=stat.median(pixelInput);
458  break;
459  case(filter::min):
460  lineOutput[imethod][x]=stat.mymin(pixelInput);
461  break;
462  case(filter::max):
463  lineOutput[imethod][x]=stat.mymax(pixelInput);
464  break;
465  case(filter::sum):
466  lineOutput[imethod][x]=stat.sum(pixelInput);
467  break;
468  case(filter::var):
469  lineOutput[imethod][x]=stat.var(pixelInput);
470  break;
471  case(filter::stdev):
472  lineOutput[imethod][x]=sqrt(stat.var(pixelInput));
473  break;
474  case(filter::mean):
475  lineOutput[imethod][x]=stat.mean(pixelInput);
476  break;
477  case(filter::percentile):{
478  assert(m_threshold.size());
479  double threshold=(ithreshold<m_threshold.size())? m_threshold[ithreshold] : m_threshold[0];
480  lineOutput[imethod][x]=stat.percentile(pixelInput,pixelInput.begin(),pixelInput.end(),threshold);
481  ++ithreshold;
482  break;
483  }
484  default:
485  std::string errorString="method not supported";
486  throw(errorString);
487  break;
488  }
489  }
490  }
491  for(int imethod=0;imethod<methods.size();++imethod){
492  try{
493  output.writeData(lineOutput[imethod],GDT_Float64,y,imethod);
494  }
495  catch(string errorstring){
496  cerr << errorstring << "in line " << y << endl;
497  }
498  }
499  progress=(1.0+y)/output.nrOfRow();
500  pfnProgress(progress,pszMessage,pProgressArg);
501  }
502 }
503 
504 void filter::Filter::filter(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim)
505 {
506  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
507  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());;
508  const char* pszMessage;
509  void* pProgressArg=NULL;
510  GDALProgressFunc pfnProgress=GDALTermProgress;
511  double progress=0;
512  pfnProgress(progress,pszMessage,pProgressArg);
513  for(int y=0;y<input.nrOfRow();++y){
514  for(int iband=0;iband<input.nrOfBand();++iband)
515  input.readData(lineInput[iband],GDT_Float64,y,iband);
516  vector<double> pixelInput(input.nrOfBand());
517  vector<double> pixelOutput;
518  for(int x=0;x<input.nrOfCol();++x){
519  pixelInput=lineInput.selectCol(x);
520  filter(pixelInput,pixelOutput,method,dim);
521  for(int iband=0;iband<pixelOutput.size();++iband){
522  lineOutput[iband][x]=pixelOutput[iband];
523  // if(pixelInput[iband]!=0)
524  // assert(pixelOutput[iband]!=0);
525  }
526  }
527  for(int iband=0;iband<input.nrOfBand();++iband){
528  try{
529  output.writeData(lineOutput[iband],GDT_Float64,y,iband);
530  }
531  catch(string errorstring){
532  cerr << errorstring << "in band " << iband << ", line " << y << endl;
533  }
534  }
535  progress=(1.0+y)/output.nrOfRow();
536  pfnProgress(progress,pszMessage,pProgressArg);
537  }
538 }
539 
540 void filter::Filter::getSavGolayCoefficients(vector<double> &tapz, int np, int nl, int nr, int ld, int m) {
541  int j, k, imj, ipj, kk, mm;
542  double d, fac, sum;
543 
544  // c.resize(nl+1+nr);
545  vector<double> tmpc(np);
546  if(np < nl + nr + 1 || nl < 0 || nr < 0 || ld > m || nl + nr < m) {
547  cerr << "bad args in savgol" << endl;
548  return;
549  }
550  vector<int> indx(m + 1, 0);
551  vector<double> a((m + 1) * (m + 1), 0.0);
552  vector<double> b(m + 1, 0.0);
553 
554  for(ipj = 0; ipj <= (m << 1); ++ipj) {
555  sum = (ipj ? 0.0 : 1.0);
556  for(k = 1; k <= nr; ++k)
557  sum += (int)pow((double)k, (double)ipj);
558  for(k = 1; k <= nl; ++k)
559  sum += (int)pow((double) - k, (double)ipj);
560  mm = (ipj < 2 * m - ipj ? ipj : 2 * m - ipj);
561  for(imj = -mm; imj <= mm; imj += 2)
562  a[(ipj + imj) / 2 * (m + 1) + (ipj - imj) / 2] = sum;
563  }
564  ludcmp(a, indx, d);
565 
566  for(j = 0; j < m + 1; ++j)
567  b[j] = 0.0;
568  b[ld] = 1.0;
569 
570  lubksb(a, indx, b);
571  // for(kk = 0; kk < np; ++kk)
572  // c[kk] = 0.0;
573  for(k = -nl; k <= nr; ++k) {
574  // for(k = -nl; k < nr; ++k) {
575  sum = b[0];
576  fac = 1.0;
577  for(mm = 1; mm <= m; ++mm)
578  sum += b[mm] * (fac *= k);
579  // store in wrap=around order
580  kk = (np - k) % np;
581  //re-order c as I need for taps
582  // kk=k+nl;
583  tmpc[kk] = sum;
584  }
585  tapz.resize(nl+1+nr);
586  // for(k=0;k<nl+1+nr)
587  tapz[tapz.size()/2]=tmpc[0];
588  //past data points
589  for(k=1;k<=tapz.size()/2;++k)
590  tapz[tapz.size()/2-k]=tmpc[k];
591  //future data points
592  for(k=1;k<=tapz.size()/2;++k)
593  tapz[tapz.size()/2+k]=tmpc[np-k];
594 }
595 
596 void filter::Filter::ludcmp(vector<double> &a, vector<int> &indx, double &d) {
597  const double TINY = 1.0e-20;
598  int i, imax = -1, j, k;
599  double big, dum, sum, temp;
600 
601  int n = indx.size();
602  vector<double> vv(n, 0.0);
603 
604  d = 1.0;
605  for(i = 0; i < n; ++i) {
606  big = 0.0;
607  for(j = 0; j < n; ++j)
608  if((temp = fabs(a[i * n + j])) > big) big = temp;
609 
610  if(big == 0.0) {
611  cerr << "Singular matrix in routine ludcmp" << endl;
612  return;
613  }
614  vv[i] = 1. / big;
615  }
616 
617  for(j = 0; j < n; ++j) {
618  for(i = 0; i < j; ++i) {
619  sum = a[i * n + j];
620  for(k = 0; k < i; ++k)
621  sum -= a[i * n + k] * a[k * n + j];
622  a[i * n + j] = sum;
623  }
624  big = 0.0;
625  for(i = j; i < n; ++i) {
626  sum = a[i * n + j];
627  for(k = 0; k < j; ++k)
628  sum -= a[i * n + k] * a[k * n + j];
629  a[i * n + j] = sum;
630  if((dum = vv[i] * fabs(sum)) >= big) {
631  big = dum;
632  imax = i;
633  }
634  }
635 
636  if(j != imax) {
637  for(k = 0; k < n; ++k) {
638  dum = a[imax * n + k];
639  a[imax * n + k] = a[j * n + k];
640  a[j * n + k] = dum;
641  }
642  d = -d;
643  vv[imax] = vv[j];
644  }
645  indx[j] = imax;
646  if(a[j * n + j] == 0.0) a[j * n + j] = TINY;
647  if(j != n - 1) {
648  dum = 1. / a[j * n + j];
649  for(i = j + 1; i < n; ++i)
650  a[i * n + j] *= dum;
651  }
652  }
653 }
654 
655 void filter::Filter::lubksb(vector<double> &a, vector<int> &indx, vector<double> &b) {
656  int i, ii = 0, ip, j;
657  double sum;
658  int n = indx.size();
659 
660  for(i = 0; i < n; ++i) {
661  ip = indx[i];
662  sum = b[ip];
663  b[ip] = b[i];
664  if(ii != 0)
665  for(j = ii - 1; j < i; ++j)
666  sum -= a[i * n + j] * b[j];
667  else if(sum != 0.0)
668  ii = i + 1;
669  b[i] = sum;
670  }
671  for(i = n - 1; i >= 0; --i) {
672  sum = b[i];
673  for(j = i + 1; j < n; ++j)
674  sum -= a[i * n + j] * b[j];
675  b[i] = sum / a[i * n + i];
676  }
677 }
678 
679 double filter::Filter::getCentreWavelength(const std::vector<double> &wavelengthIn, const Vector2d<double>& srf, const std::string& interpolationType, double delta, bool verbose)
680 {
681  assert(srf.size()==2);//[0]: wavelength, [1]: response function
682  int nband=srf[0].size();
683  double start=floor(wavelengthIn[0]);
684  double end=ceil(wavelengthIn.back());
685  if(verbose)
686  std::cout << "wavelengths in [" << start << "," << end << "]" << std::endl << std::flush;
687 
689 
690  gsl_interp_accel *acc;
691  stat.allocAcc(acc);
692  gsl_spline *spline;
693  stat.getSpline(interpolationType,nband,spline);
694  stat.initSpline(spline,&(srf[0][0]),&(srf[1][0]),nband);
695  if(verbose)
696  std::cout << "calculating norm of srf" << std::endl << std::flush;
697  double norm=0;
698  norm=gsl_spline_eval_integ(spline,srf[0].front(),srf[0].back(),acc);
699  if(verbose)
700  std::cout << "norm of srf: " << norm << std::endl << std::flush;
701  gsl_spline_free(spline);
702  gsl_interp_accel_free(acc);
703  std::vector<double> wavelength_fine;
704  for(double win=floor(wavelengthIn[0]);win<=ceil(wavelengthIn.back());win+=delta)
705  wavelength_fine.push_back(win);
706 
707  if(verbose)
708  std::cout << "interpolate wavelengths to " << wavelength_fine.size() << " entries " << std::endl;
709  std::vector<double> srf_fine;//spectral response function, interpolated for wavelength_fine
710 
711  stat.interpolateUp(srf[0],srf[1],wavelength_fine,interpolationType,srf_fine,verbose);
712  assert(srf_fine.size()==wavelength_fine.size());
713 
714  gsl_interp_accel *accOut;
715  stat.allocAcc(accOut);
716  gsl_spline *splineOut;
717  stat.getSpline(interpolationType,wavelength_fine.size(),splineOut);
718  assert(splineOut);
719 
720  std::vector<double> wavelengthOut(wavelength_fine.size());
721 
722  for(int iband=0;iband<wavelengthOut.size();++iband)
723  wavelengthOut[iband]=wavelength_fine[iband]*srf_fine[iband];
724 
725  stat.initSpline(splineOut,&(wavelength_fine[0]),&(wavelengthOut[0]),wavelength_fine.size());
726  double centreWavelength=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
727 
728  gsl_spline_free(splineOut);
729  gsl_interp_accel_free(accOut);
730 
731  return(centreWavelength);
732 }
733 
734 // 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){
735 // Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
736 // Vector2d<double> lineOutput(wavelengthOut.size(),input.nrOfCol());
737 // const char* pszMessage;
738 // void* pProgressArg=NULL;
739 // GDALProgressFunc pfnProgress=GDALTermProgress;
740 // double progress=0;
741 // pfnProgress(progress,pszMessage,pProgressArg);
742 // for(int y=0;y<input.nrOfRow();++y){
743 // for(int iband=0;iband<input.nrOfBand();++iband)
744 // input.readData(lineInput[iband],GDT_Float64,y,iband);
745 // applyFwhm<double>(wavelengthIn,lineInput,wavelengthOut,fwhm, interpolationType, lineOutput, verbose);
746 // for(int iband=0;iband<output.nrOfBand();++iband){
747 // try{
748 // output.writeData(lineOutput[iband],GDT_Float64,y,iband);
749 // }
750 // catch(string errorstring){
751 // cerr << errorstring << "in band " << iband << ", line " << y << endl;
752 // }
753 // }
754 // progress=(1.0+y)/output.nrOfRow();
755 // pfnProgress(progress,pszMessage,pProgressArg);
756 // }
757 // }
758 
759 // void filter::Filter::applySrf(const vector<double> &wavelengthIn, const ImgReaderGdal& input, const vector< Vector2d<double> > &srf, const std::string& interpolationType, ImgWriterGdal& output, bool verbose){
760 // assert(output.nrOfBand()==srf.size());
761 // double centreWavelength=0;
762 // Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
763 // const char* pszMessage;
764 // void* pProgressArg=NULL;
765 // GDALProgressFunc pfnProgress=GDALTermProgress;
766 // double progress=0;
767 // pfnProgress(progress,pszMessage,pProgressArg);
768 // for(int y=0;y<input.nrOfRow();++y){
769 // for(int iband=0;iband<input.nrOfBand();++iband)
770 // input.readData(lineInput[iband],GDT_Float64,y,iband);
771 // for(int isrf=0;isrf<srf.size();++isrf){
772 // vector<double> lineOutput(input.nrOfCol());
773 // centreWavelength=applySrf<double>(wavelengthIn,lineInput,srf[isrf], interpolationType, lineOutput, verbose);
774 // for(int iband=0;iband<output.nrOfBand();++iband){
775 // try{
776 // output.writeData(lineOutput,GDT_Float64,y,isrf);
777 // }
778 // catch(string errorstring){
779 // cerr << errorstring << "in band " << iband << ", line " << y << endl;
780 // }
781 // }
782 // }
783 // progress=(1.0+y)/output.nrOfRow();
784 // pfnProgress(progress,pszMessage,pProgressArg);
785 // }
786 // }
int nrOfBand(void) const
Get the number of bands of this dataset.
STL namespace.
Definition: Filter.h:33
int nrOfRow(void) const
Get the number of rows of this dataset.
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
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 nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98