exporting patch: # HG changeset patch # User Peter Jansen # Date 1273382390 -36000 # Node ID 74ddb0782fefebd4d72ad7c0d1527c0e365017de # Parent b0e958797e3e5b8aea4ffd5a3df68a4415f5bf29 add lasthin diff -r b0e958797e3e -r 74ddb0782fef apps/Makefile.am --- a/apps/Makefile.am Fri May 07 10:53:37 2010 -0500 +++ b/apps/Makefile.am Sun May 09 15:19:50 2010 +1000 @@ -22,6 +22,7 @@ AM_CPPFLAGS = -I../include/liblas/capi -I../include $(GDAL_CPPFLAGS) $(GEOTIFF_CPPFLAGS) $(SPATIALINDEX_CPPFLAGS) $(OCI_CPPFLAGS) $(ZLIB_CPPFLAGS) $(BOOST_FLAGS) +lasthin = lasthin.c lascommon.c lasinfo_SOURCES = lasinfo.c lascommon.c las2las_SOURCES = las2las.c lascommon.c lasmerge_SOURCES = lasmerge.c lascommon.c @@ -32,7 +33,7 @@ lasindex_SOURCES = lasindex.cpp lasindex_LDADD = ../src/liblas.la -bin_PROGRAMS = lasinfo las2las lasmerge las2txt txt2las ts2las +bin_PROGRAMS = lasthin lasinfo las2las lasmerge las2txt txt2las ts2las if SPATIALINDEX_IS_CONFIG bin_PROGRAMS += lasindex diff -r b0e958797e3e -r 74ddb0782fef apps/lasthin.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/apps/lasthin.cpp Sun May 09 15:19:50 2010 +1000 @@ -0,0 +1,360 @@ +/* +=============================================================================== + + FILE: lasthin.cpp + + CONTENTS: + + This tool is a very simplistic point thinner. It overlays a uniform grid + and keeps the lowest point with each cell. Optionally we can request to + only consider last returns (-last_only) or a particular classification + (-keep_classification 2) before doing the thinning. + + To be done: As we are essentially gridding the data we also allow this grid to be + output in form of an ESRI ARC/INFO ASCII GRID. Before outputting the grid + one can request possible empty grid cells to be flodded with values from + neighbour cells. With iterative flooding we can fill larger empty spots + (-flood 4). + + PROGRAMMERS: + + martin isenburg@cs.unc.edu + + COPYRIGHT: + + copyright (C) 2007 martin isenburg@cs.unc.edu + + This software is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + + CHANGE HISTORY: + + 12 March 2009 -- updated to ask for input if started without arguments + 19 April 2008 -- created after not going on Sheker's full moon hike + +=============================================================================== +*/ + +/* liblas */ +//#include +#include "liblas/lasreader.hpp" +#include "liblas/laswriter.hpp" + +#include +#include +#include +#include +#include +#include + +//#include "lasreader.h" +//#include "laswriter.h" + +#ifdef _WIN32 +extern "C" FILE* fopenGzipped(const char* filename, const char* mode); +#endif + +#ifdef _WIN32 +#define compare_no_case(a,b,n) _strnicmp( (a), (b), (n) ) +#else +#define compare_no_case(a,b,n) strncasecmp( (a), (b), (n) ) +#endif + +using namespace liblas; + +std::ostream* OpenOutput(std::string filename) +{ + std::ostream* ostrm; + std::ios::openmode m; + m = std::ios::out | std::ios::binary | std::ios::ate; + + if (compare_no_case(filename.c_str(),"STOUT",5) == 0) + { + ostrm = &std::cout; + } + else + { + ostrm = new std::ofstream(filename.c_str(), m); + } + + if (!ostrm->good()) + { + delete ostrm; + throw std::runtime_error("Writing stream was not able to be created"); + exit(1); + } + + return ostrm; +} + + +void usage(bool wait=false) +{ + fprintf(stderr,"usage:\n"); + fprintf(stderr,"lasthin in.las out.las\n"); + fprintf(stderr,"lasthin -i in.las -grid_spacing 0.5 -o out.las\n"); + fprintf(stderr,"lasthin -i in.las -last_only -grid_spacing 1.0 -remove_extra_header -o out.las\n"); + fprintf(stderr,"lasthin -i in.las -keep_class 2 -keep_class 3 -keep_class 4 -grid_spacing 0.5 -olas > out.las\n"); + fprintf(stderr,"lasthin -h\n"); + if (wait) + { + fprintf(stderr,"\n"); + getc(stdin); + } + exit(1); +} + +static void byebye(bool wait=false) +{ + if (wait) + { + fprintf(stderr,"\n"); + getc(stdin); + } + exit(1); +} + +void ptime(const char *const msg) +{ + float t= ((float)clock())/CLOCKS_PER_SEC; + fprintf(stderr, "cumulative CPU time thru %s = %f\n", msg, t); +} + +int main(int argc, char *argv[]) +{ + int i; + bool verbose = false; + char* file_name_in = 0; + char* file_name_out = 0; + bool olas = false; + bool olaz = false; + bool remove_extra_header = false; + int keep_classification[32] = {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1}; + bool last_only = false; + double grid_spacing = 1.0; + + for (i = 1; i < argc; i++) + { + if (strcmp(argv[i],"-verbose") == 0) + { + verbose = true; + } + else if (strcmp(argv[i],"-h") == 0) + { + usage(); + } + else if (strcmp(argv[i],"-i") == 0) + { + i++; + file_name_in = argv[i]; + } + else if (strcmp(argv[i],"-o") == 0) + { + i++; + file_name_out = argv[i]; + } + else if (strcmp(argv[i],"-olas") == 0) + { + olas = true; + } + else if (strcmp(argv[i],"-olaz") == 0) + { + olaz = true; + } + else if (strcmp(argv[i],"-grid_spacing") == 0) + { + i++; + grid_spacing = atof(argv[i]); + } + else if (strcmp(argv[i],"-keep_classification") == 0 || strcmp(argv[i],"-keep_class") == 0) + { + i++; + int j = 0; + while (keep_classification[j] != -1) j++; + keep_classification[j] = atoi(argv[i]); + } + else if (strcmp(argv[i],"-last_only") == 0) + { + last_only = true; + } + else if (strcmp(argv[i],"-remove_extra_header") == 0) + { + remove_extra_header = true; + } + else if (i == argc - 2 && file_name_in == 0 && file_name_out == 0) + { + file_name_in = argv[i]; + } + else if (i == argc - 1 && file_name_in == 0 && file_name_out == 0) + { + file_name_in = argv[i]; + } + else if (i == argc - 1 && file_name_in != 0 && file_name_out == 0) + { + file_name_out = argv[i]; + } + else + { + fprintf(stderr,"cannot understand argument '%s'\n", argv[i]); + usage(); + } + } + + // open input file + + std::ifstream ifs; + if (!liblas::Open(ifs, file_name_in)) + { + throw std::runtime_error(std::string("Can not open \'") + file_name_in + "\'"); + } + liblas::Reader lasreader(ifs); + + Header inHeader = lasreader.GetHeader(); + Header outHeader = lasreader.GetHeader(); + + // create the grid + int lowest_x = (int)((inHeader.GetMinX() / grid_spacing) + .5); + int lowest_y = (int)((inHeader.GetMinY() / grid_spacing) + .5); + int highest_x = (int)((inHeader.GetMaxX() / grid_spacing) + .5); + int highest_y = (int)((inHeader.GetMaxY() / grid_spacing) + .5); + int size_x = highest_x-lowest_x+1; + int size_y = highest_y-lowest_y+1; + int size = size_x*size_y; + + fprintf(stderr, "thinning points onto %d by %d = %d grid (grid_spacing = %.2f unit)\n", size_x, size_y, size, grid_spacing); + + int* elevation_grid = new int[size]; + if (elevation_grid == 0) + { + fprintf(stderr, "ERROR: could not allocate grid of size %d by %d\n", size_x, size_y); + byebye(argc==1); + } + for (i = 0; i < size; i++) elevation_grid[i] = 2147483647; + + // do the first pass + + unsigned int surviving_number_of_point_records = 0; + unsigned int surviving_number_of_points_by_return[] = {0,0,0,0,0,0,0,0}; + int eliminated_last_only = 0; + int eliminated_classification = 0; + int eliminated_thinning = 0; + + if (verbose) ptime("start."); + fprintf(stderr, "first pass reading %ld points ...\n", inHeader.GetPointRecordsCount()); + + // loop over points + + Point p; + while (lasreader.ReadNextPoint()) + { + p = lasreader.GetPoint(); + // put the points through two tests + if (last_only && p.GetReturnNumber() != p.GetNumberOfReturns()) + { + eliminated_last_only++; + continue; + } + if (keep_classification[0] != -1) + { + int j = 0; + while (keep_classification[j] != -1 && p.GetClassification().GetClass() != keep_classification[j]) j++; + if (keep_classification[j] == -1) + { + eliminated_classification++; + continue; + } + } + + int point_x = (int)((p.GetX() / grid_spacing) + .5); + int point_y = (int)((p.GetY() / grid_spacing) + .5); + int pos_x = point_x-lowest_x; + int pos_y = point_y-lowest_y; + int pos = pos_y*size_x+pos_x; + + if (elevation_grid[pos] == 2147483647) + { + surviving_number_of_point_records++; + surviving_number_of_points_by_return[p.GetReturnNumber()-1]++; + elevation_grid[pos] = p.GetZ(); + } + else + { + eliminated_thinning++; + if (p.GetZ() < elevation_grid[pos]) + { + elevation_grid[pos] = p.GetZ(); + } + } + } + + if (eliminated_last_only) fprintf(stderr, "eliminated based on last returns only: %d\n", eliminated_last_only); + if (eliminated_classification) fprintf(stderr, "eliminated based on classification: %d\n", eliminated_classification); + if (eliminated_thinning) fprintf(stderr, "eliminated based on thinning: %d\n", eliminated_thinning); + + fprintf(stderr, "grid saturation is %d of %d point (%.2f percent)\n", surviving_number_of_point_records, size, 100.0f*surviving_number_of_point_records/size); + + // quit if no output specified + + if (file_name_out == 0 && !olas && !olaz) + { + fprintf(stderr, "no output specified. exiting ...\n"); + exit(1); + } + + // re-open input file + + lasreader.Reset(); + + // prepare the header for the surviving points + + outHeader.SetPointRecordsCount(surviving_number_of_point_records); + for (i = 0; i < 5; i++) outHeader.SetPointRecordsByReturnCount(i, surviving_number_of_points_by_return[i]); + + fprintf(stderr, "second pass reading %ld and writing %ld points ...\n", inHeader.GetPointRecordsCount(), outHeader.GetPointRecordsCount()); + + // open output file + + std::ofstream ofs; + ofs.open(file_name_out); + liblas::Writer laswriter(ofs, outHeader); + + // loop over points + while (lasreader.ReadNextPoint()) + { + p = lasreader.GetPoint(); + + // put the points through all the tests + if (last_only && p.GetReturnNumber() != p.GetNumberOfReturns()) + { + continue; + } + if (keep_classification[0] != -1) + { + int j = 0; + while (keep_classification[j] != -1 && p.GetClassification().GetClass() != keep_classification[j]) j++; + if (keep_classification[j] == -1) continue; + } + + int point_x = (int)((p.GetX() / grid_spacing) + .5); + int point_y = (int)((p.GetY() / grid_spacing) + .5); + int pos_x = point_x-lowest_x; + int pos_y = point_y-lowest_y; + int pos = pos_y*size_x+pos_x; + + if (elevation_grid[pos] == (int)p.GetZ()) + { + elevation_grid[pos]--; + laswriter.WritePoint(p); + } + } + + ifs.close(); + ofs.close(); + + if (verbose) ptime("done."); + + byebye(argc==1); + + return 0; +}