LiDAR.h

Go to the documentation of this file.

 
#ifndef LIDARPLUGIN
#define LIDARPLUGIN
 
#include "CollisionDetection.h"
#include "Context.h"
#include "Visualizer.h"
 
#include "triangulation_cdt.h"
 
#include <functional>
 
template<class datatype>
class HitTable {
public:
    uint Ntheta, Nphi;
 
    HitTable(void) {
        Ntheta = 0;
        Nphi = 0;
    }
    HitTable(const int nx, const int ny) {
        Ntheta = nx;
        Nphi = ny;
        data.resize(Nphi);
        for (int j = 0; j < Nphi; j++) {
            data.at(j).resize(Ntheta);
        }
    }
    HitTable(const int nx, const int ny, const datatype initval) {
        Ntheta = nx;
        Nphi = ny;
        data.resize(Nphi);
        for (int j = 0; j < Nphi; j++) {
            data.at(j).resize(Ntheta, initval);
        }
    }
 
    datatype get(const int i, const int j) const {
        if (i >= 0 && i < Ntheta && j >= 0 && j < Nphi) {
            return data.at(j).at(i);
        } else {
            helios::helios_runtime_error("ERROR (hit_map.get): get index out of range. Attempting to get index map at (" + std::to_string(i) + "," + std::to_string(j) + "), but size of scan is " + std::to_string(Ntheta) + " x " +
                                         std::to_string(Nphi) + ".");
        }
        return data.at(j).at(i); // unreachable; silences control-reaches-end-of-non-void-function warning
    }
    void set(const int i, const int j, const datatype value) {
        if (i >= 0 && i < Ntheta && j >= 0 && j < Nphi) {
            data.at(j).at(i) = value;
        } else {
            helios::helios_runtime_error("ERROR (hit_map.set): set index out of range. Attempting to set index map at (" + std::to_string(i) + "," + std::to_string(j) + "), but size of scan is " + std::to_string(Ntheta) + " x " +
                                         std::to_string(Nphi) + ".");
        }
    }
    void resize(const int nx, const int ny, const datatype initval) {
        Ntheta = nx;
        Nphi = ny;
        data.resize(Nphi);
        for (int j = 0; j < Nphi; j++) {
            data.at(j).resize(Ntheta);
            for (int i = 0; i < Ntheta; i++) {
                data.at(j).at(i) = initval;
            }
        }
    }
 
private:
    std::vector<std::vector<datatype>> data;
};
 
struct HitPoint {
    helios::vec3 position;
    helios::SphericalCoord direction;
    helios::int2 row_column;
    helios::RGBcolor color;
    int gridcell;
    int scanID;
    // NOTE: per-hit scalar data (intensity, distance, timestamp, custom labels, ...) is NOT stored
    // here. It lives in cloud-level columnar storage on LiDARcloud (hit_data_columns/hit_data_present),
    // indexed by the hit's position in the `hits` vector. Storing it as N independent
    // std::map<std::string,double> trees (one per hit) made bulk field extraction O(K*N) cache-cold
    // tree descents and dominated export time. See LiDARcloud::getHitData / getHitDataColumn.
    HitPoint(void) {
        position = helios::make_vec3(0, 0, 0);
        direction = helios::make_SphericalCoord(0, 0);
        row_column = helios::make_int2(0, 0);
        color = helios::RGB::red;
        gridcell = -2;
        scanID = -1;
    }
    HitPoint(int __scanID, helios::vec3 __position, helios::SphericalCoord __direction, helios::int2 __row_column, helios::RGBcolor __color) {
        scanID = __scanID;
        position = __position;
        direction = __direction;
        row_column = __row_column;
        color = __color;
        gridcell = -2;
    }
};
 
struct Triangulation {
    helios::vec3 vertex0, vertex1, vertex2;
    int ID0, ID1, ID2;
    int scanID;
    int gridcell;
    helios::RGBcolor color;
    float area;
    Triangulation(void) {
        vertex0 = helios::make_vec3(0, 0, 0);
        vertex1 = helios::make_vec3(0, 0, 0);
        vertex2 = helios::make_vec3(0, 0, 0);
        ID0 = 0;
        ID1 = 0;
        ID2 = 0;
        scanID = -1;
        gridcell = -2;
        color = helios::RGB::green;
        area = 0;
    }
    Triangulation(int __scanID, helios::vec3 __vertex0, helios::vec3 __vertex1, helios::vec3 __vertex2, int __ID0, int __ID1, int __ID2, helios::RGBcolor __color, int __gridcell) {
        scanID = __scanID;
        vertex0 = __vertex0;
        vertex1 = __vertex1;
        vertex2 = __vertex2;
        ID0 = __ID0;
        ID1 = __ID1;
        ID2 = __ID2;
        gridcell = __gridcell;
        color = __color;
 
        // calculate area
        helios::vec3 s0 = vertex1 - vertex0;
        helios::vec3 s1 = vertex2 - vertex0;
        helios::vec3 s2 = vertex2 - vertex1;
 
        float a = s0.magnitude();
        float b = s1.magnitude();
        float c = s2.magnitude();
        float s = 0.5f * (a + b + c);
 
        area = sqrt(s * (s - a) * (s - b) * (s - c));
    }
};
 
