Raster2PointCloud.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#include <tclap/CmdLine.h>
#include <vtkPolyDataAlgorithm.h>
#include <vtkSmartPointer.h>
#include <vtkXMLPolyDataWriter.h>
 
#include <algorithm>
#include <memory>
#include <string>
 
#include "Applications/DataExplorer/VtkVis/VtkGeoImageSource.h"
#include "Applications/DataExplorer/VtkVis/VtkImageDataToPointCloudFilter.h"
#include "BaseLib/FileTools.h"
#include "BaseLib/Logging.h"
#include "BaseLib/MPI.h"
#include "BaseLib/TCLAPArguments.h"
#include "GeoLib/IO/AsciiRasterInterface.h"
#include "InfoLib/GitInfo.h"
 
int main(int argc, char* argv[])
{
    TCLAP::CmdLine cmd(
        "Batch conversion of raster files into VTK point cloud data.\n"
        "A series of raster files is converted into point clouds by creating a "
        "number of random points based on the intensity value in each pixel "
        "(the higher the value, the more points. The [x,y]-extend of these "
        "random points is limited by the extent of the pixel they represent, "
        "the elevation is limited by the max-elevation parameter. Likewise, "
        "the "
        "maximum number of points per pixel is limited by the max-points "
        "parameter. The increase of the number of points in relation to the "
        "pixel value is indicated by gamma, where 1 is a linear increase and "
        "values smaller or larger than 1 indicate a logarithmic or exponential "
        "increase, respectively.\n\n"
        "OpenGeoSys-6 software, version " +
            GitInfoLib::GitInfo::ogs_version +
            ".\n"
            "Copyright (c) 2012-2026, OpenGeoSys Community "
            "(http://www.opengeosys.org)",
        ' ', GitInfoLib::GitInfo::ogs_version);
    TCLAP::ValueArg<double> gamma_arg(
        "g", "gamma",
        "Exponental increase coefficient. A value of '1' indicates a linear "
        "increase, smaller or larger values indicate a logarithmic or "
        "exponential increase, respectively (default: 1).",
        false, 0.0, "exponential increase");
    cmd.add(gamma_arg);
    TCLAP::ValueArg<std::size_t> max_points_arg(
        "p", "max-points",
        "Maximum number of points per pixel (default: 1000), "
        "(min = 0).",
        false, 0, "maximum points per pixel");
    cmd.add(max_points_arg);
    TCLAP::ValueArg<double> max_height_arg(
        "e", "max-elevation",
        "Maximum elevation of randomly created points "
        "(default: 5000).",
        false, 0.0, "maximum point elevation");
    cmd.add(max_height_arg);
    TCLAP::ValueArg<std::string> output_arg(
        "o", "output",
        "Output (.vtp). Path to output polydata files. "
        "Specifying 'file.vtp' will "
        "make the programme write a series of files called 'file0.vtp', "
        "'file1.vtp', etc.",
        true, "", "OUTPUT_FILE");
    cmd.add(output_arg);
    TCLAP::ValueArg<std::string> input_arg(
        "i", "input",
        "Input (.asc). Path to input raster files. "
        "Specifying 'file.asc' will make "
        "the programme look for a series of files called 'file0.asc', "
        "'file1.asc', etc.",
        true, "", "INPUT_FILE");
    cmd.add(input_arg);
    auto log_level_arg = BaseLib::makeLogLevelArg();
    cmd.add(log_level_arg);
    cmd.parse(argc, argv);
 
    BaseLib::MPI::Setup mpi_setup(argc, argv);
    BaseLib::initOGSLogger(log_level_arg.getValue());
 
    std::string const input_name = input_arg.getValue().c_str();
    std::string const output_name = output_arg.getValue().c_str();
 
    double const max_height =
        (max_height_arg.isSet()) ? max_height_arg.getValue() : 5000;
    std::size_t const max_points =
        (max_points_arg.isSet()) ? max_points_arg.getValue() : 1000;
    double const gamma = (gamma_arg.isSet()) ? gamma_arg.getValue() : 1;
 
    std::string const ibase_name = BaseLib::dropFileExtension(input_name);
    std::string const ifile_ext = BaseLib::getFileExtension(input_name);
 
    std::string const obase_name = BaseLib::dropFileExtension(output_name);
    std::string const ofile_ext = BaseLib::getFileExtension(output_name);
 
    // Rainevents using point cloud filter
    std::cout << "Starting export... ";
    double global_min = std::numeric_limits<double>::max();
    double global_max = 0;
 
    std::size_t n_rasters = 0;
    for (std::size_t i = 0;; ++i)
    {
        std::string const file_name =
            ibase_name + std::to_string(i) + ifile_ext;
        std::unique_ptr<GeoLib::Raster> r(
            FileIO::AsciiRasterInterface::getRasterFromASCFile(file_name));
        if (r == nullptr)
        {
            n_rasters = i;
            break;
        }
 
        for (auto it = r->begin(); it != r->end(); ++it)
        {
            if (*it != r->getHeader().no_data)
            {
                global_min = std::min(*it, global_min);
                global_max = std::max(*it, global_max);
            }
        }
    }
 
    for (std::size_t i = 0; i < n_rasters; ++i)
    {
        std::string const file_name =
            ibase_name + std::to_string(i) + ifile_ext;
        vtkSmartPointer<VtkGeoImageSource> geo_image =
            vtkSmartPointer<VtkGeoImageSource>::New();
        geo_image->readImage(QString::fromStdString(file_name));
        VtkImageDataToPointCloudFilter* const point_cloud_filter =
            VtkImageDataToPointCloudFilter::New();
        point_cloud_filter->SetInputConnection(geo_image->GetOutputPort());
        geo_image->Update();
        point_cloud_filter->SetMinValueRange(global_min);
        point_cloud_filter->SetMaxValueRange(global_max);
        point_cloud_filter->SetMinHeight(0);
        point_cloud_filter->SetMaxHeight(max_height);
        point_cloud_filter->SetIsLinear(false);
        point_cloud_filter->SetMinNumberOfPointsPerCell(0);
        point_cloud_filter->SetMaxNumberOfPointsPerCell(max_points);
        point_cloud_filter->SetGamma(gamma);
        point_cloud_filter->Update();
        vtkSmartPointer<vtkPolyDataAlgorithm> const pd_alg =
            dynamic_cast<vtkPolyDataAlgorithm*>(point_cloud_filter);
        if (pd_alg)
        {
            std::string const output =
                obase_name + std::to_string(i) + ofile_ext;
            vtkSmartPointer<vtkXMLPolyDataWriter> pd_writer =
                vtkSmartPointer<vtkXMLPolyDataWriter>::New();
            pd_writer->SetInputData(pd_alg->GetOutputDataObject(0));
            pd_writer->SetFileName(output.c_str());
            pd_writer->Write();
        }
    }
    std::cout << "done." << std::endl;
    return EXIT_SUCCESS;
}