render.cpp

#include <filesystem>
#include <cmath>
#include <algorithm>
#include <iomanip>
#include <sstream>

#include <pdal/StageFactory.hpp>
#include <pdal/PointTable.hpp>
#include <pdal/PointView.hpp>
#include <pdal/io/BufferReader.hpp>

#include "render.hpp"
#include "utils.hpp"

namespace fs = std::filesystem;

struct Tile{
    double radius;
    Extent bounds;
    Extent bufferedBounds;
    std::string filename;
};

void render(PointSet *pset, const std::string &outDir, const std::string &outputType,
        int tileSize, 
        const std::vector<double> &radiuses, double resolution, 
        int maxTiles, bool force){

    if (outputType != "max" && outputType != "idw")
        throw std::runtime_error("Unsupported output-type: " + outputType);

    fs::path pOutDir = fs::path(outDir);
    std::vector<double> rads(radiuses);

    if (fs::exists(pOutDir)){
        if (!force) throw std::runtime_error(outDir + " exists (use --force to overwrite results)");
    }else{
        fs::create_directories(pOutDir);
    }

    // Generate tile list
    unsigned int width = static_cast<int>(std::ceil(pset->extent.width() / resolution));
    unsigned int height = static_cast<int>(std::ceil(pset->extent.height() / resolution));
    
    // Set a floor, no matter the resolution parameter
    // (sometimes a wrongly estimated scale of the model can cause the resolution
    // to be set unrealistically low, causing errors)
    const unsigned int RES_FLOOR = 64;

    if (width < RES_FLOOR && height < RES_FLOOR){
        double prev_width = width;
        double prev_height = height;
        
        if (width >= height){
            width = RES_FLOOR;
            height = static_cast<unsigned int>(std::ceil(pset->extent.height() / pset->extent.width() * RES_FLOOR));
        } else {
            width = static_cast<unsigned int>(std::ceil(pset->extent.width() / pset->extent.height() * RES_FLOOR));
            height = RES_FLOOR;
        }

        double floor_ratio = prev_width / width;
        resolution *= floor_ratio;

        for (size_t i = 0; i < rads.size(); i++){
            rads[i] *= floor_ratio;
        }

        std::cout << "Really low resolution DEM requested (" << prev_width << ", " << prev_height << ") will set floor at " << RES_FLOOR << " pixels. Resolution changed to " << resolution << ". The scale of this reconstruction might be off." << std::endl;
    }

    unsigned int numSplitsX = static_cast<int>(std::max<double>(1.0, std::ceil(width / static_cast<double>(tileSize))));
    unsigned int numSplitsY = static_cast<int>(std::max<double>(1.0, std::ceil(height / static_cast<double>(tileSize))));
    
    unsigned int numTiles = numSplitsX * numSplitsY;

    std::cout << "DEM resolution is (" << width << ", " << height << "), max tile size is " << tileSize << ", will split DEM generation into " << numTiles << " tiles" << std::endl;

    if (maxTiles > 0){
        if (numTiles > maxTiles){
            std::cerr << "Max tiles limit exceeded (" << maxTiles << "). This is a strong indicator that the reconstruction failed" << std::endl;
            exit(1);
        }
    }

    double tileBoundsWidth = pset->extent.width() / static_cast<double>(numSplitsX);
    double tileBoundsHeight = pset->extent.height() / static_cast<double>(numSplitsY);

    std::vector<Tile> tiles;

    double minx;
    double maxx;
    double miny;
    double maxy;

    for (const double &r: rads){
        minx = pset->extent.minx;
        for (unsigned int x = 0; x < numSplitsX; x++){
            miny = pset->extent.miny;
            maxx = x == numSplitsX - 1 ?
                            pset->extent.maxx : 
                            minx + tileBoundsWidth;

            for (unsigned int y = 0; y < numSplitsY; y++){
                maxy = y == numSplitsY - 1 ? 
                                pset->extent.maxy : 
                                miny + tileBoundsHeight;

                std::stringstream ss;
                ss << "r" << r << "_x" << x << "_y" << y << ".tif"; 

                Tile t;
                t.filename = (fs::absolute(pOutDir) / ss.str()).string();
                t.bounds.minx = minx;
                t.bounds.maxx = maxx;
                t.bounds.miny = miny;
                t.bounds.maxy = maxy;
                t.radius = r;

                const double buffer = r * 2;
                t.bufferedBounds.minx = t.bounds.minx - buffer;
                t.bufferedBounds.maxx = t.bounds.maxx + buffer;
                t.bufferedBounds.miny = t.bounds.miny - buffer;
                t.bufferedBounds.maxy = t.bounds.maxy + buffer;

                tiles.push_back(t);

                miny = maxy;
            }
            
            minx = maxx;
        }
    }

    // Sort tiles by decreasing radius
    std::sort(tiles.begin(), tiles.end(), 
        [](Tile const &a, Tile const &b) {
            return a.radius < b.radius; 
        });
    
    #pragma omp parallel for
    for (int i = 0; i < static_cast<int>(tiles.size()); i++){
        const Tile &t = tiles[i];

        // Build a PointView containing only points within the buffered tile bounds.
        // The buffer (2x radius) ensures the interpolation has context beyond the tile edge.
        pdal::PointTable table;
        table.layout()->registerDim(pdal::Dimension::Id::X);
        table.layout()->registerDim(pdal::Dimension::Id::Y);
        table.layout()->registerDim(pdal::Dimension::Id::Z);

        pdal::PointViewPtr view(new pdal::PointView(table));

        pdal::PointId idx = 0;
        for (size_t j = 0; j < pset->size(); j++){
            if (pset->x[j] >= t.bufferedBounds.minx && pset->x[j] <= t.bufferedBounds.maxx &&
                pset->y[j] >= t.bufferedBounds.miny && pset->y[j] <= t.bufferedBounds.maxy){
                view->setField(pdal::Dimension::Id::X, idx, pset->x[j]);
                view->setField(pdal::Dimension::Id::Y, idx, pset->y[j]);
                view->setField(pdal::Dimension::Id::Z, idx, pset->z[j]);
                idx++;
            }
        }

        bool empty = (idx == 0);

        if (!empty){
            pdal::BufferReader reader;
            reader.addView(view);

            std::ostringstream boundsStr;
            boundsStr << std::fixed << std::setprecision(12)
                      << "([" << t.bounds.minx << "," << t.bounds.maxx << "],"
                      <<  "[" << t.bounds.miny << "," << t.bounds.maxy << "])";

            pdal::Options wOpts;
            wOpts.add("filename",    t.filename);
            wOpts.add("resolution",  resolution);
            wOpts.add("radius",      t.radius);
            wOpts.add("output_type", outputType);
            wOpts.add("bounds",      boundsStr.str());
            wOpts.add("data_type",   "float");
            wOpts.add("nodata",      -9999.0);
            if (pset->srs.valid())
                wOpts.add("override_srs", pset->srs.getWKT());

            pdal::StageFactory factory;
            pdal::Stage *writer = factory.createStage("writers.gdal");
            writer->setOptions(wOpts);
            writer->setInput(reader);
            writer->prepare(table);
            writer->execute(table);
        }

        #pragma omp critical
        {
            std::cout << fs::path(t.filename).filename().string() << (empty ? " [Empty]" : "") << std::endl;
        }
    }

}