struct GridCell {
    helios::vec3 center;
    helios::vec3 global_anchor;
    helios::vec3 size;
    helios::vec3 global_size;
    helios::int3 global_ijk;
    helios::int3 global_count;
    float azimuthal_rotation;
    float leaf_area;
    float Gtheta;
    float ground_height;
    float vegetation_height;
    float maximum_height;
    // ---- LAD inversion uncertainty (sufficient statistics + sampling variance) ----
    // These quantify the statistical SAMPLING uncertainty of the leaf-area inversion,
    // conditional on the beams that entered the voxel (Pimont et al. 2018, RSE 215:343-370).
    // They do NOT capture occlusion/coverage bias (voxels shadowed so beams never penetrate).
    int beam_count = -1; 
    float I_rdi = 0.f; 
    float zbar_e = 0.f; 
    float var_path = 0.f; 
    float L1_element = -1.f; 
    float LAD_variance = -1.f; 
    bool ci_valid = false; 
    GridCell(helios::vec3 __center, helios::vec3 __global_anchor, helios::vec3 __size, helios::vec3 __global_size, float __azimuthal_rotation, helios::int3 __global_ijk, helios::int3 __global_count) {
        center = __center;
        global_anchor = __global_anchor;
        size = __size;
        global_size = __global_size;
        azimuthal_rotation = __azimuthal_rotation;
        global_ijk = __global_ijk;
        global_count = __global_count;
        leaf_area = 0;
        Gtheta = 0;
        ground_height = 0;
        vegetation_height = 0;
        maximum_height = 0;
    }
};
 
 
enum ScanPattern {
    SCAN_PATTERN_RASTER = 0,
    SCAN_PATTERN_SPINNING_MULTIBEAM = 1,
    SCAN_PATTERN_RISLEY_PRISM = 2
};
 
 
enum ScanMode {
    SCAN_MODE_STATIC_RASTER = 0,
    SCAN_MODE_MOVING_RASTER = 1,
    SCAN_MODE_SPINNING = 2,
    SCAN_MODE_RISLEY_PRISM = 3
};
 
 
enum ReturnMode {
    RETURN_MODE_MULTI = 0,
    RETURN_MODE_SINGLE = 1
};
 
 
enum SingleReturnSelection {
    SINGLE_RETURN_STRONGEST = 0,
    SINGLE_RETURN_FIRST = 1,
    SINGLE_RETURN_LAST = 2,
    SINGLE_RETURN_STRONGEST_PLUS_LAST = 3
};
 
 
struct RisleyPrism {
 
    RisleyPrism() : wedge_angle(0.0), refractive_index(1.0), rotor_rate(0.0), phase(0.0) {
    }
 
 
    RisleyPrism(double wedge_angle, double refractive_index, double rotor_rate, double phase = 0.0) : wedge_angle(wedge_angle), refractive_index(refractive_index), rotor_rate(rotor_rate), phase(phase) {
    }
 
    double wedge_angle;
    double refractive_index;
    double rotor_rate;
    double phase;
};
 
 
struct ScanMetadata {
 
    ScanMetadata();
 
 
    ScanMetadata(const helios::vec3 &origin, uint Ntheta, float thetaMin, float thetaMax, uint Nphi, float phiMin, float phiMax, float exitDiameter, float beamDivergence, float rangeNoiseStdDev, float angleNoiseStdDev,
                 const std::vector<std::string> &columnFormat, float scanTiltRoll = 0.f, float scanTiltPitch = 0.f, float scanAzimuthOffset = 0.f);
 
 
    ScanMetadata(const helios::vec3 &origin, const std::vector<float> &beamZenithAngles, uint Nphi, float phiMin, float phiMax, float exitDiameter, float beamDivergence, float rangeNoiseStdDev, float angleNoiseStdDev,
                 const std::vector<std::string> &columnFormat, float scanTiltRoll = 0.f, float scanTiltPitch = 0.f, float scanAzimuthOffset = 0.f);
 
    std::string data_file;
 
    uint Ntheta;
 
    float thetaMin;
 
    float thetaMax;
 
    uint Nphi;
 
    float phiMin;
 
    float phiMax;
 
    helios::vec3 origin;
 
 
    float exitDiameter;
 
 
    float beamDivergence;
 
 
    float rangeNoiseStdDev;
 
 
    float angleNoiseStdDev;
 
 
    ReturnMode returnMode = RETURN_MODE_MULTI;
 
 
    SingleReturnSelection singleReturnSelection = SINGLE_RETURN_STRONGEST;
 
 
    int maxReturns = 1;
 
 
    float pulseWidth = 0.f;
 
 
    float detectionThreshold = 0.05f;
 
 
    float scanTilt_roll;
 
 
    float scanTilt_pitch;
 
 
    float scanTilt_azimuth;
 
