pktools 2.6.7
Processing Kernel for geospatial data
Filter2d.h
1/**********************************************************************
2Filter2d.h: class for filtering images
3Copyright (C) 2008-2012 Pieter Kempeneers
4
5This file is part of pktools
6
7pktools is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published by
9the Free Software Foundation, either version 3 of the License, or
10(at your option) any later version.
11
12pktools is distributed in the hope that it will be useful,
13but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15GNU General Public License for more details.
16
17You should have received a copy of the GNU General Public License
18along with pktools. If not, see <http://www.gnu.org/licenses/>.
19***********************************************************************/
20#ifndef _MYFILTER2D_H_
21#define _MYFILTER2D_H_
22
23#ifndef PI
24#define PI 3.1415926535897932384626433832795
25#endif
26
27#ifndef DEG2RAD
28#define DEG2RAD(DEG) (DEG/180.0*PI)
29#endif
30
31#ifndef RAD2DEG
32#define RAD2DEG(RAD) (RAD/PI*180)
33#endif
34
35#ifdef WIN32
36#include <process.h>
37#define getpid _getpid
38#endif
39
40#include <assert.h>
41#include <math.h>
42#include <limits>
43#include <vector>
44#include <string>
45#include <map>
46extern "C" {
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>
52}
53#include "base/Vector2d.h"
54#include "Filter.h"
55#include "imageclasses/ImgReaderGdal.h"
56#include "imageclasses/ImgWriterGdal.h"
57#include "algorithms/StatFactory.h"
58
59namespace filter2d
60{
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};
62
63 enum RESAMPLE { NEAR = 0, BILINEAR = 1, BICUBIC = 2 };//bicubic not supported yet...
64
66{
67public:
68 Filter2d(void);
69 Filter2d(const Vector2d<double> &taps);
70 virtual ~Filter2d(){};
71 static FILTER_TYPE getFilterType(const std::string filterType){
72 std::map<std::string, FILTER_TYPE> m_filterMap;
73 initMap(m_filterMap);
74 return m_filterMap[filterType];
75 };
76 static const RESAMPLE getResampleType(const std::string resampleType){
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 };
86
87 void setTaps(const Vector2d<double> &taps);
88 /* void setNoValue(double noValue=0){m_noValue=noValue;}; */
89 void pushClass(short theClass=1){m_class.push_back(theClass);};
90 int pushNoDataValue(double noDataValue=0);//{m_mask.push_back(theMask);};
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;};
94 void filter(ImgReaderGdal& input, ImgWriterGdal& output, bool absolute=false, bool normalize=false, bool noData=false);
95 void smooth(ImgReaderGdal& input, ImgWriterGdal& output,int dim);
96 void smooth(ImgReaderGdal& input, ImgWriterGdal& output,int dimX, int dimY);
97 void smoothNoData(ImgReaderGdal& input, ImgWriterGdal& output,int dim);
98 void smoothNoData(ImgReaderGdal& input, ImgWriterGdal& output,int dimX, int dimY);
99 template<class T1, class T2> void filter(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector);
100 template<class T1, class T2> void smooth(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector,int dim);
101 template<class T1, class T2> void smooth(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector,int dimX, int dimY);
102 void dwtForward(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family);
103 void dwtInverse(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& wavelet_type, int family);
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>());
109 /* void homogeneousSpatial(const std::string& inputFilename, const std::string& outputFilename, int dim, bool disc=false, int noValue=0); */
110 void doit(ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down=1, bool disc=false);
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);
113 void mrf(ImgReaderGdal& input, ImgWriterGdal& output, int dimX, int dimY, Vector2d<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);
124 void shadowDsm(ImgReaderGdal& input, ImgWriterGdal& 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);
130
131private:
132 static void initMap(std::map<std::string, FILTER_TYPE>& m_filterMap){
133 //initialize selMap
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;
175 }
176
177 Vector2d<double> m_taps;
178 /* double m_noValue; */
179 std::vector<short> m_class;
180 /* std::vector<short> m_mask; */
181 std::vector<double> m_noDataValues;
182 std::vector<double> m_threshold;
183};
184
185
186 template<class T1, class T2> void Filter2d::smooth(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector,int dim)
187 {
188 smooth(inputVector,outputVector,dim,dim);
189 }
190
191 template<class T1, class T2> void Filter2d::smooth(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector,int dimX, int dimY)
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 }
201
202 template<class T1, class T2> void Filter2d::filter(const Vector2d<T1>& inputVector, Vector2d<T2>& outputVector)
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 }
258
259template<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)
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}
602
603// class Compare_mapValue{
604// public:
605// int operator() (const map<int,int>::value_type& v1, const map<int, int>::value_type& v2) const{
606// return (v1.second)>(v2.second);
607// }
608// };
609
610template<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;
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}
692
693template<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)
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}
833
834 template<class T> unsigned long int Filter2d::dsm2dtm_nwse(const Vector2d<T>& inputDSM, Vector2d<T>& outputMask, double hThreshold, int nlimit, int dim)
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}
927
928 template<class T> unsigned long int Filter2d::dsm2dtm_nesw(const Vector2d<T>& inputDSM, Vector2d<T>& outputMask, double hThreshold, int nlimit, int dim)
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}
1021
1022 template<class T> unsigned long int Filter2d::dsm2dtm_senw(const Vector2d<T>& inputDSM, Vector2d<T>& outputMask, double hThreshold, int nlimit, int dim)
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}
1111
1112 template<class T> unsigned long int Filter2d::dsm2dtm_swne(const Vector2d<T>& inputDSM, Vector2d<T>& outputMask, double hThreshold, int nlimit, int dim)
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}
1201
1202 template<class T> void Filter2d::shadowDsm(const Vector2d<T>& input, Vector2d<T>& output, double sza, double saa, double pixelSize, short shadowFlag)
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}
1242
1243template<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;
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}
1289
1290template<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;
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}
1337
1338template<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;
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}
1402
1403}
1404
1405#endif /* _MYFILTER_H_ */