| pktools 2.6.7 Processing Kernel for geospatial data

Public Member Functions | Static Public Member Functions | List of all members
filter2d::Filter2d Class Reference

## Public Member Functions
| Filter2d (const Vector2d< double > &taps)
void | setTaps (const Vector2d< double > &taps)
void | pushClass (short theClass=1)
int | pushNoDataValue (double noDataValue=0)
void | pushThreshold (double theThreshold)
void | setThresholds (const std::vector< double > &theThresholds)
void | setClasses (const std::vector< short > &theClasses)
void | filter (ImgReaderGdal &input, ImgWriterGdal &output, bool absolute=false, bool normalize=false, bool noData=false)
void | smooth (ImgReaderGdal &input, ImgWriterGdal &output, int dim)
void | smooth (ImgReaderGdal &input, ImgWriterGdal &output, int dimX, int dimY)
void | smoothNoData (ImgReaderGdal &input, ImgWriterGdal &output, int dim)
void | smoothNoData (ImgReaderGdal &input, ImgWriterGdal &output, int dimX, int dimY)
template<class T1 , class T2 >
void | filter (const Vector2d< T1 > &inputVector, Vector2d< T2 > &outputVector)
template<class T1 , class T2 >
void | smooth (const Vector2d< T1 > &inputVector, Vector2d< T2 > &outputVector, int dim)
template<class T1 , class T2 >
void | smooth (const Vector2d< T1 > &inputVector, Vector2d< T2 > &outputVector, int dimX, int dimY)
void | dwtForward (ImgReaderGdal &input, ImgWriterGdal &output, const std::string &wavelet_type, int family)
void | dwtInverse (ImgReaderGdal &input, ImgWriterGdal &output, const std::string &wavelet_type, int family)
void | dwtCut (ImgReaderGdal &input, ImgWriterGdal &output, const std::string &wavelet_type, int family, double cut, bool verbose=false)
template<class T >
void | dwtForward (Vector2d< T > &data, const std::string &wavelet_type, int family)
template<class T >
void | dwtInverse (Vector2d< T > &data, const std::string &wavelet_type, int family)
template<class T >
void | dwtCut (Vector2d< T > &data, const std::string &wavelet_type, int family, double cut)
void | majorVoting (const std::string &inputFilename, const std::string &outputFilename, int dim=0, const std::vector< int > &prior=std::vector< int >())
void | doit (ImgReaderGdal &input, ImgWriterGdal &output, const std::string &method, int dim, short down=1, bool disc=false)
void | doit (ImgReaderGdal &input, ImgWriterGdal &output, const std::string &method, int dimX, int dimY, short down=1, bool disc=false)
void | mrf (ImgReaderGdal &input, ImgWriterGdal &output, int dimX, int dimY, double beta, bool eightConnectivity=true, short down=1, bool verbose=false)
void | mrf (ImgReaderGdal &input, ImgWriterGdal &output, int dimX, int dimY, Vector2d< double > beta, bool eightConnectivity=true, short down=1, bool verbose=false)
template<class T1 , class T2 >
void | doit (const Vector2d< T1 > &inputVector, Vector2d< T2 > &outputVector, const std::string &method, int dimX, int dimY, short down=1, bool disc=false)
void | median (const std::string &inputFilename, const std::string &outputFilename, int dim, bool disc=false)
void | var (const std::string &inputFilename, const std::string &outputFilename, int dim, bool disc=false)
void | morphology (ImgReaderGdal &input, ImgWriterGdal &output, const std::string &method, int dimX, int dimY, const std::vector< double > &angle, bool disc=false)
template<class T >
unsigned long int | morphology (const Vector2d< T > &input, Vector2d< T > &output, const std::string &method, int dimX, int dimY, bool disc=false, double hThreshold=0)
template<class T >
unsigned long int | dsm2dtm_nwse (const Vector2d< T > &inputDSM, Vector2d< T > &outputMask, double hThreshold, int nlimit, int dim=3)
template<class T >
unsigned long int | dsm2dtm_nesw (const Vector2d< T > &inputDSM, Vector2d< T > &outputMask, double hThreshold, int nlimit, int dim=3)
template<class T >
unsigned long int | dsm2dtm_senw (const Vector2d< T > &inputDSM, Vector2d< T > &outputMask, double hThreshold, int nlimit, int dim=3)
template<class T >
unsigned long int | dsm2dtm_swne (const Vector2d< T > &inputDSM, Vector2d< T > &outputMask, double hThreshold, int nlimit, int dim=3)
template<class T >
void | shadowDsm (const Vector2d< T > &input, Vector2d< T > &output, double sza, double saa, double pixelSize, short shadowFlag=1)
void | shadowDsm (ImgReaderGdal &input, ImgWriterGdal &output, double sza, double saa, double pixelSize, short shadowFlag=1)
void | dwt_texture (const std::string &inputFilename, const std::string &outputFilename, int dim, int scale, int down=1, int iband=0, bool verbose=false)
void | shift (ImgReaderGdal &input, ImgWriterGdal &output, double offsetX=0, double offsetY=0, double randomSigma=0, RESAMPLE resample=BILINEAR, bool verbose=false)
template<class T >
void | shift (const Vector2d< T > &input, Vector2d< T > &output, double offsetX=0, double offsetY=0, double randomSigma=0, RESAMPLE resample=NEAR, bool verbose=false)
void | linearFeature (const Vector2d< float > &input, std::vector< Vector2d< float > > &output, float angle=361, float angleStep=1, float maxDistance=0, float eps=0, bool l1=true, bool a1=true, bool l2=true, bool a2=true, bool verbose=false)
void | linearFeature (ImgReaderGdal &input, ImgWriterGdal &output, float angle=361, float angleStep=1, float maxDistance=0, float eps=0, bool l1=true, bool a1=true, bool l2=true, bool a2=true, int band=0, bool verbose=false)

## Static Public Member Functions
static FILTER_TYPE | getFilterType (const std::string filterType)
static const RESAMPLE | getResampleType (const std::string resampleType)

## Detailed Description

Definition at line 65 of file Filter2d.h.

## Constructor & Destructor Documentation

## ◆ Filter2d() [1/2]

filter2d::Filter2d::Filter2d | ( | void |  | ) |

Definition at line 28 of file Filter2d.cc.

29{
30}

## ◆ Filter2d() [2/2]

filter2d::Filter2d::Filter2d | ( | const Vector2d< double > & | taps | ) |

Definition at line 32 of file Filter2d.cc.

33 : m_taps(taps)
34{
35}

## ◆ ~Filter2d()

virtual filter2d::Filter2d::~Filter2d | ( |  | ) |
| inline virtual

Definition at line 70 of file Filter2d.h.

70{};

## Member Function Documentation

## ◆ doit() [1/3]

template<class T1 , class T2 >

void filter2d::Filter2d::doit | ( | const Vector2d< T1 > & | inputVector,
|  | Vector2d< T2 > & | outputVector,
|  | const std::string & | method,
|  | int | dimX,
|  | int | dimY,
|  | short | down = `1`,
|  | bool | disc = `false`
| ) |  |

Definition at line 259 of file Filter2d.h.

