pktools 2.6.7
Processing Kernel for geospatial data
pkfilterdem.cc
1/**********************************************************************
2pkfilterdem.cc: Filter digital elevation model raster datasets
3Copyright (C) 2008-2014 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#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/******************************************************************************/
72using namespace std;
73/*------------------
74 Main procedure
75 ----------------*/
76int 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 omit 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 if(verbose_opt[0]){
420 std::cout << "writing data block" << std::endl;
421 std::cout << "outputData.nCols(): " << outputData.nCols() << std::endl;
422 std::cout << "outputData.nRows(): " << outputData.nRows() << std::endl;
423 std::cout << "outputWriter.nrOfCol(): " << outputWriter.nrOfCol() << std::endl;
424 std::cout << "outputWriter.nrOfRow(): " << outputWriter.nrOfRow() << std::endl;
425 std::cout << "outputWriter.getDataType(): " << outputWriter.getDataType() << std::endl;
426 }
427 //write outputData to outputWriter
428 try{
429 outputWriter.writeDataBlock(outputData,0,outputData.nCols()-1,0,outputData.nRows()-1);
430 }
431 catch(std::string errorString){
432 cout << errorString << endl;
433 exit(1);
434 }
435
436 // progress=1;
437 // pfnProgress(progress,pszMessage,pProgressArg);
438 input.close();
439 outputWriter.close();
440 return 0;
441}
int nrOfRow(void) const
Get the number of rows of this dataset.
std::string getProjection(void) const
Get the projection string (deprecated, use getProjectionRef instead)
int nrOfBand(void) const
Get the number of bands of this dataset.
void copyGeoTransform(const ImgRasterGdal &imgSrc)
Copy geotransform information from another georeferenced image.
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
double getDeltaX(void) const
Get the pixel cell spacing in x.
std::string getInterleave() const
Get the band coding (interleave)
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...
GDALDataType getDataType(int band=0) const
Get the GDAL datatype for this dataset.
CPLErr setProjection(const std::string &projection)
Set the projection for this dataset in well known text (wkt) format.
void close(void)
Set the memory (in MB) to cache a number of rows in memory.
void open(const std::string &filename, const GDALAccess &readMode=GA_ReadOnly)
Open an image.
void readDataBlock(Vector2d< T > &buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band=0)
Read pixel cell values for a range of columns and rows for a specific band (all indices start countin...
void close(void)
Close the image.
void open(const std::string &filename, const ImgReaderGdal &imgSrc, const std::vector< std::string > &options=std::vector< std::string >())
Open an image for writing, copying image attributes from a source image. Image is directly written to...
bool writeDataBlock(Vector2d< T > &buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band=0)
Write pixel cell values for a range of columns and rows for a specific band (all indices start counti...
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)