diff --git a/src/vpc.cpp b/src/vpc.cpp index dd78477..f48a291 100644 --- a/src/vpc.cpp +++ b/src/vpc.cpp @@ -10,6 +10,8 @@ * * ****************************************************************************/ +#include +#include #include #include #include @@ -152,17 +154,20 @@ bool VirtualPointCloud::read(std::string filename) statsItem["variance"])); } - // read overview file (if any, expecting at most one) - // this logic is very basic, we should be probably checking roles of assets - if (f["assets"].contains("overview")) + // read overview files (assets with "overview" role) + for (auto &asset : f["assets"]) { - vpcf.overviewFilename = f["assets"]["overview"]["href"]; + if (!asset.contains("roles") || !asset["roles"].is_array()) + continue; - if (vpcf.overviewFilename.substr(0, 2) == "./") - { - // resolve relative path - vpcf.overviewFilename = fs::weakly_canonical(filenameParent / vpcf.overviewFilename).string(); - } + const auto roles = asset["roles"]; + if (std::find(roles.cbegin(), roles.cend(), "overview") == roles.cend()) + continue; + + std::string ovFilename = asset["href"]; + if (ovFilename.substr(0, 2) == "./") + ovFilename = fs::weakly_canonical(filenameParent / ovFilename).string(); + vpcf.overviewFilenames.push_back(ovFilename); } files.push_back(vpcf); @@ -311,26 +316,20 @@ bool VirtualPointCloud::write(std::string filename) }; nlohmann::json assets = { { "data", dataAsset } }; - if (!f.overviewFilename.empty()) + for (size_t i = 0; i < f.overviewFilenames.size(); ++i) { - std::string overviewFilename; - if (pdal::Utils::isRemote(f.overviewFilename)) - { - // keep remote URLs as they are - overviewFilename = f.overviewFilename; - } - else + std::string ovFilename(f.overviewFilenames[i]); + if (!pdal::Utils::isRemote(ovFilename)) { - // turn local paths to relative - fs::path fRelative = fs::relative(f.overviewFilename, outputPath); - overviewFilename = "./" + fRelative.string(); + const fs::path fRelative = fs::relative(ovFilename, outputPath); + ovFilename = "./" + fRelative.string(); } - + const std::string key = f.overviewFilenames.size() > 1 ? ("overview-" + std::to_string(i + 1)) : "overview"; nlohmann::json overviewAsset = { - { "href", overviewFilename }, + { "href", ovFilename }, { "roles", json::array({"overview"}) }, }; - assets["overview"] = overviewAsset; + assets[key] = overviewAsset; } jFiles.push_back( @@ -418,6 +417,7 @@ void buildVpc(std::vector args) bool boundaries = false; bool stats = false; bool overview = false; + double overviewLength = 0.0; int max_threads = -1; bool verbose = false; bool help = false; @@ -430,6 +430,9 @@ void buildVpc(std::vector args) programArgs.add("boundary", "Calculate boundary polygons from data", boundaries); programArgs.add("stats", "Calculate statistics from data", stats); programArgs.add("overview", "Create overview point cloud from source data", overview); + programArgs.add("overview-length", + "Split overview into multiple tiles of specified maximum edge length in CRS units (implies --overview)", + overviewLength); pdal::Arg& argThreads = programArgs.add("threads", "Max number of concurrent threads for parallel runs", max_threads); programArgs.add("verbose", "Print extra debugging output", verbose); @@ -550,16 +553,74 @@ void buildVpc(std::vector args) // - std::string overviewFilenameBase, overviewFilenameCopc; - std::vector overviewTempFiles; - int overviewCounter = 0; + if (overviewLength > 0.0) + overview = true; + + std::string overviewFilenameBase; if (overview) { - // for /tmp/hello.vpc we will use /tmp/hello-overview.laz as overview file - fs::path outputParentDir = fs::path(outputFile).parent_path(); - fs::path outputStem = outputParentDir / fs::path(outputFile).stem(); + // for /tmp/hello.vpc we will use /tmp/hello-overview.copc.laz as overview file + // if multiple overview tiles are generated, we will use /tmp/hello-overview-1.copc.laz etc. + const fs::path outputParentDir = fs::path(outputFile).parent_path(); + const fs::path outputStem = outputParentDir / fs::path(outputFile).stem(); overviewFilenameBase = outputStem.string(); - overviewFilenameCopc = outputStem.string() + "-overview.copc.laz"; + } + + // Compute overview tiling grid + int numOverviewX = 1, numOverviewY = 1; + const BOX2D totalBox = vpc.box3d().to2d(); + + if (overview && overviewLength > 0.0) + { + const double totalW = totalBox.maxx - totalBox.minx; + const double totalH = totalBox.maxy - totalBox.miny; + numOverviewX = (std::max)(1, static_cast(std::ceil(totalW / overviewLength))); + numOverviewY = (std::max)(1, static_cast(std::ceil(totalH / overviewLength))); + } + + const int numOverviewTiles = numOverviewX * numOverviewY; + const bool multiOverview = numOverviewTiles > 1; + + struct OverviewTile { + BOX2D bbox; + std::string copcFilename; + std::vector tempFiles; + }; + std::vector overviewTiles; + + if (overview) + { + overviewTiles.reserve(numOverviewTiles); + int overviewCounter = 0; + for (int iy = 0; iy < numOverviewY; ++iy) + { + for (int ix = 0; ix < numOverviewX; ++ix) + { + if (multiOverview) + { + OverviewTile ovTile; + ovTile.bbox = BOX2D(totalBox.minx + ix * overviewLength, + totalBox.miny + iy * overviewLength, + totalBox.minx + (ix + 1) * overviewLength, + totalBox.miny + (iy + 1) * overviewLength); + + // overview tiles are on a grid, some might not overlap with vpc files + // so we skip them to preserve continuous filename numbering + if (vpc.overlappingBox2D(ovTile.bbox).empty()) + continue; + + ovTile.copcFilename = overviewFilenameBase + "-overview-" + std::to_string(++overviewCounter) + ".copc.laz"; + overviewTiles.push_back(ovTile); + } + else + { + OverviewTile ovTile; + ovTile.bbox = totalBox; + ovTile.copcFilename = overviewFilenameBase + "-overview.copc.laz"; + overviewTiles.push_back(ovTile); + } + } + } } if (boundaries || stats || overview) @@ -567,83 +628,140 @@ void buildVpc(std::vector args) std::map hexbinFilters, statsFilters; std::vector> pipelines; + int overviewCounter = 0; for (VirtualPointCloud::File &f : vpc.files) { - std::unique_ptr manager( new PipelineManager ); - - Stage* last = &makeReader(manager.get(), f.filename); - if (boundaries) + if (overview && multiOverview) { - pdal::Options hexbin_opts; - // TODO: any options? - last = &manager->makeFilter( "filters.hexbin", *last, hexbin_opts ); - hexbinFilters[f.filename] = last; - } + // For multi-tile overview: one pipeline per (source file, tile) combination + for (OverviewTile &tile : overviewTiles) + { + if (!tile.bbox.overlaps(f.bbox.to2d())) + continue; - if (stats) - { - pdal::Options stats_opts; - // TODO: any options? - last = &manager->makeFilter( "filters.stats", *last, stats_opts ); - statsFilters[f.filename] = last; - } + std::unique_ptr manager( new PipelineManager ); + Stage* last = &makeReader(manager.get(), f.filename); - if (overview) + pdal::Options crop_opts; + crop_opts.add(pdal::Option("bounds", box_to_pdal_bounds(tile.bbox))); + last = &manager->makeFilter("filters.crop", *last, crop_opts); + + pdal::Options decim_opts; + decim_opts.add(pdal::Option("step", 1000)); + last = &manager->makeFilter("filters.decimation", *last, decim_opts); + + const std::string overviewOutput = overviewFilenameBase + "-overview-tmp-" + std::to_string(++overviewCounter) + ".las"; + tile.tempFiles.push_back(overviewOutput); + makeWriter(manager.get(), overviewOutput, last); + pipelines.push_back(std::move(manager)); + } + + // Boundaries/stats for multi-overview still need a separate pipeline per file + if (boundaries || stats) + { + std::unique_ptr manager( new PipelineManager ); + Stage* last = &makeReader(manager.get(), f.filename); + if (boundaries) + { + pdal::Options hexbin_opts; + last = &manager->makeFilter("filters.hexbin", *last, hexbin_opts); + hexbinFilters[f.filename] = last; + } + if (stats) + { + pdal::Options stats_opts; + last = &manager->makeFilter("filters.stats", *last, stats_opts); + statsFilters[f.filename] = last; + } + pipelines.push_back(std::move(manager)); + } + } + else { - // TODO: configurable method and step size? - pdal::Options decim_opts; - decim_opts.add(pdal::Option("step", 1000)); - last = &manager->makeFilter( "filters.decimation", *last, decim_opts ); + std::unique_ptr manager( new PipelineManager ); + + Stage* last = &makeReader(manager.get(), f.filename); + if (boundaries) + { + pdal::Options hexbin_opts; + // TODO: any options? + last = &manager->makeFilter( "filters.hexbin", *last, hexbin_opts ); + hexbinFilters[f.filename] = last; + } + + if (stats) + { + pdal::Options stats_opts; + // TODO: any options? + last = &manager->makeFilter( "filters.stats", *last, stats_opts ); + statsFilters[f.filename] = last; + } + + if (overview) + { + // TODO: configurable method and step size? + pdal::Options decim_opts; + decim_opts.add(pdal::Option("step", 1000)); + last = &manager->makeFilter( "filters.decimation", *last, decim_opts ); - std::string overviewOutput = overviewFilenameBase + "-overview-tmp-" + std::to_string(++overviewCounter) + ".las"; - overviewTempFiles.push_back(overviewOutput); + const std::string overviewOutput = overviewFilenameBase + "-overview-tmp-" + std::to_string(++overviewCounter) + ".las"; + overviewTiles[0].tempFiles.push_back(overviewOutput); - makeWriter(manager.get(), overviewOutput, last); - } + makeWriter(manager.get(), overviewOutput, last); + } - pipelines.push_back(std::move(manager)); + pipelines.push_back(std::move(manager)); + } } runPipelineParallel(vpc.totalPoints(), true, pipelines, max_threads, verbose); if (overview) { - // When doing overviews, this is the second stage where we index overview point cloud. + // When doing overviews, this is the second stage where we index overview point cloud(s). // We do it separately because writers.copc is not streamable. We could also use // untwine instead of writers.copc... - std::unique_ptr manager( new PipelineManager ); std::vector> pipelinesCopcOverview; - // TODO: I am not really sure why we need a merge filter, but without it - // I am only getting points in output COPC from the last reader. Example - // from the documentation suggests the merge filter should not be needed: - // https://pdal.io/en/latest/stages/writers.copc.html + for (const OverviewTile &tile : overviewTiles) + { + if (tile.tempFiles.empty()) + continue; - Stage &merge = manager->makeFilter("filters.merge"); + std::unique_ptr manager( new PipelineManager ); - pdal::Options writer_opts; - //writer_opts.add(pdal::Option("forward", "all")); - Stage& writer = manager->makeWriter(overviewFilenameCopc, "writers.copc", merge, writer_opts); - (void)writer; + // TODO: I am not really sure why we need a merge filter, but without it + // I am only getting points in output COPC from the last reader. Example + // from the documentation suggests the merge filter should not be needed: + // https://pdal.io/en/latest/stages/writers.copc.html - for (const std::string &overviewTempFile : overviewTempFiles) - { - Stage& reader = makeReader(manager.get(), overviewTempFile); - merge.setInput(reader); + Stage &merge = manager->makeFilter("filters.merge"); + + pdal::Options writer_opts; + Stage& writer = manager->makeWriter(tile.copcFilename, "writers.copc", merge, writer_opts); + (void)writer; + + for (const std::string &tempFile : tile.tempFiles) + { + Stage& reader = makeReader(manager.get(), tempFile); + merge.setInput(reader); + } + + pipelinesCopcOverview.push_back(std::move(manager)); } if (verbose) { - std::cout << "Indexing overview point cloud..." << std::endl; + std::cout << "Indexing overview point cloud(s)..." << std::endl; } - pipelinesCopcOverview.push_back(std::move(manager)); runPipelineParallel(vpc.totalPoints()/1000, false, pipelinesCopcOverview, max_threads, verbose); // delete tmp overviews - for (const std::string &overviewTempFile : overviewTempFiles) + for (const OverviewTile &tile : overviewTiles) { - std::filesystem::remove(overviewTempFile); + for (const std::string &tempFile : tile.tempFiles) + std::filesystem::remove(tempFile); } } @@ -676,7 +794,17 @@ void buildVpc(std::vector args) } if (overview) { - f.overviewFilename = overviewFilenameCopc; + for (const OverviewTile &tile : overviewTiles) + { + if (tile.tempFiles.empty()) + continue; // empty tile (no source files) + // we can't use overlaps() here as we don't want to include bboxes that only touch + if (tile.bbox.minx < f.bbox.maxx && + tile.bbox.maxx > f.bbox.minx && + tile.bbox.miny < f.bbox.maxy && + tile.bbox.maxy > f.bbox.miny) + f.overviewFilenames.push_back(tile.copcFilename); + } } } } diff --git a/src/vpc.hpp b/src/vpc.hpp index 3c35633..21bf454 100644 --- a/src/vpc.hpp +++ b/src/vpc.hpp @@ -64,9 +64,9 @@ struct VirtualPointCloud std::vector schema; // we're not using it, just for STAC export std::vector stats; - // support for overview point clouds - currently we assume a file refers to at most a single overview file - // (when building VPC with overviews, we create one overview file for all source data) - std::string overviewFilename; + // support for overview point clouds - a file may refer to one or more overview files + // (when building VPC with overviews, multiple overview tiles may be created) + std::vector overviewFilenames; }; std::vector files;