260{
261 const char* pszMessage;
262 void* pProgressArg=NULL;
263 GDALProgressFunc pfnProgress=GDALTermProgress;
264 double progress=0;
265 pfnProgress(progress,pszMessage,pProgressArg);
266
267 statfactory::StatFactory stat;
268 double noDataValue=0;
269 if(m_noDataValues.size()){
270 stat.setNoDataValues(m_noDataValues);
271 noDataValue=m_noDataValues[0];
272 }
273 assert(dimX);
274 assert(dimY);
275
276 outputVector.resize((inputVector.size()+down-1)/down);
277 Vector2d<T1> inBuffer(dimY);
278 std::vector<T2> outBuffer((inputVector[0].size()+down-1)/down);
279
280 int indexI=0;
281 int indexJ=0;
282 //initialize last half of inBuffer
283 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
284 inBuffer[indexJ]=inputVector[abs(j)];
285 ++indexJ;
286 }
287 for(int y=0;y<inputVector.size();++y){
288
289 if(y){//inBuffer already initialized for y=0
290 //erase first line from inBuffer
291 inBuffer.erase(inBuffer.begin());
292 //read extra line and push back to inBuffer if not out of bounds
293 if(y+dimY/2<inputVector.size())
294 inBuffer.push_back(inputVector[y+dimY/2]);
295 else{
296 int over=y+dimY/2-inputVector.size();
297 int index=(inBuffer.size()-1)-over;
298 assert(index>=0);
299 assert(index<inBuffer.size());
300 inBuffer.push_back(inBuffer[index]);
301 }
302 }
303 if((y+1+down/2)%down)
304 continue;
305 for(int x=0;x<inputVector[0].size();++x){
306 if((x+1+down/2)%down)
307 continue;
308 outBuffer[x/down]=0;
309
310 std::vector<T1> windowBuffer;
311 std::map<int,int> occurrence;
312 int centre=dimX*(dimY-1)/2+(dimX-1)/2;
313 bool centreMasked=false;
314 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
315 for(int i=-(dimX-1)/2;i<=dimX/2;++i){
316 indexI=x+i;
317 //check if out of bounds
318 if(indexI<0)
319 indexI=-indexI;
320 else if(indexI>=inputVector[0].size())
321 indexI=inputVector[0].size()-i;
322 if(y+j<0)
323 indexJ=-j;
324 else if(y+j>=inputVector.size())
325 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
326 else
327 indexJ=(dimY-1)/2+j;
328 windowBuffer.push_back(inBuffer[indexJ][indexI]);
329 if(!stat.isNoData(inBuffer[indexJ][indexI])){
330 std::vector<short>::const_iterator vit=m_class.begin();
331 //todo: test if this works (only add occurrence if within defined classes)!
332 if(!m_class.size())
333 ++occurrence[inBuffer[indexJ][indexI]];
334 else{
335 while(vit!=m_class.end()){
336 if(inBuffer[indexJ][indexI]==*(vit++))
337 ++occurrence[inBuffer[indexJ][indexI]];
338 }
339 }
340 }
341 }
342 }
343 switch(getFilterType(method)){
344 case(filter2d::nvalid):
345 outBuffer[x/down]=stat.nvalid(windowBuffer);
346 break;
347 case(filter2d::median):
348 outBuffer[x/down]=stat.median(windowBuffer);
349 break;
350 case(filter2d::var):{
351 outBuffer[x/down]=stat.var(windowBuffer);
352 break;
353 }
354 case(filter2d::stdev):{
355 T2 varValue=stat.var(windowBuffer);
356 if(stat.isNoData(varValue))
357 outBuffer[x/down]=noDataValue;
358 else
359 outBuffer[x/down]=sqrt(varValue);
360 break;
361 }
362 case(filter2d::mean):{
363 if(windowBuffer.empty())
364 outBuffer[x/down]=noDataValue;
365 else
366 outBuffer[x/down]=stat.mean(windowBuffer);
367 break;
368 }
369 case(filter2d::min):{
370 outBuffer[x/down]=stat.mymin(windowBuffer);
371 break;
372 }
373 case(filter2d::ismin):{
374 T1 minValue=stat.mymin(windowBuffer);
375 if(stat.isNoData(minValue))
376 outBuffer[x/down]=noDataValue;
377 else
378 outBuffer[x/down]=(windowBuffer[centre]==minValue)? 1:0;
379 break;
380 }
381 case(filter2d::minmax):{
382 T1 min=0;
383 T1 max=0;
384 stat.minmax(windowBuffer,windowBuffer.begin(),windowBuffer.end(),min,max);
385 if(min!=max)
386 outBuffer[x/down]=0;
387 else
388 outBuffer[x/down]=windowBuffer[centre];//centre pixels
389 break;
390 }
391 case(filter2d::max):{
392 outBuffer[x/down]=stat.mymax(windowBuffer);
393 break;
394 }
395 case(filter2d::ismax):{
396 T1 maxValue=stat.mymax(windowBuffer);
397 if(stat.isNoData(maxValue))
398 outBuffer[x/down]=noDataValue;
399 else
400 outBuffer[x/down]=(windowBuffer[centre]==maxValue)? 1:0;
401 break;
402 }
403 case(filter2d::order):{
404 stat.eraseNoData(windowBuffer);
405 if(windowBuffer.empty())
406 outBuffer[x/down]=noDataValue;
407 else{
408 double lbound=0;
409 double ubound=dimX*dimY;
410 double theMin=stat.mymin(windowBuffer);
411 double theMax=stat.mymax(windowBuffer);
412 double scale=(ubound-lbound)/(theMax-theMin);
413 outBuffer[x/down]=static_cast<short>(scale*(windowBuffer[centre]-theMin)+lbound);
414 }
415 break;
416 }
417 case(filter2d::sum):{
418 outBuffer[x/down]=stat.sum(windowBuffer);
419 break;
420 }
421 case(filter2d::percentile):{
422 assert(m_threshold.size());
423 outBuffer[x/down]=stat.percentile(windowBuffer,windowBuffer.begin(),windowBuffer.end(),m_threshold[0]);
424 break;
425 }
426 case(filter2d::proportion):{
427 stat.eraseNoData(windowBuffer);
428 T2 sum=stat.sum(windowBuffer);
429 if(sum)
430 outBuffer[x/down]=windowBuffer[centre]/sum;
431 else
432 outBuffer[x/down]=noDataValue;
433 break;
434 }
435 case(filter2d::homog):{
436 T1 centreValue=inBuffer[(dimY-1)/2][x];
437 bool isHomog=true;
438 stat.eraseNoData(windowBuffer);
439 typename std::vector<T1>::const_iterator wit;
440 for(wit=windowBuffer.begin();wit!=windowBuffer.end();++wit){
441 if(*wit==centreValue)
442 continue;
443 else{
444 isHomog=false;
445 break;
446 }
447 }
448 if(isHomog)
449 outBuffer[x/down]=1;
450 else
451 outBuffer[x/down]=noDataValue;
452 break;
453 }
454 case(filter2d::heterog):{
455 T1 centreValue=inBuffer[(dimY-1)/2][x];
456 bool isHeterog=true;
457 stat.eraseNoData(windowBuffer);
458 typename std::vector<T1>::const_iterator wit;
459 for(wit=windowBuffer.begin();wit!=windowBuffer.end();++wit){
460 if(*wit!=centreValue)
461 continue;
462 else{
463 isHeterog=false;
464 break;
465 }
466 }
467 if(isHeterog)
468 outBuffer[x/down]=1;
469 else
470 outBuffer[x/down]=noDataValue;
471 break;
472 }
473 case(filter2d::sauvola):{
474 try{
475 double theMean=0;
476 double theStdev=0;
477 bool invalid=false;
478 T1 centreValue=inBuffer[(dimY-1)/2][x];
479 if(windowBuffer.empty()||stat.isNoData(centreValue)){
480 invalid=true;
481 throw(invalid);
482 }
483 stat.meanVar(windowBuffer,theMean,theStdev);
484 theStdev=sqrt(theStdev);
485 double kValue=0.5;
486 double rValue=128;
487 if(m_threshold.size()==2){
488 kValue=m_threshold[0];
489 rValue=m_threshold[1];
490 }
491 //from http://fiji.sc/Auto_Local_Threshold
492 //pixel = ( pixel > mean * ( 1 + k * ( standard_deviation / r - 1 ) ) ) ? object : background
493 double theThreshold=theMean*(1+kValue*(theStdev/rValue - 1));
494 //isdata value hardcoded as 1 for now
495 outBuffer[x/down]=(centreValue>theThreshold) ? 1 : noDataValue;
496 }
497 catch(bool invalid){
498 outBuffer[x/down]=noDataValue;
499 }
500 break;
501 }
502 case(filter2d::density):{
503 int nvalid=stat.nvalid(windowBuffer);
504 if(nvalid){
505 std::vector<short>::const_iterator vit=m_class.begin();
506 while(vit!=m_class.end())
507 outBuffer[x/down]+=100.0*occurrence[*(vit++)]/nvalid;
508 }
509 else
510 outBuffer[x/down]=noDataValue;
511 break;
512 }
513 case(filter2d::countid):{
514 if(occurrence.size())
515 outBuffer[x/down]=occurrence.size();
516 else
517 outBuffer[x/down]=noDataValue;
518 break;
519 }
520 case(filter2d::mode):{
521 if(occurrence.size()){
522 std::map<int,int>::const_iterator maxit=occurrence.begin();
523 for(std::map<int,int>::const_iterator mit=occurrence.begin();mit!=occurrence.end();++mit){
524 if(mit->second>maxit->second)
525 maxit=mit;
526 }
527 if(occurrence[inBuffer[(dimY-1)/2][x]]<maxit->second)//
528 outBuffer[x/down]=maxit->first;
529 else//favorize original value in case of ties
530 outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
531 }
532 else
533 outBuffer[x/down]=noDataValue;
534 break;
535 }
536 case(filter2d::threshold):{
537 assert(m_class.size()==m_threshold.size());
538 int nvalid=stat.nvalid(windowBuffer);
539 if(nvalid>0){
540 outBuffer[x/down]=inBuffer[(dimY-1)/2][x];//initialize with original value (in case thresholds not met)
541 for(int iclass=0;iclass<m_class.size();++iclass){
542 if(100.0*(occurrence[m_class[iclass]])/nvalid>m_threshold[iclass])
543 outBuffer[x/down]=m_class[iclass];
544 }
545 }
546 else
547 outBuffer[x/down]=noDataValue;
548 break;
549 }
550 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...
551 if(windowBuffer.size()){
552 int randomIndex=std::rand()%windowBuffer.size();
553 if(randomIndex>=windowBuffer.size())
554 outBuffer[x/down]=windowBuffer.back();
555 else if(randomIndex<0)
556 outBuffer[x/down]=windowBuffer[0];
557 else
558 outBuffer[x/down]=windowBuffer[randomIndex];
559 }
560 else
561 outBuffer[x/down]=noDataValue;
562 break;
563 }
564 case(filter2d::mixed):{
565 enum MixType { BF=11, CF=12, MF=13, NF=20, W=30 };
566 double nBF=occurrence[BF];
567 double nCF=occurrence[CF];
568 double nMF=occurrence[MF];
569 double nNF=occurrence[NF];
570 double nW=occurrence[W];
571 if(windowBuffer.size()){
572 if((nBF+nCF+nMF)&&(nBF+nCF+nMF>=nNF+nW)){//forest
573 if(nBF/(nBF+nCF)>=0.75)
574 outBuffer[x/down]=BF;
575 else if(nCF/(nBF+nCF)>=0.75)
576 outBuffer[x/down]=CF;
577 else
578 outBuffer[x/down]=MF;
579 }
580 else{//non-forest
581 if(nW&&(nW>=nNF))
582 outBuffer[x/down]=W;
583 else
584 outBuffer[x/down]=NF;
585 }
586 }
587 else
588 outBuffer[x/down]=inBuffer[indexJ][indexI];
589 break;
590 }
591 default:
592 break;
593 }
594 }
595 progress=(1.0+y/down);
596 progress/=outputVector.size();
597 pfnProgress(progress,pszMessage,pProgressArg);
598 //copy outBuffer to outputVector
599 outputVector[y/down]=outBuffer;
600 }
601}
Vector2d
Definition: Vector2d.h:32
statfactory::StatFactory
Definition: StatFactory.h:43

## ◆ doit() [2/3]

void filter2d::Filter2d::doit | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | const std::string & | method,
|  | int | dim,
|  | short | down = `1`,
|  | bool | disc = `false`
| ) |  |

Definition at line 352 of file Filter2d.cc.

353{
354 doit(input,output,method,dim,dim,down,disc);
355}

## ◆ doit() [3/3]

void filter2d::Filter2d::doit | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | const std::string & | method,
|  | int | dimX,
|  | int | dimY,
|  | short | down = `1`,
|  | bool | disc = `false`
| ) |  |

Definition at line 357 of file Filter2d.cc.

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}
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::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::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

## ◆ dsm2dtm_nesw()

template<class T >

unsigned long int filter2d::Filter2d::dsm2dtm_nesw | ( | const Vector2d< T > & | inputDSM,
|  | Vector2d< T > & | outputMask,
|  | double | hThreshold,
|  | int | nlimit,
|  | int | dim = `3`
| ) |  |

Definition at line 928 of file Filter2d.h.

929{
930 const char* pszMessage;
931 void* pProgressArg=NULL;
932 GDALProgressFunc pfnProgress=GDALTermProgress;
933 double progress=0;
934 pfnProgress(progress,pszMessage,pProgressArg);
935
936 Vector2d<T> tmpDSM(inputDSM);
937 double noDataValue=0;
938 if(m_noDataValues.size())
939 noDataValue=m_noDataValues[0];
940
941 unsigned long int nchange=0;
942 int dimX=dim;
943 int dimY=dim;
944 assert(dimX);
945 assert(dimY);
946 statfactory::StatFactory stat;
947 Vector2d<T> inBuffer(dimY,inputDSM.nCols());
948 if(outputMask.size()!=inputDSM.nRows())
949 outputMask.resize(inputDSM.nRows());
950 int indexI=0;
951 int indexJ=0;
952 //initialize last half of inBuffer
953 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
954 for(int i=0;i<inputDSM.nCols();++i)
955 inBuffer[indexJ][i]=tmpDSM[abs(j)][i];
956 ++indexJ;
957 }
958 for(int y=0;y<tmpDSM.nRows();++y){
959 if(y){//inBuffer already initialized for y=0
960 //erase first line from inBuffer
961 inBuffer.erase(inBuffer.begin());
962 //read extra line and push back to inBuffer if not out of bounds
963 if(y+dimY/2<tmpDSM.nRows()){
964 //allocate buffer
965 inBuffer.push_back(inBuffer.back());
966 for(int i=0;i<tmpDSM.nCols();++i)
967 inBuffer[inBuffer.size()-1][i]=tmpDSM[y+dimY/2][i];
968 }
969 else{
970 int over=y+dimY/2-tmpDSM.nRows();
971 int index=(inBuffer.size()-1)-over;
972 assert(index>=0);
973 assert(index<inBuffer.size());
974 inBuffer.push_back(inBuffer[index]);
975 }
976 }
977 for(int x=tmpDSM.nCols()-1;x>=0;--x){
978 double centerValue=inBuffer[(dimY-1)/2][x];
979 short nmasked=0;
980 std::vector<T> neighbors;
981 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
982 for(int i=-(dimX-1)/2;i<=dimX/2;++i){
983 indexI=x+i;
984 //check if out of bounds
985 if(indexI<0)
986 indexI=-indexI;
987 else if(indexI>=tmpDSM.nCols())
988 indexI=tmpDSM.nCols()-i;
989 if(y+j<0)
990 indexJ=-j;
991 else if(y+j>=tmpDSM.nRows())
992 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
993 else
994 indexJ=(dimY-1)/2+j;
995 double difference=(centerValue-inBuffer[indexJ][indexI]);
996 if(i||j)//skip centerValue
997 neighbors.push_back(inBuffer[indexJ][indexI]);
998 if(difference>hThreshold)
999 ++nmasked;
1000 }
1001 }
1002 if(nmasked<=nlimit){
1003 ++nchange;
1004 //reset pixel in outputMask
1005 outputMask[y][x]=0;
1006 }
1007 else{
1008 //reset pixel height in tmpDSM
1009 sort(neighbors.begin(),neighbors.end());
1010 assert(neighbors.size()>1);
1011 inBuffer[(dimY-1)/2][x]=neighbors[1];
1012 /* inBuffer[(dimY-1)/2][x]=stat.mymin(neighbors); */
1013 }
1014 }
1015 progress=(1.0+y);
1016 progress/=outputMask.nRows();
1017 pfnProgress(progress,pszMessage,pProgressArg);
1018 }
1019 return nchange;
1020}

## ◆ dsm2dtm_nwse()

template<class T >

unsigned long int filter2d::Filter2d::dsm2dtm_nwse | ( | const Vector2d< T > & | inputDSM,
|  | Vector2d< T > & | outputMask,
|  | double | hThreshold,
|  | int | nlimit,
|  | int | dim = `3`
| ) |  |

Definition at line 834 of file Filter2d.h.

