pktools 2.6.7
Processing Kernel for geospatial data
Filter2d.cc
1/**********************************************************************
2Filter2d.cc: class for filtering images
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 <sstream>
21#include <iomanip>
22#include <iostream>
23#include <cmath>
24#include "Filter2d.h"
25#include "StatFactory.h"
26// #include "imageclasses/ImgUtils.h"
27
28filter2d::Filter2d::Filter2d(void)
29{
30}
31
32filter2d::Filter2d::Filter2d(const Vector2d<double> &taps)
33 : m_taps(taps)
34{
35}
36
37int filter2d::Filter2d::pushNoDataValue(double noDataValue)
38{
39 if(find(m_noDataValues.begin(),m_noDataValues.end(),noDataValue)==m_noDataValues.end())
40 m_noDataValues.push_back(noDataValue);
41 return(m_noDataValues.size());
42}
43
44void filter2d::Filter2d::setTaps(const Vector2d<double> &taps)
45{
46 m_taps=taps;
47}
48
49void filter2d::Filter2d::smoothNoData(ImgReaderGdal& input, ImgWriterGdal& output, int dim)
50{
51 smoothNoData(input, output,dim,dim);
52}
53
54void filter2d::Filter2d::smooth(ImgReaderGdal& input, ImgWriterGdal& output, int dim)
55{
56 smooth(input, output,dim,dim);
57}
58
59void filter2d::Filter2d::smoothNoData(ImgReaderGdal& input, ImgWriterGdal& output, int dimX, int dimY)
60{
61 m_taps.resize(dimY);
62 for(int j=0;j<dimY;++j){
63 m_taps[j].resize(dimX);
64 for(int i=0;i<dimX;++i)
65 m_taps[j][i]=1.0;
66 }
67 filter(input,output,false,true,true);
68}
69
70void filter2d::Filter2d::smooth(ImgReaderGdal& input, ImgWriterGdal& output, int dimX, int dimY)
71{
72 m_taps.resize(dimY);
73 for(int j=0;j<dimY;++j){
74 m_taps[j].resize(dimX);
75 for(int i=0;i<dimX;++i)
76 m_taps[j][i]=1.0;
77 }
78 filter(input,output,false,true,false);
79}
80
81
82void filter2d::Filter2d::filter(ImgReaderGdal& input, ImgWriterGdal& output, bool absolute, bool normalize, bool noData)
83{
84 int dimX=m_taps[0].size();//horizontal!!!
85 int dimY=m_taps.size();//vertical!!!
86 // byte* tmpbuf=new byte[input.rowSize()];
87 const char* pszMessage;
88 void* pProgressArg=NULL;
89 GDALProgressFunc pfnProgress=GDALTermProgress;
90 double progress=0;
91 pfnProgress(progress,pszMessage,pProgressArg);
92 for(int iband=0;iband<input.nrOfBand();++iband){
93 Vector2d<double> inBuffer(dimY,input.nrOfCol());
94 std::vector<double> outBuffer(input.nrOfCol());
95 int indexI=0;
96 int indexJ=0;
97 //initialize last half of inBuffer
98 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
99 try{
100 input.readData(inBuffer[indexJ],abs(j),iband);
101 }
102 catch(std::string errorstring){
103 std::cerr << errorstring << "in line " << indexJ << std::endl;
104 }
105 ++indexJ;
106 }
107
108 for(int y=0;y<input.nrOfRow();++y){
109 if(y){//inBuffer already initialized for y=0
110 //erase first line from inBuffer
111 if(dimY>1)
112 inBuffer.erase(inBuffer.begin());
113 //read extra line and push back to inBuffer if not out of bounds
114 if(y+dimY/2<input.nrOfRow()){
115 //allocate buffer
116 if(dimY>1)
117 inBuffer.push_back(inBuffer.back());
118 try{
119 input.readData(inBuffer[inBuffer.size()-1],y+dimY/2,iband);
120 }
121 catch(std::string errorstring){
122 std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
123 }
124 }
125 else{
126 int over=y+dimY/2-input.nrOfRow();
127 int index=(inBuffer.size()-1)-over;
128 assert(index>=0);
129 assert(index<inBuffer.size());
130 inBuffer.push_back(inBuffer[index]);
131 }
132 }
133 for(int x=0;x<input.nrOfCol();++x){
134 outBuffer[x]=0;
135 double norm=0;
136 bool masked=false;
137 if(noData){//only filter noData values
138 for(int imask=0;imask<m_noDataValues.size();++imask){
139 if(inBuffer[(dimY-1)/2][x]==m_noDataValues[imask]){
140 masked=true;
141 break;
142 }
143 }
144 if(!masked){
145 outBuffer[x]=inBuffer[(dimY-1)/2][x];
146 continue;
147 }
148 }
149 assert(!noData||masked);
150 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
151 for(int i=-(dimX-1)/2;i<=dimX/2;++i){
152 indexI=x+i;
153 indexJ=(dimY-1)/2+j;
154 //check if out of bounds
155 if(x<(dimX-1)/2)
156 indexI=x+abs(i);
157 else if(x>=input.nrOfCol()-(dimX-1)/2)
158 indexI=x-abs(i);
159 if(y<(dimY-1)/2)
160 indexJ=(dimY-1)/2+abs(j);
161 else if(y>=input.nrOfRow()-(dimY-1)/2)
162 indexJ=(dimY-1)/2-abs(j);
163 //do not take masked values into account
164 masked=false;
165 for(int imask=0;imask<m_noDataValues.size();++imask){
166 if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
167 masked=true;
168 break;
169 }
170 }
171 if(!masked){
172 outBuffer[x]+=(m_taps[(dimY-1)/2+j][(dimX-1)/2+i]*inBuffer[indexJ][indexI]);
173 norm+=m_taps[(dimY-1)/2+j][(dimX-1)/2+i];
174 }
175 }
176 }
177 if(absolute)
178 outBuffer[x]=(normalize&&norm)? abs(outBuffer[x])/norm : abs(outBuffer[x]);
179 else if(normalize&&norm!=0)
180 outBuffer[x]=outBuffer[x]/norm;
181 }
182 //write outBuffer to file
183 try{
184 output.writeData(outBuffer,y,iband);
185 }
186 catch(std::string errorstring){
187 std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
188 }
189 progress=(1.0+y);
190 progress+=(output.nrOfRow()*iband);
191 progress/=output.nrOfBand()*output.nrOfRow();
192 pfnProgress(progress,pszMessage,pProgressArg);
193 }
194 }
195}
196
197
198void filter2d::Filter2d::majorVoting(const std::string& inputFilename, const std::string& outputFilename,int dim,const std::vector<int> &prior)
199{
200 const char* pszMessage;
201 void* pProgressArg=NULL;
202 GDALProgressFunc pfnProgress=GDALTermProgress;
203 double progress=0;
204 pfnProgress(progress,pszMessage,pProgressArg);
205
206 bool usePriors=true;
207 if(prior.empty()){
208 std::cout << "no prior information" << std::endl;
209 usePriors=false;
210 }
211 else{
212 std::cout << "using priors ";
213 for(int iclass=0;iclass<prior.size();++iclass)
214 std::cout << " " << static_cast<short>(prior[iclass]);
215 std::cout << std::endl;
216 }
217
218 ImgReaderGdal input;
219 ImgWriterGdal output;
220 input.open(inputFilename);
221 output.open(outputFilename,input);
222 int dimX=0;//horizontal!!!
223 int dimY=0;//vertical!!!
224 if(dim){
225 dimX=dim;
226 dimY=dim;
227 }
228 else{
229 dimX=m_taps[0].size();
230 dimY=m_taps.size();
231 }
232
233 assert(dimX);
234 assert(dimY);
235
236 Vector2d<double> inBuffer(dimY,input.nrOfCol());
237 std::vector<double> outBuffer(input.nrOfCol());
238 int indexI=0;
239 int indexJ=0;
240 //initialize last half of inBuffer
241 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
242 try{
243 input.readData(inBuffer[indexJ],abs(j));
244 }
245 catch(std::string errorstring){
246 std::cerr << errorstring << "in line " << indexJ << std::endl;
247 }
248 ++indexJ;
249 }
250
251 for(int y=0;y<input.nrOfRow();++y){
252 if(y){//inBuffer already initialized for y=0
253 //erase first line from inBuffer
254 if(dimY>1)
255 inBuffer.erase(inBuffer.begin());
256 //read extra line and push back to inBuffer if not out of bounds
257 if(y+dimY/2<input.nrOfRow()){
258 //allocate buffer
259 if(dimY>1)
260 inBuffer.push_back(inBuffer.back());
261 try{
262 input.readData(inBuffer[inBuffer.size()-1],y+dimY/2);
263 }
264 catch(std::string errorstring){
265 std::cerr << errorstring << "in line" << y << std::endl;
266 }
267 }
268 else{
269 int over=y+dimY/2-input.nrOfRow();
270 int index=(inBuffer.size()-1)-over;
271 assert(index>=0);
272 assert(index<inBuffer.size());
273 inBuffer.push_back(inBuffer[index]);
274 }
275 }
276 for(int x=0;x<input.nrOfCol();++x){
277 outBuffer[x]=0;
278 std::map<int,int> occurrence;
279 int centre=dimX*(dimY-1)/2+(dimX-1)/2;
280 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
281 for(int i=-(dimX-1)/2;i<=dimX/2;++i){
282 indexI=x+i;
283 //check if out of bounds
284 if(indexI<0)
285 indexI=-indexI;
286 else if(indexI>=input.nrOfCol())
287 indexI=input.nrOfCol()-i;
288 if(y+j<0)
289 indexJ=-j;
290 else if(y+j>=input.nrOfRow())
291 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
292 else
293 indexJ=(dimY-1)/2+j;
294
295 // if(x<dimX/2)
296 // indexI=x+abs(i);
297 // else if(x>=input.nrOfCol()-dimX/2)
298 // indexI=x-abs(i);
299 // if(y<dimY/2)
300 // indexJ=dimY/2+abs(j);
301 // else if(y>=input.nrOfRow()-dimY/2)
302 // indexJ=dimY/2-abs(j);
303 if(usePriors){
304 occurrence[inBuffer[indexJ][indexI]]+=prior[inBuffer[indexJ][indexI]-1];
305 }
306 else
307 ++occurrence[inBuffer[indexJ][indexI]];
308 }
309 }
310 std::map<int,int>::const_iterator maxit=occurrence.begin();
311 for(std::map<int,int>::const_iterator mit=occurrence.begin();mit!=occurrence.end();++mit){
312 if(mit->second>maxit->second)
313 maxit=mit;
314 }
315 if(occurrence[inBuffer[(dimY-1)/2][x]]<maxit->second)//
316 outBuffer[x]=maxit->first;
317 else//favorize original value in case of ties
318 outBuffer[x]=inBuffer[(dimY-1)/2][x];
319 }
320 //write outBuffer to file
321 try{
322 output.writeData(outBuffer,y);
323 }
324 catch(std::string errorstring){
325 std::cerr << errorstring << "in line" << y << std::endl;
326 }
327 progress=(1.0+y)/output.nrOfRow();
328 pfnProgress(progress,pszMessage,pProgressArg);
329 }
330 input.close();
331 output.close();
332}
333
334void filter2d::Filter2d::median(const std::string& inputFilename, const std::string& outputFilename,int dim, bool disc)
335{
336 ImgReaderGdal input;
337 ImgWriterGdal output;
338 input.open(inputFilename);
339 output.open(outputFilename,input);
340 doit(input,output,"median",dim,disc);
341}
342
343void filter2d::Filter2d::var(const std::string& inputFilename, const std::string& outputFilename,int dim, bool disc)
344{
345 ImgReaderGdal input;
346 ImgWriterGdal output;
347 input.open(inputFilename);
348 output.open(outputFilename,input);
349 doit(input,output,"var",dim,disc);
350}
351
352void filter2d::Filter2d::doit(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down, bool disc)
353{
354 doit(input,output,method,dim,dim,down,disc);
355}
356
357void filter2d::Filter2d::doit(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, short down, bool disc)
358{
359 const char* pszMessage;
360 void* pProgressArg=NULL;
361 GDALProgressFunc pfnProgress=GDALTermProgress;
362 double progress=0;
363 pfnProgress(progress,pszMessage,pProgressArg);
364
365 assert(dimX);
366 assert(dimY);
367
369 for(int iband=0;iband<input.nrOfBand();++iband){
370 Vector2d<double> inBuffer(dimY,input.nrOfCol());
371 std::vector<double> outBuffer((input.nrOfCol()+down-1)/down);
372 int indexI=0;
373 int indexJ=0;
374 //initialize last half of inBuffer
375 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
376 try{
377 input.readData(inBuffer[indexJ],abs(j),iband);
378 }
379 catch(std::string errorstring){
380 std::cerr << errorstring << "in line " << indexJ << std::endl;
381 }
382 ++indexJ;
383 }
384 for(int y=0;y<input.nrOfRow();++y){
385 if(y){//inBuffer already initialized for y=0
386 //erase first line from inBuffer
387 if(dimY>1)
388 inBuffer.erase(inBuffer.begin());
389 //read extra line and push back to inBuffer if not out of bounds
390 if(y+dimY/2<input.nrOfRow()){
391 //allocate buffer
392 if(dimY>1)
393 inBuffer.push_back(inBuffer.back());
394 try{
395 input.readData(inBuffer[inBuffer.size()-1],y+dimY/2,iband);
396 }
397 catch(std::string errorstring){
398 std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
399 }
400 }
401 else{
402 int over=y+dimY/2-input.nrOfRow();
403 int index=(inBuffer.size()-1)-over;
404 assert(index>=0);
405 assert(index<inBuffer.size());
406 inBuffer.push_back(inBuffer[index]);
407 }
408 }
409 if((y+1+down/2)%down)
410 continue;
411 for(int x=0;x<input.nrOfCol();++x){
412 if((x+1+down/2)%down)
413 continue;
414 outBuffer[x/down]=0;
415 std::vector<double> windowBuffer;
416 std::map<long int,int> occurrence;
417 int centre=dimX*(dimY-1)/2+(dimX-1)/2;
418 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
419 for(int i=-(dimX-1)/2;i<=dimX/2;++i){
420 double d2=i*i+j*j;//square distance
421 if(disc&&(d2>(dimX/2)*(dimY/2)))
422 continue;
423 indexI=x+i;
424 //check if out of bounds
425 if(indexI<0)
426 indexI=-indexI;
427 else if(indexI>=input.nrOfCol())
428 indexI=input.nrOfCol()-i;
429 if(y+j<0)
430 indexJ=-j;
431 else if(y+j>=input.nrOfRow())
432 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
433 else
434 indexJ=(dimY-1)/2+j;
435 bool masked=false;
436 for(int imask=0;imask<m_noDataValues.size();++imask){
437 if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
438 masked=true;
439 break;
440 }
441 }
442 if(!masked){
443 std::vector<short>::const_iterator vit=m_class.begin();
444 if(!m_class.size())
445 ++occurrence[inBuffer[indexJ][indexI]];
446 else{
447 while(vit!=m_class.end()){
448 if(inBuffer[indexJ][indexI]==*(vit++))
449 ++occurrence[inBuffer[indexJ][indexI]];
450 }
451 }
452 windowBuffer.push_back(inBuffer[indexJ][indexI]);
453 }
454 }
455 }
456 switch(getFilterType(method)){
457 case(filter2d::nvalid):
458 outBuffer[x/down]=stat.nvalid(windowBuffer);
459 break;
460 case(filter2d::median):
461 if(windowBuffer.empty())
462 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
463 else
464 outBuffer[x/down]=stat.median(windowBuffer);
465 break;
466 case(filter2d::var):{
467 if(windowBuffer.empty())
468 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
469 else
470 outBuffer[x/down]=stat.var(windowBuffer);
471 break;
472 }
473 case(filter2d::stdev):{
474 if(windowBuffer.empty())
475 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
476 else
477 outBuffer[x/down]=sqrt(stat.var(windowBuffer));
478 break;
479 }
480 case(filter2d::mean):{
481 if(windowBuffer.empty())
482 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
483 else
484 outBuffer[x/down]=stat.mean(windowBuffer);
485 break;
486 }
487 case(filter2d::min):{
488 if(windowBuffer.empty())
489 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
490 else
491 outBuffer[x/down]=stat.mymin(windowBuffer);
492 break;
493 }
494 case(filter2d::ismin):{
495 if(windowBuffer.empty())
496 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
497 else
498 outBuffer[x/down]=(stat.mymin(windowBuffer)==windowBuffer[centre])? 1:0;
499 break;
500 }
501 case(filter2d::minmax):{//is the same as homog?
502 double min=0;
503 double max=0;
504 if(windowBuffer.empty())
505 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
506 else{
507 stat.minmax(windowBuffer,windowBuffer.begin(),windowBuffer.end(),min,max);
508 if(min!=max)
509 outBuffer[x/down]=0;
510 else
511 outBuffer[x/down]=windowBuffer[centre];//centre pixels
512 }
513 break;
514 }
515 case(filter2d::max):{
516 if(windowBuffer.empty())
517 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
518 else
519 outBuffer[x/down]=stat.mymax(windowBuffer);
520 break;
521 }
522 case(filter2d::ismax):{
523 if(windowBuffer.empty())
524 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
525 else
526 outBuffer[x/down]=(stat.mymax(windowBuffer)==windowBuffer[centre])? 1:0;
527 break;
528 }
529 case(filter2d::order):{
530 if(windowBuffer.empty())
531 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
532 else{
533 double lbound=0;
534 double ubound=dimX*dimY;
535 double theMin=stat.mymin(windowBuffer);
536 double theMax=stat.mymax(windowBuffer);
537 double scale=(ubound-lbound)/(theMax-theMin);
538 outBuffer[x/down]=static_cast<short>(scale*(windowBuffer[centre]-theMin)+lbound);
539 }
540 break;
541 }
542 case(filter2d::sum):{
543 outBuffer[x/down]=stat.sum(windowBuffer);
544 break;
545 }
546 case(filter2d::percentile):{
547 assert(m_threshold.size());
548 outBuffer[x/down]=stat.percentile(windowBuffer,windowBuffer.begin(),windowBuffer.end(),m_threshold[0]);
549 break;
550 }
551 case(filter2d::proportion):{
552 if(windowBuffer.size()){
553 double sum=stat.sum(windowBuffer);
554 if(sum)
555 outBuffer[x/down]=100.0*windowBuffer[centre]/stat.sum(windowBuffer);
556 else
557 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
558 }
559 else
560 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
561 break;
562 }
563 case(filter2d::homog):
564 if(occurrence.size()==1)//all values in window are the same
565 outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
566 else
567 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
568 break;
569 case(filter2d::heterog):{
570 if(occurrence.size()==windowBuffer.size())
571 outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
572 else
573 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
574 // if(occurrence.size()==1)//all values in window are the same
575 // outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
576 // else
577 // outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
578 // break;
579 // for(std::vector<double>::const_iterator wit=windowBuffer.begin();wit!=windowBuffer.end();++wit){
580 // if(wit==windowBuffer.begin()+windowBuffer.size()/2)
581 // continue;
582 // else if(*wit!=inBuffer[(dimY-1)/2][x]){
583 // outBuffer[x/down]=1;
584 // break;
585 // }
586 // else if(*wit==inBuffer[(dimY-1)/2][x]){//todo:wit mag niet central pixel zijn
587 // outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
588 // break;
589 // }
590 // }
591 // break;
592 }
593 case(filter2d::density):{
594 if(windowBuffer.size()){
595 std::vector<short>::const_iterator vit=m_class.begin();
596 while(vit!=m_class.end())
597 outBuffer[x/down]+=100.0*occurrence[*(vit++)]/windowBuffer.size();
598 }
599 else
600 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
601 break;
602 }
603 case(filter2d::countid):{
604 if(windowBuffer.size())
605 outBuffer[x/down]=occurrence.size();
606 else
607 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
608 break;
609 }
610 case(filter2d::mode):{
611 if(occurrence.size()){
612 std::map<long int,int>::const_iterator maxit=occurrence.begin();
613 for(std::map<long int,int>::const_iterator mit=occurrence.begin();mit!=occurrence.end();++mit){
614 if(mit->second>maxit->second)
615 maxit=mit;
616 }
617 if(occurrence[inBuffer[(dimY-1)/2][x]]<maxit->second)//
618 outBuffer[x/down]=maxit->first;
619 else//favorize original value in case of ties
620 outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
621 }
622 else
623 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
624 break;
625 }
626 case(filter2d::threshold):{
627 assert(m_class.size()==m_threshold.size());
628 if(windowBuffer.size()){
629 outBuffer[x/down]=inBuffer[(dimY-1)/2][x];//initialize with original value (in case thresholds not met)
630 for(int iclass=0;iclass<m_class.size();++iclass){
631 if(100.0*(occurrence[m_class[iclass]])/windowBuffer.size()>m_threshold[iclass])
632 outBuffer[x/down]=m_class[iclass];
633 }
634 }
635 else
636 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
637 break;
638 }
639 case(filter2d::scramble):{//could be done more efficiently window by window with random shuffling entire buffer and assigning entire buffer at once to output image...
640 if(windowBuffer.size()){
641 int randomIndex=std::rand()%windowBuffer.size();
642 if(randomIndex>=windowBuffer.size())
643 outBuffer[x/down]=windowBuffer.back();
644 else if(randomIndex<0)
645 outBuffer[x/down]=windowBuffer[0];
646 else
647 outBuffer[x/down]=windowBuffer[randomIndex];
648 }
649 else
650 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
651 break;
652 }
653 case(filter2d::mixed):{
654 enum Type { BF=11, CF=12, MF=13, NF=20, W=30 };
655 double nBF=occurrence[BF];
656 double nCF=occurrence[CF];
657 double nMF=occurrence[MF];
658 double nNF=occurrence[NF];
659 double nW=occurrence[W];
660 if(windowBuffer.size()){
661 if((nBF+nCF+nMF)&&(nBF+nCF+nMF>=nNF+nW)){//forest
662 if(nBF/(nBF+nCF)>=0.75)
663 outBuffer[x/down]=BF;
664 else if(nCF/(nBF+nCF)>=0.75)
665 outBuffer[x/down]=CF;
666 else
667 outBuffer[x/down]=MF;
668 }
669 else{//non-forest
670 if(nW&&(nW>=nNF))
671 outBuffer[x/down]=W;
672 else
673 outBuffer[x/down]=NF;
674 }
675 }
676 else
677 outBuffer[x/down]=inBuffer[indexJ][indexI];
678 break;
679 }
680 default:{
681 std::ostringstream ess;
682 ess << "Error: filter method " << method << " not supported" << std::endl;
683 throw(ess.str());
684 break;
685 }
686 }
687 }
688 progress=(1.0+y/down);
689 progress+=(output.nrOfRow()*iband);
690 progress/=output.nrOfBand()*output.nrOfRow();
691 pfnProgress(progress,pszMessage,pProgressArg);
692 //write outBuffer to file
693 try{
694 output.writeData(outBuffer,y/down,iband);
695 }
696 catch(std::string errorstring){
697 std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
698 }
699 }
700 }
701 pfnProgress(1.0,pszMessage,pProgressArg);
702}
703
704void filter2d::Filter2d::mrf(ImgReaderGdal& input, ImgWriterGdal& output, int dimX, int dimY, double beta, bool eightConnectivity, short down, bool verbose){
705 assert(m_class.size()>1);
706 Vector2d<double> fullBeta(m_class.size(),m_class.size());
707 for(int iclass1=0;iclass1<m_class.size();++iclass1)
708 for(int iclass2=0;iclass2<m_class.size();++iclass2)
709 fullBeta[iclass1][iclass2]=beta;
710 mrf(input,output,dimX,dimY,fullBeta,eightConnectivity,down,verbose);
711}
712
713//beta[classTo][classFrom]
714void filter2d::Filter2d::mrf(ImgReaderGdal& input, ImgWriterGdal& output, int dimX, int dimY, Vector2d<double> beta, bool eightConnectivity, short down, bool verbose)
715{
716 const char* pszMessage;
717 void* pProgressArg=NULL;
718 GDALProgressFunc pfnProgress=GDALTermProgress;
719 double progress=0;
720 pfnProgress(progress,pszMessage,pProgressArg);
721
722 assert(dimX);
723 assert(dimY);
724
725 Vector2d<short> inBuffer(dimY,input.nrOfCol());
726 Vector2d<double> outBuffer(m_class.size(),(input.nrOfCol()+down-1)/down);
727 assert(input.nrOfBand()==1);
728 assert(output.nrOfBand()==m_class.size());
729 assert(m_class.size()>1);
730 assert(beta.size()==m_class.size());
731 int indexI=0;
732 int indexJ=0;
733 //initialize last half of inBuffer
734 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
735 try{
736 input.readData(inBuffer[indexJ],abs(j));
737 }
738 catch(std::string errorstring){
739 std::cerr << errorstring << "in line " << indexJ << std::endl;
740 }
741 ++indexJ;
742 }
743 for(int y=0;y<input.nrOfRow();++y){
744 if(y){//inBuffer already initialized for y=0
745 //erase first line from inBuffer
746 if(dimY>1)
747 inBuffer.erase(inBuffer.begin());
748 //read extra line and push back to inBuffer if not out of bounds
749 if(y+dimY/2<input.nrOfRow()){
750 //allocate buffer
751 if(dimY>1)
752 inBuffer.push_back(inBuffer.back());
753 try{
754 input.readData(inBuffer[inBuffer.size()-1],y+dimY/2);
755 }
756 catch(std::string errorstring){
757 std::cerr << errorstring << "in line " << y << std::endl;
758 }
759 }
760 else{
761 int over=y+dimY/2-input.nrOfRow();
762 int index=(inBuffer.size()-1)-over;
763 assert(index>=0);
764 assert(index<inBuffer.size());
765 inBuffer.push_back(inBuffer[index]);
766 }
767 }
768 if((y+1+down/2)%down)
769 continue;
770 for(int x=0;x<input.nrOfCol();++x){
771 if((x+1+down/2)%down)
772 continue;
773 std::vector<short> potential(m_class.size());
774 for(int iclass=0;iclass<m_class.size();++iclass){
775 potential[iclass]=0;
776 outBuffer[iclass][x/down]=0;
777 }
778 std::vector<double> windowBuffer;
779 int centre=dimX*(dimY-1)/2+(dimX-1)/2;
780 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
781 for(int i=-(dimX-1)/2;i<=dimX/2;++i){
782 if(i!=0&&j!=0&&!eightConnectivity)
783 continue;
784 if(i==0&&j==0)
785 continue;
786 indexI=x+i;
787 //check if out of bounds
788 if(indexI<0)
789 indexI=-indexI;
790 else if(indexI>=input.nrOfCol())
791 indexI=input.nrOfCol()-i;
792 if(y+j<0)
793 indexJ=-j;
794 else if(y+j>=input.nrOfRow())
795 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
796 else
797 indexJ=(dimY-1)/2+j;
798 bool masked=false;
799 for(int imask=0;imask<m_noDataValues.size();++imask){
800 if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
801 masked=true;
802 break;
803 }
804 }
805 if(!masked){
806 for(int iclass=0;iclass<m_class.size();++iclass){
807 if(inBuffer[indexJ][indexI]==m_class[iclass])
808 potential[iclass]+=1;
809 }
810 }
811 }
812 }
813 double norm=0;
814 for(int iclass1=0;iclass1<m_class.size();++iclass1){
815 assert(beta[iclass1].size()==m_class.size());
816 double pot=0;
817 for(int iclass2=0;iclass2<m_class.size();++iclass2)
818 if(iclass2!=iclass1)
819 pot+=potential[iclass2]*beta[iclass1][iclass2];
820 double prior=exp(-pot);
821 outBuffer[iclass1][x/down]=prior;
822 norm+=prior;
823 }
824 if(norm){
825 for(int iclass1=0;iclass1<m_class.size();++iclass1)
826 outBuffer[iclass1][x/down]/=norm;
827 }
828 }
829 progress=(1.0+y/down)/output.nrOfRow();
830 pfnProgress(progress,pszMessage,pProgressArg);
831 //write outBuffer to file
832 assert(outBuffer.size()==m_class.size());
833 assert(y<output.nrOfRow());
834 for(int iclass=0;iclass<m_class.size();++iclass){
835 assert(outBuffer[iclass].size()==output.nrOfCol());
836 try{
837 output.writeData(outBuffer[iclass],y/down,iclass);
838 }
839 catch(std::string errorstring){
840 std::cerr << errorstring << "in class " << iclass << ", line " << y << std::endl;
841 }
842 }
843 }
844}
845
846void filter2d::Filter2d::shift(ImgReaderGdal& input, ImgWriterGdal& output, double offsetX, double offsetY, double randomSigma, RESAMPLE resample, bool verbose)
847{
848 assert(input.nrOfCol()==output.nrOfCol());
849 assert(input.nrOfRow()==output.nrOfRow());
850 assert(input.nrOfBand()==output.nrOfBand());
851 const char* pszMessage;
852 void* pProgressArg=NULL;
853 GDALProgressFunc pfnProgress=GDALTermProgress;
854 double progress=0;
855 pfnProgress(progress,pszMessage,pProgressArg);
856 //process band per band in memory
857 Vector2d<double> inBuffer(input.nrOfRow(),output.nrOfCol());
858 Vector2d<double> outBuffer(input.nrOfRow(),output.nrOfCol());
859 for(int iband=0;iband<input.nrOfBand();++iband){
860 input.readDataBlock(inBuffer,0,inBuffer.nCols()-1,0,inBuffer.nRows()-1,iband);
861 shift(inBuffer,outBuffer,offsetX,offsetY,randomSigma,resample,verbose);
862 output.writeDataBlock(outBuffer,0,outBuffer.nCols()-1,0,outBuffer.nRows()-1,iband);
863 }
864}
865
866//todo: re-implement without dependency of CImg and reg libraries
867// void filter2d::Filter2d::dwt_texture(const std::string& inputFilename, const std::string& outputFilename,int dim, int scale, int down, int iband, bool verbose)
868// {
869// ImgReaderGdal input;
870// ImgWriterGdal output;
871// if(verbose)
872// std::cout << "opening file " << inputFilename << std::endl;
873// input.open(inputFilename);
874// double magicX=1,magicY=1;
875// output.open(outputFilename,(input.nrOfCol()+down-1)/down,(input.nrOfRow()+down-1)/down,scale*3,GDT_Float32,input.getImageType());
876// if(input.isGeoRef()){
877// output.setProjection(input.getProjection());
878// output.copyGeoTransform(input);
879// }
880// if(verbose)
881// std::cout << "Dimension texture (row x col x band) = " << (input.nrOfCol()+down-1)/down << " x " << (input.nrOfRow()+down-1)/down << " x " << scale*3 << std::endl;
882// assert(dim%2);
883// int dimX=dim;
884// int dimY=dim;
885// Vector2d<float> inBuffer(dimY,input.nrOfCol());
886// Vector2d<float> outBuffer(scale*3,(input.nrOfCol()+down-1)/down);
887// //initialize last half of inBuffer
888// int indexI=0;
889// int indexJ=0;
890// for(int j=-dimY/2;j<(dimY+1)/2;++j){
891// try{
892// if(verbose)
893// cout << "reading input line " << abs(j) << std::endl;
894// input.readData(inBuffer[indexJ],GDT_Float32,abs(j),iband);
895// ++indexJ;
896// }
897// catch(std::string errorstring){
898// std::cerr << errorstring << "in band " << iband << ", line " << indexJ << std::endl;
899// }
900// }
901// const char* pszMessage;
902// void* pProgressArg=NULL;
903// GDALProgressFunc pfnProgress=GDALTermProgress;
904// double progress=0;
905// pfnProgress(progress,pszMessage,pProgressArg);
906// for(int y=0;y<input.nrOfRow();y+=down){
907// if(verbose)
908// std::cout << "calculating line " << y/down << std::endl;
909// if(y){//inBuffer already initialized for y=0
910// //erase first line from inBuffer
911// inBuffer.erase(inBuffer.begin());
912// //read extra line and push back to inBuffer if not out of bounds
913// if(y+dimY/2<input.nrOfRow()){
914// //allocate buffer
915// inBuffer.push_back(inBuffer.back());
916// try{
917// if(verbose)
918// std::cout << "reading input line " << y+dimY/2 << std::endl;
919// input.readData(inBuffer[inBuffer.size()-1],GDT_Float32,y+dimY/2,iband);
920// }
921// catch(std::string errorstring){
922// std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
923// }
924// }
925// }
926// for(int x=0;x<input.nrOfCol();x+=down){
927// Vector2d<double> texture_feature(scale,3);
928// CImg<> texture_in(dimX,dimY);
929// int r=0;//index for row of texture_in
930// for(int j=-dimY/2;j<(dimY+1)/2;++j){
931// int c=0;
932// for(int i=-dimX/2;i<(dimX+1)/2;++i){
933// indexI=x+i;
934// //check if out of bounds
935// if(indexI<0)
936// indexI=-indexI;
937// else if(indexI>=input.nrOfCol())
938// indexI=input.nrOfCol()-i;
939// if(y+j<0)
940// indexJ=-j;
941// else if(y+j>=input.nrOfRow())
942// indexJ=dimY/2-j;//indexJ=inBuffer.size()-1-j;
943// else
944// indexJ=dimY/2+j;
945// assert(indexJ<inBuffer.size());
946// assert(indexI<inBuffer[indexJ].size());
947// texture_in(r,c)=inBuffer[indexJ][indexI];
948// c++;
949// }
950// ++r;
951// }
952// texture_in.dwt_texture(texture_feature,scale);
953// for(int v=0;v<scale*3;++v)
954// outBuffer[v][x/down]=texture_feature[v/3][v%3];
955// }
956// //write outBuffer to file
957// try{
958// if(verbose)
959// std::cout << "writing line " << y/down << std::endl;
960// for(int v=0;v<scale*3;++v)
961// output.writeData(outBuffer[v],GDT_Float32,y/down,v);
962// }
963// catch(std::string errorstring){
964// std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
965// }
966// progress=(1.0+y)/output.nrOfRow();
967// pfnProgress(progress,pszMessage,pProgressArg);
968// }
969// input.close();
970// output.close();
971// }
972
973void filter2d::Filter2d::morphology(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dimX, int dimY, const std::vector<double> &angle, bool disc)
974{
975 const char* pszMessage;
976 void* pProgressArg=NULL;
977 GDALProgressFunc pfnProgress=GDALTermProgress;
978 double progress=0;
979 pfnProgress(progress,pszMessage,pProgressArg);
980
981 assert(dimX);
982 assert(dimY);
983
985 for(int iband=0;iband<input.nrOfBand();++iband){
986 Vector2d<double> inBuffer(dimY,input.nrOfCol());
987 std::vector<double> outBuffer(input.nrOfCol());
988 int indexI=0;
989 int indexJ=0;
990 //initialize last half of inBuffer
991 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
992 try{
993 input.readData(inBuffer[indexJ],abs(j),iband);
994 ++indexJ;
995 }
996 catch(std::string errorstring){
997 std::cerr << errorstring << "in line " << indexJ << std::endl;
998 }
999 }
1000 for(int y=0;y<input.nrOfRow();++y){
1001 if(y){//inBuffer already initialized for y=0
1002 //erase first line from inBuffer
1003 if(dimY>1)
1004 inBuffer.erase(inBuffer.begin());
1005 //read extra line and push back to inBuffer if not out of bounds
1006 if(y+dimY/2<input.nrOfRow()){
1007 //allocate buffer
1008 if(dimY>1)
1009 inBuffer.push_back(inBuffer.back());
1010 try{
1011 input.readData(inBuffer[inBuffer.size()-1],y+dimY/2,iband);
1012 }
1013 catch(std::string errorstring){
1014 std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
1015 }
1016 }
1017 else{
1018 int over=y+dimY/2-input.nrOfRow();
1019 int index=(inBuffer.size()-1)-over;
1020 assert(index>=0);
1021 assert(index<inBuffer.size());
1022 inBuffer.push_back(inBuffer[index]);
1023 }
1024 }
1025 for(int x=0;x<input.nrOfCol();++x){
1026 double currentValue=inBuffer[(dimY-1)/2][x];
1027 outBuffer[x]=currentValue;
1028 std::vector<double> statBuffer;
1029 bool currentMasked=false;
1030 int centre=dimX*(dimY-1)/2+(dimX-1)/2;
1031 for(int imask=0;imask<m_noDataValues.size();++imask){
1032 if(currentValue==m_noDataValues[imask]){
1033 currentMasked=true;
1034 break;
1035 }
1036 }
1037 if(currentMasked){
1038 outBuffer[x]=currentValue;
1039 }
1040 else{
1041 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
1042 for(int i=-(dimX-1)/2;i<=dimX/2;++i){
1043 double d2=i*i+j*j;//square distance
1044 if(disc&&(d2>(dimX/2)*(dimY/2)))
1045 continue;
1046 if(angle.size()){
1047 double theta;
1048 //use polar coordinates in radians
1049 if(i>0){
1050 if(j<0)
1051 theta=atan(static_cast<double>(-j)/(static_cast<double>(i)));
1052 else
1053 theta=-atan(static_cast<double>(j)/(static_cast<double>(i)));
1054 }
1055 else if(i<0){
1056 if(j<0)
1057 theta=PI-atan(static_cast<double>(-j)/(static_cast<double>(-i)));
1058 else
1059 theta=PI+atan(static_cast<double>(j)/(static_cast<double>(-i)));
1060 }
1061 else if(j<0)
1062 theta=PI/2.0;
1063 else if(j>0)
1064 theta=3.0*PI/2.0;
1065 //convert to North (0), East (90), South (180), West (270) in degrees
1066 theta=360-(theta/PI*180)+90;
1067 if(theta<0)
1068 theta+=360;
1069 while(theta>360)
1070 theta-=360;
1071 bool alligned=false;
1072 for(int iangle=0;iangle<angle.size();++iangle){
1073 if(sqrt((theta-angle[iangle])*(theta-angle[iangle]))<10){
1074 alligned=true;
1075 break;
1076 }
1077 }
1078 if(!alligned)
1079 continue;
1080 }
1081 indexI=x+i;
1082 //check if out of bounds
1083 if(indexI<0)
1084 indexI=-indexI;
1085 else if(indexI>=input.nrOfCol())
1086 indexI=input.nrOfCol()-i;
1087 if(y+j<0)
1088 indexJ=-j;
1089 else if(y+j>=input.nrOfRow())
1090 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
1091 else
1092 indexJ=(dimY-1)/2+j;
1093 //todo: introduce novalue as this: ?
1094 // if(inBuffer[indexJ][indexI]==(m_noDataValues.size())? m_noDataValues[0] : 0)
1095 // continue;
1096 bool masked=false;
1097 for(int imask=0;imask<m_noDataValues.size();++imask){
1098 if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
1099 masked=true;
1100 break;
1101 }
1102 }
1103 if(!masked){
1104 short binValue=0;
1105 for(int iclass=0;iclass<m_class.size();++iclass){
1106 if(inBuffer[indexJ][indexI]==m_class[iclass]){
1107 binValue=1;
1108 break;
1109 }
1110 }
1111 if(m_class.size())
1112 statBuffer.push_back(binValue);
1113 else
1114 statBuffer.push_back(inBuffer[indexJ][indexI]);
1115 }
1116 }
1117 }
1118 if(statBuffer.size()){
1119 switch(getFilterType(method)){
1120 case(filter2d::dilate):
1121 outBuffer[x]=stat.mymax(statBuffer);
1122 break;
1123 case(filter2d::erode):
1124 outBuffer[x]=stat.mymin(statBuffer);
1125 break;
1126 default:
1127 std::ostringstream ess;
1128 ess << "Error: morphology method " << method << " not supported, choose " << filter2d::dilate << " (dilate) or " << filter2d::erode << " (erode)" << std::endl;
1129 throw(ess.str());
1130 break;
1131 }
1132 }
1133 if(outBuffer[x]&&m_class.size())
1134 outBuffer[x]=m_class[0];
1135 }
1136 }
1137 //write outBuffer to file
1138 try{
1139 output.writeData(outBuffer,y,iband);
1140 }
1141 catch(std::string errorstring){
1142 std::cerr << errorstring << "in band " << iband << ", line " << y << std::endl;
1143 }
1144 progress=(1.0+y);
1145 progress+=(output.nrOfRow()*iband);
1146 progress/=output.nrOfBand()*output.nrOfRow();
1147 pfnProgress(progress,pszMessage,pProgressArg);
1148 }
1149 }
1150}
1151
1152void filter2d::Filter2d::shadowDsm(ImgReaderGdal& input, ImgWriterGdal& output, double sza, double saa, double pixelSize, short shadowFlag){
1153 Vector2d<float> inputBuffer;
1154 Vector2d<float> outputBuffer;
1155 input.readDataBlock(inputBuffer, 0, input.nrOfCol()-1, 0, input.nrOfRow()-1, 0);
1156 shadowDsm(inputBuffer, outputBuffer, sza, saa, pixelSize, shadowFlag);
1157 output.writeDataBlock(outputBuffer,0,output.nrOfCol()-1,0,output.nrOfRow()-1,0);
1158}
1159
1160void filter2d::Filter2d::dwtForward(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family){
1161 Vector2d<float> theBuffer;
1162 for(int iband=0;iband<input.nrOfBand();++iband){
1163 input.readDataBlock(theBuffer, 0, input.nrOfCol()-1, 0, input.nrOfRow()-1, iband);
1164 std::cout << "filtering band " << iband << std::endl << std::flush;
1165 dwtForward(theBuffer, wavelet_type, family);
1166 output.writeDataBlock(theBuffer,0,output.nrOfCol()-1,0,output.nrOfRow()-1,iband);
1167 }
1168}
1169
1170void filter2d::Filter2d::dwtInverse(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family){
1171 Vector2d<float> theBuffer;
1172 for(int iband=0;iband<input.nrOfBand();++iband){
1173 input.readDataBlock(theBuffer, 0, input.nrOfCol()-1, 0, input.nrOfRow()-1, iband);
1174 std::cout << "filtering band " << iband << std::endl << std::flush;
1175 dwtInverse(theBuffer, wavelet_type, family);
1176 output.writeDataBlock(theBuffer,0,output.nrOfCol()-1,0,output.nrOfRow()-1,iband);
1177 }
1178}
1179
1180void filter2d::Filter2d::dwtCut(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family, double cut, bool verbose){
1181 Vector2d<float> theBuffer;
1182 for(int iband=0;iband<input.nrOfBand();++iband){
1183 input.readDataBlock(theBuffer, 0, input.nrOfCol()-1, 0, input.nrOfRow()-1, iband);
1184 std::cout << "filtering band " << iband << std::endl << std::flush;
1185 dwtCut(theBuffer, wavelet_type, family, cut);
1186 output.writeDataBlock(theBuffer,0,output.nrOfCol()-1,0,output.nrOfRow()-1,iband);
1187 }
1188}
1189
1190void filter2d::Filter2d::linearFeature(ImgReaderGdal& input, ImgWriterGdal& output, float angle, float angleStep, float maxDistance, float eps, bool l1, bool a1, bool l2, bool a2, int band, bool verbose){
1191 Vector2d<float> inputBuffer;
1192 std::vector< Vector2d<float> > outputBuffer;
1193 input.readDataBlock(inputBuffer, 0, input.nrOfCol()-1, 0, input.nrOfRow()-1, band);
1194 if(maxDistance<=0)
1195 maxDistance=sqrt(static_cast<float>(input.nrOfCol()*input.nrOfRow()));
1196 linearFeature(inputBuffer,outputBuffer,angle,angleStep,maxDistance,eps, l1, a1, l2, a2,verbose);
1197 for(int iband=0;iband<outputBuffer.size();++iband)
1198 output.writeDataBlock(outputBuffer[iband],0,output.nrOfCol()-1,0,output.nrOfRow()-1,iband);
1199}
1200
1201void filter2d::Filter2d::linearFeature(const Vector2d<float>& input, std::vector< Vector2d<float> >& output, float angle, float angleStep, float maxDistance, float eps, bool l1, bool a1, bool l2, bool a2, bool verbose)
1202{
1203 output.clear();
1204 int nband=0;//linear feature
1205 if(l1)
1206 ++nband;
1207 if(a1)
1208 ++nband;
1209 if(l2)
1210 ++nband;
1211 if(a2)
1212 ++nband;
1213 output.resize(nband);
1214 for(int iband=0;iband<output.size();++iband)
1215 output[iband].resize(input.nRows(),input.nCols());
1216 if(maxDistance<=0)
1217 maxDistance=sqrt(static_cast<float>(input.nRows()*input.nCols()));
1218 int indexI=0;
1219 int indexJ=0;
1220 const char* pszMessage;
1221 void* pProgressArg=NULL;
1222 GDALProgressFunc pfnProgress=GDALTermProgress;
1223 double progress=0;
1224 pfnProgress(progress,pszMessage,pProgressArg);
1225 for(int y=0;y<input.nRows();++y){
1226 for(int x=0;x<input.nCols();++x){
1227 float currentValue=input[y][x];
1228 //find values equal to current value with some error margin
1229 //todo: add distance for two opposite directions
1230 float lineDistance1=0;//longest line of object
1231 float lineDistance2=maxDistance;//shortest line of object
1232 float lineAngle1=0;//angle to longest line (North=0)
1233 float lineAngle2=0;//angle to shortest line (North=0)
1234 float northAngle=0;//rotating angle
1235 for(northAngle=0;northAngle<180;northAngle+=angleStep){
1236 if(angle<=360&&angle>=0&&angle!=northAngle)
1237 continue;
1238 //test
1239 if(verbose)
1240 std::cout << "northAngle: " << northAngle << std::endl;
1241 float currentDistance=0;
1242 float theDir=0;
1243 for(short side=0;side<=1;side+=1){
1244 theDir=PI/2.0-DEG2RAD(northAngle)+side*PI;//in radians
1245 //test
1246 if(verbose)
1247 std::cout << "theDir in deg: " << RAD2DEG(theDir) << std::endl;
1248 if(theDir<0)
1249 theDir+=2*PI;
1250 //test
1251 if(verbose)
1252 std::cout << "theDir in deg: " << RAD2DEG(theDir) << std::endl;
1253 float nextValue=currentValue;
1254 for(float currentRay=1;currentRay<maxDistance;++currentRay){
1255 indexI=x+currentRay*cos(theDir);
1256 indexJ=y-currentRay*sin(theDir);
1257 if(indexJ<0||indexJ>=input.size())
1258 break;
1259 if(indexI<0||indexI>=input[indexJ].size())
1260 break;
1261 nextValue=input[indexJ][indexI];
1262 if(verbose){
1263 std::cout << "x: " << x << std::endl;
1264 std::cout << "y: " << y << std::endl;
1265 std::cout << "currentValue: " << currentValue << std::endl;
1266 std::cout << "theDir in degrees: " << RAD2DEG(theDir) << std::endl;
1267 std::cout << "cos(theDir): " << cos(theDir) << std::endl;
1268 std::cout << "sin(theDir): " << sin(theDir) << std::endl;
1269 std::cout << "currentRay: " << currentRay << std::endl;
1270 std::cout << "currentDistance: " << currentDistance << std::endl;
1271 std::cout << "indexI: " << indexI << std::endl;
1272 std::cout << "indexJ: " << indexJ << std::endl;
1273 std::cout << "nextValue: " << nextValue << std::endl;
1274 }
1275 if(fabs(currentValue-nextValue)<=eps){
1276 ++currentDistance;
1277 //test
1278 if(verbose)
1279 std::cout << "currentDistance: " << currentDistance << ", continue" << std::endl;
1280 }
1281 else{
1282 if(verbose)
1283 std::cout << "currentDistance: " << currentDistance << ", break" << std::endl;
1284 break;
1285 }
1286 }
1287 }
1288 if(lineDistance1<currentDistance){
1289 lineDistance1=currentDistance;
1290 lineAngle1=northAngle;
1291 }
1292 if(lineDistance2>currentDistance){
1293 lineDistance2=currentDistance;
1294 lineAngle2=northAngle;
1295 }
1296 if(verbose){
1297 std::cout << "lineDistance1: " << lineDistance1 << std::endl;
1298 std::cout << "lineAngle1: " << lineAngle1 << std::endl;
1299 std::cout << "lineDistance2: " << lineDistance2 << std::endl;
1300 std::cout << "lineAngle2: " << lineAngle2 << std::endl;
1301 }
1302 }
1303 int iband=0;
1304 if(l1)
1305 output[iband++][y][x]=lineDistance1;
1306 if(a1)
1307 output[iband++][y][x]=lineAngle1;
1308 if(l2)
1309 output[iband++][y][x]=lineDistance2;
1310 if(a2)
1311 output[iband++][y][x]=lineAngle2;
1312 assert(iband==nband);
1313 }
1314 progress=(1.0+y);
1315 progress/=input.nRows();
1316 pfnProgress(progress,pszMessage,pProgressArg);
1317 }
1318}
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
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
void readDataBlock(Vector2d< T > &buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band=0)
Read pixel cell values for a range of columns and rows for a specific band (all indices start countin...
void close(void)
Close the image.
void open(const std::string &filename, const ImgReaderGdal &imgSrc, const std::vector< std::string > &options=std::vector< std::string >())
Open an image for writing, copying image attributes from a source image. Image is directly written to...
bool writeDataBlock(Vector2d< T > &buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band=0)
Write pixel cell values for a range of columns and rows for a specific band (all indices start counti...
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