25#include "StatFactory.h"
28filter2d::Filter2d::Filter2d(
void)
37int filter2d::Filter2d::pushNoDataValue(
double noDataValue)
39 if(find(m_noDataValues.begin(),m_noDataValues.end(),noDataValue)==m_noDataValues.end())
40 m_noDataValues.push_back(noDataValue);
41 return(m_noDataValues.size());
51 smoothNoData(input, output,dim,dim);
56 smooth(input, output,dim,dim);
62 for(
int j=0;j<dimY;++j){
63 m_taps[j].resize(dimX);
64 for(
int i=0;i<dimX;++i)
67 filter(input,output,
false,
true,
true);
73 for(
int j=0;j<dimY;++j){
74 m_taps[j].resize(dimX);
75 for(
int i=0;i<dimX;++i)
78 filter(input,output,
false,
true,
false);
84 int dimX=m_taps[0].size();
85 int dimY=m_taps.size();
87 const char* pszMessage;
88 void* pProgressArg=NULL;
89 GDALProgressFunc pfnProgress=GDALTermProgress;
91 pfnProgress(progress,pszMessage,pProgressArg);
92 for(
int iband=0;iband<input.
nrOfBand();++iband){
94 std::vector<double> outBuffer(input.
nrOfCol());
98 for(
int j=-(dimY-1)/2;j<=dimY/2;++j){
100 input.
readData(inBuffer[indexJ],abs(j),iband);
102 catch(std::string errorstring){
103 std::cerr << errorstring <<
"in line " << indexJ << std::endl;
108 for(
int y=0;y<input.
nrOfRow();++y){
112 inBuffer.erase(inBuffer.begin());
117 inBuffer.push_back(inBuffer.back());
119 input.
readData(inBuffer[inBuffer.size()-1],y+dimY/2,iband);
121 catch(std::string errorstring){
122 std::cerr << errorstring <<
"in band " << iband <<
", line " << y << std::endl;
126 int over=y+dimY/2-input.
nrOfRow();
127 int index=(inBuffer.size()-1)-over;
129 assert(index<inBuffer.size());
130 inBuffer.push_back(inBuffer[index]);
133 for(
int x=0;x<input.
nrOfCol();++x){
138 for(
int imask=0;imask<m_noDataValues.size();++imask){
139 if(inBuffer[(dimY-1)/2][x]==m_noDataValues[imask]){
145 outBuffer[x]=inBuffer[(dimY-1)/2][x];
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){
157 else if(x>=input.
nrOfCol()-(dimX-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);
165 for(
int imask=0;imask<m_noDataValues.size();++imask){
166 if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
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];
178 outBuffer[x]=(normalize&&norm)? abs(outBuffer[x])/norm : abs(outBuffer[x]);
179 else if(normalize&&norm!=0)
180 outBuffer[x]=outBuffer[x]/norm;
186 catch(std::string errorstring){
187 std::cerr << errorstring <<
"in band " << iband <<
", line " << y << std::endl;
190 progress+=(output.
nrOfRow()*iband);
192 pfnProgress(progress,pszMessage,pProgressArg);
198void filter2d::Filter2d::majorVoting(
const std::string& inputFilename,
const std::string& outputFilename,
int dim,
const std::vector<int> &prior)
200 const char* pszMessage;
201 void* pProgressArg=NULL;
202 GDALProgressFunc pfnProgress=GDALTermProgress;
204 pfnProgress(progress,pszMessage,pProgressArg);
208 std::cout <<
"no prior information" << std::endl;
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;
220 input.
open(inputFilename);
221 output.
open(outputFilename,input);
229 dimX=m_taps[0].size();
237 std::vector<double> outBuffer(input.
nrOfCol());
241 for(
int j=-(dimY-1)/2;j<=dimY/2;++j){
243 input.
readData(inBuffer[indexJ],abs(j));
245 catch(std::string errorstring){
246 std::cerr << errorstring <<
"in line " << indexJ << std::endl;
251 for(
int y=0;y<input.
nrOfRow();++y){
255 inBuffer.erase(inBuffer.begin());
260 inBuffer.push_back(inBuffer.back());
262 input.
readData(inBuffer[inBuffer.size()-1],y+dimY/2);
264 catch(std::string errorstring){
265 std::cerr << errorstring <<
"in line" << y << std::endl;
269 int over=y+dimY/2-input.
nrOfRow();
270 int index=(inBuffer.size()-1)-over;
272 assert(index<inBuffer.size());
273 inBuffer.push_back(inBuffer[index]);
276 for(
int x=0;x<input.
nrOfCol();++x){
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){
286 else if(indexI>=input.
nrOfCol())
291 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
304 occurrence[inBuffer[indexJ][indexI]]+=prior[inBuffer[indexJ][indexI]-1];
307 ++occurrence[inBuffer[indexJ][indexI]];
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)
315 if(occurrence[inBuffer[(dimY-1)/2][x]]<maxit->second)
316 outBuffer[x]=maxit->first;
318 outBuffer[x]=inBuffer[(dimY-1)/2][x];
324 catch(std::string errorstring){
325 std::cerr << errorstring <<
"in line" << y << std::endl;
327 progress=(1.0+y)/output.
nrOfRow();
328 pfnProgress(progress,pszMessage,pProgressArg);
334void filter2d::Filter2d::median(
const std::string& inputFilename,
const std::string& outputFilename,
int dim,
bool disc)
338 input.
open(inputFilename);
339 output.
open(outputFilename,input);
340 doit(input,output,
"median",dim,disc);
343void filter2d::Filter2d::var(
const std::string& inputFilename,
const std::string& outputFilename,
int dim,
bool disc)
347 input.
open(inputFilename);
348 output.
open(outputFilename,input);
349 doit(input,output,
"var",dim,disc);
352void filter2d::Filter2d::doit(
ImgReaderGdal& input,
ImgWriterGdal& output,
const std::string& method,
int dim,
short down,
bool disc)
354 doit(input,output,method,dim,dim,down,disc);
357void filter2d::Filter2d::doit(
ImgReaderGdal& input,
ImgWriterGdal& output,
const std::string& method,
int dimX,
int dimY,
short down,
bool disc)
359 const char* pszMessage;
360 void* pProgressArg=NULL;
361 GDALProgressFunc pfnProgress=GDALTermProgress;
363 pfnProgress(progress,pszMessage,pProgressArg);
369 for(
int iband=0;iband<input.
nrOfBand();++iband){
371 std::vector<double> outBuffer((input.
nrOfCol()+down-1)/down);
375 for(
int j=-(dimY-1)/2;j<=dimY/2;++j){
377 input.
readData(inBuffer[indexJ],abs(j),iband);
379 catch(std::string errorstring){
380 std::cerr << errorstring <<
"in line " << indexJ << std::endl;
384 for(
int y=0;y<input.
nrOfRow();++y){
388 inBuffer.erase(inBuffer.begin());
393 inBuffer.push_back(inBuffer.back());
395 input.
readData(inBuffer[inBuffer.size()-1],y+dimY/2,iband);
397 catch(std::string errorstring){
398 std::cerr << errorstring <<
"in band " << iband <<
", line " << y << std::endl;
402 int over=y+dimY/2-input.
nrOfRow();
403 int index=(inBuffer.size()-1)-over;
405 assert(index<inBuffer.size());
406 inBuffer.push_back(inBuffer[index]);
409 if((y+1+down/2)%down)
411 for(
int x=0;x<input.
nrOfCol();++x){
412 if((x+1+down/2)%down)
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){
421 if(disc&&(d2>(dimX/2)*(dimY/2)))
427 else if(indexI>=input.
nrOfCol())
432 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
436 for(
int imask=0;imask<m_noDataValues.size();++imask){
437 if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
443 std::vector<short>::const_iterator vit=m_class.begin();
445 ++occurrence[inBuffer[indexJ][indexI]];
447 while(vit!=m_class.end()){
448 if(inBuffer[indexJ][indexI]==*(vit++))
449 ++occurrence[inBuffer[indexJ][indexI]];
452 windowBuffer.push_back(inBuffer[indexJ][indexI]);
456 switch(getFilterType(method)){
457 case(filter2d::nvalid):
458 outBuffer[x/down]=stat.nvalid(windowBuffer);
460 case(filter2d::median):
461 if(windowBuffer.empty())
462 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
464 outBuffer[x/down]=stat.median(windowBuffer);
466 case(filter2d::var):{
467 if(windowBuffer.empty())
468 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
470 outBuffer[x/down]=stat.var(windowBuffer);
473 case(filter2d::stdev):{
474 if(windowBuffer.empty())
475 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
477 outBuffer[x/down]=sqrt(stat.var(windowBuffer));
480 case(filter2d::mean):{
481 if(windowBuffer.empty())
482 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
484 outBuffer[x/down]=stat.mean(windowBuffer);
487 case(filter2d::min):{
488 if(windowBuffer.empty())
489 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
491 outBuffer[x/down]=stat.mymin(windowBuffer);
494 case(filter2d::ismin):{
495 if(windowBuffer.empty())
496 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
498 outBuffer[x/down]=(stat.mymin(windowBuffer)==windowBuffer[centre])? 1:0;
501 case(filter2d::minmax):{
504 if(windowBuffer.empty())
505 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
507 stat.minmax(windowBuffer,windowBuffer.begin(),windowBuffer.end(),min,max);
511 outBuffer[x/down]=windowBuffer[centre];
515 case(filter2d::max):{
516 if(windowBuffer.empty())
517 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
519 outBuffer[x/down]=stat.mymax(windowBuffer);
522 case(filter2d::ismax):{
523 if(windowBuffer.empty())
524 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
526 outBuffer[x/down]=(stat.mymax(windowBuffer)==windowBuffer[centre])? 1:0;
529 case(filter2d::order):{
530 if(windowBuffer.empty())
531 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 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);
542 case(filter2d::sum):{
543 outBuffer[x/down]=stat.sum(windowBuffer);
546 case(filter2d::percentile):{
547 assert(m_threshold.size());
548 outBuffer[x/down]=stat.percentile(windowBuffer,windowBuffer.begin(),windowBuffer.end(),m_threshold[0]);
551 case(filter2d::proportion):{
552 if(windowBuffer.size()){
553 double sum=stat.sum(windowBuffer);
555 outBuffer[x/down]=100.0*windowBuffer[centre]/stat.sum(windowBuffer);
557 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
560 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
563 case(filter2d::homog):
564 if(occurrence.size()==1)
565 outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
567 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
569 case(filter2d::heterog):{
570 if(occurrence.size()==windowBuffer.size())
571 outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
573 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
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();
600 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
603 case(filter2d::countid):{
604 if(windowBuffer.size())
605 outBuffer[x/down]=occurrence.size();
607 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
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)
617 if(occurrence[inBuffer[(dimY-1)/2][x]]<maxit->second)
618 outBuffer[x/down]=maxit->first;
620 outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
623 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
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];
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];
636 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
639 case(filter2d::scramble):{
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];
647 outBuffer[x/down]=windowBuffer[randomIndex];
650 outBuffer[x/down]=(m_noDataValues.size())? m_noDataValues[0] : 0;
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)){
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;
667 outBuffer[x/down]=MF;
673 outBuffer[x/down]=NF;
677 outBuffer[x/down]=inBuffer[indexJ][indexI];
681 std::ostringstream ess;
682 ess <<
"Error: filter method " << method <<
" not supported" << std::endl;
688 progress=(1.0+y/down);
689 progress+=(output.
nrOfRow()*iband);
691 pfnProgress(progress,pszMessage,pProgressArg);
694 output.
writeData(outBuffer,y/down,iband);
696 catch(std::string errorstring){
697 std::cerr << errorstring <<
"in band " << iband <<
", line " << y << std::endl;
701 pfnProgress(1.0,pszMessage,pProgressArg);
704void filter2d::Filter2d::mrf(
ImgReaderGdal& input,
ImgWriterGdal& output,
int dimX,
int dimY,
double beta,
bool eightConnectivity,
short down,
bool verbose){
705 assert(m_class.size()>1);
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);
716 const char* pszMessage;
717 void* pProgressArg=NULL;
718 GDALProgressFunc pfnProgress=GDALTermProgress;
720 pfnProgress(progress,pszMessage,pProgressArg);
728 assert(output.
nrOfBand()==m_class.size());
729 assert(m_class.size()>1);
730 assert(beta.size()==m_class.size());
734 for(
int j=-(dimY-1)/2;j<=dimY/2;++j){
736 input.
readData(inBuffer[indexJ],abs(j));
738 catch(std::string errorstring){
739 std::cerr << errorstring <<
"in line " << indexJ << std::endl;
743 for(
int y=0;y<input.
nrOfRow();++y){
747 inBuffer.erase(inBuffer.begin());
752 inBuffer.push_back(inBuffer.back());
754 input.
readData(inBuffer[inBuffer.size()-1],y+dimY/2);
756 catch(std::string errorstring){
757 std::cerr << errorstring <<
"in line " << y << std::endl;
761 int over=y+dimY/2-input.
nrOfRow();
762 int index=(inBuffer.size()-1)-over;
764 assert(index<inBuffer.size());
765 inBuffer.push_back(inBuffer[index]);
768 if((y+1+down/2)%down)
770 for(
int x=0;x<input.
nrOfCol();++x){
771 if((x+1+down/2)%down)
773 std::vector<short> potential(m_class.size());
774 for(
int iclass=0;iclass<m_class.size();++iclass){
776 outBuffer[iclass][x/down]=0;
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)
790 else if(indexI>=input.
nrOfCol())
795 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
799 for(
int imask=0;imask<m_noDataValues.size();++imask){
800 if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
806 for(
int iclass=0;iclass<m_class.size();++iclass){
807 if(inBuffer[indexJ][indexI]==m_class[iclass])
808 potential[iclass]+=1;
814 for(
int iclass1=0;iclass1<m_class.size();++iclass1){
815 assert(beta[iclass1].size()==m_class.size());
817 for(
int iclass2=0;iclass2<m_class.size();++iclass2)
819 pot+=potential[iclass2]*beta[iclass1][iclass2];
820 double prior=exp(-pot);
821 outBuffer[iclass1][x/down]=prior;
825 for(
int iclass1=0;iclass1<m_class.size();++iclass1)
826 outBuffer[iclass1][x/down]/=norm;
829 progress=(1.0+y/down)/output.
nrOfRow();
830 pfnProgress(progress,pszMessage,pProgressArg);
832 assert(outBuffer.size()==m_class.size());
834 for(
int iclass=0;iclass<m_class.size();++iclass){
835 assert(outBuffer[iclass].size()==output.
nrOfCol());
837 output.
writeData(outBuffer[iclass],y/down,iclass);
839 catch(std::string errorstring){
840 std::cerr << errorstring <<
"in class " << iclass <<
", line " << y << std::endl;
846void filter2d::Filter2d::shift(
ImgReaderGdal& input,
ImgWriterGdal& output,
double offsetX,
double offsetY,
double randomSigma, RESAMPLE resample,
bool verbose)
851 const char* pszMessage;
852 void* pProgressArg=NULL;
853 GDALProgressFunc pfnProgress=GDALTermProgress;
855 pfnProgress(progress,pszMessage,pProgressArg);
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);
973void filter2d::Filter2d::morphology(
ImgReaderGdal& input,
ImgWriterGdal& output,
const std::string& method,
int dimX,
int dimY,
const std::vector<double> &angle,
bool disc)
975 const char* pszMessage;
976 void* pProgressArg=NULL;
977 GDALProgressFunc pfnProgress=GDALTermProgress;
979 pfnProgress(progress,pszMessage,pProgressArg);
985 for(
int iband=0;iband<input.
nrOfBand();++iband){
987 std::vector<double> outBuffer(input.
nrOfCol());
991 for(
int j=-(dimY-1)/2;j<=dimY/2;++j){
993 input.
readData(inBuffer[indexJ],abs(j),iband);
996 catch(std::string errorstring){
997 std::cerr << errorstring <<
"in line " << indexJ << std::endl;
1000 for(
int y=0;y<input.
nrOfRow();++y){
1004 inBuffer.erase(inBuffer.begin());
1009 inBuffer.push_back(inBuffer.back());
1011 input.
readData(inBuffer[inBuffer.size()-1],y+dimY/2,iband);
1013 catch(std::string errorstring){
1014 std::cerr << errorstring <<
"in band " << iband <<
", line " << y << std::endl;
1018 int over=y+dimY/2-input.
nrOfRow();
1019 int index=(inBuffer.size()-1)-over;
1021 assert(index<inBuffer.size());
1022 inBuffer.push_back(inBuffer[index]);
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]){
1038 outBuffer[x]=currentValue;
1041 for(
int j=-(dimY-1)/2;j<=dimY/2;++j){
1042 for(
int i=-(dimX-1)/2;i<=dimX/2;++i){
1044 if(disc&&(d2>(dimX/2)*(dimY/2)))
1051 theta=atan(
static_cast<double>(-j)/(
static_cast<double>(i)));
1053 theta=-atan(
static_cast<double>(j)/(
static_cast<double>(i)));
1057 theta=PI-atan(
static_cast<double>(-j)/(
static_cast<double>(-i)));
1059 theta=PI+atan(
static_cast<double>(j)/(
static_cast<double>(-i)));
1066 theta=360-(theta/PI*180)+90;
1071 bool alligned=
false;
1072 for(
int iangle=0;iangle<angle.size();++iangle){
1073 if(sqrt((theta-angle[iangle])*(theta-angle[iangle]))<10){
1085 else if(indexI>=input.
nrOfCol())
1090 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
1092 indexJ=(dimY-1)/2+j;
1097 for(
int imask=0;imask<m_noDataValues.size();++imask){
1098 if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
1105 for(
int iclass=0;iclass<m_class.size();++iclass){
1106 if(inBuffer[indexJ][indexI]==m_class[iclass]){
1112 statBuffer.push_back(binValue);
1114 statBuffer.push_back(inBuffer[indexJ][indexI]);
1118 if(statBuffer.size()){
1119 switch(getFilterType(method)){
1120 case(filter2d::dilate):
1121 outBuffer[x]=stat.mymax(statBuffer);
1123 case(filter2d::erode):
1124 outBuffer[x]=stat.mymin(statBuffer);
1127 std::ostringstream ess;
1128 ess <<
"Error: morphology method " << method <<
" not supported, choose " << filter2d::dilate <<
" (dilate) or " << filter2d::erode <<
" (erode)" << std::endl;
1133 if(outBuffer[x]&&m_class.size())
1134 outBuffer[x]=m_class[0];
1141 catch(std::string errorstring){
1142 std::cerr << errorstring <<
"in band " << iband <<
", line " << y << std::endl;
1145 progress+=(output.
nrOfRow()*iband);
1147 pfnProgress(progress,pszMessage,pProgressArg);
1152void filter2d::Filter2d::shadowDsm(
ImgReaderGdal& input,
ImgWriterGdal& output,
double sza,
double saa,
double pixelSize,
short shadowFlag){
1156 shadowDsm(inputBuffer, outputBuffer, sza, saa, pixelSize, shadowFlag);
1160void filter2d::Filter2d::dwtForward(
ImgReaderGdal& input,
ImgWriterGdal& output,
const std::string& wavelet_type,
int family){
1162 for(
int iband=0;iband<input.
nrOfBand();++iband){
1164 std::cout <<
"filtering band " << iband << std::endl << std::flush;
1165 dwtForward(theBuffer, wavelet_type, family);
1170void filter2d::Filter2d::dwtInverse(
ImgReaderGdal& input,
ImgWriterGdal& output,
const std::string& wavelet_type,
int family){
1172 for(
int iband=0;iband<input.
nrOfBand();++iband){
1174 std::cout <<
"filtering band " << iband << std::endl << std::flush;
1175 dwtInverse(theBuffer, wavelet_type, family);
1180void filter2d::Filter2d::dwtCut(
ImgReaderGdal& input,
ImgWriterGdal& output,
const std::string& wavelet_type,
int family,
double cut,
bool verbose){
1182 for(
int iband=0;iband<input.
nrOfBand();++iband){
1184 std::cout <<
"filtering band " << iband << std::endl << std::flush;
1185 dwtCut(theBuffer, wavelet_type, family, cut);
1190void filter2d::Filter2d::linearFeature(
ImgReaderGdal& input,
ImgWriterGdal& output,
float angle,
float angleStep,
float maxDistance,
float eps,
bool l1,
bool a1,
bool l2,
bool a2,
int band,
bool verbose){
1192 std::vector< Vector2d<float> > outputBuffer;
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)
1201void filter2d::Filter2d::linearFeature(
const Vector2d<float>& input, std::vector<
Vector2d<float> >& output,
float angle,
float angleStep,
float maxDistance,
float eps,
bool l1,
bool a1,
bool l2,
bool a2,
bool verbose)
1213 output.resize(nband);
1214 for(
int iband=0;iband<output.size();++iband)
1215 output[iband].resize(input.nRows(),input.nCols());
1217 maxDistance=sqrt(
static_cast<float>(input.nRows()*input.nCols()));
1220 const char* pszMessage;
1221 void* pProgressArg=NULL;
1222 GDALProgressFunc pfnProgress=GDALTermProgress;
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];
1230 float lineDistance1=0;
1231 float lineDistance2=maxDistance;
1235 for(northAngle=0;northAngle<180;northAngle+=angleStep){
1236 if(angle<=360&&angle>=0&&angle!=northAngle)
1240 std::cout <<
"northAngle: " << northAngle << std::endl;
1241 float currentDistance=0;
1243 for(
short side=0;side<=1;side+=1){
1244 theDir=PI/2.0-DEG2RAD(northAngle)+side*PI;
1247 std::cout <<
"theDir in deg: " << RAD2DEG(theDir) << std::endl;
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())
1259 if(indexI<0||indexI>=input[indexJ].size())
1261 nextValue=input[indexJ][indexI];
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;
1275 if(fabs(currentValue-nextValue)<=eps){
1279 std::cout <<
"currentDistance: " << currentDistance <<
", continue" << std::endl;
1283 std::cout <<
"currentDistance: " << currentDistance <<
", break" << std::endl;
1288 if(lineDistance1<currentDistance){
1289 lineDistance1=currentDistance;
1290 lineAngle1=northAngle;
1292 if(lineDistance2>currentDistance){
1293 lineDistance2=currentDistance;
1294 lineAngle2=northAngle;
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;
1305 output[iband++][y][x]=lineDistance1;
1307 output[iband++][y][x]=lineAngle1;
1309 output[iband++][y][x]=lineDistance2;
1311 output[iband++][y][x]=lineAngle2;
1312 assert(iband==nband);
1315 progress/=input.nRows();
1316 pfnProgress(progress,pszMessage,pProgressArg);
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.
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...
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
void readDataBlock(Vector2d< T > &buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band=0)
Read pixel cell values for a range of columns and rows for a specific band (all indices start countin...
void close(void)
Close the image.
void open(const std::string &filename, const ImgReaderGdal &imgSrc, const std::vector< std::string > &options=std::vector< std::string >())
Open an image for writing, copying image attributes from a source image. Image is directly written to...
bool writeDataBlock(Vector2d< T > &buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band=0)
Write pixel cell values for a range of columns and rows for a specific band (all indices start counti...
bool writeData(T &value, int col, int row, int band=0)
Write a single pixel cell value at a specific column and row for a specific band (all indices start c...