    std::vector<std::string> columnFormat;
 
 
    ScanPattern scanPattern;
 
 
    std::vector<float> beamZenithAngles;
 
 
    std::vector<RisleyPrism> risley_prisms;
 
 
    double risley_refractive_index_air = 1.0;
 
 
    helios::SphericalCoord rc2direction(uint row, uint column) const;
 
 
    helios::int2 direction2rc(const helios::SphericalCoord &direction) const;
 
    // ---- Moving-platform (mobile/airborne) scan support ----
    // When isMoving is true the scanner pose is driven by a dense timestamped 6-DOF trajectory and the
    // synthetic scan generates a per-pulse origin and orientation instead of a single static origin. The
    // static scanTilt_roll/pitch/azimuth fields are NOT applied in this mode: attitude comes entirely from
    // the trajectory quaternions composed with the fixed boresight misalignment. When isMoving is false the
    // scan behaves exactly as before (single static origin). Populated via \ref LiDARcloud::addScanMoving().
 
    bool isMoving = false;
 
    std::vector<double> traj_t;
 
    std::vector<helios::vec3> traj_pos;
 
 
    std::vector<helios::vec4> traj_quat;
 
    helios::vec3 lever_arm = helios::make_vec3(0, 0, 0);
 
    helios::vec3 boresight_rpy = helios::make_vec3(0, 0, 0);
 
    double pulse_period = 1.0;
 
    double t0 = 0.0;
 
    // ---- Self-describing acquisition descriptors ----
    // These make a scan introspectable ("how was it acquired / how fast did it spin / how many revolutions?")
    // without reverse-engineering the answer from scanPattern, isMoving, and phiMax. They are set by the
    // high-level scan-creation entry points (\ref LiDARcloud::addScanSpinning, \ref LiDARcloud::addScanMovingRaster,
    // etc.) and are purely descriptive; the underlying mechanism is still scanPattern + isMoving.
 
    ScanMode scanMode = SCAN_MODE_STATIC_RASTER;
 
    uint steps_per_rev = 0;
 
 
    double rotation_rate = 0.0;
 
 
    double n_revolutions = 0.0;
 
 
    void poseAt(double t, helios::vec3 &pos, helios::vec4 &quat) const;
};
 
class LiDARcloud {
private:
    size_t Nhits;
 
    std::vector<ScanMetadata> scans;
 
    std::vector<HitPoint> hits;
 
    // ---- Columnar per-hit scalar data ----
    // Per-hit scalar fields are stored column-wise (one contiguous array per label) rather than as one
    // std::map<std::string,double> per hit. This makes bulk field extraction a cache-linear pass instead
    // of N cache-cold red-black-tree lookups, matching the speed of XYZ/RGB export.
    // INVARIANT: every column in hit_data_columns and every mask in hit_data_present has length
    // hits.size() at all times. Column slot `s` holds label hit_data_labels[s]; a value is meaningful
    // only where hit_data_present[s][i] != 0 (a hit may be missing a value for a label). The values
    // for slot s of hit i are hit_data_columns[s][i]. These structures are kept in lockstep with the
    // `hits` vector everywhere it is mutated (addHitPoint, deleteHitPoint's swap-and-pop, clearHits).
    std::vector<std::string> hit_data_labels; 
    std::unordered_map<std::string, size_t> hit_data_label_index; 
    std::vector<std::vector<double>> hit_data_columns; 
    std::vector<std::vector<char>> hit_data_present; 
 
    size_t getOrCreateHitDataColumn(const std::string &label);
 
    void appendHitData(const std::map<std::string, double> &data);
 
    void clearHits();
 
    void transformHitOrigin(uint index, const std::function<helios::vec3(const helios::vec3 &)> &transform);
 
    helios::vec3 hitOriginOrFallback(uint index, const helios::vec3 &fallback) const;
 
    std::vector<GridCell> grid_cells;
 
    bool force_bruteforce_LAD = false;
 
    std::vector<Triangulation> triangles;
 
    std::vector<HitTable<int>> hit_tables;
 
    bool hitgridcellcomputed;
 
    bool triangulationcomputed;
 
    std::size_t triangulation_candidate_count;
    std::size_t triangulation_dropped_lmax;
    std::size_t triangulation_dropped_aspect;
    std::size_t triangulation_dropped_degenerate;
 
    int getContainingGridCell(const helios::vec3 &p) const;
 
    bool printmessages;
 
    std::function<void(float, const std::string &)> progress_callback;
 
    size_t synthetic_scan_memory_budget_bytes = 0; // 0 => automatic (see syntheticScan)
 
    static constexpr size_t SYNTHETIC_SCAN_DEFAULT_BUDGET_CPU = size_t(4) * 1024 * 1024 * 1024; // 4 GiB
 
    static constexpr size_t SYNTHETIC_SCAN_DEFAULT_BUDGET_GPU = size_t(8) * 1024 * 1024 * 1024; // 8 GiB
 
    CollisionDetection *collision_detection;
 
