| pktools 2.6.7 Processing Kernel for geospatial data

Filter2d.cc
1/**********************************************************************
2 Filter2d.cc: class for filtering images
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 <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
37 int 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
44 void filter2d::Filter2d::setTaps(const Vector2d<double> &taps)
45{
46 m_taps=taps;
47}
48
49 void filter2d::Filter2d::smoothNoData(ImgReaderGdal& input, ImgWriterGdal& output, int dim)
50{
51 smoothNoData(input, output,dim,dim);
52}
53
54 void filter2d::Filter2d::smooth(ImgReaderGdal& input, ImgWriterGdal& output, int dim)
55{
56 smooth(input, output,dim,dim);
57}
58
59 void 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
70 void 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
82 void 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
198 void 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
334 void 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
343 void 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
352 void 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
357 void 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
368 statfactory::StatFactory stat;
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
704 void 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]
714 void 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
846 void 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
973 void 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
984 statfactory::StatFactory stat;
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
1152 void 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
1160 void 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
1170 void 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
1180 void 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
1190 void 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
1201 void 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}
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
ImgReaderGdal::close
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
Definition: ImgReaderGdal.cc:42
ImgReaderGdal::open
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
Definition: ImgReaderGdal.cc:36
ImgReaderGdal::readDataBlock
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...
Definition: ImgReaderGdal.h:220
ImgWriterGdal
Definition: ImgWriterGdal.h:38
ImgWriterGdal::close
void close(void)
Close the image.
Definition: ImgWriterGdal.cc:74
ImgWriterGdal::open
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...
Definition: ImgWriterGdal.cc:43
ImgWriterGdal::writeDataBlock
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...
Definition: ImgWriterGdal.h:214
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