gdiam.cpp
Go to the documentation of this file.
1 /*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
2  * gdiam.C -
3  * Computes diameter, and a tight-fitting bounding box of a
4  * point-set in 3d.
5  *
6  * Copyright 2000 Sariel Har-Peled (ssaarriieell@cs.uiuc.edu)
7  *
8  * This program is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU General Public License as
11  * or (at your option) any later version.
12  *
13  * Code is based on the paper:
14  * A Practical Approach for Computing the Diameter of a
15  * Point-Set (in ACM Sym. on Computation Geometry
16  * SoCG 2001)
17  * Sariel Har-Peled (http://www.uiuc.edu/~sariel)
18  *--------------------------------------------------------------
19  * History
20  * 3/28/01 -
21  * This is a more robust version of the code. It should
22  * handlereally abnoxious inputs well (i.e., points with equal
23  * coordinates, etc.
24 \*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*/
25
26 #include <stdlib.h>
27 #include <stdio.h>
28 #include <assert.h>
29 #include <memory.h>
30 #include <math.h>
31
32 #include <vector>
33 #include <algorithm>
34
35 #include "gdiam.h"
36
37 /*--- Constants ---*/
38
39 #define HEAP_LIMIT 10000
40 #define HEAP_DELTA 10000
41 #define PI 3.14159265958979323846
42
43 /*--- Start of Code ---*/
44
45 typedef long double ldouble;
46
47 class GFSPPair;
48
49 /*void computeExtremePairDir( const ArrPoint3d & arr,
50  const Point3d & dir,
51  GPointPair & pp,
52  int start,
53  int end
54  );
55 */
56
58 {
59 private:
60  GBBox bbox;
61  //int left_ind, p_pnt_right; // range in the array
62  gdiam_point* p_pnt_left, * p_pnt_right;
63  gdiam_real bbox_diam, bbox_diam_proj;
64  GFSPTreeNode* left, * right;
65  gdiam_point_t center;
66
67 public:
68  void dump() const
69  {
70  printf("dump...\n");
71  //printf( "\tNode: Range: [%d, %d]\n", left_ind, p_pnt_right );
72  }
73  int size() const
74  {
75  return (int)(p_pnt_right - p_pnt_left);
76  }
78  {
79  return (gdiam_point)center;
80  }
81  int nodes_number() const
82  {
83  int sum;
84
85  if ((left == NULL) && (right == NULL))
86  {
87  return 1;
88  }
89
90  sum = 1;
91
92  if (left != NULL)
93  {
94  sum += left->nodes_number();
95  }
96
97  if (right != NULL)
98  {
99  sum += right->nodes_number();
100  }
101
102  return sum;
103  }
104
105
106  GFSPTreeNode(gdiam_point* _p_pnt_left,
107  gdiam_point* _p_pnt_right)
108  {
109  bbox.init();
110  left = right = NULL;
111  p_pnt_left = _p_pnt_left;
112  p_pnt_right = _p_pnt_right;
113  bbox_diam = bbox_diam_proj = -1;
114  }
115
117  {
118  return bbox_diam;
119  }
121  {
122  return bbox_diam_proj;
123  }
124
126  {
127  return left;
128  }
130  {
131  return right;
132  }
133
134  bool isSplitable() const
135  {
136  return (size() > 1);
137  }
138
140  {
141  return p_pnt_left;
142  }
144  {
145  return p_pnt_left + (rand() % (p_pnt_right - p_pnt_left + 1));
146  }
148  {
149  return p_pnt_right;
150  }
151
152
153  friend class GFSPTree;
154
155  const GBBox& getBBox() const
156  {
157  return bbox;
158  }
159 };
160
161
162 class GFSPTree
163 {
164 private:
165  gdiam_point* arr;
166  GFSPTreeNode* root;
167
168 public:
169  void init(const gdiam_point* _arr,
170  int size)
171  {
172  arr = (gdiam_point*)malloc(sizeof(gdiam_point) * size);
173  assert(arr != NULL);
174
175  for (int ind = 0; ind < size; ind++)
176  {
177  arr[ ind ] = _arr[ ind ];
178  }
179
180  root = build_node(arr, arr + size - 1);
181  }
182
183  static void terminate(GFSPTreeNode* node)
184  {
185  if (node == NULL)
186  {
187  return;
188  }
189
190  terminate(node->left);
191  terminate(node->right);
192
193  node->left = node->right = NULL;
194
195  delete node;
196  }
197  void term()
198  {
199  free(arr);
200  arr = NULL;
201
202  GFSPTree::terminate(root);
203  root = NULL;
204  }
205
206  const gdiam_point* getPoints() const
207  {
208  return arr;
209  }
210
211  int nodes_number() const
212  {
213  return root->nodes_number();
214  }
215
216  gdiam_real brute_diameter(int a_lo, int a_hi,
217  int b_lo, int b_hi,
218  GPointPair& diam) const
219  {
220  for (int ind = a_lo; ind <= a_hi; ind++)
221  for (int jnd = b_lo; jnd <= b_hi; jnd++)
222  {
223  diam.update_diam(arr[ ind ], arr[ jnd ]);
224  }
225
226  return diam.distance;
227  }
228  gdiam_real brute_diameter(int a_lo, int a_hi,
229  int b_lo, int b_hi,
230  GPointPair& diam,
231  const gdiam_point dir) const
232  {
233  for (int ind = a_lo; ind <= a_hi; ind++)
234  for (int jnd = b_lo; jnd <= b_hi; jnd++)
235  {
236  diam.update_diam(arr[ ind ], arr[ jnd ], dir);
237  }
238
239  return diam.distance;
240  }
242  const gdiam_point* a_hi,
243  const gdiam_point* b_lo,
244  const gdiam_point* b_hi,
245  GPointPair& diam) const
246  {
247  for (const gdiam_point* ind = a_lo; ind <= a_hi; ind++)
248  for (const gdiam_point* jnd = b_lo; jnd <= b_hi; jnd++)
249  {
250  diam.update_diam(*ind, *jnd);
251  }
252
253  return diam.distance;
254  }
256  const gdiam_point* a_hi,
257  const gdiam_point* b_lo,
258  const gdiam_point* b_hi,
259  GPointPair& diam,
260  const gdiam_point dir) const
261  {
262  for (const gdiam_point* ind = a_lo; ind <= a_hi; ind++)
263  for (const gdiam_point* jnd = b_lo; jnd <= b_hi; jnd++)
264  {
265  diam.update_diam(*ind, *jnd, dir);
266  }
267
268  return diam.distance;
269  }
270
271  gdiam_point getPoint(int ind) const
272  {
273  return arr[ ind ];
274  }
275
277  {
278  return root;
279  }
280
282  gdiam_point* right)
283  {
284  if (left > right)
285  {
286  printf("what!?\n");
287  fflush(stdout);
288  assert(left <= right);
289  }
290
291  while ((right > left)
292  && (pnt_isEqual(*right, *left)))
293  {
294  right--;
295  }
296
297  GFSPTreeNode* node = new GFSPTreeNode(left, right);
298
299  node->bbox.init();
300
301  for (const gdiam_point* ind = left; ind <= right; ind++)
302  {
303  node->bbox.bound(*ind);
304  }
305
306  node->bbox_diam = node->bbox.get_diam();
307
308  node->bbox.center(node->center);
309
310  return node;
311  }
312
313  // return how many elements on left side.
315  gdiam_point* right,
316  int dim, const gdiam_real& val)
317  {
318  const gdiam_point* start_left = left;
319
320  while (left < right)
321  {
322  if ((*left)[ dim ] < val)
323  {
324  left++;
325  }
326  else if ((*right)[ dim ] >= val)
327  {
328  right--;
329  }
330  else
331  {
332  gdiam_exchange(*right, *left);
333  }
334  }
335
336  return left - start_left;
337  }
338
339
341  {
342  int dim, left_size;
343  gdiam_real split_val;
344
345  if (node->left != NULL)
346  {
347  return;
348  }
349
350  GBBox& bb(node->bbox);
351
352  dim = bb.getLongestDim();
353  split_val = (bb.min_coord(dim) + bb.max_coord(dim)) / 2.0;
354
355  left_size = split_array(node->p_pnt_left, node->p_pnt_right,
356  dim, split_val);
357
358  if (left_size <= 0.0)
359  {
360  printf("bb: %g %g\n",
361  bb.min_coord(dim), bb.max_coord(dim));
362  printf("left: %p, right: %p\n",
363  node->p_pnt_left, node->p_pnt_right);
364  assert(left_size > 0);
365  }
366
367  if (left_size >= (node->p_pnt_right - node->p_pnt_left + 1))
368  {
369  printf("left size too large?\n");
370  fflush(stdout);
371  assert(left_size < (node->p_pnt_right - node->p_pnt_left + 1));
372  }
373
374  node->left = build_node(node->p_pnt_left,
375  node->p_pnt_left + left_size - 1);
376  node->right = build_node(node->p_pnt_left + left_size,
377  node->p_pnt_right);
378  }
379
381  GPointPair& max_pair_diam);
382 };
383
384 void bbox_vertex(const GBBox& bb, gdiam_point_t& ver,
385  int vert)
386 {
387  ver[ 0 ] = ((vert & 0x1) != 0) ?
388  bb.min_coord(0) : bb.max_coord(0);
389  ver[ 1 ] = ((vert & 0x2) != 0) ?
390  bb.min_coord(1) : bb.max_coord(1);
391  ver[ 2 ] = ((vert & 0x4) != 0) ?
392  bb.min_coord(2) : bb.max_coord(2);
393 }
394
395
397  const GBBox& bb2,
398  gdiam_point_cnt proj)
399 {
400  gdiam_point_t a, b;
401  gdiam_real rl;
402
403  rl = 0;
404
405  for (int ind = 0; ind < 8; ind++)
406  {
407  bbox_vertex(bb1, a, ind);
408
409  //printf( "\n\n" );
410  for (int jnd = 0; jnd < 8; jnd++)
411  {
412  bbox_vertex(bb2, b, jnd);
413  //pnt_dump( b );
414  rl = max(rl, pnt_distance(a, b, proj));
415  }
416  }
417
418  return rl;
419 }
420
421
422
423 class GFSPPair
424 {
425 private:
426  GFSPTreeNode* left, * right;
427  double max_diam;
428
429 public:
430  void init(GFSPTreeNode* _left, GFSPTreeNode* _right)
431  {
432  left = _left;
433  right = _right;
434
435  GBBox bb;
436
437  bb.init(left->getBBox(), right->getBBox());
438  //printf( "\t\tbbox: \n" );
439  //bb.dump();
440  max_diam = bb.get_diam();
441  //printf( "\t\tmax_diam: %g\n", max_diam );
442  }
443
444  void init(GFSPTreeNode* _left, GFSPTreeNode* _right,
445  gdiam_point proj_dir, gdiam_real dist)
446  {
447  left = _left;
448  right = _right;
449
450  GBBox bb;
451
452  bb.init(left->getBBox(), right->getBBox());
453  //printf( "\t\tbbox: \n" );
454  //bb.dump();
455
456  //max_diam = dist + _left->getBBox().get_diam()
457  // + _right->getBBox().get_diam();
458  max_diam = max(max(bbox_proj_dist(_left->getBBox(),
459  _right->getBBox(),
460  proj_dir),
461  bbox_proj_dist(_left->getBBox(),
462  _left->getBBox(),
463  proj_dir)),
464  bbox_proj_dist(_right->getBBox(),
465  _right->getBBox(),
466  proj_dir));
467
468  //printf( "diam: %g, smart diam: %g\n",
469  // max_diam,
470  // bbox_proj_dist( _left->getBBox(),
471  // _right->getBBox(),
472  // proj_dir ) );
473
474  //printf( "dist: %g, %g, %g:\n",
475  // dist, _left->getBBox().get_diam(),
476  // _right->getBBox().get_diam() );
477
478  //bb.get_diam_proj( proj_dir );
479  //printf( "\t\tmax_diam: %g\n", max_diam );
480  }
481
482  // error 2r/l - approximate cos of the max possible angle in the pair
483  double getProjectionError() const
484  {
485  double l, two_r;
486
487  l = pnt_distance(left->getCenter(), right->getCenter());
488  two_r = max(left->maxDiam(), right->maxDiam());
489
490  if (l == 0.0)
491  {
492  return 10;
493  }
494
495  return two_r / l;
496  }
497
499  {
500  return left;
501  }
503  {
504  return right;
505  }
506
507  void dump() const
508  {
509  printf("[\n");
510  left->dump();
511  right->dump();
512  printf("\tmax_diam; %g\n", max_diam);
513  }
514
515  bool isBelowThreshold(int threshold) const
516  {
517  if (left->size() > threshold)
518  {
519  return false;
520  }
521
522  if (right->size() > threshold)
523  {
524  return false;
525  }
526
527  return true;
528  }
529
530  double directDiameter(const GFSPTree& tree,
531  GPointPair& diam) const
532  {
533  return tree.brute_diameter(left->ptr_pnt_left(),
534  left->ptr_pnt_right(),
535  right->ptr_pnt_left(),
536  right->ptr_pnt_right(),
537  diam);
538  }
539
540  double directDiameter(const GFSPTree& tree,
541  GPointPair& diam,
542  const gdiam_point dir) const
543  {
544  return tree.brute_diameter(left->ptr_pnt_left(),
545  left->ptr_pnt_right(),
546  right->ptr_pnt_left(),
547  right->ptr_pnt_right(),
548  diam, dir);
549  }
550
551  double maxDiam() const
552  {
553  return max_diam;
554  }
555  double minAprxDiam() const
556  {
557  double sub_diam;
558
559  sub_diam = left->maxDiam() + right->maxDiam();
560
561  if (sub_diam > (max_diam / 10.0))
562  {
563  return max_diam;
564  }
565
566  return max_diam - sub_diam / 2.0;
567  }
568 };
569
570
571 class g_heap_pairs_p;
572
573 #define UDM_SIMPLE 1
574 #define UDM_BIG_SAMPLE 2
575
577 {
578 private:
579  GFSPTree tree;
580  //double diam;
581  GPointPair pair_diam;
582  int points_num;
583  int threshold_brute;
584  int update_diam_mode;
585 public:
586  GTreeDiamAlg(gdiam_point* arr, int size, int mode)
587  {
588  tree.init(arr, size);
589  //diam = 0;
590  pair_diam.init(arr[ 0 ]);
591  points_num = size;
592  threshold_brute = 40;
593  update_diam_mode = mode;
594  }
595
596  void term()
597  {
598  tree.term();
599  }
600
601  int size() const
602  {
603  return points_num;
604  }
605  const GPointPair& getDiameter() const
606  {
607  return pair_diam;
608  }
609
610  int nodes_number() const
611  {
612  return tree.nodes_number();
613  }
614
616  GFSPTreeNode* left,
617  GFSPTreeNode* right);
619  GFSPTreeNode* left,
620  GFSPTreeNode* right,
621  gdiam_point proj,
622  GFSPPair& father);
623  void split_pair(GFSPPair& pair,
624  g_heap_pairs_p& heap, double eps);
625  void split_pair_proj(GFSPPair& pair,
626  g_heap_pairs_p& heap, double eps,
627  gdiam_point proj);
628  void compute_by_heap(double eps);
629  void compute_by_heap_proj(double eps, gdiam_point proj);
630
631  double diameter()
632  {
633  return pair_diam.distance;
634  }
635 };
636
637
638 /*--- Heap implementation ---*/
639
640 typedef int (*ptrCompareFunc)(void* aPtr, void* bPtr);
641 typedef void* voidPtr_t;
642
643 struct heap_t
644 {
648 };
649
650 void heap_init(heap_t* pHeap, ptrCompareFunc _pCompFunc)
651 {
652  assert(pHeap != NULL);
653  assert(_pCompFunc != NULL);
654
655  memset(pHeap, 0, sizeof(heap_t));
656
657  pHeap->pCompFunc = _pCompFunc;
658  pHeap->max_size = 100;
659  pHeap->pArr = (voidPtr_t*)malloc(sizeof(void*)
660  * pHeap->max_size);
661  assert(pHeap->pArr != NULL);
662  pHeap->curr_size = 0;
663 }
664
665
666 static void resize(heap_t* pHeap, int size)
667 {
668  int max_sz;
669  voidPtr_t* pTmp;
670
671  if (size <= pHeap->max_size)
672  {
673  return;
674  }
675
676  max_sz = size * 2;
677  pTmp = (voidPtr_t*)malloc(max_sz * sizeof(void*));
678  assert(pTmp != NULL);
679  memset(pTmp, 0, max_sz * sizeof(void*));
680  memcpy(pTmp, pHeap->pArr, pHeap->curr_size * sizeof(void*));
681  free(pHeap->pArr);
682  pHeap->pArr = pTmp;
683 }
684
685
686 void heap_term(heap_t* pHeap)
687 {
688  assert(pHeap != NULL);
689
690  if (pHeap->pArr != NULL)
691  {
692  free(pHeap->pArr);
693  }
694
695  memset(pHeap, 0, sizeof(heap_t));
696 }
697
698
699 void heap_insert(heap_t* pHeap, void* pData)
700 {
701  int ind, father;
702  voidPtr_t pTmp;
703
704  resize(pHeap, pHeap->curr_size + 1);
705
706  ind = pHeap->curr_size;
707  pHeap->curr_size++;
708
709  pHeap->pArr[ ind ] = pData;
710
711  while (ind > 0)
712  {
713  father = (ind - 1) / 2;
714
715  if ((*pHeap->pCompFunc)(pHeap->pArr[ ind ],
716  pHeap->pArr[ father ]) <= 0)
717  {
718  break;
719  }
720
721  pTmp = pHeap->pArr[ ind ];
722  pHeap->pArr[ ind ] = pHeap->pArr[ father ];
723  pHeap->pArr[ father ] = pTmp;
724
725  ind = father;
726  }
727 }
728
729
730 void* heap_top(heap_t* pHeap)
731 {
732  return pHeap->pArr[ 0 ];
733 }
734
735
736 void* heap_delete_max(heap_t* pHeap)
737 {
738  void* res;
739  void* pTmp;
740  int ind, left, right, max_ind;
741
742  if (pHeap->curr_size <= 0)
743  {
744  return NULL;
745  }
746
747  res = pHeap->pArr[ 0 ];
748
749  //printf( "pHeap->curr_size: %d\n", pHeap->curr_size );
750  //printf( "res= %p\n", res );
751
752  pHeap->curr_size--;
753  pHeap->pArr[ 0 ] = pHeap->pArr[ pHeap->curr_size ];
754  pHeap->pArr[ pHeap->curr_size ] = NULL;
755  ind = 0;
756
757  while (ind < pHeap->curr_size)
758  {
759  left = 2 * ind + 1;
760  right = 2 * ind + 2;
761
762  if (left >= pHeap->curr_size)
763  {
764  break;
765  }
766
767  if (right >= pHeap->curr_size)
768  {
769  right = left;
770  }
771
772  if ((*pHeap->pCompFunc)(pHeap->pArr[ left ],
773  pHeap->pArr[ right ]) <= 0)
774  {
775  max_ind = right;
776  }
777  else
778  {
779  max_ind = left;
780  }
781
782  if ((*pHeap->pCompFunc)(pHeap->pArr[ ind ],
783  pHeap->pArr[ max_ind ]) >= 0)
784  {
785  break;
786  }
787
788  pTmp = pHeap->pArr[ ind ];
789  pHeap->pArr[ ind ] = pHeap->pArr[ max_ind ];
790  pHeap->pArr[ max_ind ] = pTmp;
791
792  ind = max_ind;
793  }
794
795  return res;
796 }
797
798
799 bool heap_is_empty(heap_t* pHeap)
800 {
801  assert(pHeap != NULL);
802
803  return (pHeap->curr_size <= 0);
804 }
805
806
807 int compare_pairs(void* _a, void* _b)
808 {
809  const GFSPPair& a(*(GFSPPair*)_a);
810  const GFSPPair& b(*(GFSPPair*)_b);
811
812  if (a.maxDiam() < b.maxDiam())
813  {
814  return -1;
815  }
816
817  if (a.maxDiam() > b.maxDiam())
818  {
819  return 1;
820  }
821
822  return 0;
823 }
824
825
826
828 {
829 private:
830  heap_t heap;
831
832 public:
834  {
835  heap_init(&heap, compare_pairs);
836  }
837
839  {
840  while (size() > 0)
841  {
842  pop();
843  }
844
845  heap_term(&heap);
846  }
847
848  void push(GFSPPair& pair)
849  {
850  //printf( "pushing: (%p, %p)\n", pair.get_left(),
851  // pair.get_right() );
852  GFSPPair* p_pair = new GFSPPair(pair);
853  heap_insert(&heap, (void*)p_pair);
854  }
855
856  int size() const
857  {
858  return heap.curr_size;
859  }
861  {
862  return *(GFSPPair*)heap_top(&heap);
863  }
864  void pop()
865  {
866  GFSPPair* ptr = (GFSPPair*)heap_delete_max(&heap);
867
868  delete ptr;
869  }
870 };
871
872
873 /* heap.implementation end */
874
875
876 void computeExtremePair(const gdiam_point* arr, int size,
877  int dim, GPointPair& pp)
878 {
879  pp.init(arr[ 0 ]);
880
881  for (int ind = 1; ind < size; ind++)
882  {
883  const gdiam_point pnt(arr[ ind ]);
884
885  if (pnt[ dim ] < pp.p[ dim ])
886  {
887  pp.p = pnt;
888  }
889
890  if (pnt[ dim ] > pp.q[ dim ])
891  {
892  pp.q = pnt;
893  }
894  }
895
896  pp.distance = pnt_distance(pp.p, pp.q);
897 }
898
899
900
902 {
903  g_heap_pairs_p heap;
904  int heap_limit, heap_delta;
905
906  heap_limit = HEAP_LIMIT;
907  heap_delta = HEAP_DELTA;
908
909  GFSPTreeNode* root = tree.getRoot();
910
911  computeExtremePair(tree.getPoints(), points_num,
912  root->getBBox().getLongestDim(),
913  pair_diam);
914
915  GFSPPair root_pair;
916
917  root_pair.init(root, root);
918
919  heap.push(root_pair);
920  //queue_a.pushx( root_pair );
921
922  int count = 0;
923
924  while (heap.size() > 0)
925  {
926  //printf( "heap.size: %d\n", heap.size() );
927  //printf( "%d: curr_diam: %g\n", count++, pair_diam.distance );
928  GFSPPair pair = heap.top();
929  heap.pop();
930  split_pair(pair, heap, eps);
931
932  //if ( ( count & 0x3ff ) == 0 ) {
933  // printf( "heap.size: %d\n", heap.size() );
934  //}
935  if ((count & 0x3ff) == 0)
936  {
937  if (((int)heap.size()) > heap_limit)
938  {
939  threshold_brute *= 2;
940  printf("threshold_brute: %d\n", threshold_brute);
941  heap_limit += heap_delta;
942  }
943  }
944
945  count++;
946  }
947 }
948
949
951  gdiam_point proj)
952 {
953  g_heap_pairs_p heap;
954  int heap_limit, heap_delta;
955  GPointPair pair_diam_x;
956
957
958  //pnt_dump( proj );
959  //printf( "length = %g\n", pnt_length( proj ) );
960
961  heap_limit = HEAP_LIMIT;
962  heap_delta = HEAP_DELTA;
963
964  GFSPTreeNode* root = tree.getRoot();
965
966  computeExtremePair(tree.getPoints(), points_num,
967  root->getBBox().getLongestDim(),
968  pair_diam_x);
969  pair_diam.init(pair_diam_x.p, pair_diam_x.q, proj);
970
971  GFSPPair root_pair;
972
973  root_pair.init(root, root, proj, 0);
974
975  //printf( "root pair distance: %g\n", root_pair.maxDiam() );
976  heap.push(root_pair);
977  //queue_a.pushx( root_pair );
978
979  int count = 0;
980
981  while (heap.size() > 0)
982  {
983  //printf( "heap.size: %d\n", heap.size() );
984  //printf( "%d: curr_diam: %g\n", count++, pair_diam.distance );
985  GFSPPair pair = heap.top();
986  heap.pop();
987  //printf( "pair distance: %g\n", pair.maxDiam() );
988  split_pair_proj(pair, heap, eps, proj);
989
990  //if ( ( count & 0x3ff ) == 0 ) {
991  // printf( "heap.size: %d\n", heap.size() );
992  //}
993  if ((count & 0x3ff) == 0)
994  {
995  if (((int)heap.size()) > heap_limit)
996  {
997  threshold_brute *= 2;
998  printf("threshold_brute: %d\n", threshold_brute);
999  heap_limit += heap_delta;
1000  }
1001  }
1002
1003  count++;
1004  //printf( "Poping!\n" );
1005  }
1006 }
1007
1008
1010  GFSPTreeNode* left,
1011  GFSPTreeNode* right)
1012 {
1013  if (update_diam_mode == UDM_SIMPLE)
1014  {
1015  const gdiam_point p(*(left->ptr_pnt_left()));
1016  const gdiam_point q(*(right->ptr_pnt_left()));
1017
1018  pair_diam.update_diam(p, q);
1019  }
1020  else if (update_diam_mode == UDM_BIG_SAMPLE)
1021  {
1022  if ((left->size() < 100) || (right->size() < 100))
1023  {
1024  const gdiam_point p(*(left->ptr_pnt_left()));
1025  const gdiam_point q(*(right->ptr_pnt_left()));
1026
1027  pair_diam.update_diam(p, q);
1028  }
1029  else
1030  {
1031 #define UMD_SIZE 5
1032  gdiam_point arr_left[ UMD_SIZE ],
1033  arr_right[ UMD_SIZE ];
1034
1035  for (int ind = 0; ind < UMD_SIZE; ind++)
1036  {
1037  const gdiam_point p(*(left->ptr_pnt_rand()));
1038  const gdiam_point q(*(right->ptr_pnt_rand()));
1039  arr_left[ ind ] = p;
1040  arr_right[ ind ] = q;
1041
1042  //pair_diam.update_diam( p, q );
1043  }
1044
1045  for (int ind = 0; ind < UMD_SIZE; ind++)
1046  for (int jnd = 0; jnd < UMD_SIZE; jnd++)
1047  pair_diam.update_diam(arr_left[ ind ],
1048  arr_right[ jnd ]);
1049  }
1050  }
1051  else
1052  {
1053  assert(false);
1054  }
1055
1056  //printf( "old_diam; %g\n", diam );
1057  //diam = max( diam, pnt_distance( p, q ) );
1058  //printf( "new_diam; %g\n", diam );
1059
1060  GFSPPair pair;
1061
1062  pair.init(left, right);
1063
1064  if (pair.maxDiam() <= pair_diam.distance)
1065  {
1066  return;
1067  }
1068
1069  heap.push(pair);
1070 }
1071
1072
1074  GFSPTreeNode* left,
1075  GFSPTreeNode* right,
1076  gdiam_point proj,
1077  GFSPPair& father)
1078 {
1079  const gdiam_point p(*(left->ptr_pnt_left()));
1080  const gdiam_point q(*(right->ptr_pnt_left()));
1081
1082  //printf( "addPairHeap( ?, %p, %p, ? )\n", left, right );
1083  pair_diam.update_diam(p, q, proj);
1084  //printf( "old_diam; %g\n", diam );
1085  //diam = max( diam, pnt_distance( p, q ) );
1086  //printf( "new_diam; %g\n", diam );
1087
1088  GFSPPair pair;
1089
1090  pair.init(left, right, proj,
1091  pnt_distance(p, q, proj));
1092
1093  if (pair.maxDiam() <= pair_diam.distance)
1094  {
1095  return;
1096  }
1097
1098  //printf( "pair.maxDiam: %g, father.maxDiam: %g\n",
1099  // pair.maxDiam(), father.maxDiam() );
1100  //printf( "pair_diam.distance: %g\n", pair_diam.distance );
1101  //pnt_dump( proj );
1102  heap.push(pair);
1103 }
1104
1105
1107  g_heap_pairs_p& heap, double eps,
1108  gdiam_point proj)
1109 {
1110  bool f_is_left_splitable, f_is_right_splitable;
1111
1112  //printf( "pair.maxDiam: %g, limit: %g\n",
1113  // pair.maxDiam(), ((1.0 + eps) * pair_diam.distance ));
1114
1115  if (pair.maxDiam() <= ((1.0 + eps) * pair_diam.distance))
1116  {
1117  return;
1118  }
1119
1120  if (pair.isBelowThreshold(threshold_brute))
1121  {
1122  pair.directDiameter(tree, pair_diam, proj);
1123  return;
1124  }
1125
1126  //printf( "Curr pair: (%p, %p)\n", pair.get_left(), pair.get_right() );
1127  f_is_left_splitable = pair.get_left()->isSplitable();
1128  f_is_right_splitable = pair.get_right()->isSplitable();
1129  assert(f_is_left_splitable || f_is_right_splitable);
1130
1131  if (f_is_left_splitable)
1132  {
1133  tree.split_node(pair.get_left());
1134  }
1135
1136  if (f_is_right_splitable)
1137  {
1138  tree.split_node(pair.get_right());
1139  }
1140
1141  if (f_is_left_splitable
1142  && f_is_right_splitable)
1143  {
1145  pair.get_left()->get_left(),
1146  pair.get_right()->get_left(),
1147  proj, pair);
1149  pair.get_left()->get_left(),
1150  pair.get_right()->get_right(),
1151  proj, pair);
1152
1153  // to avoid exponential blowup
1154  if (pair.get_left() != pair.get_right())
1156  pair.get_left()->get_right(),
1157  pair.get_right()->get_left(),
1158  proj, pair);
1159
1161  pair.get_left()->get_right(),
1162  pair.get_right()->get_right(),
1163  proj, pair);
1164  return;
1165  }
1166
1167  if (f_is_left_splitable)
1168  {
1170  pair.get_left()->get_left(),
1171  pair.get_right(), proj, pair);
1173  pair.get_left()->get_right(),
1174  pair.get_right(), proj, pair);
1175  return;
1176  }
1177
1178  if (f_is_right_splitable)
1179  {
1181  pair.get_left(),
1182  pair.get_right()->get_left(), proj, pair);
1184  pair.get_left(),
1185  pair.get_right()->get_right(), proj, pair);
1186  return;
1187  }
1188 }
1189
1190
1192  g_heap_pairs_p& heap,
1193  double eps)
1194 {
1195  bool f_is_left_splitable, f_is_right_splitable;
1196
1197  if (pair.maxDiam() <= ((1.0 + eps) * pair_diam.distance))
1198  {
1199  return;
1200  }
1201
1202  if (pair.isBelowThreshold(threshold_brute))
1203  {
1204  pair.directDiameter(tree, pair_diam);
1205  return;
1206  }
1207
1208  f_is_left_splitable = pair.get_left()->isSplitable();
1209  f_is_right_splitable = pair.get_right()->isSplitable();
1210  assert(f_is_left_splitable || f_is_right_splitable);
1211
1212  if (f_is_left_splitable)
1213  {
1214  tree.split_node(pair.get_left());
1215  }
1216
1217  if (f_is_right_splitable)
1218  {
1219  tree.split_node(pair.get_right());
1220  }
1221
1222  if (f_is_left_splitable
1223  && f_is_right_splitable)
1224  {
1226  pair.get_left()->get_left(),
1227  pair.get_right()->get_left());
1229  pair.get_left()->get_left(),
1230  pair.get_right()->get_right());
1231
1232  // to avoid exponential blowup
1233  if (pair.get_left() != pair.get_right())
1235  pair.get_left()->get_right(),
1236  pair.get_right()->get_left());
1237
1239  pair.get_left()->get_right(),
1240  pair.get_right()->get_right());
1241  return;
1242  }
1243
1244  if (f_is_left_splitable)
1245  {
1247  pair.get_left()->get_left(),
1248  pair.get_right());
1250  pair.get_left()->get_right(),
1251  pair.get_right());
1252  return;
1253  }
1254
1255  if (f_is_right_splitable)
1256  {
1258  pair.get_left(),
1259  pair.get_right()->get_left());
1261  pair.get_left(),
1262  pair.get_right()->get_right());
1263  return;
1264  }
1265 }
1266
1267
1269  gdiam_real eps)
1270 {
1271  //gdiam_real diam;
1272  GPointPair pair;
1273
1274  GTreeDiamAlg* pAlg = new GTreeDiamAlg(start, size, UDM_SIMPLE);
1275  pAlg->compute_by_heap(eps);
1276
1277  pair = pAlg->getDiameter();
1278
1279  delete pAlg;
1280
1281  return pair;
1282 }
1283
1284
1286  gdiam_real eps)
1287 {
1288  //gdiam_real diam;
1289  GPointPair pair;
1290
1291  GTreeDiamAlg* pAlg = new GTreeDiamAlg(start, size, UDM_BIG_SAMPLE);
1292  pAlg->compute_by_heap(eps);
1293
1294  pair = pAlg->getDiameter();
1295
1296  delete pAlg;
1297
1298  return pair;
1299 }
1300
1301
1303 {
1304  gdiam_point* p_arr;
1305
1306  assert(start != NULL);
1307  assert(size > 0);
1308
1309  p_arr = (gdiam_point*)malloc(sizeof(gdiam_point) * size);
1310  assert(p_arr != NULL);
1311
1312  for (int ind = 0; ind < size; ind++)
1313  {
1314  p_arr[ ind ] = (gdiam_point) & (start[ ind * 3 ]);
1315  }
1316
1317  return p_arr;
1318 }
1319
1320
1322  gdiam_real eps)
1323 {
1324  gdiam_point* p_arr;
1325  GPointPair pair;
1326
1327  p_arr = gdiam_convert(start, size);
1328  pair = gdiam_approx_diam(p_arr, size, eps);
1329  free(p_arr);
1330
1331  return pair;
1332 }
1333
1334
1336  gdiam_real eps)
1337 {
1338  gdiam_point* p_arr;
1339  GPointPair pair;
1340
1341  p_arr = gdiam_convert(start, size);
1342  pair = gdiam_approx_diam_UDM(p_arr, size, eps);
1343  free(p_arr);
1344
1345  return pair;
1346 }
1347
1348
1350  gdiam_real eps,
1351  GBBox* p_ap_bbox)
1352 {
1353  //gdiam_real diam;
1354  GPointPair pair, pair_2;
1355
1356  GTreeDiamAlg* pAlg = new GTreeDiamAlg(start, size, UDM_SIMPLE);
1357  pAlg->compute_by_heap(eps);
1358
1359  pair = pAlg->getDiameter();
1360
1361  /* Comput ethe resulting diameter */
1362  gdiam_point_t dir, dir_2, dir_3;
1363
1364  dir[ 0 ] = pair.p[ 0 ] - pair.q[ 0 ];
1365  dir[ 1 ] = pair.p[ 1 ] - pair.q[ 1 ];
1366  dir[ 2 ] = pair.p[ 2 ] - pair.q[ 2 ];
1367  pnt_normalize(dir);
1368
1369  // printf( "Before computing second direction!\n" );
1370  //fflush( stdout );
1371  pAlg->compute_by_heap_proj(eps, dir);
1372  //printf( "second direction computed!\n" );
1373  //fflush( stdout );
1374
1375  pair_2 = pAlg->getDiameter();
1376  dir_2[ 0 ] = pair_2.p[ 0 ] - pair_2.q[ 0 ];
1377  dir_2[ 1 ] = pair_2.p[ 1 ] - pair_2.q[ 1 ];
1378  dir_2[ 2 ] = pair_2.p[ 2 ] - pair_2.q[ 2 ];
1379
1380  gdiam_real prd;
1381
1382  prd = pnt_dot_prod(dir, dir_2);
1383  dir_2[ 0 ] = dir_2[ 0 ] - prd * dir[ 0 ];
1384  dir_2[ 1 ] = dir_2[ 1 ] - prd * dir[ 1 ];
1385  dir_2[ 2 ] = dir_2[ 2 ] - prd * dir[ 2 ];
1386
1387  pnt_normalize(dir);
1388  pnt_normalize(dir_2);
1389
1390  pnt_cross_prod(dir, dir_2, dir_3);
1391  pnt_normalize(dir_3);
1392
1393  gdiam_bbox bbox;
1394  GBBox ap_bbox;
1395
1396  if ((pnt_length(dir_2) < 1e-6)
1397  && (pnt_length(dir_3) < 1e-6))
1398  {
1399  gdiam_generate_orthonormal_base(dir, dir_2, dir_3);
1400  }
1401
1402  if ((pnt_length(dir) == 0.0)
1403  || (pnt_length(dir_2) < 1e-6)
1404  || (pnt_length(dir_3) < 1e-6))
1405  {
1406  gdiam_generate_orthonormal_base(dir, dir_2, dir_3);
1407  pnt_dump(dir);
1408  pnt_dump(dir_2);
1409  pnt_dump(dir_3);
1410
1411  fflush(stdout);
1412  fflush(stderr);
1413  assert(false);
1414  }
1415
1416  bbox.init(dir, dir_2, dir_3);
1417  ap_bbox.init();
1418
1419  for (int ind = 0; ind < size; ind++)
1420  {
1421  bbox.bound(start[ ind ]);
1422  ap_bbox.bound(start[ ind ]);
1423  }
1424
1425  //printf( "Axis parallel bounding box:\n" );
1426  //ap_bbox.dump();
1427
1428  //printf( "Arbitrarily oriented bounding box:\n" );
1429  //bbox.dump();
1430
1431  delete pAlg;
1432
1433  if (ap_bbox.volume() < bbox.volume())
1434  {
1435  bbox.init(ap_bbox);
1436  }
1437
1438  if (p_ap_bbox != NULL)
1439  {
1440  *p_ap_bbox = ap_bbox;
1441  }
1442
1443  return bbox;
1444 }
1445
1446
1448  gdiam_real eps)
1449 {
1450  gdiam_point* p_arr;
1451  gdiam_bbox gbbox;
1452
1453  p_arr = gdiam_convert(start, size);
1454  gbbox = gdiam_approx_const_mvbb(p_arr, size, eps, NULL);
1455  free(p_arr);
1456
1457  return gbbox;
1458 }
1459
1460 /*-----------------------------------------------------------------------
1461  * Code for computing best bounding box in a given drection
1462 \*-----------------------------------------------------------------------*/
1463
1465  gdiam_point out1,
1466  gdiam_point out2)
1467 {
1468  assert(pnt_length(in) > 0.0);
1469
1470  pnt_normalize(in);
1471
1472  // stupid cases..
1473  if (in[ 0 ] == 0.0)
1474  {
1475  if (in[ 1 ] == 0.0)
1476  {
1477  pnt_init_normalize(out1, 1, 0, 0);
1478  pnt_init_normalize(out2, 0, 1, 0);
1479  return;
1480  }
1481
1482  if (in[ 2 ] == 0.0)
1483  {
1484  pnt_init_normalize(out1, 1, 0, 0);
1485  pnt_init_normalize(out2, 0, 0, 1);
1486  return;
1487  }
1488
1489  pnt_init_normalize(out1, 0, -in[ 2 ], in[ 1 ]);
1490  pnt_init_normalize(out2, 1, 0, 0);
1491  return;
1492  }
1493
1494  if ((in[ 1 ] == 0.0) && (in[ 2 ] == 0.0))
1495  {
1496  pnt_init_normalize(out1, 0, 1, 0);
1497  pnt_init_normalize(out2, 0, 0, 1);
1498  return;
1499  }
1500
1501  if (in[ 1 ] == 0.0)
1502  {
1503  pnt_init_normalize(out1, -in[ 2 ], 0, in[ 0 ]);
1504  pnt_init_normalize(out2, 0, 1, 0);
1505  return;
1506  }
1507
1508  if (in[ 2 ] == 0.0)
1509  {
1510  pnt_init_normalize(out1, -in[ 1 ], in[ 0 ], 0);
1511  pnt_init_normalize(out2, 0, 0, 1);
1512  return;
1513  }
1514
1515  // all entries in the vector are not zero.
1516  pnt_init_normalize(out1, -in[ 1 ], in[ 0 ], 0);
1517  pnt_cross_prod(in, out1, out2);
1518
1519  pnt_normalize(out2);
1520 }
1521
1522
1523 class point2d
1524 {
1525 public:
1528
1529  void init(gdiam_point pnt,
1530  gdiam_point base_x,
1531  gdiam_point base_y)
1532  {
1533  src = pnt;
1534  x = pnt_dot_prod(pnt, base_x);
1535  y = pnt_dot_prod(pnt, base_y);
1536  }
1537
1539  {
1540  return sqrt((x - pnt.x) * (x - pnt.x)
1541  + (y - pnt.y) * (y - pnt.y));
1542  }
1543
1544  void dump() const
1545  {
1546  printf("(%g, %g)\n", x, y);
1547  }
1548  bool equal(const point2d& pnt) const
1549  {
1550  return ((x == pnt.x) && (y == pnt.y));
1551  }
1552  bool equal_real(const point2d& pnt) const
1553  {
1554  return ((fabs(x - pnt.x) < 1e-8)
1555  && (fabs(y - pnt.y) < 1e-8));
1556  }
1557 };
1558
1560 class vec_point_2d : public std::vector<point2d_ptr> {};
1561
1562
1563
1564 inline ldouble Area(const point2d& a,
1565  const point2d& b,
1566  const point2d& c)
1567 {
1568  double x1, y1, x2, y2, len;
1569
1570  x1 = b.x - a.x;
1571  y1 = b.y - a.y;
1572
1573  x2 = c.x - a.x;
1574  y2 = c.y - a.y;
1575
1576  if ((x1 != 0.0) || (y1 != 0.0))
1577  {
1578  len = sqrt(x1 * x1 + y1 * y1);
1579  x1 /= len;
1580  y1 /= len;
1581  }
1582
1583  if ((x2 != 0.0) || (y2 != 0.0))
1584  {
1585  len = sqrt(x2 * x2 + y2 * y2);
1586  x2 /= len;
1587  y2 /= len;
1588  }
1589
1590  ldouble area;
1591
1592  area = x1 * y2 - x2 * y1;
1593  //area = (ldouble)a.x * (ldouble)b.y - (ldouble)a.y * (ldouble)b.x +
1594  // (ldouble)a.y * (ldouble)c.x - (ldouble)a.x * (ldouble)c.y +
1595  // (ldouble)b.x * (ldouble)c.y - (ldouble)c.x * (ldouble)b.y;
1596
1597  //printf( "area: %g\n", area );
1598  return area;
1599 }
1600
1601
1602 inline int AreaSign(const point2d& a,
1603  const point2d& b,
1604  const point2d& c)
1605 {
1606  ldouble area;
1607
1608  area = a.x * b.y - a.y * b.x +
1609  a.y * c.x - a.x * c.y +
1610  b.x * c.y - c.x * b.y;
1611
1612  //printf( "area: %g\n", area );
1613  return ((area < (ldouble)0.0) ? -1 :
1614  ((area > (ldouble)0.0) ? 1 : 0));
1615 }
1616
1617 inline bool Left(const point2d& a,
1618  const point2d& b,
1619  const point2d& c)
1620 {
1621  return AreaSign(a, b, c) > 0;
1622 }
1623
1624 inline bool Left_colinear(const point2d& a,
1625  const point2d& b,
1626  const point2d& c)
1627 {
1628  return AreaSign(a, b, c) >= 0;
1629 }
1630
1631
1633 {
1634 public:
1636
1637  bool operator()(const point2d_ptr& a, const point2d_ptr& b)
1638  {
1639  int sgn;
1640  gdiam_real len1, len2;
1641
1642  if (a->equal(*b))
1643  {
1644  return false;
1645  }
1646
1647  assert(a != NULL);
1648  assert(b != NULL);
1649
1650  if (a->equal(base))
1651  {
1652  //assert( false );
1653  return true;
1654  }
1655
1656  if (b->equal(base))
1657  {
1658  // assert( false );
1659  return false;
1660  }
1661
1662  //printf( "before: %p %p...\n", a, b );
1663  //fflush( stdout );
1664  /*
1665  if ( a->equal( base ) )
1666  return true;
1667  if ( b->equal( base ) )
1668  return false;
1669  */
1670  //printf( "a: (%g, %g)\n", a->x, a->y );
1671  //printf( "b: (%g, %g)\n", b->x, b->y );
1672  //fflush( stdout );
1673
1674  sgn = AreaSign(base, *a, *b);
1675
1676  if (sgn != 0)
1677  {
1678  return (sgn > 0);
1679  }
1680
1681  len1 = base.dist(*a);
1682  len2 = base.dist(*b);
1683
1684  //printf( "bogi!\n" );
1685  //fflush( stdout );
1686
1687  return (len1 > len2);
1688  }
1689 };
1690
1691
1693  int& index)
1694 {
1695  point2d_ptr min_pnt = in[ 0 ];
1696
1697  index = 0;
1698
1699  for (int ind = 0; ind < (int)in.size(); ind++)
1700  {
1701  if (in[ ind ]->y < min_pnt->y)
1702  {
1703  min_pnt = in[ ind ];
1704  index = ind;
1705  }
1706  else if ((in[ ind ]->y == min_pnt->y)
1707  && (in[ ind ]->x < min_pnt->x))
1708  {
1709  min_pnt = in[ ind ];
1710  index = ind;
1711  }
1712  }
1713
1714  return min_pnt;
1715 }
1716
1717 void dump(vec_point_2d& vec)
1718 {
1719  for (int ind = 0; ind < (int)vec.size(); ind++)
1720  {
1721  printf("-- %11d (%-11g, %-11g)\n",
1722  ind, vec[ ind ]->x,
1723  vec[ ind ]->y);
1724  }
1725
1726  printf("\n\n");
1727 }
1728
1729
1730 static void print_pnt(point2d_ptr pnt)
1731 {
1732  printf("(%g, %g)\n", pnt->x, pnt->y);
1733 }
1734
1735
1737 {
1738  ldouble area;
1739
1740  //dump( ch );
1741  //fflush( stdout );
1742  //fflush( stderr );
1743  for (int ind = 0; ind < (int)(ch.size() - 2); ind++)
1744  assert(Left(*(ch[ ind ]), *(ch[ ind + 1 ]),
1745  *(ch[ ind + 2 ])));
1746
1747  assert(Left(*(ch[ ch.size() - 2 ]), *(ch[ ch.size() - 1 ]),
1748  *(ch[ 0 ])));
1749  assert(Left(*(ch[ ch.size() - 1 ]), *(ch[ 0 ]),
1750  *(ch[ 1 ])));
1751
1752  for (int ind = 0; ind < (int)(in.size() - 2); ind++)
1753  {
1754  for (int jnd = 0; jnd < (int)(ch.size() - 1); jnd++)
1755  {
1756  if (ch[ jnd ]->equal_real(*(in[ ind ])))
1757  {
1758  continue;
1759  }
1760
1761  if (ch[ jnd + 1 ]->equal(*(in[ ind ])))
1762  {
1763  continue;
1764  }
1765
1766  if (! Left_colinear(*(ch[ jnd ]), *(ch[ jnd + 1 ]),
1767  *(in[ ind ])))
1768  {
1769  area = Area(*(ch[ jnd ]), *(ch[ jnd + 1 ]),
1770  *(in[ ind ]));
1771
1772  if (fabs(area) < 1e-12)
1773  {
1774  continue;
1775  }
1776
1777  printf("Failure in progress!\n\n");
1778  print_pnt(ch[ jnd ]);
1779  print_pnt(ch[ jnd + 1 ]);
1780  print_pnt(in[ ind ]);
1781
1782  //dump( ch );
1783  printf("ch[ jnd ]: (%g, %g)\n",
1784  ch[ jnd ]->x, ch[ jnd ]->y);
1785  printf("ch[ jnd + 1 ]: (%g, %g)\n",
1786  ch[ jnd + 1 ]->x, ch[ jnd + 1 ]->y);
1787  printf("ch[ ind ]: (%g, %g)\n",
1788  in[ ind ]->x, in[ ind ]->y);
1789  printf("Area: %g\n", (double)Area(*(ch[ jnd ]),
1790  *(ch[ jnd + 1 ]),
1791  *(in[ ind ])));
1792  printf("jnd: %d, jnd + 1: %d, ind: %d\n",
1793  jnd, jnd + 1, ind);
1794  fflush(stdout);
1795  fflush(stderr);
1796
1797  assert(Left(*(ch[ jnd ]), *(ch[ jnd + 1 ]),
1798  *(in[ ind ])));
1799  }
1800  }
1801
1802  if (ch[ 0 ]->equal(*(in[ ind ])))
1803  {
1804  continue;
1805  }
1806
1807  if (ch[ ch.size() - 1 ]->equal(*(in[ ind ])))
1808  {
1809  continue;
1810  }
1811
1812  assert(Left(*(ch[ ch.size() - 1 ]),
1813  *(ch[ 0 ]),
1814  *(in[ ind ])));
1815  }
1816
1817  printf("Convex hull seems to be ok!\n");
1818 }
1819
1820
1821
1822 static void remove_duplicate(vec_point_2d& in,
1823  point2d_ptr ptr,
1824  int start)
1825 {
1826  int dest;
1827
1828  dest = start;
1829
1830  while (start < (int)in.size())
1831  {
1832  if (in[ start ]->equal_real(*ptr))
1833  {
1834  start++;
1835  continue;
1836  }
1837
1838  in[ dest++ ] = in[ start++ ];
1839  }
1840
1841  while ((int)in.size() > dest)
1842  {
1843  in.pop_back();
1844  }
1845 }
1846
1847
1848 static void remove_consecutive_dup(vec_point_2d& in)
1849 {
1850  int dest, pos;
1851
1852  if (in.size() < 1)
1853  {
1854  return;
1855  }
1856
1857  dest = pos = 0;
1858
1859  // copy first entry
1860  in[ dest++ ] = in[ pos++ ];
1861
1862  while (pos < ((int)in.size() - 1))
1863  {
1864  if (in[ pos - 1 ]->equal_real(*(in[ pos ])))
1865  {
1866  pos++;
1867  continue;
1868  }
1869
1870  in[ dest++ ] = in[ pos++ ];
1871  }
1872
1873  in[ dest++ ] = in[ pos++ ];
1874
1875  while ((int)in.size() > dest)
1876  {
1877  in.pop_back();
1878  }
1879 }
1880
1881
1883 {
1884
1885  CompareByAngle comp;
1886  int ind, position;
1887
1888  assert(in.size() > 1);
1889
1890  comp.base = *(get_min_point(in, position));
1891  std::swap(in[ 0 ], in[ position ]);
1892
1893  remove_duplicate(in, in[ 0 ], 1);
1894
1895  /*
1896  printf( "comp.base: (%g, %g)\n",
1897  comp.base.x,
1898  comp.base.y );
1899  */
1900  for (int ind = 0; ind < (int)in.size(); ind++)
1901  {
1902  assert(in[ ind ] != NULL);
1903  }
1904
1905  /*
1906  if ( in.size() == 24 ) {
1907  dump( in );
1908  fflush( stdout );
1909  fflush( stderr );
1910  }*/
1911
1912  //printf( "sort( %d, %d, comp )\n", 1, in.size() );
1913  sort(in.begin() + 1, in.end(), comp);
1914  remove_consecutive_dup(in);
1915
1916  //dump( in );
1917  /*
1918  for ( int ind = 0; ind < (int)in.size(); ind++ ) {
1919  double x_delta, y_delta;
1920
1921  x_delta = in[ ind ]->x - comp.base.x;
1922  y_delta = in[ ind ]->y - comp.base.y;
1923
1924  printf( "-- %11d (%-11g, %-11g) atan2: %-11g\n",
1925  ind, in[ ind ]->x,
1926  in[ ind ]->y,
1927  atan2( y_delta, x_delta ) );
1928  }*/
1929
1930  //fflush( stdout );
1931  //printf( "pushing...\n" );
1932  //fflush( stdout );
1933  out.push_back(in[ in.size() - 1 ]);
1934  out.push_back(in[ 0 ]);
1935
1936  // perform the graham scan
1937  ind = 1;
1938  int last;
1939
1940  while (ind < (int)in.size())
1941  {
1942  if (out[ out.size() - 1 ]->equal_real(*(in[ ind ])))
1943  {
1944  ind++;
1945  continue;
1946  }
1947
1948  last = out.size();
1949  assert(last > 1);
1950
1951  if (Left(*(out[ last - 2 ]),
1952  *(out[ last - 1 ]),
1953  *(in[ ind ])))
1954  {
1955  if (! out[ last - 1 ]->equal(*(in[ ind ])))
1956  {
1957  out.push_back(in[ ind ]);
1958  }
1959
1960  ind++;
1961  }
1962  else
1963  {
1964  if (out.size() < 3)
1965  {
1966  dump(out);
1967  printf("Input:\n");
1968  dump(in);
1969  printf("ind: %d\n", ind);
1970  fflush(stdout);
1971  assert(out.size() > 2);
1972  }
1973
1974  out.pop_back();
1975  }
1976  }
1977
1978  // we pushed in[ in.size() -1 ] twice in the output
1979  out.pop_back();
1980
1981  //verify_convex_hull( in, out );
1982 }
1983
1984
1986 {
1989
1990  void get_dir(int ind, gdiam_point_2d out)
1991  {
1992  out[ 0 ] = vertices[ ind ][ 0 ]
1993  - vertices[(ind + 1) % 4 ][ 0 ];
1994  out[ 1 ] = vertices[ ind ][ 1 ]
1995  - vertices[(ind + 1) % 4 ][ 1 ];
1996  }
1997
1998  void dump()
1999  {
2000  printf("--- bbox 2d ----------------------------\n");
2001
2002  for (int ind = 0; ind < 4; ind++)
2003  printf("%d: (%g, %g)\n", ind, vertices[ ind ][0],
2004  vertices[ ind ][1]);
2005
2006  for (int ind = 0; ind < 4; ind++)
2007  {
2008  gdiam_point_2d_t dir;
2009
2010  get_dir(ind, dir);
2011
2012  printf("dir %d: (%g, %g)\n", ind, dir[0],
2013  dir[1]);
2014  }
2015
2016  printf("\\----------------------------------\n");
2017  }
2018
2019 };
2020
2021 #define EPS 1e-6
2023 {
2024  return (((b - EPS) < a)
2025  && (a < (b + EPS)));
2026 }
2027
2029 {
2030 private:
2031  vec_point_2d ch;
2032  gdiam_real* angles;
2033
2034 public:
2035  void dump_ch()
2036  {
2037  for (int ind = 0; ind < (int)ch.size(); ind++)
2038  {
2039  printf("%d: (%g, %g)\n", ind, ch[ ind ]->x,
2040  ch[ ind ]->y);
2041  }
2042  }
2043
2044  void compute_edge_dir(int u)
2045  {
2046  int u2;
2047  double ang;
2048  double x1, y1, x2, y2;
2049
2050  u2 = (u + 1) % ch.size();
2051
2052  x1 = ch[ u ]->x;
2053  y1 = ch[ u ]->y;
2054  x2 = ch[ u2 ]->x;
2055  y2 = ch[ u2 ]->y;
2056
2057  ang = atan2(y2 - y1, x2 - x1);
2058
2059  if (ang < 0)
2060  {
2061  ang += 2.0 * PI;
2062  }
2063
2064  angles[ u ] = ang;
2065
2066  // floating point nonsence. A huck to solve this...
2067  if ((u > 0)
2068  && (angles[ u ] < angles[ u - 1 ])
2069  && eq_real(angles[ u ], angles[ u - 1 ]))
2070  {
2071  angles[ u ] = angles[ u - 1 ];
2072  }
2073
2074  /*
2075  printf( "angles[ %4d ]: %.20f ", u, ang );
2076  if ( u > 0 )
2077  printf( "%2d", (angles[ u ] >= angles[ u - 1 ] ) );
2078  if ( u > 1 ) {
2079  printf( " area: %.20f ", (double)Area( *( ch[ u - 2 ] ),
2080  *( ch[ u - 1 ] ),
2081  *( ch[ u ] ) ) );
2082  }
2083  printf( "\n" );
2084  */
2085  }
2086
2087  /* Finding the first vertex whose edge is over a given angle */
2088  int find_vertex(int start_vertex,
2089  double trg_angle)
2090  {
2091  double prev_angle, curr_angle;
2092  bool f_found;
2093  int ver, prev_vertex;
2094
2095  f_found = false;
2096
2097  prev_vertex = start_vertex - 1;
2098
2099  if (prev_vertex < 0)
2100  {
2101  prev_vertex = ch.size() - 1;
2102  }
2103
2104  prev_angle = angles[ prev_vertex ];
2105  ver = start_vertex;
2106
2107  while (! f_found)
2108  {
2109  curr_angle = angles[ ver ];
2110
2111  if (prev_angle <= curr_angle)
2112  f_found = ((prev_angle < trg_angle)
2113  && (trg_angle <= curr_angle));
2114  else
2115  f_found = ((prev_angle < trg_angle)
2116  || (trg_angle <= curr_angle));
2117
2118  if (f_found)
2119  {
2120  break;
2121  }
2122  else
2123  {
2124  prev_angle = curr_angle;
2125  }
2126
2127  ver = (ver + 1) % ch.size();
2128  }
2129
2130  return ver;
2131  }
2132
2133
2134  /* Computing the intersection point of two lines, each given by a
2135  point and slope */
2136  void compute_crossing(int a_ind, double a_angle,
2137  int b_ind, double b_angle,
2138  gdiam_real& x,
2139  gdiam_real& y)
2140  {
2141  double a_slope, b_slope, x_a, x_b, y_a, y_b;
2142
2143  if (eq_real(a_angle, 0.0) || eq_real(a_angle, PI))
2144  {
2145  x = ch[ b_ind ]->x;
2146  y = ch[ a_ind ]->y;
2147  return;
2148  }
2149
2150  if (eq_real(a_angle, PI / 2.0) || eq_real(a_angle, PI * 1.5))
2151  {
2152  x = ch[ a_ind ]->x;
2153  y = ch[ b_ind ]->y;
2154  return;
2155  }
2156
2157  a_slope = tan(a_angle);
2158  b_slope = tan(b_angle);
2159  x_a = ch[ a_ind ]->x;
2160  y_a = ch[ a_ind ]->y;
2161  x_b = ch[ b_ind ]->x;
2162  y_b = ch[ b_ind ]->y;
2163
2164  double a_const, b_const, x_out, y_out;
2165
2166  // Let see:
2167  // l_a = y = a_slope * x + a_const ->
2168  // y_a = a_slope * x_a + a_const (since (x_a, y_a) in l_a
2169  // a_const = y_a - a_slope * x_a
2170  a_const = y_a - a_slope * x_a;
2171  //printf( "a_const: %g, y_a: %g, a_slope: %g, x_a: %g\n",
2172  // a_const, y_a, a_slope, x_a );
2173
2174  b_const = y_b - b_slope * x_b;
2175  //printf( "b_const: %g, y_b: %g, b_slope: %g, x_b: %g\n",
2176  // b_const, y_b, b_slope, x_b );
2177
2178  x_out = - (a_const - b_const) / (a_slope - b_slope);
2179  y_out = a_slope * x_out + a_const;
2180  /*
2181  x = (y_b - y_a + x_a * a_slope - x_b * a_slope)
2182  / (a_slope - b_slope);
2183  y = ((y_a - x_a * a_slope) / a_slope - (y_b - x_b * b_slope)
2184  / b_slope) /
2185  (1.0 / a_slope - 1.0 / b_slope );
2186  printf( "gill: (%g, %g)\nMine: (%g, %g)\n",
2187  x, y, x_out, y_out );
2188  */
2189  //printf( "a_angle: %g, b_angle: %g\n",
2190  // a_angle / PI, b_angle/ PI );
2191  //printf( "a_slope: %g, b_slope: %g\n",
2192  // a_slope, b_slope );
2193  //printf( "on line_b: %g\n",
2194  // b_slope * x_out + b_const - y_out );
2195  x = x_out;
2196  y = y_out;
2197  //printf( "line_a : y = %g * x + %g\n",
2198  // a_slope, a_const );
2199  //printf( "line_b : y = %g * x + %g\n",
2200  // b_slope, b_const );
2201  //printf( "a: (%g, %g)\n"
2202  // "b: (%g, %g)\n",
2203  // x_a, y_a, x_b, y_b );
2204
2205  //printf( "out: (%g, %g)\n", x, y );
2206  }
2207
2208
2209  void get_bbox(int a_ind, gdiam_real a_angle,
2210  int b_ind, gdiam_real b_angle,
2211  int c_ind, gdiam_real c_angle,
2212  int d_ind, gdiam_real d_angle,
2213  bbox_2d_info& bbox,
2214  gdiam_real& area)
2215  {
2216  /*
2217  printf( "get_bbox: (%d, %g,\n"
2218  "%d, %g\n"
2219  "%d, %g\n"
2220  "%d, %g )\n",
2221  a_ind, a_angle / PI,
2222  b_ind, b_angle / PI,
2223  c_ind, c_angle / PI,
2224  d_ind, d_angle / PI );*/
2225
2226  compute_crossing(a_ind, a_angle, b_ind, b_angle,
2227  bbox.vertices[ 0 ][ 0 ],
2228  bbox.vertices[ 0 ][ 1 ]);
2229  compute_crossing(b_ind, b_angle, c_ind, c_angle,
2230  bbox.vertices[ 1 ][ 0 ],
2231  bbox.vertices[ 1 ][ 1 ]);
2232  compute_crossing(c_ind, c_angle, d_ind, d_angle,
2233  bbox.vertices[ 2 ][ 0 ],
2234  bbox.vertices[ 2 ][ 1 ]);
2235  compute_crossing(d_ind, d_angle, a_ind, a_angle,
2236  bbox.vertices[ 3 ][ 0 ],
2237  bbox.vertices[ 3 ][ 1 ]);
2238
2239  area = pnt_distance_2d(bbox.vertices[ 0 ],
2240  bbox.vertices[ 1 ])
2241  * pnt_distance_2d(bbox.vertices[ 0 ],
2242  bbox.vertices[ 3 ]);
2243  }
2244
2245  /* Given one angle (bbx direction), find the other three */
2246  void get_angles(int ind,
2247  gdiam_real& angle_1,
2248  gdiam_real& angle_2,
2249  gdiam_real& angle_3)
2250  {
2251  double angle;
2252
2253  angle = angles[ ind ];
2254  angle_1 = angle + PI * 0.5;
2255  //printf( "angle = %g\n", angle / PI );
2256  //printf( "angle1 = %g\n", angle_1 / PI );
2257
2258  if (angle_1 >= (2.0 * PI))
2259  {
2260  angle_1 -= 2.0 * PI;
2261  }
2262
2263  //printf( "angle1 = %g\n", angle_1 / PI );
2264
2265
2266  angle_2 = angle + PI;
2267
2268  if (angle_2 >= (2.0 * PI))
2269  {
2270  angle_2 -= 2.0 * PI;
2271  }
2272
2273  angle_3 = angle + 1.5 * PI;
2274
2275  if (angle_3 >= (2.0 * PI))
2276  {
2277  angle_3 -= 2.0 * PI;
2278  }
2279  }
2280
2282  gdiam_real& min_area)
2283  {
2284  int u, v, s, t;
2285  gdiam_real ang1, ang2, ang3, tmp_area;
2286  bbox_2d_info tmp_bbox;
2287
2288  angles = (gdiam_real*)malloc(sizeof(gdiam_real)
2289  * (int)ch.size());
2290  assert(angles != NULL);
2291
2292  /* Pre-computing all edge directions */
2293  for (u = 0; u < (int)ch.size(); u ++)
2294  {
2295  compute_edge_dir(u);
2296  }
2297
2298  /* Initializing the search */
2299  //printf( "_angles[ 0 ]: %g\n", angles[ 0 ] );
2300  get_angles(0, ang1, ang2, ang3);
2301
2302  //printf( "_angles[ 0 ]: %g\n", angles[ 0 ] );
2303  s = find_vertex(0, ang1);
2304  //printf( "_angles[ 0 ]: %g\n", angles[ 0 ] );
2305  v = find_vertex(s, ang2);
2306  t = find_vertex(v, ang3);
2307
2308  //printf( "_angles[ 0 ]: %g\n", angles[ 0 ] );
2309  get_bbox(0, angles[ 0 ],
2310  s, ang1,
2311  v, ang2,
2312  t, ang3,
2313  min_bbox, min_area);
2314  //min_bbox.dump();
2315
2316  //printf( "min_area: %g\n", min_area );
2317  /* Performing a double rotating calipers */
2318  for (u = 1; u < (int)ch.size(); u ++)
2319  {
2320  get_angles(u, ang1, ang2, ang3);
2321  s = find_vertex(s, ang1);
2322  v = find_vertex(v, ang2);
2323  t = find_vertex(t, ang3);
2324  get_bbox(u, angles[ u ],
2325  s, ang1,
2326  v, ang2,
2327  t, ang3,
2328  tmp_bbox, tmp_area);
2329
2330  //tmp_bbox.dump();
2331  //printf( "tmp_area: %g\n", tmp_area );
2332  if (min_area > tmp_area)
2333  {
2334  min_area = tmp_area;
2335  min_bbox = tmp_bbox;
2336  }
2337  }
2338
2339  free(angles);
2340  angles = NULL;
2341  }
2342
2344  bbox_2d_info& min_bbox,
2345  gdiam_real& min_area)
2346  {
2347  //printf( "before cunvex hull\n" );
2348  convex_hull(points, ch);
2349  //printf( "cunvex hull done!\n" );
2350  //dump_ch();
2351  compute_min_bbox_inner(min_bbox, min_area);
2352  //min_bbox.dump();
2353  }
2354
2355 };
2356
2358  int size)
2359 {
2360  for (int ind = 0; ind < size; ind++)
2361  {
2362  printf("%d: (%g, %g, %g)\n", ind,
2363  in_arr[ ind ][ 0 ],
2364  in_arr[ ind ][ 1 ],
2365  in_arr[ ind ][ 2 ]);
2366  }
2367 }
2368
2369 #define ZERO_EPS (1e-6)
2370
2371 static void construct_base_inner(gdiam_point pnt1,
2372  gdiam_point pnt2,
2373  gdiam_point pnt3)
2374 {
2375  pnt_normalize(pnt1);
2376  pnt_normalize(pnt2);
2377  pnt_normalize(pnt3);
2378
2379  if ((pnt_length(pnt1) < ZERO_EPS)
2380  && (pnt_length(pnt2) < ZERO_EPS)
2381  && (pnt_length(pnt3) < ZERO_EPS))
2382  {
2383  assert(false);
2384  }
2385
2386  if ((pnt_length(pnt1) < ZERO_EPS)
2387  && (pnt_length(pnt2) < ZERO_EPS))
2388  {
2389  gdiam_generate_orthonormal_base(pnt3, pnt1, pnt2);
2390  return;
2391  }
2392
2393  if ((pnt_length(pnt1) < ZERO_EPS)
2394  && (pnt_length(pnt3) < ZERO_EPS))
2395  {
2396  gdiam_generate_orthonormal_base(pnt2, pnt1, pnt3);
2397  return;
2398  }
2399
2400  if ((pnt_length(pnt2) < ZERO_EPS)
2401  && (pnt_length(pnt3) < ZERO_EPS))
2402  {
2403  gdiam_generate_orthonormal_base(pnt1, pnt2, pnt3);
2404  return;
2405  }
2406
2407  if (pnt_length(pnt1) < ZERO_EPS)
2408  {
2409  pnt_cross_prod(pnt2, pnt3, pnt1);
2410  return;
2411  }
2412
2413  if (pnt_length(pnt2) < ZERO_EPS)
2414  {
2415  pnt_cross_prod(pnt1, pnt3, pnt2);
2416  return;
2417  }
2418
2419  if (pnt_length(pnt3) < ZERO_EPS)
2420  {
2421  pnt_cross_prod(pnt1, pnt2, pnt3);
2422  return;
2423  }
2424 }
2425
2427  gdiam_point pnt2,
2428  gdiam_point pnt3)
2429 {
2430  construct_base_inner(pnt1, pnt2, pnt3);
2432  pnt1);
2433  pnt_normalize(pnt2);
2434
2436  pnt1);
2438  pnt2);
2439  pnt_normalize(pnt3);
2440 }
2441
2442
2444 {
2445 private:
2446  gdiam_point_t base_x, base_y, base_proj;
2447  point2d* arr;
2448  vec_point_2d points, ch;
2449  gdiam_bbox bbox;
2450  gdiam_point* in_arr;
2451  int size;
2452
2453 public:
2455  {
2456  term();
2457  }
2458  void init(gdiam_point dir,
2459  gdiam_point* _in_arr,
2460  int _size)
2461  {
2462  arr = NULL;
2463
2464  if (pnt_length(dir) == 0.0)
2465  {
2466  dump_points(_in_arr, _size);
2467  pnt_dump(dir);
2468  fflush(stdout);
2469  fflush(stderr);
2470  assert(pnt_length(dir) > 0.0);
2471  }
2472
2473  size = _size;
2474  in_arr = _in_arr;
2475  assert(size > 0);
2476
2477  pnt_copy(base_proj, dir);
2478  gdiam_generate_orthonormal_base(base_proj, base_x, base_y);
2479
2480  arr = (point2d*)malloc(sizeof(point2d) * size);
2481  assert(arr != 0);
2482
2483  for (int ind = 0; ind < size; ind++)
2484  {
2485  arr[ ind ].init(in_arr[ ind ], base_x, base_y);
2486  points.push_back(&(arr[ ind ]));
2487  }
2488  }
2489
2490  void compute()
2491  {
2492  MinAreaRectangle mar;
2493  bbox_2d_info min_bbox;
2494  gdiam_real min_area;
2495  double x1, y1, x2, y2;
2496  gdiam_point_t out_1, out_2;
2497
2498  mar.compute_min_bbox(points, min_bbox, min_area);
2499
2500  // We take the resulting min_area rectangle and lift it back to 3d...
2501  x1 = min_bbox.vertices[ 0 ][ 0 ] -
2502  min_bbox.vertices[ 1 ][ 0 ];
2503  y1 = min_bbox.vertices[ 0 ][ 1 ] -
2504  min_bbox.vertices[ 1 ][ 1 ];
2505  x2 = min_bbox.vertices[ 0 ][ 0 ] -
2506  min_bbox.vertices[ 3 ][ 0 ];
2507  y2 = min_bbox.vertices[ 0 ][ 1 ] -
2508  min_bbox.vertices[ 3 ][ 1 ];
2509
2510  //printf( "prod: %g\n", x1 * x2 + y1 * y2 );
2511
2512  double len;
2513
2514  len = sqrt(x1 * x1 + y1 * y1);
2515
2516  if (len > 0.0)
2517  {
2518  x1 /= len;
2519  y1 /= len;
2520  }
2521
2522  len = sqrt(x2 * x2 + y2 * y2);
2523
2524  if (len > 0.0)
2525  {
2526  x2 /= len;
2527  y2 /= len;
2528  }
2529
2530  pnt_zero(out_1);
2533  pnt_normalize(out_1);
2534
2535  pnt_zero(out_2);
2538  pnt_normalize(out_2);
2539
2540  construct_base(base_proj, out_1, out_2);
2541
2542  if ((!(fabs(pnt_dot_prod(base_proj, out_1)) < 1e-6))
2543  || (!(fabs(pnt_dot_prod(base_proj, out_2)) < 1e-6))
2544  || (!(fabs(pnt_dot_prod(out_1, out_2)) < 1e-6)))
2545  {
2546  printf("should be all close to zero: %g, %g, %g\n",
2547  pnt_dot_prod(base_proj, out_1),
2548  pnt_dot_prod(base_proj, out_2),
2549  pnt_dot_prod(out_1, out_2));
2550  pnt_dump(base_proj);
2551  pnt_dump(out_1);
2552  pnt_dump(out_2);
2553
2554
2555  printf("real points:\n");
2556  dump_points(in_arr, size);
2557
2558  printf("points on CH:\n");
2559  dump(points);
2560
2561  printf("BBox:\n");
2562  min_bbox.dump();
2563
2564  fflush(stdout);
2565  fflush(stderr);
2566  assert(fabs(pnt_dot_prod(base_proj, out_1)) < 1e-6);
2567  assert(fabs(pnt_dot_prod(base_proj, out_2)) < 1e-6);
2568  assert(fabs(pnt_dot_prod(out_1, out_2)) < 1e-6);
2569  }
2570
2571  bbox.init(base_proj, out_1, out_2);
2572
2573  for (int ind = 0; ind < size; ind++)
2574  {
2575  bbox.bound(in_arr[ ind ]);
2576  }
2577  }
2578
2580  {
2581  return bbox;
2582  }
2583
2584  void term()
2585  {
2586  if (arr != NULL)
2587  {
2588  free(arr);
2589  arr = NULL;
2590  }
2591  }
2592  void init()
2593  {
2594  }
2595 };
2596
2597
2599  gdiam_bbox bb_out, int times)
2600 {
2601  gdiam_bbox bb_tmp;
2602
2603  //printf( "gdiam_mvbb_optimize called\n" );
2604
2605  for (int ind = 0; ind < times; ind++)
2606  {
2607  ProjPointSet pps;
2608
2609  if (pnt_length(bb_out.get_dir(ind % 3)) == 0.0)
2610  {
2611  printf("Dumping!\n");
2612  bb_out.dump();
2613  fflush(stdout);
2614  continue;
2615  }
2616
2617  pps.init(bb_out.get_dir(ind % 3), start, size);
2618
2619  pps.compute();
2620
2621  //printf( "Old volume: %g\n", bb_out.volume() );
2622  bb_tmp = pps.get_bbox();
2623  //printf( "New volume: %g\n", bb_tmp.volume() );
2624  //printf( "Old bounding box:\n" );
2625  //bb2.dump();
2626
2627  if (bb_tmp.volume() < bb_out.volume())
2628  {
2629  bb_out = bb_tmp;
2630  }
2631  }
2632
2633  return bb_out;
2634 }
2635
2636
2638  gdiam_real eps)
2639 {
2640  gdiam_bbox bb, bb2;
2641
2642  bb = gdiam_approx_const_mvbb(start, size, 0.0, NULL);
2643  bb2 = gdiam_mvbb_optimize(start, size, bb, 10);
2644
2645  //bb2.dump();
2646
2647  return bb2;
2648 }
2649
2650
2651 static int gcd2(int a, int b)
2652 {
2653  if (a == 0 || a == b)
2654  {
2655  return b;
2656  }
2657
2658  if (b == 0)
2659  {
2660  return a;
2661  }
2662
2663  if (a > b)
2664  {
2665  return gcd2(a % b, b);
2666  }
2667  else
2668  {
2669  return gcd2(a, b % a);
2670  }
2671 }
2672
2673
2674 static int gcd3(int a, int b, int c)
2675 {
2676  a = abs(a);
2677  b = abs(b);
2678  c = abs(c);
2679
2680  if (a == 0)
2681  {
2682  return gcd2(b, c);
2683  }
2684
2685  if (b == 0)
2686  {
2687  return gcd2(a, c);
2688  }
2689
2690  if (c == 0)
2691  {
2692  return gcd2(a, b);
2693  }
2694
2695  return gcd2(a, gcd2(b, c));
2696 }
2697
2698
2699 static void try_direction(gdiam_bbox& bb,
2700  gdiam_bbox& bb_out,
2701  gdiam_point* start,
2702  int size,
2703  int x_coef,
2704  int y_coef,
2705  int z_coef)
2706 {
2707  gdiam_bbox bb_new;
2708
2709  if (gcd3(x_coef, y_coef, z_coef) > 1)
2710  {
2711  return;
2712  }
2713
2714  if ((x_coef == 0) && (y_coef == 0)
2715  && (z_coef == 0))
2716  {
2717  return;
2718  }
2719
2720  //printf( "trying: (%d, %d, %d)\n",
2721  // x_coef, y_coef, z_coef );
2722  //fflush( stdout );
2723  gdiam_point_t new_dir;
2724  bb.combine(new_dir, x_coef, y_coef, z_coef);
2725
2726  ProjPointSet pps;
2727
2728  pps.init(new_dir, start, size);
2729
2730  pps.compute();
2731
2732  bb_new = pps.get_bbox();
2733  bb_new = gdiam_mvbb_optimize(start, size, bb_new, 10);
2734
2735  if (bb_new.volume() < bb_out.volume())
2736  {
2737  bb_out = bb_new;
2738  }
2739 }
2740
2741
2743  int grid_size)
2744 {
2745  gdiam_bbox bb, bb_out;
2746
2747  bb_out = bb = gdiam_approx_const_mvbb(start, size, 0.0, NULL);
2748
2749  if (bb.volume() == 0)
2750  {
2751  dump_points(start, size);
2752  printf("1zero volume???\n");
2753  bb.dump();
2754  }
2755
2756  bb_out = bb = gdiam_mvbb_optimize(start, size, bb_out, 10);
2757
2758  if (bb.volume() == 0)
2759  {
2760  printf("2zero volume???\n");
2761  bb.dump();
2762  }
2763
2764  //trying: (1, -4, 4)
2765  //try_direction( bb, bb_out, start, size,
2766  // 1, -4, 4 );
2767
2768  for (int x_coef = -grid_size; x_coef <= grid_size; x_coef++)
2769  for (int y_coef = -grid_size; y_coef <= grid_size; y_coef++)
2770  for (int z_coef = 0; z_coef <= grid_size; z_coef++)
2771  try_direction(bb, bb_out, start, size,
2772  x_coef, y_coef, z_coef);
2773
2774  bb_out = gdiam_mvbb_optimize(start, size, bb_out, 10);
2775
2776  return bb_out;
2777 }
2778
2779
2780 static void register_point(gdiam_point pnt,
2781  gdiam_point* tops,
2782  gdiam_point* bottoms,
2783  int grid_size,
2784  gdiam_bbox& bb)
2785 {
2786  gdiam_point_t pnt_bb, pnt_bottom, pnt_top;
2787  int x_ind, y_ind, position;
2788
2789  bb.get_normalized_coordinates(pnt, pnt_bb);
2790
2791  // x_ind
2792  x_ind = (int)(pnt_bb[ 0 ] * grid_size);
2793  assert((-1 <= x_ind) && (x_ind <= grid_size));
2794
2795  if (x_ind < 0)
2796  {
2797  x_ind = 0;
2798  }
2799
2800  if (x_ind >= grid_size)
2801  {
2802  x_ind = grid_size - 1;
2803  }
2804
2805  // y_ind
2806  y_ind = (int)(pnt_bb[ 1 ] * grid_size);
2807  assert((-1 <= y_ind) && (y_ind <= grid_size));
2808
2809  if (y_ind < 0)
2810  {
2811  y_ind = 0;
2812  }
2813
2814  if (y_ind >= grid_size)
2815  {
2816  y_ind = grid_size - 1;
2817  }
2818
2819  position = x_ind + y_ind * grid_size;
2820
2821  if (tops[ position ] == NULL)
2822  {
2823  tops[ position ] = bottoms[ position ] = pnt;
2824  return;
2825  }
2826
2827  bb.get_normalized_coordinates(tops[ position ], pnt_top);
2828
2829  if (pnt_top[ 2 ] < pnt_bb[ 2 ])
2830  {
2831  tops[ position ] = pnt;
2832  }
2833
2834  bb.get_normalized_coordinates(bottoms[ position ], pnt_bottom);
2835
2836  if (pnt_bottom[ 2 ] > pnt_bb[ 2 ])
2837  {
2838  bottoms[ position ] = pnt;
2839  }
2840 }
2841
2842
2843
2844 // gdiam_convex_sample
2845
2846 // We are given a point set, and (hopefully) a tight fitting
2847 // bounding box. We compute a smplae of the given size that represents
2848 // the point-set. The only guarenteed is that if we use ample of size m,
2849 // we get an approximation of quality about 1/\sqrt{m}. Note that we pad
2850 // the sample if necessary to get the desired size.
2852  gdiam_bbox& bb,
2853  int sample_size)
2854 {
2855  int grid_size, mem_size, grid_entries;
2856  gdiam_point* bottoms, * tops, * out_arr;
2857  int out_count;
2858
2859  assert(sample_size > 1);
2860
2861  // we want the grid to be on the two minor dimensions.
2862  bb.set_third_dim_longest();
2863
2864  grid_size = (int)sqrt(sample_size / 2);
2865
2866  grid_entries = grid_size * grid_size;
2867  mem_size = (int)(sizeof(gdiam_point) * grid_size * grid_size);
2868
2869  bottoms = (gdiam_point*)malloc(mem_size);
2870  tops = (gdiam_point*)malloc(mem_size);
2871  out_arr = (gdiam_point*)malloc(sizeof(gdiam_point) * sample_size);
2872
2873  assert(bottoms != NULL);
2874  assert(tops != NULL);
2875  assert(out_arr != NULL);
2876
2877  for (int ind = 0; ind < grid_entries; ind++)
2878  {
2879  tops[ ind ] = bottoms[ ind ] = NULL;
2880  }
2881
2882  // No we stream the points registering them with the relevant
2883  // shaft in the grid.
2884  for (int ind = 0; ind < size; ind++)
2885  register_point(start[ ind ], tops, bottoms,
2886  grid_size, bb);
2887
2888  out_count = 0;
2889
2890  for (int ind = 0; ind < grid_entries; ind++)
2891  {
2892  if (tops[ ind ] == NULL)
2893  {
2894  continue;
2895  }
2896
2897  out_arr[ out_count++ ] = tops[ ind ];
2898
2899  if (tops[ ind ] != bottoms[ ind ])
2900  {
2901  out_arr[ out_count++ ] = bottoms[ ind ];
2902  }
2903  }
2904
2905  // We dont have neough entries in our sample - lets randomly pick
2906  // points.
2907  while (out_count < sample_size)
2908  {
2909  out_arr[ out_count++ ] = start[ rand() % size ];
2910  }
2911
2912  free(tops);
2913  free(bottoms);
2914
2915  return out_arr;
2916 }
2917
2918
2920  int grid_size, int sample_size)
2921 {
2922  gdiam_bbox bb, bb_new;
2923  gdiam_point* sample;
2924
2925  if (sample_size >= size)
2926  {
2927  return gdiam_approx_mvbb_grid(start, size, grid_size);
2928  }
2929
2930  //printf( "gdiam_approx_mvbb_grid_sample(...) called\n" );
2931  //fflush( stdout );
2932
2933  bb = gdiam_approx_const_mvbb(start, size, 0.0, NULL);
2934  //printf( "guess from const approximation:\n" );
2935
2936  sample = gdiam_convex_sample(start, size, bb, sample_size);
2937
2938  bb_new = gdiam_approx_mvbb_grid(sample, sample_size, grid_size);
2939
2940  //printf( "new volume on sample: %g\n",
2941  // bb_new.volume() );
2942  //fflush( stdout );
2943
2944  for (int ind = 0; ind < size; ind++)
2945  {
2946  bb_new.bound(start[ ind ]);
2947  }
2948
2949  //printf( "new volume on resized input: %g\n",
2950  // bb_new.volume() );
2951
2952  //bb_new = gdiam_mvbb_optimize( start, size, bb_new, 10 );
2953  //printf(s "optimized new volume: %g\n",
2954  // bb_new.volume() );
2955
2956  return bb_new;
2957 }
2958
2959
2960
2962  int grid_size, int sample_size)
2963 {
2964  gdiam_point* p_arr;
2965  gdiam_bbox bb;
2966
2967  //printf( "gdiam_approx_mvbb_grid_sample( ?, %d, ...)\n", size );
2968  //fflush( stdout );
2969
2970  p_arr = gdiam_convert(start, size);
2971  bb = gdiam_approx_mvbb_grid_sample(p_arr, size, grid_size, sample_size);
2972  free(p_arr);
2973
2974  return bb;
2975 }
2976
2977
2978 /* gdiam.C - End of File ------------------------------------------*/
g_heap_pairs_p::g_heap_pairs_p
g_heap_pairs_p()
Definition: gdiam.cpp:833
GFSPPair::get_left
GFSPTreeNode * get_left()
Definition: gdiam.cpp:498
ptrCompareFunc
int(* ptrCompareFunc)(void *aPtr, void *bPtr)
Definition: gdiam.cpp:640
Area
ldouble Area(const point2d &a, const point2d &b, const point2d &c)
Definition: gdiam.cpp:1564
gdiam_bbox::volume
gdiam_real volume() const
Definition: gdiam.h:535
GfxTL::sqrt
VectorXD< D, T > sqrt(const VectorXD< D, T > &a)
Definition: VectorXD.h:662
ProjPointSet::init
void init(gdiam_point dir, gdiam_point *_in_arr, int _size)
Definition: gdiam.cpp:2458
vec_point_2d
Definition: gdiam.cpp:1560
g_heap_pairs_p::top
GFSPPair top()
Definition: gdiam.cpp:860
MinAreaRectangle::compute_crossing
void compute_crossing(int a_ind, double a_angle, int b_ind, double b_angle, gdiam_real &x, gdiam_real &y)
Definition: gdiam.cpp:2136
ZERO_EPS
#define ZERO_EPS
Definition: gdiam.cpp:2369
point2d::y
gdiam_real y
Definition: gdiam.cpp:1526
gdiam_approx_mvbb
gdiam_bbox gdiam_approx_mvbb(gdiam_point *start, int size, gdiam_real eps)
Definition: gdiam.cpp:2637
gdiam_approx_diam_pair
GPointPair gdiam_approx_diam_pair(gdiam_real *start, int size, gdiam_real eps)
Definition: gdiam.cpp:1321
g_heap_pairs_p
Definition: gdiam.cpp:827
point2d::equal
bool equal(const point2d &pnt) const
Definition: gdiam.cpp:1548
ldouble
long double ldouble
Definition: gdiam.cpp:45
Left_colinear
bool Left_colinear(const point2d &a, const point2d &b, const point2d &c)
Definition: gdiam.cpp:1624
pnt_zero
void pnt_zero(gdiam_point dst)
Definition: gdiam.h:90
gdiam_approx_diam_UDM
GPointPair gdiam_approx_diam_UDM(gdiam_point *start, int size, gdiam_real eps)
Definition: gdiam.cpp:1285
CompareByAngle
Definition: gdiam.cpp:1632
gdiam_mvbb_optimize
gdiam_bbox gdiam_mvbb_optimize(gdiam_point *start, int size, gdiam_bbox bb_out, int times)
Definition: gdiam.cpp:2598
GFSPTreeNode::ptr_pnt_rand
const gdiam_point * ptr_pnt_rand()
Definition: gdiam.cpp:143
bbox_2d_info
Definition: gdiam.cpp:1985
gdiam_approx_mvbb_grid
gdiam_bbox gdiam_approx_mvbb_grid(gdiam_point *start, int size, int grid_size)
Definition: gdiam.cpp:2742
GTreeDiamAlg::getDiameter
const GPointPair & getDiameter() const
Definition: gdiam.cpp:605
ProjPointSet::compute
void compute()
Definition: gdiam.cpp:2490
GFSPTree::build_node
GFSPTreeNode * build_node(gdiam_point *left, gdiam_point *right)
Definition: gdiam.cpp:281
GTreeDiamAlg::split_pair
void split_pair(GFSPPair &pair, g_heap_pairs_p &heap, double eps)
Definition: gdiam.cpp:1191
GFSPTree::getPoints
const gdiam_point * getPoints() const
Definition: gdiam.cpp:206
GFSPTree::term
void term()
Definition: gdiam.cpp:197
GPointPair::init
void init(const gdiam_point _p, const gdiam_point _q)
Definition: gdiam.h:205
g_heap_pairs_p::pop
void pop()
Definition: gdiam.cpp:864
index
uint8_t index
Definition: EtherCATFrame.h:59
GFSPTreeNode::getBBox
const GBBox & getBBox() const
Definition: gdiam.cpp:155
GFSPTreeNode::get_right
GFSPTreeNode * get_right()
Definition: gdiam.cpp:129
pnt_distance_2d
gdiam_real pnt_distance_2d(gdiam_point_2d p, gdiam_point_2d q)
Definition: gdiam.h:118
GPointPair::p
gdiam_point p
Definition: gdiam.h:198
GFSPTreeNode::isSplitable
bool isSplitable() const
Definition: gdiam.cpp:134
gdiam_bbox::bound
void bound(const gdiam_point pnt)
Definition: gdiam.h:729
GFSPTreeNode::GFSPTreeNode
GFSPTreeNode(gdiam_point *_p_pnt_left, gdiam_point *_p_pnt_right)
Definition: gdiam.cpp:106
CompareByAngle::base
point2d base
Definition: gdiam.cpp:1635
GTreeDiamAlg::size
int size() const
Definition: gdiam.cpp:601
gdiam_approx_mvbb_grid_sample
gdiam_bbox gdiam_approx_mvbb_grid_sample(gdiam_point *start, int size, int grid_size, int sample_size)
Definition: gdiam.cpp:2919
GFSPTree::init
void init(const gdiam_point *_arr, int size)
Definition: gdiam.cpp:169
heap_is_empty
bool heap_is_empty(heap_t *pHeap)
Definition: gdiam.cpp:799
MinAreaRectangle
Definition: gdiam.cpp:2028
compare_pairs
int compare_pairs(void *_a, void *_b)
Definition: gdiam.cpp:807
verify_convex_hull
void verify_convex_hull(vec_point_2d &in, vec_point_2d &ch)
Definition: gdiam.cpp:1736
gdiam_bbox
Definition: gdiam.h:477
GFSPTreeNode::ptr_pnt_right
const gdiam_point * ptr_pnt_right()
Definition: gdiam.cpp:147
GFSPTreeNode::maxDiamProj
gdiam_real maxDiamProj() const
Definition: gdiam.cpp:120
bbox_proj_dist
gdiam_real bbox_proj_dist(const GBBox &bb1, const GBBox &bb2, gdiam_point_cnt proj)
Definition: gdiam.cpp:396
PI
#define PI
Definition: gdiam.cpp:41
pnt_length
gdiam_real pnt_length(const gdiam_point pnt)
Definition: gdiam.h:63
MinAreaRectangle::compute_edge_dir
void compute_edge_dir(int u)
Definition: gdiam.cpp:2044
void addPairHeap(g_heap_pairs_p &heap, GFSPTreeNode *left, GFSPTreeNode *right)
Definition: gdiam.cpp:1009
gdiam_point
gdiam_real * gdiam_point
Definition: gdiam.h:24
c
constexpr T c
Definition: UnscentedKalmanFilterTest.cpp:43
point2d::init
void init(gdiam_point pnt, gdiam_point base_x, gdiam_point base_y)
Definition: gdiam.cpp:1529
GFSPTreeNode::size
int size() const
Definition: gdiam.cpp:73
UMD_SIZE
#define UMD_SIZE
GFSPTree::brute_diameter
gdiam_real brute_diameter(int a_lo, int a_hi, int b_lo, int b_hi, GPointPair &diam, const gdiam_point dir) const
Definition: gdiam.cpp:228
MinAreaRectangle::compute_min_bbox
void compute_min_bbox(vec_point_2d &points, bbox_2d_info &min_bbox, gdiam_real &min_area)
Definition: gdiam.cpp:2343
GBBox::center
void center(gdiam_point out) const
Definition: gdiam.h:329
GBBox::max_coord
const gdiam_real & max_coord(int coord) const
Definition: gdiam.h:374
heap_delete_max
void * heap_delete_max(heap_t *pHeap)
Definition: gdiam.cpp:736
gdiam_bbox::dump
void dump() const
Definition: gdiam.h:636
GFSPPair
Definition: gdiam.cpp:423
eq_real
bool eq_real(gdiam_real a, gdiam_real b)
Definition: gdiam.cpp:2022
GFSPTree::findExtremeInProjection
void findExtremeInProjection(GFSPPair &pair, GPointPair &max_pair_diam)
GPointPair::update_diam
void update_diam(const gdiam_point _p, const gdiam_point _q)
Definition: gdiam.h:260
GFSPPair::init
void init(GFSPTreeNode *_left, GFSPTreeNode *_right)
Definition: gdiam.cpp:430
UDM_SIMPLE
#define UDM_SIMPLE
Definition: gdiam.cpp:573
gdiam_convex_sample
gdiam_point * gdiam_convex_sample(gdiam_point *start, int size, gdiam_bbox &bb, int sample_size)
Definition: gdiam.cpp:2851
HEAP_DELTA
#define HEAP_DELTA
Definition: gdiam.cpp:40
construct_base
void construct_base(gdiam_point pnt1, gdiam_point pnt2, gdiam_point pnt3)
Definition: gdiam.cpp:2426
GFSPTree::brute_diameter
gdiam_real brute_diameter(const gdiam_point *a_lo, const gdiam_point *a_hi, const gdiam_point *b_lo, const gdiam_point *b_hi, GPointPair &diam) const
Definition: gdiam.cpp:241
GFSPTreeNode::nodes_number
int nodes_number() const
Definition: gdiam.cpp:81
pnt_copy
void pnt_copy(gdiam_point_t dest, gdiam_point_t src)
Definition: gdiam.h:83
gdiam_real
double gdiam_real
Definition: gdiam.h:20
GFSPPair::init
void init(GFSPTreeNode *_left, GFSPTreeNode *_right, gdiam_point proj_dir, gdiam_real dist)
Definition: gdiam.cpp:444
GFSPPair::directDiameter
double directDiameter(const GFSPTree &tree, GPointPair &diam) const
Definition: gdiam.cpp:530
pnt_dump
void pnt_dump(gdiam_point_cnt pnt)
Definition: gdiam.h:94
bbox_2d_info::area
gdiam_real area
Definition: gdiam.cpp:1988
heap_top
void * heap_top(heap_t *pHeap)
Definition: gdiam.cpp:730
GTreeDiamAlg::split_pair_proj
void split_pair_proj(GFSPPair &pair, g_heap_pairs_p &heap, double eps, gdiam_point proj)
Definition: gdiam.cpp:1106
gdiam_bbox::combine
void combine(gdiam_point out, double a_coef, double b_coef, double c_coef)
Definition: gdiam.h:597
heap_t::pCompFunc
ptrCompareFunc pCompFunc
Definition: gdiam.cpp:647
armarx::armem::client::util::swap
void swap(SubscriptionHandle &first, SubscriptionHandle &second)
Definition: SubscriptionHandle.cpp:66
GBBox::bound
void bound(const gdiam_point pnt)
Definition: gdiam.h:379
GFSPPair::get_right
GFSPTreeNode * get_right()
Definition: gdiam.cpp:502
computeExtremePair
void computeExtremePair(const gdiam_point *arr, int size, int dim, GPointPair &pp)
Definition: gdiam.cpp:876
GFSPPair::dump
void dump() const
Definition: gdiam.cpp:507
GFSPTreeNode::ptr_pnt_left
const gdiam_point * ptr_pnt_left()
Definition: gdiam.cpp:139
GFSPTree::brute_diameter
gdiam_real brute_diameter(const gdiam_point *a_lo, const gdiam_point *a_hi, const gdiam_point *b_lo, const gdiam_point *b_hi, GPointPair &diam, const gdiam_point dir) const
Definition: gdiam.cpp:255
gdiam_convert
gdiam_point * gdiam_convert(gdiam_real *start, int size)
Definition: gdiam.cpp:1302
GFSPTree::split_array
int split_array(gdiam_point *left, gdiam_point *right, int dim, const gdiam_real &val)
Definition: gdiam.cpp:314
armarx::ctrlutil::a
double a(double t, double a0, double j)
Definition: CtrlUtil.h:45
GTreeDiamAlg::nodes_number
int nodes_number() const
Definition: gdiam.cpp:610
g_heap_pairs_p::size
int size() const
Definition: gdiam.cpp:856
armarx::abs
std::vector< T > abs(const std::vector< T > &v)
Definition: VectorHelpers.h:253
GTreeDiamAlg::term
void term()
Definition: gdiam.cpp:596
GFSPPair::maxDiam
double maxDiam() const
Definition: gdiam.cpp:551
GFSPTree::brute_diameter
gdiam_real brute_diameter(int a_lo, int a_hi, int b_lo, int b_hi, GPointPair &diam) const
Definition: gdiam.cpp:216
GTreeDiamAlg::diameter
double diameter()
Definition: gdiam.cpp:631
ProjPointSet::init
void init()
Definition: gdiam.cpp:2592
GBBox::get_diam
gdiam_real get_diam() const
Definition: gdiam.h:419
gdiam_point_2d
gdiam_real * gdiam_point_2d
Definition: gdiam.h:23
bbox_2d_info::dump
void dump()
Definition: gdiam.cpp:1998
GFSPTreeNode::get_left
GFSPTreeNode * get_left()
Definition: gdiam.cpp:125
GTreeDiamAlg::compute_by_heap_proj
void compute_by_heap_proj(double eps, gdiam_point proj)
Definition: gdiam.cpp:950
dump_points
void dump_points(gdiam_point *in_arr, int size)
Definition: gdiam.cpp:2357
GBBox::init
void init(const GBBox &a, const GBBox &b)
Definition: gdiam.h:299
heap_insert
void heap_insert(heap_t *pHeap, void *pData)
Definition: gdiam.cpp:699
GFSPTree::nodes_number
int nodes_number() const
Definition: gdiam.cpp:211
GFSPPair::directDiameter
double directDiameter(const GFSPTree &tree, GPointPair &diam, const gdiam_point dir) const
Definition: gdiam.cpp:540
pnt_distance
gdiam_real pnt_distance(gdiam_point p, gdiam_point q)
Definition: gdiam.h:128
dump
void dump(vec_point_2d &vec)
Definition: gdiam.cpp:1717
point2d::dump
void dump() const
Definition: gdiam.cpp:1544
g_heap_pairs_p::push
void push(GFSPPair &pair)
Definition: gdiam.cpp:848
MinAreaRectangle::get_angles
void get_angles(int ind, gdiam_real &angle_1, gdiam_real &angle_2, gdiam_real &angle_3)
Definition: gdiam.cpp:2246
pnt_cross_prod
void pnt_cross_prod(const gdiam_point a, const gdiam_point b, const gdiam_point out)
Definition: gdiam.h:108
bbox_2d_info::vertices
gdiam_point_2d_t vertices[4]
Definition: gdiam.cpp:1987
GBBox
Definition: gdiam.h:292
gdiam_bbox::get_dir
gdiam_point get_dir(int ind)
Definition: gdiam.h:575
EPS
#define EPS
Definition: gdiam.cpp:2021
GFSPTree::split_node
void split_node(GFSPTreeNode *node)
Definition: gdiam.cpp:340
max
T max(T t1, T t2)
Definition: gdiam.h:48
point2d::equal_real
bool equal_real(const point2d &pnt) const
Definition: gdiam.cpp:1552
gdiam_approx_const_mvbb
gdiam_bbox gdiam_approx_const_mvbb(gdiam_point *start, int size, gdiam_real eps, GBBox *p_ap_bbox)
Definition: gdiam.cpp:1349
GFSPPair::getProjectionError
double getProjectionError() const
Definition: gdiam.cpp:483
GFSPTreeNode::getCenter
gdiam_point getCenter() const
Definition: gdiam.cpp:77
bbox_vertex
void bbox_vertex(const GBBox &bb, gdiam_point_t &ver, int vert)
Definition: gdiam.cpp:384
GBBox::min_coord
const gdiam_real & min_coord(int coord) const
Definition: gdiam.h:370
gdiam_approx_diam_pair_UDM
GPointPair gdiam_approx_diam_pair_UDM(gdiam_real *start, int size, gdiam_real eps)
Definition: gdiam.cpp:1335
GFSPTree
Definition: gdiam.cpp:162
pnt_dot_prod
gdiam_real pnt_dot_prod(gdiam_point_cnt a, gdiam_point_cnt b)
Definition: gdiam.h:100
q
#define q
point2d::src
gdiam_point src
Definition: gdiam.cpp:1527
GBBox::getLongestDim
int getLongestDim() const
Definition: gdiam.h:400
GPointPair::q
gdiam_point q
Definition: gdiam.h:198
GFSPTreeNode::maxDiam
gdiam_real maxDiam() const
Definition: gdiam.cpp:116
void pnt_scale_and_add(gdiam_point dest, gdiam_real coef, gdiam_point_cnt vec)
Definition: gdiam.h:185
gdiam_point_cnt
const gdiam_real * gdiam_point_cnt
Definition: gdiam.h:25
get_min_point
point2d_ptr get_min_point(vec_point_2d &in, int &index)
Definition: gdiam.cpp:1692
gdiam_bbox::get_normalized_coordinates
void get_normalized_coordinates(gdiam_point in, gdiam_point out)
Definition: gdiam.h:779
HEAP_LIMIT
#define HEAP_LIMIT
Definition: gdiam.cpp:39
gdiam.h
GPointPair::distance
gdiam_real distance
Definition: gdiam.h:197
gdiam_bbox::set_third_dim_longest
void set_third_dim_longest()
Definition: gdiam.h:547
heap_t::max_size
int max_size
Definition: gdiam.cpp:646
convex_hull
void convex_hull(vec_point_2d &in, vec_point_2d &out)
Definition: gdiam.cpp:1882
GTreeDiamAlg
Definition: gdiam.cpp:576
point2d::x
gdiam_real x
Definition: gdiam.cpp:1526
point2d::dist
gdiam_real dist(const point2d &pnt)
Definition: gdiam.cpp:1538
MinAreaRectangle::find_vertex
int find_vertex(int start_vertex, double trg_angle)
Definition: gdiam.cpp:2088
armarx::ctrlutil::v
double v(double t, double v0, double a0, double j)
Definition: CtrlUtil.h:39
heap_t::pArr
voidPtr_t * pArr
Definition: gdiam.cpp:645
point2d
Definition: gdiam.cpp:1523
gdiam_point_2d_t
gdiam_real gdiam_point_2d_t[2]
Definition: gdiam.h:22
pnt_isEqual
bool pnt_isEqual(const gdiam_point p, const gdiam_point q)
Definition: gdiam.h:176
gdiam_generate_orthonormal_base
void gdiam_generate_orthonormal_base(gdiam_point in, gdiam_point out1, gdiam_point out2)
Definition: gdiam.cpp:1464
MinAreaRectangle::dump_ch
void dump_ch()
Definition: gdiam.cpp:2035
gdiam_exchange
void gdiam_exchange(T &a, T &b)
Definition: gdiam.h:55
heap_init
void heap_init(heap_t *pHeap, ptrCompareFunc _pCompFunc)
Definition: gdiam.cpp:650
gdiam_point_t
gdiam_real gdiam_point_t[GDIAM_DIM]
Definition: gdiam.h:21
GFSPPair::minAprxDiam
double minAprxDiam() const
Definition: gdiam.cpp:555
gdiam_approx_diam
GPointPair gdiam_approx_diam(gdiam_point *start, int size, gdiam_real eps)
Definition: gdiam.cpp:1268
GBBox::volume
gdiam_real volume() const
Definition: gdiam.h:337
angle
double angle(const Point &a, const Point &b, const Point &c)
Definition: point.hpp:100
heap_t
Definition: gdiam.cpp:643
ProjPointSet
Definition: gdiam.cpp:2443
GFSPTreeNode::dump
void dump() const
Definition: gdiam.cpp:68
pnt_normalize
void pnt_normalize(gdiam_point pnt)
Definition: gdiam.h:69
GFSPTree::getRoot
GFSPTreeNode * getRoot()
Definition: gdiam.cpp:276
UDM_BIG_SAMPLE
#define UDM_BIG_SAMPLE
Definition: gdiam.cpp:574
CompareByAngle::operator()
bool operator()(const point2d_ptr &a, const point2d_ptr &b)
Definition: gdiam.cpp:1637
GFSPTreeNode
Definition: gdiam.cpp:57
g_heap_pairs_p::~g_heap_pairs_p
~g_heap_pairs_p()
Definition: gdiam.cpp:838
GTreeDiamAlg::GTreeDiamAlg
GTreeDiamAlg(gdiam_point *arr, int size, int mode)
Definition: gdiam.cpp:586
ProjPointSet::get_bbox
gdiam_bbox & get_bbox()
Definition: gdiam.cpp:2579
MinAreaRectangle::compute_min_bbox_inner
void compute_min_bbox_inner(bbox_2d_info &min_bbox, gdiam_real &min_area)
Definition: gdiam.cpp:2281
voidPtr_t
void * voidPtr_t
Definition: gdiam.cpp:641
MinAreaRectangle::get_bbox
void get_bbox(int a_ind, gdiam_real a_angle, int b_ind, gdiam_real b_angle, int c_ind, gdiam_real c_angle, int d_ind, gdiam_real d_angle, bbox_2d_info &bbox, gdiam_real &area)
Definition: gdiam.cpp:2209
Left
bool Left(const point2d &a, const point2d &b, const point2d &c)
Definition: gdiam.cpp:1617
heap_term
void heap_term(heap_t *pHeap)
Definition: gdiam.cpp:686
armarx::ctrlutil::s
double s(double t, double s0, double v0, double a0, double j)
Definition: CtrlUtil.h:33
heap_t::curr_size
int curr_size
Definition: gdiam.cpp:646
memory.h
ProjPointSet::term
void term()
Definition: gdiam.cpp:2584
ProjPointSet::~ProjPointSet
~ProjPointSet()
Definition: gdiam.cpp:2454
point2d_ptr
point2d * point2d_ptr
Definition: gdiam.cpp:1559
GFSPTree::getPoint
gdiam_point getPoint(int ind) const
Definition: gdiam.cpp:271
GPointPair
Definition: gdiam.h:194
GFSPPair::isBelowThreshold
bool isBelowThreshold(int threshold) const
Definition: gdiam.cpp:515
AreaSign
int AreaSign(const point2d &a, const point2d &b, const point2d &c)
Definition: gdiam.cpp:1602
pnt_init_normalize
void pnt_init_normalize(gdiam_point pnt, gdiam_real x, gdiam_real y, gdiam_real z)
Definition: gdiam.h:164
bbox_2d_info::get_dir
void get_dir(int ind, gdiam_point_2d out)
Definition: gdiam.cpp:1990
gdiam_bbox::init
void init(const GBBox &bb)
Definition: gdiam.h:667
GFSPTree::terminate
static void terminate(GFSPTreeNode *node)
Definition: gdiam.cpp:183
GTreeDiamAlg::compute_by_heap
void compute_by_heap(double eps)
Definition: gdiam.cpp:901