    volatile int *cancel_flag = nullptr;
 
    volatile int *synthetic_scan_progress = nullptr;
 
 
    void prepareUnifiedRayTracing(helios::Context *context);
 
    void finishUnifiedRayTracing();
 
 
    void castRaysUnified(size_t total_rays, helios::vec3 *ray_origins, helios::vec3 *direction, float *hit_t, float *hit_fnorm, int *hit_ID, size_t packet_size = 1);
 
    // -------- I/O --------- //
 
    // -------- RECONSTRUCTION --------- //
 
    // first index: leaf group, second index: triangle #
    std::vector<std::vector<Triangulation>> reconstructed_triangles;
 
    std::vector<std::vector<Triangulation>> reconstructed_trunk_triangles;
 
    std::vector<helios::vec3> reconstructed_alphamasks_center;
    std::vector<helios::vec2> reconstructed_alphamasks_size;
    std::vector<helios::SphericalCoord> reconstructed_alphamasks_rotation;
    std::vector<uint> reconstructed_alphamasks_gridcell;
    std::string reconstructed_alphamasks_maskfile;
    std::vector<uint> reconstructed_alphamasks_direct_flag;
 
    void leafReconstructionFloodfill();
 
    void backfillLeavesAlphaMask(const std::vector<float> &leaf_size, float leaf_aspect_ratio, float solidfraction, const std::vector<bool> &group_filter_flag);
 
    void calculateLeafAngleCDF(uint Nbins, std::vector<std::vector<float>> &CDF_theta, std::vector<std::vector<float>> &CDF_phi);
 
    void floodfill(size_t t, std::vector<Triangulation> &cloud_triangles, std::vector<int> &fill_flag, std::vector<std::vector<int>> &nodes, int tag, int depth, int maxdepth);
 
    // -------- HELPERS --------- //
 
 
    void computeGtheta(uint Ncells, uint Nscans, std::vector<float> &Gtheta, std::vector<float> &Gtheta_bar);
 
 
    bool anyScanMoving() const;
 
 
    void calculateLeafArea_inner(helios::Context *context, int min_voxel_hits, float element_width, float supplied_Gtheta);
 
 
    bool invertLAD(uint voxel_index, float P, float Gtheta, const std::vector<float> &dr_samples, int min_voxel_hits, const helios::vec3 &gridsize, float &leaf_area, helios::WarningAggregator &warnings);
 
    struct LADInversionResult {
        float leaf_area = 0.f; 
        float LAD_variance = -1.f; 
        int beam_count = 0; 
        float I_rdi = 0.f; 
        float zbar_e = 0.f; 
        float var_path = 0.f; 
        float L1_element = -1.f; 
        bool converged = false; 
        bool element_size_known = false; 
    };
 
 
    LADInversionResult invertLADWithVariance(uint voxel_index, float P, float Gtheta, const std::vector<float> &dr_samples, float sum_frac_sq, float element_width, int min_voxel_hits, const helios::vec3 &gridsize,
                                             helios::WarningAggregator &warnings);
 
 
    bool ciValidPimont(float L, float L1, int N, float confidence_level) const;
 
    // -------- MULTI-RETURN HELPERS --------- //
 
 
    bool isMultiReturnData() const;
 
 
    struct BeamGrouping {
        uint Nbeams = 0;
        std::vector<uint> beam_members; 
        std::vector<uint> beam_offsets; 
 
        uint beamSize(uint k) const {
            return beam_offsets[k + 1] - beam_offsets[k];
        }
    };
 
 
    BeamGrouping groupHitsByTimestamp(const std::vector<uint> &scan_indices) const;
 
 
    struct VoxelLattice {
        bool valid = false; 
        helios::vec3 origin; 
        helios::vec3 anchor; 
        helios::vec3 cell_extent; 
        float rotation = 0.f; 
        helios::int3 count; 
        std::vector<int> ijk_to_index; 
    };
 
 
    VoxelLattice detectVoxelLattice() const;
 
 
    static void accumulateBeamCell(const uint *return_indices, size_t Nreturns, const std::vector<float> &dr, const std::vector<uint> &hit_location, float &P_equal_numerator, float &P_equal_denominator, float &P_equal_sumsq,
                                   std::vector<float> &dr_array_cell);
 
 
    std::vector<uint> loadTreeQSM_impl(helios::Context *context, const std::string &filename, uint radial_subdivisions, bool use_colormap, const std::string &colormap_or_texture);
 
 
    std::vector<helios::vec3> gapfillMisses_timestamp(uint scanID, const bool gapfill_grid_only, const bool add_flags);
 
 
    std::vector<helios::vec3> gapfillMisses_rowcolumn(uint scanID, const bool add_flags);
 
public:
    LiDARcloud();
 
    ~LiDARcloud();
 
    static int selfTest(int argc = 0, char **argv = nullptr);
 
    void validateRayDirections();
 
    void disableMessages();
 
    void enableMessages();
 
 
    void setProgressCallback(std::function<void(float, const std::string &)> callback);
 
 
    void setCancelFlag(volatile int *flag);
 
 
    void setSyntheticScanProgressPointer(volatile int *ptr);
 
    void initializeCollisionDetection(helios::Context *context);
 
    void performUnifiedRayTracing(helios::Context *context, size_t N, int Npulse, helios::vec3 *ray_origins, helios::vec3 *direction, float *hit_t, float *hit_fnorm, int *hit_ID);
 
 
    [[nodiscard]] static float applyRangeIntensityCorrection(float intensity, float distance);
 
    // ------- SCANS -------- //
 
    uint getScanCount();
 
 
    uint addScan(ScanMetadata &newscan);
 
 
    uint addScanMoving(ScanMetadata scan, const std::vector<double> &traj_t, const std::vector<helios::vec3> &traj_pos, const std::vector<helios::vec4> &traj_quat, const helios::vec3 &lever_arm, const helios::vec3 &boresight_rpy,
                       float pulse_rate_hz, double t0 = 0.0);
 
 
    uint addScanMoving(ScanMetadata scan, const std::vector<double> &traj_t, const std::vector<helios::vec3> &traj_pos, const std::vector<helios::vec3> &traj_rpy, const helios::vec3 &lever_arm, const helios::vec3 &boresight_rpy,
                       float pulse_rate_hz, double t0 = 0.0);
 
 
    uint addScanSpinning(const std::vector<float> &beamElevationAngles, float azimuthStep_rad, float pulse_rate_hz, const std::vector<double> &traj_t, const std::vector<helios::vec3> &traj_pos, const std::vector<helios::vec4> &traj_quat,
                         const helios::vec3 &lever_arm, const helios::vec3 &boresight_rpy, float exitDiameter, float beamDivergence, float rangeNoiseStdDev, float angleNoiseStdDev,
                         const std::vector<std::string> &columnFormat = {"x", "y", "z"}, double t0 = 0.0);
 
 
    uint addScanSpinning(const std::vector<float> &beamElevationAngles, float azimuthStep_rad, float pulse_rate_hz, const std::vector<double> &traj_t, const std::vector<helios::vec3> &traj_pos, const std::vector<helios::vec3> &traj_rpy,
                         const helios::vec3 &lever_arm, const helios::vec3 &boresight_rpy, float exitDiameter, float beamDivergence, float rangeNoiseStdDev, float angleNoiseStdDev,
                         const std::vector<std::string> &columnFormat = {"x", "y", "z"}, double t0 = 0.0);
 
 
    uint addScanMovingRaster(uint Ntheta, float thetaMin, float thetaMax, uint Nphi, float phiMin, float phiMax, float pulse_rate_hz, const std::vector<double> &traj_t, const std::vector<helios::vec3> &traj_pos,
                             const std::vector<helios::vec4> &traj_quat, const helios::vec3 &lever_arm, const helios::vec3 &boresight_rpy, float exitDiameter, float beamDivergence, float rangeNoiseStdDev, float angleNoiseStdDev,
                             const std::vector<std::string> &columnFormat = {"x", "y", "z"}, double t0 = 0.0);
 
 
    uint addScanRisley(const std::vector<RisleyPrism> &prisms, double refractive_index_air, float pulse_rate_hz, const std::vector<double> &traj_t, const std::vector<helios::vec3> &traj_pos, const std::vector<helios::vec4> &traj_quat,
                       const helios::vec3 &lever_arm, const helios::vec3 &boresight_rpy, float exitDiameter, float beamDivergence, float rangeNoiseStdDev, float angleNoiseStdDev, const std::vector<std::string> &columnFormat = {"x", "y", "z"},
                       double t0 = 0.0);
 
 
    uint addScanRisley(const std::vector<RisleyPrism> &prisms, double refractive_index_air, float pulse_rate_hz, const std::vector<double> &traj_t, const std::vector<helios::vec3> &traj_pos, const std::vector<helios::vec3> &traj_rpy,
                       const helios::vec3 &lever_arm, const helios::vec3 &boresight_rpy, float exitDiameter, float beamDivergence, float rangeNoiseStdDev, float angleNoiseStdDev, const std::vector<std::string> &columnFormat = {"x", "y", "z"},
                       double t0 = 0.0);
 
 
    void addHitPoint(uint scanID, const helios::vec3 &xyz, const helios::SphericalCoord &direction);
 
 
    void addHitPoint(uint scanID, const helios::vec3 &xyz, const helios::SphericalCoord &direction, const helios::RGBcolor &color);
 
 
    void addHitPoint(uint scanID, const helios::vec3 &xyz, const helios::SphericalCoord &direction, const std::map<std::string, double> &data);
 
 
    void addHitPoint(uint scanID, const helios::vec3 &xyz, const helios::SphericalCoord &direction, const helios::RGBcolor &color, const std::map<std::string, double> &data);
 
 
    void addHitPoint(uint scanID, const helios::vec3 &xyz, const helios::int2 &row_column, const helios::RGBcolor &color, const std::map<std::string, double> &data);
 
 
    void deleteHitPoint(uint index);
 
