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 ( )
inlinevirtual

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());
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}

◆ 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
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
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
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}
int nrOfRow(void) const
Get the number of rows of this dataset.
int nrOfBand(void) const
Get the number of bands of this dataset.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
void readData(T &value, int col, int row, int band=0)
Read a single pixel cell value at a specific column and row for a specific band (all indices start co...
Definition: ImgReaderGdal.h:95
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;
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
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)
inlinestatic

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)
inlinestatic

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{
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());
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());
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: