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 ( )
inlinevirtual

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

◆ 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
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}
int nrOfRow(void) const
Get the number of rows of this dataset.
int nrOfBand(void) const
Get the number of bands of this dataset.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
void readData(T &value, int col, int row, int band=0)
Read a single pixel cell value at a specific column and row for a specific band (all indices start co...
Definition: ImgReaderGdal.h:95
bool writeData(T &value, int col, int row, int band=0)
Write a single pixel cell value at a specific column and row for a specific band (all indices start c...
Definition: ImgWriterGdal.h:96

◆ 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);
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);
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);
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);
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}
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...
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...

◆ 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)
inlinestatic

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)
inlinestatic

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}
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
void close(void)
Close the image.
void open(const std::string &filename, const ImgReaderGdal &imgSrc, const std::vector< std::string > &options=std::vector< std::string >())
Open an image for writing, copying image attributes from a source image. Image is directly written to...

◆ 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);
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
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: