Skip to content

Commit

Permalink
update pdal_wrench
Browse files Browse the repository at this point in the history
  • Loading branch information
@alexbruy @wonder-sk
alexbruy authored and wonder-sk committed Apr 13, 2023
1 parent c4774d6 commit ee471bf
Show file tree
Hide file tree
Showing 16 changed files with 826 additions and 170 deletions.
17 changes: 16 additions & 1 deletion external/pdal_wrench/alg.cpp
View file
Expand Up @@ -18,6 +18,7 @@
#include <thread>

#include <pdal/QuickInfo.hpp>
#include <pdal/util/Bounds.hpp>

using namespace pdal;

Expand Down Expand Up @@ -79,14 +80,15 @@ bool runAlg(std::vector<std::string> args, Alg &alg)


bool Alg::parseArgs(std::vector<std::string> args)
{
{
pdal::Arg* argInput = nullptr;
if (hasSingleInput)
{
argInput = &programArgs.add("input,i", "Input point cloud file", inputFile);
}

(void)programArgs.add("filter,f", "Filter expression for input data", filterExpression);
(void)programArgs.add("bounds", "Filter by rectangle", filterBounds);

addArgs(); // impl in derived

Expand All @@ -112,6 +114,19 @@ bool Alg::parseArgs(std::vector<std::string> args)
return false;
}

if (!filterBounds.empty())
{
try
{
parseBounds(filterBounds);
}
catch (pdal::Bounds::error err)
{
std::cerr << "invalid bounds: " << err.what() << std::endl;
return false;
}
}

if (!checkArgs()) // impl in derived class
return false;

Expand Down
2 changes: 2 additions & 0 deletions external/pdal_wrench/alg.hpp
View file
Expand Up @@ -42,6 +42,8 @@ struct Alg

std::string filterExpression; // optional argument to limit input points

std::string filterBounds; // optional clipping rectangle for input (pdal::Bounds)

bool needsSingleCrs = true; // most algs assume that all input files in VPC are in the same CRS,
// and only few exceptions (e.g. info) tolerate mixture of multiple CRS

Expand Down
37 changes: 29 additions & 8 deletions external/pdal_wrench/boundary.cpp
View file
Expand Up @@ -68,6 +68,32 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, double r

Stage& r = manager->makeReader(tile->inputFilenames[0], "");

Stage *last = &r;

// filtering
if (!tile->filterBounds.empty())
{
Options filter_opts;
filter_opts.add(pdal::Option("bounds", tile->filterBounds));

if (readerSupportsBounds(r))
{
// Reader of the format can do the filtering - use that whenever possible!
r.addOptions(filter_opts);
}
else
{
// Reader can't do the filtering - do it with a filter
last = &manager->makeFilter( "filters.crop", *last, filter_opts);
}
}
if (!tile->filterExpression.empty())
{
Options filter_opts;
filter_opts.add(pdal::Option("expression", tile->filterExpression));
last = &manager->makeFilter( "filters.expression", *last, filter_opts);
}

// TODO: what edge size? (by default samples 5000 points if not specified
// TODO: set threshold ? (default at least 16 points to keep the cell)
// btw. if threshold=0, there are still missing points because of simplification (smooth=True)
Expand All @@ -78,13 +104,8 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, double r
hexbin_opts.add(pdal::Option("edge_size", resolution));
}
hexbin_opts.add(pdal::Option("threshold", pointsThreshold));
(void)manager->makeFilter( "filters.hexbin", *last, hexbin_opts );

if (!tile->filterExpression.empty())
{
hexbin_opts.add(pdal::Option("where", tile->filterExpression));
}

(void)manager->makeFilter( "filters.hexbin", r, hexbin_opts );
return manager;
}

Expand All @@ -99,14 +120,14 @@ void Boundary::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& p

for (const VirtualPointCloud::File& f : vpc.files)
{
ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression);
ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression, filterBounds);
tile.inputFilenames.push_back(f.filename);
pipelines.push_back(pipeline(&tile, resolution, pointsThreshold));
}
}
else
{
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression);
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression, filterBounds);
tile.inputFilenames.push_back(inputFile);
pipelines.push_back(pipeline(&tile, resolution, pointsThreshold));
}
Expand Down
37 changes: 27 additions & 10 deletions external/pdal_wrench/clip.cpp
View file
Expand Up @@ -118,21 +118,38 @@ static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, const pd