    uint getHitCount() const;
 
 
    helios::vec3 getScanOrigin(uint scanID) const;
 
 
    uint getScanSizeTheta(uint scanID) const;
 
 
    uint getScanSizePhi(uint scanID) const;
 
 
    helios::vec2 getScanRangeTheta(uint scanID) const;
 
 
    helios::vec2 getScanRangePhi(uint scanID) const;
 
 
    float getScanBeamExitDiameter(uint scanID) const;
 
 
    std::vector<std::string> getScanColumnFormat(uint scanID) const;
 
 
    ScanPattern getScanPattern(uint scanID) const;
 
 
    std::vector<float> getScanBeamZenithAngles(uint scanID) const;
 
 
    ScanMode getScanMode(uint scanID) const;
 
 
    uint getScanStepsPerRev(uint scanID) const;
 
 
    double getScanRotationRate(uint scanID) const;
 
 
    double getScanRevolutions(uint scanID) const;
 
 
    std::vector<RisleyPrism> getScanRisleyPrisms(uint scanID) const;
 
 
    double getScanRisleyRefractiveIndexAir(uint scanID) const;
 
 
    float getScanBeamDivergence(uint scanID) const;
 
 
    float getScanRangeNoiseStdDev(uint scanID) const;
 
 
    float getScanAngleNoiseStdDev(uint scanID) const;
 
 
    ReturnMode getScanReturnMode(uint scanID) const;
 
 
    void setScanReturnMode(uint scanID, ReturnMode returnMode);
 
 
    SingleReturnSelection getScanSingleReturnSelection(uint scanID) const;
 
 
    void setScanSingleReturnSelection(uint scanID, SingleReturnSelection selection);
 
 
    int getScanMaxReturns(uint scanID) const;
 
 
    void setScanMaxReturns(uint scanID, int maxReturns);
 
 
    float getScanPulseWidth(uint scanID) const;
 
 
    void setScanPulseWidth(uint scanID, float pulseWidth);
 
 
    float getScanDetectionThreshold(uint scanID) const;
 
 
    void setScanDetectionThreshold(uint scanID, float detectionThreshold);
 
 
    float getScanTiltRoll(uint scanID) const;
 
 
    float getScanTiltPitch(uint scanID) const;
 
 
    float getScanAzimuthOffset(uint scanID) const;
 
 
    helios::vec3 getHitXYZ(uint index) const;
 
 
    helios::vec3 getHitOrigin(uint index) const;
 
 
    helios::SphericalCoord getHitRaydir(uint index) const;
 
 
    double getHitData(uint index, const char *label) const;
 
 
    void setHitData(uint index, const char *label, double value);
 
 
    bool doesHitDataExist(uint index, const char *label) const;
 
 
    int getHitDataColumnIndex(const char *label) const;
 
 
    void getHitDataColumn(const char *label, std::vector<double> &data, double absent_value = -9999) const;
 
 
    static constexpr float LIDAR_MISS_DISTANCE = 20000.f;
 
 
    static constexpr float LIDAR_RAYTRACE_MISS_T = 1001.f;
 
 
    bool isHitMiss(uint index) const;
 
 
    bool hasMisses() const;
 
 
    helios::RGBcolor getHitColor(uint index) const;
 
 
    int getHitScanID(uint index) const;
 
 
    int getHitIndex(uint scanID, uint row, uint column) const;
 
 
    int getHitGridCell(uint index) const;
 
 
    void setHitGridCell(uint index, int cell);
 
 
    void coordinateShift(const helios::vec3 &shift);
 
 
    void coordinateShift(uint scanID, const helios::vec3 &shift);
 
 
    void coordinateRotation(const helios::SphericalCoord &rotation);
 
 
    void coordinateRotation(uint scanID, const helios::SphericalCoord &rotation);
 
 
    void coordinateRotation(float rotation, const helios::vec3 &line_base, const helios::vec3 &line_direction);
 
    uint getTriangleCount() const;
 
    std::size_t getTriangulationCandidateCount() const;
 
    std::size_t getTriangulationDroppedByLmax() const;
 
    std::size_t getTriangulationDroppedByAspect() const;
 
    std::size_t getTriangulationDroppedByDegenerate() const;
 
 
    Triangulation getTriangle(uint index) const;
 
    // ------- FILE I/O --------- //
 
 
    void loadXML(const char *filename);
 
 
    void loadXML(const char *filename, bool load_grid_only);
 
 
    size_t loadASCIIFile(uint scanID, const std::string &ASCII_data_file);
 
 
    void exportTriangleNormals(const char *filename);
 
 
    void exportTriangleNormals(const char *filename, int gridcell);
 
 
    void exportTriangleAreas(const char *filename);
 
 
    void exportTriangleAreas(const char *filename, int gridcell);
 
 
    void exportTriangleInclinationDistribution(const char *filename, uint Nbins);
 
 
    void exportTriangleAzimuthDistribution(const char *filename, uint Nbins);
 
 
    void exportLeafAreas(const char *filename);
 
 
    void exportLeafAreaDensities(const char *filename);
 
 
    void exportGtheta(const char *filename);
 
 
    void exportLeafAreaUncertainty(const char *filename);
 
 
    void exportPointCloud(const char *filename, bool write_header = true);
 
 
    void exportPointCloud(const char *filename, uint scanID, bool write_header = true);
 
 
    void exportPointCloudPTX(const char *filename, uint scanID);
 
 
    void exportScans(const char *filename);
 
