pktools  2.6.7
Processing Kernel for geospatial data
Filter2d.h
1 /**********************************************************************
2 Filter2d.h: class for filtering images
3 Copyright (C) 2008-2012 Pieter Kempeneers
4 
5 This file is part of pktools
6 
7 pktools is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 pktools is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along 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>
46 extern "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 
59 namespace 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 
65 class Filter2d
66 {
67 public:
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 
131 private:
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 
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)
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 
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;
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 
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)
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 
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;
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 
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;
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 
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;
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_ */
Definition: Filter.h:33