27filter::Filter::Filter(
void)
28 : m_padding(
"symmetric")
33filter::Filter::Filter(
const vector<double> &taps)
34 : m_padding(
"symmetric")
39void filter::Filter::setTaps(
const vector<double> &taps,
bool normalize)
41 m_taps.resize(taps.size());
43 for(
int itap=0;itap<taps.size();++itap)
46 for(
int itap=0;itap<taps.size();++itap)
47 m_taps[itap]=taps[itap]/norm;
51 assert(m_taps.size()%2);
54unsigned int filter::Filter::pushNoDataValue(
double noDataValue){
55 if(find(m_noDataValues.begin(),m_noDataValues.end(),noDataValue)==m_noDataValues.end())
56 m_noDataValues.push_back(noDataValue);
57 return m_noDataValues.size();
60unsigned int filter::Filter::setNoDataValues(std::vector<double> vnodata){
61 m_noDataValues=vnodata;
62 return m_noDataValues.size();
66 const char* pszMessage;
67 void* pProgressArg=NULL;
68 GDALProgressFunc pfnProgress=GDALTermProgress;
70 pfnProgress(progress,pszMessage,pProgressArg);
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];
83 for(
int iband=0;iband<input.
nrOfBand();++iband){
85 output.
writeData(lineOutput[iband],y,iband);
87 catch(
string errorstring){
88 cerr << errorstring <<
"in band " << iband <<
", line " << y << endl;
91 progress=(1.0+y)/output.
nrOfRow();
92 pfnProgress(progress,pszMessage,pProgressArg);
97 const char* pszMessage;
98 void* pProgressArg=NULL;
99 GDALProgressFunc pfnProgress=GDALTermProgress;
101 pfnProgress(progress,pszMessage,pProgressArg);
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];
114 for(
int iband=0;iband<input.
nrOfBand();++iband){
116 output.
writeData(lineOutput[iband],y,iband);
118 catch(
string errorstring){
119 cerr << errorstring <<
"in band " << iband <<
", line " << y << endl;
122 progress=(1.0+y)/output.
nrOfRow();
123 pfnProgress(progress,pszMessage,pProgressArg);
127void filter::Filter::dwtCut(
ImgReaderGdal& input,
ImgWriterGdal& output,
const std::string& wavelet_type,
int family,
double cut){
128 const char* pszMessage;
129 void* pProgressArg=NULL;
130 GDALProgressFunc pfnProgress=GDALTermProgress;
132 pfnProgress(progress,pszMessage,pProgressArg);
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];
145 for(
int iband=0;iband<input.
nrOfBand();++iband){
147 output.
writeData(lineOutput[iband],y,iband);
149 catch(
string errorstring){
150 cerr << errorstring <<
"in band " << iband <<
", line " << y << endl;
153 progress=(1.0+y)/output.
nrOfRow();
154 pfnProgress(progress,pszMessage,pProgressArg);
158void filter::Filter::dwtCutFrom(
ImgReaderGdal& input,
ImgWriterGdal& output,
const std::string& wavelet_type,
int family,
int band){
159 const char* pszMessage;
160 void* pProgressArg=NULL;
161 GDALProgressFunc pfnProgress=GDALTermProgress;
163 pfnProgress(progress,pszMessage,pProgressArg);
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){
177 dwtInverse(pixelInput,wavelet_type,family);
178 for(
int iband=0;iband<input.
nrOfBand();++iband)
179 lineOutput[iband][x]=pixelInput[iband];
181 for(
int iband=0;iband<input.
nrOfBand();++iband){
183 output.
writeData(lineOutput[iband],y,iband);
185 catch(
string errorstring){
186 cerr << errorstring <<
"in band " << iband <<
", line " << y << endl;
189 progress=(1.0+y)/output.
nrOfRow();
190 pfnProgress(progress,pszMessage,pProgressArg);
195void filter::Filter::dwtForward(std::vector<double>& data,
const std::string& wavelet_type,
int family){
196 int origsize=data.size();
198 while(data.size()&(data.size()-1))
199 data.push_back(data.back());
201 int nsize=data.size();
203 gsl_wavelet_workspace *work;
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);
214void filter::Filter::dwtInverse(std::vector<double>& data,
const std::string& wavelet_type,
int family){
215 int origsize=data.size();
217 while(data.size()&(data.size()-1))
218 data.push_back(data.back());
219 int nsize=data.size();
221 gsl_wavelet_workspace *work;
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);
232void filter::Filter::dwtCut(std::vector<double>& data,
const std::string& wavelet_type,
int family,
double cut){
233 int origsize=data.size();
235 while(data.size()&(data.size()-1))
236 data.push_back(data.back());
237 int nsize=data.size();
239 gsl_wavelet_workspace *work;
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]);
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++)
253 gsl_wavelet_transform_inverse(w,&(data[0]),1,nsize,work);
254 data.erase(data.begin()+origsize,data.end());
256 gsl_wavelet_free (w);
257 gsl_wavelet_workspace_free (work);
260void filter::Filter::morphology(
ImgReaderGdal& input,
ImgWriterGdal& output,
const std::string& method,
int dim,
short verbose)
265 const char* pszMessage;
266 void* pProgressArg=NULL;
267 GDALProgressFunc pfnProgress=GDALTermProgress;
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);
279 for(
int iband=0;iband<input.
nrOfBand();++iband)
280 lineOutput[iband][x]=pixelOutput[iband];
282 for(
int iband=0;iband<input.
nrOfBand();++iband){
284 output.
writeData(lineOutput[iband],y,iband);
286 catch(
string errorstring){
287 cerr << errorstring <<
"in band " << iband <<
", line " << y << endl;
290 progress=(1.0+y)/output.
nrOfRow();
291 pfnProgress(progress,pszMessage,pProgressArg);
299 const char* pszMessage;
300 void* pProgressArg=NULL;
301 GDALProgressFunc pfnProgress=GDALTermProgress;
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];
315 for(
int iband=0;iband<input.
nrOfBand();++iband){
317 output.
writeData(lineOutput[iband],y,iband);
319 catch(
string errorstring){
320 cerr << errorstring <<
"in band " << iband <<
", line " << y << endl;
323 progress=(1.0+y)/output.
nrOfRow();
324 pfnProgress(progress,pszMessage,pProgressArg);
332 for(
int itap=0;itap<dim;++itap)
333 m_taps[itap]=1.0/dim;
334 filter(input,output);
341 const char* pszMessage;
342 void* pProgressArg=NULL;
343 GDALProgressFunc pfnProgress=GDALTermProgress;
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];
357 for(
int iband=0;iband<input.
nrOfBand();++iband){
359 output.
writeData(lineOutput[iband],y,iband);
361 catch(
string errorstring){
362 cerr << errorstring <<
"in band " << iband <<
", line " << y << endl;
365 progress=(1.0+y)/output.
nrOfRow();
366 pfnProgress(progress,pszMessage,pProgressArg);
374 vector<double> lineOutput(output.
nrOfCol());
376 stat.setNoDataValues(m_noDataValues);
377 const char* pszMessage;
378 void* pProgressArg=NULL;
379 GDALProgressFunc pfnProgress=GDALTermProgress;
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);
393 lineOutput[x]=stat.mymin(pixelInput);
396 lineOutput[x]=stat.mymax(pixelInput);
399 lineOutput[x]=stat.sum(pixelInput);
402 lineOutput[x]=stat.var(pixelInput);
405 lineOutput[x]=sqrt(stat.var(pixelInput));
408 lineOutput[x]=stat.mean(pixelInput);
410 case(filter::percentile):
411 assert(m_threshold.size());
412 lineOutput[x]=stat.percentile(pixelInput,pixelInput.begin(),pixelInput.end(),m_threshold[0]);
415 std::string errorString=
"method not supported";
423 catch(
string errorstring){
424 cerr << errorstring <<
"in line " << y << endl;
426 progress=(1.0+y)/output.
nrOfRow();
427 pfnProgress(progress,pszMessage,pProgressArg);
433 assert(output.
nrOfBand()==methods.size());
438 stat.setNoDataValues(m_noDataValues);
439 const char* pszMessage;
440 void* pProgressArg=NULL;
441 GDALProgressFunc pfnProgress=GDALTermProgress;
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){
450 pixelInput=lineInput.selectCol(x);
453 std::cout <<
"error is caught" << std::endl;
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);
463 case(filter::median):
464 lineOutput[imethod][x]=stat.median(pixelInput);
467 lineOutput[imethod][x]=stat.mymin(pixelInput);
470 lineOutput[imethod][x]=stat.mymax(pixelInput);
473 lineOutput[imethod][x]=stat.sum(pixelInput);
476 lineOutput[imethod][x]=stat.var(pixelInput);
479 lineOutput[imethod][x]=sqrt(stat.var(pixelInput));
482 lineOutput[imethod][x]=stat.mean(pixelInput);
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);
492 std::string errorString=
"method not supported";
499 cerr <<
"An Error in line " << y << endl;
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);
513 const char* pszMessage;
514 void* pProgressArg=NULL;
515 GDALProgressFunc pfnProgress=GDALTermProgress;
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];
532 for(
int iband=0;iband<input.
nrOfBand();++iband){
534 output.
writeData(lineOutput[iband],y,iband);
536 catch(
string errorstring){
537 cerr << errorstring <<
"in band " << iband <<
", line " << y << endl;
540 progress=(1.0+y)/output.
nrOfRow();
541 pfnProgress(progress,pszMessage,pProgressArg);
545void filter::Filter::getSavGolayCoefficients(vector<double> &tapz,
int np,
int nl,
int nr,
int ld,
int m) {
546 int j, k, imj, ipj, kk, mm;
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;
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);
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;
571 for(j = 0; j < m + 1; ++j)
578 for(k = -nl; k <= nr; ++k) {
582 for(mm = 1; mm <= m; ++mm)
583 sum += b[mm] * (fac *= k);
590 tapz.resize(nl+1+nr);
592 tapz[tapz.size()/2]=tmpc[0];
594 for(k=1;k<=tapz.size()/2;++k)
595 tapz[tapz.size()/2-k]=tmpc[k];
597 for(k=1;k<=tapz.size()/2;++k)
598 tapz[tapz.size()/2+k]=tmpc[np-k];
601void filter::Filter::ludcmp(vector<double> &a, vector<int> &indx,
double &d) {
602 const double TINY = 1.0e-20;
603 int i, imax = -1, j, k;
604 double big, dum, sum, temp;
607 vector<double> vv(n, 0.0);
610 for(i = 0; i < n; ++i) {
612 for(j = 0; j < n; ++j)
613 if((temp = fabs(a[i * n + j])) > big) big = temp;
616 cerr <<
"Singular matrix in routine ludcmp" << endl;
622 for(j = 0; j < n; ++j) {
623 for(i = 0; i < j; ++i) {
625 for(k = 0; k < i; ++k)
626 sum -= a[i * n + k] * a[k * n + j];
630 for(i = j; i < n; ++i) {
632 for(k = 0; k < j; ++k)
633 sum -= a[i * n + k] * a[k * n + j];
635 if((dum = vv[i] * fabs(sum)) >= big) {
642 for(k = 0; k < n; ++k) {
643 dum = a[imax * n + k];
644 a[imax * n + k] = a[j * n + k];
651 if(a[j * n + j] == 0.0) a[j * n + j] = TINY;
653 dum = 1. / a[j * n + j];
654 for(i = j + 1; i < n; ++i)
660void filter::Filter::lubksb(vector<double> &a, vector<int> &indx, vector<double> &b) {
661 int i, ii = 0, ip, j;
665 for(i = 0; i < n; ++i) {
670 for(j = ii - 1; j < i; ++j)
671 sum -= a[i * n + j] * b[j];
676 for(i = n - 1; i >= 0; --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];
684double filter::Filter::getCentreWavelength(
const std::vector<double> &wavelengthIn,
const Vector2d<double>& srf,
const std::string& interpolationType,
double delta,
bool verbose)
686 assert(srf.size()==2);
687 int nband=srf[0].size();
688 double start=floor(wavelengthIn[0]);
689 double end=ceil(wavelengthIn.back());
691 std::cout <<
"wavelengths in [" << start <<
"," << end <<
"]" << std::endl << std::flush;
695 gsl_interp_accel *acc;
698 stat.getSpline(interpolationType,nband,spline);
699 stat.initSpline(spline,&(srf[0][0]),&(srf[1][0]),nband);
701 std::cout <<
"calculating norm of srf" << std::endl << std::flush;
703 norm=gsl_spline_eval_integ(spline,srf[0].front(),srf[0].back(),acc);
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);
713 std::cout <<
"interpolate wavelengths to " << wavelength_fine.size() <<
" entries " << std::endl;
714 std::vector<double> srf_fine;
716 stat.interpolateUp(srf[0],srf[1],wavelength_fine,interpolationType,srf_fine,verbose);
717 assert(srf_fine.size()==wavelength_fine.size());
719 gsl_interp_accel *accOut;
720 stat.allocAcc(accOut);
721 gsl_spline *splineOut;
722 stat.getSpline(interpolationType,wavelength_fine.size(),splineOut);
725 std::vector<double> wavelengthOut(wavelength_fine.size());
727 for(
int iband=0;iband<wavelengthOut.size();++iband)
728 wavelengthOut[iband]=wavelength_fine[iband]*srf_fine[iband];
730 stat.initSpline(splineOut,&(wavelength_fine[0]),&(wavelengthOut[0]),wavelength_fine.size());
731 double centreWavelength=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
733 gsl_spline_free(splineOut);
734 gsl_interp_accel_free(accOut);
736 return(centreWavelength);
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.
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...
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...