    // ------- VISUALIZER --------- //
 
 
    void addHitsToVisualizer(Visualizer *visualizer, uint pointsize) const;
 
 
    void addHitsToVisualizer(Visualizer *visualizer, uint pointsize, const helios::RGBcolor &point_color) const;
 
 
    void addHitsToVisualizer(Visualizer *visualizer, uint pointsize, const char *color_value) const;
 
 
    void addGridToVisualizer(Visualizer *visualizer) const;
 
 
    void addGridWireFrametoVisualizer(Visualizer *visualizer, float linewidth_pixels = 1.0f) const;
 
 
    void addGrid(const helios::vec3 &center, const helios::vec3 &size, const helios::int3 &ndiv, float rotation);
 
 
    void addTrianglesToVisualizer(Visualizer *visualizer) const;
 
 
    void addTrianglesToVisualizer(Visualizer *visualizer, uint gridcell) const;
 
 
    void addLeafReconstructionToVisualizer(Visualizer *visualizer) const;
 
 
    void addTrunkReconstructionToVisualizer(Visualizer *visualizer) const;
 
 
    void addTrunkReconstructionToVisualizer(Visualizer *visualizer, const helios::RGBcolor &trunk_color) const;
 
 
    std::vector<uint> addLeafReconstructionToContext(helios::Context *context) const;
 
 
    std::vector<uint> addLeafReconstructionToContext(helios::Context *context, const helios::int2 &subpatches) const;
 
 
    std::vector<uint> addReconstructedTriangleGroupsToContext(helios::Context *context) const;
 
 
    std::vector<uint> addTrunkReconstructionToContext(helios::Context *context) const;
 
 
    void getHitBoundingBox(helios::vec3 &boxmin, helios::vec3 &boxmax) const;
 
 
    void getGridBoundingBox(helios::vec3 &boxmin, helios::vec3 &boxmax) const;
 
    // --------- POINT FILTERING ----------- //
 
 
    void distanceFilter(float maxdistance);
 
 
    void xyzFilter(float xmin, float xmax, float ymin, float ymax, float zmin, float zmax);
 
 
    void xyzFilter(float xmin, float xmax, float ymin, float ymax, float zmin, float zmax, bool deleteOutside);
 
 
 
    void reflectanceFilter(float minreflectance);
 
 
    void scalarFilter(const char *scalar_field, float threshold, const char *comparator);
 
 
    void maxPulseFilter(const char *scalar);
 
 
    void minPulseFilter(const char *scalar);
 
 
    void firstHitFilter();
 
 
    void lastHitFilter();
 
    // ------- TRIANGULATION --------- //
 
 
    void triangulateHitPoints(float Lmax, float max_aspect_ratio);
 
    // ERK
 
    void triangulateHitPoints(float Lmax, float max_aspect_ratio, const char *scalar_field, float threshold, const char *comparator);
 
 
    void setExternalTriangulation(const std::vector<helios::vec3> &triangle_vertices, const std::vector<int> &scanIDs);
 
 
 
    void addTrianglesToContext(helios::Context *context) const;
 
    // -------- GRID ----------- //
 
    uint getGridCellCount() const;
 
 
    void addGridCell(const helios::vec3 &center, const helios::vec3 &size, float rotation);
 
 
    void addGridCell(const helios::vec3 &center, const helios::vec3 &global_anchor, const helios::vec3 &size, const helios::vec3 &global_size, float rotation, const helios::int3 &global_ijk, const helios::int3 &global_count);
 
 
    helios::vec3 getCellCenter(uint index) const;
 
 
    helios::vec3 getCellGlobalAnchor(uint index) const;
 
 
    helios::vec3 getCellSize(uint index) const;
 
 
    float getCellRotation(uint index) const;
 
