26using namespace helios;
31 bool isStandardColumnToken(
const std::string &label) {
32 static const std::set<std::string> standard_tokens = {
"x",
"y",
"z",
"r",
"g",
"b",
"r255",
"g255",
"b255",
"row",
"column",
"zenith",
"azimuth",
"zenith_rad",
"azimuth_rad",
36 return standard_tokens.find(label) != standard_tokens.end();
47 if (
context->doesPrimitiveDataExist(UUID, label.c_str())) {
48 switch (
context->getPrimitiveDataType(label.c_str())) {
51 context->getPrimitiveData(UUID, label.c_str(), v);
56 context->getPrimitiveData(UUID, label.c_str(), value);
61 context->getPrimitiveData(UUID, label.c_str(), v);
67 context->getPrimitiveData(UUID, label.c_str(), v);
78 uint objID =
context->getPrimitiveParentObjectID(UUID);
79 if (
context->doesObjectExist(objID) &&
context->doesObjectDataExist(objID, label.c_str())) {
80 switch (
context->getObjectDataType(label.c_str())) {
83 context->getObjectData(objID, label.c_str(), v);
88 context->getObjectData(objID, label.c_str(), value);
93 context->getObjectData(objID, label.c_str(), v);
99 context->getObjectData(objID, label.c_str(), v);
114 double normalQuantile(
double p) {
115 if (p <= 0.0 || p >= 1.0) {
119 static const double a[] = {-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02, 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00};
120 static const double b[] = {-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02, 6.680131188771972e+01, -1.328068155288572e+01};
121 static const double c[] = {-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00, -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00};
122 static const double d[] = {7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00, 3.754408661907416e+00};
123 const double p_low = 0.02425;
124 const double p_high = 1.0 - p_low;
127 q = std::sqrt(-2.0 * std::log(p));
128 return (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0);
129 }
else if (p <= p_high) {
132 return (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q / (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1.0);
134 q = std::sqrt(-2.0 * std::log(1.0 - p));
135 return -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0);
158 double dot = double(q0.
x) * q1.
x + double(q0.
y) * q1.
y + double(q0.
z) * q1.
z + double(q0.
w) * q1.
w;
170 const double theta_0 = std::acos(dot);
171 const double theta = theta_0 * u;
172 const double sin_theta_0 = std::sin(theta_0);
173 const double s0 = std::sin(theta_0 - theta) / sin_theta_0;
174 const double s1 = std::sin(theta) / sin_theta_0;
181 helios::vec4 quat_from_rpy(
float roll,
float pitch,
float yaw) {
182 const float cr = std::cos(roll * 0.5f), sr = std::sin(roll * 0.5f);
183 const float cp = std::cos(pitch * 0.5f), sp = std::sin(pitch * 0.5f);
184 const float cy = std::cos(yaw * 0.5f), sy = std::sin(yaw * 0.5f);
187 q.
w = cr * cp * cy + sr * sp * sy;
188 q.
x = sr * cp * cy - cr * sp * sy;
189 q.
y = cr * sp * cy + sr * cp * sy;
190 q.
z = cr * cp * sy - sr * sp * cy;
198 inline bool triangulationCancelled(
volatile int *cancel_flag) {
199 return cancel_flag !=
nullptr && *cancel_flag != 0;
225ScanMetadata::ScanMetadata(
const vec3 &a_origin,
uint a_Ntheta,
float a_thetaMin,
float a_thetaMax,
uint a_Nphi,
float a_phiMin,
float a_phiMax,
float a_exitDiameter,
float a_beamDivergence,
float a_rangeNoiseStdDev,
float a_angleNoiseStdDev,
226 const vector<string> &a_columnFormat,
float a_scanTiltRoll,
float a_scanTiltPitch,
float a_scanAzimuthOffset) {
249ScanMetadata::ScanMetadata(
const vec3 &a_origin,
const std::vector<float> &a_beamZenithAngles,
uint a_Nphi,
float a_phiMin,
float a_phiMax,
float a_exitDiameter,
float a_beamDivergence,
float a_rangeNoiseStdDev,
float a_angleNoiseStdDev,
250 const vector<string> &a_columnFormat,
float a_scanTiltRoll,
float a_scanTiltPitch,
float a_scanAzimuthOffset) {
252 if (a_beamZenithAngles.empty()) {
253 helios_runtime_error(
"ERROR (ScanMetadata): A spinning multibeam scan requires at least one beam (channel) zenith angle, but the provided beamZenithAngles vector is empty.");
261 thetaMin = *std::min_element(a_beamZenithAngles.begin(), a_beamZenithAngles.end());
262 thetaMax = *std::max_element(a_beamZenithAngles.begin(), a_beamZenithAngles.end());
285 return {-v.x, -v.y, -v.z};
287 inline double dot(
const dvec3 &a,
const dvec3 &b) {
288 return a.x * b.x + a.y * b.y + a.z * b.z;
290 inline dvec3 normalize(
const dvec3 &v) {
291 double m = std::sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
292 return {v.x / m, v.y / m, v.z / m};
310 bool refractRay(
const dvec3 &incident,
const dvec3 &normal,
double n_from,
double n_to, dvec3 &refracted) {
311 const double r = n_from / n_to;
312 const double c = -dot(normal, incident);
313 const double radicand = 1.0 - r * r * (1.0 - c * c);
314 if (radicand < 0.0) {
317 const double k = r * c - std::sqrt(radicand);
318 refracted = normalize({r * incident.x + k * normal.x, r * incident.y + k * normal.y, r * incident.z + k * normal.z});
349 const double t = double(pulse_index) * scan.
pulse_period;
353 dvec3 beam = {0.0, 1.0, 0.0};
354 const dvec3 flat_normal = {0.0, -1.0, 0.0};
356 for (
const RisleyPrism &prism : scan.risley_prisms) {
357 const double phi = prism.phase + prism.rotor_rate * t;
358 const double sinW = std::sin(prism.wedge_angle);
359 const double cosW = std::cos(prism.wedge_angle);
365 const dvec3 tilted_normal = normalize({sinW * std::cos(phi), cosW, sinW * std::sin(phi)});
370 if (!refractRay(beam, flat_normal, n_air, prism.refractive_index, in_glass)) {
374 if (!refractRay(in_glass, -tilted_normal, prism.refractive_index, n_air, out_glass)) {
380 beam = normalize(beam);
393 return cart2sphere(risleyBodyDirection(*
this, column));
404 float elevation = 0.5f *
M_PI - zenith;
419 float theta = direction.
zenith;
442 }
else if (row >=
Ntheta) {
447 }
else if (column >=
Nphi) {
456 const size_t M =
traj_t.size();
458 helios_runtime_error(
"ERROR (ScanMetadata::poseAt): the scan has no trajectory samples. This scan was not created as a moving-platform scan (see LiDARcloud::addScanMoving).");
461 helios_runtime_error(
"ERROR (ScanMetadata::poseAt): trajectory arrays have inconsistent lengths (traj_t=" + std::to_string(M) +
", traj_pos=" + std::to_string(
traj_pos.size()) +
", traj_quat=" + std::to_string(
traj_quat.size()) +
462 "). All three must have the same number of samples.");
473 if (t <=
traj_t.front()) {
486 const auto upper = std::upper_bound(
traj_t.begin(),
traj_t.end(), t);
487 const size_t i1 = size_t(upper -
traj_t.begin());
488 const size_t i0 = i1 - 1;
490 const double t0_s =
traj_t.at(i0);
491 const double t1_s =
traj_t.at(i1);
492 const double denom = t1_s - t0_s;
494 helios_runtime_error(
"ERROR (ScanMetadata::poseAt): trajectory times are not strictly increasing (traj_t[" + std::to_string(i0) +
"]=" + std::to_string(t0_s) +
" >= traj_t[" + std::to_string(i1) +
"]=" + std::to_string(t1_s) +
").");
496 const double u = (t - t0_s) / denom;
506 hitgridcellcomputed =
false;
507 triangulationcomputed =
false;
508 triangulation_candidate_count = 0;
509 triangulation_dropped_lmax = 0;
510 triangulation_dropped_aspect = 0;
511 triangulation_dropped_degenerate = 0;
512 printmessages =
true;
513 collision_detection =
nullptr;
517 delete collision_detection;
521 printmessages =
false;
525 if (collision_detection ==
nullptr) {
537 if (collision_detection !=
nullptr) {
543 synthetic_scan_progress = ptr;
554void LiDARcloud::finishUnifiedRayTracing() {
558void LiDARcloud::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) {
560 constexpr uint MISS_UUID = 0xFFFFFFFFu;
562 if (total_rays == 0) {
569 std::vector<uint> uuid(total_rays);
570 std::vector<helios::vec3> normal(total_rays);
574 if (packet_size > 1) {
575 collision_detection->
castRaysSoA_packets(ray_origins, direction, total_rays, packet_size, miss_distance, hit_t, normal.data(), uuid.data());
577 collision_detection->
castRaysSoA(ray_origins, direction, total_rays, miss_distance, hit_t, normal.data(), uuid.data());
581 for (
size_t i = 0; i < total_rays; i++) {
582 if (uuid[i] != MISS_UUID) {
583 hit_ID[i] =
static_cast<int>(uuid[i]);
586 hit_fnorm[i] = ray_dir.
x * normal[i].x + ray_dir.
y * normal[i].y + ray_dir.
z * normal[i].z;
588 hit_t[i] = miss_distance;
599 prepareUnifiedRayTracing(
context);
600 castRaysUnified(N *
size_t(Npulse), ray_origins, direction, hit_t, hit_fnorm, hit_ID,
size_t(Npulse));
601 finishUnifiedRayTracing();
622 printmessages =
true;
626 progress_callback = std::move(callback);
631 helios_runtime_error(
"ERROR (LiDARcloud::setSyntheticScanMemoryBudget): the memory budget must be greater than zero.");
633 synthetic_scan_memory_budget_bytes = bytes;
637 return synthetic_scan_memory_budget_bytes;
640bool LiDARcloud::anyScanMoving()
const {
641 for (
const auto &scan: scans) {
649void LiDARcloud::validateRayDirections() {
653 if (anyScanMoving()) {
654 helios_runtime_error(
"ERROR (LiDARcloud::validateRayDirections): ray-direction validation is not supported for moving-platform scans (see addScanMoving), because the per-pulse origins make a single-origin direction check meaningless.");
671 if (err_theta > 1e-6 || err_phi > 1e-6) {
672 helios_runtime_error(
"ERROR (LiDARcloud::validateRayDirections): validation of ray directions failed.");
686 float epsilon = 1e-5;
689 std::cerr <<
"WARNING (LiDARcloud::addScan): Specified scan minimum zenith angle of " << newscan.
thetaMin <<
" is less than 0. Truncating to 0." << std::endl;
693 std::cerr <<
"WARNING (LiDARcloud::addScan): Specified scan minimum azimuth angle of " << newscan.
phiMin <<
" is less than 0. Truncating to 0." << std::endl;
697 std::cerr <<
"WARNING (LiDARcloud::addScan): Specified scan maximum zenith angle of " << newscan.
thetaMax <<
" is greater than pi. Setting thetaMin to 0 and truncating thetaMax to pi. Did you mistakenly use degrees instead of radians?"
706 std::cerr <<
"WARNING (LiDARcloud::addScan): Specified scan maximum azimuth angle of " << newscan.
phiMax <<
" is greater than 2pi. Did you mistakenly use degrees instead of radians?" << std::endl;
711 table.resize(newscan.
Ntheta, newscan.
Nphi, -1);
712 hit_tables.push_back(table);
714 scans.emplace_back(newscan);
716 return scans.size() - 1;
721 const size_t M = traj_t.size();
723 helios_runtime_error(
"ERROR (LiDARcloud::addScanMoving): the trajectory is empty. At least one pose sample is required.");
725 if (traj_pos.size() != M || traj_quat.size() != M) {
726 helios_runtime_error(
"ERROR (LiDARcloud::addScanMoving): trajectory arrays have inconsistent lengths (traj_t=" + std::to_string(M) +
", traj_pos=" + std::to_string(traj_pos.size()) +
", traj_quat=" + std::to_string(traj_quat.size()) +
727 "). All three must have the same number of samples.");
729 for (
size_t k = 1; k < M; k++) {
730 if (traj_t.at(k) <= traj_t.at(k - 1)) {
731 helios_runtime_error(
"ERROR (LiDARcloud::addScanMoving): trajectory times must be strictly increasing, but traj_t[" + std::to_string(k - 1) +
"]=" + std::to_string(traj_t.at(k - 1)) +
" >= traj_t[" + std::to_string(k) +
"]=" +
732 std::to_string(traj_t.at(k)) +
".");
738 for (
size_t k = 0; k < M; k++) {
739 if (!std::isfinite(traj_t.at(k))) {
740 helios_runtime_error(
"ERROR (LiDARcloud::addScanMoving): trajectory time traj_t[" + std::to_string(k) +
"] is not finite (NaN or infinity).");
743 if (!std::isfinite(p.
x) || !std::isfinite(p.
y) || !std::isfinite(p.
z)) {
744 helios_runtime_error(
"ERROR (LiDARcloud::addScanMoving): trajectory position traj_pos[" + std::to_string(k) +
"] is not finite (NaN or infinity).");
747 if (!std::isfinite(q.
x) || !std::isfinite(q.
y) || !std::isfinite(q.
z) || !std::isfinite(q.
w)) {
748 helios_runtime_error(
"ERROR (LiDARcloud::addScanMoving): trajectory quaternion traj_quat[" + std::to_string(k) +
"] is not finite (NaN or infinity).");
751 helios_runtime_error(
"ERROR (LiDARcloud::addScanMoving): trajectory quaternion traj_quat[" + std::to_string(k) +
"] has near-zero magnitude and cannot be normalized to a valid rotation.");
754 if (!std::isfinite(t0)) {
755 helios_runtime_error(
"ERROR (LiDARcloud::addScanMoving): t0 must be finite, but a non-finite value was provided.");
757 if (pulse_rate_hz <= 0.f) {
758 helios_runtime_error(
"ERROR (LiDARcloud::addScanMoving): pulse_rate_hz must be greater than 0, but " + std::to_string(pulse_rate_hz) +
" was provided.");
762 helios_runtime_error(
"ERROR (LiDARcloud::addScanMoving): the scan specifies a non-zero static tilt (scanTilt_roll/pitch/azimuth), which is not applied for moving-platform scans. Platform attitude must be supplied entirely "
763 "through the trajectory quaternions and the boresight; leave the static tilt at zero.");
775 scan.
origin = traj_pos.front();
785 if (traj_rpy.size() != traj_t.size()) {
786 helios_runtime_error(
"ERROR (LiDARcloud::addScanMoving): trajectory arrays have inconsistent lengths (traj_t=" + std::to_string(traj_t.size()) +
", traj_rpy=" + std::to_string(traj_rpy.size()) +
"). All trajectory arrays must have the same number of samples.");
789 std::vector<vec4> traj_quat;
790 traj_quat.reserve(traj_rpy.size());
791 for (
const vec3 &rpy: traj_rpy) {
792 traj_quat.push_back(quat_from_rpy(rpy.x, rpy.y, rpy.z));
795 return addScanMoving(scan, traj_t, traj_pos, traj_quat, lever_arm, boresight_rpy, pulse_rate_hz, t0);
798uint LiDARcloud::addScanSpinning(
const std::vector<float> &beamElevationAngles,
float azimuthStep_rad,
float pulse_rate_hz,
const std::vector<double> &traj_t,
const std::vector<vec3> &traj_pos,
const std::vector<vec4> &traj_quat,
799 const vec3 &lever_arm,
const vec3 &boresight_rpy,
float exitDiameter,
float beamDivergence,
float rangeNoiseStdDev,
float angleNoiseStdDev,
const std::vector<std::string> &columnFormat,
double t0) {
806 if (beamElevationAngles.empty()) {
807 helios_runtime_error(
"ERROR (LiDARcloud::addScanSpinning): beamElevationAngles is empty. A spinning multibeam sensor requires at least one channel.");
809 if (azimuthStep_rad <= 0.f) {
810 helios_runtime_error(
"ERROR (LiDARcloud::addScanSpinning): azimuthStep_rad must be greater than 0, but " + std::to_string(azimuthStep_rad) +
" was provided.");
812 if (pulse_rate_hz <= 0.f) {
813 helios_runtime_error(
"ERROR (LiDARcloud::addScanSpinning): pulse_rate_hz must be greater than 0, but " + std::to_string(pulse_rate_hz) +
" was provided.");
815 const size_t M = traj_t.size();
817 helios_runtime_error(
"ERROR (LiDARcloud::addScanSpinning): the trajectory is empty. At least one pose sample is required (a stationary capture is expressed as two coincident poses separated by the acquisition duration).");
822 if (traj_pos.size() != M || traj_quat.size() != M) {
823 helios_runtime_error(
"ERROR (LiDARcloud::addScanSpinning): trajectory arrays have inconsistent lengths (traj_t=" + std::to_string(M) +
", traj_pos=" + std::to_string(traj_pos.size()) +
", traj_quat=" + std::to_string(traj_quat.size()) +
824 "). All trajectory arrays must have the same number of samples.");
828 std::vector<float> beamZenithAngles;
829 beamZenithAngles.reserve(beamElevationAngles.size());
830 for (
float elevation: beamElevationAngles) {
831 beamZenithAngles.push_back(0.5f *
float(
M_PI) - elevation);
834 const uint channels =
uint(beamZenithAngles.size());
837 const uint steps_per_rev =
uint(std::lround(2.0 *
M_PI /
double(azimuthStep_rad)));
838 if (steps_per_rev == 0) {
839 helios_runtime_error(
"ERROR (LiDARcloud::addScanSpinning): azimuthStep_rad=" + std::to_string(azimuthStep_rad) +
" is larger than 2pi, which yields zero azimuth steps per revolution. Use a finer azimuth resolution.");
842 const double duration = traj_t.back() - traj_t.front();
843 if (duration <= 0.0) {
844 helios_runtime_error(
"ERROR (LiDARcloud::addScanSpinning): the trajectory duration (traj_t.back() - traj_t.front() = " + std::to_string(duration) +
") must be greater than 0.");
847 const double rotation_rate = double(pulse_rate_hz) / (double(channels) * double(steps_per_rev));
848 const double n_revolutions = rotation_rate * duration;
849 const uint Nphi =
uint(std::lround(
double(steps_per_rev) * n_revolutions));
851 helios_runtime_error(
"ERROR (LiDARcloud::addScanSpinning): the derived azimuth-step count is zero (PRF=" + std::to_string(pulse_rate_hz) +
" Hz over a " + std::to_string(duration) +
852 " s trajectory yields less than one azimuth step). Increase the PRF, the trajectory duration, or the azimuth resolution.");
857 ScanMetadata scan(traj_pos.front(), beamZenithAngles, Nphi, 0.f,
float(n_revolutions * 2.0 *
M_PI), exitDiameter, beamDivergence, rangeNoiseStdDev, angleNoiseStdDev, columnFormat);
863 return addScanMoving(scan, traj_t, traj_pos, traj_quat, lever_arm, boresight_rpy, pulse_rate_hz, t0);
866uint LiDARcloud::addScanSpinning(
const std::vector<float> &beamElevationAngles,
float azimuthStep_rad,
float pulse_rate_hz,
const std::vector<double> &traj_t,
const std::vector<vec3> &traj_pos,
const std::vector<vec3> &traj_rpy,
867 const vec3 &lever_arm,
const vec3 &boresight_rpy,
float exitDiameter,
float beamDivergence,
float rangeNoiseStdDev,
float angleNoiseStdDev,
const std::vector<std::string> &columnFormat,
double t0) {
870 if (traj_rpy.size() != traj_t.size()) {
871 helios_runtime_error(
"ERROR (LiDARcloud::addScanSpinning): trajectory arrays have inconsistent lengths (traj_t=" + std::to_string(traj_t.size()) +
", traj_rpy=" + std::to_string(traj_rpy.size()) +
"). All trajectory arrays must have the same number of samples.");
874 std::vector<vec4> traj_quat;
875 traj_quat.reserve(traj_rpy.size());
876 for (
const vec3 &rpy: traj_rpy) {
877 traj_quat.push_back(quat_from_rpy(rpy.x, rpy.y, rpy.z));
880 return addScanSpinning(beamElevationAngles, azimuthStep_rad, pulse_rate_hz, traj_t, traj_pos, traj_quat, lever_arm, boresight_rpy, exitDiameter, beamDivergence, rangeNoiseStdDev, angleNoiseStdDev, columnFormat, t0);
884 const std::vector<vec4> &traj_quat,
const vec3 &lever_arm,
const vec3 &boresight_rpy,
float exitDiameter,
float beamDivergence,
float rangeNoiseStdDev,
float angleNoiseStdDev,
885 const std::vector<std::string> &columnFormat,
double t0) {
890 ScanMetadata scan(traj_pos.empty() ?
make_vec3(0, 0, 0) : traj_pos.front(), Ntheta, thetaMin, thetaMax, Nphi, phiMin, phiMax, exitDiameter, beamDivergence, rangeNoiseStdDev, angleNoiseStdDev, columnFormat);
893 return addScanMoving(scan, traj_t, traj_pos, traj_quat, lever_arm, boresight_rpy, pulse_rate_hz, t0);
896uint LiDARcloud::addScanRisley(
const std::vector<RisleyPrism> &prisms,
double refractive_index_air,
float pulse_rate_hz,
const std::vector<double> &traj_t,
const std::vector<vec3> &traj_pos,
const std::vector<vec4> &traj_quat,
897 const vec3 &lever_arm,
const vec3 &boresight_rpy,
float exitDiameter,
float beamDivergence,
float rangeNoiseStdDev,
float angleNoiseStdDev,
const std::vector<std::string> &columnFormat,
double t0) {
904 if (prisms.empty()) {
905 helios_runtime_error(
"ERROR (LiDARcloud::addScanRisley): the prism stack is empty. A Risley-prism scanner requires at least one rotating wedge prism (a Livox-style sensor uses two counter-rotating prisms).");
907 if (pulse_rate_hz <= 0.f) {
908 helios_runtime_error(
"ERROR (LiDARcloud::addScanRisley): pulse_rate_hz must be greater than 0, but " + std::to_string(pulse_rate_hz) +
" was provided.");
910 if (refractive_index_air <= 0.0) {
911 helios_runtime_error(
"ERROR (LiDARcloud::addScanRisley): refractive_index_air must be greater than 0, but " + std::to_string(refractive_index_air) +
" was provided.");
913 for (
size_t k = 0; k < prisms.size(); k++) {
914 if (prisms.at(k).refractive_index <= 0.0) {
915 helios_runtime_error(
"ERROR (LiDARcloud::addScanRisley): prism " + std::to_string(k) +
" has a non-positive refractive index (" + std::to_string(prisms.at(k).refractive_index) +
").");
918 const size_t M = traj_t.size();
920 helios_runtime_error(
"ERROR (LiDARcloud::addScanRisley): the trajectory is empty. At least one pose sample is required (a stationary capture is expressed as two coincident poses separated by the acquisition duration).");
925 if (traj_pos.size() != M || traj_quat.size() != M) {
926 helios_runtime_error(
"ERROR (LiDARcloud::addScanRisley): trajectory arrays have inconsistent lengths (traj_t=" + std::to_string(M) +
", traj_pos=" + std::to_string(traj_pos.size()) +
", traj_quat=" + std::to_string(traj_quat.size()) +
927 "). All trajectory arrays must have the same number of samples.");
930 const double duration = traj_t.back() - traj_t.front();
931 if (duration <= 0.0) {
932 helios_runtime_error(
"ERROR (LiDARcloud::addScanRisley): the trajectory duration (traj_t.back() - traj_t.front() = " + std::to_string(duration) +
") must be greater than 0.");
938 const uint Npulses =
uint(std::lround(
double(pulse_rate_hz) * duration));
940 helios_runtime_error(
"ERROR (LiDARcloud::addScanRisley): the derived pulse count is zero (PRF=" + std::to_string(pulse_rate_hz) +
" Hz over a " + std::to_string(duration) +
941 " s trajectory yields less than one pulse). Increase the PRF or the trajectory duration.");
948 ScanMetadata scan(traj_pos.front(), 1u, 0.f,
float(
M_PI), Npulses, 0.f,
float(2.0 *
M_PI), exitDiameter, beamDivergence, rangeNoiseStdDev, angleNoiseStdDev, columnFormat);
957 const size_t Nsample = std::min<size_t>(Npulses, 20000);
958 float max_halfangle = 0.f;
959 for (
size_t k = 0; k < Nsample; k++) {
961 float halfangle = std::acos(std::max(-1.f, std::min(1.f, dir.
y)));
962 if (halfangle > max_halfangle) {
963 max_halfangle = halfangle;
966 scan.
thetaMin = std::max(0.f, 0.5f *
float(
M_PI) - max_halfangle);
969 return addScanMoving(scan, traj_t, traj_pos, traj_quat, lever_arm, boresight_rpy, pulse_rate_hz, t0);
972uint LiDARcloud::addScanRisley(
const std::vector<RisleyPrism> &prisms,
double refractive_index_air,
float pulse_rate_hz,
const std::vector<double> &traj_t,
const std::vector<vec3> &traj_pos,
const std::vector<vec3> &traj_rpy,
973 const vec3 &lever_arm,
const vec3 &boresight_rpy,
float exitDiameter,
float beamDivergence,
float rangeNoiseStdDev,
float angleNoiseStdDev,
const std::vector<std::string> &columnFormat,
double t0) {
976 if (traj_rpy.size() != traj_t.size()) {
977 helios_runtime_error(
"ERROR (LiDARcloud::addScanRisley): trajectory arrays have inconsistent lengths (traj_t=" + std::to_string(traj_t.size()) +
", traj_rpy=" + std::to_string(traj_rpy.size()) +
"). All trajectory arrays must have the same number of samples.");
980 std::vector<vec4> traj_quat;
981 traj_quat.reserve(traj_rpy.size());
982 for (
const vec3 &rpy: traj_rpy) {
983 traj_quat.push_back(quat_from_rpy(rpy.x, rpy.y, rpy.z));
986 return addScanRisley(prisms, refractive_index_air, pulse_rate_hz, traj_t, traj_pos, traj_quat, lever_arm, boresight_rpy, exitDiameter, beamDivergence, rangeNoiseStdDev, angleNoiseStdDev, columnFormat, t0);
995 std::map<std::string, double> data;
1011 std::map<std::string, double> data;
1016size_t LiDARcloud::getOrCreateHitDataColumn(
const std::string &label) {
1017 auto it = hit_data_label_index.find(label);
1018 if (it != hit_data_label_index.end()) {
1025 const size_t slot = hit_data_labels.size();
1026 hit_data_labels.push_back(label);
1027 hit_data_label_index[label] = slot;
1028 hit_data_columns.emplace_back(hits.size(), 0.0);
1029 hit_data_present.emplace_back(hits.size(),
char(0));
1035void LiDARcloud::appendHitData(
const std::map<std::string, double> &data) {
1036 const size_t i = hits.size() - 1;
1039 for (
size_t s = 0; s < hit_data_columns.size(); s++) {
1040 hit_data_columns[s].push_back(0.0);
1041 hit_data_present[s].push_back(
char(0));
1047 for (
const auto &kv: data) {
1048 const size_t s = getOrCreateHitDataColumn(kv.first);
1049 hit_data_columns[s][i] = kv.second;
1050 hit_data_present[s][i] = char(1);
1057 if (scanID >= scans.size()) {
1058 helios_runtime_error(
"ERROR (LiDARcloud::addHitPoint): Hit point cannot be added to scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1064 HitPoint hit(scanID, xyz, direction, row_column, color);
1066 hits.push_back(hit);
1067 appendHitData(data);
1075 HitPoint hit(scanID, xyz, direction, row_column, color);
1077 hits.push_back(hit);
1078 appendHitData(data);
1083 if (index >= hits.size()) {
1084 cerr <<
"WARNING (deleteHitPoint): Hit point #" << index <<
" cannot be deleted from the scan because there have only been " << hits.size() <<
" hit points added." << endl;
1091 const size_t last = hits.size() - 1;
1092 for (
size_t s = 0; s < hit_data_columns.size(); s++) {
1093 std::swap(hit_data_columns[s][index], hit_data_columns[s][last]);
1094 std::swap(hit_data_present[s][index], hit_data_present[s][last]);
1095 hit_data_columns[s].pop_back();
1096 hit_data_present[s].pop_back();
1099 std::swap(hits.at(index), hits.back());
1108 if (scanID >= scans.size()) {
1109 helios_runtime_error(
"ERROR (LiDARcloud::getScanOrigin): Cannot get origin of scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1111 return scans.at(scanID).origin;
1115 if (scanID >= scans.size()) {
1116 helios_runtime_error(
"ERROR (LiDARcloud::getScanSizeTheta): Cannot get theta size for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1118 return scans.at(scanID).Ntheta;
1122 if (scanID >= scans.size()) {
1123 helios_runtime_error(
"ERROR (LiDARcloud::getScanSizePhi): Cannot get phi size for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1125 return scans.at(scanID).Nphi;
1129 if (scanID >= scans.size()) {
1130 helios_runtime_error(
"ERROR (LiDARcloud::getScanRangeTheta): Cannot get theta range for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1136 if (scanID >= scans.size()) {
1137 helios_runtime_error(
"ERROR (LiDARcloud::getScanRangePhi): Cannot get phi range for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1143 if (scanID >= scans.size()) {
1144 helios_runtime_error(
"ERROR (LiDARcloud::getScanBeamExitDiameter): Cannot get exit diameter for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1146 return scans.at(scanID).exitDiameter;
1150 if (scanID >= scans.size()) {
1151 helios_runtime_error(
"ERROR (LiDARcloud::getScanBeamDivergence): Cannot get beam divergence for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1153 return scans.at(scanID).beamDivergence;
1157 if (scanID >= scans.size()) {
1158 helios_runtime_error(
"ERROR (LiDARcloud::getScanRangeNoiseStdDev): Cannot get range noise standard deviation for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1160 return scans.at(scanID).rangeNoiseStdDev;
1164 if (scanID >= scans.size()) {
1165 helios_runtime_error(
"ERROR (LiDARcloud::getScanAngleNoiseStdDev): Cannot get angular noise standard deviation for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1167 return scans.at(scanID).angleNoiseStdDev;
1171 if (scanID >= scans.size()) {
1172 helios_runtime_error(
"ERROR (LiDARcloud::getScanReturnMode): Cannot get return mode for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1174 return scans.at(scanID).returnMode;
1178 if (scanID >= scans.size()) {
1179 helios_runtime_error(
"ERROR (LiDARcloud::setScanReturnMode): Cannot set return mode for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1181 scans.at(scanID).returnMode = returnMode;
1185 if (scanID >= scans.size()) {
1186 helios_runtime_error(
"ERROR (LiDARcloud::getScanSingleReturnSelection): Cannot get single-return selection for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1188 return scans.at(scanID).singleReturnSelection;
1192 if (scanID >= scans.size()) {
1193 helios_runtime_error(
"ERROR (LiDARcloud::setScanSingleReturnSelection): Cannot set single-return selection for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1195 scans.at(scanID).singleReturnSelection = selection;
1199 if (scanID >= scans.size()) {
1200 helios_runtime_error(
"ERROR (LiDARcloud::getScanMaxReturns): Cannot get maximum returns for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1202 return scans.at(scanID).maxReturns;
1206 if (scanID >= scans.size()) {
1207 helios_runtime_error(
"ERROR (LiDARcloud::setScanMaxReturns): Cannot set maximum returns for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1209 if (maxReturns < 1) {
1210 helios_runtime_error(
"ERROR (LiDARcloud::setScanMaxReturns): Maximum returns must be at least 1, but " + std::to_string(maxReturns) +
" was given.");
1212 scans.at(scanID).maxReturns = maxReturns;
1216 if (scanID >= scans.size()) {
1217 helios_runtime_error(
"ERROR (LiDARcloud::getScanPulseWidth): Cannot get pulse width for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1219 return scans.at(scanID).pulseWidth;
1223 if (scanID >= scans.size()) {
1224 helios_runtime_error(
"ERROR (LiDARcloud::setScanPulseWidth): Cannot set pulse width for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1226 if (pulseWidth < 0.f) {
1227 helios_runtime_error(
"ERROR (LiDARcloud::setScanPulseWidth): Pulse width must be non-negative, but " + std::to_string(pulseWidth) +
" was given.");
1229 scans.at(scanID).pulseWidth = pulseWidth;
1233 if (scanID >= scans.size()) {
1234 helios_runtime_error(
"ERROR (LiDARcloud::getScanDetectionThreshold): Cannot get detection threshold for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1236 return scans.at(scanID).detectionThreshold;
1240 if (scanID >= scans.size()) {
1241 helios_runtime_error(
"ERROR (LiDARcloud::setScanDetectionThreshold): Cannot set detection threshold for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1243 if (detectionThreshold < 0.f) {
1244 helios_runtime_error(
"ERROR (LiDARcloud::setScanDetectionThreshold): Detection threshold must be non-negative, but " + std::to_string(detectionThreshold) +
" was given.");
1246 scans.at(scanID).detectionThreshold = detectionThreshold;
1250 if (scanID >= scans.size()) {
1251 helios_runtime_error(
"ERROR (LiDARcloud::getScanTiltRoll): Cannot get scanner tilt roll for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1253 return scans.at(scanID).scanTilt_roll;
1257 if (scanID >= scans.size()) {
1258 helios_runtime_error(
"ERROR (LiDARcloud::getScanTiltPitch): Cannot get scanner tilt pitch for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1260 return scans.at(scanID).scanTilt_pitch;
1264 if (scanID >= scans.size()) {
1265 helios_runtime_error(
"ERROR (LiDARcloud::getScanAzimuthOffset): Cannot get scanner azimuth offset for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1267 return scans.at(scanID).scanTilt_azimuth;
1271 if (scanID >= scans.size()) {
1272 helios_runtime_error(
"ERROR (LiDARcloud::getScanColumnFormat): Cannot get column format for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1274 return scans.at(scanID).columnFormat;
1278 if (scanID >= scans.size()) {
1279 helios_runtime_error(
"ERROR (LiDARcloud::getScanPattern): Cannot get scan pattern for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1281 return scans.at(scanID).scanPattern;
1285 if (scanID >= scans.size()) {
1286 helios_runtime_error(
"ERROR (LiDARcloud::getScanBeamZenithAngles): Cannot get beam zenith angles for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1288 return scans.at(scanID).beamZenithAngles;
1292 if (scanID >= scans.size()) {
1293 helios_runtime_error(
"ERROR (LiDARcloud::getScanMode): Cannot get scan mode for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1295 return scans.at(scanID).scanMode;
1299 if (scanID >= scans.size()) {
1300 helios_runtime_error(
"ERROR (LiDARcloud::getScanStepsPerRev): Cannot get steps per revolution for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1302 return scans.at(scanID).steps_per_rev;
1306 if (scanID >= scans.size()) {
1307 helios_runtime_error(
"ERROR (LiDARcloud::getScanRotationRate): Cannot get rotation rate for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1309 return scans.at(scanID).rotation_rate;
1313 if (scanID >= scans.size()) {
1314 helios_runtime_error(
"ERROR (LiDARcloud::getScanRevolutions): Cannot get revolution count for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1316 return scans.at(scanID).n_revolutions;
1320 if (scanID >= scans.size()) {
1321 helios_runtime_error(
"ERROR (LiDARcloud::getScanRisleyPrisms): Cannot get Risley prisms for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1323 return scans.at(scanID).risley_prisms;
1327 if (scanID >= scans.size()) {
1328 helios_runtime_error(
"ERROR (LiDARcloud::getScanRisleyRefractiveIndexAir): Cannot get the Risley air refractive index for scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1330 return scans.at(scanID).risley_refractive_index_air;
1335 if (index >= hits.size()) {
1336 helios_runtime_error(
"ERROR (LiDARcloud::getHitXYZ): Hit point index out of bounds. Requesting hit #" + std::to_string(index) +
" but scan only has " + std::to_string(hits.size()) +
" hits.");
1339 return hits.at(index).position;
1344 if (index >= hits.size()) {
1345 helios_runtime_error(
"ERROR (LiDARcloud::getHitOrigin): Hit point index out of bounds. Requesting hit #" + std::to_string(index) +
" but scan only has " + std::to_string(hits.size()) +
" hits.");
1359 if (index >= hits.size()) {
1360 helios_runtime_error(
"ERROR (LiDARcloud::getHitRaydir): Hit point index out of bounds. Requesting hit #" + std::to_string(index) +
" but scan only has " + std::to_string(hits.size()) +
" hits.");
1371 if (index >= hits.size()) {
1372 helios_runtime_error(
"ERROR (LiDARcloud::setHitScalarData): Hit point index out of bounds. Tried to set hit #" + std::to_string(index) +
" but scan only has " + std::to_string(hits.size()) +
" hits.");
1377 const size_t slot = getOrCreateHitDataColumn(label);
1378 hit_data_columns[slot][index] = value;
1379 hit_data_present[slot][index] = char(1);
1384 if (index >= hits.size()) {
1385 helios_runtime_error(
"ERROR (LiDARcloud::getHitData): Hit point index out of bounds. Requesting hit #" + std::to_string(index) +
" but scan only has " + std::to_string(hits.size()) +
" hits.");
1390 auto it = hit_data_label_index.find(label);
1391 if (it == hit_data_label_index.end() || hit_data_present[it->second][index] ==
char(0)) {
1392 helios_runtime_error(
"ERROR (LiDARcloud::getHitData): Data value ``" + std::string(label) +
"'' does not exist.");
1395 return hit_data_columns[it->second][index];
1400 if (index >= hits.size()) {
1404 auto it = hit_data_label_index.find(label);
1405 return it != hit_data_label_index.end() && hit_data_present[it->second][index] != char(0);
1409 auto it = hit_data_label_index.find(label);
1410 if (it == hit_data_label_index.end()) {
1413 return int(it->second);
1417 const size_t N = hits.size();
1420 auto it = hit_data_label_index.find(label);
1421 if (it == hit_data_label_index.end()) {
1423 std::fill(data.begin(), data.end(), absent_value);
1428 const std::vector<double> &column = hit_data_columns[it->second];
1429 const std::vector<char> &present = hit_data_present[it->second];
1430 for (
size_t i = 0; i < N; i++) {
1431 data[i] = present[i] != char(0) ? column[i] : absent_value;
1435void LiDARcloud::clearHits() {
1437 hit_data_labels.clear();
1438 hit_data_label_index.clear();
1439 hit_data_columns.clear();
1440 hit_data_present.clear();
1445 if (index >= hits.size()) {
1446 helios_runtime_error(
"ERROR (LiDARcloud::getHitColor): Hit point index out of bounds. Requesting hit #" + std::to_string(index) +
" but scan only has " + std::to_string(hits.size()) +
" hits.");
1449 return hits.at(index).color;
1454 if (index >= hits.size()) {
1455 helios_runtime_error(
"ERROR (LiDARcloud::getHitColor): Hit point index out of bounds. Requesting hit #" + std::to_string(index) +
" but scan only has " + std::to_string(hits.size()) +
" hits.");
1458 return hits.at(index).scanID;
1463 if (scanID >= scans.size()) {
1464 helios_runtime_error(
"ERROR (LiDARcloud::deleteHitPoint): Hit point cannot be deleted from scan #" + std::to_string(scanID) +
" because there have only been " + std::to_string(scans.size()) +
" scans added.");
1469 helios_runtime_error(
"ERROR (LiDARcloud::getHitIndex): Column in scan data table out of range.");
1472 int hit = hit_tables.at(scanID).get(row, column);
1481 if (index >= hits.size()) {
1482 helios_runtime_error(
"ERROR (LiDARcloud::getHitGridCell): Hit point index out of bounds. Requesting hit #" + std::to_string(index) +
" but scan only has " + std::to_string(hits.size()) +
" hits.");
1483 }
else if (hits.at(index).gridcell == -2) {
1484 cerr <<
"WARNING (LiDARcloud::getHitGridCell): hit grid cell for point #" << index <<
" was never set. Returning a value of `-1'. Did you forget to call calculateHitGridCell[*] first?" << endl;
1488 return hits.at(index).gridcell;
1493 if (index >= hits.size()) {
1494 helios_runtime_error(
"ERROR (LiDARcloud::setHitGridCell): Hit point index out of bounds. Tried to set hit #" + std::to_string(index) +
" but scan only has " + std::to_string(hits.size()) +
" hits.");
1497 hits.at(index).gridcell = cell;
1522 for (
auto &scan: scans) {
1526 for (
size_t i = 0; i < hits.size(); i++) {
1527 hits[i].position = hits[i].position + shift;
1528 transformHitOrigin(
uint(i), [&](
const vec3 &o) {
return o + shift; });
1534 if (scanID >= scans.size()) {
1535 helios_runtime_error(
"ERROR (LiDARcloud::coordinateShift): Cannot apply coordinate shift to scan " + std::to_string(scanID) +
" because it does not exist.");
1538 scans.at(scanID).origin = scans.at(scanID).origin + shift;
1540 for (
size_t i = 0; i < hits.size(); i++) {
1541 if (hits[i].scanID == scanID) {
1542 hits[i].position = hits[i].position + shift;
1543 transformHitOrigin(
uint(i), [&](
const vec3 &o) {
return o + shift; });
1550 for (
auto &scan: scans) {
1554 for (
size_t i = 0; i < hits.size(); i++) {
1555 hits[i].position =
rotatePoint(hits[i].position, rotation);
1558 hits[i].direction =
cart2sphere(hits[i].position - hitOriginOrFallback(
uint(i), scans.at(hits[i].scanID).origin));
1564 if (scanID >= scans.size()) {
1565 helios_runtime_error(
"ERROR (LiDARcloud::coordinateRotation): Cannot apply rotation to scan " + std::to_string(scanID) +
" because it does not exist.");
1568 scans.at(scanID).origin =
rotatePoint(scans.at(scanID).origin, rotation);
1570 for (
size_t i = 0; i < hits.size(); i++) {
1571 if (hits[i].scanID == scanID) {
1572 hits[i].position =
rotatePoint(hits[i].position, rotation);
1574 hits[i].direction =
cart2sphere(hits[i].position - hitOriginOrFallback(
uint(i), scans.at(scanID).origin));
1581 for (
auto &scan: scans) {
1585 for (
size_t i = 0; i < hits.size(); i++) {
1586 hits[i].position =
rotatePointAboutLine(hits[i].position, line_base, line_direction, rotation);
1588 hits[i].direction =
cart2sphere(hits[i].position - hitOriginOrFallback(
uint(i), scans.at(hits[i].scanID).origin));
1593 return triangles.size();
1597 return triangulation_candidate_count;
1601 return triangulation_dropped_lmax;
1605 return triangulation_dropped_aspect;
1609 return triangulation_dropped_degenerate;
1613 if (index >= triangles.size()) {
1614 helios_runtime_error(
"ERROR (LiDARcloud::getTriangle): Triangle index out of bounds. Tried to get triangle #" + std::to_string(index) +
" but point cloud only has " + std::to_string(triangles.size()) +
" triangles.");
1617 return triangles.at(index);
1626 if (printmessages && scans.size() == 0) {
1627 std::cout <<
"WARNING (LiDARcloud::addHitsToVisualizer): There are no scans in the point cloud, and thus there is no geometry to add...skipping." << std::endl;
1640 if (printmessages && scans.size() == 0) {
1641 std::cout <<
"WARNING (LiDARcloud::addHitsToVisualizer): There are no scans in the point cloud, and thus there is no geometry to add...skipping." << std::endl;
1647 float maxval = -1e9;
1648 if (strcmp(color_value,
"gridcell") == 0) {
1651 }
else if (strcmp(color_value,
"") != 0) {
1654 float data = float(
getHitData(i, color_value));
1655 if (data < minval) {
1658 if (data > maxval) {
1667 if (minval != 1e9 && maxval != -1e9) {
1668 cmap.setRange(minval, maxval);
1673 if (strcmp(color_value,
"") == 0) {
1675 }
else if (strcmp(color_value,
"gridcell") == 0) {
1685 float data = float(
getHitData(i, color_value));
1686 color = cmap.query(data);
1698 if (printmessages && scans.size() == 0) {
1699 std::cout <<
"WARNING (LiDARcloud::addGridToVisualizer): There are no scans in the point cloud, and thus there is no geometry to add...skipping." << std::endl;
1704 float maxval = -1e9;
1707 if (data < minval) {
1710 if (data > maxval) {
1716 if (minval != 1e9 && maxval != -1e9) {
1717 cmap.setRange(minval, maxval);
1743 vec3 boxmin, boxmax;
1746 float R = 2.f * sqrt(pow(boxmax.
x - boxmin.
x, 2) + pow(boxmax.
y - boxmin.
y, 2) + pow(boxmax.
z - boxmin.
z, 2));
1751 if (printmessages && scans.size() == 0) {
1752 std::cout <<
"WARNING (LiDARcloud::addGeometryToVisualizer): There are no scans in the point cloud, and thus there is no geometry to add...skipping." << std::endl;
1756 for (
uint i = 0; i < triangles.size(); i++) {
1766 if (printmessages && scans.size() == 0) {
1767 std::cout <<
"WARNING (LiDARcloud::addTrianglesToVisualizer): There are no scans in the point cloud, and thus there is no geometry to add...skipping." << std::endl;
1771 for (
uint i = 0; i < triangles.size(); i++) {
1775 if (tri.gridcell == gridcell) {
1782 if (size.
x <= 0 || size.
y <= 0 || size.z <= 0) {
1783 cerr <<
"failed.\n";
1787 if (ndiv.
x <= 0 || ndiv.
y <= 0 || ndiv.
z <= 0) {
1788 cerr <<
"failed.\n";
1789 helios_runtime_error(
"ERROR (LiDARcloud::addGrid): The number of grid cells in each direction must be positive.");
1793 vec3 gsubsize =
make_vec3(
float(size.
x) /
float(ndiv.
x),
float(size.
y) /
float(ndiv.
y),
float(size.z) /
float(ndiv.
z));
1797 for (
int k = 0; k < ndiv.
z; k++) {
1798 z = -0.5f * float(size.z) + (float(k) + 0.5f) *
float(gsubsize.
z);
1799 for (
int j = 0; j < ndiv.
y; j++) {
1800 y = -0.5f * float(size.
y) + (float(j) + 0.5f) *
float(gsubsize.
y);
1801 for (
int i = 0; i < ndiv.
x; i++) {
1802 x = -0.5f * float(size.
x) + (float(i) + 0.5f) *
float(gsubsize.
x);
1808 if (printmessages) {
1809 cout <<
"Adding grid cell #" << count <<
" with center " << subcenter_rot.
x + center.
x <<
"," << subcenter_rot.
y + center.
y <<
"," << subcenter.
z + center.
z <<
" and size " << gsubsize.
x <<
" x " << gsubsize.
y <<
" x "
1810 << gsubsize.
z << endl;
1829 boxmin =
make_vec3(center.
x - 0.5 * size.
x, center.
y - 0.5 * size.
y, center.
z - 0.5 * size.z);
1830 boxmax =
make_vec3(center.
x + 0.5 * size.
x, center.
y + 0.5 * size.
y, center.
z + 0.5 * size.z);
1854 size_t Ngroups = reconstructed_triangles.size();
1856 std::vector<helios::RGBcolor> ctable;
1857 std::vector<float> clocs;
1859 ctable.push_back(RGB::violet);
1860 ctable.push_back(RGB::blue);
1861 ctable.push_back(RGB::green);
1862 ctable.push_back(RGB::yellow);
1863 ctable.push_back(RGB::orange);
1864 ctable.push_back(RGB::red);
1866 clocs.push_back(0.f);
1867 clocs.push_back(0.2f);
1868 clocs.push_back(0.4f);
1869 clocs.push_back(0.6f);
1870 clocs.push_back(0.8f);
1871 clocs.push_back(1.f);
1873 Colormap colormap(ctable, clocs, 100, 0, Ngroups - 1);
1875 for (
size_t g = 0; g < Ngroups; g++) {
1877 float randi =
randu() * (Ngroups - 1);
1878 RGBcolor color = colormap.query(randi);
1880 for (
size_t t = 0; t < reconstructed_triangles.at(g).size(); t++) {
1882 helios::vec3 v0 = reconstructed_triangles.at(g).at(t).vertex0;
1883 helios::vec3 v1 = reconstructed_triangles.at(g).at(t).vertex1;
1884 helios::vec3 v2 = reconstructed_triangles.at(g).at(t).vertex2;
1892 Ngroups = reconstructed_alphamasks_center.size();
1894 for (
size_t g = 0; g < Ngroups; g++) {
1902 size_t Ngroups = reconstructed_trunk_triangles.size();
1904 for (
size_t g = 0; g < Ngroups; g++) {
1906 for (
size_t t = 0; t < reconstructed_trunk_triangles.at(g).size(); t++) {
1908 helios::vec3 v0 = reconstructed_trunk_triangles.at(g).at(t).vertex0;
1909 helios::vec3 v1 = reconstructed_trunk_triangles.at(g).at(t).vertex1;
1910 helios::vec3 v2 = reconstructed_trunk_triangles.at(g).at(t).vertex2;
1912 RGBcolor color = reconstructed_trunk_triangles.at(g).at(t).color;
1921 size_t Ngroups = reconstructed_trunk_triangles.size();
1923 for (
size_t g = 0; g < Ngroups; g++) {
1925 for (
size_t t = 0; t < reconstructed_trunk_triangles.at(g).size(); t++) {
1927 helios::vec3 v0 = reconstructed_trunk_triangles.at(g).at(t).vertex0;
1928 helios::vec3 v1 = reconstructed_trunk_triangles.at(g).at(t).vertex1;
1929 helios::vec3 v2 = reconstructed_trunk_triangles.at(g).at(t).vertex2;
1942 std::vector<uint> UUIDs;
1944 std::vector<uint> UUID_leaf_template;
1945 if (subpatches.
x > 1 || subpatches.
y > 1) {
1949 size_t Ngroups = reconstructed_alphamasks_center.size();
1951 for (
size_t g = 0; g < Ngroups; g++) {
1955 uint zone = reconstructed_alphamasks_gridcell.at(g);
1957 if (reconstructed_alphamasks_size.at(g).x > 0 && reconstructed_alphamasks_size.at(g).y > 0) {
1958 std::vector<uint> UUIDs_leaf;
1959 if (subpatches.
x == 1 && subpatches.
y == 1) {
1960 UUIDs_leaf.push_back(
context->addPatch(reconstructed_alphamasks_center.at(g), reconstructed_alphamasks_size.at(g), reconstructed_alphamasks_rotation.at(g), reconstructed_alphamasks_maskfile.c_str()));
1962 UUIDs_leaf =
context->copyPrimitive(UUID_leaf_template);
1963 context->scalePrimitive(UUIDs_leaf,
make_vec3(reconstructed_alphamasks_size.at(g).x, reconstructed_alphamasks_size.at(g).y, 1));
1964 context->rotatePrimitive(UUIDs_leaf, -reconstructed_alphamasks_rotation.at(g).elevation,
"x");
1965 context->rotatePrimitive(UUIDs_leaf, -reconstructed_alphamasks_rotation.at(g).azimuth,
"z");
1966 context->translatePrimitive(UUIDs_leaf, reconstructed_alphamasks_center.at(g));
1968 context->setPrimitiveData(UUIDs_leaf,
"gridCell", zone);
1969 uint flag = reconstructed_alphamasks_direct_flag.at(g);
1970 context->setPrimitiveData(UUIDs_leaf,
"directFlag", flag);
1971 UUIDs.insert(UUIDs.end(), UUIDs_leaf.begin(), UUIDs_leaf.end());
1975 context->deletePrimitive(UUID_leaf_template);
1982 std::vector<uint> UUIDs;
1984 size_t Ngroups = reconstructed_triangles.size();
1986 for (
size_t g = 0; g < Ngroups; g++) {
1988 int leafGroup = round(
context->randu() * (Ngroups - 1));
1990 for (
size_t t = 0; t < reconstructed_triangles.at(g).size(); t++) {
1992 helios::vec3 v0 = reconstructed_triangles.at(g).at(t).vertex0;
1993 helios::vec3 v1 = reconstructed_triangles.at(g).at(t).vertex1;
1994 helios::vec3 v2 = reconstructed_triangles.at(g).at(t).vertex2;
1996 RGBcolor color = reconstructed_triangles.at(g).at(t).color;
1998 UUIDs.push_back(
context->addTriangle(v0, v1, v2, color));
2000 uint zone = reconstructed_triangles.at(g).at(t).gridcell;
2001 context->setPrimitiveData(UUIDs.back(),
"gridCell", zone);
2003 context->setPrimitiveData(UUIDs.back(),
"leafGroup", leafGroup);
2012 std::vector<uint> UUIDs;
2014 size_t Ngroups = reconstructed_trunk_triangles.size();
2016 for (
size_t g = 0; g < Ngroups; g++) {
2018 for (
size_t t = 0; t < reconstructed_trunk_triangles.at(g).size(); t++) {
2020 helios::vec3 v0 = reconstructed_trunk_triangles.at(g).at(t).vertex0;
2021 helios::vec3 v1 = reconstructed_trunk_triangles.at(g).at(t).vertex1;
2022 helios::vec3 v2 = reconstructed_trunk_triangles.at(g).at(t).vertex2;
2024 RGBcolor color = reconstructed_trunk_triangles.at(g).at(t).color;
2026 UUIDs.push_back(
context->addTriangle(v0, v1, v2, color));
2035 if (printmessages && hits.size() == 0) {
2036 std::cout <<
"WARNING (getHitBoundingBox): There are no hit points in the point cloud, cannot determine bounding box...skipping." << std::endl;
2043 for (std::size_t i = 0; i < hits.size(); i++) {
2047 if (xyz.
x < boxmin.
x) {
2050 if (xyz.
x > boxmax.
x) {
2053 if (xyz.
y < boxmin.
y) {
2056 if (xyz.
y > boxmax.
y) {
2059 if (xyz.
z < boxmin.
z) {
2062 if (xyz.
z > boxmax.
z) {
2071 std::cout <<
"WARNING (getGridBoundingBox): There are no grid cells in the point cloud, cannot determine bounding box...skipping." << std::endl;
2078 std::size_t count = 0;
2086 vec3 xyz_min = center - 0.5f * size;
2088 vec3 xyz_max = center + 0.5f * size;
2091 if (xyz_min.
x < boxmin.
x) {
2092 boxmin.
x = xyz_min.
x;
2094 if (xyz_max.
x > boxmax.
x) {
2095 boxmax.
x = xyz_max.
x;
2097 if (xyz_min.
y < boxmin.
y) {
2098 boxmin.
y = xyz_min.
y;
2100 if (xyz_max.
y > boxmax.
y) {
2101 boxmax.
y = xyz_max.
y;
2103 if (xyz_min.
z < boxmin.
z) {
2104 boxmin.
z = xyz_min.
z;
2106 if (xyz_max.
z > boxmax.
z) {
2107 boxmax.
z = xyz_max.
z;
2114 std::size_t delete_count = 0;
2127 if (printmessages) {
2128 std::cout <<
"Removed " << delete_count <<
" hit points based on distance filter." << std::endl;
2134 std::size_t delete_count = 0;
2138 if (
R < minreflectance) {
2145 if (printmessages) {
2146 std::cout <<
"Removed " << delete_count <<
" hit points based on reflectance filter." << std::endl;
2152 std::size_t delete_count = 0;
2156 if (strcmp(comparator,
"<") == 0) {
2157 if (
R < threshold) {
2161 }
else if (strcmp(comparator,
">") == 0) {
2162 if (
R > threshold) {
2166 }
else if (strcmp(comparator,
"=") == 0) {
2167 if (
R == threshold) {
2175 if (printmessages) {
2176 std::cout <<
"Removed " << delete_count <<
" hit points based on scalar filter." << std::endl;
2182 xyzFilter(xmin, xmax, ymin, ymax, zmin, zmax,
true);
2187 if (xmin > xmax || ymin > ymax || zmin > zmax) {
2188 std::cout <<
"WARNING: at least one minimum value provided is greater than one maximum value. " << std::endl;
2191 std::size_t delete_count = 0;
2193 if (deleteOutside) {
2197 if (xyz.
x < xmin || xyz.
x > xmax || xyz.
y < ymin || xyz.
y > ymax || xyz.
z < zmin || xyz.
z > zmax) {
2206 if (xyz.
x >= xmin && xyz.
x < xmax && xyz.
y > ymin && xyz.
y < ymax && xyz.
z > zmin && xyz.
z < zmax) {
2214 if (printmessages) {
2215 std::cout <<
"Removed " << delete_count <<
" hit points based on provided bounding box." << std::endl;
2227bool sortcol0(
const std::vector<double> &v0,
const std::vector<double> &v1) {
2228 return v0.at(0) < v1.at(0);
2231bool sortcol1(
const std::vector<double> &v0,
const std::vector<double> &v1) {
2232 return v0.at(1) < v1.at(1);
2238 double median_double(std::vector<double> &v) {
2239 const size_t n = v.size();
2240 const size_t mid = n / 2;
2241 std::nth_element(v.begin(), v.begin() + mid, v.end());
2242 const double hi = v.at(mid);
2247 const double lo = *std::max_element(v.begin(), v.begin() + mid);
2248 return 0.5 * (lo + hi);
2265 bool theilSenFit(
const std::vector<double> &x,
const std::vector<double> &y,
double &slope,
double &intercept) {
2266 const size_t n = x.size();
2273 const size_t pair_cap = 1000;
2274 const size_t stride = (n > pair_cap) ? (n / pair_cap) : 1;
2276 const size_t n_strided = n / stride + 1;
2277 std::vector<double> slopes;
2278 slopes.reserve(n_strided * n_strided / 2);
2279 for (
size_t i = 0; i < n; i += stride) {
2280 for (
size_t j = i + 1; j < n; j += stride) {
2281 const double dx = x.at(j) - x.at(i);
2285 slopes.push_back((y.at(j) - y.at(i)) / dx);
2289 if (slopes.empty()) {
2293 slope = median_double(slopes);
2295 std::vector<double> intercepts;
2296 intercepts.reserve(n);
2297 for (
size_t i = 0; i < n; i++) {
2298 intercepts.push_back(y.at(i) - slope * x.at(i));
2300 intercept = median_double(intercepts);
2309 if (printmessages) {
2310 std::cout <<
"Filtering point cloud by maximum " << scalar <<
" per pulse..." << std::flush;
2313 std::vector<std::vector<double>> timestamps;
2316 std::size_t delete_count = 0;
2320 helios_runtime_error(
"ERROR (LiDARcloud::maxPulseFilter): Hit point " + std::to_string(r) +
" does not have scalar data 'timestamp', which is required for max pulse filtering.");
2322 helios_runtime_error(
"ERROR (LiDARcloud::maxPulseFilter): Hit point " + std::to_string(r) +
" does not have scalar data '" + scalar +
"', which is required for max pulse filtering.");
2328 timestamps.at(r) = v;
2331 std::sort(timestamps.begin(), timestamps.end(), sortcol0);
2333 std::vector<std::vector<double>> isort;
2334 std::vector<int> to_delete;
2335 double time_old = timestamps.at(0).at(0);
2336 for (std::size_t r = 0; r < timestamps.size(); r++) {
2338 if (timestamps.at(r).at(0) != time_old) {
2340 if (isort.size() > 1) {
2342 std::sort(isort.begin(), isort.end(), sortcol1);
2344 for (
int i = 0; i < isort.size() - 1; i++) {
2345 to_delete.push_back(
int(isort.at(i).at(2)));
2350 time_old = timestamps.at(r).at(0);
2353 isort.push_back(timestamps.at(r));
2356 std::sort(to_delete.begin(), to_delete.end());
2358 for (
int i = to_delete.size() - 1; i >= 0; i--) {
2362 if (printmessages) {
2363 std::cout <<
"done." << std::endl;
2369 if (printmessages) {
2370 std::cout <<
"Filtering point cloud by minimum " << scalar <<
" per pulse..." << std::flush;
2373 std::vector<std::vector<double>> timestamps;
2376 std::size_t delete_count = 0;
2380 helios_runtime_error(
"ERROR (LiDARcloud::minPulseFilter): Hit point " + std::to_string(r) +
" does not have scalar data 'timestamp', which is required for min pulse filtering.");
2382 helios_runtime_error(
"ERROR (LiDARcloud::minPulseFilter): Hit point " + std::to_string(r) +
" does not have scalar data '" + scalar +
"', which is required for min pulse filtering.");
2390 timestamps.at(r) = v;
2393 std::sort(timestamps.begin(), timestamps.end(), sortcol0);
2395 std::vector<std::vector<double>> isort;
2396 std::vector<int> to_delete;
2397 double time_old = timestamps.at(0).at(0);
2398 for (std::size_t r = 0; r < timestamps.size(); r++) {
2400 if (timestamps.at(r).at(0) != time_old) {
2402 if (isort.size() > 1) {
2404 std::sort(isort.begin(), isort.end(), sortcol1);
2406 for (
int i = 1; i < isort.size(); i++) {
2407 to_delete.push_back(
int(isort.at(i).at(2)));
2412 time_old = timestamps.at(r).at(0);
2415 isort.push_back(timestamps.at(r));
2418 std::sort(to_delete.begin(), to_delete.end());
2420 for (
int i = to_delete.size() - 1; i >= 0; i--) {
2424 if (printmessages) {
2425 std::cout <<
"done." << std::endl;
2431 if (printmessages) {
2432 std::cout <<
"Filtering point cloud to only first hits per pulse..." << std::flush;
2435 std::vector<float> target_index;
2439 for (std::size_t r = 0; r < target_index.size(); r++) {
2442 std::cerr <<
"failed\nERROR (LiDARcloud::firstHitFilter): Hit point " << r
2443 <<
" does not have scalar data "
2445 ". No filtering will be performed."
2450 target_index.at(r) =
getHitData(r,
"target_index");
2452 if (target_index.at(r) == 0) {
2457 for (
int r = target_index.size() - 1; r >= 0; r--) {
2459 if (target_index.at(r) != min_tindex) {
2464 if (printmessages) {
2465 std::cout <<
"done." << std::endl;
2471 if (printmessages) {
2472 std::cout <<
"Filtering point cloud to only last hits per pulse..." << std::flush;
2475 std::vector<float> target_index;
2479 for (std::size_t r = 0; r < target_index.size(); r++) {
2482 std::cout <<
"failed\n";
2483 std::cerr <<
"ERROR (LiDARcloud::lastHitFilter): Hit point " << r
2484 <<
" does not have scalar data "
2486 ". No filtering will be performed."
2490 std::cout <<
"failed\n";
2491 std::cerr <<
"ERROR (LiDARcloud::lastHitFilter): Hit point " << r
2492 <<
" does not have scalar data "
2494 ". No filtering will be performed."
2499 target_index.at(r) =
getHitData(r,
"target_index");
2501 if (target_index.at(r) == 0) {
2506 for (
int r = target_index.size() - 1; r >= 0; r--) {
2508 float target_count =
getHitData(r,
"target_count");
2510 if (target_index.at(r) == target_count - 1 + min_tindex) {
2515 if (printmessages) {
2516 std::cout <<
"done." << std::endl;
2520std::vector<helios::vec3> LiDARcloud::gapfillMisses_rowcolumn(
uint scanID,
const bool add_flags) {
2522 if (printmessages) {
2523 std::cout <<
"Gap filling complete misses in scan " << scanID <<
" using row/column indices..." << std::flush;
2530 if (scans.at(scanID).isMoving) {
2531 helios_runtime_error(
"ERROR (LiDARcloud::gapfillMisses): the row/column gap-filling path does not support moving-platform scans (see addScanMoving). A moving scan should be gap-filled via its per-pulse timestamps; ensure the scan "
2532 "data carries 'timestamp' (and not 'row'/'column') so the timestamp-based path is used.");
2537 const int Ntheta = (int) scans.at(scanID).Ntheta;
2538 const int Nphi = (int) scans.at(scanID).Nphi;
2540 std::vector<helios::vec3> xyz_filled;
2546 std::vector<std::vector<double>> row_cols(Ntheta);
2547 std::vector<std::vector<double>> row_zeniths(Ntheta);
2548 std::vector<std::vector<double>> row_azimuths(Ntheta);
2549 std::set<std::pair<int, int>> occupied;
2568 const int row = (int) std::lround(
getHitData(r,
"row"));
2569 const int col = (int) std::lround(
getHitData(r,
"column"));
2570 if (row < 0 || row >= Ntheta || col < 0 || col >= Nphi) {
2575 row_cols.at(row).push_back((
double) col);
2576 row_zeniths.at(row).push_back(raydir.
zenith);
2577 row_azimuths.at(row).push_back(raydir.
azimuth);
2578 occupied.insert(std::make_pair(row, col));
2585 const int min_returns_for_fit = 4;
2586 std::vector<double> zenith_lut(Ntheta, 0.0);
2587 std::vector<double> az_intercept_lut(Ntheta, 0.0);
2588 std::vector<double> az_slope_lut(Ntheta, 0.0);
2589 std::vector<bool> row_fitted(Ntheta,
false);
2591 for (
int row = 0; row < Ntheta; row++) {
2592 if ((
int) row_cols.at(row).size() < min_returns_for_fit) {
2597 std::vector<double> zeniths_copy = row_zeniths.at(row);
2598 const double zen = median_double(zeniths_copy);
2601 std::vector<double> az_center_copy = row_azimuths.at(row);
2602 const double az_center = median_double(az_center_copy);
2603 std::vector<double> az_unwrapped = row_azimuths.at(row);
2604 for (
double &a: az_unwrapped) {
2605 while (a - az_center > M_PI) {
2608 while (a - az_center < -M_PI) {
2613 double slope = 0.0, intercept = 0.0;
2614 if (!theilSenFit(row_cols.at(row), az_unwrapped, slope, intercept)) {
2618 zenith_lut.at(row) = zen;
2619 az_slope_lut.at(row) = slope;
2620 az_intercept_lut.at(row) = intercept;
2621 row_fitted.at(row) =
true;
2629 std::vector<double> fitted_row_idx, fitted_zenith, fitted_intercept, fitted_slope;
2630 for (
int row = 0; row < Ntheta; row++) {
2631 if (row_fitted.at(row)) {
2632 fitted_row_idx.push_back((
double) row);
2633 fitted_zenith.push_back(zenith_lut.at(row));
2634 fitted_intercept.push_back(az_intercept_lut.at(row));
2635 fitted_slope.push_back(az_slope_lut.at(row));
2639 if ((
int) fitted_row_idx.size() < 2) {
2640 helios_runtime_error(
"ERROR (LiDARcloud::gapfillMisses): scan " + std::to_string(scanID) +
" has too few populated scan rows (" + std::to_string(fitted_row_idx.size()) +
2641 ") to robustly reconstruct the row/column scan-grid model. At least 2 rows with >= " + std::to_string(min_returns_for_fit) +
" returns are required.");
2644 double zen_slope = 0.0, zen_intercept = 0.0;
2645 double int_slope = 0.0, int_intercept = 0.0;
2646 const bool zen_ok = theilSenFit(fitted_row_idx, fitted_zenith, zen_slope, zen_intercept);
2647 const bool int_ok = theilSenFit(fitted_row_idx, fitted_intercept, int_slope, int_intercept);
2648 const double median_slope = median_double(fitted_slope);
2650 for (
int row = 0; row < Ntheta; row++) {
2651 if (row_fitted.at(row)) {
2654 zenith_lut.at(row) = zen_ok ? (zen_intercept + zen_slope * (double) row) : fitted_zenith.front();
2655 az_intercept_lut.at(row) = int_ok ? (int_intercept + int_slope * (double) row) : fitted_intercept.front();
2656 az_slope_lut.at(row) = median_slope;
2660 uint npoints_interior = 0;
2661 uint npoints_extrapolated = 0;
2662 for (
int row = 0; row < Ntheta; row++) {
2663 for (
int col = 0; col < Nphi; col++) {
2665 if (occupied.find(std::make_pair(row, col)) != occupied.end()) {
2669 const double zenith = zenith_lut.at(row);
2670 double azimuth = az_intercept_lut.at(row) + az_slope_lut.at(row) * (double) col;
2672 azimuth = std::fmod(azimuth, 2.0 * M_PI);
2673 if (azimuth < 0.0) {
2674 azimuth += 2.0 *
M_PI;
2679 xyz_filled.push_back(xyz);
2681 std::map<std::string, double> data;
2682 data.insert(std::make_pair(
"is_miss", 1.0));
2683 data.insert(std::make_pair(
"row", (
double) row));
2684 data.insert(std::make_pair(
"column", (
double) col));
2685 data.insert(std::make_pair(
"nRaysHit", 0.0));
2688 data.insert(std::make_pair(
"gapfillMisses_code", row_fitted.at(row) ? 1.0 : 4.0));
2692 if (row_fitted.at(row)) {
2695 npoints_extrapolated++;
2700 if (printmessages) {
2701 std::cout <<
"filled " << xyz_filled.size() <<
" points (" << npoints_interior <<
" interior, " << npoints_extrapolated <<
" extrapolated-row)." << std::endl;
2708 std::vector<helios::vec3> xyz_filled;
2710 std::vector<helios::vec3> filled_this_scan =
gapfillMisses(scanID,
false,
false);
2711 xyz_filled.insert(xyz_filled.end(), filled_this_scan.begin(), filled_this_scan.end());
2724 helios_runtime_error(
"ERROR (LiDARcloud::gapfillMisses): Invalid scanID " + std::to_string(scanID) +
". Only " + std::to_string(
getScanCount()) +
" scans exist.");
2732 bool has_rowcolumn =
false;
2733 bool has_timestamp =
false;
2734 size_t scan_hit_count = 0;
2741 has_rowcolumn =
true;
2744 has_timestamp =
true;
2746 if (has_rowcolumn && has_timestamp) {
2754 if (scan_hit_count == 0) {
2755 if (printmessages) {
2756 std::cout <<
"Gap filling complete misses in scan " << scanID <<
"...scan has no hits. Skipping gap fill." << std::endl;
2761 if (has_rowcolumn) {
2762 return gapfillMisses_rowcolumn(scanID, add_flags);
2763 }
else if (has_timestamp) {
2764 return gapfillMisses_timestamp(scanID, gapfill_grid_only, add_flags);
2767 " has neither 'timestamp' nor 'row'/'column' hit data; cannot reconstruct miss directions. "
2768 "Provide either per-hit timestamps or scan row/column indices.");
2773std::vector<helios::vec3> LiDARcloud::gapfillMisses_timestamp(
uint scanID,
const bool gapfill_grid_only,
const bool add_flags) {
2775 if (printmessages) {
2776 std::cout <<
"Gap filling complete misses in scan " << scanID <<
"..." << std::flush;
2782 std::vector<helios::vec3> xyz_filled;
2789 const bool gapfill_is_moving = gapfill_scan.
isMoving;
2790 auto originAtTime = [&](
double timestep) ->
helios::vec3 {
2791 if (!gapfill_is_moving) {
2796 gapfill_scan.
poseAt(timestep, pos, quat);
2797 return pos + quat_rotate(quat, gapfill_scan.
lever_arm);
2802 std::vector<std::vector<double>> hit_table;
2820 helios_runtime_error(
"ERROR (LiDARcloud::gapfillMisses): Hit " + std::to_string(r) +
" is missing required 'timestamp' data. Cannot perform gap filling.");
2823 double timestamp =
getHitData(r,
"timestamp");
2824 std::vector<double> data;
2826 data.at(0) = float(r);
2827 data.at(1) = timestamp;
2828 data.at(2) = raydir.
zenith;
2830 hit_table.push_back(data);
2835 if (hit_table.empty()) {
2836 if (printmessages) {
2837 std::cout <<
"scan has no hits. Skipping gap fill." << std::endl;
2845 std::sort(hit_table.begin(), hit_table.end(), sortcol1);
2848 for (
size_t r = 0; r < hit_table.size(); r++) {
2852 if (
getHitData(hit_table.at(r).at(0),
"target_index") == 0) {
2860 int ndup_target = 0;
2865 std::vector<std::vector<double>> hit_table_semiclean;
2866 for (
size_t r = 0; r < hit_table.size(); r++) {
2870 if (
getHitData(hit_table.at(r).at(0),
"target_index") > min_tindex) {
2876 hit_table_semiclean.push_back(hit_table.at(r));
2881 std::vector<double> dt_semiclean;
2882 dt_semiclean.resize(hit_table_semiclean.size(), 0.0);
2883 for (
size_t r = 0; r + 1 < hit_table_semiclean.size(); r++) {
2885 dt_semiclean.at(r) = hit_table_semiclean.at(r + 1).at(1) - hit_table_semiclean.at(r).at(1);
2887 hit_table_semiclean.at(r).at(0) = r;
2894 std::vector<std::vector<double>> hit_table_clean;
2895 for (
size_t r = 0; r + 1 < hit_table_semiclean.size(); r++) {
2899 if (dt_semiclean.at(r) == 0) {
2904 hit_table_clean.push_back(hit_table_semiclean.at(r));
2912 if (hit_table_clean.size() < 2) {
2913 if (printmessages) {
2914 std::cout <<
"insufficient hits to reconstruct scan grid. Skipping gap fill." << std::endl;
2921 std::vector<double> dt_clean;
2922 std::vector<float> dtheta_clean;
2923 dt_clean.resize(hit_table_clean.size(), 0.0);
2924 dtheta_clean.resize(hit_table_clean.size(), 0.f);
2926 double dt_clean_min = 1e6;
2927 for (
size_t r = 0; r + 1 < hit_table_clean.size(); r++) {
2929 dt_clean.at(r) = hit_table_clean.at(r + 1).at(1) - hit_table_clean.at(r).at(1);
2930 dtheta_clean.at(r) = hit_table_clean.at(r + 1).at(2) - hit_table_clean.at(r).at(2);
2932 hit_table_clean.at(r).at(0) = r;
2934 if (dt_clean.at(r) < dt_clean_min) {
2935 dt_clean_min = dt_clean.at(r);
2941 std::vector<std::vector<std::vector<double>>> hit_table2D;
2944 hit_table2D.resize(1);
2945 for (
size_t r = 0; r + 1 < hit_table_clean.size(); r++) {
2947 hit_table2D.at(column).push_back(hit_table_clean.at(r));
2951 if (dtheta_clean.at(r) < -0.1745329f) {
2953 hit_table2D.resize(column + 1);
2964 float dtheta_avg = 0;
2967 for (
int j = 0; j < hit_table2D.size(); j++) {
2968 for (
int i = 0; i < hit_table2D.at(j).size(); i++) {
2969 int r = int(hit_table2D.at(j).at(i).at(0));
2970 if (dt_clean.at(r) >= dt_clean_min && dt_clean.at(r) < 1.5 * dt_clean_min) {
2971 dt_avg += dt_clean.at(r);
2975 dtheta_avg += dtheta_clean.at(r);
2981 if (dt_sum == 0 || dtheta_sum == 0) {
2982 if (printmessages) {
2983 std::cout <<
"insufficient valid hit pairs. Skipping gap fill." << std::endl;
2988 dt_avg = dt_avg / float(dt_sum);
2990 dtheta_avg = dtheta_avg / float(dtheta_sum);
2996 if (!std::isfinite(dt_avg) || dt_avg <= 0.f || !std::isfinite(dt_clean_min) || dt_clean_min <= 0.0) {
2997 if (printmessages) {
2998 std::cout <<
"degenerate timestamp spacing. Skipping gap fill." << std::endl;
3007 std::set<std::pair<int, int>> filled_positions;
3013 helios::int2 rc = scans.at(scanID).direction2rc(raydir);
3014 filled_positions.insert(std::make_pair(rc.
x, rc.
y));
3019 for (
int j = 0; j < hit_table2D.size(); j++) {
3021 if (hit_table2D.at(j).size() > 0) {
3022 for (
size_t i = 0; i + 1 < hit_table2D.at(j).size(); i++) {
3024 double dt = hit_table2D.at(j).at(i + 1).at(1) - hit_table2D.at(j).at(i).at(1);
3026 if (dt > 1.5f * dt_clean_min) {
3029 int Ngap = round(dt / dt_avg) - 1;
3034 int Ngap_max = (int) scans.at(scanID).Ntheta;
3035 if (Ngap > Ngap_max) {
3040 for (
int k = 1; k <= Ngap; k++) {
3042 float timestep = hit_table2D.at(j).at(i).at(1) + dt_avg * float(k);
3045 float theta = hit_table2D.at(j).at(i).at(2) + (hit_table2D.at(j).at(i + 1).at(2) - hit_table2D.at(j).at(i).at(2)) *
float(k) / float(Ngap + 1);
3046 float phi = hit_table2D.at(j).at(i).at(3) + (hit_table2D.at(j).at(i + 1).at(3) - hit_table2D.at(j).at(i).at(3)) *
float(k) / float(Ngap + 1);
3048 if (phi > 2.f * M_PI) {
3049 phi = phi - 2.f *
M_PI;
3050 }
else if (phi < 0.f) {
3051 phi = phi + 2.f *
M_PI;
3056 helios::int2 rc = scans.at(scanID).direction2rc(dir_to_check);
3060 if (rc.
x < 0 || rc.
x >= (
int) scans.at(scanID).Ntheta || rc.
y < 0 || rc.
y >= (
int) scans.at(scanID).Nphi) {
3064 auto grid_key = std::make_pair(rc.
x, rc.
y);
3067 if (filled_positions.find(grid_key) == filled_positions.end()) {
3072 xyz_filled.push_back(xyz);
3074 std::map<std::string, double> data;
3075 data.insert(std::pair<std::string, double>(
"timestamp", timestep));
3076 data.insert(std::pair<std::string, double>(
"target_index", min_tindex));
3077 data.insert(std::pair<std::string, double>(
"nRaysHit", 0.0));
3078 data.insert(std::pair<std::string, double>(
"is_miss", 1.0));
3079 if (gapfill_is_moving) {
3081 data.insert(std::pair<std::string, double>(
"origin_x", fill_origin.
x));
3082 data.insert(std::pair<std::string, double>(
"origin_y", fill_origin.
y));
3083 data.insert(std::pair<std::string, double>(
"origin_z", fill_origin.
z));
3087 data.insert(std::pair<std::string, double>(
"gapfillMisses_code", 1.0));
3090 filled_positions.insert(grid_key);
3097 uint npointsfilled = xyz_filled.size();
3100 float grid_dtheta = (theta_range.
y - theta_range.
x) /
float(scans.at(scanID).Ntheta - 1);
3101 float grid_dphi = (scans.at(scanID).phiMax - scans.at(scanID).phiMin) /
float(scans.at(scanID).Nphi - 1);
3103 if (gapfill_grid_only ==
true) {
3106 std::vector<helios::vec3> grid_vertices;
3109 grid_vertices.push_back(boxmin);
3110 grid_vertices.push_back(boxmax);
3111 grid_vertices.push_back(helios::make_vec3(boxmin.
x, boxmin.
y, boxmax.
z));
3112 grid_vertices.push_back(helios::make_vec3(boxmax.
x, boxmax.
y, boxmin.
z));
3113 grid_vertices.push_back(helios::make_vec3(boxmin.
x, boxmax.
y, boxmin.
z));
3114 grid_vertices.push_back(helios::make_vec3(boxmin.
x, boxmax.
y, boxmax.
z));
3115 grid_vertices.push_back(helios::make_vec3(boxmax.
x, boxmin.
y, boxmin.
z));
3116 grid_vertices.push_back(helios::make_vec3(boxmax.
x, boxmin.
y, boxmax.
z));
3118 float max_theta = 0;
3119 float min_theta =
M_PI;
3121 float min_phi = 2 *
M_PI;
3122 for (
uint gg = 0; gg < grid_vertices.size(); gg++) {
3133 if (sc.
zenith < min_theta) {
3137 if (sc.
zenith > max_theta) {
3143 if (min_theta < theta_range.
x) {
3144 min_theta = theta_range.
x;
3147 if (max_theta > theta_range.
y) {
3148 max_theta = theta_range.
y;
3155 for (
int j = 0; j < hit_table2D.size(); j++) {
3157 if (hit_table2D.at(j).size() > 0) {
3160 if (hit_table2D.at(j).front().at(2) > theta_range.
x) {
3162 float dtheta = dtheta_avg;
3163 float theta = hit_table2D.at(j).at(0).at(2) - dtheta;
3165 float phi = hit_table2D.at(j).at(0).at(3);
3166 float timestep = hit_table2D.at(j).at(0).at(1) - dt_avg;
3171 while (theta > theta_range.
x) {
3175 helios::int2 rc = scans.at(scanID).direction2rc(dir_to_check);
3178 if (rc.
x >= 0 && rc.
x < (
int) scans.at(scanID).Ntheta && rc.
y >= 0 && rc.
y < (
int) scans.at(scanID).Nphi) {
3181 auto grid_key = std::make_pair(rc.
x, rc.
y);
3182 if (filled_positions.find(grid_key) == filled_positions.end()) {
3187 xyz_filled.push_back(xyz);
3189 std::map<std::string, double> data;
3190 data.insert(std::pair<std::string, double>(
"timestamp", timestep));
3191 data.insert(std::pair<std::string, double>(
"target_index", min_tindex));
3192 data.insert(std::pair<std::string, double>(
"nRaysHit", 0.0));
3193 data.insert(std::pair<std::string, double>(
"is_miss", 1.0));
3194 if (gapfill_is_moving) {
3195 data.insert(std::pair<std::string, double>(
"origin_x", fill_origin.
x));
3196 data.insert(std::pair<std::string, double>(
"origin_y", fill_origin.
y));
3197 data.insert(std::pair<std::string, double>(
"origin_z", fill_origin.
z));
3201 data.insert(std::pair<std::string, double>(
"gapfillMisses_code", 3.0));
3205 filled_positions.insert(grid_key);
3209 theta = theta - dtheta;
3210 timestep = timestep - dt_avg;
3215 if (hit_table2D.at(j).back().at(2) < theta_range.
y) {
3217 int sz = hit_table2D.at(j).size();
3219 float dtheta = dtheta_avg;
3220 float theta = hit_table2D.at(j).at(sz - 1).at(2) + dtheta;
3221 float phi = hit_table2D.at(j).at(sz - 1).at(3);
3222 float timestep = hit_table2D.at(j).at(sz - 1).at(1) + dt_avg;
3223 while (theta < theta_range.
y) {
3227 helios::int2 rc = scans.at(scanID).direction2rc(dir_to_check);
3230 if (rc.
x >= 0 && rc.
x < (
int) scans.at(scanID).Ntheta && rc.
y >= 0 && rc.
y < (
int) scans.at(scanID).Nphi) {
3233 auto grid_key = std::make_pair(rc.
x, rc.
y);
3234 if (filled_positions.find(grid_key) == filled_positions.end()) {
3239 xyz_filled.push_back(xyz);
3241 std::map<std::string, double> data;
3242 data.insert(std::pair<std::string, double>(
"timestamp", timestep));
3243 data.insert(std::pair<std::string, double>(
"target_index", min_tindex));
3244 data.insert(std::pair<std::string, double>(
"nRaysHit", 0.0));
3245 data.insert(std::pair<std::string, double>(
"is_miss", 1.0));
3246 if (gapfill_is_moving) {
3247 data.insert(std::pair<std::string, double>(
"origin_x", fill_origin.
x));
3248 data.insert(std::pair<std::string, double>(
"origin_y", fill_origin.
y));
3249 data.insert(std::pair<std::string, double>(
"origin_z", fill_origin.
z));
3253 data.insert(std::pair<std::string, double>(
"gapfillMisses_code", 2.0));
3257 filled_positions.insert(grid_key);
3261 theta = theta + dtheta;
3262 timestep = timestep + dt_avg;
3268 uint npointsextrapolated = xyz_filled.size() - npointsfilled;
3270 if (printmessages) {
3271 std::cout <<
"filled " << xyz_filled.size() <<
" points (" << npointsfilled <<
" interior, " << npointsextrapolated <<
" edge)." << std::endl;
3272 std::cout <<
" Processed " << hit_table2D.size() <<
" scan columns" << std::endl;
3283 if (anyScanMoving()) {
3284 helios_runtime_error(
"ERROR (LiDARcloud::triangulateHitPoints): triangulation is not supported for moving-platform scans (see addScanMoving), which have no fixed theta-phi scan grid to triangulate. For leaf-area inversion of a "
3285 "moving scan, call the calculateLeafArea overload that takes a G(theta) argument (it does not require triangulation).");
3289 cout <<
"WARNING (triangulateHitPoints): No scans have been added to the point cloud. Skipping triangulation..." << endl;
3292 cout <<
"WARNING (triangulateHitPoints): No hit points have been added to the point cloud. Skipping triangulation..." << endl;
3296 if (!hitgridcellcomputed) {
3303 triangulation_candidate_count = 0;
3304 triangulation_dropped_lmax = 0;
3305 triangulation_dropped_aspect = 0;
3306 triangulation_dropped_degenerate = 0;
3309 bool use_adaptive_threshold = isMultiReturnData();
3310 float adaptive_sep_threshold = 0.0f;
3312 if (use_adaptive_threshold) {
3313 if (printmessages) {
3314 std::cout <<
"Multi-return data detected - calculating adaptive separation ratio threshold..." << std::endl;
3318 std::vector<float> all_separation_ratios;
3321 std::vector<int> Delaunay_inds_pass1;
3322 std::vector<Shx> pts_pass1, pts_copy_pass1;
3323 int count_pass1 = 0;
3328 if (use_adaptive_threshold) {
3339 pt.id = count_pass1;
3342 pts_pass1.push_back(pt);
3343 Delaunay_inds_pass1.push_back(r);
3348 if (pts_pass1.size() == 0)
3352 float h[2] = {0, 0};
3353 for (
int r = 0; r < pts_pass1.size(); r++) {
3354 if (pts_pass1.at(r).c < 0.5 *
M_PI)
3356 else if (pts_pass1.at(r).c > 1.5 *
M_PI)
3359 h[0] /= float(pts_pass1.size());
3360 h[1] /= float(pts_pass1.size());
3361 if (h[0] + h[1] > 0.4) {
3362 for (
int r = 0; r < pts_pass1.size(); r++) {
3363 pts_pass1.at(r).c +=
M_PI;
3364 if (pts_pass1.at(r).c > 2.f *
M_PI)
3365 pts_pass1.at(r).c -= 2.f *
M_PI;
3369 std::vector<int> dupes_pass1;
3370 de_duplicate(pts_pass1, dupes_pass1);
3372 std::vector<Triad> triads_pass1;
3376 int success = triangulate_CDT(pts_pass1, triads_pass1);
3382 for (
int t = 0; t < triads_pass1.size(); t++) {
3383 int ID0 = Delaunay_inds_pass1.at(triads_pass1.at(t).a);
3384 int ID1 = Delaunay_inds_pass1.at(triads_pass1.at(t).b);
3385 int ID2 = Delaunay_inds_pass1.at(triads_pass1.at(t).c);
3394 float L0 = (v0 - v1).magnitude();
3395 float L1 = (v0 - v2).magnitude();
3396 float L2 = (v1 - v2).magnitude();
3399 if (L0 > Lmax || L1 > Lmax || L2 > Lmax) {
3407 float ratio01 = L0 / (ang01 + 1e-6);
3408 float ratio02 = L1 / (ang02 + 1e-6);
3409 float ratio12 = L2 / (ang12 + 1e-6);
3410 float max_sep_ratio =
max(
max(ratio01, ratio02), ratio12);
3412 all_separation_ratios.push_back(max_sep_ratio);
3417 if (!all_separation_ratios.empty()) {
3418 std::sort(all_separation_ratios.begin(), all_separation_ratios.end());
3419 size_t idx_25 = all_separation_ratios.size() / 4;
3420 float percentile_25 = all_separation_ratios[idx_25];
3421 adaptive_sep_threshold = 8.5f * percentile_25;
3423 if (printmessages) {
3424 std::cout <<
" 25th percentile separation ratio: " << percentile_25 << std::endl;
3425 std::cout <<
" Adaptive threshold: " << adaptive_sep_threshold << std::endl;
3435 if (triangulationCancelled(cancel_flag)) {
3440 std::vector<int> Delaunay_inds;
3442 std::vector<Shx> pts, pts_copy;
3449 if ((r & 0x1FFFF) == 0 && triangulationCancelled(cancel_flag)) {
3456 if (use_adaptive_threshold) {
3474 Delaunay_inds.push_back(r);
3480 if (pts.size() == 0) {
3481 if (printmessages) {
3482 std::cout <<
"Scan " << s <<
" contains no triangles. Skipping this scan..." << std::endl;
3487 float h[2] = {0, 0};
3488 for (
int r = 0; r < pts.size(); r++) {
3489 if (pts.at(r).c < 0.5 *
M_PI) {
3491 }
else if (pts.at(r).c > 1.5 *
M_PI) {
3495 h[0] /= float(pts.size());
3496 h[1] /= float(pts.size());
3497 if (h[0] + h[1] > 0.4) {
3498 if (printmessages) {
3499 std::cout <<
"Shifting scan " << s << std::endl;
3501 for (
int r = 0; r < pts.size(); r++) {
3502 pts.at(r).c +=
M_PI;
3503 if (pts.at(r).c > 2.f *
M_PI) {
3504 pts.at(r).c -= 2.f *
M_PI;
3515 const float COORD_SNAP_PRECISION = 1e-6f;
3516 for (
auto &pt: pts) {
3517 pt.r = std::round(pt.r / COORD_SNAP_PRECISION) * COORD_SNAP_PRECISION;
3518 pt.c = std::round(pt.c / COORD_SNAP_PRECISION) * COORD_SNAP_PRECISION;
3521 std::vector<int> dupes;
3522 int nx = de_duplicate(pts, dupes);
3525 std::vector<Triad> triads;
3527 if (printmessages) {
3528 std::cout <<
"starting triangulation for scan " << s <<
"..." << std::endl;
3533 if (triangulationCancelled(cancel_flag)) {
3541 int success = triangulate_CDT(pts, triads);
3545 if (triangulationCancelled(cancel_flag)) {
3551 if (printmessages) {
3552 std::cout <<
"FAILED: could not triangulate scan " << s <<
". Skipping this scan." << std::endl;
3555 }
else if (printmessages) {
3556 std::cout <<
"finished triangulation" << std::endl;
3559 triangulation_candidate_count += triads.size();
3561 for (
int t = 0; t < triads.size(); t++) {
3564 if ((t & 0x1FFFF) == 0 && triangulationCancelled(cancel_flag)) {
3569 int ID0 = Delaunay_inds.at(triads.at(t).a);
3570 int ID1 = Delaunay_inds.at(triads.at(t).b);
3571 int ID2 = Delaunay_inds.at(triads.at(t).c);
3583 v = vertex0 - vertex1;
3585 v = vertex0 - vertex2;
3587 v = vertex1 - vertex2;
3590 float aspect_ratio =
max(
max(L0, L1), L2) /
min(
min(L0, L1), L2);
3596 bool dropped_lmax = (L0 > Lmax || L1 > Lmax || L2 > Lmax);
3597 bool dropped_aspect =
false;
3599 if (use_adaptive_threshold) {
3605 float ratio01 = L0 / (ang01 + 1e-6);
3606 float ratio02 = L1 / (ang02 + 1e-6);
3607 float ratio12 = L2 / (ang12 + 1e-6);
3608 float max_sep_ratio =
max(
max(ratio01, ratio02), ratio12);
3610 dropped_aspect = (max_sep_ratio > adaptive_sep_threshold) || (aspect_ratio > max_aspect_ratio);
3613 dropped_aspect = (aspect_ratio > max_aspect_ratio);
3617 triangulation_dropped_lmax++;
3620 if (dropped_aspect) {
3621 triangulation_dropped_aspect++;
3627 if (printmessages && gridcell == -2) {
3628 cout <<
"WARNING (triangulateHitPoints): You typically want to define the hit grid cell for all hit points before performing triangulation." << endl;
3632 color.
r = (hits.at(ID0).color.r + hits.at(ID1).color.r + hits.at(ID2).color.r) / 3.f;
3633 color.
g = (hits.at(ID0).color.g + hits.at(ID1).color.g + hits.at(ID2).color.g) / 3.f;
3634 color.
b = (hits.at(ID0).color.b + hits.at(ID1).color.b + hits.at(ID2).color.b) / 3.f;
3636 Triangulation tri(s, vertex0, vertex1, vertex2, ID0, ID1, ID2, color, gridcell);
3638 if (tri.area != tri.area) {
3639 triangulation_dropped_degenerate++;
3643 triangles.push_back(tri);
3649 triangulationcomputed =
true;
3651 if (printmessages) {
3653 cout <<
"\rTriangulating...formed " << Ntriangles <<
" total triangles." << endl;
3657int LiDARcloud::getContainingGridCell(
const helios::vec3 &p)
const {
3666 for (
uint c = 0; c < Ncells; c++) {
3675 if (fabs(rotation) > 1e-6f) {
3677 q =
rotatePointAboutLine(p - anchor, helios::make_vec3(0, 0, 0), helios::make_vec3(0, 0, 1), -rotation) + anchor;
3680 if (q.
x >= lo.
x && q.
x <= hi.
x && q.
y >= lo.
y && q.
y <= hi.
y && q.
z >= lo.
z && q.
z <= hi.
z) {
3681 return static_cast<int>(c);
3689 if (triangle_vertices.size() % 3 != 0) {
3690 helios_runtime_error(
"ERROR (LiDARcloud::setExternalTriangulation): triangle_vertices size (" + std::to_string(triangle_vertices.size()) +
") must be a multiple of 3 (three vertices per triangle).");
3693 const size_t Ntri = triangle_vertices.size() / 3;
3695 if (scanIDs.size() != Ntri) {
3696 helios_runtime_error(
"ERROR (LiDARcloud::setExternalTriangulation): scanIDs size (" + std::to_string(scanIDs.size()) +
") must equal the triangle count (" + std::to_string(Ntri) +
3697 "). Each triangle requires a source scan for the G(theta) ray direction.");
3701 helios_runtime_error(
"ERROR (LiDARcloud::setExternalTriangulation): a grid must be defined (see addGrid()) before supplying an external triangulation, so each triangle can be assigned to a grid cell.");
3705 for (
size_t t = 0; t < Ntri; t++) {
3706 if (scanIDs.at(t) < 0 ||
static_cast<uint>(scanIDs.at(t)) >= Nscans) {
3707 helios_runtime_error(
"ERROR (LiDARcloud::setExternalTriangulation): triangle " + std::to_string(t) +
" has scanID " + std::to_string(scanIDs.at(t)) +
", which is not a valid scan index in [0, " + std::to_string(Nscans) +
3708 "). Per-scan provenance is required; a merged mesh with no scan association is not a valid input.");
3714 triangulation_candidate_count = Ntri;
3715 triangulation_dropped_lmax = 0;
3716 triangulation_dropped_aspect = 0;
3717 triangulation_dropped_degenerate = 0;
3719 for (
size_t t = 0; t < Ntri; t++) {
3721 const helios::vec3 &v0 = triangle_vertices.at(3 * t + 0);
3722 const helios::vec3 &v1 = triangle_vertices.at(3 * t + 1);
3723 const helios::vec3 &v2 = triangle_vertices.at(3 * t + 2);
3727 int gridcell = getContainingGridCell(centroid);
3730 Triangulation tri(scanIDs.at(t), v0, v1, v2, -1, -1, -1, helios::RGB::green, gridcell);
3734 if (tri.area != tri.area || tri.area <= 0.f) {
3735 triangulation_dropped_degenerate++;
3739 triangles.push_back(tri);
3742 triangulationcomputed =
true;
3744 if (printmessages) {
3745 std::cout <<
"Set external triangulation: " << triangles.size() <<
" triangles (" << triangulation_dropped_degenerate <<
" degenerate dropped)." << std::endl;
3752 if (anyScanMoving()) {
3753 helios_runtime_error(
"ERROR (LiDARcloud::triangulateHitPoints): triangulation is not supported for moving-platform scans (see addScanMoving), which have no fixed theta-phi scan grid to triangulate. For leaf-area inversion of a "
3754 "moving scan, call the calculateLeafArea overload that takes a G(theta) argument (it does not require triangulation).");
3758 cout <<
"WARNING (triangulateHitPoints): No scans have been added to the point cloud. Skipping triangulation..." << endl;
3761 cout <<
"WARNING (triangulateHitPoints): No hit points have been added to the point cloud. Skipping triangulation..." << endl;
3765 if (!hitgridcellcomputed) {
3772 triangulation_candidate_count = 0;
3773 triangulation_dropped_lmax = 0;
3774 triangulation_dropped_aspect = 0;
3775 triangulation_dropped_degenerate = 0;
3778 bool use_adaptive_threshold = isMultiReturnData();
3779 float adaptive_sep_threshold = 0.0f;
3781 if (use_adaptive_threshold) {
3782 if (printmessages) {
3783 std::cout <<
"Multi-return data detected - calculating adaptive separation ratio threshold..." << std::endl;
3787 std::vector<float> all_separation_ratios;
3790 std::vector<int> Delaunay_inds_pass1;
3791 std::vector<Shx> pts_pass1, pts_copy_pass1;
3792 int count_pass1 = 0;
3801 pt.id = count_pass1;
3804 pts_pass1.push_back(pt);
3805 Delaunay_inds_pass1.push_back(r);
3810 if (pts_pass1.size() == 0)
3814 float h[2] = {0, 0};
3815 for (
int r = 0; r < pts_pass1.size(); r++) {
3816 if (pts_pass1.at(r).c < 0.5 *
M_PI)
3818 else if (pts_pass1.at(r).c > 1.5 *
M_PI)
3821 h[0] /= float(pts_pass1.size());
3822 h[1] /= float(pts_pass1.size());
3823 if (h[0] + h[1] > 0.4) {
3824 for (
int r = 0; r < pts_pass1.size(); r++) {
3825 pts_pass1.at(r).c +=
M_PI;
3826 if (pts_pass1.at(r).c > 2.f *
M_PI)
3827 pts_pass1.at(r).c -= 2.f *
M_PI;
3831 std::vector<int> dupes_pass1;
3832 de_duplicate(pts_pass1, dupes_pass1);
3834 std::vector<Triad> triads_pass1;
3838 int success = triangulate_CDT(pts_pass1, triads_pass1);
3844 for (
int t = 0; t < triads_pass1.size(); t++) {
3845 int ID0 = Delaunay_inds_pass1.at(triads_pass1.at(t).a);
3846 int ID1 = Delaunay_inds_pass1.at(triads_pass1.at(t).b);
3847 int ID2 = Delaunay_inds_pass1.at(triads_pass1.at(t).c);
3856 float L0 = (v0 - v1).magnitude();
3857 float L1 = (v0 - v2).magnitude();
3858 float L2 = (v1 - v2).magnitude();
3861 if (L0 > Lmax || L1 > Lmax || L2 > Lmax) {
3869 float ratio01 = L0 / (ang01 + 1e-6);
3870 float ratio02 = L1 / (ang02 + 1e-6);
3871 float ratio12 = L2 / (ang12 + 1e-6);
3872 float max_sep_ratio =
max(
max(ratio01, ratio02), ratio12);
3874 all_separation_ratios.push_back(max_sep_ratio);
3879 if (!all_separation_ratios.empty()) {
3880 std::sort(all_separation_ratios.begin(), all_separation_ratios.end());
3881 size_t idx_25 = all_separation_ratios.size() / 4;
3882 float percentile_25 = all_separation_ratios[idx_25];
3883 adaptive_sep_threshold = 8.5f * percentile_25;
3885 if (printmessages) {
3886 std::cout <<
" 25th percentile separation ratio: " << percentile_25 << std::endl;
3887 std::cout <<
" Adaptive threshold: " << adaptive_sep_threshold << std::endl;
3897 if (triangulationCancelled(cancel_flag)) {
3902 std::vector<int> Delaunay_inds;
3904 std::vector<Shx> pts, pts_copy;
3906 std::size_t delete_count = 0;
3913 if ((r & 0x1FFFF) == 0 && triangulationCancelled(cancel_flag)) {
3922 if (strcmp(comparator,
"<") == 0) {
3923 if (
R < threshold) {
3927 }
else if (strcmp(comparator,
">") == 0) {
3928 if (
R > threshold) {
3932 }
else if (strcmp(comparator,
"=") == 0) {
3933 if (
R == threshold) {
3953 Delaunay_inds.push_back(r);
3959 if (printmessages) {
3960 std::cout <<
"Scan " << s <<
" triangulation: " << count <<
" points used, " << delete_count <<
" points filtered out";
3961 if (strlen(scalar_field) > 0) {
3962 std::cout <<
" (filter: " << scalar_field <<
" " << comparator <<
" " << threshold <<
")";
3964 std::cout << std::endl;
3967 if (pts.size() == 0) {
3968 if (printmessages) {
3969 std::cout <<
"Scan " << s <<
" contains no triangles. Skipping this scan..." << std::endl;
3974 float h[2] = {0, 0};
3975 for (
int r = 0; r < pts.size(); r++) {
3976 if (pts.at(r).c < 0.5 *
M_PI) {
3978 }
else if (pts.at(r).c > 1.5 *
M_PI) {
3982 h[0] /= float(pts.size());
3983 h[1] /= float(pts.size());
3984 if (h[0] + h[1] > 0.4) {
3985 if (printmessages) {
3986 std::cout <<
"Shifting scan " << s << std::endl;
3988 for (
int r = 0; r < pts.size(); r++) {
3989 pts.at(r).c +=
M_PI;
3990 if (pts.at(r).c > 2.f *
M_PI) {
3991 pts.at(r).c -= 2.f *
M_PI;
4002 const float COORD_SNAP_PRECISION = 1e-6f;
4003 for (
auto &pt: pts) {
4004 pt.r = std::round(pt.r / COORD_SNAP_PRECISION) * COORD_SNAP_PRECISION;
4005 pt.c = std::round(pt.c / COORD_SNAP_PRECISION) * COORD_SNAP_PRECISION;
4008 std::vector<int> dupes;
4009 int nx = de_duplicate(pts, dupes);
4012 std::vector<Triad> triads;
4014 if (printmessages) {
4015 std::cout <<
"starting triangulation for scan " << s <<
"..." << std::endl;
4020 if (triangulationCancelled(cancel_flag)) {
4028 int success = triangulate_CDT(pts, triads);
4032 if (triangulationCancelled(cancel_flag)) {
4038 if (printmessages) {
4039 std::cout <<
"FAILED: could not triangulate scan " << s <<
". Skipping this scan." << std::endl;
4042 }
else if (printmessages) {
4043 std::cout <<
"finished triangulation" << std::endl;
4046 triangulation_candidate_count += triads.size();
4048 for (
int t = 0; t < triads.size(); t++) {
4051 if ((t & 0x1FFFF) == 0 && triangulationCancelled(cancel_flag)) {
4056 int ID0 = Delaunay_inds.at(triads.at(t).a);
4057 int ID1 = Delaunay_inds.at(triads.at(t).b);
4058 int ID2 = Delaunay_inds.at(triads.at(t).c);
4070 v = vertex0 - vertex1;
4072 v = vertex0 - vertex2;
4074 v = vertex1 - vertex2;
4077 float aspect_ratio =
max(
max(L0, L1), L2) /
min(
min(L0, L1), L2);
4083 bool dropped_lmax = (L0 > Lmax || L1 > Lmax || L2 > Lmax);
4084 bool dropped_aspect =
false;
4086 if (use_adaptive_threshold) {
4092 float ratio01 = L0 / (ang01 + 1e-6);
4093 float ratio02 = L1 / (ang02 + 1e-6);
4094 float ratio12 = L2 / (ang12 + 1e-6);
4095 float max_sep_ratio =
max(
max(ratio01, ratio02), ratio12);
4097 dropped_aspect = (max_sep_ratio > adaptive_sep_threshold) || (aspect_ratio > max_aspect_ratio);
4100 dropped_aspect = (aspect_ratio > max_aspect_ratio);
4104 triangulation_dropped_lmax++;
4107 if (dropped_aspect) {
4108 triangulation_dropped_aspect++;
4114 if (printmessages && gridcell == -2) {
4115 cout <<
"WARNING (triangulateHitPoints): You typically want to define the hit grid cell for all hit points before performing triangulation." << endl;
4119 color.
r = (hits.at(ID0).color.r + hits.at(ID1).color.r + hits.at(ID2).color.r) / 3.f;
4120 color.
g = (hits.at(ID0).color.g + hits.at(ID1).color.g + hits.at(ID2).color.g) / 3.f;
4121 color.
b = (hits.at(ID0).color.b + hits.at(ID1).color.b + hits.at(ID2).color.b) / 3.f;
4123 Triangulation tri(s, vertex0, vertex1, vertex2, ID0, ID1, ID2, color, gridcell);
4125 if (tri.area != tri.area) {
4126 triangulation_dropped_degenerate++;
4130 triangles.push_back(tri);
4136 triangulationcomputed =
true;
4138 if (printmessages) {
4140 cout <<
"\rTriangulating...formed " << Ntriangles <<
" total triangles." << endl;
4147 if (scans.size() == 0) {
4148 if (printmessages) {
4149 std::cout <<
"WARNING (addTrianglesToContext): There are no scans in the point cloud, and thus there are no triangles to add...skipping." << std::endl;
4158 context->addTriangle(tri.vertex0, tri.vertex1, tri.vertex2, tri.color);
4163 return grid_cells.size();
4167 addGridCell(center, center, size, size, rotation,
make_int3(1, 1, 1),
make_int3(1, 1, 1));
4172 GridCell newcell(center, global_anchor, size, global_size, rotation, global_ijk, global_count);
4174 grid_cells.push_back(newcell);
4180 helios_runtime_error(
"ERROR (LiDARcloud::getCellCenter): grid cell index out of range. Requested center of cell #" + std::to_string(index) +
" but there are only " + std::to_string(
getGridCellCount()) +
" cells in the grid.");
4183 return grid_cells.at(index).center;
4189 helios_runtime_error(
"ERROR (LiDARcloud::getCellGlobalAnchor): grid cell index out of range. Requested anchor of cell #" + std::to_string(index) +
" but there are only " + std::to_string(
getGridCellCount()) +
" cells in the grid.");
4192 return grid_cells.at(index).global_anchor;
4198 helios_runtime_error(
"ERROR (LiDARcloud::getCellCenter): grid cell index out of range. Requested size of cell #" + std::to_string(index) +
" but there are only " + std::to_string(
getGridCellCount()) +
" cells in the grid.");
4201 return grid_cells.at(index).size;
4207 helios_runtime_error(
"ERROR (LiDARcloud::getCellRotation): grid cell index out of range. Requested rotation of cell #" + std::to_string(index) +
" but there are only " + std::to_string(
getGridCellCount()) +
" cells in the grid.");
4210 return grid_cells.at(index).azimuthal_rotation;
4215 size_t Nprims =
context->getPrimitiveCount();
4221 std::vector<float> Gtheta;
4222 Gtheta.resize(Ncells);
4224 std::vector<float> area_sum;
4225 area_sum.resize(Ncells, 0.f);
4226 std::vector<uint> cell_tri_count;
4227 cell_tri_count.resize(Ncells, 0);
4229 std::vector<uint> UUIDs =
context->getAllUUIDs();
4230 for (
int p = 0; p < UUIDs.size(); p++) {
4232 uint UUID = UUIDs.at(p);
4234 if (
context->doesPrimitiveDataExist(UUID,
"gridCell")) {
4237 context->getPrimitiveData(UUID,
"gridCell", gridCell);
4239 std::vector<vec3> vertices =
context->getPrimitiveVertices(UUID);
4240 float area =
context->getPrimitiveArea(UUID);
4243 for (
int s = 0; s < Nscans; s++) {
4245 vec3 raydir = vertices.front() - origin;
4250 Gtheta.at(gridCell) += fabs(normal * raydir) * area;
4252 area_sum.at(gridCell) += area;
4253 cell_tri_count.at(gridCell) += 1;
4259 for (
uint v = 0; v < Ncells; v++) {
4260 if (cell_tri_count[v] > 0) {
4261 Gtheta[v] *= float(cell_tri_count[v]) / (area_sum[v]);
4266 std::vector<float> output_Gtheta;
4267 output_Gtheta.resize(Ncells, 0.f);
4269 for (
int v = 0; v < Ncells; v++) {
4270 output_Gtheta.at(v) = Gtheta.at(v);
4271 if (
context->doesPrimitiveDataExist(UUIDs.at(v),
"gridCell")) {
4272 context->setPrimitiveData(UUIDs.at(v),
"synthetic_Gtheta", Gtheta.at(v));
4276 return output_Gtheta;
4285 grid_cells.at(index).leaf_area = area;
4291 helios_runtime_error(
"ERROR (LiDARcloud::getCellLeafArea): grid cell index out of range. Requested leaf area of cell #" + std::to_string(index) +
" but there are only " + std::to_string(
getGridCellCount()) +
" cells in the grid.");
4294 return grid_cells.at(index).leaf_area;
4300 helios_runtime_error(
"ERROR (LiDARcloud::getCellLeafAreaDensity): grid cell index out of range. Requested leaf area density of cell #" + std::to_string(index) +
" but there are only " + std::to_string(
getGridCellCount()) +
4301 " cells in the grid.");
4305 return grid_cells.at(index).leaf_area / (gridsize.
x * gridsize.
y * gridsize.
z);
4314 grid_cells.at(index).Gtheta = Gtheta;
4320 helios_runtime_error(
"ERROR (LiDARcloud::getCellGtheta): grid cell index out of range. Requested leaf area of cell #" + std::to_string(index) +
" but there are only " + std::to_string(
getGridCellCount()) +
" cells in the grid.");
4323 return grid_cells.at(index).Gtheta;
4329 helios_runtime_error(
"ERROR (LiDARcloud::getCellBeamCount): grid cell index out of range. Requested beam count of cell #" + std::to_string(index) +
" but there are only " + std::to_string(
getGridCellCount()) +
" cells in the grid.");
4332 return grid_cells.at(index).beam_count;
4338 helios_runtime_error(
"ERROR (LiDARcloud::getCellRelativeDensityIndex): grid cell index out of range. Requested RDI of cell #" + std::to_string(index) +
" but there are only " + std::to_string(
getGridCellCount()) +
" cells in the grid.");
4341 return grid_cells.at(index).I_rdi;
4347 helios_runtime_error(
"ERROR (LiDARcloud::getCellMeanPathLength): grid cell index out of range. Requested mean path length of cell #" + std::to_string(index) +
" but there are only " + std::to_string(
getGridCellCount()) +
4348 " cells in the grid.");
4351 return grid_cells.at(index).zbar_e;
4357 helios_runtime_error(
"ERROR (LiDARcloud::getCellLADVariance): grid cell index out of range. Requested LAD variance of cell #" + std::to_string(index) +
" but there are only " + std::to_string(
getGridCellCount()) +
" cells in the grid.");
4360 return grid_cells.at(index).LAD_variance;
4366 helios_runtime_error(
"ERROR (LiDARcloud::getCellLeafAreaConfidenceInterval): grid cell index out of range. Requested cell #" + std::to_string(index) +
" but there are only " + std::to_string(
getGridCellCount()) +
" cells in the grid.");
4368 if (confidence_level <= 0.f || confidence_level >= 1.f) {
4369 helios_runtime_error(
"ERROR (LiDARcloud::getCellLeafAreaConfidenceInterval): confidence_level must be strictly between 0 and 1.");
4372 const GridCell &cell = grid_cells.at(index);
4377 const float volume = cell.size.
x * cell.size.
y * cell.size.
z;
4378 const float a = (volume > 0.f) ? cell.leaf_area / volume : 0.f;
4379 const float L = a * cell.Gtheta * cell.
zbar_e;
4385 const double z = normalQuantile(1.0 - (1.0 - (
double) confidence_level) / 2.0);
4387 const float half_width = (float) z * volume * lad_se;
4388 lower = cell.leaf_area - half_width;
4392 upper = cell.leaf_area + half_width;
4398 if (confidence_level <= 0.f || confidence_level >= 1.f) {
4399 helios_runtime_error(
"ERROR (LiDARcloud::getGroupLADConfidenceInterval): confidence_level must be strictly between 0 and 1.");
4405 double sum_lad = 0.0;
4406 double sum_variance = 0.0;
4408 for (
uint index: indices) {
4410 helios_runtime_error(
"ERROR (LiDARcloud::getGroupLADConfidenceInterval): grid cell index out of range. Requested cell #" + std::to_string(index) +
" but there are only " + std::to_string(
getGridCellCount()) +
" cells in the grid.");
4412 const GridCell &cell = grid_cells.at(index);
4416 const float volume = cell.size.
x * cell.size.
y * cell.size.
z;
4417 const float a = (volume > 0.f) ? cell.leaf_area / volume : 0.f;
4418 const float L = a * cell.Gtheta * cell.
zbar_e;
4431 mean_lad = (float) (sum_lad / (
double) n_valid);
4432 const double z = normalQuantile(1.0 - (1.0 - (
double) confidence_level) / 2.0);
4433 const float half_width = (float) (z * std::sqrt(sum_variance) / (double) n_valid);
4434 lower = mean_lad - half_width;
4438 upper = mean_lad + half_width;
4442void LiDARcloud::leafReconstructionFloodfill() {
4444 size_t group_count = 0;
4445 int current_group = 0;
4447 vector<vector<int>> nodes;
4455 if (tri.gridcell >= 0) {
4457 nodes.at(tri.ID0).push_back(t);
4458 nodes.at(tri.ID1).push_back(t);
4459 nodes.at(tri.ID2).push_back(t);
4465 std::vector<int> fill_flag;
4466 fill_flag.resize(Ntri);
4467 for (
size_t t = 0; t < Ntri; t++) {
4468 fill_flag.at(t) = -1;
4471 for (
size_t t = 0; t < Ntri; t++) {
4473 if (fill_flag.at(t) < 0) {
4475 floodfill(t, triangles, fill_flag, nodes, current_group, 0, 1e3);
4481 for (
size_t t = 0; t < Ntri; t++) {
4483 if (fill_flag.at(t) >= 0) {
4484 int fill_group = fill_flag.at(t);
4486 if (fill_group >= reconstructed_triangles.size()) {
4487 reconstructed_triangles.resize(fill_group + 1);
4490 reconstructed_triangles.at(fill_group).push_back(triangles.at(t));
4495void LiDARcloud::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) {
4499 int verts[3] = {tri.ID0, tri.ID1, tri.ID2};
4501 std::vector<int> connection_list;
4503 for (
int i = 0; i < 3; i++) {
4504 std::vector<int> connected_tris = nodes.at(verts[i]);
4505 connection_list.insert(connection_list.begin(), connected_tris.begin(), connected_tris.end());
4508 std::sort(connection_list.begin(), connection_list.end());
4511 for (
int tt = 1; tt < connection_list.size(); tt++) {
4512 if (connection_list.at(tt - 1) != connection_list.at(tt)) {
4516 int index = connection_list.at(tt - 1);
4518 if (fill_flag.at(index) == -1 && index != t) {
4520 fill_flag.at(index) = tag;
4522 if (depth < maxdepth) {
4523 floodfill(index, cloud_triangles, fill_flag, nodes, tag, depth + 1, maxdepth);
4541 if (printmessages) {
4542 cout <<
"Performing alphamask leaf reconstruction..." << flush;
4545 if (triangles.size() == 0) {
4546 std::cout <<
"failed." << std::endl;
4547 helios_runtime_error(
"ERROR (LiDARcloud::leafReconstructionAlphamask): There are no triangulated points. Either the triangulation failed or 'triangulateHitPoints()' was not called.");
4550 std::string file = mask_file;
4551 if (file.substr(file.find_last_of(
".") + 1) !=
"png") {
4552 std::cout <<
"failed." << std::endl;
4553 helios_runtime_error(
"ERROR (LiDARcloud::leafReconstructionAlphaMask): Mask data file " + std::string(mask_file) +
" must be PNG image format.");
4555 std::vector<std::vector<bool>> maskdata =
readPNGAlpha(mask_file);
4556 if (maskdata.size() == 0) {
4557 std::cout <<
"failed." << std::endl;
4558 helios_runtime_error(
"ERROR (LiDARcloud::leafReconstructionAlphaMask): Could not load mask file " + std::string(mask_file) +
". It contains no data.");
4560 int ix = maskdata.front().size();
4561 int jy = maskdata.size();
4565 for (
uint j = 0; j < masksize.
y; j++) {
4566 for (
uint i = 0; i < masksize.
x; i++) {
4568 if (maskdata.at(j).at(i)) {
4574 float solidfraction = float(Asolid) / float(Atotal);
4576 float total_area = 0.f;
4578 std::vector<std::vector<float>> group_areas;
4581 reconstructed_alphamasks_maskfile = mask_file;
4583 leafReconstructionFloodfill();
4587 uint group_count = reconstructed_triangles.size();
4589 float group_area_max = 0;
4591 std::vector<bool> group_filter_flag;
4592 group_filter_flag.resize(reconstructed_triangles.size());
4594 for (
int group = group_count - 1; group >= 0; group--) {
4598 for (
size_t t = 0; t < reconstructed_triangles.at(group).size(); t++) {
4600 float triangle_area = reconstructed_triangles.at(group).at(t).area;
4602 garea += triangle_area;
4605 if (garea < minimum_leaf_group_area || garea > maximum_leaf_group_area) {
4606 group_filter_flag.at(group) =
false;
4609 group_filter_flag.at(group) =
true;
4610 int cell = reconstructed_triangles.at(group).front().gridcell;
4611 group_areas.at(cell).push_back(garea);
4622 std::sort(group_areas.at(v).begin(), group_areas.at(v).end());
4625 if (group_areas.at(v).size() > Navg) {
4626 for (
int i = group_areas.at(v).size() - 1; i >= group_areas.at(v).size() - Navg; i--) {
4627 Lavg.at(v) += sqrtf(group_areas.at(v).at(i)) / float(Navg);
4629 }
else if (group_areas.at(v).size() == 0) {
4632 for (
int i = 0; i < group_areas.at(v).size(); i++) {
4633 Lavg.at(v) += sqrtf(group_areas.at(v).at(i)) / float(group_areas.at(v).size());
4637 if (printmessages) {
4638 std::cout <<
"Average leaf length for volume #" << v <<
" : " << Lavg.at(v) << endl;
4644 for (
int group = 0; group < reconstructed_triangles.size(); group++) {
4646 if (!group_filter_flag.at(group)) {
4650 int cell = reconstructed_triangles.at(group).front().gridcell;
4653 for (
int t = 0; t < reconstructed_triangles.at(group).size(); t++) {
4654 position = position + reconstructed_triangles.at(group).at(t).vertex0 / float(reconstructed_triangles.at(group).size());
4657 int gind = round(
randu() * (reconstructed_triangles.at(group).size() - 1));
4659 reconstructed_alphamasks_center.push_back(position);
4660 float l = Lavg.at(reconstructed_triangles.at(group).front().gridcell) * sqrt(leaf_aspect_ratio / solidfraction);
4661 float w = l / leaf_aspect_ratio;
4663 helios::vec3 normal =
cross(reconstructed_triangles.at(group).at(gind).vertex1 - reconstructed_triangles.at(group).at(gind).vertex0, reconstructed_triangles.at(group).at(gind).vertex2 - reconstructed_triangles.at(group).at(gind).vertex0);
4665 reconstructed_alphamasks_gridcell.push_back(reconstructed_triangles.at(group).front().gridcell);
4666 reconstructed_alphamasks_direct_flag.push_back(1);
4669 if (printmessages) {
4670 cout <<
"done." << endl;
4671 cout <<
"Directly reconstructed " << reconstructed_alphamasks_center.size() <<
" leaf groups." << endl;
4674 backfillLeavesAlphaMask(Lavg, leaf_aspect_ratio, solidfraction, group_filter_flag);
4676 for (
int group = 0; group < reconstructed_triangles.size(); group++) {
4678 if (!group_filter_flag.at(group)) {
4679 std::swap(reconstructed_triangles.at(group), reconstructed_triangles.back());
4680 reconstructed_triangles.pop_back();
4688void LiDARcloud::backfillLeavesAlphaMask(
const vector<float> &leaf_size,
float leaf_aspect_ratio,
float solidfraction,
const vector<bool> &group_filter_flag) {
4690 if (printmessages) {
4691 cout <<
"Backfilling leaves..." << endl;
4694 unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
4695 std::minstd_rand0 generator;
4696 generator.seed(seed);
4697 std::normal_distribution<float> randn;
4699 uint Ngroups = reconstructed_triangles.size();
4703 std::vector<std::vector<uint>> group_gridcell;
4704 group_gridcell.resize(Ncells);
4707 std::vector<float> leaf_area_current;
4708 leaf_area_current.resize(Ncells);
4712 for (
uint g = 0; g < Ngroups; g++) {
4713 if (group_filter_flag.at(g)) {
4714 if (reconstructed_triangles.at(g).size() > 0) {
4715 cell = reconstructed_triangles.at(g).front().gridcell;
4716 leaf_area_current.at(cell) += leaf_size.at(cell) * leaf_size.at(cell) * solidfraction;
4717 group_gridcell.at(cell).push_back(count);
4723 std::vector<int> deleted_groups;
4724 int backfill_count = 0;
4730 for (
uint v = 0; v < Ncells; v++) {
4734 float reconstruct_frac = (leaf_area_total - leaf_area_current.at(v)) / leaf_area_total;
4736 if (leaf_area_total == 0 || reconstructed_alphamasks_size.size() == 0) {
4737 backfill_warnings.
addWarning(
"volume_no_measured_leaf_area",
"skipping volume #" + std::to_string(v) +
" because it has no measured leaf area.");
4740 backfill_warnings.
addWarning(
"volume_no_triangles",
"skipping volume #" + std::to_string(v) +
" because it has no triangles.");
4742 }
else if (leaf_area_current.at(v) == 0) {
4744 std::vector<SphericalCoord> tri_rots;
4749 if (tri.gridcell == v) {
4750 helios::vec3 normal =
cross(tri.vertex1 - tri.vertex0, tri.vertex2 - tri.vertex0);
4755 while (leaf_area_current.at(v) < leaf_area_total) {
4757 int randi = round(
randu() * (tri_rots.size() - 1));
4765 reconstructed_alphamasks_center.push_back(shift);
4766 reconstructed_alphamasks_size.push_back(reconstructed_alphamasks_size.front());
4767 reconstructed_alphamasks_rotation.push_back(tri_rots.at(randi));
4768 reconstructed_alphamasks_gridcell.push_back(v);
4769 reconstructed_alphamasks_direct_flag.push_back(0);
4771 leaf_area_current.at(v) += reconstructed_alphamasks_size.back().x * reconstructed_alphamasks_size.back().y * solidfraction;
4775 }
else if (leaf_area_current.at(v) > leaf_area_total) {
4777 while (leaf_area_current.at(v) > leaf_area_total) {
4779 int randi = round(
randu() * (group_gridcell.at(v).size() - 1));
4781 int group_index = group_gridcell.at(v).at(randi);
4783 deleted_groups.push_back(group_index);
4785 leaf_area_current.at(v) -= reconstructed_alphamasks_size.at(group_index).x * reconstructed_alphamasks_size.at(group_index).y * solidfraction;
4790 while (leaf_area_current.at(v) < leaf_area_total) {
4792 int randi = round(
randu() * (group_gridcell.at(v).size() - 1));
4794 int group_index = group_gridcell.at(v).at(randi);
4802 helios::vec3 shift = reconstructed_alphamasks_center.at(group_index) +
helios::make_vec3(0.25 * randn(generator) * cellsize.
x, 0.25 * randn(generator) * cellsize.
y, 0.25 * randn(generator) * cellsize.
z);
4806 if (group_index >= reconstructed_alphamasks_center.size()) {
4807 helios_runtime_error(
"FAILED: " + std::to_string(group_index) +
" " + std::to_string(reconstructed_alphamasks_center.size()) +
" " + std::to_string(randi));
4808 }
else if (reconstructed_alphamasks_gridcell.at(group_index) != v) {
4812 reconstructed_alphamasks_center.push_back(shift);
4813 reconstructed_alphamasks_size.push_back(reconstructed_alphamasks_size.at(group_index));
4814 reconstructed_alphamasks_rotation.push_back(reconstructed_alphamasks_rotation.at(group_index));
4815 reconstructed_alphamasks_gridcell.push_back(v);
4816 reconstructed_alphamasks_direct_flag.push_back(0);
4818 leaf_area_current.at(v) += reconstructed_alphamasks_size.at(group_index).x * reconstructed_alphamasks_size.at(group_index).y * solidfraction;
4825 for (
uint v = 0; v < Ncells; v++) {
4829 float current_area = 0;
4830 for (
uint i = 0; i < reconstructed_alphamasks_size.size(); i++) {
4831 if (reconstructed_alphamasks_gridcell.at(i) == v) {
4832 current_area += reconstructed_alphamasks_size.at(i).x * reconstructed_alphamasks_size.at(i).y * solidfraction;
4837 if (printmessages) {
4838 cout <<
"Backfilled " << backfill_count <<
" total leaf groups." << endl;
4839 cout <<
"Deleted " << deleted_groups.size() <<
" total leaf groups." << endl;
4842 for (
int i = deleted_groups.size() - 1; i >= 0; i--) {
4843 int group_index = deleted_groups.at(i);
4844 if (group_index >= 0 && group_index < reconstructed_alphamasks_center.size()) {
4846 std::swap(reconstructed_alphamasks_center.at(group_index), reconstructed_alphamasks_center.back());
4847 reconstructed_alphamasks_center.pop_back();
4848 std::swap(reconstructed_alphamasks_size.at(group_index), reconstructed_alphamasks_size.back());
4849 reconstructed_alphamasks_size.pop_back();
4850 std::swap(reconstructed_alphamasks_rotation.at(group_index), reconstructed_alphamasks_rotation.back());
4851 reconstructed_alphamasks_rotation.pop_back();
4852 std::swap(reconstructed_alphamasks_gridcell.at(group_index), reconstructed_alphamasks_gridcell.back());
4853 reconstructed_alphamasks_gridcell.pop_back();
4854 std::swap(reconstructed_alphamasks_direct_flag.at(group_index), reconstructed_alphamasks_direct_flag.back());
4855 reconstructed_alphamasks_direct_flag.pop_back();
4859 backfill_warnings.
report(std::cerr);
4861 if (printmessages) {
4862 cout <<
"done." << endl;
4866void LiDARcloud::calculateLeafAngleCDF(
uint Nbins, std::vector<std::vector<float>> &CDF_theta, std::vector<std::vector<float>> &CDF_phi) {
4870 std::vector<std::vector<float>> PDF_theta, PDF_phi;
4871 CDF_theta.resize(Ncells);
4872 PDF_theta.resize(Ncells);
4873 CDF_phi.resize(Ncells);
4874 PDF_phi.resize(Ncells);
4875 for (
uint v = 0; v < Ncells; v++) {
4876 CDF_theta.at(v).resize(Nbins, 0.f);
4877 PDF_theta.at(v).resize(Nbins, 0.f);
4878 CDF_phi.at(v).resize(Nbins, 0.f);
4879 PDF_phi.at(v).resize(Nbins, 0.f);
4881 float db_theta = 0.5 *
M_PI / Nbins;
4882 float db_phi = 2.f *
M_PI / Nbins;
4885 for (
size_t t = 0; t < triangles.size(); t++) {
4886 float triangle_area = triangles.at(t).area;
4887 int gridcell = triangles.at(t).gridcell;
4889 if (gridcell >= 0 && gridcell < (
int) Ncells) {
4890 helios::vec3 normal =
cross(triangles.at(t).vertex1 - triangles.at(t).vertex0, triangles.at(t).vertex2 - triangles.at(t).vertex0);
4891 normal.z = fabs(normal.z);
4895 int bin_theta = floor(normal_dir.
zenith / db_theta);
4896 if (bin_theta >= Nbins) {
4897 bin_theta = Nbins - 1;
4900 int bin_phi = floor(normal_dir.
azimuth / db_phi);
4901 if (bin_phi >= Nbins) {
4902 bin_phi = Nbins - 1;
4905 PDF_theta.at(gridcell).at(bin_theta) += triangle_area;
4906 PDF_phi.at(gridcell).at(bin_phi) += triangle_area;
4911 for (
uint v = 0; v < Ncells; v++) {
4912 for (
uint i = 0; i < Nbins; i++) {
4913 for (
uint j = 0; j <= i; j++) {
4914 CDF_theta.at(v).at(i) += PDF_theta.at(v).at(j);
4915 CDF_phi.at(v).at(i) += PDF_phi.at(v).at(j);
4939 std::vector<helios::vec3> grid_vertices;
4942 grid_vertices.push_back(boxmin);
4943 grid_vertices.push_back(boxmax);
4951 float max_theta = 0;
4952 float min_theta =
M_PI;
4954 float min_phi = 2 *
M_PI;
4955 for (
uint gg = 0; gg < grid_vertices.size(); gg++) {
4967 if (sc.
zenith < min_theta) {
4971 if (sc.
zenith > max_theta) {
4982 float this_theta = raydir.
zenith;
4983 float this_phi = raydir.
azimuth;
4984 double this_phi_d = double(this_phi);
4986 if (this_phi < phi_range.x || this_phi > phi_range.
y || this_theta < theta_range.x || this_theta > theta_range.
y) {
4995void LiDARcloud::computeGtheta(
uint Ncells,
uint Nscans, std::vector<float> &Gtheta, std::vector<float> &Gtheta_bar) {
4998 Gtheta.resize(Ncells, 0.f);
4999 Gtheta_bar.resize(Ncells, 0.f);
5003 std::vector<float> denom_sum;
5004 denom_sum.resize(Ncells, 0.f);
5005 std::vector<uint> cell_tri_count;
5006 cell_tri_count.resize(Ncells, 0);
5009 for (
size_t t = 0; t < Ntri; t++) {
5012 int cell = tri.gridcell;
5014 if (cell >= 0 && cell < Ncells) {
5029 float S = 0.5f * (L0 + L1 + L2);
5030 float area = sqrt(S * (S - L0) * (S - L1) * (S - L2));
5044 if (std::isfinite(area) && area > 0.f) {
5045 float normal_dot_ray = fabs(normal * raydir);
5046 Gtheta.at(cell) += normal_dot_ray * area * fabs(sin(theta));
5047 denom_sum.at(cell) += fabs(sin(theta)) * area;
5048 cell_tri_count.at(cell) += 1;
5054 for (
uint v = 0; v < Ncells; v++) {
5055 if (cell_tri_count[v] > 0) {
5056 Gtheta[v] = Gtheta[v] / denom_sum[v];
5057 Gtheta_bar[v] += Gtheta[v] / float(Nscans);
5062bool LiDARcloud::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) {
5065 if (Gtheta == 0 || Gtheta != Gtheta) {
5070 if (dr_samples.size() < min_voxel_hits) {
5084 const float P_error_denom = fmax(P, 1e-6f);
5092 for (
size_t j = 0; j < dr_samples.size(); j++) {
5093 mean += exp(-a * dr_samples[j] * Gtheta);
5095 mean /= float(dr_samples.size());
5096 float error = fabs(
mean - P) / P_error_denom;
5104 while (error > etol && iter < maxiter) {
5110 for (
size_t j = 0; j < dr_samples.size(); j++) {
5111 mean += exp(-a * dr_samples[j] * Gtheta);
5113 mean /= float(dr_samples.size());
5114 error = fabs(
mean - P) / P_error_denom;
5118 if (error == eold) {
5123 a = fabs((aold * error - a * eold) / (error - eold));
5128 float dr_bar = 0.0f;
5129 for (
size_t i = 0; i < dr_samples.size(); i++) {
5130 dr_bar += dr_samples[i];
5132 dr_bar /= float(dr_samples.size());
5140 bool converged = (error <= etol && a == a && a <= 100);
5141 bool used_fallback =
false;
5144 warnings.
addWarning(
"invertLAD_did_not_converge",
"LAD inversion failed for volume #" + std::to_string(voxel_index) +
". Using average dr formulation.");
5145 a = (1.f - P) / (dr_bar * Gtheta);
5146 used_fallback =
true;
5151 a = fmin((1.f - P) / dr_bar / Gtheta, -log(P_error_denom) / dr_bar / Gtheta);
5155 leaf_area = a * gridsize.
x * gridsize.
y * gridsize.
z;
5160LiDARcloud::LADInversionResult LiDARcloud::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,
5174 LADInversionResult result;
5175 result.beam_count = (int) dr_samples.size();
5176 result.I_rdi = 1.f - P;
5179 float leaf_area = 0.f;
5180 bool ok = invertLAD(voxel_index, P, Gtheta, dr_samples, min_voxel_hits, gridsize, leaf_area, warnings);
5181 result.leaf_area = leaf_area;
5182 result.converged = ok;
5185 const int N = result.beam_count;
5188 for (
float d: dr_samples) {
5191 result.zbar_e =
sum / float(N);
5193 for (
float d: dr_samples) {
5194 ss += (d - result.zbar_e) * (d - result.zbar_e);
5196 result.var_path = ss / float(N);
5200 const float dr_bar = result.zbar_e;
5201 if (!ok || N < min_voxel_hits || Gtheta <= 0.f || Gtheta != Gtheta || dr_bar <= 0.f) {
5202 result.LAD_variance = -1.f;
5206 const float I = result.I_rdi;
5208 const float I_b = std::min(I, 1.f - 1.f / (2.f *
float(N) + 2.f));
5217 const float binomial_varI = I_b * (1.f - I_b) /
float(N);
5218 float empirical_varI = (sum_frac_sq / float(N) - P * P) /
float(N);
5219 if (empirical_varI < 0.f) {
5220 empirical_varI = 0.f;
5222 const float var_I_a = std::max(binomial_varI, empirical_varI);
5229 float var_I_b = 0.f;
5230 if (element_width > 0.f) {
5231 const float volume = gridsize.
x * gridsize.
y * gridsize.
z;
5232 const float delta = std::cbrt(volume);
5233 const float L1 = (float) (M_PI) *element_width * element_width / (8.f * delta * delta);
5234 result.L1_element = L1;
5235 result.element_size_known =
true;
5236 if (L1 < 0.3f && I > 0.f && I < 1.f) {
5237 var_I_b = 0.23f * L1 * (1.f - I) * std::pow(I, 1.9f - 2.3f * L1);
5238 if (var_I_b < 0.f) {
5245 const float denom = (1.f - I_b) * (1.f - I_b) * dr_bar * dr_bar * Gtheta * Gtheta;
5246 result.LAD_variance = (var_I_a + var_I_b) / denom;
5251bool LiDARcloud::ciValidPimont(
float L,
float L1,
int N,
float confidence_level)
const {
5257 const float l1 = (L1 < 0.f) ? 0.f : L1;
5258 if (confidence_level > 0.925f) {
5259 return (L <= 2.0f && l1 <= 0.05f && N >= 30) || (L <= 2.5f && l1 <= 0.01f && N >= 150) || (l1 <= 0.05f && N >= 150);
5262 return (L >= 0.5f && L <= 2.0f && l1 <= 0.05f && N >= 40) || (l1 <= 0.01f && N >= 100) || (l1 <= 0.05f && N >= 200);
5278 calculateLeafArea_inner(
context, min_voxel_hits, element_width, -1.f);
5283 if (!(Gtheta > 0.f) || Gtheta > 1.f) {
5284 helios_runtime_error(
"ERROR (LiDARcloud::calculateLeafArea): The supplied G(theta) must be in the range (0,1], but " + std::to_string(Gtheta) +
" was provided. Use 0.5 for a spherical (random) leaf-angle distribution.");
5286 calculateLeafArea_inner(
context, min_voxel_hits, element_width, Gtheta);
5289void LiDARcloud::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,
5290 std::vector<float> &dr_array_cell) {
5292 float E_before = 0, E_inside = 0, E_after = 0;
5298 for (
size_t r = 0; r < Nreturns; r++) {
5299 uint local_index = return_indices[r];
5300 if (dr[local_index] > 0) {
5301 drr += dr[local_index];
5305 if (hit_location[local_index] == 1)
5307 else if (hit_location[local_index] == 2)
5309 else if (hit_location[local_index] == 3)
5317 if (E_inside != 0 || E_after != 0) {
5318 float frac = E_after / (E_inside + E_after);
5319 P_equal_numerator += frac;
5320 P_equal_sumsq += frac * frac;
5321 P_equal_denominator += 1;
5326 float drrx = drr / float(dr_count);
5327 dr_array_cell.push_back(drrx);
5336 float tx_min = (voxel_min.
x - origin.
x) / direction.
x;
5337 float tx_max = (voxel_max.x - origin.
x) / direction.x;
5338 if (tx_min > tx_max)
5339 std::swap(tx_min, tx_max);
5341 float ty_min = (voxel_min.y - origin.
y) / direction.y;
5342 float ty_max = (voxel_max.y - origin.
y) / direction.y;
5343 if (ty_min > ty_max)
5344 std::swap(ty_min, ty_max);
5346 float tz_min = (voxel_min.z - origin.
z) / direction.z;
5347 float tz_max = (voxel_max.z - origin.
z) / direction.z;
5348 if (tz_min > tz_max)
5349 std::swap(tz_min, tz_max);
5351 t0 = std::max({tx_min, ty_min, tz_min});
5352 t1 = std::min({tx_max, ty_max, tz_max});
5359 const float o[3] = {origin.
x, origin.
y, origin.
z};
5360 const float d[3] = {direction.
x, direction.
y, direction.
z};
5361 const float lo[3] = {grid_min.
x, grid_min.
y, grid_min.
z};
5362 const float hi[3] = {grid_max.
x, grid_max.
y, grid_max.
z};
5364 float tmin = -std::numeric_limits<float>::max();
5365 float tmax = std::numeric_limits<float>::max();
5366 for (
int ax = 0; ax < 3; ax++) {
5367 if (fabs(d[ax]) < 1e-9f) {
5368 if (o[ax] < lo[ax] || o[ax] > hi[ax])
5371 float t1 = (lo[ax] - o[ax]) / d[ax];
5372 float t2 = (hi[ax] - o[ax]) / d[ax];
5375 tmin = std::max(tmin, t1);
5376 tmax = std::min(tmax, t2);
5386LiDARcloud::VoxelLattice LiDARcloud::detectVoxelLattice()
const {
5388 VoxelLattice lattice;
5396 const GridCell &ref = grid_cells.front();
5398 if (count.x <= 0 || count.y <= 0 || count.z <= 0) {
5404 const size_t expected_cells = (size_t) count.x * (
size_t) count.y * (size_t) count.z;
5405 if ((
size_t) Ncells != expected_cells) {
5409 const helios::vec3 cell_extent =
make_vec3(ref.global_size.
x /
float(count.x), ref.global_size.
y /
float(count.y), ref.global_size.
z /
float(count.z));
5410 const helios::vec3 lattice_origin = ref.global_anchor - ref.global_size * 0.5f;
5413 const float pos_tol = 1e-4f * std::max({cell_extent.
x, cell_extent.
y, cell_extent.
z, 1e-6f});
5414 const float rot_tol = 1e-6f;
5416 std::vector<int> ijk_to_index(expected_cells, -1);
5418 for (
uint c = 0; c < Ncells; c++) {
5419 const GridCell &cell = grid_cells.at(c);
5422 if (cell.global_count.
x != count.x || cell.global_count.
y != count.y || cell.global_count.
z != count.z) {
5425 if ((cell.global_anchor - ref.global_anchor).magnitude() > pos_tol || (cell.global_size - ref.global_size).magnitude() > pos_tol || fabs(cell.
azimuthal_rotation - ref.
azimuthal_rotation) > rot_tol) {
5430 if (fabs(cell.size.
x - cell_extent.
x) > pos_tol || fabs(cell.size.
y - cell_extent.
y) > pos_tol || fabs(cell.size.
z - cell_extent.
z) > pos_tol) {
5436 if (ijk.
x < 0 || ijk.
x >= count.x || ijk.
y < 0 || ijk.
y >= count.y || ijk.
z < 0 || ijk.
z >= count.z) {
5439 const size_t flat = ((size_t) ijk.
z * count.y + ijk.
y) * count.x + ijk.
x;
5440 if (ijk_to_index[flat] != -1) {
5443 ijk_to_index[flat] = (int) c;
5447 const helios::vec3 expected_center = lattice_origin +
make_vec3((
float(ijk.
x) + 0.5f) * cell_extent.
x, (
float(ijk.
y) + 0.5f) * cell_extent.
y, (
float(ijk.
z) + 0.5f) * cell_extent.
z);
5448 if ((cell.center - expected_center).magnitude() > pos_tol) {
5453 lattice.valid =
true;
5454 lattice.origin = lattice_origin;
5455 lattice.anchor = ref.global_anchor;
5456 lattice.cell_extent = cell_extent;
5458 lattice.count = count;
5459 lattice.ijk_to_index = std::move(ijk_to_index);
5464void LiDARcloud::calculateLeafArea_inner(
helios::Context *
context,
int min_voxel_hits,
float element_width,
float supplied_Gtheta) {
5466 const bool use_supplied_Gtheta = (supplied_Gtheta > 0.f);
5468 if (printmessages) {
5469 std::cout <<
"Calculating leaf area (CollisionDetection)..." << std::endl;
5474 if (!use_supplied_Gtheta && !triangulationcomputed) {
5475 helios_runtime_error(
"ERROR (LiDARcloud::calculateLeafAreaCD): Triangulation must be performed prior to leaf area calculation. See triangulateHitPoints(). For scans that cannot be triangulated (e.g. moving-platform scans), use the "
5476 "calculateLeafArea overload that takes a G(theta) argument.");
5479 if (!hitgridcellcomputed) {
5500 "ERROR (LiDARcloud::calculateLeafArea): No miss points found in the point cloud. Leaf area inversion requires fired pulses that returned nothing (misses) in order to count transmitted beams. Provide a scan format that retains misses, or call gapfillMisses() to synthesize them, before calling calculateLeafArea().");
5503 if (printmessages) {
5504 if (isMultiReturnData()) {
5505 std::cout <<
"Multi-return data detected - using beam-based equal weighting algorithm (CD)" << std::endl;
5507 std::cout <<
"Single-return data with misses - using beam-based equal weighting algorithm (CD)" << std::endl;
5515 std::vector<std::vector<float>> P_equal_numerator_array(Ncells);
5516 std::vector<std::vector<float>> P_equal_denominator_array(Ncells);
5520 std::vector<std::vector<float>> P_equal_sumsq_array(Ncells);
5521 std::vector<std::vector<float>> dr_array(Ncells);
5524 std::vector<std::vector<float>> dr_agg;
5525 dr_agg.resize(Ncells);
5526 std::vector<float> Gtheta_bar;
5527 Gtheta_bar.resize(Ncells, 0.f);
5533 const VoxelLattice lattice = detectVoxelLattice();
5536 for (
uint s = 0; s < Nscans; s++) {
5539 std::vector<helios::vec3> this_scan_xyz;
5540 std::vector<uint> this_scan_index;
5545 std::vector<uint> global_to_local(
getHitCount(), 0);
5548 global_to_local[(
uint) r] = (
uint) this_scan_xyz.size();
5550 this_scan_index.push_back(r);
5553 size_t Nhits = this_scan_xyz.size();
5560 BeamGrouping beams = groupHitsByTimestamp(this_scan_index);
5561 uint Nbeams = beams.Nbeams;
5567 std::vector<helios::vec3> this_scan_origin(Nhits);
5568 for (
size_t i = 0; i < Nhits; i++) {
5569 this_scan_origin[i] =
getHitOrigin(this_scan_index[i]);
5575 std::vector<uint> beam_members_local(beams.beam_members.size());
5576 for (
size_t m = 0; m < beams.beam_members.size(); m++) {
5577 beam_members_local[m] = global_to_local[beams.beam_members[m]];
5580 if (lattice.valid && !force_bruteforce_LAD) {
5589 auto process_beam = [&](
uint k, std::vector<float> &P_num, std::vector<float> &P_denom, std::vector<float> &P_sumsq, std::vector<std::vector<float>> &dr_cells, std::vector<float> &ret_dist, std::vector<float> &dr_cell,
5590 std::vector<uint> &hl_cell, std::vector<uint> &local_seq) {
5591 const uint beam_start = beams.beam_offsets[k];
5592 const size_t Nret = beams.beam_offsets[k + 1] - beam_start;
5597 helios::vec3 origin = this_scan_origin[beam_members_local[beam_start]];
5603 if (fabs(lattice.rotation) > 1e-6f) {
5604 origin_L =
rotatePointAboutLine(origin - lattice.anchor, helios::make_vec3(0, 0, 0), helios::make_vec3(0, 0, 1), -lattice.rotation) + lattice.anchor;
5608 float max_hit_distance = -1.f;
5610 ret_dist.resize(Nret);
5611 for (
size_t j = 0; j < Nret; j++) {
5612 uint i = beam_members_local[beam_start + j];
5614 if (fabs(lattice.rotation) > 1e-6f) {
5615 hit_L =
rotatePointAboutLine(hit_L - lattice.anchor, helios::make_vec3(0, 0, 0), helios::make_vec3(0, 0, 1), -lattice.rotation) + lattice.anchor;
5620 if (dist > max_hit_distance) {
5621 max_hit_distance = dist;
5625 if (max_hit_distance <= 0.f)
5631 helios::vec3 grid_max = lattice.origin +
make_vec3(lattice.cell_extent.x * lattice.count.x, lattice.cell_extent.y * lattice.count.y, lattice.cell_extent.z * lattice.count.z);
5633 float t_enter, t_exit;
5634 if (!rayGridIntersect(origin_L, direction, grid_min, grid_max, t_enter, t_exit))
5637 float t_start = std::max(t_enter, 0.f);
5638 if (t_exit <= 1e-6f)
5647 const float origin_arr[3] = {origin_L.
x, origin_L.y, origin_L.z};
5648 const float dir_arr[3] = {direction.
x, direction.
y, direction.
z};
5649 const float entry_arr[3] = {entry.
x, entry.
y, entry.
z};
5650 const float gmin_arr[3] = {grid_min.
x, grid_min.
y, grid_min.
z};
5651 const float extent_arr[3] = {lattice.cell_extent.x, lattice.cell_extent.y, lattice.cell_extent.z};
5652 const int count_arr[3] = {lattice.count.x, lattice.count.y, lattice.count.z};
5653 for (
int ax = 0; ax < 3; ax++) {
5654 int idx = (int) std::floor((entry_arr[ax] - gmin_arr[ax]) / extent_arr[ax]);
5657 if (idx >= count_arr[ax])
5658 idx = count_arr[ax] - 1;
5660 if (fabs(dir_arr[ax]) < 1e-9f) {
5663 tMax[ax] = std::numeric_limits<float>::max();
5664 tDelta[ax] = std::numeric_limits<float>::max();
5666 step[ax] = (dir_arr[ax] > 0) ? 1 : -1;
5667 float next_boundary = gmin_arr[ax] + float(idx + (step[ax] > 0 ? 1 : 0)) * extent_arr[ax];
5668 tMax[ax] = (next_boundary - origin_arr[ax]) / dir_arr[ax];
5669 tDelta[ax] = extent_arr[ax] / fabs(dir_arr[ax]);
5675 dr_cell.assign(Nret, 0.f);
5676 hl_cell.assign(Nret, 0);
5677 local_seq.resize(Nret);
5678 for (
uint j = 0; j < Nret; j++)
5681 const size_t max_steps = (size_t) count_arr[0] + count_arr[1] + count_arr[2] + 3;
5682 for (
size_t stepcount = 0; stepcount <= max_steps; stepcount++) {
5683 if (ijk[0] < 0 || ijk[0] >= count_arr[0] || ijk[1] < 0 || ijk[1] >= count_arr[1] || ijk[2] < 0 || ijk[2] >= count_arr[2])
5686 const size_t flat = ((size_t) ijk[2] * count_arr[1] + ijk[1]) * count_arr[0] + ijk[0];
5687 int cell_index = lattice.ijk_to_index[flat];
5688 if (cell_index >= 0) {
5691 helios::vec3 cmin = lattice.origin +
make_vec3(ijk[0] * lattice.cell_extent.x, ijk[1] * lattice.cell_extent.y, ijk[2] * lattice.cell_extent.z);
5695 if (cellSlab(origin_L, direction, cmin, cmax, ct0, ct1) && ct1 > 1e-6f) {
5696 float drval = fabs(ct1 - ct0);
5698 for (
size_t j = 0; j < Nret; j++) {
5699 float hd = ret_dist[j];
5701 if (hd >= ct0 && hd <= ct1)
5708 accumulateBeamCell(local_seq.data(), Nret, dr_cell, hl_cell, P_num[cell_index], P_denom[cell_index], P_sumsq[cell_index], dr_cells[cell_index]);
5714 if (tMax[1] < tMax[axis])
5716 if (tMax[2] < tMax[axis])
5718 if (tMax[axis] > t_exit)
5720 if (step[axis] == 0)
5722 ijk[axis] += step[axis];
5723 tMax[axis] += tDelta[axis];
5729 std::vector<float> P_num_scratch(Ncells, 0.f);
5730 std::vector<float> P_denom_scratch(Ncells, 0.f);
5731 std::vector<float> P_sumsq_scratch(Ncells, 0.f);
5732 std::vector<std::vector<float>> dr_scratch(Ncells);
5740 int num_threads = omp_get_max_threads();
5742 int num_threads = 1;
5744 std::vector<std::vector<float>> P_num_thread(num_threads, std::vector<float>(Ncells, 0.f));
5745 std::vector<std::vector<float>> P_denom_thread(num_threads, std::vector<float>(Ncells, 0.f));
5746 std::vector<std::vector<float>> P_sumsq_thread(num_threads, std::vector<float>(Ncells, 0.f));
5747 std::vector<std::vector<std::vector<float>>> dr_thread(num_threads, std::vector<std::vector<float>>(Ncells));
5752 int tid = omp_get_thread_num();
5757 std::vector<float> ret_dist, dr_cell;
5758 std::vector<uint> hl_cell, local_seq;
5759#pragma omp for schedule(dynamic, 256)
5760 for (
int k = 0; k < static_cast<int>(Nbeams); k++) {
5761 process_beam((
uint) k, P_num_thread[tid], P_denom_thread[tid], P_sumsq_thread[tid], dr_thread[tid], ret_dist, dr_cell, hl_cell, local_seq);
5766 for (
int t = 0; t < num_threads; t++) {
5767 for (
uint c = 0; c < Ncells; c++) {
5768 P_num_scratch[c] += P_num_thread[t][c];
5769 P_denom_scratch[c] += P_denom_thread[t][c];
5770 P_sumsq_scratch[c] += P_sumsq_thread[t][c];
5771 for (
float v: dr_thread[t][c])
5772 dr_scratch[c].push_back(v);
5777 for (
uint c = 0; c < Ncells; c++) {
5778 P_equal_numerator_array.at(c).push_back(P_num_scratch[c]);
5779 P_equal_denominator_array.at(c).push_back(P_denom_scratch[c]);
5780 P_equal_sumsq_array.at(c).push_back(P_sumsq_scratch[c]);
5781 for (
float v: dr_scratch[c])
5782 dr_array.at(c).push_back(v);
5790 std::vector<float> dr(Nhits, 0.0f);
5791 std::vector<uint> hit_location(Nhits, 0);
5794 for (
uint c = 0; c < Ncells; c++) {
5802 std::fill(dr.begin(), dr.end(), 0.0f);
5803 std::fill(hit_location.begin(), hit_location.end(), 0);
5806#pragma omp parallel for
5807 for (
int i = 0; i < static_cast<int>(Nhits); i++) {
5813 if (fabs(rotation) > 1e-6f) {
5814 hit_xyz =
rotatePointAboutLine(hit_xyz - anchor, helios::make_vec3(0, 0, 0), helios::make_vec3(0, 0, 1), -rotation) + anchor;
5815 origin =
rotatePointAboutLine(origin - anchor, helios::make_vec3(0, 0, 0), helios::make_vec3(0, 0, 1), -rotation) + anchor;
5820 float hit_distance = direction.
magnitude();
5828 float tx_min = (voxel_min.
x - origin.x) / direction.
x;
5829 float tx_max = (voxel_max.x - origin.x) / direction.
x;
5830 if (tx_min > tx_max)
5831 std::swap(tx_min, tx_max);
5833 float ty_min = (voxel_min.y - origin.y) / direction.
y;
5834 float ty_max = (voxel_max.y - origin.y) / direction.
y;
5835 if (ty_min > ty_max)
5836 std::swap(ty_min, ty_max);
5838 float tz_min = (voxel_min.z - origin.
z) / direction.
z;
5839 float tz_max = (voxel_max.z - origin.
z) / direction.
z;
5840 if (tz_min > tz_max)
5841 std::swap(tz_min, tz_max);
5843 float t0 = std::max({tx_min, ty_min, tz_min});
5844 float t1 = std::min({tx_max, ty_max, tz_max});
5860 if (t0 < t1 && t1 > 1e-6f) {
5861 dr[i] = fabs(t1 - t0);
5863 if (hit_distance >= t0 && hit_distance <= t1) {
5864 hit_location[i] = 2;
5865 }
else if (hit_distance > t1) {
5866 hit_location[i] = 3;
5867 }
else if (hit_distance < t0) {
5868 hit_location[i] = 1;
5874 float P_equal_numerator = 0;
5875 float P_equal_denominator = 0;
5876 float P_equal_sumsq = 0;
5878 for (
uint k = 0; k < Nbeams; k++) {
5879 accumulateBeamCell(&beam_members_local[beams.beam_offsets[k]], beams.beamSize(k), dr, hit_location, P_equal_numerator, P_equal_denominator, P_equal_sumsq, dr_array.at(c));
5882 P_equal_numerator_array.at(c).push_back(P_equal_numerator);
5883 P_equal_denominator_array.at(c).push_back(P_equal_denominator);
5884 P_equal_sumsq_array.at(c).push_back(P_equal_sumsq);
5891 std::vector<float> Gtheta;
5892 if (use_supplied_Gtheta) {
5893 Gtheta.assign(Ncells, supplied_Gtheta);
5895 computeGtheta(Ncells, Nscans, Gtheta, Gtheta_bar);
5899 if (printmessages) {
5900 std::cout <<
"Inverting to find LAD..." << std::flush;
5904 invertLAD_warnings.
setEnabled(printmessages);
5906 for (
uint v = 0; v < Ncells; v++) {
5909 float P_num_sum = 0.0f, P_denom_sum = 0.0f, P_sumsq_sum = 0.0f;
5910 for (
uint s = 0; s < P_equal_numerator_array[v].size(); s++) {
5911 P_num_sum += P_equal_numerator_array[v][s];
5912 P_denom_sum += P_equal_denominator_array[v][s];
5913 P_sumsq_sum += P_equal_sumsq_array[v][s];
5915 if (P_denom_sum > 0) {
5916 P = P_num_sum / P_denom_sum;
5920 for (
uint s = 0; s < dr_array[v].size(); s++) {
5921 if (dr_array[v][s] > 0) {
5922 dr_agg[v].push_back(dr_array[v][s]);
5927 if (dr_agg[v].size() < min_voxel_hits) {
5930 grid_cells.at(v).beam_count = (int) dr_agg[v].size();
5931 grid_cells.at(v).LAD_variance = -1.f;
5932 grid_cells.at(v).ci_valid =
false;
5940 LADInversionResult inv = invertLADWithVariance(v, P, Gtheta[v], dr_agg[v], P_sumsq_sum, element_width, min_voxel_hits, gridsize, invertLAD_warnings);
5947 cell.
I_rdi = inv.I_rdi;
5948 cell.
zbar_e = inv.zbar_e;
5956 const float volume = gridsize.
x * gridsize.
y * gridsize.
z;
5957 const float a = (volume > 0.f) ? inv.leaf_area / volume : 0.f;
5958 const float L = a * Gtheta[v] * inv.zbar_e;
5959 cell.
ci_valid = (inv.LAD_variance >= 0.f) && ciValidPimont(L, inv.L1_element, inv.beam_count, 0.95f);
5962 invertLAD_warnings.
report(std::cerr);
5964 if (printmessages) {
5965 std::cout <<
"done." << std::endl;
5980 if (collision_detection !=
nullptr) {
5986 if (collision_detection !=
nullptr) {
6003 if (printmessages) {
6004 std::cout <<
"Grouping hit points by grid cell (CPU)..." << std::flush;
6010 if (total_hits == 0) {
6011 std::cout <<
"WARNING (calculateHitGridCellCD): There are no hits currently in the point cloud. Skipping grid cell binning calculation." << std::endl;
6019 std::vector<helios::vec3> cell_min(Ncells), cell_max(Ncells), cell_anchor(Ncells);
6020 std::vector<float> cell_rotation(Ncells);
6021 std::vector<bool> cell_rotated(Ncells);
6022 for (
uint c = 0; c < Ncells; c++) {
6025 cell_min[c] = center - size * 0.5f;
6026 cell_max[c] = center + size * 0.5f;
6029 cell_rotated[c] = (fabs(cell_rotation[c]) > 1e-6f);
6033#pragma omp parallel for schedule(dynamic, 1000)
6034 for (
int r = 0; r < static_cast<int>(total_hits); r++) {
6037 int assigned_cell = -1;
6044 for (
uint c = 0; c < Ncells; c++) {
6048 if (cell_rotated[c]) {
6054 if (p.
x >= lo.
x && p.
x <= hi.
x && p.
y >= lo.
y && p.
y <= hi.
y && p.
z >= lo.
z && p.
z <= hi.
z) {
6064 if (printmessages) {
6065 std::cout <<
"done." << std::endl;
6068 hitgridcellcomputed =
true;
6071bool LiDARcloud::isMultiReturnData()
const {
6078 helios_runtime_error(
"ERROR (isMultiReturnData): Multi-return data detected (target_count > 1) but 'timestamp' field is missing. Cannot group hits into beams.");
6082 helios_runtime_error(
"ERROR (isMultiReturnData): Multi-return data detected (target_count > 1) but 'target_index' field is missing. Cannot filter first returns for triangulation.");
6120LiDARcloud::BeamGrouping LiDARcloud::groupHitsByTimestamp(
const std::vector<uint> &scan_indices)
const {
6122 BeamGrouping result;
6124 if (scan_indices.empty()) {
6129 const size_t N = scan_indices.size();
6134 bool has_timestamp =
true;
6135 for (
uint idx: scan_indices) {
6137 has_timestamp =
false;
6141 if (!has_timestamp) {
6143 result.Nbeams = (
uint) N;
6144 result.beam_members = scan_indices;
6145 result.beam_offsets.resize(N + 1);
6146 for (
size_t i = 0; i <= N; i++) {
6147 result.beam_offsets[i] = (
uint) i;
6154 std::vector<double> timestamps(N);
6155 for (
size_t i = 0; i < N; i++) {
6156 timestamps[i] =
getHitData(scan_indices[i],
"timestamp");
6158 std::vector<uint> order(N);
6159 for (
size_t i = 0; i < N; i++) {
6160 order[i] = (
uint) i;
6162 std::sort(order.begin(), order.end(), [&](
uint a,
uint b) { return timestamps[a] < timestamps[b]; });
6166 result.beam_members.resize(N);
6167 result.beam_offsets.clear();
6168 result.beam_offsets.push_back(0);
6169 double previous_time = 0.0;
6170 for (
size_t i = 0; i < N; i++) {
6172 result.beam_members[i] = scan_indices[si];
6174 previous_time = timestamps[si];
6175 }
else if (timestamps[si] != previous_time) {
6176 result.beam_offsets.push_back((
uint) i);
6177 previous_time = timestamps[si];
6180 result.beam_offsets.push_back((
uint) N);
6181 result.Nbeams = (
uint) result.beam_offsets.size() - 1;
6188 std::vector<uint> UUIDs_all =
context->getAllUUIDs();
6189 const uint N = UUIDs_all.size();
6193 std::vector<int> prim_vol(N, -1);
6196#pragma omp parallel for
6197 for (
int p = 0; p < static_cast<int>(N); p++) {
6198 std::vector<helios::vec3> verts =
context->getPrimitiveVertices(UUIDs_all[p]);
6202 for (
uint c = 0; c < Ncells; c++) {
6210 if (fabs(rotation) > 1e-6f) {
6220 float x0 = center.
x - 0.5f * size.
x;
6221 float x1 = center.
x + 0.5f * size.
x;
6222 float y0 = center.
y - 0.5f * size.
y;
6223 float y1 = center.
y + 0.5f * size.
y;
6224 float z0 = center.
z - 0.5f * size.z;
6225 float z1 = center.
z + 0.5f * size.z;
6228 float tx_min = (x0 - origin_pt.
x) / direction.
x;
6229 float tx_max = (x1 - origin_pt.
x) / direction.
x;
6230 if (tx_min > tx_max)
6231 std::swap(tx_min, tx_max);
6233 float ty_min = (y0 - origin_pt.
y) / direction.
y;
6234 float ty_max = (y1 - origin_pt.
y) / direction.
y;
6235 if (ty_min > ty_max)
6236 std::swap(ty_min, ty_max);
6238 float tz_min = (z0 - origin_pt.
z) / direction.
z;
6239 float tz_max = (z1 - origin_pt.
z) / direction.
z;
6240 if (tz_min > tz_max)
6241 std::swap(tz_min, tz_max);
6243 float t0 = std::max({tx_min, ty_min, tz_min});
6244 float t1 = std::min({tx_max, ty_max, tz_max});
6247 if (t0 < t1 && t1 > 1e-6f) {
6248 float T = (prim_xyz_rot - origin_pt).magnitude();
6249 if (T >= t0 && T <= t1) {
6258 std::vector<float> total_area(Ncells, 0.f);
6259 for (
size_t p = 0; p < N; p++) {
6260 if (prim_vol[p] >= 0) {
6261 uint gridcell = prim_vol[p];
6262 total_area[gridcell] +=
context->getPrimitiveArea(UUIDs_all[p]);
6263 context->setPrimitiveData(UUIDs_all[p],
"gridCell", gridcell);
6268 std::vector<float> output_LeafArea(Ncells);
6269 for (
uint v = 0; v < Ncells; v++) {
6270 output_LeafArea[v] = total_area[v];
6277 for (
size_t p = 0; p < N; p++) {
6278 if (prim_vol[p] >= 0) {
6279 context->setPrimitiveData(UUIDs_all[p],
"synthetic_leaf_area", total_area[prim_vol[p]]);
6283 return output_LeafArea;
6309 std::vector<std::vector<float>> detectReturnsFromSubrays(std::vector<std::vector<float>> &t_pulse,
float total_pulse_weight,
int Npulse,
float range_resolution,
float detection_threshold,
int max_returns,
6312 std::vector<std::vector<float>> t_hit;
6313 if (t_pulse.empty()) {
6317 std::sort(t_pulse.begin(), t_pulse.end(), [](
const std::vector<float> &a,
const std::vector<float> &b) { return a[0] < b[0]; });
6321 const float denom = (total_pulse_weight > 0.f) ? total_pulse_weight : float(Npulse);
6327 while (i < t_pulse.size()) {
6328 const float t0 = t_pulse[i][0];
6329 double sum_w = 0.0, sum_wt = 0.0, sum_wt2 = 0.0, sum_wcos = 0.0;
6331 float lastID = t_pulse[i][2];
6333 while (j < t_pulse.size() && (t_pulse[j][0] - t0) <= range_resolution) {
6334 const float t = t_pulse[j][0];
6335 const float cosval = t_pulse[j][1];
6336 const float w = t_pulse[j][3];
6338 sum_wt += double(w) * t;
6339 sum_wt2 += double(w) * double(t) * t;
6340 sum_wcos += double(w) * cosval;
6341 lastID = t_pulse[j][2];
6347 const bool is_miss = (t0 >= 0.98f * miss_distance);
6348 const float distance = (sum_w > 0.0) ?
float(sum_wt / sum_w) : t0;
6354 intensity = (count == Npulse) ? 1.0f : 0.0f;
6357 intensity = float(sum_wcos / denom);
6358 double var = (sum_w > 0.0) ? (sum_wt2 / sum_w -
double(distance) * distance) : 0.0;
6362 echo_width = sqrtf(range_resolution * range_resolution +
float(var));
6365 t_hit.push_back({distance, intensity, float(count), lastID, echo_width});
6370 if (detection_threshold > 0.f) {
6371 std::vector<std::vector<float>> kept;
6372 kept.reserve(t_hit.size());
6373 for (
const auto &h: t_hit) {
6374 const bool h_is_miss = (h[0] >= 0.98f * miss_distance);
6375 if (h_is_miss || fabsf(h[1]) >= detection_threshold) {
6390 if (max_returns > 0) {
6391 std::vector<size_t> real_idx;
6392 real_idx.reserve(t_hit.size());
6393 for (
size_t k = 0; k < t_hit.size(); k++) {
6394 const bool h_is_miss = (t_hit[k][0] >= 0.98f * miss_distance);
6396 real_idx.push_back(k);
6400 if (real_idx.empty()) {
6402 if (t_hit.size() > 1) {
6403 std::vector<std::vector<float>> only_miss{t_hit[0]};
6404 t_hit.swap(only_miss);
6411 const size_t last_idx = real_idx.back();
6412 size_t strongest_idx = real_idx.front();
6413 for (
const size_t k: real_idx) {
6414 if (fabsf(t_hit[k][1]) > fabsf(t_hit[strongest_idx][1])) {
6418 std::vector<size_t> ranked;
6419 if (strongest_idx == last_idx) {
6420 ranked = {last_idx};
6422 ranked = {strongest_idx, last_idx};
6424 std::sort(ranked.begin(), ranked.end(), [&](
size_t a,
size_t b) { return t_hit[a][0] < t_hit[b][0]; });
6426 std::vector<std::vector<float>> kept;
6427 kept.reserve(ranked.size());
6428 for (
const size_t k: ranked) {
6429 kept.push_back(t_hit[k]);
6433 std::vector<size_t> ranked = real_idx;
6434 if (
int(ranked.size()) > max_returns) {
6436 auto better = [&](
size_t a,
size_t b) {
6438 return fabsf(t_hit[a][1]) > fabsf(t_hit[b][1]);
6440 return t_hit[a][0] < t_hit[b][0];
6442 return t_hit[a][0] > t_hit[b][0];
6445 std::partial_sort(ranked.begin(), ranked.begin() + max_returns, ranked.end(), better);
6446 ranked.resize(
size_t(max_returns));
6448 std::sort(ranked.begin(), ranked.end(), [&](
size_t a,
size_t b) { return t_hit[a][0] < t_hit[b][0]; });
6451 std::vector<std::vector<float>> kept;
6452 kept.reserve(ranked.size());
6453 for (
const size_t k: ranked) {
6454 kept.push_back(t_hit[k]);
6489 syntheticScan(
context, rays_per_pulse, pulse_distance_threshold, scan_grid_only, record_misses,
true);
6496 std::vector<ReturnMode> saved_modes(scans.size());
6497 for (
size_t s = 0; s < scans.size(); s++) {
6498 saved_modes[s] = scans[s].returnMode;
6499 scans[s].returnMode = return_mode;
6501 syntheticScan(
context, rays_per_pulse, pulse_distance_threshold, scan_grid_only, record_misses, append);
6502 for (
size_t s = 0; s < scans.size(); s++) {
6503 scans[s].returnMode = saved_modes[s];
6514 for (
auto &hit_table: hit_tables) {
6515 hit_table.resize(hit_table.Ntheta, hit_table.Nphi, -1);
6517 hitgridcellcomputed =
false;
6518 triangulationcomputed =
false;
6522 if (rays_per_pulse < 1) {
6525 Npulse = rays_per_pulse;
6528 if (printmessages) {
6530 std::cout <<
"Performing multi-return synthetic LiDAR scan..." << std::endl;
6532 std::cout <<
"Performing single-return synthetic LiDAR scan..." << std::endl;
6537 std::cout <<
"WARNING (syntheticScan): No scans added to the point cloud. Exiting.." << std::endl;
6549 if (scan_grid_only ==
false) {
6553 context->getDomainBoundingBox(xbounds, ybounds, zbounds);
6554 bb_center =
helios::make_vec3(xbounds.
x + 0.5 * (xbounds.
y - xbounds.
x), ybounds.
x + 0.5 * (ybounds.
y - ybounds.
x), zbounds.
x + 0.5 * (zbounds.
y - zbounds.
x));
6562 const float bb_pad = 1e-4f;
6563 if (bb_size.
x < bb_pad)
6565 if (bb_size.
y < bb_pad)
6567 if (bb_size.
z < bb_pad)
6575 bb_center =
helios::make_vec3(boxmin.
x + 0.5 * (boxmax.
x - boxmin.
x), boxmin.
y + 0.5 * (boxmax.
y - boxmin.
y), boxmin.
z + 0.5 * (boxmax.
z - boxmin.
z));
6580 const float bb_pad = 1e-4f;
6581 if (bb_size.
x < bb_pad)
6583 if (bb_size.
y < bb_pad)
6585 if (bb_size.
z < bb_pad)
6593 std::map<std::string, int> textures;
6594 std::map<std::string, helios::int2> texture_size;
6595 std::map<std::string, std::vector<std::vector<bool>>> texture_data;
6598 std::vector<uint> UUIDs_all =
context->getAllUUIDs();
6600 std::vector<uint> ID_mapping;
6605 size_t Npatches = 0;
6606 for (
int p = 0; p < UUIDs_all.size(); p++) {
6612 ID_mapping.resize(Npatches);
6615 int *patch_textureID = (
int *) malloc(Npatches *
sizeof(
int));
6619 for (
int p = 0; p < UUIDs_all.size(); p++) {
6620 uint UUID = UUIDs_all.at(p);
6622 std::vector<helios::vec3> verts =
context->getPrimitiveVertices(UUID);
6623 patch_vertex[4 * c] = verts.at(0);
6624 patch_vertex[4 * c + 1] = verts.at(1);
6625 patch_vertex[4 * c + 2] = verts.at(2);
6626 patch_vertex[4 * c + 3] = verts.at(3);
6628 ID_mapping.at(c) = UUIDs_all.at(p);
6630 if (!
context->getPrimitiveTextureFile(UUID).empty() &&
context->primitiveTextureHasTransparencyChannel(UUID)) {
6631 std::string tex =
context->getPrimitiveTextureFile(UUID);
6632 std::map<std::string, int>::iterator it = textures.find(tex);
6633 if (it != textures.end()) {
6634 patch_textureID[c] = textures.at(tex);
6636 patch_textureID[c] = tID;
6637 textures[tex] = tID;
6640 texture_data[tex] = *
context->getPrimitiveTextureTransparencyData(UUID);
6644 std::vector<helios::vec2> uv =
context->getPrimitiveTextureUV(UUID);
6645 if (uv.size() == 4) {
6646 patch_uv[2 * c] = uv.at(1);
6647 patch_uv[2 * c + 1] = uv.at(3);
6654 patch_textureID[c] = -1;
6666 size_t Ntriangles = 0;
6667 for (
int p = 0; p < UUIDs_all.size(); p++) {
6673 ID_mapping.resize(Npatches + Ntriangles);
6676 int *tri_textureID = (
int *) malloc(Ntriangles *
sizeof(
int));
6680 for (
int p = 0; p < UUIDs_all.size(); p++) {
6681 uint UUID = UUIDs_all.at(p);
6683 std::vector<helios::vec3> verts =
context->getPrimitiveVertices(UUID);
6684 tri_vertex[3 * c] = verts.at(0);
6685 tri_vertex[3 * c + 1] = verts.at(1);
6686 tri_vertex[3 * c + 2] = verts.at(2);
6688 ID_mapping.at(Npatches + c) = UUIDs_all.at(p);
6690 if (!
context->getPrimitiveTextureFile(UUID).empty() &&
context->primitiveTextureHasTransparencyChannel(UUID)) {
6691 std::string tex =
context->getPrimitiveTextureFile(UUID);
6692 std::map<std::string, int>::iterator it = textures.find(tex);
6693 if (it != textures.end()) {
6694 tri_textureID[c] = textures.at(tex);
6696 tri_textureID[c] = tID;
6697 textures[tex] = tID;
6700 texture_data[tex] = *
context->getPrimitiveTextureTransparencyData(UUID);
6704 std::vector<helios::vec2> uv =
context->getPrimitiveTextureUV(UUID);
6705 assert(uv.size() == 3);
6706 tri_uv[3 * c] = uv.at(0);
6707 tri_uv[3 * c + 1] = uv.at(1);
6708 tri_uv[3 * c + 2] = uv.at(2);
6711 tri_textureID[c] = -1;
6721 const int Ntextures = textures.size();
6724 for (std::map<std::string, helios::int2>::iterator it = texture_size.begin(); it != texture_size.end(); ++it) {
6725 if (it->second.x > masksize_max.
x) {
6726 masksize_max.
x = it->second.x;
6728 if (it->second.y > masksize_max.
y) {
6729 masksize_max.
y = it->second.y;
6733 bool *maskdata = (
bool *) malloc(Ntextures * masksize_max.
x * masksize_max.
y *
sizeof(
bool));
6736 for (std::map<std::string, helios::int2>::iterator it = texture_size.begin(); it != texture_size.end(); ++it) {
6737 std::string texture_file = it->first;
6739 int ID = textures.at(texture_file);
6741 masksize[ID] = it->second;
6744 for (
int j = 0; j < masksize_max.
y; j++) {
6745 for (
int i = 0; i < masksize_max.
x; i++) {
6747 if (i < texture_size.at(texture_file).x && j < texture_size.at(texture_file).y) {
6748 maskdata[ID * masksize_max.
x * masksize_max.
y + ind] = texture_data.at(texture_file).at(j).at(i);
6750 maskdata[ID * masksize_max.
x * masksize_max.
y + ind] =
false;
6761 struct TextureColorMap {
6764 std::vector<helios::RGBcolor> pixels;
6766 std::map<std::string, TextureColorMap> texture_color_cache;
6768 auto load_texture_colors = [&](
const std::string &filename) ->
const TextureColorMap & {
6769 auto it = texture_color_cache.find(filename);
6770 if (it != texture_color_cache.end()) {
6773 TextureColorMap entry;
6775 size_t dot = filename.find_last_of(
'.');
6776 if (dot != std::string::npos) {
6777 ext = filename.substr(dot);
6778 for (
char &ch: ext) {
6779 ch =
static_cast<char>(std::tolower(
static_cast<unsigned char>(ch)));
6782 if (ext ==
".png") {
6783 std::vector<helios::RGBAcolor> rgba;
6785 entry.pixels.resize(rgba.size());
6786 for (
size_t i = 0; i < rgba.size(); i++) {
6789 }
else if (ext ==
".jpg" || ext ==
".jpeg") {
6792 return texture_color_cache.emplace(filename, std::move(entry)).first->second;
6796 const std::string tex_file =
context->getPrimitiveTextureFile(UUID);
6797 if (tex_file.empty() ||
context->isPrimitiveTextureColorOverridden(UUID)) {
6798 return context->getPrimitiveColor(UUID);
6800 const TextureColorMap &tex = load_texture_colors(tex_file);
6801 if (tex.pixels.empty()) {
6802 return context->getPrimitiveColor(UUID);
6805 std::vector<helios::vec3> verts =
context->getPrimitiveVertices(UUID);
6806 std::vector<helios::vec2> uvs =
context->getPrimitiveTextureUV(UUID);
6815 float e1_sq = e1 * e1;
6816 float e2_sq = e2 * e2;
6817 float s_param = (e1_sq > 0.f) ? (d * e1) / e1_sq : 0.f;
6818 float t_param = (e2_sq > 0.f) ? (d * e2) / e2_sq : 0.f;
6821 else if (s_param > 1.f)
6825 else if (t_param > 1.f)
6827 if (uvs.size() == 4) {
6828 uv = (1.f - s_param) * (1.f - t_param) * uvs[0] + s_param * (1.f - t_param) * uvs[1] + s_param * t_param * uvs[2] + (1.f - s_param) * t_param * uvs[3];
6836 float dot11 = e1 * e1;
6837 float dot12 = e1 * e2;
6838 float dot22 = e2 * e2;
6839 float dot1d = e1 * d;
6840 float dot2d = e2 * d;
6841 float denom = dot11 * dot22 - dot12 * dot12;
6842 if (std::fabs(denom) < 1e-20f) {
6843 return context->getPrimitiveColor(UUID);
6845 float inv_denom = 1.f / denom;
6846 float beta = (dot22 * dot1d - dot12 * dot2d) * inv_denom;
6847 float gamma = (dot11 * dot2d - dot12 * dot1d) * inv_denom;
6848 uv = uvs[0] + beta * (uvs[1] - uvs[0]) + gamma * (uvs[2] - uvs[0]);
6850 return context->getPrimitiveColor(UUID);
6854 uv.
x -= std::floor(uv.
x);
6855 uv.
y -= std::floor(uv.
y);
6857 int px =
static_cast<int>(uv.
x *
static_cast<float>(tex.width));
6860 if (px >=
static_cast<int>(tex.width))
6861 px =
static_cast<int>(tex.width) - 1;
6863 int py =
static_cast<int>((1.f - uv.
y) *
static_cast<float>(tex.height));
6866 if (py >=
static_cast<int>(tex.height))
6867 py =
static_cast<int>(tex.height) - 1;
6868 return tex.pixels[
static_cast<size_t>(py) * tex.width +
static_cast<size_t>(px)];
6875 if (progress_callback) {
6879 if (synthetic_scan_progress !=
nullptr) {
6880 *synthetic_scan_progress = 0;
6887 if (synthetic_scan_progress !=
nullptr) {
6888 *synthetic_scan_progress = s;
6894 progress_bar.
update(
static_cast<size_t>(s));
6902 float thetamin = thetarange.
x;
6903 float thetamax = thetarange.
y;
6905 float phimin = phirange.
x;
6906 float phimax = phirange.
y;
6928 bool apply_azimuth = (scanAzimuthOffset != 0.f);
6929 bool apply_tilt = (scanTiltRoll != 0.f || scanTiltPitch != 0.f);
6938 const float heading = phimin - scanAzimuthOffset;
6957 const bool is_moving = scan.
isMoving;
6959 const double pulse_t0 = scan.
t0;
6963 std::vector<helios::vec3> raydir;
6964 raydir.resize(Ntheta * Nphi);
6968 std::vector<helios::vec3> raygrid_origin;
6969 raygrid_origin.resize(Ntheta * Nphi, scan_origin);
6979 const float dphi = spinning_multibeam ? (phimax - phimin) /
float(Nphi) : ((Nphi > 1) ? (phimax - phimin) / float(Nphi - 1) : 0.f);
6980 const float dtheta = (Ntheta > 1) ? (thetamax - thetamin) / float(Ntheta - 1) : 0.f;
6991 const float dphi_per_row = (!spinning_multibeam && !risley && Ntheta > 0) ? dphi /
float(Ntheta) : 0.f;
6993 for (
uint j = 0; j < Nphi; j++) {
6994 float phi = phimin + float(j) * dphi;
6995 for (
uint i = 0; i < Ntheta; i++) {
7000 dir = risleyBodyDirection(scan, j);
7002 float theta_z = spinning_multibeam ? scan.
beamZenithAngles.at(i) : (thetamin + float(i) * dtheta);
7003 float theta_elev = 0.5f *
M_PI - theta_z;
7005 float phi_skew = phi + float(i) * dphi_per_row;
7011 const size_t ordinal = size_t(Ntheta) * j + i;
7012 const double t = pulse_t0 + double(ordinal) * pulse_period;
7015 scan.
poseAt(t, pos, quat);
7016 helios::vec3 dir_body = quat_rotate(boresight_quat, dir);
7017 dir = quat_rotate(quat, dir_body);
7019 raygrid_origin.at(Ntheta * j + i) = pos + quat_rotate(quat, scan.
lever_arm);
7021 if (apply_azimuth) {
7029 raydir.at(Ntheta * j + i) = dir;
7033 size_t N = Ntheta * Nphi;
7036 std::vector<uint> bb_hit(N, 0);
7044 bool origin_inside_bb = !is_moving && (scan_origin.
x >= bb_min.
x && scan_origin.
x <= bb_max.
x && scan_origin.
y >= bb_min.
y && scan_origin.
y <= bb_max.
y && scan_origin.
z >= bb_min.
z && scan_origin.
z <= bb_max.
z);
7046 for (
size_t r = 0; r < N; r++) {
7048 if (origin_inside_bb) {
7057 if (is_moving && cell_origin.
x >= bb_min.
x && cell_origin.
x <= bb_max.
x && cell_origin.
y >= bb_min.
y && cell_origin.
y <= bb_max.
y && cell_origin.
z >= bb_min.
z && cell_origin.
z <= bb_max.
z) {
7065 float tx_min, tx_max, ty_min, ty_max, tz_min, tz_max;
7067 float a = 1.0f / ray_dir.
x;
7069 tx_min = (bb_min.
x - cell_origin.
x) * a;
7070 tx_max = (bb_max.
x - cell_origin.
x) * a;
7072 tx_min = (bb_max.
x - cell_origin.
x) * a;
7073 tx_max = (bb_min.
x - cell_origin.
x) * a;
7076 float b = 1.0f / ray_dir.
y;
7078 ty_min = (bb_min.
y - cell_origin.
y) * b;
7079 ty_max = (bb_max.
y - cell_origin.
y) * b;
7081 ty_min = (bb_max.
y - cell_origin.
y) * b;
7082 ty_max = (bb_min.
y - cell_origin.
y) * b;
7085 float c = 1.0f / ray_dir.
z;
7087 tz_min = (bb_min.
z - cell_origin.
z) * c;
7088 tz_max = (bb_max.
z - cell_origin.
z) * c;
7090 tz_min = (bb_max.
z - cell_origin.
z) * c;
7091 tz_max = (bb_min.
z - cell_origin.
z) * c;
7108 if (t0 < t1 && t1 > 1e-6f) {
7114 size_t total_scan_rays = Ntheta * Nphi;
7117 for (
int i = 0; i < total_scan_rays; i++) {
7118 if (bb_hit[i] == 1) {
7121 hit_out += sin(dir.
zenith);
7126 std::vector<helios::vec3> base_directions;
7127 base_directions.reserve(N);
7128 std::vector<helios::int2> pulse_scangrid_ij(N);
7130 std::vector<helios::vec3> pulse_origin(N);
7133 for (
int i = 0; i < Ntheta * Nphi; i++) {
7134 if (bb_hit[i] == 1) {
7136 base_directions.push_back(raydir.at(i));
7138 int jj = floor(i / Ntheta);
7139 int ii = i - jj * Ntheta;
7141 pulse_origin[count] = raygrid_origin.at(i);
7152 if (record_misses) {
7159 for (
int i = 0; i < Ntheta * Nphi; i++) {
7160 std::map<std::string, double> data;
7161 data[
"target_index"] = 0;
7162 data[
"target_count"] = 1;
7163 data[
"deviation"] = 0.0;
7166 data[
"timestamp"] = pulse_t0 + double(i) * pulse_period;
7167 data[
"intensity"] = 1.0;
7168 data[
"distance"] = miss_dist;
7169 data[
"nRaysHit"] = Npulse;
7170 data[
"echo_width"] = 0.0;
7171 data[
"is_miss"] = 1.0;
7172 if (spinning_multibeam) {
7173 data[
"channel"] = double(i % Ntheta);
7175 if (std::find(column_format.begin(), column_format.end(),
"reflectance") != column_format.end()) {
7176 data[
"reflectance"] = 0.0;
7182 data[
"pulse_id"] = double(i);
7183 data[
"origin_x"] = cell_origin.
x;
7184 data[
"origin_y"] = cell_origin.
y;
7185 data[
"origin_z"] = cell_origin.
z;
7191 scan_warnings.
addWarning(
"synthetic_rays_no_hit",
"Synthetic rays did not hit any primitives.");
7204 const ReturnMode return_mode = scans.at(s).returnMode;
7208 const int max_returns = (return_mode ==
RETURN_MODE_MULTI) ? 0 : scans.at(s).maxReturns;
7209 const float detection_threshold = scans.at(s).detectionThreshold;
7210 const float range_resolution = (scans.at(s).pulseWidth > 0.f) ? scans.at(s).pulseWidth : pulse_distance_threshold;
7217 std::vector<helios::vec3> nominal_directions(N);
7218 for (
size_t beam = 0; beam < N; beam++) {
7220 if (angle_noise_stddev > 0) {
7227 float a =
context->randn(0.f, angle_noise_stddev);
7228 float b =
context->randn(0.f, angle_noise_stddev);
7231 nominal_directions[beam] = jittered;
7233 nominal_directions[beam] = base_dir;
7243 constexpr size_t BYTES_PER_SUBRAY = 40;
7244 constexpr size_t MIN_RAYS_PER_CHUNK = 1050000;
7250 size_t effective_budget_bytes = synthetic_scan_memory_budget_bytes;
7251 if (effective_budget_bytes == 0) {
7252 effective_budget_bytes = collision_detection->
isGPUAccelerationEnabled() ? SYNTHETIC_SCAN_DEFAULT_BUDGET_GPU : SYNTHETIC_SCAN_DEFAULT_BUDGET_CPU;
7254 size_t target_subrays = effective_budget_bytes / BYTES_PER_SUBRAY;
7255 size_t chunk_beams = target_subrays / size_t(Npulse);
7256 size_t min_chunk_beams = (MIN_RAYS_PER_CHUNK + size_t(Npulse) - 1) /
size_t(Npulse);
7257 if (chunk_beams < min_chunk_beams) {
7258 chunk_beams = min_chunk_beams;
7260 if (chunk_beams < 1) {
7263 if (chunk_beams > N) {
7269 size_t beams_with_zero_hits = 0;
7270 size_t beams_with_one_hit = 0;
7271 size_t beams_with_multi_hits = 0;
7275 prepareUnifiedRayTracing(
context);
7279 const size_t chunk_capacity = chunk_beams * size_t(Npulse);
7282 float *hit_t = (
float *) malloc(chunk_capacity *
sizeof(
float));
7283 float *hit_fnorm = (
float *) malloc(chunk_capacity *
sizeof(
float));
7284 int *hit_ID = (
int *) malloc(chunk_capacity *
sizeof(
int));
7285 std::vector<float> subray_weight(chunk_capacity, 1.0f);
7286 if (direction ==
nullptr || ray_origins ==
nullptr || hit_t ==
nullptr || hit_fnorm ==
nullptr || hit_ID ==
nullptr) {
7287 helios_runtime_error(
"ERROR (LiDARcloud::syntheticScan): failed to allocate ray-tracing scratch buffers for a beam chunk of " + std::to_string(chunk_capacity) +
" sub-rays. Lower the synthetic-scan memory budget (setSyntheticScanMemoryBudget) or reduce rays_per_pulse.");
7295 std::set<std::string> warm_tex_files;
7296 std::vector<uint> all_uuids =
context->getAllUUIDs();
7297 for (
uint warm_uuid: all_uuids) {
7298 std::string tf =
context->getPrimitiveTextureFile(warm_uuid);
7300 warm_tex_files.insert(tf);
7303 for (
const std::string &tf: warm_tex_files) {
7304 load_texture_colors(tf);
7311 const unsigned int range_noise_seed_salt = (range_noise_stddev > 0.f) ?
static_cast<unsigned int>(
context->randu() * 4294967000.0f) : 0u;
7314 struct SyntheticBeamHit {
7317 std::map<std::string, double> data;
7319 struct SyntheticBeamOutput {
7320 std::vector<SyntheticBeamHit> hits;
7322 int npulse_hits = 0;
7328 constexpr float GOLDEN_ANGLE = 2.39996322972865332f;
7332 auto radicalInverse2 = [](uint32_t i) ->
float {
7333 i = (i << 16) | (i >> 16);
7334 i = ((i & 0x55555555u) << 1) | ((i & 0xAAAAAAAAu) >> 1);
7335 i = ((i & 0x33333333u) << 2) | ((i & 0xCCCCCCCCu) >> 2);
7336 i = ((i & 0x0F0F0F0Fu) << 4) | ((i & 0xF0F0F0F0u) >> 4);
7337 i = ((i & 0x00FF00FFu) << 8) | ((i & 0xFF00FF00u) >> 8);
7338 return float(i) * 2.3283064365386963e-10f;
7341 for (
size_t chunk_begin = 0; chunk_begin < N; chunk_begin += chunk_beams) {
7347 if (cancel_flag !=
nullptr && *cancel_flag != 0) {
7351 const size_t chunk_end = std::min(chunk_begin + chunk_beams, N);
7352 const size_t chunk_N = chunk_end - chunk_begin;
7376 const bool draw_divergence = (beam_divergence != 0.0f && Npulse > 1);
7377 const bool draw_aperture = (exit_diameter > 0.0f);
7383 std::vector<float> divergence_rand, aperture_rand;
7384 if (draw_divergence) {
7385 divergence_rand.resize(chunk_N * 2, 0.f);
7386 for (
size_t local = 0; local < chunk_N; local++) {
7387 divergence_rand[local * 2] =
context->randu();
7388 divergence_rand[local * 2 + 1] = 2.0f * float(
M_PI) *
context->randu();
7391 if (draw_aperture) {
7392 aperture_rand.resize(chunk_N * 2, 0.f);
7393 for (
size_t local = 0; local < chunk_N; local++) {
7394 aperture_rand[local * 2] =
context->randu();
7395 aperture_rand[local * 2 + 1] = 2.0f * float(
M_PI) *
context->randu();
7407 const float footprint_trim = fminf(0.5f, (detection_threshold > 0.f) ? detection_threshold : 1.0e-3f);
7408 const float u_max_foot = 1.0f - footprint_trim;
7410 const float aperture_radius = 0.5f * exit_diameter;
7411#pragma omp parallel for schedule(dynamic, 256)
7412 for (
int local = 0; local < static_cast<int>(chunk_N); local++) {
7413 const size_t global_r = chunk_begin + local;
7414 const helios::vec3 base_dir = nominal_directions[global_r];
7415 const helios::vec3 beam_origin = pulse_origin[global_r];
7419 if (draw_aperture) {
7428 for (
int p = 0; p < Npulse; p++) {
7429 const size_t idx = local * size_t(Npulse) + size_t(p);
7435 if (p == 0 || beam_divergence == 0.0f) {
7436 direction[idx] = base_dir;
7438 const float xi_r = divergence_rand[local * 2];
7439 const float phi0 = divergence_rand[local * 2 + 1];
7440 const int j = p - 1;
7441 const int M = Npulse - 1;
7442 float uu = (float(j) + 0.5f) /
float(M) + xi_r;
7448 float theta_offset = beam_divergence * sqrtf(-0.5f * logf(1.0f - uu));
7449 const float phi_offset = phi0 + float(j) * GOLDEN_ANGLE;
7451 const float new_elevation = base_spherical.
elevation + theta_offset * cosf(phi_offset);
7452 const float new_azimuth = base_spherical.
azimuth + theta_offset * sinf(phi_offset) / fmaxf(cosf(base_spherical.
elevation), 1e-6f);
7455 direction[idx] = perturbed_dir;
7462 if (draw_aperture) {
7464 ray_origins[idx] = beam_origin;
7466 const float xi_ap = aperture_rand[local * 2];
7467 const float phi0_ap = aperture_rand[local * 2 + 1];
7468 const int j = p - 1;
7469 float s = radicalInverse2((uint32_t)(j + 1)) + xi_ap;
7473 float r_sample = aperture_radius * sqrtf(-0.5f * logf(1.0f - s));
7474 const float theta = phi0_ap + float(j) * GOLDEN_ANGLE;
7475 const float x_disk = r_sample * cosf(theta);
7476 const float y_disk = r_sample * sinf(theta);
7478 ray_origins[idx] = beam_origin + offset;
7481 ray_origins[idx] = beam_origin;
7484 subray_weight[idx] = w;
7491 castRaysUnified(chunk_N *
size_t(Npulse), ray_origins, direction, hit_t, hit_fnorm, hit_ID,
size_t(Npulse));
7497 std::vector<SyntheticBeamOutput> beam_outputs(chunk_N);
7498#pragma omp parallel for schedule(dynamic, 256)
7499 for (
int local = 0; local < static_cast<int>(chunk_N); local++) {
7500 const size_t r = local;
7501 const size_t global_r = chunk_begin + local;
7502 SyntheticBeamOutput &beam_out = beam_outputs[local];
7508 std::vector<std::vector<float>> t_pulse;
7509 float total_pulse_weight = 0.f;
7512 for (
size_t p = 0; p < Npulse; p++) {
7514 float t = hit_t[r * Npulse + p];
7515 float i = hit_fnorm[r * Npulse + p];
7516 float ID = float(hit_ID[r * Npulse + p]);
7517 float w = subray_weight[r * Npulse + p];
7519 total_pulse_weight += w;
7521 if (record_misses || (!record_misses && t < miss_distance)) {
7522 std::vector<float> v{t, i, ID, w};
7523 t_pulse.push_back(v);
7528 beam_out.npulse_hits = int(t_pulse.size());
7532 std::vector<std::vector<float>> t_hit = detectReturnsFromSubrays(t_pulse, total_pulse_weight, Npulse, range_resolution, detection_threshold, max_returns, single_return_selection, miss_distance);
7535 int non_miss_count = 0;
7536 for (
size_t hit = 0; hit < t_hit.size(); hit++) {
7537 if (t_hit.at(hit).at(0) < 0.98f * miss_distance) {
7545 std::mt19937 beam_rng;
7546 std::normal_distribution<float> range_noise_dist;
7547 if (range_noise_stddev > 0.f) {
7548 beam_rng.seed(range_noise_seed_salt +
static_cast<unsigned int>(global_r));
7549 range_noise_dist = std::normal_distribution<float>(0.f, range_noise_stddev);
7552 int real_hit_index = 0;
7553 for (
size_t hit = 0; hit < t_hit.size(); hit++) {
7555 std::map<std::string, double> data;
7558 bool is_miss = (t_hit.at(hit).at(0) >= 0.98f * miss_distance);
7564 float measured_distance = t_hit.at(hit).at(0);
7565 if (!is_miss && range_noise_stddev > 0.f) {
7566 measured_distance += range_noise_dist(beam_rng);
7571 data[
"target_index"] = 99;
7573 data[
"target_index"] = real_hit_index;
7577 data[
"is_miss"] = is_miss ? 1.0 : 0.0;
7578 data[
"target_count"] = t_hit.size();
7585 const float echo_width = t_hit.at(hit).at(4);
7586 float deviation = 0.f;
7587 if (!is_miss && range_resolution > 0.f) {
7588 const float excess_var = echo_width * echo_width - range_resolution * range_resolution;
7589 deviation = (excess_var > 0.f) ? sqrtf(excess_var) / range_resolution : 0.f;
7591 data[
"deviation"] = deviation;
7596 const size_t pulse_ordinal = size_t(pulse_scangrid_ij.at(global_r).y) * Ntheta + size_t(pulse_scangrid_ij.at(global_r).x);
7597 data[
"timestamp"] = pulse_t0 + double(pulse_ordinal) * pulse_period;
7601 data[
"distance"] = measured_distance;
7602 data[
"nRaysHit"] = t_hit.at(hit).at(2);
7605 data[
"echo_width"] = t_hit.at(hit).at(4);
7606 if (spinning_multibeam) {
7607 data[
"channel"] = double(pulse_scangrid_ij.at(global_r).x);
7610 float UUID = t_hit.at(hit).at(3);
7616 const helios::vec3 beam_origin = pulse_origin[global_r];
7617 helios::vec3 p = beam_origin + dir * measured_distance;
7619 data[
"pulse_id"] = double(pulse_ordinal);
7620 data[
"origin_x"] = beam_origin.
x;
7621 data[
"origin_y"] = beam_origin.
y;
7622 data[
"origin_z"] = beam_origin.
z;
7627 if (UUID >= 0 &&
context->doesPrimitiveExist(
uint(UUID))) {
7629 color = sample_hit_color(
uint(UUID), p);
7636 for (
const std::string &label: column_format) {
7637 if (isStandardColumnToken(label)) {
7642 if (!resolveScalarHitData(
context,
uint(UUID), label, value)) {
7646 if (label ==
"reflectivity_lidar") {
7648 data.at(
"intensity") *= value;
7649 }
else if (label ==
"reflectance") {
7655 data[label] = value;
7670 if (std::find(column_format.begin(), column_format.end(),
"reflectance") != column_format.end()) {
7671 constexpr double REFLECTANCE_FLOOR_DB = -999.0;
7672 double abs_intensity = fabs(data.at(
"intensity"));
7673 data[
"reflectance"] = (abs_intensity > 0.0) ? 10.0 * log10(abs_intensity) : REFLECTANCE_FLOOR_DB;
7678 beam_out.hits.push_back(SyntheticBeamHit{p, color, std::move(data)});
7695 std::vector<std::set<std::string>> thread_label_sets(
static_cast<size_t>(omp_get_max_threads()));
7696#pragma omp parallel for schedule(dynamic, 256)
7697 for (
int local = 0; local < static_cast<int>(chunk_N); local++) {
7698 std::set<std::string> &my_labels = thread_label_sets[
static_cast<size_t>(omp_get_thread_num())];
7699 for (
const SyntheticBeamHit &bh: beam_outputs[local].hits) {
7700 for (
const auto &kv: bh.data) {
7701 my_labels.insert(kv.first);
7705 std::set<std::string> all_labels;
7706 for (
const std::set<std::string> &ts: thread_label_sets) {
7707 all_labels.insert(ts.begin(), ts.end());
7711 std::unordered_map<std::string, size_t> slot_of;
7712 for (
const std::string &lbl: all_labels) {
7713 slot_of[lbl] = getOrCreateHitDataColumn(lbl);
7717 std::vector<size_t> row_offset(chunk_N + 1, 0);
7718 for (
size_t local = 0; local < chunk_N; local++) {
7719 const SyntheticBeamOutput &beam_out = beam_outputs[local];
7720 if (beam_out.npulse_hits == 0) {
7721 beams_with_zero_hits++;
7722 }
else if (beam_out.npulse_hits == 1) {
7723 beams_with_one_hit++;
7725 beams_with_multi_hits++;
7727 row_offset[local + 1] = row_offset[local] + beam_out.hits.size();
7729 const size_t chunk_hits = row_offset[chunk_N];
7730 const size_t old_n = hits.size();
7731 const size_t new_n = old_n + chunk_hits;
7734 for (
size_t sl = 0; sl < hit_data_columns.size(); sl++) {
7735 hit_data_columns[sl].resize(new_n, 0.0);
7736 hit_data_present[sl].resize(new_n,
char(0));
7741#pragma omp parallel for schedule(dynamic, 256)
7742 for (
int local = 0; local < static_cast<int>(chunk_N); local++) {
7743 SyntheticBeamOutput &beam_out = beam_outputs[local];
7744 if (beam_out.hits.empty()) {
7748 const size_t base = old_n + row_offset[local];
7749 for (
size_t h = 0; h < beam_out.hits.size(); h++) {
7750 SyntheticBeamHit &bh = beam_out.hits[h];
7751 const size_t row = base + h;
7753 hp.position = bh.xyz;
7754 hp.direction = beam_out.dir_sph;
7756 hp.color = bh.color;
7758 for (
const auto &kv: bh.data) {
7759 const size_t sl = slot_of.at(kv.first);
7760 hit_data_columns[sl][row] = kv.second;
7761 hit_data_present[sl][row] = char(1);
7765 Nhits += chunk_hits;
7770 finishUnifiedRayTracing();
7779 if (printmessages) {
7780 std::cout <<
"Created synthetic scan #" << s <<
" with " << Nhits <<
" hit points." << std::endl;
7785 if (cancel_flag !=
nullptr && *cancel_flag != 0) {
7791 if (synthetic_scan_progress !=
nullptr) {
7797 scan_warnings.
report(std::cerr);
7801 free(patch_textureID);
7804 free(tri_textureID);