Stage& r = manager->makeReader( tile->inputFilenames[0], "");

Stage& f = manager->makeFilter( "filters.crop", r, crop_opts );
Stage *last = &r;

pdal::Options writer_opts;
writer_opts.add(pdal::Option("forward", "all"));

Stage& w = manager->makeWriter( tile->outputFilename, "", f, writer_opts);
// filtering
if (!tile->filterBounds.empty())
{
Options filter_opts;
filter_opts.add(pdal::Option("bounds", tile->filterBounds));

if (readerSupportsBounds(r))
{
// Reader of the format can do the filtering - use that whenever possible!
r.addOptions(filter_opts);
}
else
{
// Reader can't do the filtering - do it with a filter
last = &manager->makeFilter( "filters.crop", *last, filter_opts);
}
}
if (!tile->filterExpression.empty())
{
Options filter_opts;
filter_opts.add(pdal::Option("where", tile->filterExpression));
f.addOptions(filter_opts);
w.addOptions(filter_opts);
filter_opts.add(pdal::Option("expression", tile->filterExpression));
last = &manager->makeFilter( "filters.expression", *last, filter_opts);
}

last = &manager->makeFilter( "filters.crop", *last, crop_opts );

pdal::Options writer_opts;
writer_opts.add(pdal::Option("forward", "all"));
manager->makeWriter( tile->outputFilename, "", *last, writer_opts);

return manager;
}

Expand Down Expand Up @@ -175,7 +192,7 @@ void Clip::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pipel
std::cout << "using " << f.filename << std::endl;
}

ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression);
ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression, filterBounds);
tile.inputFilenames.push_back(f.filename);

// for input file /x/y/z.las that goes to /tmp/hello.vpc,
Expand All @@ -190,7 +207,7 @@ void Clip::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pipel
}
else
{
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression);
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression, filterBounds);
tile.inputFilenames.push_back(inputFile);
tile.outputFilename = outputFile;
pipelines.push_back(pipeline(&tile, crop_opts));
Expand Down
85 changes: 65 additions & 20 deletions external/pdal_wrench/density.cpp
View file
Expand Up @@ -73,20 +73,53 @@ std::unique_ptr<PipelineManager> Density::pipeline(ParallelJobInfo *tile) const
readers.push_back(&manager->makeReader(f, ""));
}

if (tile->mode == ParallelJobInfo::Spatial)
std::vector<Stage*> last = readers;

// find out what will be the bounding box for this job
// (there could be also no bbox if there's no "bounds" filter and no tiling)
BOX2D filterBox = !tile->filterBounds.empty() ? parseBounds(tile->filterBounds).to2d() : BOX2D();
BOX2D box = intersectTileBoxWithFilterBox(tile->box, filterBox);

if (box.valid())
{
// We are going to do filtering of points based on 2D box. Ideally we want to do
// the filtering in the reader (if the reader can do it efficiently like copc/ept),
// otherwise we have to add filters.crop stage to filter points after they were read

for (Stage* reader : readers)
{
// with COPC files, we can also specify bounds at the reader
// that will only read the required parts of the file
if (reader->getName() == "readers.copc")
if (readerSupportsBounds(*reader))
{
// add "bounds" option to reader
pdal::Options copc_opts;
copc_opts.add(pdal::Option("threads", 1));
copc_opts.add(pdal::Option("bounds", box_to_pdal_bounds(tile->box)));
copc_opts.add(pdal::Option("bounds", box_to_pdal_bounds(box)));
reader->addOptions(copc_opts);
}
}

if (!allReadersSupportBounds(readers) && !tile->filterBounds.empty())
{
// At least some readers can't do the filtering - do it with a filter
Options filter_opts;
filter_opts.add(pdal::Option("bounds", tile->filterBounds));
Stage *filterCrop = &manager->makeFilter( "filters.crop", filter_opts);
for (Stage *s : last)
filterCrop->setInput(*s);
last.clear();
last.push_back(filterCrop);
}
}

