Add support for multiple tiled vpc overviews

src/vpc.cpp

@@ -10,6 +10,8 @@
  *                                                                           *
  ****************************************************************************/
 
+#include <algorithm>
+#include <cmath>
 #include <iostream>
 #include <fstream>
 #include <filesystem>
@@ -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<std::string> 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<std::string> 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,100 +553,215 @@ void buildVpc(std::vector<std::string> args)
 
     //
 
-    std::string overviewFilenameBase, overviewFilenameCopc;
-    std::vector<std::string> 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<int>(std::ceil(totalW / overviewLength)));
+        numOverviewY = (std::max)(1, static_cast<int>(std::ceil(totalH / overviewLength)));
+    }
+
+    const int numOverviewTiles = numOverviewX * numOverviewY;
+    const bool multiOverview = numOverviewTiles > 1;
+
+    struct OverviewTile {
+        BOX2D bbox;
+        std::string copcFilename;
+        std::vector<std::string> tempFiles;
+    };
+    std::vector<OverviewTile> 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)
     {
         std::map<std::string, Stage*> hexbinFilters, statsFilters;
         std::vector<std::unique_ptr<PipelineManager>> pipelines;
 
+        int overviewCounter = 0;
         for (VirtualPointCloud::File &f : vpc.files)
         {
-            std::unique_ptr<PipelineManager> 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<PipelineManager> 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<PipelineManager> 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<PipelineManager> 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<PipelineManager> manager( new PipelineManager );
             std::vector<std::unique_ptr<PipelineManager>> 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<PipelineManager> 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<std::string> 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);
+                }
             }
         }
     }

src/vpc.hpp

@@ -64,9 +64,9 @@ struct VirtualPointCloud
         std::vector<SchemaItem> schema;  // we're not using it, just for STAC export
         std::vector<StatsItem> 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<std::string> overviewFilenames;
     };
 
     std::vector<File> files;