835{
836 const char* pszMessage;
837 void* pProgressArg=NULL;
838 GDALProgressFunc pfnProgress=GDALTermProgress;
839 double progress=0;
840 pfnProgress(progress,pszMessage,pProgressArg);
841
842 Vector2d<T> tmpDSM(inputDSM);
843 double noDataValue=0;
844 if(m_noDataValues.size())
845 noDataValue=m_noDataValues[0];
846
847 unsigned long int nchange=0;
848 int dimX=dim;
849 int dimY=dim;
850 assert(dimX);
851 assert(dimY);
852 statfactory::StatFactory stat;
853 Vector2d<T> inBuffer(dimY,inputDSM.nCols());
854 if(outputMask.size()!=inputDSM.nRows())
855 outputMask.resize(inputDSM.nRows());
856 int indexI=0;
857 int indexJ=0;
858 //initialize last half of inBuffer
859 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
860 for(int i=0;i<inputDSM.nCols();++i)
861 inBuffer[indexJ][i]=tmpDSM[abs(j)][i];
862 ++indexJ;
863 }
864 for(int y=0;y<tmpDSM.nRows();++y){
865 if(y){//inBuffer already initialized for y=0
866 //erase first line from inBuffer
867 inBuffer.erase(inBuffer.begin());
868 //read extra line and push back to inBuffer if not out of bounds
869 if(y+dimY/2<tmpDSM.nRows()){
870 //allocate buffer
871 inBuffer.push_back(inBuffer.back());
872 for(int i=0;i<tmpDSM.nCols();++i)
873 inBuffer[inBuffer.size()-1][i]=tmpDSM[y+dimY/2][i];
874 }
875 else{
876 int over=y+dimY/2-tmpDSM.nRows();
877 int index=(inBuffer.size()-1)-over;
878 assert(index>=0);
879 assert(index<inBuffer.size());
880 inBuffer.push_back(inBuffer[index]);
881 }
882 }
883 for(int x=0;x<tmpDSM.nCols();++x){
884 double centerValue=inBuffer[(dimY-1)/2][x];
885 short nmasked=0;
886 std::vector<T> neighbors;
887 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
888 for(int i=-(dimX-1)/2;i<=dimX/2;++i){
889 indexI=x+i;
890 //check if out of bounds
891 if(indexI<0)
892 indexI=-indexI;
893 else if(indexI>=tmpDSM.nCols())
894 indexI=tmpDSM.nCols()-i;
895 if(y+j<0)
896 indexJ=-j;
897 else if(y+j>=tmpDSM.nRows())
898 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
899 else
900 indexJ=(dimY-1)/2+j;
901 double difference=(centerValue-inBuffer[indexJ][indexI]);
902 if(i||j)//skip centerValue
903 neighbors.push_back(inBuffer[indexJ][indexI]);
904 if(difference>hThreshold)
905 ++nmasked;
906 }
907 }
908 if(nmasked<=nlimit){
909 ++nchange;
910 //reset pixel in outputMask
911 outputMask[y][x]=0;
912 }
913 else{
914 //reset pixel height in tmpDSM
915 sort(neighbors.begin(),neighbors.end());
916 assert(neighbors.size()>1);
917 inBuffer[(dimY-1)/2][x]=neighbors[1];
918 /* inBuffer[(dimY-1)/2][x]=stat.mymin(neighbors); */
919 }
920 }
921 progress=(1.0+y);
922 progress/=outputMask.nRows();
923 pfnProgress(progress,pszMessage,pProgressArg);
924 }
925 return nchange;
926}

## ◆ dsm2dtm_senw()

template<class T >

unsigned long int filter2d::Filter2d::dsm2dtm_senw | ( | const Vector2d< T > & | inputDSM,
|  | Vector2d< T > & | outputMask,
|  | double | hThreshold,
|  | int | nlimit,
|  | int | dim = `3`
| ) |  |

Definition at line 1022 of file Filter2d.h.

1023{
1024 const char* pszMessage;
1025 void* pProgressArg=NULL;
1026 GDALProgressFunc pfnProgress=GDALTermProgress;
1027 double progress=0;
1028 pfnProgress(progress,pszMessage,pProgressArg);
1029
1030 Vector2d<T> tmpDSM(inputDSM);
1031 double noDataValue=0;
1032 if(m_noDataValues.size())
1033 noDataValue=m_noDataValues[0];
1034
1035 unsigned long int nchange=0;
1036 int dimX=dim;
1037 int dimY=dim;
1038 assert(dimX);
1039 assert(dimY);
1040 statfactory::StatFactory stat;
1041 Vector2d<T> inBuffer(dimY,inputDSM.nCols());
1042 if(outputMask.size()!=inputDSM.nRows())
1043 outputMask.resize(inputDSM.nRows());
1044 int indexI=0;
1045 int indexJ=inputDSM.nRows()-1;
1046 //initialize first half of inBuffer
1047 for(int j=inputDSM.nRows()-dimY/2;j<inputDSM.nRows();--j){
1048 for(int i=0;i<inputDSM.nCols();++i)
1049 inBuffer[indexJ][i]=tmpDSM[abs(j)][i];
1050 ++indexJ;
1051 }
1052 for(int y=tmpDSM.nRows()-1;y>=0;--y){
1053 if(y<tmpDSM.nRows()-1){//inBuffer already initialized for y=tmpDSM.nRows()-1
1054 //erase last line from inBuffer
1055 inBuffer.erase(inBuffer.end()-1);
1056 //read extra line and insert to inBuffer if not out of bounds
1057 if(y-dimY/2>0){
1058 //allocate buffer
1059 inBuffer.insert(inBuffer.begin(),inBuffer.back());
1060 for(int i=0;i<tmpDSM.nCols();++i)
1061 inBuffer[0][i]=tmpDSM[y-dimY/2][i];
1062 }
1063 else{
1064 inBuffer.insert(inBuffer.begin(),inBuffer[abs(y-dimY/2)]);
1065 }
1066 }
1067 for(int x=tmpDSM.nCols()-1;x>=0;--x){
1068 double centerValue=inBuffer[(dimY-1)/2][x];
1069 short nmasked=0;
1070 std::vector<T> neighbors;
1071 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
1072 for(int i=-(dimX-1)/2;i<=dimX/2;++i){
1073 indexI=x+i;
1074 //check if out of bounds
1075 if(indexI<0)
1076 indexI=-indexI;
1077 else if(indexI>=tmpDSM.nCols())
1078 indexI=tmpDSM.nCols()-i;
1079 if(y+j<0)
1080 indexJ=-j;
1081 else if(y+j>=tmpDSM.nRows())
1082 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
1083 else
1084 indexJ=(dimY-1)/2+j;
1085 double difference=(centerValue-inBuffer[indexJ][indexI]);
1086 if(i||j)//skip centerValue
1087 neighbors.push_back(inBuffer[indexJ][indexI]);
1088 if(difference>hThreshold)
1089 ++nmasked;
1090 }
1091 }
1092 if(nmasked<=nlimit){
1093 ++nchange;
1094 //reset pixel in outputMask
1095 outputMask[y][x]=0;
1096 }
1097 else{
1098 //reset pixel height in tmpDSM
1099 sort(neighbors.begin(),neighbors.end());
1100 assert(neighbors.size()>1);
1101 inBuffer[(dimY-1)/2][x]=neighbors[1];
1102 /* inBuffer[(dimY-1)/2][x]=stat.mymin(neighbors); */
1103 }
1104 }
1105 progress=(1.0+y);
1106 progress/=outputMask.nRows();
1107 pfnProgress(progress,pszMessage,pProgressArg);
1108 }
1109 return nchange;
1110}

## ◆ dsm2dtm_swne()

template<class T >

unsigned long int filter2d::Filter2d::dsm2dtm_swne | ( | const Vector2d< T > & | inputDSM,
|  | Vector2d< T > & | outputMask,
|  | double | hThreshold,
|  | int | nlimit,
|  | int | dim = `3`
| ) |  |

Definition at line 1112 of file Filter2d.h.

