| pktools 2.6.7 Processing Kernel for geospatial data

Filter.cc
1/**********************************************************************
2 Filter.cc: class for filtering
3 Copyright (C) 2008-2012 Pieter Kempeneers
4
5 This file is part of pktools
6
7 pktools is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10(at your option) any later version.
11
12 pktools is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along 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
25 using 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
39 void 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
54 unsigned 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
60 unsigned int filter::Filter::setNoDataValues(std::vector<double> vnodata){
61 m_noDataValues=vnodata;
62 return m_noDataValues.size();
63};
64
65 void 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
96 void 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
127 void 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
158 void 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
195 void 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
214 void 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
232 void 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
260 void 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
295 void 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
328 void 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
337 void 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
370 void 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());
375 statfactory::StatFactory stat;
376 stat.setNoDataValues(m_noDataValues);
377 const char* pszMessage;
378 void* pProgressArg=NULL;
379 GDALProgressFunc pfnProgress=GDALTermProgress;
380 double progress=0;
381 pfnProgress(progress,pszMessage,pProgressArg);
382 for(int y=0;y<input.nrOfRow();++y){
383 for(int iband=0;iband<input.nrOfBand();++iband)
384 input.readData(lineInput[iband],y,iband);
385 vector<double> pixelInput(input.nrOfBand());
386 for(int x=0;x<input.nrOfCol();++x){
387 pixelInput=lineInput.selectCol(x);
388 switch(getFilterType(method)){
389 case(filter::median):
390 lineOutput[x]=stat.median(pixelInput);
391 break;
392 case(filter::min):
393 lineOutput[x]=stat.mymin(pixelInput);
394 break;
395 case(filter::max):
396 lineOutput[x]=stat.mymax(pixelInput);
397 break;
398 case(filter::sum):
399 lineOutput[x]=stat.sum(pixelInput);
400 break;
401 case(filter::var):
402 lineOutput[x]=stat.var(pixelInput);
403 break;
404 case(filter::stdev):
405 lineOutput[x]=sqrt(stat.var(pixelInput));
406 break;
407 case(filter::mean):
408 lineOutput[x]=stat.mean(pixelInput);
409 break;
410 case(filter::percentile):
411 assert(m_threshold.size());
412 lineOutput[x]=stat.percentile(pixelInput,pixelInput.begin(),pixelInput.end(),m_threshold[0]);
413 break;
414 default:
415 std::string errorString="method not supported";
416 throw(errorString);
417 break;
418 }
419 }
420 try{
421 output.writeData(lineOutput,y);
422 }
423 catch(string errorstring){
424 cerr << errorstring << "in line " << y << endl;
425 }
426 progress=(1.0+y)/output.nrOfRow();
427 pfnProgress(progress,pszMessage,pProgressArg);
428 }
429}
430
431 void 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());
437 statfactory::StatFactory stat;
438 stat.setNoDataValues(m_noDataValues);
439 const char* pszMessage;
440 void* pProgressArg=NULL;
441 GDALProgressFunc pfnProgress=GDALTermProgress;
442 double progress=0;
443 pfnProgress(progress,pszMessage,pProgressArg);
444 for(int y=0;y<input.nrOfRow();++y){
445 for(int iband=0;iband<input.nrOfBand();++iband)
446 input.readData(lineInput[iband],y,iband);
447 vector<double> pixelInput(input.nrOfBand());
448 for(int x=0;x<input.nrOfCol();++x){
449 try{
450 pixelInput=lineInput.selectCol(x);
451 }
452 catch(...){
453 std::cout << "error is caught" << std::endl;
454 exit(1);
455 }
456 int ithreshold=0;//threshold to use for percentiles
457 try{
458 for(int imethod=0;imethod<methods.size();++imethod){
459 switch(getFilterType(methods[imethod])){
460 case(filter::nvalid):
461 lineOutput[imethod][x]=stat.nvalid(pixelInput);
462 break;
463 case(filter::median):
464 lineOutput[imethod][x]=stat.median(pixelInput);
465 break;
466 case(filter::min):
467 lineOutput[imethod][x]=stat.mymin(pixelInput);
468 break;
469 case(filter::max):
470 lineOutput[imethod][x]=stat.mymax(pixelInput);
471 break;
472 case(filter::sum):
473 lineOutput[imethod][x]=stat.sum(pixelInput);
474 break;
475 case(filter::var):
476 lineOutput[imethod][x]=stat.var(pixelInput);
477 break;
478 case(filter::stdev):
479 lineOutput[imethod][x]=sqrt(stat.var(pixelInput));
480 break;
481 case(filter::mean):
482 lineOutput[imethod][x]=stat.mean(pixelInput);
483 break;
484 case(filter::percentile):{
485 assert(m_threshold.size());
486 double threshold=(ithreshold<m_threshold.size())? m_threshold[ithreshold] : m_threshold[0];
487 lineOutput[imethod][x]=stat.percentile(pixelInput,pixelInput.begin(),pixelInput.end(),threshold);
488 ++ithreshold;
489 break;
490 }
491 default:
492 std::string errorString="method not supported";
493 throw(errorString);
494 break;
495 }
496 }
497 }
498 catch(...){
499 cerr << "An Error in line " << y << endl;
500 }
501 }
502 for(int imethod=0;imethod<methods.size();++imethod)
503 output.writeData(lineOutput[imethod],y,imethod);
504 progress=(1.0+y)/output.nrOfRow();
505 pfnProgress(progress,pszMessage,pProgressArg);
506 }
507}
508
509 void 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
545 void 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
601 void 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
660 void 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
684 double 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
693 statfactory::StatFactory stat;
694
695 gsl_interp_accel *acc;
696 stat.allocAcc(acc);
697 gsl_spline *spline;
698 stat.getSpline(interpolationType,nband,spline);
699 stat.initSpline(spline,&(srf[0][0]),&(srf[1][0]),nband);
700 if(verbose)
701 std::cout << "calculating norm of srf" << std::endl << std::flush;
702 double norm=0;
703 norm=gsl_spline_eval_integ(spline,srf[0].front(),srf[0].back(),acc);
704 if(verbose)
705 std::cout << "norm of srf: " << norm << std::endl << std::flush;
706 gsl_spline_free(spline);
707 gsl_interp_accel_free(acc);
708 std::vector<double> wavelength_fine;
709 for(double win=floor(wavelengthIn[0]);win<=ceil(wavelengthIn.back());win+=delta)
710 wavelength_fine.push_back(win);
711
712 if(verbose)
713 std::cout << "interpolate wavelengths to " << wavelength_fine.size() << " entries " << std::endl;
714 std::vector<double> srf_fine;//spectral response function, interpolated for wavelength_fine
715
716 stat.interpolateUp(srf[0],srf[1],wavelength_fine,interpolationType,srf_fine,verbose);
717 assert(srf_fine.size()==wavelength_fine.size());
718
719 gsl_interp_accel *accOut;
720 stat.allocAcc(accOut);
721 gsl_spline *splineOut;
722 stat.getSpline(interpolationType,wavelength_fine.size(),splineOut);
723 assert(splineOut);
724
725 std::vector<double> wavelengthOut(wavelength_fine.size());
726
727 for(int iband=0;iband<wavelengthOut.size();++iband)
728 wavelengthOut[iband]=wavelength_fine[iband]*srf_fine[iband];
729
730 stat.initSpline(splineOut,&(wavelength_fine[0]),&(wavelengthOut[0]),wavelength_fine.size());
731 double centreWavelength=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
732
733 gsl_spline_free(splineOut);
734 gsl_interp_accel_free(accOut);
735
736 return(centreWavelength);
737}
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// }
ImgRasterGdal::nrOfRow
int nrOfRow(void) const
Get the number of rows of this dataset.
Definition: ImgRasterGdal.h:100
ImgRasterGdal::nrOfBand
int nrOfBand(void) const
Get the number of bands of this dataset.
Definition: ImgRasterGdal.h:102
ImgRasterGdal::nrOfCol
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
ImgReaderGdal
Definition: ImgReaderGdal.h:37
ImgReaderGdal::readData
void readData(T &value, int col, int row, int band=0)
Read a single pixel cell value at a specific column and row for a specific band (all indices start co...
Definition: ImgReaderGdal.h:95
ImgWriterGdal
Definition: ImgWriterGdal.h:38
ImgWriterGdal::writeData
bool writeData(T &value, int col, int row, int band=0)
Write a single pixel cell value at a specific column and row for a specific band (all indices start c...
Definition: ImgWriterGdal.h:96
Vector2d< double >
statfactory::StatFactory
Definition: StatFactory.h:43