1.3.64
 
Loading...
Searching...
No Matches
LiDAR.h
Go to the documentation of this file.
1
16#ifndef LIDARPLUGIN
17#define LIDARPLUGIN
18
19#include "CollisionDetection.h"
20#include "Context.h"
21#include "Visualizer.h"
22
23#include "s_hull_pro.h"
24
25template<class datatype>
26class HitTable {
27public:
28 uint Ntheta, Nphi;
29
30 HitTable(void) {
31 Ntheta = 0;
32 Nphi = 0;
33 }
34 HitTable(const int nx, const int ny) {
35 Ntheta = nx;
36 Nphi = ny;
37 data.resize(Nphi);
38 for (int j = 0; j < Nphi; j++) {
39 data.at(j).resize(Ntheta);
40 }
41 }
42 HitTable(const int nx, const int ny, const datatype initval) {
43 Ntheta = nx;
44 Nphi = ny;
45 data.resize(Nphi);
46 for (int j = 0; j < Nphi; j++) {
47 data.at(j).resize(Ntheta, initval);
48 }
49 }
50
51 datatype get(const int i, const int j) const {
52 if (i >= 0 && i < Ntheta && j >= 0 && j < Nphi) {
53 return data.at(j).at(i);
54 } else {
55 std::cerr << "ERROR (hit_map.get): get index out of range. Attempting to get index map at (" << i << "," << j << "), but size of scan is " << Ntheta << " x " << Nphi << "." << std::endl;
56 exit(EXIT_FAILURE);
57 }
58 }
59 void set(const int i, const int j, const datatype value) {
60 if (i >= 0 && i < Ntheta && j >= 0 && j < Nphi) {
61 data.at(j).at(i) = value;
62 } else {
63 std::cerr << "ERROR (hit_map.set): set index out of range. Attempting to set index map at (" << i << "," << j << "), but size of scan is " << Ntheta << " x " << Nphi << "." << std::endl;
64 exit(EXIT_FAILURE);
65 }
66 }
67 void resize(const int nx, const int ny, const datatype initval) {
68 Ntheta = nx;
69 Nphi = ny;
70 data.resize(Nphi);
71 for (int j = 0; j < Nphi; j++) {
72 data.at(j).resize(Ntheta);
73 for (int i = 0; i < Ntheta; i++) {
74 data.at(j).at(i) = initval;
75 }
76 }
77 }
78
79private:
80 std::vector<std::vector<datatype>> data;
81};
82
83struct HitPoint {
84 helios::vec3 position;
85 helios::SphericalCoord direction;
86 helios::int2 row_column;
87 helios::RGBcolor color;
88 std::map<std::string, double> data;
89 int gridcell;
90 int scanID;
91 HitPoint(void) {
92 position = helios::make_vec3(0, 0, 0);
93 direction = helios::make_SphericalCoord(0, 0);
94 row_column = helios::make_int2(0, 0);
95 color = helios::RGB::red;
96 gridcell = -2;
97 scanID = -1;
98 }
99 HitPoint(int __scanID, helios::vec3 __position, helios::SphericalCoord __direction, helios::int2 __row_column, helios::RGBcolor __color, std::map<std::string, double> __data) {
100 scanID = __scanID;
101 position = __position;
102 direction = __direction;
103 row_column = __row_column;
104 color = __color;
105 data = __data;
106 }
107};
108
109struct Triangulation {
110 helios::vec3 vertex0, vertex1, vertex2;
111 int ID0, ID1, ID2;
112 int scanID;
113 int gridcell;
114 helios::RGBcolor color;
115 float area;
116 Triangulation(void) {
117 vertex0 = helios::make_vec3(0, 0, 0);
118 vertex1 = helios::make_vec3(0, 0, 0);
119 vertex2 = helios::make_vec3(0, 0, 0);
120 ID0 = 0;
121 ID1 = 0;
122 ID2 = 0;
123 scanID = -1;
124 gridcell = -2;
125 color = helios::RGB::green;
126 area = 0;
127 }
128 Triangulation(int __scanID, helios::vec3 __vertex0, helios::vec3 __vertex1, helios::vec3 __vertex2, int __ID0, int __ID1, int __ID2, helios::RGBcolor __color, int __gridcell) {
129 scanID = __scanID;
130 vertex0 = __vertex0;
131 vertex1 = __vertex1;
132 vertex2 = __vertex2;
133 ID0 = __ID0;
134 ID1 = __ID1;
135 ID2 = __ID2;
136 gridcell = __gridcell;
137 color = __color;
138
139 // calculate area
140 helios::vec3 s0 = vertex1 - vertex0;
141 helios::vec3 s1 = vertex2 - vertex0;
142 helios::vec3 s2 = vertex2 - vertex1;
143
144 float a = s0.magnitude();
145 float b = s1.magnitude();
146 float c = s2.magnitude();
147 float s = 0.5f * (a + b + c);
148
149 area = sqrt(s * (s - a) * (s - b) * (s - c));
150 }
151};
152
153struct GridCell {
154 helios::vec3 center;
155 helios::vec3 global_anchor;
156 helios::vec3 size;
157 helios::vec3 global_size;
158 helios::int3 global_ijk;
159 helios::int3 global_count;
160 float azimuthal_rotation;
161 float leaf_area;
162 float Gtheta;
163 float ground_height;
164 float vegetation_height;
165 float maximum_height;
166 GridCell(helios::vec3 __center, helios::vec3 __global_anchor, helios::vec3 __size, helios::vec3 __global_size, float __azimuthal_rotation, helios::int3 __global_ijk, helios::int3 __global_count) {
167 center = __center;
168 global_anchor = __global_anchor;
169 size = __size;
170 global_size = __global_size;
171 azimuthal_rotation = __azimuthal_rotation;
172 global_ijk = __global_ijk;
173 global_count = __global_count;
174 leaf_area = 0;
175 Gtheta = 0;
176 ground_height = 0;
177 vegetation_height = 0;
178 maximum_height = 0;
179 }
180};
181
183
187struct ScanMetadata {
188
190 ScanMetadata();
191
193
202 ScanMetadata(const helios::vec3 &origin, uint Ntheta, float thetaMin, float thetaMax, uint Nphi, float phiMin, float phiMax, float exitDiameter, float beamDivergence, const std::vector<std::string> &columnFormat);
203
205 std::string data_file;
206
208 uint Ntheta;
210
213 float thetaMin;
215
218 float thetaMax;
219
221 uint Nphi;
223
226 float phiMin;
228
231 float phiMax;
232
234 helios::vec3 origin;
235
237
240 float exitDiameter;
241
243
246 float beamDivergence;
247
249 std::vector<std::string> columnFormat;
250
252
257 helios::SphericalCoord rc2direction(uint row, uint column) const;
258
260
264 helios::int2 direction2rc(const helios::SphericalCoord &direction) const;
265};
266
268class LiDARcloud {
269private:
270 size_t Nhits;
271
272 std::vector<ScanMetadata> scans;
273
274 std::vector<HitPoint> hits;
275
276 std::vector<GridCell> grid_cells;
277
278 std::vector<Triangulation> triangles;
279
281 std::vector<HitTable<int>> hit_tables;
282
284 bool hitgridcellcomputed;
285
287 bool triangulationcomputed;
288
290 bool printmessages;
291
293 CollisionDetection *collision_detection;
294
295 // -------- I/O --------- //
296
297 // -------- RECONSTRUCTION --------- //
298
299 // first index: leaf group, second index: triangle #
300 std::vector<std::vector<Triangulation>> reconstructed_triangles;
301
302 std::vector<std::vector<Triangulation>> reconstructed_trunk_triangles;
303
304 std::vector<helios::vec3> reconstructed_alphamasks_center;
305 std::vector<helios::vec2> reconstructed_alphamasks_size;
306 std::vector<helios::SphericalCoord> reconstructed_alphamasks_rotation;
307 std::vector<uint> reconstructed_alphamasks_gridcell;
308 std::string reconstructed_alphamasks_maskfile;
309 std::vector<uint> reconstructed_alphamasks_direct_flag;
310
311 void leafReconstructionFloodfill();
312
313 void backfillLeavesAlphaMask(const std::vector<float> &leaf_size, float leaf_aspect_ratio, float solidfraction, const std::vector<bool> &group_filter_flag);
314
315 void calculateLeafAngleCDF(uint Nbins, std::vector<std::vector<float>> &CDF_theta, std::vector<std::vector<float>> &CDF_phi);
316
317 void floodfill(size_t t, std::vector<Triangulation> &cloud_triangles, std::vector<int> &fill_flag, std::vector<std::vector<int>> &nodes, int tag, int depth, int maxdepth);
318
320
326 std::vector<float> LAD_inversion(std::vector<float> &P, std::vector<float> &Gtheta, std::vector<std::vector<float>> &dr_array, bool fillAnalytic);
327
328 // -------- HELPERS --------- //
329
331
337 void computeGtheta(uint Ncells, uint Nscans,
338 std::vector<float> &Gtheta, std::vector<float> &Gtheta_bar);
339
341
351 bool invertLAD(uint voxel_index, float P, float Gtheta,
352 const std::vector<float> &dr_samples, int min_voxel_hits,
353 const helios::vec3 &gridsize, float &leaf_area);
354
356
364 uint filterRaysByBoundingBox(const helios::vec3 &scan_origin,
365 const std::vector<helios::vec3> &ray_endpoints,
366 const helios::vec3 &bb_center,
367 const helios::vec3 &bb_size,
368 std::vector<uint> &filtered_indices);
369
371
381 void calculateVoxelPathLengths(const helios::vec3 &scan_origin,
382 const std::vector<helios::vec3> &ray_directions,
383 const std::vector<helios::vec3> &voxel_centers,
384 const std::vector<helios::vec3> &voxel_sizes,
385 const std::vector<float> &voxel_rotations,
386 std::vector<std::vector<float>> &dr_agg,
387 std::vector<float> &hit_before_agg,
388 std::vector<float> &hit_after_agg);
389
390 // -------- MULTI-RETURN HELPERS --------- //
391
393
396 bool isMultiReturnData() const;
397
399 struct BeamGrouping {
400 uint Nbeams;
401 std::vector<std::vector<uint>> beam_array;
402 };
403
405
409 BeamGrouping groupHitsByTimestamp(const std::vector<uint>& scan_indices) const;
410
412
420 std::vector<uint> loadTreeQSM_impl(helios::Context *context, const std::string &filename, uint radial_subdivisions, bool use_colormap, const std::string &colormap_or_texture);
421
422public:
424 LiDARcloud();
425
427 ~LiDARcloud();
428
430 static int selfTest(int argc = 0, char **argv = nullptr);
431
432 void validateRayDirections();
433
435 void disableMessages();
436
438 void enableMessages();
439
441 void initializeCollisionDetection(helios::Context *context);
442
444 void performUnifiedRayTracing(helios::Context *context, size_t N, int Npulse, helios::vec3 *ray_origins, helios::vec3 *direction, float *hit_t, float *hit_fnorm, int *hit_ID);
445
446 // ------- SCANS -------- //
447
449 uint getScanCount();
450
452
456 uint addScan(ScanMetadata &newscan);
457
459
465 void addHitPoint(uint scanID, const helios::vec3 &xyz, const helios::SphericalCoord &direction);
466
468
475 void addHitPoint(uint scanID, const helios::vec3 &xyz, const helios::SphericalCoord &direction, const helios::RGBcolor &color);
476
478
484 void addHitPoint(uint scanID, const helios::vec3 &xyz, const helios::SphericalCoord &direction, const std::map<std::string, double> &data);
485
487
494 void addHitPoint(uint scanID, const helios::vec3 &xyz, const helios::SphericalCoord &direction, const helios::RGBcolor &color, const std::map<std::string, double> &data);
495
497
504 void addHitPoint(uint scanID, const helios::vec3 &xyz, const helios::int2 &row_column, const helios::RGBcolor &color, const std::map<std::string, double> &data);
505
507
510 void deleteHitPoint(uint index);
511
513 uint getHitCount() const;
514
516
519 helios::vec3 getScanOrigin(uint scanID) const;
520
522
525 uint getScanSizeTheta(uint scanID) const;
526
528
531 uint getScanSizePhi(uint scanID) const;
532
534
538 helios::vec2 getScanRangeTheta(uint scanID) const;
539
541
545 helios::vec2 getScanRangePhi(uint scanID) const;
546
548
552 float getScanBeamExitDiameter(uint scanID) const;
553
555
558 std::vector<std::string> getScanColumnFormat(uint scanID) const;
559
561
565 float getScanBeamDivergence(uint scanID) const;
566
568
571 helios::vec3 getHitXYZ(uint index) const;
572
574
577 helios::SphericalCoord getHitRaydir(uint index) const;
578
580
585 double getHitData(uint index, const char *label) const;
586
588
593 void setHitData(uint index, const char *label, double value);
594
596
600 bool doesHitDataExist(uint index, const char *label) const;
601
603
606 helios::RGBcolor getHitColor(uint index) const;
607
609
612 int getHitScanID(uint index) const;
613
615
621 int getHitIndex(uint scanID, uint row, uint column) const;
622
624
629 int getHitGridCell(uint index) const;
630
632
636 void setHitGridCell(uint index, int cell);
637
639
642 void coordinateShift(const helios::vec3 &shift);
643
645
649 void coordinateShift(uint scanID, const helios::vec3 &shift);
650
652
655 void coordinateRotation(const helios::SphericalCoord &rotation);
656
658
662 void coordinateRotation(uint scanID, const helios::SphericalCoord &rotation);
663
665
670 void coordinateRotation(float rotation, const helios::vec3 &line_base, const helios::vec3 &line_direction);
671
673 uint getTriangleCount() const;
674
676
680 Triangulation getTriangle(uint index) const;
681
682 // ------- FILE I/O --------- //
683
685
688 void loadXML(const char *filename);
689
691
695 void loadXML(const char *filename, bool load_grid_only);
696
698
703 size_t loadASCIIFile(uint scanID, const std::string &ASCII_data_file);
704
706
709 void exportTriangleNormals(const char *filename);
710
712
716 void exportTriangleNormals(const char *filename, int gridcell);
717
719
722 void exportTriangleAreas(const char *filename);
723
725
729 void exportTriangleAreas(const char *filename, int gridcell);
730
733
737 void exportTriangleInclinationDistribution(const char *filename, uint Nbins);
738
741
745 void exportTriangleAzimuthDistribution(const char *filename, uint Nbins);
746
748
751 void exportLeafAreas(const char *filename);
752
754
757 void exportLeafAreaDensities(const char *filename);
758
760
763 void exportGtheta(const char *filename);
764
766
770 void exportPointCloud(const char *filename);
771
773
777 void exportPointCloud(const char *filename, uint scanID);
778
780
784 void exportPointCloudPTX(const char *filename, uint scanID);
785
786 // ------- VISUALIZER --------- //
787
789
793 void addHitsToVisualizer(Visualizer *visualizer, uint pointsize) const;
794
796
801 void addHitsToVisualizer(Visualizer *visualizer, uint pointsize, const helios::RGBcolor &point_color) const;
802
804
809 void addHitsToVisualizer(Visualizer *visualizer, uint pointsize, const char *color_value) const;
810
812
815 void addGridToVisualizer(Visualizer *visualizer) const;
816
818
822 void addGridWireFrametoVisualizer(Visualizer *visualizer, float linewidth_pixels = 1.0f) const;
823
825
831 void addGrid(const helios::vec3 &center, const helios::vec3 &size, const helios::int3 &ndiv, float rotation);
832
834
837 void addTrianglesToVisualizer(Visualizer *visualizer) const;
838
840
844 void addTrianglesToVisualizer(Visualizer *visualizer, uint gridcell) const;
845
847
850 void addLeafReconstructionToVisualizer(Visualizer *visualizer) const;
851
853
856 void addTrunkReconstructionToVisualizer(Visualizer *visualizer) const;
857
859
863 void addTrunkReconstructionToVisualizer(Visualizer *visualizer, const helios::RGBcolor &trunk_color) const;
864
866
871 std::vector<uint> addLeafReconstructionToContext(helios::Context *context) const;
872
874
880 std::vector<uint> addLeafReconstructionToContext(helios::Context *context, const helios::int2 &subpatches) const;
881
883
887 std::vector<uint> addReconstructedTriangleGroupsToContext(helios::Context *context) const;
888
890
893 std::vector<uint> addTrunkReconstructionToContext(helios::Context *context) const;
894
896
900 void getHitBoundingBox(helios::vec3 &boxmin, helios::vec3 &boxmax) const;
901
903
907 void getGridBoundingBox(helios::vec3 &boxmin, helios::vec3 &boxmax) const;
908
909 // --------- POINT FILTERING ----------- //
910
912
915 void distanceFilter(float maxdistance);
916
918
927 void xyzFilter(float xmin, float xmax, float ymin, float ymax, float zmin, float zmax);
928
930
940 void xyzFilter(float xmin, float xmax, float ymin, float ymax, float zmin, float zmax, bool deleteOutside);
941
942
944
948 void reflectanceFilter(float minreflectance);
949
951
957 void scalarFilter(const char *scalar_field, float threshold, const char *comparator);
958
960
964 void maxPulseFilter(const char *scalar);
965
967
971 void minPulseFilter(const char *scalar);
972
974
978 void firstHitFilter();
979
981
985 void lastHitFilter();
986
987 // ------- TRIANGULATION --------- //
988
990
994 void triangulateHitPoints(float Lmax, float max_aspect_ratio);
995
996 // ERK
998
1006 void triangulateHitPoints(float Lmax, float max_aspect_ratio, const char *scalar_field, float threshold, const char *comparator);
1007
1008
1010
1013 void addTrianglesToContext(helios::Context *context) const;
1014
1015 // -------- GRID ----------- //
1016
1018 uint getGridCellCount() const;
1019
1021
1026 void addGridCell(const helios::vec3 &center, const helios::vec3 &size, float rotation);
1027
1029
1038 void addGridCell(const helios::vec3 &center, const helios::vec3 &global_anchor, const helios::vec3 &size, const helios::vec3 &global_size, float rotation, const helios::int3 &global_ijk, const helios::int3 &global_count);
1039
1041
1044 helios::vec3 getCellCenter(uint index) const;
1045
1047
1050 helios::vec3 getCellGlobalAnchor(uint index) const;
1051
1053
1057 helios::vec3 getCellSize(uint index) const;
1058
1060
1063 float getCellRotation(uint index) const;
1064
1065 // ------- SYNTHETIC SCAN ------ //
1066
1068
1071 void syntheticScan(helios::Context *context);
1072
1074
1078 void syntheticScan(helios::Context *context, bool append);
1079
1081
1087 void syntheticScan(helios::Context *context, bool scan_grid_only, bool record_misses);
1088
1090
1097 void syntheticScan(helios::Context *context, bool scan_grid_only, bool record_misses, bool append);
1098
1100
1106 void syntheticScan(helios::Context *context, int rays_per_pulse, float pulse_distance_threshold);
1107
1109
1116 void syntheticScan(helios::Context *context, int rays_per_pulse, float pulse_distance_threshold, bool append);
1117
1119
1127 void syntheticScan(helios::Context *context, int rays_per_pulse, float pulse_distance_threshold, bool scan_grid_only, bool record_misses);
1128
1130
1139 void syntheticScan(helios::Context *context, int rays_per_pulse, float pulse_distance_threshold, bool scan_grid_only, bool record_misses, bool append);
1140
1142
1145 std::vector<float> calculateSyntheticLeafArea(helios::Context *context);
1146
1148
1151 std::vector<float> calculateSyntheticGtheta(helios::Context *context);
1152
1153 // -------- LEAF AREA -------- //
1154
1156
1160 void setCellLeafArea(float area, uint index);
1161
1163
1166 float getCellLeafArea(uint index) const;
1167
1169
1172 float getCellLeafAreaDensity(uint index) const;
1173
1175
1179 void setCellGtheta(float Gtheta, uint index);
1180
1182
1185 float getCellGtheta(uint index) const;
1186
1188
1191 std::vector<helios::vec3> gapfillMisses();
1192
1194
1198 std::vector<helios::vec3> gapfillMisses(uint scanID);
1199
1201
1207 std::vector<helios::vec3> gapfillMisses(uint scanID, const bool gapfill_grid_only, const bool add_flags);
1208
1209
1211
1214 void calculateLeafArea(helios::Context *context);
1215
1217
1221 void calculateLeafArea(helios::Context *context, int min_voxel_hits);
1222
1224
1228 [[deprecated("Use calculateLeafArea() instead. GPU functionality is now provided by the CollisionDetection plugin.")]]
1229 void calculateLeafAreaGPU(helios::Context *context);
1230
1232
1237 [[deprecated("Use calculateLeafArea(context, min_voxel_hits) instead. GPU functionality is now provided by the CollisionDetection plugin.")]]
1238 void calculateLeafAreaGPU(helios::Context *context, int min_voxel_hits);
1239
1241 void enableGPUAcceleration();
1242
1244 void disableGPUAcceleration();
1245
1247 void calculateHitGridCell();
1248
1249 // -------- RECONSTRUCTION --------- //
1250
1252
1258 void leafReconstructionAlphaMask(float minimum_leaf_group_area, float maximum_leaf_group_area, float leaf_aspect_ratio, const char *mask_file);
1259
1261
1268 void leafReconstructionAlphaMask(float minimum_leaf_group_area, float maximum_leaf_group_area, float leaf_aspect_ratio, float leaf_length_constant, const char *mask_file);
1269
1272
1278 void trunkReconstruction(const helios::vec3 &box_center, const helios::vec3 &box_size, float Lmax, float max_aspect_ratio);
1279
1281
1288 std::vector<uint> loadTreeQSM(helios::Context *context, const std::string &filename, uint radial_subdivisions, const std::string &texture_file = "");
1289
1291
1299 std::vector<uint> loadTreeQSMColormap(helios::Context *context, const std::string &filename, uint radial_subdivisions, const std::string &colormap_name);
1300
1302
1305 void cropBeamsToGridAngleRange(uint source);
1306
1308
1311 std::vector<uint> peakFinder(std::vector<float> signal);
1312};
1313
1314bool sortcol0(const std::vector<double> &v0, const std::vector<double> &v1);
1315
1316bool sortcol1(const std::vector<double> &v0, const std::vector<double> &v1);
1317
1318#endif