1113{
1114 const char* pszMessage;
1115 void* pProgressArg=NULL;
1116 GDALProgressFunc pfnProgress=GDALTermProgress;
1117 double progress=0;
1118 pfnProgress(progress,pszMessage,pProgressArg);
1119
1120 Vector2d<T> tmpDSM(inputDSM);
1121 double noDataValue=0;
1122 if(m_noDataValues.size())
1123 noDataValue=m_noDataValues[0];
1124
1125 unsigned long int nchange=0;
1126 int dimX=dim;
1127 int dimY=dim;
1128 assert(dimX);
1129 assert(dimY);
1130 statfactory::StatFactory stat;
1131 Vector2d<T> inBuffer(dimY,inputDSM.nCols());
1132 if(outputMask.size()!=inputDSM.nRows())
1133 outputMask.resize(inputDSM.nRows());
1134 int indexI=0;
1135 int indexJ=0;
1136 //initialize first half of inBuffer
1137 for(int j=inputDSM.nRows()-dimY/2;j<inputDSM.nRows();--j){
1138 for(int i=0;i<inputDSM.nCols();++i)
1139 inBuffer[indexJ][i]=tmpDSM[abs(j)][i];
1140 ++indexJ;
1141 }
1142 for(int y=tmpDSM.nRows()-1;y>=0;--y){
1143 if(y<tmpDSM.nRows()-1){//inBuffer already initialized for y=0
1144 //erase last line from inBuffer
1145 inBuffer.erase(inBuffer.end()-1);
1146 //read extra line and insert to inBuffer if not out of bounds
1147 if(y-dimY/2>0){
1148 //allocate buffer
1149 inBuffer.insert(inBuffer.begin(),inBuffer.back());
1150 for(int i=0;i<tmpDSM.nCols();++i)
1151 inBuffer[0][i]=tmpDSM[y-dimY/2][i];
1152 }
1153 else{
1154 inBuffer.insert(inBuffer.begin(),inBuffer[abs(y-dimY/2)]);
1155 }
1156 }
1157 for(int x=0;x<tmpDSM.nCols();++x){
1158 double centerValue=inBuffer[(dimY-1)/2][x];
1159 short nmasked=0;
1160 std::vector<T> neighbors;
1161 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
1162 for(int i=-(dimX-1)/2;i<=dimX/2;++i){
1163 indexI=x+i;
1164 //check if out of bounds
1165 if(indexI<0)
1166 indexI=-indexI;
1167 else if(indexI>=tmpDSM.nCols())
1168 indexI=tmpDSM.nCols()-i;
1169 if(y+j<0)
1170 indexJ=-j;
1171 else if(y+j>=tmpDSM.nRows())
1172 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
1173 else
1174 indexJ=(dimY-1)/2+j;
1175 double difference=(centerValue-inBuffer[indexJ][indexI]);
1176 if(i||j)//skip centerValue
1177 neighbors.push_back(inBuffer[indexJ][indexI]);
1178 if(difference>hThreshold)
1179 ++nmasked;
1180 }
1181 }
1182 if(nmasked<=nlimit){
1183 ++nchange;
1184 //reset pixel in outputMask
1185 outputMask[y][x]=0;
1186 }
1187 else{
1188 //reset pixel height in tmpDSM
1189 sort(neighbors.begin(),neighbors.end());
1190 assert(neighbors.size()>1);
1191 inBuffer[(dimY-1)/2][x]=neighbors[1];
1192 /* inBuffer[(dimY-1)/2][x]=stat.mymin(neighbors); */
1193 }
1194 }
1195 progress=(1.0+y);
1196 progress/=outputMask.nRows();
1197 pfnProgress(progress,pszMessage,pProgressArg);
1198 }
1199 return nchange;
1200}

## ◆ dwtCut() [1/2]

void filter2d::Filter2d::dwtCut | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | const std::string & | wavelet_type,
|  | int | family,
|  | double | cut,
|  | bool | verbose = `false`
| ) |  |

Definition at line 1180 of file Filter2d.cc.

1180 {
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}
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::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

## ◆ dwtCut() [2/2]

template<class T >

void filter2d::Filter2d::dwtCut | ( | Vector2d< T > & | data,
|  | const std::string & | wavelet_type,
|  | int | family,
|  | double | cut
| ) |  |

Definition at line 1338 of file Filter2d.h.

1338 {
1339 const char* pszMessage;
1340 void* pProgressArg=NULL;
1341 GDALProgressFunc pfnProgress=GDALTermProgress;
1342 double progress=0;
1343 pfnProgress(progress,pszMessage,pProgressArg);
1344
1345 int nRow=theBuffer.size();
1346 assert(nRow);
1347 int nCol=theBuffer[0].size();
1348 assert(nCol);
1349 //make sure data size if power of 2
1350 while(theBuffer.size()&(theBuffer.size()-1))
1351 theBuffer.push_back(theBuffer.back());
1352 for(int irow=0;irow<theBuffer.size();++irow)
1353 while(theBuffer[irow].size()&(theBuffer[irow].size()-1))
1354 theBuffer[irow].push_back(theBuffer[irow].back());
1355 double* data=new double[theBuffer.size()*theBuffer[0].size()];
1356 double* abscoeff=new double[theBuffer.size()*theBuffer[0].size()];
1357 size_t* p=new size_t[theBuffer.size()*theBuffer[0].size()];
1358 for(int irow=0;irow<theBuffer.size();++irow){
1359 for(int icol=0;icol<theBuffer[0].size();++icol){
1360 int index=irow*theBuffer[0].size()+icol;
1361 assert(index<theBuffer.size()*theBuffer[0].size());
1362 data[index]=theBuffer[irow][icol];
1363 }
1364 }
1365 int nsize=theBuffer.size()*theBuffer[0].size();
1366 gsl_wavelet *w;
1367 gsl_wavelet_workspace *work;
1368 assert(nsize);
1369 w=gsl_wavelet_alloc(filter::Filter::getWaveletType(wavelet_type),family);
1370 work=gsl_wavelet_workspace_alloc(nsize);
1371 gsl_wavelet2d_nstransform_forward (w, data, theBuffer.size(), theBuffer[0].size(),theBuffer[0].size(), work);
1372 for(int irow=0;irow<theBuffer.size();++irow){
1373 for(int icol=0;icol<theBuffer[0].size();++icol){
1374 int index=irow*theBuffer[0].size()+icol;
1375 abscoeff[index]=fabs(data[index]);
1376 }
1377 }
1378 int nc=(100-cut)/100.0*nsize;
1379 gsl_sort_index(p,abscoeff,1,nsize);
1380 for(int i=0;(i+nc)<nsize;i++)
1381 data[p[i]]=0;
1382 gsl_wavelet2d_nstransform_inverse (w, data, theBuffer.size(), theBuffer[0].size(),theBuffer[0].size(), work);
1383 for(int irow=0;irow<theBuffer.size();++irow){
1384 for(int icol=0;icol<theBuffer[irow].size();++icol){
1385 int index=irow*theBuffer[irow].size()+icol;
1386 theBuffer[irow][icol]=data[index];
1387 }
1388 progress=(1.0+irow);
1389 progress/=theBuffer.nRows();
1390 pfnProgress(progress,pszMessage,pProgressArg);
1391 }
1392 theBuffer.erase(theBuffer.begin()+nRow,theBuffer.end());
1393 for(int irow=0;irow<theBuffer.size();++irow)
1394 theBuffer[irow].erase(theBuffer[irow].begin()+nCol,theBuffer[irow].end());
1395 delete[] data;
1396 delete[] abscoeff;
1397 delete[] p;
1398 gsl_wavelet_free (w);
1399 gsl_wavelet_workspace_free (work);
1400
1401}

## ◆ dwtForward() [1/2]

void filter2d::Filter2d::dwtForward | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | const std::string & | wavelet_type,
|  | int | family
| ) |  |

Definition at line 1160 of file Filter2d.cc.

1160 {
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}

## ◆ dwtForward() [2/2]

template<class T >

void filter2d::Filter2d::dwtForward | ( | Vector2d< T > & | data,
|  | const std::string & | wavelet_type,
|  | int | family
| ) |  |

Definition at line 1243 of file Filter2d.h.

1243 {
1244 const char* pszMessage;
1245 void* pProgressArg=NULL;
1246 GDALProgressFunc pfnProgress=GDALTermProgress;
1247 double progress=0;
1248 pfnProgress(progress,pszMessage,pProgressArg);
1249
1250 int nRow=theBuffer.size();
1251 assert(nRow);
1252 int nCol=theBuffer[0].size();
1253 assert(nCol);
1254 //make sure data size if power of 2
1255 while(theBuffer.size()&(theBuffer.size()-1))
1256 theBuffer.push_back(theBuffer.back());
1257 for(int irow=0;irow<theBuffer.size();++irow)
1258 while(theBuffer[irow].size()&(theBuffer[irow].size()-1))
1259 theBuffer[irow].push_back(theBuffer[irow].back());
1260 std::vector<double> vdata(theBuffer.size()*theBuffer[0].size());
1261 double* data=&(vdata[0]);
1262 for(int irow=0;irow<theBuffer.size();++irow){
1263 for(int icol=0;icol<theBuffer[0].size();++icol){
1264 int index=irow*theBuffer[0].size()+icol;
1265 data[index]=theBuffer[irow][icol];
1266 }
1267 }
1268 int nsize=theBuffer.size()*theBuffer[0].size();
1269 gsl_wavelet *w;
1270 gsl_wavelet_workspace *work;
1271 assert(nsize);
1272 w=gsl_wavelet_alloc(filter::Filter::getWaveletType(wavelet_type),family);
1273 work=gsl_wavelet_workspace_alloc(nsize);
1274 gsl_wavelet2d_nstransform_forward (w, data, theBuffer.size(), theBuffer.size(),theBuffer[0].size(), work);
1275 theBuffer.erase(theBuffer.begin()+nRow,theBuffer.end());
1276 for(int irow=0;irow<theBuffer.size();++irow){
1277 theBuffer[irow].erase(theBuffer[irow].begin()+nCol,theBuffer[irow].end());
1278 for(int icol=0;icol<theBuffer[irow].size();++icol){
1279 int index=irow*theBuffer[irow].size()+icol;
1280 theBuffer[irow][icol]=data[index];
1281 }
1282 progress=(1.0+irow);
1283 progress/=theBuffer.nRows();
1284 pfnProgress(progress,pszMessage,pProgressArg);
1285 }
1286 gsl_wavelet_free (w);
1287 gsl_wavelet_workspace_free (work);
1288}

