20 #ifndef _MYFILTER2D_H_ 21 #define _MYFILTER2D_H_ 24 #define PI 3.1415926535897932384626433832795 28 #define DEG2RAD(DEG) (DEG/180.0*PI) 32 #define RAD2DEG(RAD) (RAD/PI*180) 37 #define getpid _getpid 47 #include <gsl/gsl_sort.h> 48 #include <gsl/gsl_wavelet.h> 49 #include <gsl/gsl_wavelet2d.h> 50 #include <gsl/gsl_rng.h> 51 #include <gsl/gsl_randist.h> 53 #include "base/Vector2d.h" 55 #include "imageclasses/ImgReaderGdal.h" 56 #include "imageclasses/ImgWriterGdal.h" 57 #include "algorithms/StatFactory.h" 61 enum FILTER_TYPE { median=100, var=101 , min=102, max=103, sum=104, mean=105, minmax=106, dilate=107, erode=108, close=109, open=110, homog=111, sobelx=112, sobely=113, sobelxy=114, sobelyx=115, smooth=116, density=117, mode=118, mixed=119, threshold=120, ismin=121, ismax=122, heterog=123, order=124, stdev=125, mrf=126, dwt=127, dwti=128, dwt_cut=129, scramble=130, shift=131, linearfeature=132, smoothnodata=133, countid=134, dwt_cut_from=135, savgolay=136, percentile=137, proportion=138, nvalid=139, sauvola=140};
63 enum RESAMPLE { NEAR = 0, BILINEAR = 1, BICUBIC = 2 };
71 static FILTER_TYPE getFilterType(
const std::string filterType){
72 std::map<std::string, FILTER_TYPE> m_filterMap;
74 return m_filterMap[filterType];
76 static const RESAMPLE getResampleType(
const std::string resampleType){
77 if(resampleType==
"near")
return(NEAR);
78 else if(resampleType==
"bilinear")
return(BILINEAR);
80 std::string errorString=
"resampling type not supported: ";
81 errorString+=resampleType;
82 errorString+=
" use near or bilinear";
89 void pushClass(
short theClass=1){m_class.push_back(theClass);};
90 int pushNoDataValue(
double noDataValue=0);
91 void pushThreshold(
double theThreshold){m_threshold.push_back(theThreshold);};
92 void setThresholds(
const std::vector<double>& theThresholds){m_threshold=theThresholds;};
93 void setClasses(
const std::vector<short>& theClasses){m_class=theClasses;};
101 template<
class T1,
class T2>
void smooth(
const Vector2d<T1>& inputVector,
Vector2d<T2>& outputVector,
int dimX,
int dimY);
104 void dwtCut(
ImgReaderGdal& input,
ImgWriterGdal& output,
const std::string& wavelet_type,
int family,
double cut,
bool verbose=
false);
105 template<
class T>
void dwtForward(
Vector2d<T>& data,
const std::string& wavelet_type,
int family);
106 template<
class T>
void dwtInverse(
Vector2d<T>& data,
const std::string& wavelet_type,
int family);
107 template<
class T>
void dwtCut(
Vector2d<T>& data,
const std::string& wavelet_type,
int family,
double cut);
108 void majorVoting(
const std::string& inputFilename,
const std::string& outputFilename,
int dim=0,
const std::vector<int> &prior=std::vector<int>());
111 void doit(
ImgReaderGdal& input,
ImgWriterGdal& output,
const std::string& method,
int dimX,
int dimY,
short down=1,
bool disc=
false);
112 void mrf(
ImgReaderGdal& input,
ImgWriterGdal& output,
int dimX,
int dimY,
double beta,
bool eightConnectivity=
true,
short down=1,
bool verbose=
false);
114 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);
115 void median(
const std::string& inputFilename,
const std::string& outputFilename,
int dim,
bool disc=
false);
116 void var(
const std::string& inputFilename,
const std::string& outputFilename,
int dim,
bool disc=
false);
117 void morphology(
ImgReaderGdal& input,
ImgWriterGdal& output,
const std::string& method,
int dimX,
int dimY,
const std::vector<double> &angle,
bool disc=
false);
118 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);
119 template<
class T>
unsigned long int dsm2dtm_nwse(
const Vector2d<T>& inputDSM,
Vector2d<T>& outputMask,
double hThreshold,
int nlimit,
int dim=3);
120 template<
class T>
unsigned long int dsm2dtm_nesw(
const Vector2d<T>& inputDSM,
Vector2d<T>& outputMask,
double hThreshold,
int nlimit,
int dim=3);
121 template<
class T>
unsigned long int dsm2dtm_senw(
const Vector2d<T>& inputDSM,
Vector2d<T>& outputMask,
double hThreshold,
int nlimit,
int dim=3);
122 template<
class T>
unsigned long int dsm2dtm_swne(
const Vector2d<T>& inputDSM,
Vector2d<T>& outputMask,
double hThreshold,
int nlimit,
int dim=3);
123 template<
class T>
void shadowDsm(
const Vector2d<T>& input,
Vector2d<T>& output,
double sza,
double saa,
double pixelSize,
short shadowFlag=1);
125 void dwt_texture(
const std::string& inputFilename,
const std::string& outputFilename,
int dim,
int scale,
int down=1,
int iband=0,
bool verbose=
false);
126 void shift(
ImgReaderGdal& input,
ImgWriterGdal& output,
double offsetX=0,
double offsetY=0,
double randomSigma=0, RESAMPLE resample=BILINEAR,
bool verbose=
false);
127 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);
128 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);
129 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);
132 static void initMap(std::map<std::string, FILTER_TYPE>& m_filterMap){
134 m_filterMap[
"median"]=filter2d::median;
135 m_filterMap[
"nvalid"]=filter2d::nvalid;
136 m_filterMap[
"var"]=filter2d::var;
137 m_filterMap[
"min"]=filter2d::min;
138 m_filterMap[
"max"]=filter2d::max;
139 m_filterMap[
"sum"]=filter2d::sum;
140 m_filterMap[
"mean"]=filter2d::mean;
141 m_filterMap[
"minmax"]=filter2d::minmax;
142 m_filterMap[
"dilate"]=filter2d::dilate;
143 m_filterMap[
"erode"]=filter2d::erode;
144 m_filterMap[
"close"]=filter2d::close;
145 m_filterMap[
"open"]=filter2d::open;
146 m_filterMap[
"homog"]=filter2d::homog;
147 m_filterMap[
"sobelx"]=filter2d::sobelx;
148 m_filterMap[
"sobely"]=filter2d::sobely;
149 m_filterMap[
"sobelxy"]=filter2d::sobelxy;
150 m_filterMap[
"sobelyx"]=filter2d::sobelyx;
151 m_filterMap[
"smooth"]=filter2d::smooth;
152 m_filterMap[
"density"]=filter2d::density;
153 m_filterMap[
"mode"]=filter2d::mode;
154 m_filterMap[
"mixed"]=filter2d::mixed;
155 m_filterMap[
"smoothnodata"]=filter2d::smoothnodata;
156 m_filterMap[
"threshold"]=filter2d::threshold;
157 m_filterMap[
"ismin"]=filter2d::ismin;
158 m_filterMap[
"ismax"]=filter2d::ismax;
159 m_filterMap[
"heterog"]=filter2d::heterog;
160 m_filterMap[
"sauvola"]=filter2d::sauvola;
161 m_filterMap[
"order"]=filter2d::order;
162 m_filterMap[
"stdev"]=filter2d::stdev;
163 m_filterMap[
"mrf"]=filter2d::mrf;
164 m_filterMap[
"dwt"]=filter2d::dwt;
165 m_filterMap[
"dwti"]=filter2d::dwti;
166 m_filterMap[
"dwt_cut"]=filter2d::dwt_cut;
167 m_filterMap[
"dwt_cut_from"]=filter2d::dwt_cut_from;
168 m_filterMap[
"scramble"]=filter2d::scramble;
169 m_filterMap[
"shift"]=filter2d::shift;
170 m_filterMap[
"linearfeature"]=filter2d::linearfeature;
171 m_filterMap[
"countid"]=filter2d::countid;
172 m_filterMap[
"savgolay"]=filter2d::savgolay;
173 m_filterMap[
"percentile"]=filter2d::percentile;
174 m_filterMap[
"proportion"]=filter2d::proportion;
179 std::vector<short> m_class;
181 std::vector<double> m_noDataValues;
182 std::vector<double> m_threshold;
186 template<
class T1,
class T2>
void Filter2d::smooth(
const Vector2d<T1>& inputVector,
Vector2d<T2>& outputVector,
int dim)
188 smooth(inputVector,outputVector,dim,dim);
191 template<
class T1,
class T2>
void Filter2d::smooth(
const Vector2d<T1>& inputVector,
Vector2d<T2>& outputVector,
int dimX,
int 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;
199 filter(inputVector,outputVector);
204 outputVector.resize(inputVector.size());
205 int dimX=m_taps[0].size();
206 int dimY=m_taps.size();
208 std::vector<T2> outBuffer(inputVector[0].size());
213 for(
int j=-(dimY-1)/2;j<=dimY/2;++j){
214 inBuffer[indexJ]=inputVector[abs(j)];
218 for(
int y=0;y<inputVector.size();++y){
221 inBuffer.erase(inBuffer.begin());
223 if(y+dimY/2<inputVector.size()){
225 inBuffer.push_back(inputVector[y+dimY/2]);
228 int over=y+dimY/2-inputVector.nRows();
229 int index=(inBuffer.size()-1)-over;
231 assert(index<inBuffer.size());
232 inBuffer.push_back(inBuffer[index]);
235 for(
int x=0;x<inputVector.nCols();++x){
237 for(
int j=-(dimY-1)/2;j<=dimY/2;++j){
238 for(
int i=-(dimX-1)/2;i<=dimX/2;++i){
244 else if(x>=inputVector.nCols()-(dimX-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]);
255 outputVector[y]=outBuffer;
259 template<
class T1,
class T2>
void Filter2d::doit(
const Vector2d<T1>& inputVector,
Vector2d<T2>& outputVector,
const std::string& method,
int dimX,
int dimY,
short down,
bool disc)
261 const char* pszMessage;
262 void* pProgressArg=NULL;
263 GDALProgressFunc pfnProgress=GDALTermProgress;
265 pfnProgress(progress,pszMessage,pProgressArg);
268 double noDataValue=0;
269 if(m_noDataValues.size()){
270 stat.setNoDataValues(m_noDataValues);
271 noDataValue=m_noDataValues[0];
276 outputVector.resize((inputVector.size()+down-1)/down);
278 std::vector<T2> outBuffer((inputVector[0].size()+down-1)/down);
283 for(
int j=-(dimY-1)/2;j<=dimY/2;++j){
284 inBuffer[indexJ]=inputVector[abs(j)];
287 for(
int y=0;y<inputVector.size();++y){
291 inBuffer.erase(inBuffer.begin());
293 if(y+dimY/2<inputVector.size())
294 inBuffer.push_back(inputVector[y+dimY/2]);
296 int over=y+dimY/2-inputVector.size();
297 int index=(inBuffer.size()-1)-over;
299 assert(index<inBuffer.size());
300 inBuffer.push_back(inBuffer[index]);
303 if((y+1+down/2)%down)
305 for(
int x=0;x<inputVector[0].size();++x){
306 if((x+1+down/2)%down)
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){
320 else if(indexI>=inputVector[0].size())
321 indexI=inputVector[0].size()-i;
324 else if(y+j>=inputVector.size())
325 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
328 windowBuffer.push_back(inBuffer[indexJ][indexI]);
329 if(!stat.isNoData(inBuffer[indexJ][indexI])){
330 std::vector<short>::const_iterator vit=m_class.begin();
333 ++occurrence[inBuffer[indexJ][indexI]];
335 while(vit!=m_class.end()){
336 if(inBuffer[indexJ][indexI]==*(vit++))
337 ++occurrence[inBuffer[indexJ][indexI]];
343 switch(getFilterType(method)){
344 case(filter2d::nvalid):
345 outBuffer[x/down]=stat.nvalid(windowBuffer);
347 case(filter2d::median):
348 outBuffer[x/down]=stat.median(windowBuffer);
350 case(filter2d::var):{
351 outBuffer[x/down]=stat.var(windowBuffer);
354 case(filter2d::stdev):{
355 T2 varValue=stat.var(windowBuffer);
356 if(stat.isNoData(varValue))
357 outBuffer[x/down]=noDataValue;
359 outBuffer[x/down]=sqrt(varValue);
362 case(filter2d::mean):{
363 if(windowBuffer.empty())
364 outBuffer[x/down]=noDataValue;
366 outBuffer[x/down]=stat.mean(windowBuffer);
369 case(filter2d::min):{
370 outBuffer[x/down]=stat.mymin(windowBuffer);
373 case(filter2d::ismin):{
374 T1 minValue=stat.mymin(windowBuffer);
375 if(stat.isNoData(minValue))
376 outBuffer[x/down]=noDataValue;
378 outBuffer[x/down]=(windowBuffer[centre]==minValue)? 1:0;
381 case(filter2d::minmax):{
384 stat.minmax(windowBuffer,windowBuffer.begin(),windowBuffer.end(),min,max);
388 outBuffer[x/down]=windowBuffer[centre];
391 case(filter2d::max):{
392 outBuffer[x/down]=stat.mymax(windowBuffer);
395 case(filter2d::ismax):{
396 T1 maxValue=stat.mymax(windowBuffer);
397 if(stat.isNoData(maxValue))
398 outBuffer[x/down]=noDataValue;
400 outBuffer[x/down]=(windowBuffer[centre]==maxValue)? 1:0;
403 case(filter2d::order):{
404 stat.eraseNoData(windowBuffer);
405 if(windowBuffer.empty())
406 outBuffer[x/down]=noDataValue;
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);
417 case(filter2d::sum):{
418 outBuffer[x/down]=stat.sum(windowBuffer);
421 case(filter2d::percentile):{
422 assert(m_threshold.size());
423 outBuffer[x/down]=stat.percentile(windowBuffer,windowBuffer.begin(),windowBuffer.end(),m_threshold[0]);
426 case(filter2d::proportion):{
427 stat.eraseNoData(windowBuffer);
428 T2 sum=stat.sum(windowBuffer);
430 outBuffer[x/down]=windowBuffer[centre]/sum;
432 outBuffer[x/down]=noDataValue;
435 case(filter2d::homog):{
436 T1 centreValue=inBuffer[(dimY-1)/2][x];
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)
451 outBuffer[x/down]=noDataValue;
454 case(filter2d::heterog):{
455 T1 centreValue=inBuffer[(dimY-1)/2][x];
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)
470 outBuffer[x/down]=noDataValue;
473 case(filter2d::sauvola):{
478 T1 centreValue=inBuffer[(dimY-1)/2][x];
479 if(windowBuffer.empty()||stat.isNoData(centreValue)){
483 stat.meanVar(windowBuffer,theMean,theStdev);
484 theStdev=sqrt(theStdev);
487 if(m_threshold.size()==2){
488 kValue=m_threshold[0];
489 rValue=m_threshold[1];
493 double theThreshold=theMean*(1+kValue*(theStdev/rValue - 1));
495 outBuffer[x/down]=(centreValue>theThreshold) ? 1 : noDataValue;
498 outBuffer[x/down]=noDataValue;
502 case(filter2d::density):{
503 int nvalid=stat.nvalid(windowBuffer);
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;
510 outBuffer[x/down]=noDataValue;
513 case(filter2d::countid):{
514 if(occurrence.size())
515 outBuffer[x/down]=occurrence.size();
517 outBuffer[x/down]=noDataValue;
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)
527 if(occurrence[inBuffer[(dimY-1)/2][x]]<maxit->second)
528 outBuffer[x/down]=maxit->first;
530 outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
533 outBuffer[x/down]=noDataValue;
536 case(filter2d::threshold):{
537 assert(m_class.size()==m_threshold.size());
538 int nvalid=stat.nvalid(windowBuffer);
540 outBuffer[x/down]=inBuffer[(dimY-1)/2][x];
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];
547 outBuffer[x/down]=noDataValue;
550 case(filter2d::scramble):{
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];
558 outBuffer[x/down]=windowBuffer[randomIndex];
561 outBuffer[x/down]=noDataValue;
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)){
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;
578 outBuffer[x/down]=MF;
584 outBuffer[x/down]=NF;
588 outBuffer[x/down]=inBuffer[indexJ][indexI];
595 progress=(1.0+y/down);
596 progress/=outputVector.size();
597 pfnProgress(progress,pszMessage,pProgressArg);
599 outputVector[y/down]=outBuffer;
610 template<
class T>
void Filter2d::shift(
const Vector2d<T>& input,
Vector2d<T>& output,
double offsetX,
double offsetY,
double randomSigma, RESAMPLE resample,
bool verbose){
611 output.resize(input.nRows(),input.nCols());
612 const gsl_rng_type *rangenType;
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;
623 pfnProgress(progress,pszMessage,pProgressArg);
624 for(
int j=0;j<input.nRows();++j){
625 for(
int i=0;i<input.nCols();++i){
630 randomX=gsl_ran_gaussian(rangen,randomSigma);
631 randomY=gsl_ran_gaussian(rangen,randomSigma);
633 double readCol=i+offsetX+randomX;
634 double readRow=j+offsetY+randomY;
637 if(readRow>input.nRows()-1)
638 readRow=input.nRows()-1;
641 if(readCol>input.nCols()-1)
642 readCol=input.nCols()-1;
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);
654 assert(lowerRow<input.nRows());
656 assert(lowerCol<input.nCols());
658 assert(upperRow<input.nRows());
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;
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;
677 theValue=input[
static_cast<int>(readRow)][static_cast<int>(readCol)];
681 assert(j<output.nRows());
683 assert(i<output.nCols());
684 output[j][i]=theValue;
687 progress/=output.nRows();
688 pfnProgress(progress,pszMessage,pProgressArg);
690 gsl_rng_free(rangen);
693 template<
class T>
unsigned long int Filter2d::morphology(
const Vector2d<T>& input,
Vector2d<T>& output,
const std::string& method,
int dimX,
int dimY,
bool disc,
double hThreshold)
695 const char* pszMessage;
696 void* pProgressArg=NULL;
697 GDALProgressFunc pfnProgress=GDALTermProgress;
699 pfnProgress(progress,pszMessage,pProgressArg);
701 double noDataValue=0;
702 if(m_noDataValues.size())
703 noDataValue=m_noDataValues[0];
705 unsigned long int nchange=0;
711 output.resize(input.nRows(),input.nCols());
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];
720 for(
int y=0;y<input.nRows();++y){
723 inBuffer.erase(inBuffer.begin());
725 if(y+dimY/2<input.nRows()){
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];
732 int over=y+dimY/2-input.nRows();
733 int index=(inBuffer.size()-1)-over;
735 assert(index<inBuffer.size());
736 inBuffer.push_back(inBuffer[index]);
739 for(
int x=0;x<input.nCols();++x){
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]){
750 output[y][x]=currentValue;
752 output[y][x]=currentValue;
755 for(
int j=-(dimY-1)/2;j<=dimY/2;++j){
756 for(
int i=-(dimX-1)/2;i<=dimX/2;++i){
758 if(disc&&(d2>(dimX/2)*(dimY/2)))
764 else if(indexI>=input.nCols())
765 indexI=input.nCols()-i;
768 else if(y+j>=input.nRows())
769 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
772 if(inBuffer[indexJ][indexI]==noDataValue)
775 for(
int imask=0;imask<m_noDataValues.size();++imask){
776 if(inBuffer[indexJ][indexI]==m_noDataValues[imask]){
783 for(
int iclass=0;iclass<m_class.size();++iclass){
784 if(inBuffer[indexJ][indexI]==m_class[iclass]){
790 statBuffer.push_back(binValue);
792 statBuffer.push_back(inBuffer[indexJ][indexI]);
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);
804 case(filter2d::erode):
805 if(output[y][x]>stat.mymin(statBuffer)+hThreshold){
806 output[y][x]=stat.mymin(statBuffer);
811 std::ostringstream ess;
812 ess <<
"Error: morphology method " << method <<
" not supported, choose " << filter2d::dilate <<
" (dilate) or " << filter2d::erode <<
" (erode)" << std::endl;
816 if(output[y][x]&&m_class.size())
817 output[y][x]=m_class[0];
824 output[y][x]=noDataValue;
828 progress/=output.nRows();
829 pfnProgress(progress,pszMessage,pProgressArg);
834 template<
class T>
unsigned long int Filter2d::dsm2dtm_nwse(
const Vector2d<T>& inputDSM,
Vector2d<T>& outputMask,
double hThreshold,
int nlimit,
int dim)
836 const char* pszMessage;
837 void* pProgressArg=NULL;
838 GDALProgressFunc pfnProgress=GDALTermProgress;
840 pfnProgress(progress,pszMessage,pProgressArg);
843 double noDataValue=0;
844 if(m_noDataValues.size())
845 noDataValue=m_noDataValues[0];
847 unsigned long int nchange=0;
854 if(outputMask.size()!=inputDSM.nRows())
855 outputMask.resize(inputDSM.nRows());
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];
864 for(
int y=0;y<tmpDSM.nRows();++y){
867 inBuffer.erase(inBuffer.begin());
869 if(y+dimY/2<tmpDSM.nRows()){
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];
876 int over=y+dimY/2-tmpDSM.nRows();
877 int index=(inBuffer.size()-1)-over;
879 assert(index<inBuffer.size());
880 inBuffer.push_back(inBuffer[index]);
883 for(
int x=0;x<tmpDSM.nCols();++x){
884 double centerValue=inBuffer[(dimY-1)/2][x];
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){
893 else if(indexI>=tmpDSM.nCols())
894 indexI=tmpDSM.nCols()-i;
897 else if(y+j>=tmpDSM.nRows())
898 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
901 double difference=(centerValue-inBuffer[indexJ][indexI]);
903 neighbors.push_back(inBuffer[indexJ][indexI]);
904 if(difference>hThreshold)
915 sort(neighbors.begin(),neighbors.end());
916 assert(neighbors.size()>1);
917 inBuffer[(dimY-1)/2][x]=neighbors[1];
922 progress/=outputMask.nRows();
923 pfnProgress(progress,pszMessage,pProgressArg);
928 template<
class T>
unsigned long int Filter2d::dsm2dtm_nesw(
const Vector2d<T>& inputDSM,
Vector2d<T>& outputMask,
double hThreshold,
int nlimit,
int dim)
930 const char* pszMessage;
931 void* pProgressArg=NULL;
932 GDALProgressFunc pfnProgress=GDALTermProgress;
934 pfnProgress(progress,pszMessage,pProgressArg);
937 double noDataValue=0;
938 if(m_noDataValues.size())
939 noDataValue=m_noDataValues[0];
941 unsigned long int nchange=0;
948 if(outputMask.size()!=inputDSM.nRows())
949 outputMask.resize(inputDSM.nRows());
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];
958 for(
int y=0;y<tmpDSM.nRows();++y){
961 inBuffer.erase(inBuffer.begin());
963 if(y+dimY/2<tmpDSM.nRows()){
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];
970 int over=y+dimY/2-tmpDSM.nRows();
971 int index=(inBuffer.size()-1)-over;
973 assert(index<inBuffer.size());
974 inBuffer.push_back(inBuffer[index]);
977 for(
int x=tmpDSM.nCols()-1;x>=0;--x){
978 double centerValue=inBuffer[(dimY-1)/2][x];
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){
987 else if(indexI>=tmpDSM.nCols())
988 indexI=tmpDSM.nCols()-i;
991 else if(y+j>=tmpDSM.nRows())
992 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
995 double difference=(centerValue-inBuffer[indexJ][indexI]);
997 neighbors.push_back(inBuffer[indexJ][indexI]);
998 if(difference>hThreshold)
1002 if(nmasked<=nlimit){
1009 sort(neighbors.begin(),neighbors.end());
1010 assert(neighbors.size()>1);
1011 inBuffer[(dimY-1)/2][x]=neighbors[1];
1016 progress/=outputMask.nRows();
1017 pfnProgress(progress,pszMessage,pProgressArg);
1022 template<
class T>
unsigned long int Filter2d::dsm2dtm_senw(
const Vector2d<T>& inputDSM,
Vector2d<T>& outputMask,
double hThreshold,
int nlimit,
int dim)
1024 const char* pszMessage;
1025 void* pProgressArg=NULL;
1026 GDALProgressFunc pfnProgress=GDALTermProgress;
1028 pfnProgress(progress,pszMessage,pProgressArg);
1031 double noDataValue=0;
1032 if(m_noDataValues.size())
1033 noDataValue=m_noDataValues[0];
1035 unsigned long int nchange=0;
1042 if(outputMask.size()!=inputDSM.nRows())
1043 outputMask.resize(inputDSM.nRows());
1045 int indexJ=inputDSM.nRows()-1;
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];
1052 for(
int y=tmpDSM.nRows()-1;y>=0;--y){
1053 if(y<tmpDSM.nRows()-1){
1055 inBuffer.erase(inBuffer.end()-1);
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];
1064 inBuffer.insert(inBuffer.begin(),inBuffer[abs(y-dimY/2)]);
1067 for(
int x=tmpDSM.nCols()-1;x>=0;--x){
1068 double centerValue=inBuffer[(dimY-1)/2][x];
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){
1077 else if(indexI>=tmpDSM.nCols())
1078 indexI=tmpDSM.nCols()-i;
1081 else if(y+j>=tmpDSM.nRows())
1082 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
1084 indexJ=(dimY-1)/2+j;
1085 double difference=(centerValue-inBuffer[indexJ][indexI]);
1087 neighbors.push_back(inBuffer[indexJ][indexI]);
1088 if(difference>hThreshold)
1092 if(nmasked<=nlimit){
1099 sort(neighbors.begin(),neighbors.end());
1100 assert(neighbors.size()>1);
1101 inBuffer[(dimY-1)/2][x]=neighbors[1];
1106 progress/=outputMask.nRows();
1107 pfnProgress(progress,pszMessage,pProgressArg);
1112 template<
class T>
unsigned long int Filter2d::dsm2dtm_swne(
const Vector2d<T>& inputDSM,
Vector2d<T>& outputMask,
double hThreshold,
int nlimit,
int dim)
1114 const char* pszMessage;
1115 void* pProgressArg=NULL;
1116 GDALProgressFunc pfnProgress=GDALTermProgress;
1118 pfnProgress(progress,pszMessage,pProgressArg);
1121 double noDataValue=0;
1122 if(m_noDataValues.size())
1123 noDataValue=m_noDataValues[0];
1125 unsigned long int nchange=0;
1132 if(outputMask.size()!=inputDSM.nRows())
1133 outputMask.resize(inputDSM.nRows());
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];
1142 for(
int y=tmpDSM.nRows()-1;y>=0;--y){
1143 if(y<tmpDSM.nRows()-1){
1145 inBuffer.erase(inBuffer.end()-1);
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];
1154 inBuffer.insert(inBuffer.begin(),inBuffer[abs(y-dimY/2)]);
1157 for(
int x=0;x<tmpDSM.nCols();++x){
1158 double centerValue=inBuffer[(dimY-1)/2][x];
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){
1167 else if(indexI>=tmpDSM.nCols())
1168 indexI=tmpDSM.nCols()-i;
1171 else if(y+j>=tmpDSM.nRows())
1172 indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
1174 indexJ=(dimY-1)/2+j;
1175 double difference=(centerValue-inBuffer[indexJ][indexI]);
1177 neighbors.push_back(inBuffer[indexJ][indexI]);
1178 if(difference>hThreshold)
1182 if(nmasked<=nlimit){
1189 sort(neighbors.begin(),neighbors.end());
1190 assert(neighbors.size()>1);
1191 inBuffer[(dimY-1)/2][x]=neighbors[1];
1196 progress/=outputMask.nRows();
1197 pfnProgress(progress,pszMessage,pProgressArg);
1202 template<
class T>
void Filter2d::shadowDsm(
const Vector2d<T>& input,
Vector2d<T>& output,
double sza,
double saa,
double pixelSize,
short shadowFlag)
1204 unsigned int ncols=input.nCols();
1206 output.resize(input.nRows(),ncols);
1213 const char* pszMessage;
1214 void* pProgressArg=NULL;
1215 GDALProgressFunc pfnProgress=GDALTermProgress;
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)));
1222 double theDir=DEG2RAD(saa)+PI/2.0;
1225 for(
int d=0;d<theDist;++d){
1226 indexI=x+d*cos(theDir);
1227 indexJ=y+d*sin(theDir);
1228 if(indexJ<0||indexJ>=input.size())
1230 if(indexI<0||indexI>=input[indexJ].size())
1232 if(input[indexJ][indexI]<currentValue-d*pixelSize/tan(DEG2RAD(sza))){
1233 output[indexJ][indexI]=shadowFlag;
1238 progress/=output.nRows();
1239 pfnProgress(progress,pszMessage,pProgressArg);
1243 template<
class T>
void Filter2d::dwtForward(
Vector2d<T>& theBuffer,
const std::string& wavelet_type,
int family){
1244 const char* pszMessage;
1245 void* pProgressArg=NULL;
1246 GDALProgressFunc pfnProgress=GDALTermProgress;
1248 pfnProgress(progress,pszMessage,pProgressArg);
1250 int nRow=theBuffer.size();
1252 int nCol=theBuffer[0].size();
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];
1268 int nsize=theBuffer.size()*theBuffer[0].size();
1270 gsl_wavelet_workspace *work;
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];
1282 progress=(1.0+irow);
1283 progress/=theBuffer.nRows();
1284 pfnProgress(progress,pszMessage,pProgressArg);
1286 gsl_wavelet_free (w);
1287 gsl_wavelet_workspace_free (work);
1290 template<
class T>
void Filter2d::dwtInverse(
Vector2d<T>& theBuffer,
const std::string& wavelet_type,
int family){
1291 const char* pszMessage;
1292 void* pProgressArg=NULL;
1293 GDALProgressFunc pfnProgress=GDALTermProgress;
1295 pfnProgress(progress,pszMessage,pProgressArg);
1297 int nRow=theBuffer.size();
1299 int nCol=theBuffer[0].size();
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]);
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];
1316 int nsize=theBuffer.size()*theBuffer[0].size();
1318 gsl_wavelet_workspace *work;
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];
1330 progress=(1.0+irow);
1331 progress/=theBuffer.nRows();
1332 pfnProgress(progress,pszMessage,pProgressArg);
1334 gsl_wavelet_free (w);
1335 gsl_wavelet_workspace_free (work);
1338 template<
class T>
void Filter2d::dwtCut(
Vector2d<T>& theBuffer,
const std::string& wavelet_type,
int family,
double cut){
1339 const char* pszMessage;
1340 void* pProgressArg=NULL;
1341 GDALProgressFunc pfnProgress=GDALTermProgress;
1343 pfnProgress(progress,pszMessage,pProgressArg);
1345 int nRow=theBuffer.size();
1347 int nCol=theBuffer[0].size();
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];
1365 int nsize=theBuffer.size()*theBuffer[0].size();
1367 gsl_wavelet_workspace *work;
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]);
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++)
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];
1388 progress=(1.0+irow);
1389 progress/=theBuffer.nRows();
1390 pfnProgress(progress,pszMessage,pProgressArg);
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());
1398 gsl_wavelet_free (w);
1399 gsl_wavelet_workspace_free (work);