lasutility.cpp
Go to the documentation of this file.
1 /*
2 ===============================================================================
3 
4  FILE: lasutility.cpp
5 
6  CONTENTS:
7 
8  see corresponding header file
9 
10  PROGRAMMERS:
11 
12  martin.isenburg@gmail.com
13 
14  COPYRIGHT:
15 
16  (c) 2010-2011, Martin Isenburg, LASSO - tools to catch reality
17 
18  This is free software; you can redistribute and/or modify it under the
19  terms of the GNU Lesser General Licence as published by the Free Software
20  Foundation. See the COPYING file for more information.
21 
22  This software is distributed WITHOUT ANY WARRANTY and without even the
23  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
24 
25  CHANGE HISTORY:
26 
27  see corresponding header file
28 
29 ===============================================================================
30 */
31 #include "lasutility.hpp"
32 
33 #include <stdio.h>
34 #include <stdlib.h>
35 #include <string.h>
36 
38 {
39  U32 i;
41  for (i = 0; i < 8; i++) number_of_points_by_return[i] = 0;
42  raw_max_x = raw_min_x = 0;
43  raw_max_y = raw_min_y = 0;
44  raw_max_z = raw_min_z = 0;
45  first = true;
46 }
47 
49 {
52  if (first)
53  {
54  raw_min_x = raw_max_x = point->x;
55  raw_min_y = raw_max_y = point->y;
56  raw_min_z = raw_max_z = point->z;
57  first = false;
58  }
59  else
60  {
61  if (point->x < raw_min_x) raw_min_x = point->x;
62  else if (point->x > raw_max_x) raw_max_x = point->x;
63  if (point->y < raw_min_y) raw_min_y = point->y;
64  else if (point->y > raw_max_y) raw_max_y = point->y;
65  if (point->z < raw_min_z) raw_min_z = point->z;
66  else if (point->z > raw_max_z) raw_max_z = point->z;
67  }
68  return TRUE;
69 }
70 
72 {
73  U32 i;
75  for (i = 0; i < 8; i++) number_of_points_by_return[i] = 0;
76  for (i = 0; i < 8; i++) number_of_returns_of_given_pulse[i] = 0;
77  for (i = 0; i < 32; i++) classification[i] = 0;
81  first = true;
82 }
83 
85 {
89  classification[point->classification&31]++;
90  if (point->classification & 32) classification_synthetic++;
91  if (point->classification & 64) classification_keypoint++;
92  if (point->classification & 128) classification_withheld++;
93  if (first)
94  {
95  min = *point;
96  max = *point;
97  first = false;
98  }
99  else
100  {
101  if (point->x < min.x) min.x = point->x;
102  else if (point->x > max.x) max.x = point->x;
103  if (point->y < min.y) min.y = point->y;
104  else if (point->y > max.y) max.y = point->y;
105  if (point->z < min.z) min.z = point->z;
106  else if (point->z > max.z) max.z = point->z;
107  if (point->intensity < min.intensity) min.intensity = point->intensity;
108  else if (point->intensity > max.intensity) max.intensity = point->intensity;
116  else if (point->return_number > max.return_number) max.return_number = point->return_number;
121  if (point->user_data < min.user_data) min.user_data = point->user_data;
122  else if (point->user_data > max.user_data) max.user_data = point->user_data;
127  if (point->have_gps_time)
128  {
129  if (point->gps_time < min.gps_time) min.gps_time = point->gps_time;
130  else if (point->gps_time > max.gps_time) max.gps_time = point->gps_time;
131  }
132  if (point->have_rgb)
133  {
134  if (point->rgb[0] < min.rgb[0]) min.rgb[0] = point->rgb[0];
135  else if (point->rgb[0] > max.rgb[0]) max.rgb[0] = point->rgb[0];
136  if (point->rgb[1] < min.rgb[1]) min.rgb[1] = point->rgb[1];
137  else if (point->rgb[1] > max.rgb[1]) max.rgb[1] = point->rgb[1];
138  if (point->rgb[2] < min.rgb[2]) min.rgb[2] = point->rgb[2];
139  else if (point->rgb[2] > max.rgb[2]) max.rgb[2] = point->rgb[2];
140  }
141  if (point->have_wavepacket)
142  {
148  else if (point->wavepacket.getSize() > max.wavepacket.getSize()) max.wavepacket.setSize(point->wavepacket.getSize());
151  if (point->wavepacket.getXt() < min.wavepacket.getXt()) min.wavepacket.setXt(point->wavepacket.getXt());
152  else if (point->wavepacket.getXt() > max.wavepacket.getXt()) max.wavepacket.setXt(point->wavepacket.getXt());
153  if (point->wavepacket.getYt() < min.wavepacket.getYt()) min.wavepacket.setYt(point->wavepacket.getYt());
154  else if (point->wavepacket.getYt() > max.wavepacket.getYt()) max.wavepacket.setYt(point->wavepacket.getYt());
155  if (point->wavepacket.getZt() < min.wavepacket.getZt()) min.wavepacket.setZt(point->wavepacket.getZt());
156  else if (point->wavepacket.getZt() > max.wavepacket.getZt()) max.wavepacket.setZt(point->wavepacket.getZt());
157  }
158  }
159  return TRUE;
160 }
161 
163 {
164  total = 0;
165  count = 0;
166  this->one_over_step = 1.0f/step;
167  first = TRUE;
168  size_pos = 0;
169  size_neg = 0;
170  bins_pos = 0;
171  bins_neg = 0;
172  values_pos = 0;
173  values_neg = 0;
174 }
175 
177 {
178  if (bins_pos) free(bins_pos);
179  if (bins_neg) free(bins_neg);
180  if (values_pos) free(values_pos);
181  if (values_neg) free(values_neg);
182 }
183 
184 void LASbin::add(I32 item)
185 {
186  total += item;
187  count++;
188  I32 bin = I32_FLOOR(one_over_step*item);
189  add_to_bin(bin);
190 }
191 
192 void LASbin::add(F64 item)
193 {
194  total += item;
195  count++;
196  I32 bin = I32_FLOOR(one_over_step*item);
197  add_to_bin(bin);
198 }
199 
200 void LASbin::add(I64 item)
201 {
202  total += item;
203  count++;
204  I32 bin = I32_FLOOR(one_over_step*item);
205  add_to_bin(bin);
206 }
207 
209 {
210  if (first)
211  {
212  anker = bin;
213  first = FALSE;
214  }
215  bin = bin - anker;
216  if (bin >= 0)
217  {
218  if (bin >= size_pos)
219  {
220  I32 i;
221  if (size_pos == 0)
222  {
223  size_pos = 1024;
224  bins_pos = (U32*)malloc(sizeof(U32)*size_pos);
225  for (i = 0; i < size_pos; i++) bins_pos[i] = 0;
226  }
227  else
228  {
229  I32 new_size = bin + 1024;
230  bins_pos = (U32*)realloc(bins_pos, sizeof(U32)*new_size);
231  for (i = size_pos; i < new_size; i++) bins_pos[i] = 0;
232  size_pos = new_size;
233  }
234  }
235  bins_pos[bin]++;
236  }
237  else
238  {
239  bin = -(bin+1);
240  if (bin >= size_neg)
241  {
242  I32 i;
243  if (size_neg == 0)
244  {
245  size_neg = 1024;
246  bins_neg = (U32*)malloc(sizeof(U32)*size_neg);
247  for (i = 0; i < size_neg; i++) bins_neg[i] = 0;
248  }
249  else
250  {
251  I32 new_size = bin + 1024;
252  bins_neg = (U32*)realloc(bins_neg, sizeof(U32)*new_size);
253  for (i = size_neg; i < new_size; i++) bins_neg[i] = 0;
254  size_neg = new_size;
255  }
256  }
257  bins_neg[bin]++;
258  }
259 }
260 
261 void LASbin::add(I32 item, I32 value)
262 {
263  total += item;
264  count++;
265  I32 bin = I32_FLOOR(one_over_step*item);
266  if (first)
267  {
268  anker = bin;
269  first = FALSE;
270  }
271  bin = bin - anker;
272  if (bin >= 0)
273  {
274  if (bin >= size_pos)
275  {
276  I32 i;
277  if (size_pos == 0)
278  {
279  size_pos = 1024;
280  bins_pos = (U32*)malloc(sizeof(U32)*size_pos);
281  values_pos = (F64*)malloc(sizeof(F64)*size_pos);
282  for (i = 0; i < size_pos; i++) { bins_pos[i] = 0; values_pos[i] = 0; }
283  }
284  else
285  {
286  I32 new_size = bin + 1024;
287  bins_pos = (U32*)realloc(bins_pos, sizeof(U32)*new_size);
288  values_pos = (F64*)realloc(values_pos, sizeof(F64)*new_size);
289  for (i = size_pos; i < new_size; i++) { bins_pos[i] = 0; values_pos[i] = 0; }
290  size_pos = new_size;
291  }
292  }
293  bins_pos[bin]++;
294  values_pos[bin] += value;
295  }
296  else
297  {
298  bin = -(bin+1);
299  if (bin >= size_neg)
300  {
301  I32 i;
302  if (size_neg == 0)
303  {
304  size_neg = 1024;
305  bins_neg = (U32*)malloc(sizeof(U32)*size_neg);
306  values_neg = (F64*)malloc(sizeof(F64)*size_pos);
307  for (i = 0; i < size_neg; i++) { bins_neg[i] = 0; values_neg[i] = 0; }
308  }
309  else
310  {
311  I32 new_size = bin + 1024;
312  bins_neg = (U32*)realloc(bins_neg, sizeof(U32)*new_size);
313  values_neg = (F64*)realloc(values_neg, sizeof(F64)*new_size);
314  for (i = size_neg; i < new_size; i++) { bins_neg[i] = 0; values_neg[i] = 0; }
315  size_neg = new_size;
316  }
317  }
318  bins_neg[bin]++;
319  values_neg[bin] += value;
320  }
321 }
322 
323 void LASbin::report(FILE* file, const char* name, const char* name_avg) const
324 {
325  I32 i, bin;
326  if (name)
327  {
328  if (values_pos)
329  {
330  if (name_avg)
331  fprintf(file, "%s histogram of %s averages with bin size %g\012", name, name_avg, 1.0f/one_over_step);
332  else
333  fprintf(file, "%s histogram of averages with bin size %g\012", name, 1.0f/one_over_step);
334  }
335  else
336  fprintf(file, "%s histogram with bin size %g\012", name, 1.0f/one_over_step);
337  }
338  if (size_neg)
339  {
340  for (i = size_neg-1; i >= 0; i--)
341  {
342  if (bins_neg[i])
343  {
344  bin = -(i+1) + anker;
345  if (one_over_step == 1)
346  {
347  if (values_neg)
348  fprintf(file, " bin %d has average %g (of %d)\012", bin, values_neg[i]/bins_neg[i], bins_neg[i]);
349  else
350  fprintf(file, " bin %d has %d\012", bin, bins_neg[i]);
351  }
352  else
353  {
354  if (values_neg)
355  fprintf(file, " bin [%g,%g) has average %g (of %d)\012", ((F32)bin)/one_over_step, ((F32)(bin+1))/one_over_step, values_neg[i]/bins_neg[i], bins_neg[i]);
356  else
357  fprintf(file, " bin [%g,%g) has %d\012", ((F32)bin)/one_over_step, ((F32)(bin+1))/one_over_step, bins_neg[i]);
358  }
359  }
360  }
361  }
362  if (size_pos)
363  {
364  for (i = 0; i < size_pos; i++)
365  {
366  if (bins_pos[i])
367  {
368  bin = i + anker;
369  if (one_over_step == 1)
370  {
371  if (values_pos)
372  fprintf(file, " bin %d has average %g (of %d)\012", bin, values_pos[i]/bins_pos[i], bins_pos[i]);
373  else
374  fprintf(file, " bin %d has %d\012", bin, bins_pos[i]);
375  }
376  else
377  {
378  if (values_pos)
379  fprintf(file, " bin [%g,%g) average has %g (of %d)\012", ((F32)bin)/one_over_step, ((F32)(bin+1))/one_over_step, values_pos[i]/bins_pos[i], bins_pos[i]);
380  else
381  fprintf(file, " bin [%g,%g) has %d\012", ((F32)bin)/one_over_step, ((F32)(bin+1))/one_over_step, bins_pos[i]);
382  }
383  }
384  }
385  }
386  if (name)
387  fprintf(file, " average %s %g\012", name, total/count);
388  else
389  fprintf(file, " average %g\012", total/count);
390 }
391 
393 {
394  is_active = FALSE;
395  // counter bins
396  x_bin = 0;
397  y_bin = 0;
398  z_bin = 0;
399  intensity_bin = 0;
400  classification_bin = 0;
401  scan_angle_bin = 0;
403  gps_time_bin = 0;
408  // averages bins
411  scan_angle_bin_z = 0;
415 }
416 
418 {
419  // counter bins
420  if (x_bin) delete x_bin;
421  if (y_bin) delete y_bin;
422  if (z_bin) delete z_bin;
423  if (intensity_bin) delete intensity_bin;
425  if (scan_angle_bin) delete scan_angle_bin;
427  if (gps_time_bin) delete gps_time_bin;
432  // averages bins
439 }
440 
442 {
443  int i;
444  for (i = 1; i < argc; i++)
445  {
446  if (argv[i][0] == '\0')
447  {
448  continue;
449  }
450  else if (strcmp(argv[i],"-h") == 0 || strcmp(argv[i],"-help") == 0)
451  {
452  return TRUE;
453  }
454  else if (strcmp(argv[i],"-histo") == 0)
455  {
456  if ((i+2) >= argc)
457  {
458  fprintf(stderr,"ERROR: '%s' needs 2 arguments: name step\n", argv[i]);
459  return FALSE;
460  }
461  if (!histo(argv[i+1], (F32)atof(argv[i+2]))) return FALSE;
462  *argv[i]='\0'; *argv[i+1]='\0'; *argv[i+2]='\0'; i+=2;
463  }
464  else if (strcmp(argv[i],"-histo_avg") == 0)
465  {
466  if ((i+3) >= argc)
467  {
468  fprintf(stderr,"ERROR: '%s' needs 3 arguments: name step name_avg\n", argv[i]);
469  return FALSE;
470  }
471  if (!histo_avg(argv[i+1], (F32)atof(argv[i+2]), argv[i+3])) return FALSE;
472  *argv[i]='\0'; *argv[i+1]='\0'; *argv[i+2]='\0'; *argv[i+3]='\0'; i+=3;
473  }
474  }
475  return TRUE;
476 }
477 
478 BOOL LAShistogram::histo(const char* name, F32 step)
479 {
480  if (strcmp(name, "x") == 0)
481  x_bin = new LASbin(step);
482  else if (strcmp(name, "y") == 0)
483  y_bin = new LASbin(step);
484  else if (strcmp(name, "z") == 0)
485  z_bin = new LASbin(step);
486  else if (strcmp(name, "intensity") == 0)
487  intensity_bin = new LASbin(step);
488  else if (strcmp(name, "classification") == 0)
490  else if (strstr(name, "scan_angle") != 0)
491  scan_angle_bin = new LASbin(step);
492  else if (strstr(name, "point_source") != 0)
494  else if (strstr(name, "gps_time") != 0)
495  gps_time_bin = new LASbin(step);
496  else if (strstr(name, "wavepacket_index") != 0)
498  else if (strstr(name, "wavepacket_offset") != 0)
500  else if (strstr(name, "wavepacket_size") != 0)
502  else if (strstr(name, "wavepacket_location") != 0)
504  else
505  {
506  fprintf(stderr,"ERROR: histogram of '%s' not implemented\n", name);
507  return FALSE;
508  }
509  is_active = TRUE;
510  return TRUE;
511 }
512 
513 BOOL LAShistogram::histo_avg(const char* name, F32 step, const char* name_avg)
514 {
515  if (strcmp(name, "classification") == 0)
516  {
517  if (strcmp(name_avg, "intensity") == 0)
519  else if (strstr(name_avg, "scan_angle") != 0)
521  else
522  {
523  fprintf(stderr,"ERROR: histogram of '%s' with '%s' averages not implemented\n", name, name_avg);
524  return FALSE;
525  }
526  }
527  else if (strcmp(name, "scan_angle") == 0)
528  {
529  if (strcmp(name_avg, "z") == 0)
531  else if (strcmp(name_avg, "number_of_returns") == 0)
533  else if (strcmp(name_avg, "intensity") == 0)
535  else
536  {
537  fprintf(stderr,"ERROR: histogram of '%s' with '%s' averages not implemented\n", name, name_avg);
538  return FALSE;
539  }
540  }
541  else if (strcmp(name, "return_map") == 0)
542  {
543  if (strcmp(name_avg, "intensity") == 0)
545  else
546  {
547  fprintf(stderr,"ERROR: histogram of '%s' with '%s' averages not implemented\n", name, name_avg);
548  return FALSE;
549  }
550  }
551  else
552  {
553  fprintf(stderr,"ERROR: histogram of '%s' not implemented\n", name);
554  return FALSE;
555  }
556  is_active = TRUE;
557  return TRUE;
558 }
559 
560 void LAShistogram::add(const LASpoint* point)
561 {
562  // counter bins
563  if (x_bin) x_bin->add(point->x);
564  if (y_bin) y_bin->add(point->y);
565  if (z_bin) z_bin->add(point->z);
570  if (gps_time_bin) gps_time_bin->add(point->gps_time);
575  // averages bins
578  if (scan_angle_bin_z) scan_angle_bin_z->add(point->scan_angle_rank, point->z);
582  {
583  int n = point->number_of_returns_of_given_pulse;
584  int r = point->return_number;
585  return_map_bin_intensity->add((n == 1 ? 0 : (n == 2 ? r : (n == 3 ? r+2 : (n == 4 ? r+5 : (n == 5 ? r+9 : 15))))), point->intensity);
586  }
587 }
588 
589 void LAShistogram::report(FILE* file) const
590 {
591  // counter bins
592  if (x_bin) x_bin->report(file, "x coordinate");
593  if (y_bin) y_bin->report(file, "y coordinate");
594  if (z_bin) z_bin->report(file, "z coordinate");
595  if (intensity_bin) intensity_bin->report(file, "intensity");
596  if (classification_bin) classification_bin->report(file, "classification");
597  if (scan_angle_bin) scan_angle_bin->report(file, "scan angle");
598  if (point_source_id_bin) point_source_id_bin->report(file, "point source id");
599  if (gps_time_bin) gps_time_bin->report(file, "gps_time");
600  if (wavepacket_index_bin) wavepacket_index_bin->report(file, "wavepacket_index");
601  if (wavepacket_offset_bin) wavepacket_offset_bin->report(file, "wavepacket_offset");
602  if (wavepacket_size_bin) wavepacket_size_bin->report(file, "wavepacket_size");
603  if (wavepacket_location_bin) wavepacket_location_bin->report(file, "wavepacket_location");
604  // averages bins
605  if (classification_bin_intensity) classification_bin_intensity->report(file, "classification", "intensity");
606  if (classification_bin_scan_angle) classification_bin_scan_angle->report(file, "classification", "scan_angle");
607  if (scan_angle_bin_z) scan_angle_bin_z->report(file, "scan angle", "z coordinate");
608  if (scan_angle_bin_number_of_returns) scan_angle_bin_number_of_returns->report(file, "scan_angle", "number_of_returns");
609  if (scan_angle_bin_intensity) scan_angle_bin_intensity->report(file, "scan angle", "intensity");
610  if (return_map_bin_intensity) return_map_bin_intensity->report(file, "return map", "intensity");
611 }
612 
614 {
615  I32 pos_x, pos_y;
616  if (grid_spacing < 0)
617  {
619  pos_x = I32_FLOOR(point->get_x() / grid_spacing);
620  pos_y = I32_FLOOR(point->get_y() / grid_spacing);
621  anker = pos_y;
622  min_x = max_x = pos_x;
623  min_y = max_y = pos_y;
624  }
625  else
626  {
627  pos_x = I32_FLOOR(point->get_x() / grid_spacing);
628  pos_y = I32_FLOOR(point->get_y() / grid_spacing);
629  if (pos_x < min_x) min_x = pos_x; else if (pos_x > max_x) max_x = pos_x;
630  if (pos_y < min_y) min_y = pos_y; else if (pos_y > max_y) max_y = pos_y;
631  }
632  return add_internal(pos_x, pos_y);
633 }
634 
636 {
637  if (grid_spacing < 0)
638  {
640  anker = pos_y;
641  min_x = max_x = pos_x;
642  min_y = max_y = pos_y;
643  }
644  else
645  {
646  if (pos_x < min_x) min_x = pos_x; else if (pos_x > max_x) max_x = pos_x;
647  if (pos_y < min_y) min_y = pos_y; else if (pos_y > max_y) max_y = pos_y;
648  }
649  return add_internal(pos_x, pos_y);
650 }
651 
653 {
654  pos_y = pos_y - anker;
655  BOOL no_x_anker = FALSE;
656  U32* array_size;
657  I32** ankers;
658  U32*** array;
659  U16** array_sizes;
660  if (pos_y < 0)
661  {
662  pos_y = -pos_y - 1;
663  ankers = &minus_ankers;
664  if ((U32)pos_y < minus_plus_size && minus_plus_sizes[pos_y])
665  {
666  pos_x -= minus_ankers[pos_y];
667  if (pos_x < 0)
668  {
669  pos_x = -pos_x - 1;
670  array_size = &minus_minus_size;
671  array = &minus_minus;
672  array_sizes = &minus_minus_sizes;
673  }
674  else
675  {
676  array_size = &minus_plus_size;
677  array = &minus_plus;
678  array_sizes = &minus_plus_sizes;
679  }
680  }
681  else
682  {
683  no_x_anker = TRUE;
684  array_size = &minus_plus_size;
685  array = &minus_plus;
686  array_sizes = &minus_plus_sizes;
687  }
688  }
689  else
690  {
691  ankers = &plus_ankers;
692  if ((U32)pos_y < plus_plus_size && plus_plus_sizes[pos_y])
693  {
694  pos_x -= plus_ankers[pos_y];
695  if (pos_x < 0)
696  {
697  pos_x = -pos_x - 1;
698  array_size = &plus_minus_size;
699  array = &plus_minus;
700  array_sizes = &plus_minus_sizes;
701  }
702  else
703  {
704  array_size = &plus_plus_size;
705  array = &plus_plus;
706  array_sizes = &plus_plus_sizes;
707  }
708  }
709  else
710  {
711  no_x_anker = TRUE;
712  array_size = &plus_plus_size;
713  array = &plus_plus;
714  array_sizes = &plus_plus_sizes;
715  }
716  }
717  // maybe grow banded grid in y direction
718  if ((U32)pos_y >= *array_size)
719  {
720  U32 array_size_new = ((pos_y/1024)+1)*1024;
721  if (*array_size)
722  {
723  if (array == &minus_plus || array == &plus_plus) *ankers = (I32*)realloc(*ankers, array_size_new*sizeof(I32));
724  *array = (U32**)realloc(*array, array_size_new*sizeof(U32*));
725  *array_sizes = (U16*)realloc(*array_sizes, array_size_new*sizeof(U16));
726  }
727  else
728  {
729  if (array == &minus_plus || array == &plus_plus) *ankers = (I32*)malloc(array_size_new*sizeof(I32));
730  *array = (U32**)malloc(array_size_new*sizeof(U32*));
731  *array_sizes = (U16*)malloc(array_size_new*sizeof(U16));
732  }
733  for (U32 i = *array_size; i < array_size_new; i++)
734  {
735  (*array)[i] = 0;
736  (*array_sizes)[i] = 0;
737  }
738  *array_size = array_size_new;
739  }
740  // is this the first x anker for this y pos?
741  if (no_x_anker)
742  {
743  (*ankers)[pos_y] = pos_x;
744  pos_x = 0;
745  }
746  // maybe grow banded grid in x direction
747  U32 pos_x_pos = pos_x/32;
748  if (pos_x_pos >= (*array_sizes)[pos_y])
749  {
750  U32 array_sizes_new = ((pos_x_pos/256)+1)*256;
751  if ((*array_sizes)[pos_y])
752  {
753  (*array)[pos_y] = (U32*)realloc((*array)[pos_y], array_sizes_new*sizeof(U32));
754  }
755  else
756  {
757  (*array)[pos_y] = (U32*)malloc(array_sizes_new*sizeof(U32));
758  }
759  for (U16 i = (*array_sizes)[pos_y]; i < array_sizes_new; i++)
760  {
761  (*array)[pos_y][i] = 0;
762  }
763  (*array_sizes)[pos_y] = array_sizes_new;
764  }
765  U32 pos_x_bit = 1 << (pos_x%32);
766  if ((*array)[pos_y][pos_x_pos] & pos_x_bit) return FALSE;
767  (*array)[pos_y][pos_x_pos] |= pos_x_bit;
768  num_occupied++;
769  return TRUE;
770 }
771 
773 {
774  I32 pos_x = I32_FLOOR(point->get_x() / grid_spacing);
775  I32 pos_y = I32_FLOOR(point->get_y() / grid_spacing);
776  return occupied(pos_x, pos_y);
777 }
778 
780 {
781  if (grid_spacing < 0)
782  {
783  return FALSE;
784  }
785  pos_y = pos_y - anker;
786  U32 array_size;
787  const I32* ankers;
788  const U32* const * array;
789  const U16* array_sizes;
790  if (pos_y < 0)
791  {
792  pos_y = -pos_y - 1;
793  ankers = minus_ankers;
794  if ((U32)pos_y < minus_plus_size && minus_plus_sizes[pos_y])
795  {
796  pos_x -= minus_ankers[pos_y];
797  if (pos_x < 0)
798  {
799  pos_x = -pos_x - 1;
800  array_size = minus_minus_size;
801  array = minus_minus;
802  array_sizes = minus_minus_sizes;
803  }
804  else
805  {
806  array_size = minus_plus_size;
807  array = minus_plus;
808  array_sizes = minus_plus_sizes;
809  }
810  }
811  else
812  {
813  return FALSE;
814  }
815  }
816  else
817  {
818  ankers = plus_ankers;
819  if ((U32)pos_y < plus_plus_size && plus_plus_sizes[pos_y])
820  {
821  pos_x -= plus_ankers[pos_y];
822  if (pos_x < 0)
823  {
824  pos_x = -pos_x - 1;
825  array_size = plus_minus_size;
826  array = plus_minus;
827  array_sizes = plus_minus_sizes;
828  }
829  else
830  {
831  array_size = plus_plus_size;
832  array = plus_plus;
833  array_sizes = plus_plus_sizes;
834  }
835  }
836  else
837  {
838  return FALSE;
839  }
840  }
841  // maybe out of bounds in y direction
842  if ((U32)pos_y >= array_size)
843  {
844  return FALSE;
845  }
846  // maybe out of bounds in x direction
847  U32 pos_x_pos = pos_x/32;
848  if (pos_x_pos >= array_sizes[pos_y])
849  {
850  return FALSE;
851  }
852  U32 pos_x_bit = 1 << (pos_x%32);
853  if (array[pos_y][pos_x_pos] & pos_x_bit) return TRUE;
854  return FALSE;
855 }
856 
858 {
859  if (grid_spacing < 0) return FALSE;
860  return TRUE;
861 }
862 
864 {
865  min_x = min_y = max_x = max_y = 0;
867  if (minus_minus_size)
868  {
869  for (U32 i = 0; i < minus_minus_size; i++) if (minus_minus[i]) free(minus_minus[i]);
870  free(minus_minus);
871  minus_minus = 0;
872  free(minus_minus_sizes);
873  minus_minus_sizes = 0;
874  minus_minus_size = 0;
875  }
876  if (minus_plus_size)
877  {
878  free(minus_ankers);
879  minus_ankers = 0;
880  for (U32 i = 0; i < minus_plus_size; i++) if (minus_plus[i]) free(minus_plus[i]);
881  free(minus_plus);
882  minus_plus = 0;
883  free(minus_plus_sizes);
884  minus_plus_sizes = 0;
885  minus_plus_size = 0;
886  }
887  if (plus_minus_size)
888  {
889  for (U32 i = 0; i < plus_minus_size; i++) if (plus_minus[i]) free(plus_minus[i]);
890  free(plus_minus);
891  plus_minus = 0;
892  free(plus_minus_sizes);
893  plus_minus_sizes = 0;
894  plus_minus_size = 0;
895  }
896  if (plus_plus_size)
897  {
898  free(plus_ankers);
899  plus_ankers = 0;
900  for (U32 i = 0; i < plus_plus_size; i++) if (plus_plus[i]) free(plus_plus[i]);
901  free(plus_plus);
902  plus_plus = 0;
903  free(plus_plus_sizes);
904  plus_plus_sizes = 0;
905  plus_plus_size = 0;
906  }
907  num_occupied = 0;
908 }
909 
910 BOOL LASoccupancyGrid::write_asc_grid(const char* file_name) const
911 {
912  FILE* file = fopen(file_name, "w");
913  if (file == 0) return FALSE;
914  fprintf(file, "ncols %d\012", max_x-min_x+1);
915  fprintf(file, "nrows %d\012", max_y-min_y+1);
916  fprintf(file, "xllcorner %f\012", grid_spacing*min_x);
917  fprintf(file, "yllcorner %f\012", grid_spacing*min_y);
918  fprintf(file, "cellsize %lf\012", grid_spacing);
919  fprintf(file, "NODATA_value %d\012", 0);
920  fprintf(file, "\012");
921  I32 pos_x, pos_y;
922  for (pos_y = min_y; pos_y <= max_y; pos_y++)
923  {
924  for (pos_x = min_x; pos_x <= max_x; pos_x++)
925  {
926  if (occupied(pos_x, pos_y))
927  {
928  fprintf(file, "1 ");
929  }
930  else
931  {
932  fprintf(file, "0 ");
933  }
934  }
935  fprintf(file, "\012");
936  }
937  fclose(file);
938  return TRUE;
939 }
940 
942 {
943  min_x = min_y = max_x = max_y = 0;
944  this->grid_spacing = -grid_spacing;
945  minus_ankers = 0;
946  minus_minus_size = 0;
947  minus_minus = 0;
948  minus_minus_sizes = 0;
949  minus_plus_size = 0;
950  minus_plus = 0;
951  minus_plus_sizes = 0;
952  plus_ankers = 0;
953  plus_minus_size = 0;
954  plus_minus = 0;
955  plus_minus_sizes = 0;
956  plus_plus_size = 0;
957  plus_plus = 0;
958  plus_plus_sizes = 0;
959  num_occupied = 0;
960 }
961 
963 {
964  reset();
965 }
966 
LASpoint::edge_of_flight_line
U8 edge_of_flight_line
Definition: lasdefinitions.hpp:485
LASsummary::number_of_returns_of_given_pulse
U32 number_of_returns_of_given_pulse[8]
Definition: lasutility.hpp:60
LASbin::bins_neg
U32 * bins_neg
Definition: lasutility.hpp:93
LAShistogram::classification_bin_intensity
LASbin * classification_bin_intensity
Definition: lasutility.hpp:125
LAShistogram::intensity_bin
LASbin * intensity_bin
Definition: lasutility.hpp:115
LASpoint::user_data
U8 user_data
Definition: lasdefinitions.hpp:488
LASoccupancyGrid::minus_minus_sizes
U16 * minus_minus_sizes
Definition: lasutility.hpp:159
LASbin::one_over_step
F32 one_over_step
Definition: lasutility.hpp:87
LAShistogram::scan_angle_bin_z
LASbin * scan_angle_bin_z
Definition: lasutility.hpp:127
LASwavepacket::setOffset
void setOffset(U64 offset)
Definition: lasdefinitions.hpp:78
LAShistogram::return_map_bin_intensity
LASbin * return_map_bin_intensity
Definition: lasutility.hpp:130
LASpoint::gps_time
F64 gps_time
Definition: lasdefinitions.hpp:497
LAShistogram::gps_time_bin
LASbin * gps_time_bin
Definition: lasutility.hpp:119
LAShistogram::scan_angle_bin_intensity
LASbin * scan_angle_bin_intensity
Definition: lasutility.hpp:129
LASpoint::scan_angle_rank
I8 scan_angle_rank
Definition: lasdefinitions.hpp:487
LASsummary::min
LASpoint min
Definition: lasutility.hpp:65
LASoccupancyGrid::occupied
BOOL occupied(const LASpoint *point) const
Definition: lasutility.cpp:772
LASpoint::scan_direction_flag
U8 scan_direction_flag
Definition: lasdefinitions.hpp:484
LASpoint
Definition: lasdefinitions.hpp:472
LASpoint::point_source_ID
U16 point_source_ID
Definition: lasdefinitions.hpp:489
LASwavepacket::setIndex
void setIndex(U8 index)
Definition: lasdefinitions.hpp:77
LASoccupancyGrid::reset
void reset()
Definition: lasutility.cpp:863
LASwavepacket::getZt
F32 getZt() const
Definition: lasdefinitions.hpp:76
LASwavepacket::setYt
void setYt(F32 yt)
Definition: lasdefinitions.hpp:82
LAShistogram::~LAShistogram
~LAShistogram()
Definition: lasutility.cpp:417
I64
long long I64
Definition: mydefs.hpp:48
LASbin::add_to_bin
void add_to_bin(I32 bin)
Definition: lasutility.cpp:208
F64
double F64
Definition: mydefs.hpp:52
LASinventory::LASinventory
LASinventory()
Definition: lasutility.cpp:37
lasutility.hpp
LASinventory::raw_max_x
I32 raw_max_x
Definition: lasutility.hpp:42
LASoccupancyGrid::write_asc_grid
BOOL write_asc_grid(const char *file_name) const
Definition: lasutility.cpp:910
I32
int I32
Definition: mydefs.hpp:35
LASinventory::add
BOOL add(const LASpoint *point)
Definition: lasutility.cpp:48
LASpoint::have_rgb
BOOL have_rgb
Definition: lasdefinitions.hpp:518
LASpoint::wavepacket
LASwavepacket wavepacket
Definition: lasdefinitions.hpp:499
TRUE
#define TRUE
Definition: mydefs.hpp:137
LASbin
Definition: lasutility.hpp:73
LASsummary::first
BOOL first
Definition: lasutility.hpp:70
LAShistogram::point_source_id_bin
LASbin * point_source_id_bin
Definition: lasutility.hpp:118
LASsummary::number_of_point_records
U32 number_of_point_records
Definition: lasutility.hpp:57
LAShistogram::wavepacket_index_bin
LASbin * wavepacket_index_bin
Definition: lasutility.hpp:120
LASwavepacket::getSize
U32 getSize() const
Definition: lasdefinitions.hpp:72
LASsummary::classification
U32 classification[32]
Definition: lasutility.hpp:61
LAShistogram::histo
BOOL histo(const char *name, F32 step)
Definition: lasutility.cpp:478
LASoccupancyGrid::num_occupied
U32 num_occupied
Definition: lasutility.hpp:170
LAShistogram::wavepacket_offset_bin
LASbin * wavepacket_offset_bin
Definition: lasutility.hpp:121
LASoccupancyGrid::active
BOOL active() const
Definition: lasutility.cpp:857
LASoccupancyGrid::anker
I32 anker
Definition: lasutility.hpp:155
I32_FLOOR
#define I32_FLOOR(n)
Definition: mydefs.hpp:118
LAShistogram::classification_bin_scan_angle
LASbin * classification_bin_scan_angle
Definition: lasutility.hpp:126
LAShistogram::z_bin
LASbin * z_bin
Definition: lasutility.hpp:114
LASpoint::y
I32 y
Definition: lasdefinitions.hpp:479
LASwavepacket::getYt
F32 getYt() const
Definition: lasdefinitions.hpp:75
LASpoint::classification
U8 classification
Definition: lasdefinitions.hpp:486
LASpoint::rgb
U16 rgb[4]
Definition: lasdefinitions.hpp:498
LASoccupancyGrid::max_y
I32 max_y
Definition: lasutility.hpp:151
LAShistogram::histo_avg
BOOL histo_avg(const char *name, F32 step, const char *name_avg)
Definition: lasutility.cpp:513
LASpoint::return_number
U8 return_number
Definition: lasdefinitions.hpp:482
LASwavepacket::setZt
void setZt(F32 zt)
Definition: lasdefinitions.hpp:83
LAShistogram::LAShistogram
LAShistogram()
Definition: lasutility.cpp:392
LASinventory::raw_max_z
I32 raw_max_z
Definition: lasutility.hpp:46
LAShistogram::add
void add(const LASpoint *point)
Definition: lasutility.cpp:560
LASoccupancyGrid::minus_minus_size
U32 minus_minus_size
Definition: lasutility.hpp:157
LASinventory::raw_min_z
I32 raw_min_z
Definition: lasutility.hpp:47
LASbin::count
I64 count
Definition: lasutility.hpp:86
LASsummary::max
LASpoint max
Definition: lasutility.hpp:66
LAShistogram::is_active
BOOL is_active
Definition: lasutility.hpp:110
LASoccupancyGrid::add_internal
BOOL add_internal(I32 pos_x, I32 pos_y)
Definition: lasutility.cpp:652
LASwavepacket::getXt
F32 getXt() const
Definition: lasdefinitions.hpp:74
LASpoint::get_y
F64 get_y() const
Definition: lasdefinitions.hpp:812
LASwavepacket::setSize
void setSize(U32 size)
Definition: lasdefinitions.hpp:79
LASbin::add
void add(I32 item)
Definition: lasutility.cpp:184
LASbin::size_pos
I32 size_pos
Definition: lasutility.hpp:90
LASoccupancyGrid::plus_plus_size
U32 plus_plus_size
Definition: lasutility.hpp:167
LASoccupancyGrid::plus_minus
U32 ** plus_minus
Definition: lasutility.hpp:165
LASsummary::classification_keypoint
U32 classification_keypoint
Definition: lasutility.hpp:63
LAShistogram::x_bin
LASbin * x_bin
Definition: lasutility.hpp:112
U16
unsigned short U16
Definition: mydefs.hpp:40
LASsummary::classification_withheld
U32 classification_withheld
Definition: lasutility.hpp:64
LASbin::report
void report(FILE *file, const char *name=0, const char *name_avg=0) const
Definition: lasutility.cpp:323
LASoccupancyGrid::max_x
I32 max_x
Definition: lasutility.hpp:151
LASoccupancyGrid::~LASoccupancyGrid
~LASoccupancyGrid()
Definition: lasutility.cpp:962
LASpoint::have_gps_time
BOOL have_gps_time
Definition: lasdefinitions.hpp:517
LAShistogram::wavepacket_size_bin
LASbin * wavepacket_size_bin
Definition: lasutility.hpp:122
LASsummary::LASsummary
LASsummary()
Definition: lasutility.cpp:71
BOOL
int BOOL
Definition: mydefs.hpp:57
LASpoint::z
I32 z
Definition: lasdefinitions.hpp:480
LASbin::bins_pos
U32 * bins_pos
Definition: lasutility.hpp:92
LASoccupancyGrid::minus_plus_sizes
U16 * minus_plus_sizes
Definition: lasutility.hpp:162
LASwavepacket::setLocation
void setLocation(F32 location)
Definition: lasdefinitions.hpp:80
argc
int argc
Definition: tests_high_five_parallel.cpp:27
LAShistogram::wavepacket_location_bin
LASbin * wavepacket_location_bin
Definition: lasutility.hpp:123
LAShistogram::y_bin
LASbin * y_bin
Definition: lasutility.hpp:113
LASbin::values_pos
F64 * values_pos
Definition: lasutility.hpp:94
FALSE
#define FALSE
Definition: mydefs.hpp:133
file
FILE * file
Definition: arithmeticencoder.cpp:77
LASoccupancyGrid::LASoccupancyGrid
LASoccupancyGrid(F32 grid_spacing)
Definition: lasutility.cpp:941
F32
float F32
Definition: mydefs.hpp:51
LASoccupancyGrid::min_x
I32 min_x
Definition: lasutility.hpp:151
LASsummary::add
BOOL add(const LASpoint *point)
Definition: lasutility.cpp:84
LASbin::total
F64 total
Definition: lasutility.hpp:85
LASpoint::have_wavepacket
BOOL have_wavepacket
Definition: lasdefinitions.hpp:519
LASoccupancyGrid::plus_plus
U32 ** plus_plus
Definition: lasutility.hpp:168
LASbin::values_neg
F64 * values_neg
Definition: lasutility.hpp:95
LASwavepacket::getLocation
F32 getLocation() const
Definition: lasdefinitions.hpp:73
LASoccupancyGrid::plus_ankers
I32 * plus_ankers
Definition: lasutility.hpp:163
LASbin::~LASbin
~LASbin()
Definition: lasutility.cpp:176
LASoccupancyGrid::min_y
I32 min_y
Definition: lasutility.hpp:151
LASoccupancyGrid::plus_minus_sizes
U16 * plus_minus_sizes
Definition: lasutility.hpp:166
LASpoint::number_of_returns_of_given_pulse
U8 number_of_returns_of_given_pulse
Definition: lasdefinitions.hpp:483
LAShistogram::scan_angle_bin_number_of_returns
LASbin * scan_angle_bin_number_of_returns
Definition: lasutility.hpp:128
LASoccupancyGrid::minus_ankers
I32 * minus_ankers
Definition: lasutility.hpp:156
LASbin::anker
I32 anker
Definition: lasutility.hpp:89
LASinventory::first
BOOL first
Definition: lasutility.hpp:51
LASoccupancyGrid::grid_spacing
F32 grid_spacing
Definition: lasutility.hpp:154
LASsummary::classification_synthetic
U32 classification_synthetic
Definition: lasutility.hpp:62
LASinventory::number_of_point_records
U32 number_of_point_records
Definition: lasutility.hpp:39
U32
unsigned int U32
Definition: mydefs.hpp:39
LASoccupancyGrid::add
BOOL add(const LASpoint *point)
Definition: lasutility.cpp:613
LASbin::LASbin
LASbin(F32 step)
Definition: lasutility.cpp:162
LASwavepacket::getIndex
U8 getIndex() const
Definition: lasdefinitions.hpp:70
LASbin::size_neg
I32 size_neg
Definition: lasutility.hpp:91
LASoccupancyGrid::minus_plus
U32 ** minus_plus
Definition: lasutility.hpp:161
LASpoint::x
I32 x
Definition: lasdefinitions.hpp:478
LAShistogram::report
void report(FILE *file) const
Definition: lasutility.cpp:589
LAShistogram::scan_angle_bin
LASbin * scan_angle_bin
Definition: lasutility.hpp:117
LASpoint::get_x
F64 get_x() const
Definition: lasdefinitions.hpp:811
LASwavepacket::setXt
void setXt(F32 xt)
Definition: lasdefinitions.hpp:81
LASsummary::number_of_points_by_return
U32 number_of_points_by_return[8]
Definition: lasutility.hpp:59
LAShistogram::parse
BOOL parse(int argc, char *argv[])
Definition: lasutility.cpp:441
LASbin::first
BOOL first
Definition: lasutility.hpp:88
LASoccupancyGrid::minus_plus_size
U32 minus_plus_size
Definition: lasutility.hpp:160
argv
char ** argv
Definition: tests_high_five_parallel.cpp:28
LASoccupancyGrid::minus_minus
U32 ** minus_minus
Definition: lasutility.hpp:158
LAShistogram::classification_bin
LASbin * classification_bin
Definition: lasutility.hpp:116
LASinventory::raw_max_y
I32 raw_max_y
Definition: lasutility.hpp:44
scripts.create_png.step
step
Definition: create_png.py:42
LASpoint::intensity
U16 intensity
Definition: lasdefinitions.hpp:481
LASinventory::raw_min_x
I32 raw_min_x
Definition: lasutility.hpp:43
LASinventory::raw_min_y
I32 raw_min_y
Definition: lasutility.hpp:45
LASwavepacket::getOffset
U64 getOffset() const
Definition: lasdefinitions.hpp:71
LASoccupancyGrid::plus_minus_size
U32 plus_minus_size
Definition: lasutility.hpp:164
LASoccupancyGrid::plus_plus_sizes
U16 * plus_plus_sizes
Definition: lasutility.hpp:169
LASinventory::number_of_points_by_return
U32 number_of_points_by_return[8]
Definition: lasutility.hpp:41


lvr2
Author(s): Thomas Wiemann , Sebastian Pütz , Alexander Mock , Lars Kiesow , Lukas Kalbertodt , Tristan Igelbrink , Johan M. von Behren , Dominik Feldschnieders , Alexander Löhr
autogenerated on Wed Mar 2 2022 00:37:23