## ◆ dwtInverse() [1/2]

void filter2d::Filter2d::dwtInverse | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | const std::string & | wavelet_type,
|  | int | family
| ) |  |

Definition at line 1170 of file Filter2d.cc.

1170 {
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}

## ◆ dwtInverse() [2/2]

template<class T >

void filter2d::Filter2d::dwtInverse | ( | Vector2d< T > & | data,
|  | const std::string & | wavelet_type,
|  | int | family
| ) |  |

Definition at line 1290 of file Filter2d.h.

1290 {
1291 const char* pszMessage;
1292 void* pProgressArg=NULL;
1293 GDALProgressFunc pfnProgress=GDALTermProgress;
1294 double progress=0;
1295 pfnProgress(progress,pszMessage,pProgressArg);
1296
1297 int nRow=theBuffer.size();
1298 assert(nRow);
1299 int nCol=theBuffer[0].size();
1300 assert(nCol);
1301 //make sure data size if power of 2
1302 while(theBuffer.size()&(theBuffer.size()-1))
1303 theBuffer.push_back(theBuffer.back());
1304 for(int irow=0;irow<theBuffer.size();++irow)
1305 while(theBuffer[irow].size()&(theBuffer[irow].size()-1))
1306 theBuffer[irow].push_back(theBuffer[irow].back());
1307 std::vector<double> vdata(theBuffer.size()*theBuffer[0].size());
1308 double* data=&(vdata[0]);
1309 //double data[theBuffer.size()*theBuffer[0].size()];
1310 for(int irow=0;irow<theBuffer.size();++irow){
1311 for(int icol=0;icol<theBuffer[0].size();++icol){
1312 int index=irow*theBuffer[0].size()+icol;
1313 data[index]=theBuffer[irow][icol];
1314 }
1315 }
1316 int nsize=theBuffer.size()*theBuffer[0].size();
1317 gsl_wavelet *w;
1318 gsl_wavelet_workspace *work;
1319 assert(nsize);
1320 w=gsl_wavelet_alloc(filter::Filter::getWaveletType(wavelet_type),family);
1321 work=gsl_wavelet_workspace_alloc(nsize);
1322 gsl_wavelet2d_nstransform_inverse (w, data, theBuffer.size(), theBuffer.size(),theBuffer[0].size(), work);
1323 theBuffer.erase(theBuffer.begin()+nRow,theBuffer.end());
1324 for(int irow=0;irow<theBuffer.size();++irow){
1325 theBuffer[irow].erase(theBuffer[irow].begin()+nCol,theBuffer[irow].end());
1326 for(int icol=0;icol<theBuffer[irow].size();++icol){
1327 int index=irow*theBuffer[irow].size()+icol;
1328 theBuffer[irow][icol]=data[index];
1329 }
1330 progress=(1.0+irow);
1331 progress/=theBuffer.nRows();
1332 pfnProgress(progress,pszMessage,pProgressArg);
1333 }
1334 gsl_wavelet_free (w);
1335 gsl_wavelet_workspace_free (work);
1336}

## ◆ filter() [1/2]

template<class T1 , class T2 >

void filter2d::Filter2d::filter | ( | const Vector2d< T1 > & | inputVector,
|  | Vector2d< T2 > & | outputVector
| ) |  |

Definition at line 202 of file Filter2d.h.

203 {
204 outputVector.resize(inputVector.size());
205 int dimX=m_taps[0].size();//horizontal!!!
206 int dimY=m_taps.size();//vertical!!!
207 Vector2d<T1> inBuffer(dimY);
208 std::vector<T2> outBuffer(inputVector[0].size());
209 //initialize last half of inBuffer
210 int indexI=0;
211 int indexJ=0;
212 //initialize last half of inBuffer
213 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
214 inBuffer[indexJ]=inputVector[abs(j)];
215 ++indexJ;
216 }
217
218 for(int y=0;y<inputVector.size();++y){
219 if(y){//inBuffer already initialized for y=0
220 //erase first line from inBuffer
221 inBuffer.erase(inBuffer.begin());
222 //read extra line and push back to inBuffer if not out of bounds
223 if(y+dimY/2<inputVector.size()){
224 //allocate buffer
225 inBuffer.push_back(inputVector[y+dimY/2]);
226 }
227 else{
228 int over=y+dimY/2-inputVector.nRows();
229 int index=(inBuffer.size()-1)-over;
230 assert(index>=0);
231 assert(index<inBuffer.size());
232 inBuffer.push_back(inBuffer[index]);
233 }
234 }
235 for(int x=0;x<inputVector.nCols();++x){
236 outBuffer[x]=0;
237 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
238 for(int i=-(dimX-1)/2;i<=dimX/2;++i){
239 indexI=x+i;
240 indexJ=(dimY-1)/2+j;
241 //check if out of bounds
242 if(x<(dimX-1)/2)
243 indexI=x+abs(i);
244 else if(x>=inputVector.nCols()-(dimX-1)/2)
245 indexI=x-abs(i);
246 if(y<(dimY-1)/2)
247 indexJ=(dimY-1)/2+abs(j);
248 else if(y>=inputVector.nRows()-(dimY-1)/2)
249 indexJ=(dimY-1)/2-abs(j);
250 outBuffer[x]+=(m_taps[(dimY-1)/2+j][(dimX-1)/2+i]*inBuffer[indexJ][indexI]);
251 }
252 }
253 }
254 //copy outBuffer to outputVector
255 outputVector[y]=outBuffer;
256 }
257 }

## ◆ filter() [2/2]

void filter2d::Filter2d::filter | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | bool | absolute = `false`,
|  | bool | normalize = `false`,
|  | bool | noData = `false`
| ) |  |

Definition at line 82 of file Filter2d.cc.

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}

## ◆ getFilterType()

static FILTER_TYPE filter2d::Filter2d::getFilterType | ( | const std::string | filterType | ) |
| inline static

Definition at line 71 of file Filter2d.h.

71 {
72 std::map<std::string, FILTER_TYPE> m_filterMap;
73 initMap(m_filterMap);
74 return m_filterMap[filterType];
75 };

## ◆ getResampleType()

static const RESAMPLE filter2d::Filter2d::getResampleType | ( | const std::string | resampleType | ) |
| inline static

Definition at line 76 of file Filter2d.h.

76 {
77 if(resampleType=="near") return(NEAR);
78 else if(resampleType=="bilinear") return(BILINEAR);
79 else{
80 std::string errorString="resampling type not supported: ";
81 errorString+=resampleType;
82 errorString+=" use near or bilinear";
83 throw(errorString);
84 }
85 };

## ◆ linearFeature() [1/2]

void filter2d::Filter2d::linearFeature | ( | const Vector2d< float > & | input,
|  | std::vector< Vector2d< float > > & | output,
|  | float | angle = `361`,
|  | float | angleStep = `1`,
|  | float | maxDistance = `0`,
|  | float | eps = `0`,
|  | bool | l1 = `true`,
|  | bool | a1 = `true`,
|  | bool | l2 = `true`,
|  | bool | a2 = `true`,
|  | bool | verbose = `false`
| ) |  |

Definition at line 1201 of file Filter2d.cc.

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}

## ◆ linearFeature() [2/2]

void filter2d::Filter2d::linearFeature | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | float | angle = `361`,
|  | float | angleStep = `1`,
|  | float | maxDistance = `0`,
|  | float | eps = `0`,
|  | bool | l1 = `true`,
|  | bool | a1 = `true`,
|  | bool | l2 = `true`,
|  | bool | a2 = `true`,
|  | int | band = `0`,
|  | bool | verbose = `false`
| ) |  |

Definition at line 1190 of file Filter2d.cc.

1190 {
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}

## ◆ majorVoting()

void filter2d::Filter2d::majorVoting | ( | const std::string & | inputFilename,
|  | const std::string & | outputFilename,
|  | int | dim = `0`,
|  | const std::vector< int > & | prior = `std::vector<int>()`
| ) |  |

Definition at line 198 of file Filter2d.cc.

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}
ImgReaderGdal
Definition: ImgReaderGdal.h:37
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
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

## ◆ median()

void filter2d::Filter2d::median | ( | const std::string & | inputFilename,
|  | const std::string & | outputFilename,
|  | int | dim,
|  | bool | disc = `false`
| ) |  |

