pktools 2.6.7
Processing Kernel for geospatial data
Filter.h
1/**********************************************************************
2Filter.h: class for filtering
3Copyright (C) 2008-2012 Pieter Kempeneers
4
5This file is part of pktools
6
7pktools is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published by
9the Free Software Foundation, either version 3 of the License, or
10(at your option) any later version.
11
12pktools is distributed in the hope that it will be useful,
13but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15GNU General Public License for more details.
16
17You should have received a copy of the GNU General Public License
18along with pktools. If not, see <http://www.gnu.org/licenses/>.
19***********************************************************************/
20#ifndef _MYFILTER_H_
21#define _MYFILTER_H_
22
23#include <vector>
24#include <iostream>
25extern "C" {
26#include <gsl/gsl_sort.h>
27#include <gsl/gsl_wavelet.h>
28}
29#include "StatFactory.h"
30#include "imageclasses/ImgReaderGdal.h"
31#include "imageclasses/ImgWriterGdal.h"
32
33namespace filter
34{
35
36 enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6, dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13, sobelxy=14, sobelyx=-14, smooth=15, density=16, mode=17, mixed=18, smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24, stdev=25, dwt=26, dwti=27, dwt_cut=28, dwt_cut_from=29, savgolay=30, percentile=31, nvalid=32};
37
38 enum PADDING { symmetric=0, replicate=1, circular=2, zero=3};
39
40class Filter
41{
42public:
43 Filter(void);
44 Filter(const std::vector<double> &taps);
45 virtual ~Filter(){};
46
47 void setPadding(const std::string& padString){
48 m_padding=padString;
49 };
50
51 static const gsl_wavelet_type* getWaveletType(const std::string waveletType){
52 if(waveletType=="daubechies") return(gsl_wavelet_daubechies);
53 if(waveletType=="daubechies_centered") return(gsl_wavelet_daubechies_centered);
54 if(waveletType=="haar") return(gsl_wavelet_haar);
55 if(waveletType=="haar_centered") return(gsl_wavelet_haar_centered);
56 if(waveletType=="bspline") return(gsl_wavelet_bspline);
57 if(waveletType=="bspline_centered") return(gsl_wavelet_bspline_centered);
58 }
59 static FILTER_TYPE getFilterType(const std::string filterType){
60 std::map<std::string, FILTER_TYPE> m_filterMap;
61 initFilterMap(m_filterMap);
62 return m_filterMap[filterType];
63 };
64
65 void setTaps(const std::vector<double> &taps, bool normalize=true);
66 void pushClass(short theClass=1){m_class.push_back(theClass);};
67 void pushMask(short theMask=0){m_mask.push_back(theMask);};
68 unsigned int pushNoDataValue(double noDataValue);
69 unsigned int setNoDataValues(std::vector<double> vnodata);
70 void pushThreshold(double theThreshold){m_threshold.push_back(theThreshold);};
71 void setThresholds(const std::vector<double>& theThresholds){m_threshold=theThresholds;};
72 template<class T> void filter(const std::vector<T>& input, std::vector<T>& output);
73 template<class T> void filter(const std::vector<T>& input, std::vector<T>& output, const std::string& method, int dim);
74 template<class T> void smooth(const std::vector<T>& input, std::vector<T>& output, short dim);
75 template<class T> void smoothNoData(const std::vector<T>& input, const std::string& interpolationType, std::vector<T>& output);
76 template<class T> void filter(T* input, int inputSize, std::vector<T>& output);
77 template<class T> void smooth(T* input, int inputSize, std::vector<T>& output, short dim);
78 //template<class T> void morphology(const std::vector<T>& input, std::vector<T>& output, const std::string& method, int dim, bool verbose=false);
79 void morphology(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short verbose=0);
80 void filter(ImgReaderGdal& input, ImgWriterGdal& output);
81 void stat(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method);
82 void stats(ImgReaderGdal& input, ImgWriterGdal& output, const std::vector<std::string >& methods);
83 void filter(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim);
84 void getSavGolayCoefficients(std::vector<double> &c, int np, int nl, int nr, int ld, int m);
85 void ludcmp(std::vector<double> &a, std::vector<int> &indx, double &d);
86 void lubksb(std::vector<double> &a, std::vector<int> &indx, std::vector<double> &b);
87 /* void savgolay(const ImgReaderGdal& input, ImgWriterGdal& output, int np, int nl, int nr, int m); */
88 void smooth(ImgReaderGdal& input, ImgWriterGdal& output, short dim);
89 void smoothNoData(ImgReaderGdal& input, const std::string& interpolationType, ImgWriterGdal& output);
90 double getCentreWavelength(const std::vector<double> &wavelengthIn, const Vector2d<double>& srf, const std::string& interpolationType, double delta=1.0, bool verbose=false);
91 template<class T> double applySrf(const std::vector<double> &wavelengthIn, const std::vector<T>& input, const Vector2d<double>& srf, const std::string& interpolationType, T& output, double delta=1.0, bool normalize=false, bool verbose=false);
92 template<class T> double applySrf(const std::vector<double> &wavelengthIn, const Vector2d<T>& input, const Vector2d<double>& srf, const std::string& interpolationType, std::vector<T>& output, double delta=1.0, bool normalize=false, int down=1, bool transposeInput=false, bool verbose=false);
93
94 template<class T> void 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=false);
95 template<class T> void applyFwhm(const std::vector<double> &wavelengthIn, const Vector2d<T>& input, const std::vector<double> &wavelengthOut, const std::vector<double> &fwhm, const std::string& interpolationType, Vector2d<T>& output, int down=1, bool verbose=false);
96 void dwtForward(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family);
97 void dwtInverse(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family);
98 void dwtCut(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, double cut);
99 void dwtForward(std::vector<double>& data, const std::string& wavelet_type, int family);
100 void dwtInverse(std::vector<double>& data, const std::string& wavelet_type, int family);
101 void dwtCut(std::vector<double>& data, const std::string& wavelet_type, int family, double cut);
102 void dwtCutFrom(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, int band);
103
104private:
105
106 static void initFilterMap(std::map<std::string, FILTER_TYPE>& m_filterMap){
107 //initialize Map
108 m_filterMap["dwt"]=filter::dwt;
109 m_filterMap["dwti"]=filter::dwti;
110 m_filterMap["dwt_cut"]=filter::dwt_cut;
111 m_filterMap["dwt_cut_from"]=filter::dwt_cut_from;
112 m_filterMap["stdev"]=filter::stdev;
113 m_filterMap["var"]=filter::var;
114 m_filterMap["min"]=filter::min;
115 m_filterMap["max"]=filter::max;
116 m_filterMap["sum"]=filter::sum;
117 m_filterMap["mean"]=filter::mean;
118 m_filterMap["minmax"]=filter::minmax;
119 m_filterMap["dilate"]=filter::dilate;
120 m_filterMap["erode"]=filter::erode;
121 m_filterMap["close"]=filter::close;
122 m_filterMap["open"]=filter::open;
123 m_filterMap["homog"]=filter::homog;
124 m_filterMap["sobelx"]=filter::sobelx;
125 m_filterMap["sobely"]=filter::sobely;
126 m_filterMap["sobelxy"]=filter::sobelxy;
127 m_filterMap["sobelyx"]=filter::sobelyx;
128 m_filterMap["smooth"]=filter::smooth;
129 m_filterMap["density"]=filter::density;
130 m_filterMap["mode"]=filter::mode;
131 m_filterMap["mixed"]=filter::mixed;
132 m_filterMap["smoothnodata"]=filter::smoothnodata;
133 m_filterMap["threshold"]=filter::threshold;
134 m_filterMap["ismin"]=filter::ismin;
135 m_filterMap["ismax"]=filter::ismax;
136 m_filterMap["heterog"]=filter::heterog;
137 m_filterMap["order"]=filter::order;
138 m_filterMap["nvalid"]=filter::nvalid;
139 m_filterMap["median"]=filter::median;
140 m_filterMap["savgolay"]=filter::savgolay;
141 m_filterMap["percentile"]=filter::percentile;
142 }
143
144
145 static PADDING getPadding(const std::string& padString){
146 std::map<std::string, PADDING> padMap;
147 padMap["zero"]=filter::zero;
148 padMap["symmetric"]=filter::symmetric;
149 padMap["replicate"]=filter::replicate;
150 padMap["circular"]=filter::circular;
151 return(padMap[padString]);
152 };
153
154 std::vector<double> m_taps;
155 std::vector<short> m_class;
156 std::vector<short> m_mask;
157 std::string m_padding;
158 std::vector<double> m_noDataValues;
159 std::vector<double> m_threshold;
160};
161
162
163//input[band], output
164//returns wavelength for which srf is maximum
165 template<class T> double Filter::applySrf(const std::vector<double> &wavelengthIn, const std::vector<T>& input, const Vector2d<double>& srf, const std::string& interpolationType, T& output, double delta, bool normalize, bool verbose)
166{
167 assert(srf.size()==2);//[0]: wavelength, [1]: response function
168 int nband=srf[0].size();
169 double start=floor(wavelengthIn[0]);
170 double end=ceil(wavelengthIn.back());
171 if(verbose)
172 std::cout << "wavelengths in [" << start << "," << end << "]" << std::endl << std::flush;
173
175
176 gsl_interp_accel *acc;
177 stat.allocAcc(acc);
178 gsl_spline *spline;
179 stat.getSpline(interpolationType,nband,spline);
180 stat.initSpline(spline,&(srf[0][0]),&(srf[1][0]),nband);
181 if(verbose)
182 std::cout << "calculating norm of srf" << std::endl << std::flush;
183 double norm=0;
184 norm=gsl_spline_eval_integ(spline,srf[0].front(),srf[0].back(),acc);
185 if(verbose)
186 std::cout << "norm of srf: " << norm << std::endl << std::flush;
187 gsl_spline_free(spline);
188 gsl_interp_accel_free(acc);
189 //interpolate input and srf to delta
190
191 std::vector<double> wavelength_fine;
192 for(double win=floor(wavelengthIn[0]);win<=ceil(wavelengthIn.back());win+=delta)
193 wavelength_fine.push_back(win);
194
195 if(verbose)
196 std::cout << "interpolate wavelengths to " << wavelength_fine.size() << " entries " << std::endl;
197 std::vector<double> srf_fine;//spectral response function, interpolated for wavelength_fine
198
199 stat.interpolateUp(srf[0],srf[1],wavelength_fine,interpolationType,srf_fine,verbose);
200 assert(srf_fine.size()==wavelength_fine.size());
201
202 gsl_interp_accel *accOut;
203 stat.allocAcc(accOut);
204 gsl_spline *splineOut;
205 stat.getSpline(interpolationType,wavelength_fine.size(),splineOut);
206 assert(splineOut);
207
208 assert(wavelengthIn.size()==input.size());
209 std::vector<double> input_fine;
210 std::vector<double> product(wavelength_fine.size());
211 std::vector<double> wavelengthOut(wavelength_fine.size());
212 stat.interpolateUp(wavelengthIn,input,wavelength_fine,interpolationType,input_fine,verbose);
213
214 if(verbose)
215 std::cout << "input_fine.size(): " << input_fine.size() << std::endl;
216 for(int iband=0;iband<input_fine.size();++iband){
217 product[iband]=input_fine[iband]*srf_fine[iband];
218 wavelengthOut[iband]=wavelength_fine[iband]*srf_fine[iband];
219 }
220
221 assert(input_fine.size()==srf_fine.size());
222 assert(input_fine.size()==wavelength_fine.size());
223 stat.initSpline(splineOut,&(wavelength_fine[0]),&(product[0]),wavelength_fine.size());
224 if(normalize)
225 output=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
226 else
227 output=gsl_spline_eval_integ(splineOut,start,end,accOut);
228
229 stat.initSpline(splineOut,&(wavelength_fine[0]),&(wavelengthOut[0]),wavelength_fine.size());
230 double centreWavelength=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
231
232 gsl_spline_free(splineOut);
233 gsl_interp_accel_free(accOut);
234
235 return(centreWavelength);
236}
237
238//input[band][sample], output[sample] (if !transposeInput)
239//returns wavelength for which srf is maximum
240 template<class T> double Filter::applySrf(const std::vector<double> &wavelengthIn, const Vector2d<T>& input, const Vector2d<double>& srf, const std::string& interpolationType, std::vector<T>& output, double delta, bool normalize, int down, bool transposeInput, bool verbose)
241{
242 assert(srf.size()==2);//[0]: wavelength, [1]: response function
243 int nband=srf[0].size();
244 unsigned int nsample=(transposeInput)? input.size():input[0].size();
245 output.resize((nsample+down-1)/down);
246 double start=floor(wavelengthIn[0]);
247 double end=ceil(wavelengthIn.back());
248 if(verbose)
249 std::cout << "wavelengths in [" << start << "," << end << "]" << std::endl << std::flush;
250
252
253 gsl_interp_accel *acc;
254 stat.allocAcc(acc);
255 gsl_spline *spline;
256 stat.getSpline(interpolationType,nband,spline);
257 stat.initSpline(spline,&(srf[0][0]),&(srf[1][0]),nband);
258 if(verbose)
259 std::cout << "calculating norm of srf" << std::endl << std::flush;
260 double norm=0;
261 norm=gsl_spline_eval_integ(spline,srf[0].front(),srf[0].back(),acc);
262 if(verbose)
263 std::cout << "norm of srf: " << norm << std::endl << std::flush;
264 gsl_spline_free(spline);
265 gsl_interp_accel_free(acc);
266 //interpolate input and srf to delta
267
268 std::vector<double> wavelength_fine;
269 for(double win=floor(wavelengthIn[0]);win<=ceil(wavelengthIn.back());win+=delta)
270 wavelength_fine.push_back(win);
271
272 if(verbose)
273 std::cout << "interpolate wavelengths to " << wavelength_fine.size() << " entries " << std::endl;
274 std::vector<double> srf_fine;//spectral response function, interpolated for wavelength_fine
275
276 stat.interpolateUp(srf[0],srf[1],wavelength_fine,interpolationType,srf_fine,verbose);
277 assert(srf_fine.size()==wavelength_fine.size());
278
279 gsl_interp_accel *accOut;
280 stat.allocAcc(accOut);
281 gsl_spline *splineOut;
282 stat.getSpline(interpolationType,wavelength_fine.size(),splineOut);
283 assert(splineOut);
284
285 std::vector<double> wavelengthOut;
286 double centreWavelength=0;
287 for(int isample=0;isample<nsample;++isample){
288 if((isample+1+down/2)%down)
289 continue;
290 std::vector<T> inputValues;
291 if(transposeInput)
292 inputValues=input[isample];
293 else
294 input.selectCol(isample,inputValues);
295 assert(wavelengthIn.size()==inputValues.size());
296 std::vector<double> input_fine;
297 std::vector<double> product(wavelength_fine.size());
298 stat.interpolateUp(wavelengthIn,inputValues,wavelength_fine,interpolationType,input_fine,verbose);
299
300 for(int iband=0;iband<input_fine.size();++iband){
301 product[iband]=input_fine[iband]*srf_fine[iband];
302 if(wavelengthOut.size()<input_fine.size())
303 wavelengthOut.push_back(wavelength_fine[iband]*srf_fine[iband]);
304 }
305
306 assert(input_fine.size()==srf_fine.size());
307 assert(input_fine.size()==wavelength_fine.size());
308 stat.initSpline(splineOut,&(wavelength_fine[0]),&(product[0]),wavelength_fine.size());
309 if(normalize)
310 output[isample/down]=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
311 else
312 output[isample/down]=gsl_spline_eval_integ(splineOut,start,end,accOut);
313
314 stat.initSpline(splineOut,&(wavelength_fine[0]),&(wavelengthOut[0]),wavelength_fine.size());
315 if(centreWavelength>0);
316 else
317 centreWavelength=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
318 }
319 gsl_spline_free(splineOut);
320 gsl_interp_accel_free(accOut);
321
322 return(centreWavelength);
323}
324
325template<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){
326 double delta=1;//1 nm resolution
327 std::vector<double> stddev(fwhm.size());
328 for(int index=0;index<fwhm.size();++index)
329 stddev[index]=fwhm[index]/2.0/sqrt(2*log(2.0));//http://mathworld.wolfram.com/FullWidthatHalfMaximum.html
330 assert(wavelengthOut.size()==fwhm.size());
331 assert(wavelengthIn.size()==input.size());
332 assert(wavelengthIn[0]<=wavelengthOut[0]);
333 assert(wavelengthIn.back()>=wavelengthOut.back());
335 std::vector<double> input_fine;
336 std::vector<double> wavelength_fine;
337 for(double win=floor(wavelengthIn[0]);win<=ceil(wavelengthIn.back());win+=delta)
338 wavelength_fine.push_back(win);
339 if(verbose){
340 for(int index=0;index<wavelength_fine.size();++index)
341 std::cout << " " << wavelength_fine[index];
342 std::cout << std::endl;
343 std::cout << "interpolate input wavelength to " << delta << " nm resolution (size=" << wavelength_fine.size() << ")" << std::endl;
344 }
345 stat.interpolateUp(wavelengthIn,input,wavelength_fine,interpolationType,input_fine,verbose);
346 int nbandIn=wavelength_fine.size();
347
348 int nbandOut=wavelengthOut.size();
349 output.resize(nbandOut);
350 Vector2d<double> tf(nbandIn,nbandOut);
351 for(int indexOut=0;indexOut<nbandOut;++indexOut){
352 double norm=0;
353 for(int indexIn=0;indexIn<nbandIn;++indexIn){
354 // tf(indexIn,indexOut)=
355 tf[indexIn][indexOut]=
356 exp((wavelengthOut[indexOut]-wavelength_fine[indexIn])
357 *(wavelength_fine[indexIn]-wavelengthOut[indexOut])
358 /2.0/stddev[indexOut]
359 /stddev[indexOut]);
360 tf[indexIn][indexOut]/=sqrt(2.0*M_PI);
361 tf[indexIn][indexOut]/=stddev[indexOut];
362 norm+=tf[indexIn][indexOut];
363 }
364 output[indexOut]=0;
365 for(int indexIn=0;indexIn<nbandIn;++indexIn)
366 output[indexOut]+=input_fine[indexIn]*tf[indexIn][indexOut]/norm;
367 }
368}
369
370
371 //input[inBand][sample], output[outBand][sample]
372 template<class T> void Filter::applyFwhm(const std::vector<double> &wavelengthIn, const Vector2d<T>& input, const std::vector<double> &wavelengthOut, const std::vector<double> &fwhm, const std::string& interpolationType, Vector2d<T>& output, int down, bool verbose){
373 double delta=1;//1 nm resolution
374 std::vector<double> stddev(fwhm.size());
375 for(int index=0;index<fwhm.size();++index)
376 stddev[index]=fwhm[index]/2.0/sqrt(2*log(2.0));//http://mathworld.wolfram.com/FullWidthatHalfMaximum.html
378 std::vector<double> wavelength_fine;
379 for(double win=floor(wavelengthIn[0]);win<=ceil(wavelengthIn.back());win+=delta)
380 wavelength_fine.push_back(win);
381 assert(wavelengthOut.size()==fwhm.size());
382 assert(wavelengthIn[0]<=wavelengthOut[0]);
383 assert(wavelengthIn.back()>=wavelengthOut.back());
384 if(verbose){
385 for(int index=0;index<wavelength_fine.size();++index)
386 std::cout << " " << wavelength_fine[index];
387 std::cout << std::endl;
388 std::cout << "interpolate input wavelength to " << delta << " nm resolution (size=" << wavelength_fine.size() << ")" << std::endl;
389 }
390 int nbandIn=wavelength_fine.size();
391 int nbandOut=wavelengthOut.size();
392 output.resize(nbandOut,(input[0].size()+down-1)/down);
393
394 Vector2d<double> tf(nbandIn,nbandOut);
395 std::vector<double> norm(nbandOut);
396 for(int indexOut=0;indexOut<nbandOut;++indexOut){
397 norm[indexOut]=0;
398 for(int indexIn=0;indexIn<nbandIn;++indexIn){
399 tf[indexIn][indexOut]=
400 exp((wavelengthOut[indexOut]-wavelength_fine[indexIn])
401 *(wavelength_fine[indexIn]-wavelengthOut[indexOut])
402 /2.0/stddev[indexOut]
403 /stddev[indexOut]);
404 tf[indexIn][indexOut]/=sqrt(2.0*M_PI);
405 tf[indexIn][indexOut]/=stddev[indexOut];
406 norm[indexOut]+=tf[indexIn][indexOut];
407 }
408 }
409
410 for(int isample=0;isample<input[0].size();++isample){
411 if((isample+1+down/2)%down)
412 continue;
413 std::vector<T> inputValues;
414 input.selectCol(isample,inputValues);
415 assert(wavelengthIn.size()==inputValues.size());
416 for(int indexOut=0;indexOut<nbandOut;++indexOut){
417 std::vector<double> input_fine;
418 stat.interpolateUp(wavelengthIn,inputValues,wavelength_fine,interpolationType,input_fine,verbose);
419 output[indexOut][(isample+down-1)/down]=0;
420 for(int indexIn=0;indexIn<nbandIn;++indexIn){
421 output[indexOut][(isample+down-1)/down]+=input_fine[indexIn]*tf[indexIn][indexOut]/norm[indexOut];
422 }
423 }
424 }
425}
426
427 template<class T> void Filter::smooth(const std::vector<T>& input, std::vector<T>& output, short dim)
428{
429 assert(dim>0);
430 m_taps.resize(dim);
431 for(int itap=0;itap<dim;++itap)
432 m_taps[itap]=1.0/dim;
433 filter(input,output);
434 }
435
436 template<class T> void Filter::smoothNoData(const std::vector<T>& input, const std::string& interpolationType, std::vector<T>& output)
437{
439 stat.setNoDataValues(m_noDataValues);
440 std::vector<double> abscis(input.size());
441 for(int i=0;i<abscis.size();++i)
442 abscis[i]=i;
443 stat.interpolateNoData(abscis,input,interpolationType,output);
444 }
445
446template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T>& output)
447{
448 assert(input.size()>=m_taps.size());
449 output.resize(input.size());
450 int i=0;
451 //start: extend input by padding
452 for(i=0;i<m_taps.size()/2;++i){
453 //todo:introduce nodata?
454 output[i]=m_taps[m_taps.size()/2]*input[i];
455 for(int t=1;t<=m_taps.size()/2;++t){
456 output[i]+=m_taps[m_taps.size()/2+t]*input[i+t];
457 if(i>=t)
458 output[i]+=m_taps[m_taps.size()/2-t]*input[i-t];
459 else{
460 switch(getPadding(m_padding)){
461 case(replicate):
462 output[i]+=m_taps[m_taps.size()/2-t]*input[0];
463 break;
464 case(circular):
465 output[i]+=m_taps[m_taps.size()/2-t]*input[input.size()+i-t];
466 break;
467 case(zero):
468 output[i]+=m_taps[m_taps.size()/2-t]*0;
469 break;
470 case(symmetric):
471 default:
472 output[i]+=m_taps[m_taps.size()/2-t]*input[t-i];
473 break;
474 }
475 }
476 }
477 }
478 //main
479 for(i=m_taps.size()/2;i<input.size()-m_taps.size()/2;++i){
480 //todo:introduce nodata
481 T leaveOut=(*(m_taps.begin()))*input[i-m_taps.size()/2];
482 T include=(m_taps.back())*input[i+m_taps.size()/2];
483 output[i]=0;
484 for(int t=0;t<m_taps.size();++t)
485 output[i]+=input[i-m_taps.size()/2+t]*m_taps[t];
486 }
487 //end: extend input by padding
488 for(i=input.size()-m_taps.size()/2;i<input.size();++i){
489 //todo:introduce nodata?
490 output[i]=m_taps[m_taps.size()/2]*input[i];
491 //todo:introduce nodata?
492 for(int t=1;t<=m_taps.size()/2;++t){
493 output[i]+=m_taps[m_taps.size()/2-t]*input[i-t];
494 if(i+t<input.size())
495 output[i]+=m_taps[m_taps.size()/2+t]*input[i+t];
496 else{
497 switch(getPadding(m_padding)){
498 case(replicate):
499 output[i]+=m_taps[m_taps.size()/2+t]*input.back();
500 break;
501 case(circular):
502 output[i]+=m_taps[m_taps.size()/2+t]*input[t-1];
503 break;
504 case(zero):
505 output[i]+=m_taps[m_taps.size()/2+t]*0;
506 break;
507 case(symmetric):
508 default:
509 output[i]+=m_taps[m_taps.size()/2+t]*input[i-t];
510 break;
511 }
512 }
513 //output[i]+=(m_taps[m_taps.size()/2+t]+m_taps[m_taps.size()/2-t])*input[i-t];
514 }
515 }
516}
517
518//todo: filling statBuffer can be optimized (no need to clear and fill entire buffer, just push back new value...)
519 template<class T> void Filter::filter(const std::vector<T>& input, std::vector<T>& output, const std::string& method, int dim)
520{
521 bool verbose=false;
522 assert(dim);
523 output.resize(input.size());
524 int i=0;
526 stat.setNoDataValues(m_noDataValues);
527 std::vector<T> statBuffer;
528 short binValue=0;
529 //start: extend input by padding
530 for(i=0;i<dim/2;++i){
531 binValue=0;
532 for(int iclass=0;iclass<m_class.size();++iclass){
533 if(input[i]==m_class[iclass]){
534 binValue=m_class[0];
535 break;
536 }
537 }
538 if(m_class.size())
539 statBuffer.push_back(binValue);
540 else
541 statBuffer.push_back(input[i]);
542
543 for(int t=1;t<=dim/2;++t){
544 T theValue=input[i+t];
545 for(int iclass=0;iclass<m_class.size();++iclass){
546 if(theValue==m_class[iclass]){
547 binValue=m_class[0];
548 break;
549 }
550 }
551 if(m_class.size())
552 statBuffer.push_back(binValue);
553 else
554 statBuffer.push_back(theValue);
555
556 if(i>=t){
557 theValue=input[i-t];
558 }
559 else{
560 switch(getPadding(m_padding)){
561 case(replicate):
562 theValue=input[0];
563 break;
564 case(circular):
565 theValue=input[input.size()+i-t];
566 break;
567 case(zero):
568 theValue=0;
569 break;
570 case(symmetric):
571 default:
572 theValue=input[t-i];
573 break;
574 }
575 }
576 for(int iclass=0;iclass<m_class.size();++iclass){
577 if(theValue==m_class[iclass]){
578 binValue=m_class[0];
579 break;
580 }
581 }
582 if(m_class.size())
583 statBuffer.push_back(binValue);
584 else
585 statBuffer.push_back(theValue);
586 }
587
588 switch(getFilterType(method)){
589 case(filter::nvalid):
590 output[i]=stat.nvalid(statBuffer);
591 break;
592 case(filter::median):
593 output[i]=stat.median(statBuffer);
594 break;
595 case(filter::min):
596 case(filter::erode):
597 output[i]=stat.mymin(statBuffer);
598 break;
599 case(filter::max):
600 case(filter::dilate):
601 output[i]=stat.mymax(statBuffer);
602 break;
603 case(filter::sum):
604 output[i]=sqrt(stat.sum(statBuffer));
605 break;
606 case(filter::var):
607 output[i]=stat.var(statBuffer);
608 break;
609 case(filter::stdev):
610 output[i]=sqrt(stat.var(statBuffer));
611 break;
612 case(filter::mean):
613 output[i]=stat.mean(statBuffer);
614 break;
615 case(filter::percentile):
616 assert(m_threshold.size());
617 output[i]=stat.percentile(statBuffer,statBuffer.begin(),statBuffer.end(),m_threshold[0]);
618 break;
619 default:{
620 std::ostringstream ess;
621 ess << "method " << method << " (" << getFilterType(method) << ") not supported";
622 throw(ess.str());
623 break;
624 }
625 }
626 }
627 //main
628 statBuffer.clear();
629 for(i=dim/2;i<input.size()-dim/2;++i){
630 binValue=0;
631 for(int t=0;t<dim;++t){
632 for(int iclass=0;iclass<m_class.size();++iclass){
633 if(input[i-dim/2+t]==m_class[iclass]){
634 binValue=m_class[0];
635 break;
636 }
637 }
638 if(m_class.size())
639 statBuffer.push_back(binValue);
640 else
641 statBuffer.push_back(input[i-dim/2+t]);
642 }
643 switch(getFilterType(method)){
644 case(filter::nvalid):
645 output[i]=stat.nvalid(statBuffer);
646 break;
647 case(filter::median):
648 output[i]=stat.median(statBuffer);
649 break;
650 case(filter::min):
651 case(filter::erode):
652 output[i]=stat.mymin(statBuffer);
653 break;
654 case(filter::max):
655 case(filter::dilate):
656 output[i]=stat.mymax(statBuffer);
657 break;
658 case(filter::sum):
659 output[i]=sqrt(stat.sum(statBuffer));
660 break;
661 case(filter::var):
662 output[i]=stat.var(statBuffer);
663 break;
664 case(filter::mean):
665 output[i]=stat.mean(statBuffer);
666 break;
667 case(filter::percentile):
668 assert(m_threshold.size());
669 output[i]=stat.percentile(statBuffer,statBuffer.begin(),statBuffer.end(),m_threshold[0]);
670 break;
671 default:
672 std::string errorString="method not supported";
673 throw(errorString);
674 break;
675 }
676 statBuffer.clear();
677 }
678 //end: extend input by padding
679 for(i=input.size()-dim/2;i<input.size();++i){
680 binValue=0;
681 for(int iclass=0;iclass<m_class.size();++iclass){
682 if(input[i]==m_class[iclass]){
683 binValue=m_class[0];
684 break;
685 }
686 }
687 if(m_class.size())
688 statBuffer.push_back(binValue);
689 else
690 statBuffer.push_back(input[i]);
691
692 for(int t=1;t<=dim/2;++t){
693 T theValue=input[i-t];
694 for(int iclass=0;iclass<m_class.size();++iclass){
695 if(theValue==m_class[iclass]){
696 binValue=m_class[0];
697 break;
698 }
699 }
700 if(m_class.size())
701 statBuffer.push_back(binValue);
702 else
703 statBuffer.push_back(theValue);
704 if(i+t<input.size())
705 theValue=input[i+t];
706 else{
707 switch(getPadding(m_padding)){
708 case(replicate):
709 theValue=input.back();
710 break;
711 case(circular):
712 theValue=input[t-1];
713 break;
714 case(zero):
715 theValue=0;
716 break;
717 case(symmetric):
718 default:
719 theValue=input[i-t];
720 break;
721 }
722 }
723 for(int iclass=0;iclass<m_class.size();++iclass){
724 if(theValue==m_class[iclass]){
725 binValue=m_class[0];
726 break;
727 }
728 }
729 if(m_class.size())
730 statBuffer.push_back(binValue);
731 else
732 statBuffer.push_back(theValue);
733 }
734 switch(getFilterType(method)){
735 case(filter::nvalid):
736 output[i]=stat.nvalid(statBuffer);
737 break;
738 case(filter::median):
739 output[i]=stat.median(statBuffer);
740 break;
741 case(filter::min):
742 case(filter::erode):
743 output[i]=stat.mymin(statBuffer);
744 break;
745 case(filter::max):
746 case(filter::dilate):
747 output[i]=stat.mymax(statBuffer);
748 break;
749 case(filter::sum):
750 output[i]=sqrt(stat.sum(statBuffer));
751 break;
752 case(filter::var):
753 output[i]=stat.var(statBuffer);
754 break;
755 case(filter::mean):
756 output[i]=stat.mean(statBuffer);
757 break;
758 case(filter::percentile):
759 assert(m_threshold.size());
760 output[i]=stat.percentile(statBuffer,statBuffer.begin(),statBuffer.end(),m_threshold[0]);
761 break;
762 default:
763 std::string errorString="method not supported";
764 throw(errorString);
765 break;
766 }
767 }
768 }
769
770 template<class T> void Filter::smooth(T* input, int inputSize, std::vector<T>& output, short dim)
771{
772 assert(dim>0);
773 m_taps.resize(dim);
774 for(int itap=0;itap<dim;++itap)
775 m_taps[itap]=1.0/dim;
776 filter(input,output);
777 }
778
779template<class T> void Filter::filter(T* input, int inputSize, std::vector<T>& output)
780{
781 assert(inputSize>=m_taps.size());
782 output.resize(inputSize);
783 int i=0;
784
785 //start: extend input by padding
786 for(i=0;i<m_taps.size()/2;++i){
787 //todo:introduce nodata
788 output[i]=m_taps[m_taps.size()/2]*input[i];
789
790 for(int t=1;t<=m_taps.size()/2;++t){
791 output[i]+=m_taps[m_taps.size()/2+t]*input[i+t];
792 if(i>=t)
793 output[i]+=m_taps[m_taps.size()/2-t]*input[i-t];
794 else{
795 switch(getPadding(m_padding)){
796 case(replicate):
797 output[i]+=m_taps[m_taps.size()/2-t]*input[0];
798 break;
799 case(circular):
800 output[i]+=m_taps[m_taps.size()/2-t]*input[input.size()+i-t];
801 break;
802 case(zero):
803 output[i]+=m_taps[m_taps.size()/2-t]*0;
804 break;
805 case(symmetric):
806 default:
807 output[i]+=m_taps[m_taps.size()/2-t]*input[t-i];
808 break;
809 }
810 }
811 }
812 }
813 //main
814 for(i=m_taps.size()/2;i<input.size()-m_taps.size()/2;++i){
815 //todo:introduce nodata
816 T leaveOut=(*(m_taps.begin()))*input[i-m_taps.size()/2];
817 T include=(m_taps.back())*input[i+m_taps.size()/2];
818 output[i]=0;
819 for(int t=0;t<m_taps.size();++t)
820 output[i]+=input[i-m_taps.size()/2+t]*m_taps[t];
821 }
822 //end: extend input by padding
823 for(i=input.size()-m_taps.size()/2;i<input.size();++i){
824 //todo:introduce nodata
825 output[i]=m_taps[m_taps.size()/2]*input[i];
826 //todo:introduce nodata
827 for(int t=1;t<=m_taps.size()/2;++t){
828 output[i]+=m_taps[m_taps.size()/2-t]*input[i-t];
829 if(i+t<input.size())
830 output[i]+=m_taps[m_taps.size()/2+t]*input[i+t];
831 else{
832 switch(getPadding(m_padding)){
833 case(replicate):
834 output[i]+=m_taps[m_taps.size()/2+t]*input.back();
835 break;
836 case(circular):
837 output[i]+=m_taps[m_taps.size()/2+t]*input[t-1];
838 break;
839 case(zero):
840 output[i]+=m_taps[m_taps.size()/2+t]*0;
841 break;
842 case(symmetric):
843 default:
844 output[i]+=m_taps[m_taps.size()/2+t]*input[i-t];
845 break;
846 }
847 }
848 }
849 }
850}
851
852}
853
854#endif /* _MYFILTER_H_ */