pktools 2.6.7
Processing Kernel for geospatial data
pklas2img.cc
1/**********************************************************************
2pklas2img.cc: Rasterize LAS/LAZ point clouds with filtering/compositing options
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 <iostream>
21#include "base/Optionpk.h"
22#include "imageclasses/ImgReaderGdal.h"
23#include "imageclasses/ImgWriterGdal.h"
24#include "imageclasses/ImgReaderOgr.h"
25#include "lasclasses/FileReaderLas.h"
26#include "algorithms/StatFactory.h"
27#include "algorithms/Filter2d.h"
28
29/******************************************************************************/
83using namespace std;
84
85int main(int argc,char **argv) {
86 Optionpk<string> input_opt("i", "input", "Input las file");
87 Optionpk<string> attribute_opt("n", "name", "names of the point attribute to select: intensity, angle, return, nreturn, spacing, z", "z");
88 // Optionpk<bool> disc_opt("circ", "circular", "circular disc kernel for dilation and erosion", false);
89 // Optionpk<double> maxSlope_opt("s", "maxSlope", "Maximum slope used for morphological filtering", 0.0);
90 // Optionpk<double> hThreshold_opt("ht", "maxHeight", "initial and maximum height threshold for progressive morphological filtering (e.g., -ht 0.2 -ht 2.5)", 0.2);
91 // Optionpk<short> maxIter_opt("maxit", "maxit", "Maximum number of iterations in post filter", 5);
92 Optionpk<unsigned short> returns_opt("ret", "ret", "number(s) of returns to include");
93 Optionpk<unsigned short> classes_opt("class", "class", "classes to keep: 0 (created, never classified), 1 (unclassified), 2 (ground), 3 (low vegetation), 4 (medium vegetation), 5 (high vegetation), 6 (building), 7 (low point, noise), 8 (model key-point), 9 (water), 10 (reserved), 11 (reserved), 12 (overlap)");
94 Optionpk<string> composite_opt("comp", "comp", "composite for multiple points in cell (min, max, absmin, absmax, median, mean, sum, first, last, profile (percentile height values), percentile, number (point density)). Last: overwrite cells with latest point", "last");
95 Optionpk<string> filter_opt("fir", "filter", "filter las points (first,last,single,multiple,all).", "all");
96 Optionpk<short> angle_min_opt("angle_min", "angle_min", "Minimum scan angle to read points.");
97 Optionpk<short> angle_max_opt("angle_max", "angle_max", "Maximum scan angle to read points.");
98 // Optionpk<string> postFilter_opt("pf", "pfilter", "post processing filter (etew_min,promorph (progressive morphological filter),bunting (adapted promorph),open,close,none).", "none");
99 // Optionpk<short> dimx_opt("dimx", "dimx", "Dimension X of postFilter", 3);
100 // Optionpk<short> dimy_opt("dimy", "dimy", "Dimension Y of postFilter", 3);
101 Optionpk<string> output_opt("o", "output", "Output image file");
102 Optionpk<string> projection_opt("a_srs", "a_srs", "assign the projection for the output file in epsg code, e.g., epsg:3035 for European LAEA projection");
103 Optionpk<double> ulx_opt("ulx", "ulx", "Upper left x value bounding box (in geocoordinates if georef is true). 0 is read from input file", 0.0);
104 Optionpk<double> uly_opt("uly", "uly", "Upper left y value bounding box (in geocoordinates if georef is true). 0 is read from input file", 0.0);
105 Optionpk<double> lrx_opt("lrx", "lrx", "Lower right x value bounding box (in geocoordinates if georef is true). 0 is read from input file", 0.0);
106 Optionpk<double> lry_opt("lry", "lry", "Lower right y value bounding box (in geocoordinates if georef is true). 0 is read from input file", 0.0);
107 Optionpk<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", "Byte");
108 Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate).", "GTiff");
109 Optionpk<double> dx_opt("dx", "dx", "Output resolution in x (in meter)", 1.0);
110 Optionpk<double> dy_opt("dy", "dy", "Output resolution in y (in meter)", 1.0);
111 Optionpk<short> nbin_opt("nbin", "nbin", "Number of percentile bins for calculating percentile height value profile (=number of output bands)", 10.0);
112 Optionpk<double> percentile_opt("perc","percentile","Percentile value used for rule percentile",95);
113 Optionpk<short> nodata_opt("nodata", "nodata", "nodata value to put in image", 0);
114 Optionpk<string> option_opt("co", "co", "Creation option for output file. Multiple options can be specified.");
115 Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)");
116 Optionpk<short> verbose_opt("v", "verbose", "verbose mode", 0,2);
117
118 nbin_opt.setHide(1);
119 percentile_opt.setHide(1);
120 nodata_opt.setHide(1);
121 option_opt.setHide(1);
122 colorTable_opt.setHide(1);
123
124 bool doProcess;//stop process when program was invoked with help option (-h --help)
125 try{
126 doProcess=input_opt.retrieveOption(argc,argv);
127 attribute_opt.retrieveOption(argc,argv);
128 returns_opt.retrieveOption(argc,argv);
129 classes_opt.retrieveOption(argc,argv);
130 composite_opt.retrieveOption(argc,argv);
131 filter_opt.retrieveOption(argc,argv);
132 angle_min_opt.retrieveOption(argc,argv);
133 angle_max_opt.retrieveOption(argc,argv);
134 output_opt.retrieveOption(argc,argv);
135 projection_opt.retrieveOption(argc,argv);
136 ulx_opt.retrieveOption(argc,argv);
137 uly_opt.retrieveOption(argc,argv);
138 lrx_opt.retrieveOption(argc,argv);
139 lry_opt.retrieveOption(argc,argv);
140 otype_opt.retrieveOption(argc,argv);
141 oformat_opt.retrieveOption(argc,argv);
142 dx_opt.retrieveOption(argc,argv);
143 dy_opt.retrieveOption(argc,argv);
144 nbin_opt.retrieveOption(argc,argv);
145 percentile_opt.retrieveOption(argc,argv);
146 nodata_opt.retrieveOption(argc,argv);
147 option_opt.retrieveOption(argc,argv);
148 colorTable_opt.retrieveOption(argc,argv);
149 verbose_opt.retrieveOption(argc,argv);
150 }
151 catch(string predefinedString){
152 std::cout << predefinedString << std::endl;
153 exit(0);
154 }
155
156 if(!doProcess){
157 cout << endl;
158 cout << "pklas2img -i lasfile -o output" << endl;
159 cout << endl;
160 std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
161 exit(0);//help was invoked, stop processing
162 }
163 //todo: is this needed?
164 GDALAllRegister();
165
166 double dfComplete=0.0;
167 const char* pszMessage;
168 void* pProgressArg=NULL;
169 GDALProgressFunc pfnProgress=GDALTermProgress;
170 double progress=0;
171
172 Vector2d<vector<double> > inputData;//row,col,point
173
174
175 ImgReaderGdal maskReader;
176 ImgWriterGdal outputWriter;
177 GDALDataType theType=GDT_Unknown;
178 if(verbose_opt[0])
179 cout << "possible output data types: ";
180 for(int iType = 0; iType < GDT_TypeCount; ++iType){
181 if(verbose_opt[0])
182 cout << " " << GDALGetDataTypeName((GDALDataType)iType);
183 if( GDALGetDataTypeName((GDALDataType)iType) != NULL
184 && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
185 otype_opt[0].c_str()))
186 theType=(GDALDataType) iType;
187 }
188 if(attribute_opt[0]=="spacing"){
189 if(theType!=GDT_Float32||theType!=GDT_Float64)
190 theType=GDT_Float32;
191 }
192 if(verbose_opt[0]){
193 if(theType==GDT_Unknown)
194 cout << "Unknown output pixel type: " << otype_opt[0] << endl;
195 else
196 cout << "Output pixel type: " << GDALGetDataTypeName(theType) << endl;
197 }
198
199 double maxLRX=0;
200 double maxULY=0;
201 double minULX=0;
202 double minLRY=0;
203
204 unsigned long int totalPoints=0;
205 unsigned long int nPoints=0;
206 unsigned long int ipoint=0;
207 for(int iinput=0;iinput<input_opt.size();++iinput){
208 // assert(input_opt[iinput].find(".las")!=string::npos);
209 FileReaderLas lasReader;
210 try{
211 lasReader.open(input_opt[iinput]);
212 }
213 catch(string errorString){
214 cerr << errorString << endl;
215 exit(1);
216 }
217 catch(...){
218 cerr << "Error opening input " << input_opt[iinput] << endl;
219 exit(2);
220 }
221 nPoints=lasReader.getPointCount();
222 totalPoints+=nPoints;
223
224 if(ulx_opt[0]>=lrx_opt[0]||uly_opt[0]<=lry_opt[0]){
225 double ulx,uly,lrx,lry;
226 lasReader.getExtent(ulx,uly,lrx,lry);
227 lrx+=dx_opt[0];//pixel coordinates are referenced to upper left corner (las coordinates are centres)
228 lry-=dy_opt[0];//pixel coordinates are referenced to upper left corner (las coordinates are centres)
229 if(ulx>=lrx){
230 ulx=ulx-dx_opt[0]/2.0;
231 lrx=ulx+dx_opt[0]/2.0;
232 }
233 if(uly<=lry){
234 uly=lry+dy_opt[0]/2.0;
235 lry=lry-dy_opt[0]/2.0;
236 }
237 if(maxLRX>minULX){
238 maxLRX=(lrx>maxLRX)?lrx:maxLRX;
239 maxULY=(uly>maxULY)?uly:maxULY;
240 minULX=(ulx<minULX)?ulx:minULX;
241 minLRY=(lry<minLRY)?lry:minLRY;
242 }
243 else{//initialize
244 maxLRX=lrx;
245 maxULY=uly;
246 minULX=ulx;
247 minLRY=lry;
248 }
249 }
250 else{
251 maxLRX=lrx_opt[0];
252 maxULY=uly_opt[0];
253 minULX=ulx_opt[0];
254 minLRY=lry_opt[0];
255 }
256 lasReader.close();
257 }
258 if(verbose_opt[0]){
259 std::cout << setprecision(12) << "--ulx=" << minULX << " --uly=" << maxULY << " --lrx=" << maxLRX << " --lry=" << minLRY << std::endl;
260 std::cout << "total number of points before filtering: " << totalPoints << std::endl;
261 std::cout << "filter set to " << filter_opt[0] << std::endl;
262 if(angle_min_opt.size())
263 std::cout << "minimum scan angle: " << angle_min_opt[0] << std::endl;
264 if(angle_max_opt.size())
265 std::cout << "maximum scan angle: " << angle_max_opt[0] << std::endl;
266 // std::cout << "postFilter set to " << postFilter_opt[0] << std::endl;
267 }
268 int ncol=ceil(maxLRX-minULX)/dx_opt[0];//number of columns in outputGrid
269 int nrow=ceil(maxULY-minLRY)/dy_opt[0];//number of rows in outputGrid
270 //todo: multiple bands
271 int nband=(composite_opt[0]=="profile")? nbin_opt[0] : 1;
272 if(!output_opt.size()){
273 cerr << "Error: no output file defined" << endl;
274 exit(1);
275 }
276 if(verbose_opt[0])
277 cout << "opening output file " << output_opt[0] << endl;
278 outputWriter.open(output_opt[0],ncol,nrow,nband,theType,oformat_opt[0],option_opt);
279 outputWriter.GDALSetNoDataValue(nodata_opt[0]);
280 //set projection
281 double gt[6];
282 gt[0]=minULX;
283 gt[1]=dx_opt[0];
284 gt[2]=0;
285 gt[3]=maxULY;
286 gt[4]=0;
287 gt[5]=-dy_opt[0];
288 outputWriter.setGeoTransform(gt);
289 if(projection_opt.size()){
290 outputWriter.setProjectionProj4(projection_opt[0]);
291 }
292 if(!outputWriter.isGeoRef())
293 cout << "Warning: output image " << output_opt[0] << " is not georeferenced!" << endl;
294 if(colorTable_opt.size())
295 outputWriter.setColorTable(colorTable_opt[0]);
296
297 inputData.clear();
298 inputData.resize(nrow,ncol);
299 Vector2d<double> outputData(nrow,ncol);
300 for(int irow=0;irow<nrow;++irow)
301 for(int icol=0;icol<ncol;++icol)
302 outputData[irow][icol]=0;
303
304 cout << "Reading " << input_opt.size() << " point cloud files" << endl;
305 pfnProgress(progress,pszMessage,pProgressArg);
306 for(int iinput=0;iinput<input_opt.size();++iinput){
307 FileReaderLas lasReader;
308 try{
309 lasReader.open(input_opt[iinput]);
310 }
311 catch(string errorString){
312 cout << errorString << endl;
313 exit(1);
314 }
315 if(verbose_opt[0]){
316 if(lasReader.isCompressed())
317 cout << "Reading compressed point cloud " << input_opt[iinput]<< endl;
318 else
319 cout << "Reading uncompressed point cloud " << input_opt[iinput] << endl;
320 }
321 //set bounding filter
322 // lasReader.addBoundsFilter(minULX,maxULY,maxLRX,minLRY);
323 //set returns filter
324 if(returns_opt.size())
325 lasReader.addReturnsFilter(returns_opt);
326 if(classes_opt.size())
327 lasReader.addClassFilter(classes_opt);
328 lasReader.setFilters();
329
330 if(attribute_opt[0]!="z"){
331 vector<boost::uint16_t> returnsVector;
332 vector<string>::iterator ait=attribute_opt.begin();
333 while(ait!=attribute_opt.end()){
334 if(*ait=="intensity"){
335 if(verbose_opt[0])
336 std::cout << "writing intensity" << std::endl;
337 ++ait;
338 }
339 else if(*ait=="angle"){
340 if(verbose_opt[0])
341 std::cout << "writing angle" << std::endl;
342 ++ait;
343 }
344 else if(*ait=="return"){
345 if(verbose_opt[0])
346 std::cout << "writing return number" << std::endl;
347 ++ait;
348 }
349 else if(*ait=="nreturn"){
350 if(verbose_opt[0])
351 std::cout << "writing number of returns" << std::endl;
352 ++ait;
353 }
354 else if(*ait=="spacing"){
355 if(verbose_opt[0])
356 std::cout << "writing spacing" << std::endl;
357 ++ait;
358 }
359 else
360 attribute_opt.erase(ait);
361 }
362 }
363
364 // liblas::Point thePoint(&(lasReader.getHeader()));
365 // while(lasReader.readNextPoint(thePoint)){
366 OGRSpatialReference projectionRef(outputWriter.getProjectionRef().c_str());
367 OGRPoint ogrPoint;
368 OGRPoint ogrCenter;
369 if(attribute_opt[0]=="spacing"){
370 ogrPoint.assignSpatialReference(&projectionRef);
371 ogrCenter.assignSpatialReference(&projectionRef);
372 }
373 while(lasReader.getReader()->ReadNextPoint()){
374 liblas::Point const& thePoint = lasReader.getReader()->GetPoint();
375 // liblas::Point const& thePoint=lasReader.getPoint();
376 progress=static_cast<float>(ipoint)/totalPoints;
377 pfnProgress(progress,pszMessage,pProgressArg);
378 double theX=thePoint.GetX();
379 double theY=thePoint.GetY();
380 if(verbose_opt[0]>1)
381 cout << "reading point " << ipoint << endl;
382 if(theX<minULX||theX>=maxLRX||theY>=maxULY||theY<minLRY)
383 continue;
384 if((filter_opt[0]=="single")&&(thePoint.GetNumberOfReturns()!=1))
385 continue;
386 if((filter_opt[0]=="multiple")&&(thePoint.GetNumberOfReturns()<2))
387 continue;
388 if((filter_opt[0]=="last")&&(thePoint.GetReturnNumber()!=thePoint.GetNumberOfReturns()))
389 continue;
390 if((filter_opt[0]=="first")&&(thePoint.GetReturnNumber()!=1))
391 continue;
392 if(angle_min_opt.size()){
393 if(angle_min_opt[0]>thePoint.GetScanAngleRank())
394 continue;
395 }
396 if(angle_max_opt.size()){
397 if(angle_max_opt[0]<thePoint.GetScanAngleRank())
398 continue;
399 }
400 double dcol,drow;
401 outputWriter.geo2image(theX,theY,dcol,drow);
402 int icol=static_cast<int>(dcol);
403 int irow=static_cast<int>(drow);
404 if(irow<0||irow>=nrow){
405 // //test
406 // cout << "Error: thePoint.GetX(),thePoint.GetY(),dcol,drow" << thePoint.GetX() << ", " << thePoint.GetY() << ", " << dcol << ", " << drow << endl;
407 continue;
408 }
409 if(icol<0||icol>=ncol){
410 // //test
411 // cout << "Error: thePoint.GetX(),thePoint.GetY(),dcol,drow" << thePoint.GetX() << ", " << thePoint.GetY() << ", " << dcol << ", " << drow << endl;
412 continue;
413 }
414 assert(irow>=0);
415 assert(irow<nrow);
416 assert(icol>=0);
417 assert(icol<ncol);
418 if(composite_opt[0]=="number")
419 outputData[irow][icol]+=1;
420 else if(attribute_opt[0]=="z")
421 inputData[irow][icol].push_back(thePoint.GetZ());
422 else if(attribute_opt[0]=="intensity")
423 inputData[irow][icol].push_back(thePoint.GetIntensity());
424 else if(attribute_opt[0]=="angle")
425 inputData[irow][icol].push_back(thePoint.GetScanAngleRank());
426 else if(attribute_opt[0]=="return")
427 inputData[irow][icol].push_back(thePoint.GetReturnNumber());
428 else if(attribute_opt[0]=="nreturn")
429 inputData[irow][icol].push_back(thePoint.GetNumberOfReturns());
430 else if(attribute_opt[0]=="spacing"){
431 ogrPoint.setX(theX);
432 ogrPoint.setY(theY);
433 double centerX;
434 double centerY;
435 outputWriter.image2geo(icol,irow,centerX,centerY);
436 ogrCenter.setX(centerX);
437 ogrCenter.setY(centerY);
438 inputData[irow][icol].push_back(ogrPoint.Distance(&ogrCenter));
439 }
440 else{
441 std::string errorString="attribute not supported";
442 throw(errorString);
443 }
444 ++ipoint;
445 }
446 if(verbose_opt[0])
447 std::cout << "number of points: " << ipoint << std::endl;
448 lasReader.close();
449 }
450 progress=1;
451 pfnProgress(progress,pszMessage,pProgressArg);
452
453 std::cout << "processing LiDAR points" << std::endl;
454 progress=0;
455 pfnProgress(progress,pszMessage,pProgressArg);
457 //fill inputData in outputData
458 // if(composite_opt[0]=="profile"){
459 // assert(postFilter_opt[0]=="none");
460 // for(int iband=0;iband<nband;++iband)
461 // outputProfile[iband].resize(nrow,ncol);
462 // }
463 for(int irow=0;irow<nrow;++irow){
464 if(composite_opt[0]=="number")
465 continue;//outputData already set
466 Vector2d<double> outputProfile(nband,ncol);
467 for(int icol=0;icol<ncol;++icol){
468 std::vector<double> profile;
469 if(!inputData[irow][icol].size())
470 outputData[irow][icol]=(static_cast<double>((nodata_opt[0])));
471 else{
473 if(composite_opt[0]=="min")
474 outputData[irow][icol]=stat.mymin(inputData[irow][icol]);
475 else if(composite_opt[0]=="max")
476 outputData[irow][icol]=stat.mymax(inputData[irow][icol]);
477 else if(composite_opt[0]=="absmin")
478 outputData[irow][icol]=stat.absmin(inputData[irow][icol]);
479 else if(composite_opt[0]=="absmax")
480 outputData[irow][icol]=stat.absmax(inputData[irow][icol]);
481 else if(composite_opt[0]=="median")
482 outputData[irow][icol]=stat.median(inputData[irow][icol]);
483 else if(composite_opt[0]=="percentile")
484 outputData[irow][icol]=stat.percentile(inputData[irow][icol],inputData[irow][icol].begin(),inputData[irow][icol].end(),percentile_opt[0]);
485 else if(composite_opt[0]=="mean")
486 outputData[irow][icol]=stat.mean(inputData[irow][icol]);
487 else if(composite_opt[0]=="sum")
488 outputData[irow][icol]=stat.sum(inputData[irow][icol]);
489 else if(composite_opt[0]=="first")
490 outputData[irow][icol]=inputData[irow][icol][0];
491 else if(composite_opt[0]=="last")
492 outputData[irow][icol]=inputData[irow][icol].back();
493 else if(composite_opt[0]=="profile"){
494 if(inputData[irow][icol].size()<2){
495 for(int iband=0;iband<nband;++iband)
496 outputProfile[iband][icol]=static_cast<double>(nodata_opt[0]);
497 continue;
498 }
499 double min=0;
500 double max=0;
501 stat.minmax(inputData[irow][icol],inputData[irow][icol].begin(),inputData[irow][icol].end(),min,max);
502 if(verbose_opt[0])
503 std::cout << "min,max: " << min << "," << max << std::endl;
504 if(max>min){
505 stat.percentiles(inputData[irow][icol],inputData[irow][icol].begin(),inputData[irow][icol].end(),profile,nband,min,max);
506 assert(profile.size()==nband);
507 for(int iband=0;iband<nband;++iband)
508 outputProfile[iband][icol]=profile[iband];
509 }
510 else{
511 for(int iband=0;iband<nband;++iband)
512 outputProfile[iband][icol]=max;
513 }
514 }
515 else{
516 std::cout << "Error: composite_opt " << composite_opt[0] << " not supported" << std::endl;
517 exit(2);
518 }
519 }
520 }
521 if(composite_opt[0]=="profile"){
522 for(int iband=0;iband<nband;++iband){
523 // assert(outputProfile[iband].size()==outputWriter.nrOfRow());
524 assert(outputProfile[iband].size()==outputWriter.nrOfCol());
525 try{
526 outputWriter.writeData(outputProfile[iband],irow,iband);
527 }
528 catch(std::string errorString){
529 cout << errorString << endl;
530 exit(1);
531 }
532 }
533 }
534 progress=static_cast<float>(irow)/outputWriter.nrOfRow();
535 pfnProgress(progress,pszMessage,pProgressArg);
536 }
537 progress=1;
538 pfnProgress(progress,pszMessage,pProgressArg);
539 inputData.clear();//clean up memory
540 //apply post filter
541 // std::cout << "Applying post processing filter: " << postFilter_opt[0] << std::endl;
542 // if(postFilter_opt[0]=="etew_min"){
543 // if(composite_opt[0]!="min")
544 // std::cout << "Warning: composite option is not set to min!" << std::endl;
545 // //Elevation Threshold with Expand Window (ETEW) Filter (p.73 frmo Airborne LIDAR Data Processing and Analysis Tools ALDPAT 1.0)
546 // //first iteration is performed assuming only minima are selected using options -fir all -comp min
547 // unsigned long int nchange=1;
548 // //increase cells and thresholds until no points from the previous iteration are discarded.
549 // int dimx=dimx_opt[0];
550 // int dimy=dimy_opt[0];
551 // filter2d::Filter2d morphFilter;
552 // // morphFilter.setNoValue(0);
553 // Vector2d<float> currentOutput=outputData;
554 // int iteration=1;
555 // while(nchange&&iteration<=maxIter_opt[0]){
556 // double hThreshold=maxSlope_opt[0]*dimx;
557 // Vector2d<float> newOutput;
558 // nchange=morphFilter.morphology(currentOutput,newOutput,"erode",dimx,dimy,disc_opt[0],hThreshold);
559 // currentOutput=newOutput;
560 // dimx+=2;//change from theory: originally double cellCize
561 // dimy+=2;//change from theory: originally double cellCize
562 // std::cout << "iteration " << iteration << ": " << nchange << " pixels changed" << std::endl;
563 // ++iteration;
564 // }
565 // outputData=currentOutput;
566 // }
567 // else if(postFilter_opt[0]=="promorph"||postFilter_opt[0]=="bunting"){
568 // if(composite_opt[0]!="min")
569 // std::cout << "Warning: composite option is not set to min!" << std::endl;
570 // assert(hThreshold_opt.size()>1);
571 // //Progressive morphological filter tgrs2003_zhang vol41 pp 872-882
572 // //first iteration is performed assuming only minima are selected using options -fir all -comp min
573 // //increase cells and thresholds until no points from the previous iteration are discarded.
574 // int dimx=dimx_opt[0];
575 // int dimy=dimy_opt[0];
576 // filter2d::Filter2d theFilter;
577 // // theFilter.setNoValue(0);
578 // Vector2d<float> currentOutput=outputData;
579 // double hThreshold=hThreshold_opt[0];
580 // int iteration=1;
581 // while(iteration<=maxIter_opt[0]){
582 // std::cout << "iteration " << iteration << " with window size " << dimx << " and dh_max: " << hThreshold << std::endl;
583 // Vector2d<float> newOutput;
584 // try{
585 // theFilter.morphology(outputData,currentOutput,"erode",dimx,dimy,disc_opt[0],hThreshold);
586 // theFilter.morphology(currentOutput,outputData,"dilate",dimx,dimy,disc_opt[0],hThreshold);
587 // if(postFilter_opt[0]=="bunting"){
588 // theFilter.doit(outputData,currentOutput,"median",dimx,dimy,1,disc_opt[0]);
589 // outputData=currentOutput;
590 // }
591 // }
592 // catch(std::string errorString){
593 // cout << errorString << endl;
594 // exit(1);
595 // }
596 // int newdimx=(dimx==1)? 3: 2*(dimx-1)+1;
597 // int newdimy=(dimx==1)? 3: 2*(dimy-1)+1;//from PE&RS vol 71 pp313-324
598 // hThreshold=hThreshold_opt[0]+maxSlope_opt[0]*(newdimx-dimx)*dx_opt[0];
599 // dimx=newdimx;
600 // dimy=newdimy;
601 // if(hThreshold>hThreshold_opt[1])
602 // hThreshold=hThreshold_opt[1];
603 // ++iteration;
604 // }
605 // outputData=currentOutput;
606 // }
607 // else if(postFilter_opt[0]=="open"){
608 // if(composite_opt[0]!="min")
609 // std::cout << "Warning: composite option is not set to min!" << std::endl;
610 // filter2d::Filter2d morphFilter;
611 // // morphFilter.setNoValue(0);
612 // Vector2d<float> filterInput=outputData;
613 // try{
614 // morphFilter.morphology(outputData,filterInput,"erode",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
615 // morphFilter.morphology(filterInput,outputData,"dilate",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
616 // }
617 // catch(std::string errorString){
618 // cout << errorString << endl;
619 // exit(1);
620 // }
621 // }
622 // else if(postFilter_opt[0]=="close"){
623 // if(composite_opt[0]!="max")
624 // std::cout << "Warning: composite option is not set to max!" << std::endl;
625 // filter2d::Filter2d morphFilter;
626 // // morphFilter.setNoValue(0);
627 // Vector2d<float> filterInput=outputData;
628 // try{
629 // morphFilter.morphology(outputData,filterInput,"dilate",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
630 // morphFilter.morphology(filterInput,outputData,"erode",dimx_opt[0],dimy_opt[0],disc_opt[0],maxSlope_opt[0]);
631 // }
632 // catch(std::string errorString){
633 // cout << errorString << endl;
634 // exit(1);
635 // }
636 // }
637 if(composite_opt[0]!="profile"){
638 //write output file
639 std::cout << "writing output raster file" << std::endl;
640 progress=0;
641 pfnProgress(progress,pszMessage,pProgressArg);
642 for(int irow=0;irow<nrow;++irow){
643 try{
644 assert(outputData.size()==outputWriter.nrOfRow());
645 assert(outputData[0].size()==outputWriter.nrOfCol());
646 outputWriter.writeData(outputData[irow],irow,0);
647 }
648 catch(std::string errorString){
649 cout << errorString << endl;
650 exit(1);
651 }
652 progress=static_cast<float>(irow)/outputWriter.nrOfRow();
653 pfnProgress(progress,pszMessage,pProgressArg);
654 }
655 }
656 progress=1;
657 pfnProgress(progress,pszMessage,pProgressArg);
658 if(verbose_opt[0])
659 std::cout << "closing lasReader" << std::endl;
660 outputWriter.close();
661}
int nrOfRow(void) const
Get the number of rows of this dataset.
bool isGeoRef() const
Is this dataset georeferenced (pixel size in y must be negative) ?
bool geo2image(double x, double y, double &i, double &j) const
Convert georeferenced coordinates (x and y) to image coordinates (column and row)
CPLErr setGeoTransform(double *gt)
Set the geotransform data for this dataset.
std::string getProjectionRef(void) const
Get the projection reference.
CPLErr setProjectionProj4(const std::string &projection)
Set the projection for this dataset from user input (supports epsg:<number> format)
int nrOfCol(void) const
Get the number of columns of this dataset.
Definition: ImgRasterGdal.h:98
bool image2geo(double i, double j, double &x, double &y) const
Convert image coordinates (column and row) to georeferenced coordinates (x and y)
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 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...
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)
bool writeData(T &value, int col, int row, int band=0)
Write a single pixel cell value at a specific column and row for a specific band (all indices start c...
Definition: ImgWriterGdal.h:96