Definition at line 334 of file Filter2d.cc.

335{
336 ImgReaderGdal input;
337 ImgWriterGdal output;
338 input.open(inputFilename);
339 output.open(outputFilename,input);
340 doit(input,output,"median",dim,disc);
341}

## ◆ morphology() [1/2]

template<class T >

unsigned long int filter2d::Filter2d::morphology | ( | const Vector2d< T > & | input,
|  | Vector2d< T > & | output,
|  | const std::string & | method,
|  | int | dimX,
|  | int | dimY,
|  | bool | disc = `false`,
|  | double | hThreshold = `0`
| ) |  |

Definition at line 693 of file Filter2d.h.

694{
695 const char* pszMessage;
696 void* pProgressArg=NULL;
697 GDALProgressFunc pfnProgress=GDALTermProgress;
698 double progress=0;
699 pfnProgress(progress,pszMessage,pProgressArg);
700
701 double noDataValue=0;
702 if(m_noDataValues.size())
703 noDataValue=m_noDataValues[0];
704
705 unsigned long int nchange=0;
706 assert(dimX);
707 assert(dimY);
708 statfactory::StatFactory stat;
709 Vector2d<T> inBuffer(dimY,input.nCols());
710 output.clear();
711 output.resize(input.nRows(),input.nCols());
712 int indexI=0;
713 int indexJ=0;
714 //initialize last half of inBuffer
715 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
716 for(int i=0;i<input.nCols();++i)
717 inBuffer[indexJ][i]=input[abs(j)][i];
718 ++indexJ;
719 }
720 for(int y=0;y<input.nRows();++y){
721 if(y){//inBuffer already initialized for y=0
722 //erase first line from inBuffer
723 inBuffer.erase(inBuffer.begin());
724 //read extra line and push back to inBuffer if not out of bounds
725 if(y+dimY/2<input.nRows()){
726 //allocate buffer
727 inBuffer.push_back(inBuffer.back());
728 for(int i=0;i<input.nCols();++i)
729 inBuffer[inBuffer.size()-1][i]=input[y+dimY/2][i];
730 }
731 else{
732 int over=y+dimY/2-input.nRows();
733 int index=(inBuffer.size()-1)-over;
734 assert(index>=0);
735 assert(index<inBuffer.size());
736 inBuffer.push_back(inBuffer[index]);
737 }
738 }
739 for(int x=0;x<input.nCols();++x){
740 output[y][x]=0;
741 double currentValue=inBuffer[(dimY-1)/2][x];
742 std::vector<double> statBuffer;
743 bool currentMasked=false;
744 for(int imask=0;imask<m_noDataValues.size();++imask){
745 if(currentValue==m_noDataValues[imask]){
746 currentMasked=true;
747 break;
748 }
749 }
750 output[y][x]=currentValue;//introduced due to hThreshold
751 if(currentMasked){
752 output[y][x]=currentValue;
753 }
754 else{
755 for(int j=-(dimY-1)/2;j<=dimY/2;++j){
756 for(int i=-(dimX-1)/2;i<=dimX/2;++i){
757 double d2=i*i+j*j;//square distance
758 if(disc&&(d2>(dimX/2)*(dimY/2)))
759 continue;
760 indexI=x+i;
761 //check if out of bounds
762 if(indexI<0)
763 indexI=-indexI;
764 else if(indexI>=input.nCols())
765 indexI=input.nCols()-i;
766 if(y+j<0)
767 indexJ=-j;
768 else if(y+j>=input.nRows())
769 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
770 else
771 indexJ=(dimY-1)/2+j;
772 if(inBuffer[indexJ][indexI]==noDataValue)
773 continue;
774 bool masked=false;
775 for(int imask=0;imask<m_noDataValues.size();++imask){
776 if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
777 masked=true;
778 break;
779 }
780 }
781 if(!masked){
782 short binValue=0;
783 for(int iclass=0;iclass<m_class.size();++iclass){
784 if(inBuffer[indexJ][indexI]==m_class[iclass]){
785 binValue=1;
786 break;
787 }
788 }
789 if(m_class.size())
790 statBuffer.push_back(binValue);
791 else
792 statBuffer.push_back(inBuffer[indexJ][indexI]);
793 }
794 }
795 }
796 if(statBuffer.size()){
797 switch(getFilterType(method)){
798 case(filter2d::dilate):
799 if(output[y][x]<stat.mymax(statBuffer)-hThreshold){
800 output[y][x]=stat.mymax(statBuffer);
801 ++nchange;
802 }
803 break;
804 case(filter2d::erode):
805 if(output[y][x]>stat.mymin(statBuffer)+hThreshold){
806 output[y][x]=stat.mymin(statBuffer);
807 ++nchange;
808 }
809 break;
810 default:
811 std::ostringstream ess;
812 ess << "Error: morphology method " << method << " not supported, choose " << filter2d::dilate << " (dilate) or " << filter2d::erode << " (erode)" << std::endl;
813 throw(ess.str());
814 break;
815 }
816 if(output[y][x]&&m_class.size())
817 output[y][x]=m_class[0];
818 // else{
819 // assert(m_noDataValues.size());
820 // output[x]=m_noDataValues[0];
821 // }
822 }
823 else
824 output[y][x]=noDataValue;
825 }
826 }
827 progress=(1.0+y);
828 progress/=output.nRows();
829 pfnProgress(progress,pszMessage,pProgressArg);
830 }
831 return nchange;
832}

## ◆ morphology() [2/2]

void filter2d::Filter2d::morphology | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | const std::string & | method,
|  | int | dimX,
|  | int | dimY,
|  | const std::vector< double > & | angle,
|  | bool | disc = `false`
| ) |  |

Definition at line 973 of file Filter2d.cc.

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}

## ◆ mrf() [1/2]

void filter2d::Filter2d::mrf | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | int | dimX,
|  | int | dimY,
|  | double | beta,
|  | bool | eightConnectivity = `true`,
|  | short | down = `1`,
|  | bool | verbose = `false`
| ) |  |

Definition at line 704 of file Filter2d.cc.

704 {
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}

## ◆ mrf() [2/2]

void filter2d::Filter2d::mrf | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | int | dimX,
|  | int | dimY,
|  | Vector2d< double > | beta,
|  | bool | eightConnectivity = `true`,
|  | short | down = `1`,
|  | bool | verbose = `false`
| ) |  |

Definition at line 714 of file Filter2d.cc.

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}

## ◆ pushClass()

void filter2d::Filter2d::pushClass | ( | short | theClass = `1` | ) |
| inline

Definition at line 89 of file Filter2d.h.

89{m_class.push_back(theClass);};

## ◆ pushNoDataValue()

int filter2d::Filter2d::pushNoDataValue | ( | double | noDataValue = `0` | ) |

Definition at line 37 of file Filter2d.cc.

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}

## ◆ pushThreshold()

void filter2d::Filter2d::pushThreshold | ( | double | theThreshold | ) |
| inline

Definition at line 91 of file Filter2d.h.

91{m_threshold.push_back(theThreshold);};

## ◆ setClasses()

void filter2d::Filter2d::setClasses | ( | const std::vector< short > & | theClasses | ) |
| inline

Definition at line 93 of file Filter2d.h.

93{m_class=theClasses;};

## ◆ setTaps()

void filter2d::Filter2d::setTaps | ( | const Vector2d< double > & | taps | ) |

Definition at line 44 of file Filter2d.cc.

45{
46 m_taps=taps;
47}

## ◆ setThresholds()

void filter2d::Filter2d::setThresholds | ( | const std::vector< double > & | theThresholds | ) |
| inline

Definition at line 92 of file Filter2d.h.

92{m_threshold=theThresholds;};

## ◆ shadowDsm() [1/2]

template<class T >

void filter2d::Filter2d::shadowDsm | ( | const Vector2d< T > & | input,
|  | Vector2d< T > & | output,
|  | double | sza,
|  | double | saa,
|  | double | pixelSize,
|  | short | shadowFlag = `1`
| ) |  |

Definition at line 1202 of file Filter2d.h.