if (!tile->filterExpression.empty())
{
Options filter_opts;
filter_opts.add(pdal::Option("expression", tile->filterExpression));
Stage *filterExpr = &manager->makeFilter( "filters.expression", filter_opts);
for (Stage *s : last)
filterExpr->setInput(*s);
last.clear();
last.push_back(filterExpr);
}

pdal::Options writer_opts;
Expand All @@ -98,9 +131,9 @@ std::unique_ptr<PipelineManager> Density::pipeline(ParallelJobInfo *tile) const
writer_opts.add(pdal::Option("gdalopts", "TILED=YES"));
writer_opts.add(pdal::Option("gdalopts", "COMPRESS=DEFLATE"));

if (tile->box.valid())
if (box.valid())
{
BOX2D box2 = tile->box;
BOX2D box2 = box;
// fix tile size - PDAL's writers.gdal adds one pixel (see GDALWriter::createGrid()),
// because it probably expects that that the bounds and resolution do not perfectly match
box2.maxx -= resolution;
Expand All @@ -109,20 +142,13 @@ std::unique_ptr<PipelineManager> Density::pipeline(ParallelJobInfo *tile) const
writer_opts.add(pdal::Option("bounds", box_to_pdal_bounds(box2)));
}

if (!tile->filterExpression.empty())
{
writer_opts.add(pdal::Option("where", tile->filterExpression));
}

// TODO: "writers.gdal: Requested driver 'COG' does not support file creation.""
// writer_opts.add(pdal::Option("gdaldriver", "COG"));

pdal::StageCreationOptions opts{ tile->outputFilename, "", nullptr, writer_opts, "" };
Stage& w = manager->makeWriter( opts );
for (Stage *stage : readers)
{
w.setInput(*stage); // connect all readers to the writer
}
for (Stage *stage : last)
w.setInput(*stage);

return manager;
}
Expand Down Expand Up @@ -172,7 +198,14 @@ void Density::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pi
// to avoid empty areas in resulting rasters
tileBox.clip(gridBounds);

ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression);
if (!filterBounds.empty() && !intersectionBox2D(tileBox, parseBounds(filterBounds).to2d()).valid())
{
if (verbose)
std::cout << "skipping tile " << iy << " " << ix << " -- " << tileBox.toBox() << std::endl;
continue;
}

ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression, filterBounds);
for (const VirtualPointCloud::File & f: vpc.overlappingBox2D(tileBox))
{
tile.inputFilenames.push_back(f.filename);
Expand Down Expand Up @@ -226,7 +259,14 @@ void Density::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pi
{
BOX2D tileBox = t.boxAt(ix, iy);

ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression);
if (!filterBounds.empty() && !intersectionBox2D(tileBox, parseBounds(filterBounds).to2d()).valid())
{
if (verbose)
std::cout << "skipping tile " << iy << " " << ix << " -- " << tileBox.toBox() << std::endl;
continue;
}

ParallelJobInfo tile(ParallelJobInfo::Spatial, tileBox, filterExpression, filterBounds);
tile.inputFilenames.push_back(inputFile);

// create temp output file names
Expand All @@ -245,7 +285,7 @@ void Density::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pi
{
// single input LAS/LAZ - no parallelism

ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression);
ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression, filterBounds);
tile.inputFilenames.push_back(inputFile);
tile.outputFilename = outputFile;
pipelines.push_back(pipeline(&tile));
Expand All @@ -256,8 +296,13 @@ void Density::preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pi

void Density::finalize(std::vector<std::unique_ptr<PipelineManager>>& pipelines)
{
if (pipelines.size() > 1)
if (!tileOutputFiles.empty())
{
rasterTilesToCog(tileOutputFiles, outputFile);

// clean up the temporary directory
fs::path outputParentDir = fs::path(outputFile).parent_path();
fs::path outputSubdir = outputParentDir / fs::path(outputFile).stem();
fs::remove_all(outputSubdir);
}
}

0 comments on commit ee471bf

Please sign in to comment.