pktools  2.6.7
Processing Kernel for geospatial data
pkfilterdem.cc
1 /**********************************************************************
2 pkfilterdem.cc: Filter digital elevation model raster datasets
3 Copyright (C) 2008-2014 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 #include <assert.h>
21 #include <iostream>
22 #include <string>
23 #include "base/Optionpk.h"
24 #include "base/Vector2d.h"
25 #include "algorithms/Filter2d.h"
26 #include "imageclasses/ImgReaderGdal.h"
27 #include "imageclasses/ImgWriterGdal.h"
28 
29 /******************************************************************************/
72 using namespace std;
73 /*------------------
74  Main procedure
75  ----------------*/
76 int main(int argc,char **argv) {
77  Optionpk<std::string> input_opt("i","input","input image file");
78  Optionpk<std::string> output_opt("o", "output", "Output image file");
79  // Optionpk<std::string> tmpdir_opt("tmp", "tmp", "Temporary directory","/tmp",2);
80  Optionpk<bool> disc_opt("circ", "circular", "circular disc kernel for dilation and erosion", false);
81  Optionpk<string> postFilter_opt("f", "filter", "post processing filter: vito, etew_min, promorph (progressive morphological filter),open,close).");
82  Optionpk<double> dim_opt("dim", "dim", "maximum filter kernel size", 17);
83  Optionpk<double> maxSlope_opt("st", "st", "slope threshold used for morphological filtering. Use a low values to remove more height objects in flat terrains", 0.0);
84  Optionpk<double> hThreshold_opt("ht", "ht", "initial height threshold for progressive morphological filtering. Use low values to remove more height objects. Optionally, a maximum height threshold can be set via a second argument (e.g., -ht 0.2 -ht 2.5 sets an initial threshold at 0.2 m and caps the threshold at 2.5 m).", 0.2);
85  Optionpk<short> minChange_opt("minchange", "minchange", "Stop iterations when no more pixels are changed than this threshold.", 0);
86  Optionpk<std::string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image","");
87  Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate).","GTiff");
88  Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid). Use none to ommit color table");
89  Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
90  Optionpk<short> nodata_opt("nodata", "nodata", "nodata value");
91  Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0,2);
92 
93  disc_opt.setHide(1);
94  maxSlope_opt.setHide(1);
95  hThreshold_opt.setHide(1);
96  minChange_opt.setHide(1);
97  otype_opt.setHide(1);
98  oformat_opt.setHide(1);
99  colorTable_opt.setHide(1);
100  nodata_opt.setHide(1);
101 
102  bool doProcess;//stop process when program was invoked with help option (-h --help)
103  try{
104  doProcess=input_opt.retrieveOption(argc,argv);
105  output_opt.retrieveOption(argc,argv);
106  // tmpdir_opt.retrieveOption(argc,argv);
107  postFilter_opt.retrieveOption(argc,argv);
108  dim_opt.retrieveOption(argc,argv);
109  disc_opt.retrieveOption(argc,argv);
110  maxSlope_opt.retrieveOption(argc,argv);
111  hThreshold_opt.retrieveOption(argc,argv);
112  minChange_opt.retrieveOption(argc,argv);
113  otype_opt.retrieveOption(argc,argv);
114  oformat_opt.retrieveOption(argc,argv);
115  colorTable_opt.retrieveOption(argc,argv);
116  nodata_opt.retrieveOption(argc,argv);
117  verbose_opt.retrieveOption(argc,argv);
118  }
119  catch(string predefinedString){
120  std::cout << predefinedString << std::endl;
121  exit(0);
122  }
123  if(!doProcess){
124  cout << endl;
125  cout << "Usage: pkfilterdem -i input.txt -o output" << endl;
126  cout << endl;
127  std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
128  exit(0);//help was invoked, stop processing
129  }
130 
131  ImgReaderGdal input;
132  ImgWriterGdal outputWriter;
133  if(input_opt.empty()){
134  cerr << "Error: no input file selected, use option -i" << endl;
135  exit(1);
136  }
137  if(output_opt.empty()){
138  cerr << "Error: no outputWriter file selected, use option -o" << endl;
139  exit(1);
140  }
141  if(postFilter_opt.empty()){
142  cerr << "Error: no filter selected, use option -f" << endl;
143  exit(1);
144  }
145  input.open(input_opt[0]);
146  GDALDataType theType=GDT_Unknown;
147  if(verbose_opt[0])
148  cout << "possible output data types: ";
149  for(int iType = 0; iType < GDT_TypeCount; ++iType){
150  if(verbose_opt[0])
151  cout << " " << GDALGetDataTypeName((GDALDataType)iType);
152  if( GDALGetDataTypeName((GDALDataType)iType) != NULL
153  && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
154  otype_opt[0].c_str()))
155  theType=(GDALDataType) iType;
156  }
157  if(theType==GDT_Unknown)
158  theType=input.getDataType();
159 
160  if(verbose_opt[0])
161  std::cout << std::endl << "Output pixel type: " << GDALGetDataTypeName(theType) << endl;
162 
163  string imageType;//=input.getImageType();
164  if(oformat_opt.size())
165  imageType=oformat_opt[0];
166 
167  if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
168  string theInterleave="INTERLEAVE=";
169  theInterleave+=input.getInterleave();
170  option_opt.push_back(theInterleave);
171  }
172 
173  if(verbose_opt[0])
174  cout << "opening output file " << output_opt[0] << endl;
175  outputWriter.open(output_opt[0],input.nrOfCol(),input.nrOfRow(),1,theType,imageType,option_opt);
176  //set projection
177  outputWriter.setProjection(input.getProjection());
178  outputWriter.copyGeoTransform(input);
179  if(colorTable_opt.size())
180  outputWriter.setColorTable(colorTable_opt[0]);
181 
182  //set nodata value
183  if(nodata_opt.size()){
184  for(int iband=0;iband<outputWriter.nrOfBand();++iband)
185  outputWriter.GDALSetNoDataValue(nodata_opt[0],iband);
186  }
187 
188  Vector2d<double> inputData(input.nrOfRow(),input.nrOfCol());
189  Vector2d<double> outputData(outputWriter.nrOfRow(),outputWriter.nrOfCol());
190  Vector2d<double> tmpData(outputWriter.nrOfRow(),outputWriter.nrOfCol());
191  input.readDataBlock(inputData,0,inputData.nCols()-1,0,inputData.nRows()-1);
192 
193  //apply post filter
194  std::cout << "Applying post processing filter: " << postFilter_opt[0] << std::endl;
195 
196  // const char* pszMessage;
197  // void* pProgressArg=NULL;
198  // GDALProgressFunc pfnProgress=GDALTermProgress;
199  // double progress=0;
200  // pfnProgress(progress,pszMessage,pProgressArg);
201 
202  //make sure dim_opt contains initial [0] and maximum [1] kernel sizes in this order
203  if(dim_opt.size()<2)
204  dim_opt.insert(dim_opt.begin(),3);
205  if(dim_opt[0]>dim_opt[1]){
206  dim_opt.insert(dim_opt.begin(),dim_opt[1]);
207  dim_opt.erase(dim_opt.begin()+2);
208  }
209 
210  filter2d::Filter2d theFilter;
211  if(nodata_opt.size()){
212  for(int inodata=0;inodata<nodata_opt.size();++inodata)
213  theFilter.pushNoDataValue(nodata_opt[inodata]);
214  }
215 
216  unsigned long int nchange=1;
217  if(postFilter_opt[0]=="vito"){
218  //todo: fill empty pixels
219  // hThreshold_opt.resize(4);
220  // hThreshold_opt[0]=0.7;
221  // hThreshold_opt[1]=0.3;
222  // hThreshold_opt[2]=0.1;
223  // hThreshold_opt[2]=-0.2;
224  vector<int> nlimit(4);
225  nlimit[0]=2;
226  nlimit[1]=3;
227  nlimit[2]=4;
228  nlimit[2]=2;
229  //init finalMask
230  for(int irow=0;irow<tmpData.nRows();++irow)
231  for(int icol=0;icol<tmpData.nCols();++icol)
232  tmpData[irow][icol]=1;
233  for(int iheight=0;iheight<hThreshold_opt.size();++iheight){
234  if(verbose_opt[0])
235  cout << "height: " << hThreshold_opt[iheight] << endl;
236  //todo:replace with binary mask (or short) -> adapt template with T1,T2 in Filter2d
237  Vector2d<double> tmpMask(input.nrOfRow(),input.nrOfCol());
238  for(int irow=0;irow<tmpMask.nRows();++irow)
239  for(int icol=0;icol<tmpMask.nCols();++icol)
240  tmpMask[irow][icol]=1;//1=surface, 0=terrain
241  if(verbose_opt[0])
242  cout << "filtering NWSE" << endl;
243  //from here
244  // Vector2d<double> tmpDSM(inputData);
245  // int dimX=dim_opt[0];
246  // int dimY=dim_opt[0];
247  // assert(dimX);
248  // assert(dimY);
249  // statfactory::StatFactory stat;
250  // Vector2d<double> inBuffer(dimY,inputData.nCols());
251  // if(tmpData.size()!=inputData.nRows())
252  // tmpData.resize(inputData.nRows());
253  // int indexI=0;
254  // int indexJ=inputData.nRows()-1;
255  // // int indexJ=0;
256  // //initialize last half of inBuffer
257  // for(int j=-(dimY-1)/2;j<=dimY/2;++j){
258  // for(int i=0;i<inputData.nCols();++i)
259  // inBuffer[indexJ][i]=tmpDSM[abs(j)][i];
260  // --indexJ;
261  // // ++indexJ;
262  // }
263  // for(int y=tmpDSM.nRows()-1;y>=0;--y){
264  // if(y){//inBuffer already initialized for y=0
265  // //erase first line from inBuffer
266  // inBuffer.erase(inBuffer.end()-1);
267  // // inBuffer.erase(inBuffer.begin());
268  // //read extra line and push back to inBuffer if not out of bounds
269  // if(y+dimY/2<tmpDSM.nRows()){
270  // //allocate buffer
271  // // inBuffer.push_back(inBuffer.back());
272  // inBuffer.insert(inBuffer.begin(),*(inBuffer.begin()));
273  // for(int i=0;i<tmpDSM.nCols();++i)
274  // inBuffer[0][i]=tmpDSM[y-dimY/2][i];
275  // // inBuffer[inBuffer.size()-1][i]=tmpDSM[y+dimY/2][i];
276  // }
277  // else{
278  // int over=y+dimY/2-tmpDSM.nRows();
279  // int index=(inBuffer.size()-1)-over;
280  // assert(index>=0);
281  // assert(index<inBuffer.size());
282  // inBuffer.push_back(inBuffer[index]);
283  // }
284  // }
285  // for(int x=tmpDSM.nCols()-1;x>=0;--x){
286  // double centerValue=inBuffer[(dimY-1)/2][x];
287  // //test
288  // cout << "pixel (" << x << "," << y << "): " << centerValue << endl;
289  // short nmasked=0;
290  // std::vector<double> neighbors;
291  // for(int j=-(dimY-1)/2;j<=dimY/2;++j){
292  // for(int i=-(dimX-1)/2;i<=dimX/2;++i){
293  // indexI=x+i;
294  // //check if out of bounds
295  // if(indexI<0)
296  // indexI=-indexI;
297  // else if(indexI>=tmpDSM.nCols())
298  // indexI=tmpDSM.nCols()-i;
299  // if(y+j<0)
300  // indexJ=-j;
301  // else if(y+j>=tmpDSM.nRows())
302  // indexJ=(dimY>2) ? (dimY-1)/2-j : 0;
303  // else
304  // indexJ=(dimY-1)/2+j;
305  // double difference=(centerValue-inBuffer[indexJ][indexI]);
306  // //test
307  // cout << "centerValue-inBuffer[" << indexJ << "][" << indexI << "]=" << centerValue << " - " << inBuffer[indexJ][indexI] << " = " << difference << endl;
308  // if(i||j)//skip centerValue
309  // neighbors.push_back(inBuffer[indexJ][indexI]);
310  // if(difference>hThreshold_opt[iheight])
311  // ++nmasked;
312  // }
313  // }
314  // //test
315  // cout << "pixel " << x << ", " << y << ": nmasked is " << nmasked << endl;
316  // if(nmasked<=nlimit[iheight]){
317  // ++nchange;
318  // //reset pixel in outputMask
319  // tmpData[y][x]=0;
320  // //test
321  // cout << "pixel " << x << ", " << y << " is ground" << endl;
322  // }
323  // else{
324  // //reset pixel height in tmpDSM
325  // sort(neighbors.begin(),neighbors.end());
326  // assert(neighbors.size()>1);
327  // inBuffer[(dimY-1)/2][x]=neighbors[1];
328  // //test
329  // cout << "pixel " << x << ", " << y << " is surface, reset DSM to " << neighbors[1] << endl;
330  // /* inBuffer[(dimY-1)/2][x]=stat.mymin(neighbors); */
331  // }
332  // }
333  // }
334  //to here
335 
336  theFilter.dsm2dtm_nwse(inputData,tmpData,hThreshold_opt[iheight],nlimit[iheight],dim_opt[0]);
337  if(verbose_opt[0])
338  cout << "filtering NESW" << endl;
339  theFilter.dsm2dtm_nesw(inputData,tmpData,hThreshold_opt[iheight],nlimit[iheight],dim_opt[0]);
340  if(verbose_opt[0])
341  cout << "filtering SENW" << endl;
342  theFilter.dsm2dtm_senw(inputData,tmpData,hThreshold_opt[iheight],nlimit[iheight],dim_opt[0]);
343  if(verbose_opt[0])
344  cout << "filtering SWNE" << endl;
345  theFilter.dsm2dtm_swne(inputData,tmpData,hThreshold_opt[iheight],nlimit[iheight],dim_opt[0]);
346  }
347  outputData=tmpData;
348  //todo: interpolate
349  //outputData.setMask(tmpData,1,0);
350  }
351  else if(postFilter_opt[0]=="etew_min"){
352  //Elevation Threshold with Expand Window (ETEW) Filter (p.73 from Airborne LIDAR Data Processing and Analysis Tools ALDPAT 1.0)
353  //first iteration is performed assuming only minima are selected using options -fir all -comp min
354  //increase cells and thresholds until no points from the previous iteration are discarded.
355  int dim=dim_opt[0];
356  // theFilter.setNoValue(0);
357  int iteration=1;
358  while(nchange>minChange_opt[0]&&dim<=dim_opt[1]){
359  double hThreshold=maxSlope_opt[0]*dim;
360  nchange=theFilter.morphology(inputData,outputData,"erode",dim,dim,disc_opt[0],hThreshold);
361  inputData=outputData;
362  dim+=2;//change from theory: originally double cellCize
363  std::cout << "iteration " << iteration << ": " << nchange << " pixels changed" << std::endl;
364  ++iteration;
365  }
366  }
367  else if(postFilter_opt[0]=="promorph"){
368  //Progressive morphological filter tgrs2003_zhang vol41 pp 872-882
369  //first iteration is performed assuming only minima are selected using options -fir all -comp min
370  //increase cells and thresholds until no points from the previous iteration are discarded.
371  int dim=dim_opt[0];
372  double hThreshold=hThreshold_opt[0];
373  int iteration=1;
374  while(nchange>minChange_opt[0]&&dim<=dim_opt[1]){
375  std::cout << "iteration " << iteration << " with window size " << dim << " and dh_max: " << hThreshold << std::endl;
376  try{
377  nchange=theFilter.morphology(inputData,outputData,"erode",dim,dim,disc_opt[0],hThreshold);
378  theFilter.morphology(outputData,inputData,"dilate",dim,dim,disc_opt[0],hThreshold);
379  theFilter.doit(inputData,outputData,"median",dim,dim,1,disc_opt[0]);
380  inputData=outputData;
381  }
382  catch(std::string errorString){
383  cout << errorString << endl;
384  exit(1);
385  }
386  int newdim=(dim==1)? 3: 2*(dim-1)+1;
387  hThreshold=hThreshold_opt[0]+maxSlope_opt[0]*(newdim-dim)*input.getDeltaX();
388  dim=newdim;
389  if(hThreshold_opt.size()>1){
390  if(hThreshold>hThreshold_opt[1]){
391  hThreshold=hThreshold_opt[1];
392  }
393  }
394  std::cout << "iteration " << iteration << ": " << nchange << " pixels changed" << std::endl;
395  ++iteration;
396  }
397  }
398  else if(postFilter_opt[0]=="open"){
399  try{
400  theFilter.morphology(inputData,tmpData,"erode",dim_opt[0],dim_opt[0],disc_opt[0],hThreshold_opt[0]);
401  theFilter.morphology(tmpData,outputData,"dilate",dim_opt[0],dim_opt[0],disc_opt[0],hThreshold_opt[0]);
402  outputData=inputData;
403  }
404  catch(std::string errorString){
405  cout << errorString << endl;
406  exit(1);
407  }
408  }
409  else if(postFilter_opt[0]=="close"){
410  try{
411  theFilter.morphology(inputData,tmpData,"dilate",dim_opt[0],dim_opt[0],disc_opt[0],hThreshold_opt[0]);
412  theFilter.morphology(tmpData,outputData,"erode",dim_opt[0],dim_opt[0],disc_opt[0],hThreshold_opt[0]);
413  }
414  catch(std::string errorString){
415  cout << errorString << endl;
416  exit(1);
417  }
418  }
419  //write outputData to outputWriter
420  outputWriter.writeDataBlock(outputData,GDT_Float64,0,outputData.nCols()-1,0,outputData.nRows()-1);
421 
422  // progress=1;
423  // pfnProgress(progress,pszMessage,pProgressArg);
424  input.close();
425  outputWriter.close();
426  return 0;
427 }
int nrOfBand(void) const
Get the number of bands of this dataset.
STL namespace.
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
int nrOfRow(void) const
Get the number of rows of this dataset.
void setColorTable(const std::string &filename, int band=0)
Set the color table using an (ASCII) file with 5 columns (value R G B alpha)
void copyGeoTransform(const ImgRasterGdal &imgSrc)
Copy geotransform information from another georeferenced image.
bool writeDataBlock(Vector2d< T > &buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band=0)
Write pixel cell values for a range of columns and rows for a specific band (all indices start counti...
void open(const std::string &filename, const ImgReaderGdal &imgSrc, const std::vector< std::string > &options=std::vector< std::string >())
Open an image for writing, copying image attributes from a source image. Image is directly written to...
double getDeltaX(void) const
Get the pixel cell spacing in x.
CPLErr setProjection(const std::string &projection)
Set the projection for this dataset in well known text (wkt) format.
std::string getInterleave() const
Get the band coding (interleave)
void close(void)
Close the image.
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
void readDataBlock(Vector2d< T > &buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band=0)
Read pixel cell values for a range of columns and rows for a specific band (all indices start countin...
CPLErr GDALSetNoDataValue(double noDataValue, int band=0)
Set the GDAL (internal) no data value for this data set. Only a single no data value per band is supp...
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98