1203{
1204 unsigned int ncols=input.nCols();
1205 output.clear();
1206 output.resize(input.nRows(),ncols);
1207 //do we need to initialize output?
1208 // for(int y=0;y<output.nRows();++y)
1209 // for(int x=0;x<output.nCols();++x)
1210 // output[y][x]=0;
1211 int indexI=0;
1212 int indexJ=0;
1213 const char* pszMessage;
1214 void* pProgressArg=NULL;
1215 GDALProgressFunc pfnProgress=GDALTermProgress;
1216 double progress=0;
1217 pfnProgress(progress,pszMessage,pProgressArg);
1218 for(int y=0;y<input.nRows();++y){
1219 for(int x=0;x<input.nCols();++x){
1220 double currentValue=input[y][x];
1221 int theDist=static_cast<int>(sqrt((currentValue*tan(DEG2RAD(sza))/pixelSize)*(currentValue*tan(DEG2RAD(sza))/pixelSize)));//in pixels
1222 double theDir=DEG2RAD(saa)+PI/2.0;
1223 if(theDir<0)
1224 theDir+=2*PI;
1225 for(int d=0;d<theDist;++d){//d in pixels
1226 indexI=x+d*cos(theDir);//in pixels
1227 indexJ=y+d*sin(theDir);//in pixels
1228 if(indexJ<0||indexJ>=input.size())
1229 continue;
1230 if(indexI<0||indexI>=input[indexJ].size())
1231 continue;
1232 if(input[indexJ][indexI]<currentValue-d*pixelSize/tan(DEG2RAD(sza))){//in m
1233 output[indexJ][indexI]=shadowFlag;
1234 }
1235 }
1236 }
1237 progress=(1.0+y);
1238 progress/=output.nRows();
1239 pfnProgress(progress,pszMessage,pProgressArg);
1240 }
1241}

## ◆ shadowDsm() [2/2]

void filter2d::Filter2d::shadowDsm | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | double | sza,
|  | double | saa,
|  | double | pixelSize,
|  | short | shadowFlag = `1`
| ) |  |

Definition at line 1152 of file Filter2d.cc.

1152 {
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}

## ◆ shift() [1/2]

template<class T >

void filter2d::Filter2d::shift | ( | const Vector2d< T > & | input,
|  | Vector2d< T > & | output,
|  | double | offsetX = `0`,
|  | double | offsetY = `0`,
|  | double | randomSigma = `0`,
|  | RESAMPLE | resample = `NEAR`,
|  | bool | verbose = `false`
| ) |  |

Definition at line 610 of file Filter2d.h.

610 {
611 output.resize(input.nRows(),input.nCols());
612 const gsl_rng_type *rangenType;
613 gsl_rng *rangen;
614 gsl_rng_env_setup();
615 rangenType=gsl_rng_default;
616 rangen=gsl_rng_alloc(rangenType);
617 long seed=time(NULL)*getpid();
618 gsl_rng_set(rangen,seed);
619 const char* pszMessage;
620 void* pProgressArg=NULL;
621 GDALProgressFunc pfnProgress=GDALTermProgress;
622 double progress=0;
623 pfnProgress(progress,pszMessage,pProgressArg);
624 for(int j=0;j<input.nRows();++j){
625 for(int i=0;i<input.nCols();++i){
626 T theValue=0;
627 double randomX=0;
628 double randomY=0;
629 if(randomSigma>0){
630 randomX=gsl_ran_gaussian(rangen,randomSigma);
631 randomY=gsl_ran_gaussian(rangen,randomSigma);
632 }
633 double readCol=i+offsetX+randomX;
634 double readRow=j+offsetY+randomY;
635 if(readRow<0)
636 readRow=0;
637 if(readRow>input.nRows()-1)
638 readRow=input.nRows()-1;
639 if(readCol<0)
640 readCol=0;
641 if(readCol>input.nCols()-1)
642 readCol=input.nCols()-1;
643 switch(resample){
644 case(BILINEAR):{
645 double lowerRow=readRow-0.5;
646 double upperRow=readRow+0.5;
647 lowerRow=static_cast<int>(lowerRow);
648 upperRow=static_cast<int>(upperRow);
649 double lowerCol=readCol-0.5;
650 double upperCol=readCol+0.5;
651 lowerCol=static_cast<int>(lowerCol);
652 upperCol=static_cast<int>(upperCol);
653 assert(lowerRow>=0);
654 assert(lowerRow<input.nRows());
655 assert(lowerCol>=0);
656 assert(lowerCol<input.nCols());
657 assert(upperRow>=0);
658 assert(upperRow<input.nRows());
659 assert(upperCol>=0);
660 if(upperCol>=input.nCols()){
661 std::cout << "upperCol: " << upperCol << std::endl;
662 std::cout << "readCol: " << readCol << std::endl;
663 std::cout << "readCol+0.5: " << readCol+0.5 << std::endl;
664 std::cout << "static_cast<int>(readCol+0.5): " << static_cast<int>(readCol+0.5) << std::endl;
665 }
666 assert(upperCol<input.nCols());
667 double c00=input[lowerRow][lowerCol];
668 double c11=input[upperRow][upperCol];
669 double c01=input[lowerRow][upperCol];
670 double c10=input[upperRow][lowerCol];
671 double a=(upperCol-readCol)*c00+(readCol-lowerCol)*c01;
672 double b=(upperCol-readCol)*c10+(readCol-lowerCol)*c11;
673 theValue=(upperRow-readRow)*a+(readRow-lowerRow)*b;
674 break;
675 }
676 default:
677 theValue=input[static_cast<int>(readRow)][static_cast<int>(readCol)];
678 break;
679 }
680 assert(j>=0);
681 assert(j<output.nRows());
682 assert(i>=0);
683 assert(i<output.nCols());
684 output[j][i]=theValue;
685 }
686 progress=(1.0+j);
687 progress/=output.nRows();
688 pfnProgress(progress,pszMessage,pProgressArg);
689 }
690 gsl_rng_free(rangen);
691}

## ◆ shift() [2/2]

void filter2d::Filter2d::shift | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | double | offsetX = `0`,
|  | double | offsetY = `0`,
|  | double | randomSigma = `0`,
|  | RESAMPLE | resample = `BILINEAR`,
|  | bool | verbose = `false`
| ) |  |

Definition at line 846 of file Filter2d.cc.

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}

## ◆ smooth() [1/4]

template<class T1 , class T2 >

void filter2d::Filter2d::smooth | ( | const Vector2d< T1 > & | inputVector,
|  | Vector2d< T2 > & | outputVector,
|  | int | dim
| ) |  |

Definition at line 186 of file Filter2d.h.

187 {
188 smooth(inputVector,outputVector,dim,dim);
189 }

## ◆ smooth() [2/4]

template<class T1 , class T2 >

void filter2d::Filter2d::smooth | ( | const Vector2d< T1 > & | inputVector,
|  | Vector2d< T2 > & | outputVector,
|  | int | dimX,
|  | int | dimY
| ) |  |

Definition at line 191 of file Filter2d.h.

192 {
193 m_taps.resize(dimY);
194 for(int j=0;j<dimY;++j){
195 m_taps[j].resize(dimX);
196 for(int i=0;i<dimX;++i)
197 m_taps[j][i]=1.0/dimX/dimY;
198 }
199 filter(inputVector,outputVector);
200 }

## ◆ smooth() [3/4]

void filter2d::Filter2d::smooth | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | int | dim
| ) |  |

Definition at line 54 of file Filter2d.cc.

55{
56 smooth(input, output,dim,dim);
57}

## ◆ smooth() [4/4]

void filter2d::Filter2d::smooth | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | int | dimX,
|  | int | dimY
| ) |  |

Definition at line 70 of file Filter2d.cc.

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}

## ◆ smoothNoData() [1/2]

void filter2d::Filter2d::smoothNoData | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | int | dim
| ) |  |

Definition at line 49 of file Filter2d.cc.

50{
51 smoothNoData(input, output,dim,dim);
52}

## ◆ smoothNoData() [2/2]

void filter2d::Filter2d::smoothNoData | ( | ImgReaderGdal & | input,
|  | ImgWriterGdal & | output,
|  | int | dimX,
|  | int | dimY
| ) |  |

Definition at line 59 of file Filter2d.cc.

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}

## ◆ var()

void filter2d::Filter2d::var | ( | const std::string & | inputFilename,
|  | const std::string & | outputFilename,
|  | int | dim,
|  | bool | disc = `false`
| ) |  |

Definition at line 343 of file Filter2d.cc.

344{
345 ImgReaderGdal input;
346 ImgWriterGdal output;
347 input.open(inputFilename);
348 output.open(outputFilename,input);
349 doit(input,output,"var",dim,disc);
350}

The documentation for this class was generated from the following files:

- /tmp/slackdce/slackrepo.NvCr9x/build_pktools/PKTOOLS-2.6.7.6/src/algorithms/Filter2d.h
- /tmp/slackdce/slackrepo.NvCr9x/build_pktools/PKTOOLS-2.6.7.6/src/algorithms/Filter2d.cc