| pktools 2.6.7 Processing Kernel for geospatial data

Public Member Functions | Static Public Member Functions | List of all members
filter::Filter Class Reference

## Public Member Functions
| Filter (const std::vector< double > &taps)
void | setPadding (const std::string &padString)
void | setTaps (const std::vector< double > &taps, bool normalize=true)
void | pushClass (short theClass=1)
void | pushMask (short theMask=0)
unsigned int | pushNoDataValue (double noDataValue)
unsigned int | setNoDataValues (std::vector< double > vnodata)
void | pushThreshold (double theThreshold)
void | setThresholds (const std::vector< double > &theThresholds)
template<class T >
void | filter (const std::vector< T > &input, std::vector< T > &output)
template<class T >
void | filter (const std::vector< T > &input, std::vector< T > &output, const std::string &method, int dim)
template<class T >
void | smooth (const std::vector< T > &input, std::vector< T > &output, short dim)
template<class T >
void | smoothNoData (const std::vector< T > &input, const std::string &interpolationType, std::vector< T > &output)
template<class T >
void | filter (T *input, int inputSize, std::vector< T > &output)
template<class T >
void | smooth (T *input, int inputSize, std::vector< T > &output, short dim)
void | morphology (ImgReaderGdal &input, ImgWriterGdal &output, const std::string &method, int dim, short verbose=0)
void | filter (ImgReaderGdal &input, ImgWriterGdal &output)
void | stat (ImgReaderGdal &input, ImgWriterGdal &output, const std::string &method)
void | stats (ImgReaderGdal &input, ImgWriterGdal &output, const std::vector< std::string > &methods)
void | filter (ImgReaderGdal &input, ImgWriterGdal &output, const std::string &method, int dim)
void | getSavGolayCoefficients (std::vector< double > &c, int np, int nl, int nr, int ld, int m)
void | ludcmp (std::vector< double > &a, std::vector< int > &indx, double &d)
void | lubksb (std::vector< double > &a, std::vector< int > &indx, std::vector< double > &b)
void | smooth (ImgReaderGdal &input, ImgWriterGdal &output, short dim)
void | smoothNoData (ImgReaderGdal &input, const std::string &interpolationType, ImgWriterGdal &output)
double | getCentreWavelength (const std::vector< double > &wavelengthIn, const Vector2d< double > &srf, const std::string &interpolationType, double delta=1.0, bool verbose=false)
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)
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)
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)
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)
void | dwtForward (ImgReaderGdal &input, ImgWriterGdal &output, const std::string &wavelet_type, int family)
void | dwtInverse (ImgReaderGdal &input, ImgWriterGdal &output, const std::string &wavelet_type, int family)
void | dwtCut (ImgReaderGdal &input, ImgWriterGdal &output, const std::string &wavelet_type, int family, double cut)
void | dwtForward (std::vector< double > &data, const std::string &wavelet_type, int family)
void | dwtInverse (std::vector< double > &data, const std::string &wavelet_type, int family)
void | dwtCut (std::vector< double > &data, const std::string &wavelet_type, int family, double cut)
void | dwtCutFrom (ImgReaderGdal &input, ImgWriterGdal &output, const std::string &wavelet_type, int family, int band)

## Static Public Member Functions
static const gsl_wavelet_type * | getWaveletType (const std::string waveletType)
static FILTER_TYPE | getFilterType (const std::string filterType)

## Detailed Description

Definition at line 40 of file Filter.h.

## Constructor & Destructor Documentation

## ◆ Filter() [1/2]

filter::Filter::Filter | ( | void |  | ) |

Definition at line 27 of file Filter.cc.

28 : m_padding("symmetric")
29{
30}

## ◆ Filter() [2/2]

filter::Filter::Filter | ( | const std::vector< double > & | taps | ) |

Definition at line 33 of file Filter.cc.

34 : m_padding("symmetric")
35{
36 setTaps(taps);
37}

## ◆ ~Filter()

virtual filter::Filter::~Filter | ( |  | ) |
| inline virtual

Definition at line 45 of file Filter.h.

45{};

## Member Function Documentation

## ◆ applyFwhm() [1/2]

template<class T >

void filter::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 = `false`
| ) |  |

Definition at line 325 of file Filter.h.

325 {
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());
334 statfactory::StatFactory stat;
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}
Vector2d< double >
statfactory::StatFactory
Definition: StatFactory.h:43

## ◆ applyFwhm() [2/2]

template<class T >

void filter::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 = `1`,
|  | bool | verbose = `false`
| ) |  |

Definition at line 372 of file Filter.h.

372 {
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
377 statfactory::StatFactory stat;
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}

## ◆ applySrf() [1/2]

template<class T >

