pktools 2.6.7
Processing Kernel for geospatial data
Filter.cc
1/**********************************************************************
2Filter.cc: class for filtering
3Copyright (C) 2008-2012 Pieter Kempeneers
4
5This file is part of pktools
6
7pktools is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published by
9the Free Software Foundation, either version 3 of the License, or
10(at your option) any later version.
11
12pktools is distributed in the hope that it will be useful,
13but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15GNU General Public License for more details.
16
17You should have received a copy of the GNU General Public License
18along with pktools. If not, see <http://www.gnu.org/licenses/>.
19***********************************************************************/
20#include "Filter.h"
21#include <assert.h>
22#include <math.h>
23#include <iostream>
24
25using namespace std;
26
27filter::Filter::Filter(void)
28 : m_padding("symmetric")
29{
30}
31
32
33filter::Filter::Filter(const vector<double> &taps)
34 : m_padding("symmetric")
35{
36 setTaps(taps);
37}
38
39void filter::Filter::setTaps(const vector<double> &taps, bool normalize)
40{
41 m_taps.resize(taps.size());
42 double norm=0;
43 for(int itap=0;itap<taps.size();++itap)
44 norm+=taps[itap];
45 if(norm){
46 for(int itap=0;itap<taps.size();++itap)
47 m_taps[itap]=taps[itap]/norm;
48 }
49 else
50 m_taps=taps;
51 assert(m_taps.size()%2);
52}
53
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();
58};
59
60unsigned int filter::Filter::setNoDataValues(std::vector<double> vnodata){
61 m_noDataValues=vnodata;
62 return m_noDataValues.size();
63};
64
65void filter::Filter::dwtForward(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family){
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}
95
96void filter::Filter::dwtInverse(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family){
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}
126
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;
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}
157
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;
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}
193
194//todo: support different padding strategies
195void filter::Filter::dwtForward(std::vector<double>& data, const std::string& wavelet_type, int family){
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}
212
213//todo: support different padding strategies
214void filter::Filter::dwtInverse(std::vector<double>& data, const std::string& wavelet_type, int family){
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}
230
231//todo: support different padding strategies
232void filter::Filter::dwtCut(std::vector<double>& data, const std::string& wavelet_type, int family, double cut){
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}
259
260void filter::Filter::morphology(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short verbose)
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}
294
295void filter::Filter::smoothNoData(ImgReaderGdal& input, const std::string& interpolationType, ImgWriterGdal& output)
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}
327
328void filter::Filter::smooth(ImgReaderGdal& input, ImgWriterGdal& output, short dim)
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}
336
337void filter::Filter::filter(ImgReaderGdal& input, ImgWriterGdal& output)
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}
369
370void filter::Filter::stat(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method)
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}
430
431void filter::Filter::stats(ImgReaderGdal& input, ImgWriterGdal& output, const vector<std::string>& methods)
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}
508
509void filter::Filter::filter(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim)
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}
544
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;
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}
600
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;
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}
659
660void filter::Filter::lubksb(vector<double> &a, vector<int> &indx, vector<double> &b) {
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}
683
684double filter::Filter::getCentreWavelength(const std::vector<double> &wavelengthIn, const Vector2d<double>& srf, const std::string& interpolationType, double delta, bool verbose)
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}
738
739// void filter::Filter::applyFwhm(const vector<double> &wavelengthIn, const ImgReaderGdal& input, const vector<double> &wavelengthOut, const vector<double> &fwhm, const std::string& interpolationType, ImgWriterGdal& output, bool verbose){
740// Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
741// Vector2d<double> lineOutput(wavelengthOut.size(),input.nrOfCol());
742// const char* pszMessage;
743// void* pProgressArg=NULL;
744// GDALProgressFunc pfnProgress=GDALTermProgress;
745// double progress=0;
746// pfnProgress(progress,pszMessage,pProgressArg);
747// for(int y=0;y<input.nrOfRow();++y){
748// for(int iband=0;iband<input.nrOfBand();++iband)
749// input.readData(lineInput[iband],GDT_Float64,y,iband);
750// applyFwhm<double>(wavelengthIn,lineInput,wavelengthOut,fwhm, interpolationType, lineOutput, verbose);
751// for(int iband=0;iband<output.nrOfBand();++iband){
752// try{
753// output.writeData(lineOutput[iband],GDT_Float64,y,iband);
754// }
755// catch(string errorstring){
756// cerr << errorstring << "in band " << iband << ", line " << y << endl;
757// }
758// }
759// progress=(1.0+y)/output.nrOfRow();
760// pfnProgress(progress,pszMessage,pProgressArg);
761// }
762// }
763
764// void filter::Filter::applySrf(const vector<double> &wavelengthIn, const ImgReaderGdal& input, const vector< Vector2d<double> > &srf, const std::string& interpolationType, ImgWriterGdal& output, bool verbose){
765// assert(output.nrOfBand()==srf.size());
766// double centreWavelength=0;
767// Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
768// const char* pszMessage;
769// void* pProgressArg=NULL;
770// GDALProgressFunc pfnProgress=GDALTermProgress;
771// double progress=0;
772// pfnProgress(progress,pszMessage,pProgressArg);
773// for(int y=0;y<input.nrOfRow();++y){
774// for(int iband=0;iband<input.nrOfBand();++iband)
775// input.readData(lineInput[iband],GDT_Float64,y,iband);
776// for(int isrf=0;isrf<srf.size();++isrf){
777// vector<double> lineOutput(input.nrOfCol());
778// centreWavelength=applySrf<double>(wavelengthIn,lineInput,srf[isrf], interpolationType, lineOutput, verbose);
779// for(int iband=0;iband<output.nrOfBand();++iband){
780// try{
781// output.writeData(lineOutput,GDT_Float64,y,isrf);
782// }
783// catch(string errorstring){
784// cerr << errorstring << "in band " << iband << ", line " << y << endl;
785// }
786// }
787// }
788// progress=(1.0+y)/output.nrOfRow();
789// pfnProgress(progress,pszMessage,pProgressArg);
790// }
791// }
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