    // ------- SYNTHETIC SCAN ------ //
 
 
    void syntheticScan(helios::Context *context);
 
 
    void syntheticScan(helios::Context *context, bool append);
 
 
    void syntheticScan(helios::Context *context, bool scan_grid_only, bool record_misses);
 
 
    void syntheticScan(helios::Context *context, bool scan_grid_only, bool record_misses, bool append);
 
 
    void syntheticScan(helios::Context *context, int rays_per_pulse, float pulse_distance_threshold);
 
 
    void syntheticScan(helios::Context *context, int rays_per_pulse, float pulse_distance_threshold, bool append);
 
 
    void syntheticScan(helios::Context *context, int rays_per_pulse, float pulse_distance_threshold, bool scan_grid_only, bool record_misses);
 
 
    void syntheticScan(helios::Context *context, int rays_per_pulse, float pulse_distance_threshold, bool scan_grid_only, bool record_misses, bool append);
 
 
    void syntheticScan(helios::Context *context, int rays_per_pulse, float pulse_distance_threshold, ReturnMode return_mode, bool scan_grid_only = false, bool record_misses = false, bool append = true);
 
 
    void setSyntheticScanMemoryBudget(size_t bytes);
 
 
    [[nodiscard]] size_t getSyntheticScanMemoryBudget() const;
 
 
    std::vector<float> calculateSyntheticLeafArea(helios::Context *context);
 
 
    std::vector<float> calculateSyntheticGtheta(helios::Context *context);
 
    // -------- LEAF AREA -------- //
 
 
    void setCellLeafArea(float area, uint index);
 
 
    float getCellLeafArea(uint index) const;
 
 
    float getCellLeafAreaDensity(uint index) const;
 
 
    void setCellGtheta(float Gtheta, uint index);
 
 
    float getCellGtheta(uint index) const;
 
    // -------- LEAF AREA INVERSION UNCERTAINTY -------- //
    //
    // The following accessors expose the per-voxel statistical SAMPLING uncertainty of the
    // leaf-area inversion (Pimont et al. 2018, RSE 215:343-370). This is the uncertainty owing
    // to the finite number of beams that sampled the voxel and to vegetation-element position
    // variability. It is CONDITIONAL on the beams that entered the voxel and does NOT capture
    // occlusion/coverage bias (voxels shadowed so that beams never penetrate). The group form
    // (\ref getGroupLADConfidenceInterval()) is the recommended path: single-voxel intervals are
    // routinely +-50-100%, whereas group intervals (a vertical slice, a whole plant) are +-5-10%.
 
 
    int getCellBeamCount(uint index) const;
 
 
    float getCellRelativeDensityIndex(uint index) const;
 
 
    float getCellMeanPathLength(uint index) const;
 
 
    float getCellLADVariance(uint index) const;
 
 
    bool getCellLeafAreaConfidenceInterval(uint index, float confidence_level, float &lower, float &upper) const;
 
 
    bool getGroupLADConfidenceInterval(const std::vector<uint> &indices, float confidence_level, float &mean_lad, float &lower, float &upper) const;
 
 
    std::vector<helios::vec3> gapfillMisses();
 
 
    std::vector<helios::vec3> gapfillMisses(uint scanID);
 
 
    std::vector<helios::vec3> gapfillMisses(uint scanID, const bool gapfill_grid_only, const bool add_flags);
 
 
 
    void forceBruteForceLeafArea(bool force) {
        force_bruteforce_LAD = force;
    }
 
 
    void calculateLeafArea(helios::Context *context);
 
 
    void calculateLeafArea(helios::Context *context, int min_voxel_hits);
 
 
    void calculateLeafArea(helios::Context *context, int min_voxel_hits, float element_width);
 
 
    void calculateLeafArea(helios::Context *context, float Gtheta, int min_voxel_hits, float element_width);
 
 
    [[deprecated("Use calculateLeafArea() instead. GPU functionality is now provided by the CollisionDetection plugin.")]]
    void calculateLeafAreaGPU(helios::Context *context);
 
 
    [[deprecated("Use calculateLeafArea(context, min_voxel_hits) instead. GPU functionality is now provided by the CollisionDetection plugin.")]]
    void calculateLeafAreaGPU(helios::Context *context, int min_voxel_hits);
 
    void enableGPUAcceleration();
 
    void disableGPUAcceleration();
 
 
    [[nodiscard]] bool isGPUAvailable() const;
 
    [[nodiscard]] bool isGPUAccelerationEnabled() const;
 
    void calculateHitGridCell();
 
    // -------- RECONSTRUCTION --------- //
 
 
    void leafReconstructionAlphaMask(float minimum_leaf_group_area, float maximum_leaf_group_area, float leaf_aspect_ratio, const char *mask_file);
 
 
    void leafReconstructionAlphaMask(float minimum_leaf_group_area, float maximum_leaf_group_area, float leaf_aspect_ratio, float leaf_length_constant, const char *mask_file);
 
 
    void trunkReconstruction(const helios::vec3 &box_center, const helios::vec3 &box_size, float Lmax, float max_aspect_ratio);
 
 
    std::vector<uint> loadTreeQSM(helios::Context *context, const std::string &filename, uint radial_subdivisions, const std::string &texture_file = "");
 
 
    std::vector<uint> loadTreeQSMColormap(helios::Context *context, const std::string &filename, uint radial_subdivisions, const std::string &colormap_name);
 
 
    void cropBeamsToGridAngleRange(uint source);
 
 
    std::vector<uint> peakFinder(std::vector<float> signal);
};
 
bool sortcol0(const std::vector<double> &v0, const std::vector<double> &v1);
 
bool sortcol1(const std::vector<double> &v0, const std::vector<double> &v1);
 
#endif