double filter::Filter::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`
| ) |  |

Definition at line 165 of file Filter.h.

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
174 statfactory::StatFactory stat;
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}

## ◆ applySrf() [2/2]

template<class T >

double filter::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 = `1.0`,
|  | bool | normalize = `false`,
|  | int | down = `1`,
|  | bool | transposeInput = `false`,
|  | bool | verbose = `false`
| ) |  |

Definition at line 240 of file Filter.h.

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
251 statfactory::StatFactory stat;
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}

## ◆ dwtCut() [1/2]

void filter::Filter::dwtCut | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | const std::string & | wavelet_type,
|  | int | family,
|  | double | cut
| ) |  |

Definition at line 127 of file Filter.cc.

127 {
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}
ImgRasterGdal::nrOfRow
int nrOfRow(void) const
Get the number of rows of this dataset.
Definition: ImgRasterGdal.h:100
ImgRasterGdal::nrOfBand
int nrOfBand(void) const
Get the number of bands of this dataset.
Definition: ImgRasterGdal.h:102
ImgRasterGdal::nrOfCol
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
ImgReaderGdal::readData
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
ImgWriterGdal::writeData
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

## ◆ dwtCut() [2/2]

void filter::Filter::dwtCut | ( | std::vector< double > & | data,
|  | const std::string & | wavelet_type,
|  | int | family,
|  | double | cut
| ) |  |

Definition at line 232 of file Filter.cc.

232 {
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}

## ◆ dwtCutFrom()

void filter::Filter::dwtCutFrom | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | const std::string & | wavelet_type,
|  | int | family,
|  | int | band
| ) |  |

Definition at line 158 of file Filter.cc.

158 {
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}

## ◆ dwtForward() [1/2]

void filter::Filter::dwtForward | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | const std::string & | wavelet_type,
|  | int | family
| ) |  |

Definition at line 65 of file Filter.cc.

65 {
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}

## ◆ dwtForward() [2/2]

void filter::Filter::dwtForward | ( | std::vector< double > & | data,
|  | const std::string & | wavelet_type,
|  | int | family
| ) |  |

Definition at line 195 of file Filter.cc.

195 {
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}

## ◆ dwtInverse() [1/2]

void filter::Filter::dwtInverse | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | const std::string & | wavelet_type,
|  | int | family
| ) |  |

Definition at line 96 of file Filter.cc.

96 {
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}

## ◆ dwtInverse() [2/2]

void filter::Filter::dwtInverse | ( | std::vector< double > & | data,
|  | const std::string & | wavelet_type,
|  | int | family
| ) |  |

Definition at line 214 of file Filter.cc.

214 {
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}

## ◆ filter() [1/5]

template<class T >

void filter::Filter::filter | ( | const std::vector< T > & | input,
|  | std::vector< T > & | output
| ) |  |

Definition at line 446 of file Filter.h.

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}

## ◆ filter() [2/5]

template<class T >

void filter::Filter::filter | ( | const std::vector< T > & | input,
|  | std::vector< T > & | output,
|  | const std::string & | method,
|  | int | dim
| ) |  |

Definition at line 519 of file Filter.h.

520{
521 bool verbose=false;
522 assert(dim);
523 output.resize(input.size());
524 int i=0;
525 statfactory::StatFactory stat;
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 }

## ◆ filter() [3/5]

void filter::Filter::filter | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output
| ) |  |

Definition at line 337 of file Filter.cc.

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}

## ◆ filter() [4/5]

void filter::Filter::filter | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | const std::string & | method,
|  | int | dim
| ) |  |

Definition at line 509 of file Filter.cc.

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}

## ◆ filter() [5/5]

template<class T >

void filter::Filter::filter | ( | T * | input,
|  | int | inputSize,
|  | std::vector< T > & | output
| ) |  |

Definition at line 779 of file Filter.h.

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}

## ◆ getCentreWavelength()

double filter::Filter::getCentreWavelength | ( | const std::vector< double > & | wavelengthIn,
|  | const Vector2d< double > & | srf,
|  | const std::string & | interpolationType,
|  | double | delta = `1.0`,
|  | bool | verbose = `false`
| ) |  |

Definition at line 684 of file Filter.cc.

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
693 statfactory::StatFactory stat;
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}

## ◆ getFilterType()

static FILTER_TYPE filter::Filter::getFilterType | ( | const std::string | filterType | ) |
| inline static

Definition at line 59 of file Filter.h.

59 {
60 std::map<std::string, FILTER_TYPE> m_filterMap;
61 initFilterMap(m_filterMap);
62 return m_filterMap[filterType];
63 };

## ◆ getSavGolayCoefficients()

void filter::Filter::getSavGolayCoefficients | ( | std::vector< double > & | c,
|  | int | np,
|  | int | nl,
|  | int | nr,
|  | int | ld,
|  | int | m
| ) |  |

Definition at line 545 of file Filter.cc.

545 {
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}

## ◆ getWaveletType()

static const gsl_wavelet_type * filter::Filter::getWaveletType | ( | const std::string | waveletType | ) |
| inline static

Definition at line 51 of file Filter.h.

51 {
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 }

## ◆ lubksb()

void filter::Filter::lubksb | ( | std::vector< double > & | a,
|  | std::vector< int > & | indx,
|  | std::vector< double > & | b
| ) |  |

Definition at line 660 of file Filter.cc.

660 {
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}

## ◆ ludcmp()

void filter::Filter::ludcmp | ( | std::vector< double > & | a,
|  | std::vector< int > & | indx,
|  | double & | d
| ) |  |

Definition at line 601 of file Filter.cc.

601 {
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}

## ◆ morphology()

void filter::Filter::morphology | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | const std::string & | method,
|  | int | dim,
|  | short | verbose = `0`
| ) |  |

Definition at line 260 of file Filter.cc.

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}

## ◆ pushClass()

void filter::Filter::pushClass | ( | short | theClass = `1` | ) |
| inline

Definition at line 66 of file Filter.h.

66{m_class.push_back(theClass);};

## ◆ pushMask()

void filter::Filter::pushMask | ( | short | theMask = `0` | ) |
| inline

Definition at line 67 of file Filter.h.

67{m_mask.push_back(theMask);};

## ◆ pushNoDataValue()

unsigned int filter::Filter::pushNoDataValue | ( | double | noDataValue | ) |

Definition at line 54 of file Filter.cc.

54 {
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};

## ◆ pushThreshold()

void filter::Filter::pushThreshold | ( | double | theThreshold | ) |
| inline

Definition at line 70 of file Filter.h.

70{m_threshold.push_back(theThreshold);};

## ◆ setNoDataValues()

unsigned int filter::Filter::setNoDataValues | ( | std::vector< double > | vnodata | ) |

Definition at line 60 of file Filter.cc.

60 {
61 m_noDataValues=vnodata;
62 return m_noDataValues.size();
63};

## ◆ setPadding()

void filter::Filter::setPadding | ( | const std::string & | padString | ) |
| inline

Definition at line 47 of file Filter.h.

47 {
48 m_padding=padString;
49 };

## ◆ setTaps()

void filter::Filter::setTaps | ( | const std::vector< double > & | taps,
|  | bool | normalize = `true`
| ) |  |

Definition at line 39 of file Filter.cc.

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}

## ◆ setThresholds()

void filter::Filter::setThresholds | ( | const std::vector< double > & | theThresholds | ) |
| inline

Definition at line 71 of file Filter.h.

71{m_threshold=theThresholds;};

## ◆ smooth() [1/3]

template<class T >

void filter::Filter::smooth | ( | const std::vector< T > & | input,
|  | std::vector< T > & | output,
|  | short | dim
| ) |  |

Definition at line 427 of file Filter.h.

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 }

## ◆ smooth() [2/3]

void filter::Filter::smooth | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | short | dim
| ) |  |

Definition at line 328 of file Filter.cc.

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}

## ◆ smooth() [3/3]

template<class T >

void filter::Filter::smooth | ( | T * | input,
|  | int | inputSize,
|  | std::vector< T > & | output,
|  | short | dim
| ) |  |

Definition at line 770 of file Filter.h.

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 }

## ◆ smoothNoData() [1/2]

template<class T >

void filter::Filter::smoothNoData | ( | const std::vector< T > & | input,
|  | const std::string & | interpolationType,
|  | std::vector< T > & | output
| ) |  |

Definition at line 436 of file Filter.h.

437{
438 statfactory::StatFactory stat;
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 }

## ◆ smoothNoData() [2/2]

void filter::Filter::smoothNoData | ( | ImgReaderGdal & | input,
|  | const std::string & | interpolationType,
|  | ImgWriterGdal & | output
| ) |  |

Definition at line 295 of file Filter.cc.

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}

## ◆ stat()

void filter::Filter::stat | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | const std::string & | method
| ) |  |

Definition at line 370 of file Filter.cc.

371{
372 Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
373 assert(output.nrOfCol()==input.nrOfCol());
374 vector<double> lineOutput(output.nrOfCol());
375 statfactory::StatFactory stat;
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}

## ◆ stats()

void filter::Filter::stats | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | const std::vector< std::string > & | methods
| ) |  |

Definition at line 431 of file Filter.cc.

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());
437 statfactory::StatFactory stat;
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}

The documentation for this class was generated from the following files:

- /tmp/slackdce/slackrepo.NvCr9x/build_pktools/PKTOOLS-2.6.7.6/src/algorithms/Filter.h
- /tmp/slackdce/slackrepo.NvCr9x/build_pktools/PKTOOLS-2.6.7.6/src/algorithms/Filter.cc