Processing Kernel for remote sensing 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 void filter::Filter::dwtForward(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family){
55  const char* pszMessage;
56  void* pProgressArg=NULL;
57  GDALProgressFunc pfnProgress=GDALTermProgress;
58  double progress=0;
59  pfnProgress(progress,pszMessage,pProgressArg);
60  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
61  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
62  for(int y=0;y<input.nrOfRow();++y){
63  for(int iband=0;iband<input.nrOfBand();++iband)
64  input.readData(lineInput[iband],GDT_Float64,y,iband);
65  vector<double> pixelInput(input.nrOfBand());
66  for(int x=0;x<input.nrOfCol();++x){
67  pixelInput=lineInput.selectCol(x);
68  dwtForward(pixelInput,wavelet_type,family);
69  for(int iband=0;iband<input.nrOfBand();++iband)
70  lineOutput[iband][x]=pixelInput[iband];
71  }
72  for(int iband=0;iband<input.nrOfBand();++iband){
73  try{
74  output.writeData(lineOutput[iband],GDT_Float64,y,iband);
75  }
76  catch(string errorstring){
77  cerr << errorstring << "in band " << iband << ", line " << y << endl;
78  }
79  }
80  progress=(1.0+y)/output.nrOfRow();
81  pfnProgress(progress,pszMessage,pProgressArg);
82  }
83 }
84 
85 void filter::Filter::dwtInverse(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family){
86  const char* pszMessage;
87  void* pProgressArg=NULL;
88  GDALProgressFunc pfnProgress=GDALTermProgress;
89  double progress=0;
90  pfnProgress(progress,pszMessage,pProgressArg);
91  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
92  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
93  for(int y=0;y<input.nrOfRow();++y){
94  for(int iband=0;iband<input.nrOfBand();++iband)
95  input.readData(lineInput[iband],GDT_Float64,y,iband);
96  vector<double> pixelInput(input.nrOfBand());
97  for(int x=0;x<input.nrOfCol();++x){
98  pixelInput=lineInput.selectCol(x);
99  dwtInverse(pixelInput,wavelet_type,family);
100  for(int iband=0;iband<input.nrOfBand();++iband)
101  lineOutput[iband][x]=pixelInput[iband];
102  }
103  for(int iband=0;iband<input.nrOfBand();++iband){
104  try{
105  output.writeData(lineOutput[iband],GDT_Float64,y,iband);
106  }
107  catch(string errorstring){
108  cerr << errorstring << "in band " << iband << ", line " << y << endl;
109  }
110  }
111  progress=(1.0+y)/output.nrOfRow();
112  pfnProgress(progress,pszMessage,pProgressArg);
113  }
114 }
115 
116 void filter::Filter::dwtCut(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, double cut){
117  const char* pszMessage;
118  void* pProgressArg=NULL;
119  GDALProgressFunc pfnProgress=GDALTermProgress;
120  double progress=0;
121  pfnProgress(progress,pszMessage,pProgressArg);
122  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
123  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
124  for(int y=0;y<input.nrOfRow();++y){
125  for(int iband=0;iband<input.nrOfBand();++iband)
126  input.readData(lineInput[iband],GDT_Float64,y,iband);
127  vector<double> pixelInput(input.nrOfBand());
128  for(int x=0;x<input.nrOfCol();++x){
129  pixelInput=lineInput.selectCol(x);
130  dwtCut(pixelInput,wavelet_type,family,cut);
131  for(int iband=0;iband<input.nrOfBand();++iband)
132  lineOutput[iband][x]=pixelInput[iband];
133  }
134  for(int iband=0;iband<input.nrOfBand();++iband){
135  try{
136  output.writeData(lineOutput[iband],GDT_Float64,y,iband);
137  }
138  catch(string errorstring){
139  cerr << errorstring << "in band " << iband << ", line " << y << endl;
140  }
141  }
142  progress=(1.0+y)/output.nrOfRow();
143  pfnProgress(progress,pszMessage,pProgressArg);
144  }
145 }
146 
147 void filter::Filter::dwtCutFrom(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, int band){
148  const char* pszMessage;
149  void* pProgressArg=NULL;
150  GDALProgressFunc pfnProgress=GDALTermProgress;
151  double progress=0;
152  pfnProgress(progress,pszMessage,pProgressArg);
153  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
154  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
155  for(int y=0;y<input.nrOfRow();++y){
156  for(int iband=0;iband<input.nrOfBand();++iband)
157  input.readData(lineInput[iband],GDT_Float64,y,iband);
158  vector<double> pixelInput(input.nrOfBand());
159  for(int x=0;x<input.nrOfCol();++x){
160  pixelInput=lineInput.selectCol(x);
161  dwtForward(pixelInput,wavelet_type,family);
162  for(int iband=0;iband<input.nrOfBand();++iband){
163  if(iband>=band)
164  pixelInput[iband]=0;
165  }
166  dwtInverse(pixelInput,wavelet_type,family);
167  for(int iband=0;iband<input.nrOfBand();++iband)
168  lineOutput[iband][x]=pixelInput[iband];
169  }
170  for(int iband=0;iband<input.nrOfBand();++iband){
171  try{
172  output.writeData(lineOutput[iband],GDT_Float64,y,iband);
173  }
174  catch(string errorstring){
175  cerr << errorstring << "in band " << iband << ", line " << y << endl;
176  }
177  }
178  progress=(1.0+y)/output.nrOfRow();
179  pfnProgress(progress,pszMessage,pProgressArg);
180  }
181 }
182 
183 //todo: support different padding strategies
184 void filter::Filter::dwtForward(std::vector<double>& data, const std::string& wavelet_type, int family){
185  int origsize=data.size();
186  //make sure data size if power of 2
187  while(data.size()&(data.size()-1))
188  data.push_back(data.back());
189 
190  int nsize=data.size();
191  gsl_wavelet *w;
192  gsl_wavelet_workspace *work;
193  assert(nsize);
194  w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
195  work=gsl_wavelet_workspace_alloc(nsize);
196  gsl_wavelet_transform_forward(w,&(data[0]),1,nsize,work);
197  data.erase(data.begin()+origsize,data.end());
198  gsl_wavelet_free (w);
199  gsl_wavelet_workspace_free (work);
200 }
201 
202 //todo: support different padding strategies
203 void filter::Filter::dwtInverse(std::vector<double>& data, const std::string& wavelet_type, int family){
204  int origsize=data.size();
205  //make sure data size if power of 2
206  while(data.size()&(data.size()-1))
207  data.push_back(data.back());
208  int nsize=data.size();
209  gsl_wavelet *w;
210  gsl_wavelet_workspace *work;
211  assert(nsize);
212  w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
213  work=gsl_wavelet_workspace_alloc(nsize);
214  gsl_wavelet_transform_inverse(w,&(data[0]),1,nsize,work);
215  data.erase(data.begin()+origsize,data.end());
216  gsl_wavelet_free (w);
217  gsl_wavelet_workspace_free (work);
218 }
219 
220 //todo: support different padding strategies
221 void filter::Filter::dwtCut(std::vector<double>& data, const std::string& wavelet_type, int family, double cut){
222  int origsize=data.size();
223  //make sure data size if power of 2
224  while(data.size()&(data.size()-1))
225  data.push_back(data.back());
226  int nsize=data.size();
227  gsl_wavelet *w;
228  gsl_wavelet_workspace *work;
229  assert(nsize);
230  w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
231  work=gsl_wavelet_workspace_alloc(nsize);
232  gsl_wavelet_transform_forward(w,&(data[0]),1,nsize,work);
233  std::vector<double> abscoeff(data.size());
234  size_t* p=new size_t[data.size()];
235  for(int index=0;index<data.size();++index){
236  abscoeff[index]=fabs(data[index]);
237  }
238  int nc=(100-cut)/100.0*nsize;
239  gsl_sort_index(p,&(abscoeff[0]),1,nsize);
240  for(int i=0;(i+nc)<nsize;i++)
241  data[p[i]]=0;
242  gsl_wavelet_transform_inverse(w,&(data[0]),1,nsize,work);
243  data.erase(data.begin()+origsize,data.end());
244  delete[] p;
245  gsl_wavelet_free (w);
246  gsl_wavelet_workspace_free (work);
247 }
248 
249 void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short verbose)
250 {
251  // bool bverbose=(verbose>1)? true:false;
252  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
253  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
254  const char* pszMessage;
255  void* pProgressArg=NULL;
256  GDALProgressFunc pfnProgress=GDALTermProgress;
257  double progress=0;
258  pfnProgress(progress,pszMessage,pProgressArg);
259  for(int y=0;y<input.nrOfRow();++y){
260  for(int iband=0;iband<input.nrOfBand();++iband)
261  input.readData(lineInput[iband],GDT_Float64,y,iband);
262  vector<double> pixelInput(input.nrOfBand());
263  vector<double> pixelOutput(input.nrOfBand());
264  for(int x=0;x<input.nrOfCol();++x){
265  pixelInput=lineInput.selectCol(x);
266  filter(pixelInput,pixelOutput,method,dim);
267  // morphology(pixelInput,pixelOutput,method,dim,bverbose);
268  for(int iband=0;iband<input.nrOfBand();++iband)
269  lineOutput[iband][x]=pixelOutput[iband];
270  }
271  for(int iband=0;iband<input.nrOfBand();++iband){
272  try{
273  output.writeData(lineOutput[iband],GDT_Float64,y,iband);
274  }
275  catch(string errorstring){
276  cerr << errorstring << "in band " << iband << ", line " << y << endl;
277  }
278  }
279  progress=(1.0+y)/output.nrOfRow();
280  pfnProgress(progress,pszMessage,pProgressArg);
281  }
282 }
283 
284 void filter::Filter::smooth(const ImgReaderGdal& input, ImgWriterGdal& output, short dim)
285 {
286  assert(dim>0);
287  m_taps.resize(dim);
288  for(int itap=0;itap<dim;++itap)
289  m_taps[itap]=1.0/dim;
290  filter(input,output);
291 }
292 
293 // void filter::Filter::smoothnodata(const ImgReaderGdal& input, ImgWriterGdal& output, short dim, short down, int offset)
294 // {
295 // assert(dim>0);
296 // m_taps.resize(dim);
297 // for(int itap=0;itap<dim;++itap)
298 // m_taps[itap]=1.0/dim;
299 // filter(input,output,down,offset);
300 // }
301 
302 void filter::Filter::filter(const ImgReaderGdal& input, ImgWriterGdal& output)
303 {
304  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
305  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());
306  const char* pszMessage;
307  void* pProgressArg=NULL;
308  GDALProgressFunc pfnProgress=GDALTermProgress;
309  double progress=0;
310  pfnProgress(progress,pszMessage,pProgressArg);
311  for(int y=0;y<input.nrOfRow();++y){
312  for(int iband=0;iband<input.nrOfBand();++iband)
313  input.readData(lineInput[iband],GDT_Float64,y,iband);
314  vector<double> pixelInput(input.nrOfBand());
315  vector<double> pixelOutput(input.nrOfBand());
316  for(int x=0;x<input.nrOfCol();++x){
317  pixelInput=lineInput.selectCol(x);
318  filter(pixelInput,pixelOutput);
319  for(int iband=0;iband<input.nrOfBand();++iband)
320  lineOutput[iband][x]=pixelOutput[iband];
321  }
322  for(int iband=0;iband<input.nrOfBand();++iband){
323  try{
324  output.writeData(lineOutput[iband],GDT_Float64,y,iband);
325  }
326  catch(string errorstring){
327  cerr << errorstring << "in band " << iband << ", line " << y << endl;
328  }
329  }
330  progress=(1.0+y)/output.nrOfRow();
331  pfnProgress(progress,pszMessage,pProgressArg);
332  }
333 }
334 
335 void filter::Filter::stat(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method)
336 {
337  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
338  assert(output.nrOfCol()==input.nrOfCol());
339  vector<double> lineOutput(output.nrOfCol());
340  statfactory::StatFactory stat;
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  for(int x=0;x<input.nrOfCol();++x){
351  pixelInput=lineInput.selectCol(x);
352  switch(getFilterType(method)){
353  case(filter::median):
354  lineOutput[x]=stat.median(pixelInput);
355  break;
356  case(filter::min):
357  lineOutput[x]=stat.mymin(pixelInput);
358  break;
359  case(filter::max):
360  lineOutput[x]=stat.mymax(pixelInput);
361  break;
362  case(filter::sum):
363  lineOutput[x]=stat.sum(pixelInput);
364  break;
365  case(filter::var):
366  lineOutput[x]=stat.var(pixelInput);
367  break;
368  case(filter::stdev):
369  lineOutput[x]=sqrt(stat.var(pixelInput));
370  break;
371  case(filter::mean):
372  lineOutput[x]=stat.mean(pixelInput);
373  break;
374  default:
375  std::string errorString="method not supported";
376  throw(errorString);
377  break;
378  }
379  }
380  try{
381  output.writeData(lineOutput,GDT_Float64,y);
382  }
383  catch(string errorstring){
384  cerr << errorstring << "in line " << y << endl;
385  }
386  progress=(1.0+y)/output.nrOfRow();
387  pfnProgress(progress,pszMessage,pProgressArg);
388  }
389 }
390 
391 void filter::Filter::filter(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim)
392 {
393  Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
394  Vector2d<double> lineOutput(input.nrOfBand(),input.nrOfCol());;
395  const char* pszMessage;
396  void* pProgressArg=NULL;
397  GDALProgressFunc pfnProgress=GDALTermProgress;
398  double progress=0;
399  pfnProgress(progress,pszMessage,pProgressArg);
400  for(int y=0;y<input.nrOfRow();++y){
401  for(int iband=0;iband<input.nrOfBand();++iband)
402  input.readData(lineInput[iband],GDT_Float64,y,iband);
403  vector<double> pixelInput(input.nrOfBand());
404  vector<double> pixelOutput;
405  for(int x=0;x<input.nrOfCol();++x){
406  pixelInput=lineInput.selectCol(x);
407  filter(pixelInput,pixelOutput,method,dim);
408  for(int iband=0;iband<pixelOutput.size();++iband){
409  lineOutput[iband][x]=pixelOutput[iband];
410  // if(pixelInput[iband]!=0)
411  // assert(pixelOutput[iband]!=0);
412  }
413  }
414  for(int iband=0;iband<input.nrOfBand();++iband){
415  try{
416  output.writeData(lineOutput[iband],GDT_Float64,y,iband);
417  }
418  catch(string errorstring){
419  cerr << errorstring << "in band " << iband << ", line " << y << endl;
420  }
421  }
422  progress=(1.0+y)/output.nrOfRow();
423  pfnProgress(progress,pszMessage,pProgressArg);
424  }
425 }
426 
427 double filter::Filter::getCentreWavelength(const std::vector<double> &wavelengthIn, const Vector2d<double>& srf, const std::string& interpolationType, double delta, bool verbose)
428 {
429  assert(srf.size()==2);//[0]: wavelength, [1]: response function
430  int nband=srf[0].size();
431  double start=floor(wavelengthIn[0]);
432  double end=ceil(wavelengthIn.back());
433  if(verbose)
434  std::cout << "wavelengths in [" << start << "," << end << "]" << std::endl << std::flush;
435 
436  statfactory::StatFactory stat;
437 
438  gsl_interp_accel *acc;
439  stat.allocAcc(acc);
440  gsl_spline *spline;
441  stat.getSpline(interpolationType,nband,spline);
442  stat.initSpline(spline,&(srf[0][0]),&(srf[1][0]),nband);
443  if(verbose)
444  std::cout << "calculating norm of srf" << std::endl << std::flush;
445  double norm=0;
446  norm=gsl_spline_eval_integ(spline,srf[0].front(),srf[0].back(),acc);
447  if(verbose)
448  std::cout << "norm of srf: " << norm << std::endl << std::flush;
449  gsl_spline_free(spline);
450  gsl_interp_accel_free(acc);
451  std::vector<double> wavelength_fine;
452  for(double win=floor(wavelengthIn[0]);win<=ceil(wavelengthIn.back());win+=delta)
453  wavelength_fine.push_back(win);
454 
455  if(verbose)
456  std::cout << "interpolate wavelengths to " << wavelength_fine.size() << " entries " << std::endl;
457  std::vector<double> srf_fine;//spectral response function, interpolated for wavelength_fine
458 
459  stat.interpolateUp(srf[0],srf[1],wavelength_fine,interpolationType,srf_fine,verbose);
460  assert(srf_fine.size()==wavelength_fine.size());
461 
462  gsl_interp_accel *accOut;
463  stat.allocAcc(accOut);
464  gsl_spline *splineOut;
465  stat.getSpline(interpolationType,wavelength_fine.size(),splineOut);
466  assert(splineOut);
467 
468  std::vector<double> wavelengthOut(wavelength_fine.size());
469 
470  for(int iband=0;iband<wavelengthOut.size();++iband)
471  wavelengthOut[iband]=wavelength_fine[iband]*srf_fine[iband];
472 
473  stat.initSpline(splineOut,&(wavelength_fine[0]),&(wavelengthOut[0]),wavelength_fine.size());
474  double centreWavelength=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
475 
476  gsl_spline_free(splineOut);
477  gsl_interp_accel_free(accOut);
478 
479  return(centreWavelength);
480 }
481 
482 // 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){
483 // Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
484 // Vector2d<double> lineOutput(wavelengthOut.size(),input.nrOfCol());
485 // const char* pszMessage;
486 // void* pProgressArg=NULL;
487 // GDALProgressFunc pfnProgress=GDALTermProgress;
488 // double progress=0;
489 // pfnProgress(progress,pszMessage,pProgressArg);
490 // for(int y=0;y<input.nrOfRow();++y){
491 // for(int iband=0;iband<input.nrOfBand();++iband)
492 // input.readData(lineInput[iband],GDT_Float64,y,iband);
493 // applyFwhm<double>(wavelengthIn,lineInput,wavelengthOut,fwhm, interpolationType, lineOutput, verbose);
494 // for(int iband=0;iband<output.nrOfBand();++iband){
495 // try{
496 // output.writeData(lineOutput[iband],GDT_Float64,y,iband);
497 // }
498 // catch(string errorstring){
499 // cerr << errorstring << "in band " << iband << ", line " << y << endl;
500 // }
501 // }
502 // progress=(1.0+y)/output.nrOfRow();
503 // pfnProgress(progress,pszMessage,pProgressArg);
504 // }
505 // }
506 
507 // void filter::Filter::applySrf(const vector<double> &wavelengthIn, const ImgReaderGdal& input, const vector< Vector2d<double> > &srf, const std::string& interpolationType, ImgWriterGdal& output, bool verbose){
508 // assert(output.nrOfBand()==srf.size());
509 // double centreWavelength=0;
510 // Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
511 // const char* pszMessage;
512 // void* pProgressArg=NULL;
513 // GDALProgressFunc pfnProgress=GDALTermProgress;
514 // double progress=0;
515 // pfnProgress(progress,pszMessage,pProgressArg);
516 // for(int y=0;y<input.nrOfRow();++y){
517 // for(int iband=0;iband<input.nrOfBand();++iband)
518 // input.readData(lineInput[iband],GDT_Float64,y,iband);
519 // for(int isrf=0;isrf<srf.size();++isrf){
520 // vector<double> lineOutput(input.nrOfCol());
521 // centreWavelength=applySrf<double>(wavelengthIn,lineInput,srf[isrf], interpolationType, lineOutput, verbose);
522 // for(int iband=0;iband<output.nrOfBand();++iband){
523 // try{
524 // output.writeData(lineOutput,GDT_Float64,y,isrf);
525 // }
526 // catch(string errorstring){
527 // cerr << errorstring << "in band " << iband << ", line " << y << endl;
528 // }
529 // }
530 // }
531 // progress=(1.0+y)/output.nrOfRow();
532 // pfnProgress(progress,pszMessage,pProgressArg);
533 // }
534 // }
Generated on Wed Nov 19 2014 12:26:58 by   doxygen 1.8.6