BRL-CAD
chull3d.cpp
Go to the documentation of this file.
1 /*
2  * Ken Clarkson wrote this. Copyright (c) 1995 by AT&T..
3  * Permission to use, copy, modify, and distribute this software for any
4  * purpose without fee is hereby granted, provided that this entire notice
5  * is included in all copies of any software which is or includes a copy
6  * or modification of this software and in all copies of the supporting
7  * documentation for such software.
8  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
9  * WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
10  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
11  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
12  */
13 
14 #include "common.h"
15 
16 #include <map>
17 
18 #include <float.h>
19 #include <locale.h>
20 #include <math.h>
21 #include <stdarg.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <string.h>
25 #include <ctype.h>
26 #include <time.h>
27 
28 /*#include "chull3d.h"*/
29 #include "bu/log.h"
30 #include "bu/ptbl.h"
31 #include "bu/malloc.h"
32 #include "bn/chull.h"
33 #include "bn/randmt.h"
34 
35 #define BLOCKSIZE 100000
36 #define MAXDIM 3
37 
38 /* It is faster, when dealing with large point sets, to
39  * repeatedly calculate the convex hull for subsets of
40  * the total set and then fine the hull of that set of
41  * hulls. These parameters control that process - see
42  * chull3d_intermediate_set and bn-3d_chull. In the
43  * worst case this approach will slow the calculation,
44  * but so long as some significant percentage of the
45  * bot's vertices are not part of the hull it should
46  * help */
47 #define CHULL3D_COUNT_THRESHOLD 1000
48 #define CHULL3D_DELTA_THRESHOLD 200
49 #define CHULL3D_SUBSET_SIZE 275
50 
51 typedef double Coord;
52 typedef Coord* point;
53 typedef point site;
54 typedef Coord* normalp;
56 
57 typedef struct basis_s {
58  struct basis_s *next; /* free list */
59  int ref_count; /* storage management */
60  int lscale; /* the log base 2 of total scaling of vector */
61  Coord sqa, sqb; /* sums of squared norms of a part and b part */
62  Coord vecs[1]; /* the actual vectors, extended by malloc'ing bigger */
63 } basis_s;
64 
65 typedef struct neighbor {
66  site vert; /* vertex of simplex */
67  struct simplex *simp; /* neighbor sharing all vertices but vert */
68  basis_s *basis; /* derived vectors */
69 } neighbor;
70 
71 typedef struct simplex {
72  struct simplex *next; /* free list */
73  long visit; /* number of last site visiting this simplex */
74  short mark;
75  basis_s* normal; /* normal vector pointing inward */
76  neighbor peak; /* if null, remaining vertices give facet */
77  neighbor neigh[1]; /* neighbors of simplex */
78 } simplex;
79 
80 typedef struct fg_node fg;
81 typedef struct tree_node Tree;
82 struct tree_node {
85  int size; /* maintained to be the number of nodes rooted here */
86  fg *fgs;
87  Tree *next; /* freelist */
88 };
89 
90 typedef struct fg_node {
92  double dist, vol; /* of Voronoi face dual to this */
93  fg *next; /* freelist */
94  short mark;
95  int ref_count;
96 } fg_node;
97 
98 
99 typedef site gsitef(void *);
100 typedef long site_n(void *, site);
101 
102 /* We use this instead of globals */
103 struct chull3d_data {
104  size_t simplex_size;
106  size_t basis_s_size;
108  size_t Tree_size;
110  int pdim; /* point dimension */
114  double Huge;
117  float b_err_min;
122  int *B;
126  double mult_up;
129  site p; /* the current site */
130  long pnum;
131  int rdim; /* region dimension: (max) number of sites specifying region */
132  int cdim; /* number of sites currently specifying region */
133  int site_size; /* size of malloc needed for a site */
134  int point_size; /* size of malloc needed for a point */
135  short *mi;
136  short *mo;
137  long *shufmat;
138 
139  /* these were static variables in functions */
146  int lscale;
147  double max_scale;
148  double ldetbound;
149  double Sb;
152  int fg_vn;
153  long vnum;
156  long ss;
157  long s_num;
158  std::map<long *, int> *point_lookup;
160 
161  /* libbn style inputs and outputs */
162  const point_t *input_vert_array;
167  point_t center_pnt;
168 
169  /* Output containers */
170  point_t *vert_array;
171  int *vert_cnt;
172  int *faces;
173  int *face_cnt;
174 };
175 
176 
177 #define CHULL3D_VA(x) ((x)->vecs+cdata->rdim)
178 #define CHULL3D_VB(x) ((x)->vecs)
179 #define CHULL3D_trans(z,p,q) {int ti; for (ti=0;ti<cdata->pdim;ti++) z[ti+cdata->rdim] = z[ti] = p[ti] - q[ti];}
180 #define CHULL3D_DELIFT 0
181 
182 #define CHULL3D_push(x) *(cdata->st+tms++) = x;
183 #define CHULL3D_pop(x) x = *(cdata->st + --tms);
184 
185 #define CHULL3D_NEWL(cdata, list, X, p) \
186 { \
187  p = list ? list : chull3d_new_block_##X(cdata, 1); \
188  /*assert(p);*/ \
189  list = p->next; \
190 }
191 
192 #define CHULL3D_NEWLRC(cdata, list, X, p) \
193 { \
194  p = list ? list : chull3d_new_block_##X(cdata, 1); \
195  /*assert(p);*/ \
196  list = p->next; \
197  p->ref_count = 1; \
198 }
199 
200 #define CHULL3D_NULLIFY(cdata, X,v) { \
201  if ((v) && --(v)->ref_count == 0) { \
202  memset((v),0,cdata->X##_size); \
203  (v)->next = cdata->X##_list; \
204  cdata->X##_list = v; \
205  } \
206  v = NULL;\
207 }
208 
209 #define CHULL3D_copy_simp(cdata, cnew, s) \
210 { \
211  int mi; \
212  neighbor *mrsn; \
213  CHULL3D_NEWL(cdata, cdata->simplex_list, simplex, cnew); \
214  memcpy(cnew,s,cdata->simplex_size); \
215  for (mi=-1,mrsn=s->neigh-1;mi<(cdata->cdim);mi++,mrsn++) \
216  if (mrsn->basis) mrsn->basis->ref_count++; \
217 }
218 
219 HIDDEN simplex *
220 chull3d_new_block_simplex(struct chull3d_data *cdata, int make_blocks)
221 {
222  int i;
223  simplex *xlm, *xbt;
224  if (make_blocks && cdata) {
225  xbt = cdata->simplex_block_table[(cdata->num_simplex_blocks)++] = (simplex*)malloc(cdata->input_vert_cnt * cdata->simplex_size);
226  memset(xbt,0,cdata->input_vert_cnt * cdata->simplex_size);
227  xlm = ((simplex*) ( (char*)xbt + (cdata->input_vert_cnt) * cdata->simplex_size));
228  for (i=0; i<cdata->input_vert_cnt; i++) {
229  xlm = ((simplex*) ( (char*)xlm + ((-1)) * cdata->simplex_size));
230  xlm->next = cdata->simplex_list;
231  cdata->simplex_list = xlm;
232  }
233  return cdata->simplex_list;
234  }
235  for (i=0; i<cdata->num_simplex_blocks; i++) free(cdata->simplex_block_table[i]);
236  cdata->num_simplex_blocks = 0;
237  if (cdata) cdata->simplex_list = 0;
238  return 0;
239 }
240 
241 HIDDEN void
243 {
244  chull3d_new_block_simplex(cdata, 0);
245 }
246 
247 HIDDEN basis_s *
248 chull3d_new_block_basis_s(struct chull3d_data *cdata, int make_blocks)
249 {
250  int i;
251  basis_s *xlm, *xbt;
252  if (make_blocks && cdata) {
253  xbt = cdata->basis_s_block_table[(cdata->num_basis_s_blocks)++] = (basis_s*)malloc(cdata->input_vert_cnt * cdata->basis_s_size);
254  memset(xbt,0,cdata->input_vert_cnt * cdata->basis_s_size);
255  xlm = ((basis_s*) ( (char*)xbt + (cdata->input_vert_cnt) * cdata->basis_s_size));
256  for (i=0; i<cdata->input_vert_cnt; i++) {
257  xlm = ((basis_s*) ( (char*)xlm + ((-1)) * cdata->basis_s_size));
258  xlm->next = cdata->basis_s_list;
259  cdata->basis_s_list = xlm;
260  }
261  return cdata->basis_s_list;
262  }
263  for (i=0; i<cdata->num_basis_s_blocks; i++) free(cdata->basis_s_block_table[i]);
264  cdata->num_basis_s_blocks = 0;
265  if (cdata) cdata->basis_s_list = 0;
266  return 0;
267 }
268 
269 HIDDEN void
271 {
272  chull3d_new_block_basis_s(cdata, 0);
273 }
274 
275 
276 HIDDEN Tree *
277 chull3d_new_block_Tree(struct chull3d_data *cdata, int make_blocks) {
278  int i;
279  Tree *xlm, *xbt;
280  if (make_blocks && cdata) {
281  xbt = cdata->Tree_block_table[cdata->num_Tree_blocks++] = (Tree*)malloc(cdata->input_vert_cnt * cdata->Tree_size);
282  memset(xbt,0,cdata->input_vert_cnt * cdata->Tree_size);
283  xlm = ((Tree*) ( (char*)xbt + (cdata->input_vert_cnt) * cdata->Tree_size));
284  for (i=0; i<cdata->input_vert_cnt; i++) {
285  xlm = ((Tree*) ( (char*)xlm + ((-1)) * cdata->Tree_size));
286  xlm->next = cdata->Tree_list;
287  cdata->Tree_list = xlm;
288  }
289  return cdata->Tree_list;
290  }
291  for (i=0; i<cdata->num_Tree_blocks; i++) free(cdata->Tree_block_table[i]);
292  cdata->num_Tree_blocks = 0;
293  if (cdata) cdata->Tree_list = 0;
294  return 0;
295 }
296 
297 HIDDEN void
299 {
300  chull3d_new_block_Tree(cdata, 0);
301 }
302 
303 /* ch.c : numerical functions for hull computation */
304 
307 {
308  int i;
309  Coord sum = 0;
310  for (i=0;i<cdata->rdim;i++) sum += x[i] * y[i];
311  return sum;
312 }
313 
314  HIDDEN Coord
316 {
317  int i;
318  Coord sum = 0;
319  for (i=0;i<cdata->pdim;i++) sum += x[i] * y[i];
320  return sum;
321 }
322 
323  HIDDEN Coord
325 {
326  int i;
327  Coord sum = 0;
328  for (i=0;i<cdata->rdim;i++) sum += x[i] * x[i];
329  return sum;
330 }
331 
332  HIDDEN void
334 {
335  int i;
336  for (i=0;i<cdata->rdim;i++) {
337  *y++ += a * *x++;
338  }
339 }
340 
341  HIDDEN void
342 chull3d_Vec_scale(struct chull3d_data *UNUSED(cdata), int n, Coord a, Coord *x)
343 {
344  register Coord *xx = x;
345  register Coord *xend = xx + n;
346 
347  while (xx!=xend) {
348  *xx *= a;
349  xx++;
350  }
351 }
352 
353 
354 /* amount by which to scale up vector, for reduce_inner */
355  HIDDEN double
356 chull3d_sc(struct chull3d_data *cdata, basis_s *v,simplex *s, int k, int j)
357 {
358  double labound;
359 
360  if (j<10) {
361  labound = logb(v->sqa)/2;
362  cdata->max_scale = cdata->exact_bits - labound - 0.66*(k-2)-1 -CHULL3D_DELIFT;
363  if (cdata->max_scale<1) {
364  cdata->max_scale = 1;
365  }
366 
367  if (j==0) {
368  int i;
369  neighbor *sni;
370  basis_s *snib;
371 
372  cdata->ldetbound = CHULL3D_DELIFT;
373 
374  cdata->Sb = 0;
375  for (i=k-1,sni=s->neigh+k-1;i>0;i--,sni--) {
376  snib = sni->basis;
377  cdata->Sb += snib->sqb;
378  cdata->ldetbound += logb(snib->sqb)/2 + 1;
379  cdata->ldetbound -= snib->lscale;
380  }
381  }
382  }
383  if (cdata->ldetbound - v->lscale + logb(v->sqb)/2 + 1 < 0) {
384  return 0;
385  } else {
386  double dlscale = floor(logb(2*cdata->Sb/(v->sqb + v->sqa*cdata->b_err_min)))/2;
387  cdata->lscale = (int)dlscale;
388  if (cdata->lscale > cdata->max_scale) {
389  dlscale = floor(cdata->max_scale);
390  cdata->lscale = (int)dlscale;
391  } else if (cdata->lscale<0) cdata->lscale = 0;
392  v->lscale += cdata->lscale;
393  return ( ((cdata->lscale)<20) ? 1<<(cdata->lscale) : ldexp(1,(cdata->lscale)) );
394  }
395 }
396 
397 
398  HIDDEN int
400 {
401  point va = CHULL3D_VA(v);
402  point vb = CHULL3D_VB(v);
403  int i,j;
404  double dd;
405  double scale;
406  basis_s *snibv;
407  neighbor *sni;
408 
409  v->sqa = v->sqb = chull3d_Norm2(cdata, vb);
410  if (k<=1) {
411  memcpy(vb,va,cdata->basis_vec_size);
412  return 1;
413  }
414  for (j=0;j<250;j++) {
415 
416  memcpy(vb,va,cdata->basis_vec_size);
417  for (i=k-1,sni=s->neigh+k-1;i>0;i--,sni--) {
418  snibv = sni->basis;
419  dd = -chull3d_Vec_dot(cdata, CHULL3D_VB(snibv),vb)/ snibv->sqb;
420  chull3d_Ax_plus_y(cdata, dd, CHULL3D_VA(snibv), vb);
421  }
422  v->sqb = chull3d_Norm2(cdata, vb);
423  v->sqa = chull3d_Norm2(cdata, va);
424 
425  if (2*v->sqb >= v->sqa) {cdata->B[j]++; return 1;}
426 
427  chull3d_Vec_scale(cdata, cdata->rdim,scale = chull3d_sc(cdata,v,s,k,j),va);
428 
429  for (i=k-1,sni=s->neigh+k-1;i>0;i--,sni--) {
430  snibv = sni->basis;
431  dd = -chull3d_Vec_dot(cdata, CHULL3D_VB(snibv),va)/snibv->sqb;
432  dd = floor(dd+0.5);
433  chull3d_Ax_plus_y(cdata, dd, CHULL3D_VA(snibv), va);
434  }
435  }
436  return 0;
437 }
438 
439 HIDDEN int
440 chull3d_reduce(struct chull3d_data *cdata, basis_s **v, point p, simplex *s, int k)
441 {
442  point z;
443  point tt = s->neigh[0].vert;
444 
445  if (!*v) CHULL3D_NEWLRC(cdata, cdata->basis_s_list,basis_s,(*v))
446  else (*v)->lscale = 0;
447  z = CHULL3D_VB(*v);
448  CHULL3D_trans(z,p,tt);
449  return chull3d_reduce_inner(cdata,*v,s,k);
450 }
451 
452  HIDDEN void
454 {
455  int k=1;
456  neighbor *sn = s->neigh+1,
457  *sn0 = s->neigh;
458  if (!sn0->basis) {
459  sn0->basis = cdata->tt_basis;
460  cdata->tt_basis->ref_count++;
461  } else while (k < (cdata->cdim) && sn->basis) {k++;sn++;}
462  while (k<(cdata->cdim)) {
463  CHULL3D_NULLIFY(cdata, basis_s,sn->basis);
464  chull3d_reduce(cdata, &sn->basis,sn->vert,s,k);
465  k++; sn++;
466  }
467 }
468 
469  HIDDEN int
471 {
472  if (!cdata->p_neigh.basis) cdata->p_neigh.basis = (basis_s*) malloc(cdata->basis_s_size);
473  cdata->p_neigh.vert = p;
474  (cdata->cdim)++;
475  root->neigh[(cdata->cdim)-1].vert = root->peak.vert;
476  CHULL3D_NULLIFY(cdata, basis_s,root->neigh[(cdata->cdim)-1].basis);
477  chull3d_get_basis_sede(cdata, root);
478  chull3d_reduce(cdata, &cdata->p_neigh.basis,p,root,(cdata->cdim));
479  /*if (cdata->p_neigh.basis->sqa != 0) return 1;*/
480  if (!ZERO(cdata->p_neigh.basis->sqa)) return 1;
481  (cdata->cdim)--;
482  return 0;
483 }
484 
485 HIDDEN int chull3d_sees(struct chull3d_data *, site, simplex *);
486 
487 HIDDEN void
489 {
490  neighbor *rn;
491  int i,j;
492 
493  chull3d_get_basis_sede(cdata, s);
494  if (cdata->rdim==3 && (cdata->cdim)==3) {
495  point c;
496  point a = CHULL3D_VB(s->neigh[1].basis);
497  point b = CHULL3D_VB(s->neigh[2].basis);
498  CHULL3D_NEWLRC(cdata, cdata->basis_s_list, basis_s,s->normal);
499  c = CHULL3D_VB(s->normal);
500  c[0] = a[1]*b[2] - a[2]*b[1];
501  c[1] = a[2]*b[0] - a[0]*b[2];
502  c[2] = a[0]*b[1] - a[1]*b[0];
503  s->normal->sqb = chull3d_Norm2(cdata, c);
504  for (i=(cdata->cdim)+1,rn = cdata->ch_root->neigh+(cdata->cdim)-1; i; i--, rn--) {
505  for (j = 0; j<(cdata->cdim) && rn->vert != s->neigh[j].vert;j++);
506  if (j<(cdata->cdim)) continue;
507  if (rn->vert==cdata->hull_infinity) {
508  if (c[2] > -cdata->b_err_min) continue;
509  } else if (!chull3d_sees(cdata, rn->vert,s)) continue;
510  c[0] = -c[0]; c[1] = -c[1]; c[2] = -c[2];
511  break;
512  }
513  return;
514  }
515 
516  for (i=(cdata->cdim)+1,rn = cdata->ch_root->neigh+(cdata->cdim)-1; i; i--, rn--) {
517  for (j = 0; j<(cdata->cdim) && rn->vert != s->neigh[j].vert;j++);
518  if (j<(cdata->cdim)) continue;
519  chull3d_reduce(cdata, &s->normal,rn->vert,s,(cdata->cdim));
520  if (!ZERO(s->normal->sqb)) break;
521  }
522 }
523 
524 
525 HIDDEN void chull3d_get_normal(struct chull3d_data *cdata, simplex *s) {chull3d_get_normal_sede(cdata,s); return;}
526 
527 HIDDEN int
528 chull3d_sees(struct chull3d_data *cdata, site p, simplex *s) {
529 
530  point tt,zz;
531  double dd,dds;
532  int i;
533 
534  if (!cdata->b) cdata->b = (basis_s*)malloc(cdata->basis_s_size);
535  else cdata->b->lscale = 0;
536  zz = CHULL3D_VB(cdata->b);
537  if ((cdata->cdim)==0) return 0;
538  if (!s->normal) {
539  chull3d_get_normal_sede(cdata, s);
540  for (i=0;i<(cdata->cdim);i++) CHULL3D_NULLIFY(cdata, basis_s,s->neigh[i].basis);
541  }
542  tt = s->neigh[0].vert;
543  CHULL3D_trans(zz,p,tt);
544  for (i=0;i<3;i++) {
545  dd = chull3d_Vec_dot(cdata, zz,s->normal->vecs);
546  if (ZERO(dd)) {
547  return 0;
548  }
549  dds = dd*dd/s->normal->sqb/chull3d_Norm2(cdata, zz);
550  if (dds > cdata->b_err_min_sq) return (dd<0);
551  chull3d_get_basis_sede(cdata, s);
552  chull3d_reduce_inner(cdata,cdata->b,s,(cdata->cdim));
553  }
554  return 0;
555 }
556 
557 HIDDEN void
559 {
563 }
564 
565 HIDDEN void
566 chull3d_visit_fg_i(struct chull3d_data *cdata, void (*v_fg)(struct chull3d_data *,Tree *, int, int),
567  Tree *t, int depth, int vn, int boundary)
568 {
569  int boundaryc = boundary;
570 
571  if (!t) return;
572 
573  //assert(t->fgs);
574  if (t->fgs->mark!=vn) {
575  t->fgs->mark = vn;
576  if (t->key!=cdata->hull_infinity && !cdata->mo[cdata->site_num((void *)cdata, t->key)]) boundaryc = 0;
577  v_fg(cdata, t,depth, boundaryc);
578  chull3d_visit_fg_i(cdata, v_fg, t->fgs->facets,depth+1, vn, boundaryc);
579  }
580  chull3d_visit_fg_i(cdata, v_fg, t->left,depth,vn, boundary);
581  chull3d_visit_fg_i(cdata, v_fg, t->right,depth,vn,boundary);
582 }
583 
584 HIDDEN void
585 chull3d_visit_fg(struct chull3d_data *cdata, fg *faces_gr, void (*v_fg)(struct chull3d_data *, Tree *, int, int))
586 {
587  chull3d_visit_fg_i(cdata, v_fg, faces_gr->facets, 0, ++(cdata->fg_vn), 1);
588 }
589 
590 HIDDEN int
591 chull3d_visit_fg_i_far(struct chull3d_data *cdata, void (*v_fg)(struct chull3d_data *, Tree *, int),
592  Tree *t, int depth, int vn)
593 {
594  int nb = 0;
595 
596  if (!t) return 0;
597 
598  //assert(t->fgs);
599  if (t->fgs->mark!=vn) {
600  t->fgs->mark = vn;
601  nb = (t->key==cdata->hull_infinity) || cdata->mo[cdata->site_num((void *)cdata, t->key)];
602  if (!nb && !chull3d_visit_fg_i_far(cdata, v_fg, t->fgs->facets,depth+1,vn))
603  v_fg(cdata, t,depth);
604  }
605  nb = chull3d_visit_fg_i_far(cdata, v_fg, t->left,depth,vn) || nb;
606  nb = chull3d_visit_fg_i_far(cdata, v_fg, t->right,depth,vn) || nb;
607  return nb;
608 }
609 
610 HIDDEN void
611 chull3d_visit_fg_far(struct chull3d_data *cdata, fg *faces_gr, void (*v_fg)(struct chull3d_data *, Tree *, int))
612 {
613  chull3d_visit_fg_i_far(cdata, v_fg,faces_gr->facets, 0, --(cdata->fg_vn));
614 }
615 
616 /* hull.c : "combinatorial" functions for hull computation */
617 
618 typedef int chull3d_test_func(struct chull3d_data *, simplex *, int, void *);
619 typedef void* chull3d_visit_func(struct chull3d_data *, simplex *, void *);
620 
621 /* starting at s, visit simplices t such that test(s,i,0) is true,
622  * and t is the i'th neighbor of s;
623  * apply visit function to all visited simplices;
624  * when visit returns nonNULL, exit and return its value */
625 HIDDEN void *
627 {
628  neighbor *sn;
629  void *v;
630  simplex *t;
631  int i;
632  long tms = 0;
633 
634  (cdata->vnum)--;
635  if (s) CHULL3D_push(s);
636  while (tms) {
637  if (tms>cdata->ss) {
638  cdata->st=(simplex**)realloc(cdata->st,((cdata->ss+=cdata->ss)+MAXDIM+1)*sizeof(simplex*));
639  }
640  CHULL3D_pop(t);
641  if (!t || t->visit == cdata->vnum) continue;
642  t->visit = cdata->vnum;
643  if ((v=(*visit)(cdata,t,0))) {return v;}
644  for (i=-1,sn = t->neigh-1;i<(cdata->cdim);i++,sn++)
645  if ((sn->simp->visit != cdata->vnum) && sn->simp && test(cdata,t,i,0))
646  CHULL3D_push(sn->simp);
647  }
648  return NULL;
649 }
650 
651 /* visit the whole triangulation */
652 HIDDEN int chull3d_truet(struct chull3d_data *UNUSED(cdata), simplex *UNUSED(s), int UNUSED(i), void *UNUSED(dum)) {return 1;}
653 HIDDEN void *
655 {
656  return chull3d_visit_triang_gen(cdata, root, visit, chull3d_truet);
657 }
658 
659 
660 HIDDEN int chull3d_hullt(struct chull3d_data *UNUSED(cdata), simplex *UNUSED(s), int i, void *UNUSED(dummy)) {return i>-1;}
661 HIDDEN void *chull3d_facet_test(struct chull3d_data *UNUSED(cdata), simplex *s, void *UNUSED(dummy)) {return (!s->peak.vert) ? s : NULL;}
662 HIDDEN void *
664 {
665  return chull3d_visit_triang_gen(cdata, (struct simplex *)chull3d_visit_triang(cdata, root, &chull3d_facet_test), visit, chull3d_hullt);
666 }
667 
668 HIDDEN void *
669 chull3d_check_simplex(struct chull3d_data *cdata, simplex *s, void *UNUSED(dum))
670 {
671  int i,j,k,l;
672  neighbor *sn, *snn, *sn2;
673  simplex *sns;
674  site vn;
675 
676  for (i=-1,sn=s->neigh-1;i<(cdata->cdim);i++,sn++) {
677  sns = sn->simp;
678  if (!sns) {
679  return s;
680  }
681  if (!s->peak.vert && sns->peak.vert && i!=-1) {
682  exit(1);
683  }
684  for (j=-1,snn=sns->neigh-1; j<(cdata->cdim) && snn->simp!=s; j++,snn++);
685  if (j==(cdata->cdim)) {
686  exit(1);
687  }
688  for (k=-1,snn=sns->neigh-1; k<(cdata->cdim); k++,snn++){
689  vn = snn->vert;
690  if (k!=j) {
691  for (l=-1,sn2 = s->neigh-1;
692  l<(cdata->cdim) && sn2->vert != vn;
693  l++,sn2++);
694  if (l==(cdata->cdim)) {
695  exit(1);
696  }
697  }
698  }
699  }
700  return NULL;
701 }
702 
703 HIDDEN void *
704 chull3d_collect_hull_pnts(struct chull3d_data *cdata, simplex *s, void *UNUSED(p)) {
705 
706  point v[MAXDIM];
707  int j;
708  long int ip[3];
709  /*point_t p1, p2, p3;*/
710  const point_t *pp[3];
711  int f;
712 
713  if (!s) return NULL;
714 
715  for (j=0;j<(cdata->cdim);j++) v[j] = s->neigh[j].vert;
716  for (j=0;j<(cdata->cdim);j++) ip[j] = (cdata->site_num)((void *)cdata, v[j]);
717  for (j=0;j<(cdata->cdim);j++) pp[j] = &(cdata->input_vert_array[ip[j]]);
718 
719  /* Don't add a point if it's already added */
720  for (j=0;j<(cdata->cdim);j++){
721  std::map<long *, int>::iterator pli;
722  pli = cdata->point_lookup->find((long *)pp[j]);
723  if (pli == cdata->point_lookup->end()) {
724  f = bu_ptbl_ins(cdata->output_pnts, (long *)pp[j]);
725  cdata->point_lookup->insert(std::make_pair((long *)pp[j], f));
726  VMOVE(cdata->vert_array[f],*pp[j]);
727  (*cdata->vert_cnt)++;
728  }
729  }
730  return NULL;
731 }
732 
733 
734 HIDDEN void *
735 chull3d_collect_faces(struct chull3d_data *cdata, simplex *s, void *UNUSED(p)) {
736 
737  point v[MAXDIM];
738  int j;
739  long int ip[3];
740  /*point_t p1, p2, p3;*/
741  const point_t *pp[3];
742  int f[3];
743  vect_t a, b, normal, center_to_edge;
744 
745  if (!s) return NULL;
746 
747  /* If we're not 3d, it doesn't make sense to do this */
748  if (cdata->cdim != 3) return NULL;
749 
750  for (j=0;j<(cdata->cdim);j++) v[j] = s->neigh[j].vert;
751  for (j=0;j<(cdata->cdim);j++) ip[j] = (cdata->site_num)((void *)cdata, v[j]);
752  for (j=0;j<(cdata->cdim);j++) pp[j] = &(cdata->input_vert_array[ip[j]]);
753 
754  for (j=0;j<(cdata->cdim);j++){
755  std::map<long *, int>::iterator pli;
756  pli = cdata->point_lookup->find((long *)pp[j]);
757  f[j] = pli->second;
758  }
759 
760  /* Use cdata->center_pnt and the center pnt of the triangle to construct
761  * a vector, and find the normal of the proposed face. If the face normal
762  * does not point in the same direction (dot product is positive) then
763  * the face needs to be reversed. Do this so the eventual bot primitive
764  * won't need a bot sync. Since by definition the chull is convex, the
765  * center point should be "inside" relative to all faces and this test
766  * should be reliable. */
767  VSUB2(a, cdata->vert_array[f[1]], cdata->vert_array[f[0]]);
768  VSUB2(b, cdata->vert_array[f[2]], cdata->vert_array[f[0]]);
769  VCROSS(normal, a, b);
770  VUNITIZE(normal);
771 
772  VSUB2(center_to_edge, cdata->vert_array[f[0]], cdata->center_pnt);
773  VUNITIZE(center_to_edge);
774 
775  cdata->faces[(*cdata->face_cnt)*3] = f[0];
776  if(VDOT(normal, center_to_edge) < 0) {
777  cdata->faces[(*cdata->face_cnt)*3+1] = f[2];
778  cdata->faces[(*cdata->face_cnt)*3+2] = f[1];
779  } else {
780  cdata->faces[(*cdata->face_cnt)*3+1] = f[1];
781  cdata->faces[(*cdata->face_cnt)*3+2] = f[2];
782  }
783  (*cdata->face_cnt)++;
784  return NULL;
785 }
786 
787 
788 #ifndef HAVE_LOGB
789 double logb(double x) {
790  if (x<=0) return -1e305;
791  return log(x)/log(2);
792 }
793 #endif
794 
796 chull3d_maxdist(int dim, point p1, point p2)
797 {
798  Coord x = 0;
799  Coord y = 0;
800  Coord d = 0;
801  int i = dim;
802  while (i--) {
803  x = *p1++;
804  y = *p2++;
805  d += (x<y) ? y-x : x-y ;
806  }
807  return d;
808 }
809 
810 /* the neighbor entry of a containing b */
813 {
814  int i;
815  neighbor *x;
816  for (i=0, x = a->neigh; (x->simp != b) && (i<(cdata->cdim)) ; i++, x++);
817  if (i<(cdata->cdim))
818  return x;
819  else {
820  return 0;
821  }
822 }
823 
826 {
827  int i;
828  neighbor *x;
829  for (i=0, x = a->neigh; (x->vert != b) && (i<(cdata->cdim)) ; i++, x++);
830  if (i<(cdata->cdim))
831  return x;
832  else {
833  return 0;
834  }
835 }
836 
837 
838 /* make neighbor connections between newly created simplices incident to p */
839 HIDDEN void
841 {
842  site xf,xb,xfi;
843  simplex *sb, *sf, *seen;
844  int i;
845  neighbor *sn;
846 
847  if (!s) return;
848  //assert(!s->peak.vert && s->peak.simp->peak.vert==cdata->p && !chull3d_op_vert(cdata,s,cdata->p)->simp->peak.vert);
849  if (s->visit==cdata->pnum) return;
850  s->visit = cdata->pnum;
851  seen = s->peak.simp;
852  xfi = chull3d_op_simp(cdata,seen,s)->vert;
853  for (i=0, sn = s->neigh; i<(cdata->cdim); i++,sn++) {
854  xb = sn->vert;
855  if (cdata->p == xb) continue;
856  sb = seen;
857  sf = sn->simp;
858  xf = xfi;
859  if (!sf->peak.vert) { /* are we done already? */
860  sf = chull3d_op_vert(cdata,seen,xb)->simp;
861  if (sf->peak.vert) continue;
862  } else do {
863  xb = xf;
864  xf = chull3d_op_simp(cdata,sf,sb)->vert;
865  sb = sf;
866  sf = chull3d_op_vert(cdata,sb,xb)->simp;
867  } while (sf->peak.vert);
868 
869  sn->simp = sf;
870  chull3d_op_vert(cdata,sf,xf)->simp = s;
871 
872  chull3d_connect(cdata, sf);
873  }
874 }
875 
876 
877 /* visit simplices s with sees(p,s), and make a facet for every neighbor
878  * of s not seen by p */
879 HIDDEN simplex *
881 {
882  simplex *n;
883  neighbor *bn;
884  int i;
885 
886 
887  if (!seen) return NULL;
888  seen->peak.vert = cdata->p;
889 
890  for (i=0,bn = seen->neigh; i<(cdata->cdim); i++,bn++) {
891  neighbor *ns;
892  n = bn->simp;
893  if (cdata->pnum != n->visit) {
894  n->visit = cdata->pnum;
895  if (chull3d_sees(cdata, cdata->p,n)) chull3d_make_facets(cdata, n);
896  }
897  if (n->peak.vert) continue;
898  CHULL3D_copy_simp(cdata, cdata->ns, seen);
899  cdata->ns->visit = 0;
900  cdata->ns->peak.vert = 0;
901  cdata->ns->normal = 0;
902  cdata->ns->peak.simp = seen;
903  CHULL3D_NULLIFY(cdata, basis_s,cdata->ns->neigh[i].basis);
904  cdata->ns->neigh[i].vert = cdata->p;
905  bn->simp = cdata->ns;
906  ns = chull3d_op_simp(cdata,n,seen);
907  ns->simp = cdata->ns;
908  }
909  return cdata->ns;
910 }
911 
912 
913 /* p lies outside flat containing previous sites;
914  * make p a vertex of every current simplex, and create some new simplices */
915 HIDDEN simplex *
917 {
918  int i;
919  int ocdim=(cdata->cdim)-1;
920  simplex *ns = NULL;
921  neighbor *nsn = NULL;
922 
923 
924  if (s->visit == cdata->pnum) return s->peak.vert ? s->neigh[ocdim].simp : s;
925  s->visit = cdata->pnum;
926  s->neigh[ocdim].vert = cdata->p;
927  CHULL3D_NULLIFY(cdata, basis_s,s->normal);
928  CHULL3D_NULLIFY(cdata, basis_s,s->neigh[0].basis);
929  if (!s->peak.vert) {
930  s->neigh[ocdim].simp = chull3d_extend_simplices(cdata, s->peak.simp);
931  return s;
932  } else {
933  CHULL3D_copy_simp(cdata, ns, s);
934  s->neigh[ocdim].simp = ns;
935  ns->peak.vert = NULL;
936  ns->peak.simp = s;
937  ns->neigh[ocdim] = s->peak;
938  if (s->peak.basis) s->peak.basis->ref_count++;
939  for (i=0,nsn=ns->neigh;i<(cdata->cdim);i++,nsn++)
940  nsn->simp = chull3d_extend_simplices(cdata, nsn->simp);
941  }
942  return ns;
943 }
944 
945 
946 /* return a simplex s that corresponds to a facet of the
947  * current hull, and sees(p, s) */
948 HIDDEN simplex *
949 chull3d_search(struct chull3d_data *cdata, simplex *root)
950 {
951  simplex *s;
952  neighbor *sn;
953  int i;
954  long tms = 0;
955 
956 
957  CHULL3D_push(root->peak.simp);
958  root->visit = cdata->pnum;
959  if (!chull3d_sees(cdata, cdata->p,root))
960  for (i=0,sn=root->neigh;i<(cdata->cdim);i++,sn++) CHULL3D_push(sn->simp);
961  while (tms) {
962  if (tms>cdata->ss)
963  cdata->st=(simplex**)realloc(cdata->st, ((cdata->ss+=cdata->ss)+MAXDIM+1)*sizeof(simplex*));
964  CHULL3D_pop(s);
965  if (s->visit == cdata->pnum) continue;
966  s->visit = cdata->pnum;
967  if (!chull3d_sees(cdata, cdata->p,s)) continue;
968  if (!s->peak.vert) return s;
969  for (i=0, sn=s->neigh; i<(cdata->cdim); i++,sn++) CHULL3D_push(sn->simp);
970  }
971  return NULL;
972 }
973 
974 
977 {
978  point pnext;
979  pnext = (*(cdata->get_site))((void *)cdata);
980  if (!pnext) return NULL;
981  cdata->pnum = cdata->site_num((void *)cdata, pnext)+2;
982  return pnext;
983 }
984 
985 
986 
987 HIDDEN void
989 {
990  while ((cdata->cdim) < cdata->rdim) {
991  cdata->p = chull3d_get_another_site(cdata);
992  if (!cdata->p) return;
993  if (chull3d_out_of_flat(cdata, root,cdata->p))
994  chull3d_extend_simplices(cdata, root);
995  else
996  chull3d_connect(cdata, chull3d_make_facets(cdata, chull3d_search(cdata, root)));
997  }
998  while ((cdata->p = chull3d_get_another_site(cdata)))
999  chull3d_connect(cdata, chull3d_make_facets(cdata, chull3d_search(cdata, root)));
1000 }
1001 
1002 
1003 /*
1004  get_s returns next site each call;
1005  hull construction stops when NULL returned;
1006  site_numm returns number of site when given site;
1007  dim dimension of point set;
1008  */
1009 HIDDEN simplex *
1010 chull3d_build_convex_hull(struct chull3d_data *cdata, gsitef *get_s, site_n *site_numm, short dim)
1011 {
1012  double febits;
1013  simplex *s = NULL;
1014  simplex *root = NULL;
1015 
1016  if (ZERO(cdata->Huge)) cdata->Huge = DBL_MAX*DBL_MAX;
1017 
1018  (cdata->cdim) = 0;
1019  cdata->get_site = get_s;
1020  cdata->site_num = site_numm;
1021  cdata->pdim = dim;
1022 
1023  febits = floor(DBL_MANT_DIG*log(FLT_RADIX)/log(2));
1024  cdata->exact_bits = (int)febits;
1025  cdata->b_err_min = (float)(DBL_EPSILON*MAXDIM*(1<<MAXDIM)*MAXDIM*3.01);
1026  cdata->b_err_min_sq = cdata->b_err_min * cdata->b_err_min;
1027 
1028  //assert(cdata->get_site!=NULL); assert(cdata->site_num!=NULL);
1029 
1030  cdata->rdim = cdata->pdim;
1031 
1032  cdata->point_size = cdata->site_size = sizeof(Coord)*cdata->pdim;
1033  cdata->basis_vec_size = sizeof(Coord)*cdata->rdim;
1034  cdata->basis_s_size = sizeof(basis_s)+ (2*cdata->rdim-1)*sizeof(Coord);
1035  cdata->simplex_size = sizeof(simplex) + (cdata->rdim-1)*sizeof(neighbor);
1036  cdata->Tree_size = sizeof(Tree);
1037 
1038  root = NULL;
1039  if (!(cdata->p = (*(cdata->get_site))((void *)cdata))) return 0;
1040 
1041  CHULL3D_NEWL(cdata, cdata->simplex_list, simplex,root);
1042 
1043  cdata->ch_root = root;
1044 
1045  CHULL3D_copy_simp(cdata, s, root);
1046  root->peak.vert = cdata->p;
1047  root->peak.simp = s;
1048  s->peak.simp = root;
1049 
1050  chull3d_buildhull(cdata, root);
1051  return root;
1052 }
1053 
1054 
1055 /* hullmain.c */
1056 
1057 HIDDEN long
1059 {
1060  long i,j;
1061  struct chull3d_data *cdata = (struct chull3d_data *)data;
1062 
1063  if (!p) return -2;
1064  for (i=0; i<cdata->num_blocks; i++)
1065  if ((j=p-cdata->site_blocks[i])>=0 && j < BLOCKSIZE*cdata->pdim)
1066  return j/cdata->pdim+BLOCKSIZE*i;
1067  return -3;
1068 }
1069 
1070 HIDDEN site
1071 chull3d_new_site(struct chull3d_data *cdata, site p, long j)
1072 {
1073  //assert(cdata->num_blocks+1<cdata->input_vert_cnt);
1074  if (0==(j%BLOCKSIZE)) {
1075  //assert(cdata->num_blocks < cdata->input_vert_cnt);
1076  return(cdata->site_blocks[cdata->num_blocks++]=(site)malloc(BLOCKSIZE*cdata->site_size));
1077  } else
1078  return p+cdata->pdim;
1079 }
1080 
1081 HIDDEN site
1082 chull3d_read_next_site(struct chull3d_data *cdata, long j)
1083 {
1084  int i=0;
1085 
1086  cdata->p = chull3d_new_site(cdata, cdata->p,j);
1087 
1088  if (cdata->next_vert == cdata->input_vert_cnt) return 0;
1089 
1090  cdata->p[0] = cdata->input_vert_array[cdata->next_vert][0];
1091  cdata->p[1] = cdata->input_vert_array[cdata->next_vert][1];
1092  cdata->p[2] = cdata->input_vert_array[cdata->next_vert][2];
1093  (cdata->next_vert)++;
1094  for(i = 0; i < 3; i++) {
1095  (cdata->p)[i] = floor(cdata->mult_up*(cdata->p)[i]+0.5);
1096  cdata->mins[i] = (cdata->mins[i]<(cdata->p)[i]) ? cdata->mins[i] : (cdata->p)[i];
1097  cdata->maxs[i] = (cdata->maxs[i]>(cdata->p)[i]) ? cdata->maxs[i] : (cdata->p)[i];
1098  }
1099 
1100  return cdata->p;
1101 }
1102 
1103 HIDDEN void
1105 {
1106  long i, j, t;
1107  for (i=0;i<cdata->input_vert_cnt;i++) cdata->shufmat[i] = i;
1108  for (i=0;i<cdata->input_vert_cnt;i++) {
1109  double random_num = bn_randmt();
1110  t = cdata->shufmat[i];
1111  j = i + (cdata->input_vert_cnt - i) * (long)random_num;
1112  cdata->shufmat[i] = cdata->shufmat[j];
1113  cdata->shufmat[j] = t;
1114  }
1115 }
1116 
1117 HIDDEN site
1119  struct chull3d_data *cdata = (struct chull3d_data *)data;
1120  return chull3d_read_next_site(cdata, cdata->shufmat[(cdata->s_num)++]);
1121 }
1122 
1123 HIDDEN void
1125 {
1126  int i = 0;
1127  data->simplex_list = 0;
1128  data->basis_s_list = 0;
1129  data->Tree_list = 0;
1130  data->site_blocks = (point *)bu_calloc(vert_cnt * 3, sizeof(point), "site_blocks");
1131  BU_GET(data->tt_basis, basis_s);
1132  data->tt_basis->next = 0;
1133  data->tt_basis->ref_count = 1;
1134  data->tt_basis->lscale = -1;
1135  data->tt_basis->sqa = 0;
1136  data->tt_basis->sqb = 0;
1137  data->tt_basis->vecs[0] = 0;
1138  /* point at infinity for Delaunay triang; value not used */
1139  data->hull_infinity = (Coord *)bu_calloc(10, sizeof(Coord), "hull_infinity");
1140  data->hull_infinity[0] = 57.2;
1141  data->B = (int *)bu_calloc(vert_cnt, sizeof(int), "A");
1142  data->mins = (Coord *)bu_calloc(MAXDIM, sizeof(Coord), "mins");
1143  for(i = 0; i < MAXDIM; i++) data->mins[i] = DBL_MAX;
1144  data->maxs = (Coord *)bu_calloc(MAXDIM, sizeof(Coord), "maxs");
1145  for(i = 0; i < MAXDIM; i++) data->maxs[i] = -DBL_MAX;
1146  data->mult_up = 1000.0; /* we'll need to multiply based on a tolerance parameter at some point... */
1147  data->mi = (short *)bu_calloc(vert_cnt, sizeof(short), "mi");
1148  data->mo = (short *)bu_calloc(vert_cnt, sizeof(short), "mo");
1149  data->pdim = 3;
1150  data->shufmat = (long*)bu_calloc(vert_cnt+1, sizeof(long), "shufmat");
1151 
1152  /* These were static variables in functions */
1153  data->simplex_block_table = (simplex **)bu_calloc(vert_cnt, sizeof(simplex *), "simplex_block_table");
1154  data->num_simplex_blocks = 0;
1155  data->basis_s_block_table = (basis_s **)bu_calloc(vert_cnt, sizeof(basis_s *), "basis_s_block_table");
1156  data->num_basis_s_blocks = 0;
1157  data->Tree_block_table = (Tree **)bu_calloc(vert_cnt, sizeof(Tree *), "Tree_block_table");
1158  data->num_Tree_blocks = 0;
1159  data->p_neigh.vert = 0;
1160  data->p_neigh.simp = NULL;
1161  data->p_neigh.basis = NULL;
1162  data->b = NULL;
1163  data->vnum = -1;
1164  data->ss = vert_cnt + MAXDIM;
1165  data->s_num = 0;
1166  data->fg_vn = 0;
1167  data->st=(simplex**)malloc((data->ss+MAXDIM+1)*sizeof(simplex*));
1168  data->ns = NULL;
1169  data->free_point_lookup = 0;
1170  data->free_output_pnts = 0;
1171  if (!data->point_lookup) {
1172  data->point_lookup = new std::map<long *, int>;
1173  data->free_point_lookup = 1;
1174  }
1175  if (!data->output_pnts) {
1176  BU_GET(data->output_pnts, struct bu_ptbl);
1177  bu_ptbl_init(data->output_pnts, vert_cnt, "output pnts container");
1178  data->free_output_pnts = 1;
1179  }
1180  data->input_vert_cnt = vert_cnt;
1181 
1182  VSET(data->center_pnt,0,0,0);
1183 
1184  /* Output containers */
1185  if (!data->vert_array) data->vert_array = (point_t *)bu_calloc(vert_cnt, sizeof(point_t), "vertex array");
1186  data->next_vert = 0;
1187  if (!data->faces) data->faces = (int *)bu_calloc(3*vert_cnt, sizeof(int), "face array");
1188 }
1189 
1190  HIDDEN void
1192 {
1193  bu_free(data->site_blocks, "site_blocks");
1194  bu_free(data->hull_infinity, "hull_infinity");
1195  bu_free(data->B, "B");
1196  bu_free(data->mins, "mins");
1197  bu_free(data->maxs, "maxs");
1198  bu_free(data->mi, "mi");
1199  bu_free(data->mo, "mo");
1200  bu_free(data->simplex_block_table, "simplex_block_table");
1201  bu_free(data->basis_s_block_table, "basis_s_block_table");
1202  bu_free(data->Tree_block_table, "Tree_block_table");
1203  BU_PUT(data->tt_basis, basis_s);
1204  if (data->free_output_pnts) {
1205  bu_ptbl_free(data->output_pnts);
1206  BU_PUT(data->output_pnts, struct bu_ptbl);
1207  }
1208  if (data->free_point_lookup) delete data->point_lookup;
1209 }
1210 
1211 HIDDEN
1212 int
1213 chull3d_intermediate_set(point_t **vertices, int *num_vertices, const point_t *input_points_3d, int num_input_pnts)
1214 {
1215  int nv = 0;
1216  int nf = 0;
1217  int curr_stp = 0;
1218  int pnt_stp = 0;
1219  int last_stp = 0;
1220  int pnts_per_stp = CHULL3D_SUBSET_SIZE;
1221  const point_t *curr_inputs;
1222  struct bu_ptbl *opnts;
1223  std::map<long *, int> *pl = new std::map<long *, int>;
1224  point_t *intermediate_vert_array = (point_t *)bu_calloc(num_input_pnts, sizeof(point_t), "vertex array");
1225  simplex *root = NULL;
1226  BU_GET(opnts, struct bu_ptbl);
1227  bu_ptbl_init(opnts, num_input_pnts, "output pnts container");
1228  pnt_stp = (num_input_pnts < pnts_per_stp) ? num_input_pnts : pnts_per_stp;
1229  last_stp = num_input_pnts / pnts_per_stp;
1230  //bu_log("input point cnt: %d\n", num_input_pnts);
1231  while (curr_stp <= last_stp) {
1232  struct chull3d_data *cdata;
1233  //bu_log("step %d of %d\n", curr_stp, last_stp);
1234  curr_inputs = input_points_3d + curr_stp * pnt_stp;
1235  if (curr_stp == last_stp) pnt_stp = num_input_pnts - curr_stp * pnt_stp;
1236  BU_GET(cdata, struct chull3d_data);
1237  cdata->output_pnts = opnts;
1238  cdata->point_lookup = pl;
1239  cdata->vert_array = intermediate_vert_array;
1240  chull3d_data_init(cdata, pnt_stp);
1241  cdata->input_vert_array = curr_inputs;
1242  cdata->vert_cnt = &nv;
1243  cdata->face_cnt = &nf;
1244  chull3d_make_shuffle(cdata);
1246  if (!root) return -1;
1249  chull3d_data_free(cdata);
1250  BU_PUT(cdata, struct chull3d_data);
1251  curr_stp++;
1252  //bu_log("points accumulated: %d\n", nv);
1253  }
1254  //bu_log("output pnt cnt: %d\n", nv);
1255  bu_ptbl_free(opnts);
1256  BU_PUT(opnts, struct bu_ptbl);
1257  delete pl;
1258 
1259  (*vertices) = intermediate_vert_array;
1260  (*num_vertices) = nv;
1261  return 0;
1262 }
1263 
1264 int
1265 bn_3d_chull(int **faces, int *num_faces, point_t **vertices, int *num_vertices,
1266  const point_t *input_points_3d, int num_input_pnts)
1267 {
1268  int i;
1269  struct chull3d_data *cdata;
1270  simplex *root = NULL;
1271  int output_dim = 0;
1272  point_t* hull_2d;
1273 
1274  const point_t *iv_1;
1275  point_t *iv_2, *iva;
1276  int nv_1, nv_2, nva;
1277  iv_1 = input_points_3d;
1278  nv_1 = num_input_pnts;
1279 
1280  (void)chull3d_intermediate_set(&iv_2, &nv_2, iv_1, nv_1);
1281  iva = iv_2;
1282  nva = nv_2;
1283  while (nva > CHULL3D_COUNT_THRESHOLD && abs(nv_1 - nv_2) > CHULL3D_DELTA_THRESHOLD) {
1284  nv_1 = nv_2;
1285  iv_1 = (const point_t *)iv_2;
1286  (void)chull3d_intermediate_set(&iv_2, &nv_2, iv_1, nv_1);
1287  bu_free((point_t *)iv_1, "free old set");
1288  iva = iv_2;
1289  nva = nv_2;
1290  }
1291 
1292  BU_GET(cdata, struct chull3d_data);
1293  chull3d_data_init(cdata, nva);
1294  cdata->input_vert_array = (const point_t *)iva;
1295  cdata->vert_cnt = num_vertices;
1296  cdata->face_cnt = num_faces;
1297  (*cdata->vert_cnt) = 0;
1298  (*cdata->face_cnt) = 0;
1299  chull3d_make_shuffle(cdata);
1301  if (!root) return -1;
1303  switch(cdata->cdim) {
1304  case 0:
1305  bu_log("point?\n");
1306  break;
1307  case 1:
1308  bu_log("line?\n");
1309  break;
1310  case 2:
1311  /* We already have the hull points, but we need to assemble them into a CCW hull */
1312  bn_3d_coplanar_chull(&hull_2d, (const point_t *)cdata->vert_array, (*cdata->vert_cnt));
1313  (*vertices) = hull_2d;
1314  bu_free(cdata->vert_array, "using 2d vertex list from coplanar_chull");
1315  (*faces) = NULL; /* Should we tessellate the hull and make faces? */
1316  bu_free(cdata->faces, "free face list - not used in 2d");
1317  break;
1318  case 3:
1319  /* 3d hull - calculate the center point so we can correctly add triangle faces
1320  * with normals outward from the convex hull */
1321  for(i = 0; i < (int)BU_PTBL_LEN(cdata->output_pnts); i++) {
1322  point_t p1;
1323  VMOVE(p1,cdata->vert_array[i]);
1324  cdata->center_pnt[0] += p1[0];
1325  cdata->center_pnt[1] += p1[1];
1326  cdata->center_pnt[2] += p1[2];
1327  }
1328  cdata->center_pnt[0] = cdata->center_pnt[0]/(*cdata->vert_cnt);
1329  cdata->center_pnt[1] = cdata->center_pnt[1]/(*cdata->vert_cnt);
1330  cdata->center_pnt[2] = cdata->center_pnt[2]/(*cdata->vert_cnt);
1331  /* Build the faces list and assign the results */
1333  (*vertices) = cdata->vert_array;
1334  (*faces) = cdata->faces;
1335  break;
1336  default:
1337  break;
1338  }
1339  output_dim = cdata->cdim;
1341  chull3d_data_free(cdata);
1342  bu_free(iva, "intermediate_vert_array");
1343  BU_PUT(cdata, struct chull3d_data);
1344 
1345  return output_dim;
1346 }
1347 
1348 /*
1349  * Local Variables:
1350  * mode: C
1351  * tab-width: 8
1352  * indent-tabs-mode: t
1353  * c-file-style: "stroustrup"
1354  * End:
1355  * ex: shiftwidth=4 tabstop=8
1356  */
1357 
neighbor neigh[1]
Definition: chull3d.cpp:77
HIDDEN simplex * chull3d_new_block_simplex(struct chull3d_data *cdata, int make_blocks)
Definition: chull3d.cpp:220
HIDDEN double chull3d_sc(struct chull3d_data *cdata, basis_s *v, simplex *s, int k, int j)
Definition: chull3d.cpp:356
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define CHULL3D_SUBSET_SIZE
Definition: chull3d.cpp:49
size_t Tree_size
Definition: chull3d.cpp:108
size_t simplex_size
Definition: chull3d.cpp:104
Tree * next
Definition: chull3d.cpp:87
int num_simplex_blocks
Definition: chull3d.cpp:141
struct tree_node Tree
Definition: chull3d.cpp:81
struct simplex * next
Definition: chull3d.cpp:72
int point_size
Definition: chull3d.cpp:134
#define MAXDIM
Definition: chull3d.cpp:36
HIDDEN void * chull3d_check_simplex(struct chull3d_data *cdata, simplex *s, void *dum)
Definition: chull3d.cpp:669
HIDDEN simplex * chull3d_search(struct chull3d_data *cdata, simplex *root)
Definition: chull3d.cpp:949
int basis_vec_size
Definition: chull3d.cpp:115
double logb(double x)
Definition: chull3d.cpp:789
#define CHULL3D_DELTA_THRESHOLD
Definition: chull3d.cpp:48
#define CHULL3D_trans(z, p, q)
Definition: chull3d.cpp:179
Coord * maxs
Definition: chull3d.cpp:128
double Coord
Definition: chull3d.cpp:51
if lu s
Definition: nmg_mod.c:3860
HIDDEN Coord chull3d_Norm2(struct chull3d_data *cdata, point x)
Definition: chull3d.cpp:324
#define CHULL3D_NEWLRC(cdata, list, X, p)
Definition: chull3d.cpp:192
void bu_ptbl_init(struct bu_ptbl *b, size_t len, const char *str)
Definition: ptbl.c:32
#define VSET(a, b, c, d)
Definition: color.c:53
gsitef * get_site
Definition: chull3d.cpp:123
HIDDEN Coord chull3d_Vec_dot(struct chull3d_data *cdata, point x, point y)
Definition: chull3d.cpp:306
int * faces
Definition: chull3d.cpp:172
fg * next
Definition: chull3d.cpp:93
const point_t * input_vert_array
Definition: chull3d.cpp:162
#define CHULL3D_NULLIFY(cdata, X, v)
Definition: chull3d.cpp:200
int free_output_pnts
Definition: chull3d.cpp:165
simplex * simplex_list
Definition: chull3d.cpp:105
Coord sqa
Definition: chull3d.cpp:61
Coord vecs[1]
Definition: chull3d.cpp:62
double ldetbound
Definition: chull3d.cpp:148
HIDDEN basis_s * chull3d_new_block_basis_s(struct chull3d_data *cdata, int make_blocks)
Definition: chull3d.cpp:248
site vert
Definition: chull3d.cpp:66
HIDDEN int chull3d_visit_fg_i_far(struct chull3d_data *cdata, void(*v_fg)(struct chull3d_data *, Tree *, int), Tree *t, int depth, int vn)
Definition: chull3d.cpp:591
HIDDEN site chull3d_get_next_site(void *data)
Definition: chull3d.cpp:1118
HIDDEN void * chull3d_visit_triang_gen(struct chull3d_data *cdata, simplex *s, chull3d_visit_func *visit, chull3d_test_func *test)
Definition: chull3d.cpp:626
HIDDEN int chull3d_sees(struct chull3d_data *, site, simplex *)
Definition: chull3d.cpp:528
fg * faces_gr_t
Definition: chull3d.cpp:125
Header file for the BRL-CAD common definitions.
fg * fgs
Definition: chull3d.cpp:86
HIDDEN point chull3d_get_another_site(struct chull3d_data *cdata)
Definition: chull3d.cpp:976
HIDDEN int chull3d_truet(struct chull3d_data *cdata, simplex *s, int i, void *dum)
Definition: chull3d.cpp:652
point site
Definition: chull3d.cpp:53
Coord site_struct
Definition: chull3d.cpp:55
HIDDEN void chull3d_Ax_plus_y(struct chull3d_data *cdata, Coord a, point x, point y)
Definition: chull3d.cpp:333
int bu_ptbl_ins(struct bu_ptbl *b, long *p)
Tree * right
Definition: chull3d.cpp:83
HIDDEN Coord chull3d_maxdist(int dim, point p1, point p2)
Definition: chull3d.cpp:796
HIDDEN void chull3d_free_simplex_storage(struct chull3d_data *cdata)
Definition: chull3d.cpp:242
#define HIDDEN
Definition: common.h:86
point_t center_pnt
Definition: chull3d.cpp:167
int num_Tree_blocks
Definition: chull3d.cpp:145
HIDDEN void chull3d_get_basis_sede(struct chull3d_data *cdata, simplex *s)
Definition: chull3d.cpp:453
point_t * vert_array
Definition: chull3d.cpp:170
HIDDEN simplex * chull3d_make_facets(struct chull3d_data *cdata, simplex *seen)
Definition: chull3d.cpp:880
basis_s * normal
Definition: chull3d.cpp:75
HIDDEN int chull3d_out_of_flat(struct chull3d_data *cdata, simplex *root, point p)
Definition: chull3d.cpp:470
float b_err_min
Definition: chull3d.cpp:117
HIDDEN void chull3d_data_free(struct chull3d_data *data)
Definition: chull3d.cpp:1191
void * chull3d_visit_func(struct chull3d_data *, simplex *, void *)
Definition: chull3d.cpp:619
Definition: ptbl.h:62
HIDDEN void * chull3d_visit_triang(struct chull3d_data *cdata, simplex *root, chull3d_visit_func *visit)
Definition: chull3d.cpp:654
Coord * hull_infinity
Definition: chull3d.cpp:121
HIDDEN neighbor * chull3d_op_simp(struct chull3d_data *cdata, simplex *a, simplex *b)
Definition: chull3d.cpp:812
#define CHULL3D_push(x)
Definition: chull3d.cpp:182
HIDDEN void chull3d_Vec_scale(struct chull3d_data *cdata, int n, Coord a, Coord *x)
Definition: chull3d.cpp:342
COMPLEX data[64]
Definition: fftest.c:34
short * mi
Definition: chull3d.cpp:135
HIDDEN void chull3d_make_shuffle(struct chull3d_data *cdata)
Definition: chull3d.cpp:1104
void * memset(void *s, int c, size_t n)
struct basis_s * next
Definition: chull3d.cpp:58
HIDDEN int chull3d_hullt(struct chull3d_data *cdata, simplex *s, int i, void *dummy)
Definition: chull3d.cpp:660
HIDDEN void chull3d_visit_fg_far(struct chull3d_data *cdata, fg *faces_gr, void(*v_fg)(struct chull3d_data *, Tree *, int))
Definition: chull3d.cpp:611
HIDDEN site chull3d_read_next_site(struct chull3d_data *cdata, long j)
Definition: chull3d.cpp:1082
HIDDEN void chull3d_visit_fg(struct chull3d_data *cdata, fg *faces_gr, void(*v_fg)(struct chull3d_data *, Tree *, int, int))
Definition: chull3d.cpp:585
#define CHULL3D_VA(x)
Definition: chull3d.cpp:177
long visit
Definition: chull3d.cpp:73
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
long site_n(void *, site)
Definition: chull3d.cpp:100
int bn_3d_chull(int **faces, int *num_faces, point_t **vertices, int *num_vertices, const point_t *input_points_3d, int num_input_pnts)
Find 3D point convex hull for unordered point sets.
Definition: chull3d.cpp:1265
int * face_cnt
Definition: chull3d.cpp:173
int num_basis_s_blocks
Definition: chull3d.cpp:143
HIDDEN Tree * chull3d_new_block_Tree(struct chull3d_data *cdata, int make_blocks)
Definition: chull3d.cpp:277
Tree * facets
Definition: chull3d.cpp:91
struct basis_s basis_s
HIDDEN void chull3d_get_normal(struct chull3d_data *cdata, simplex *s)
Definition: chull3d.cpp:525
#define CHULL3D_copy_simp(cdata, cnew, s)
Definition: chull3d.cpp:209
#define BU_GET(_ptr, _type)
Definition: malloc.h:201
#define CHULL3D_COUNT_THRESHOLD
Definition: chull3d.cpp:47
HIDDEN void chull3d_data_init(struct chull3d_data *data, int vert_cnt)
Definition: chull3d.cpp:1124
HIDDEN void * chull3d_facet_test(struct chull3d_data *cdata, simplex *s, void *dummy)
Definition: chull3d.cpp:661
HIDDEN simplex * chull3d_extend_simplices(struct chull3d_data *cdata, simplex *s)
Definition: chull3d.cpp:916
int size
Definition: chull3d.cpp:85
std::map< long *, int > * point_lookup
Definition: chull3d.cpp:158
Coord * point
Definition: chull3d.cpp:52
struct fg_node fg_node
float b_err_min_sq
Definition: chull3d.cpp:118
double Sb
Definition: chull3d.cpp:149
HIDDEN void * chull3d_collect_hull_pnts(struct chull3d_data *cdata, simplex *s, void *p)
Definition: chull3d.cpp:704
struct bu_list l
Definition: ptbl.h:63
Tree * Tree_list
Definition: chull3d.cpp:109
HIDDEN void * chull3d_collect_faces(struct chull3d_data *cdata, simplex *s, void *p)
Definition: chull3d.cpp:735
#define UNUSED(parameter)
Definition: common.h:239
#define BU_PUT(_ptr, _type)
Definition: malloc.h:215
HIDDEN void chull3d_free_Tree_storage(struct chull3d_data *cdata)
Definition: chull3d.cpp:298
short * mo
Definition: chull3d.cpp:136
site gsitef(void *)
Definition: chull3d.cpp:99
HIDDEN Coord chull3d_Vec_dot_pdim(struct chull3d_data *cdata, point x, point y)
Definition: chull3d.cpp:315
neighbor p_neigh
Definition: chull3d.cpp:150
HIDDEN simplex * chull3d_build_convex_hull(struct chull3d_data *cdata, gsitef *get_s, site_n *site_numm, short dim)
Definition: chull3d.cpp:1010
#define BU_PTBL_LEN(ptbl)
Definition: ptbl.h:107
void bu_ptbl_free(struct bu_ptbl *b)
Definition: ptbl.c:226
double Huge
Definition: chull3d.cpp:114
HIDDEN neighbor * chull3d_op_vert(struct chull3d_data *cdata, simplex *a, site b)
Definition: chull3d.cpp:825
#define ZERO(val)
Definition: units.c:38
int bn_3d_coplanar_chull(point_t **hull, const point_t *points_3d, int n)
Find 3D coplanar point convex hull for unordered co-planar point sets.
Definition: chull.c:151
basis_s * tt_basis
Definition: chull3d.cpp:119
double mult_up
Definition: chull3d.cpp:126
int num_blocks
Definition: chull3d.cpp:112
long * shufmat
Definition: chull3d.cpp:137
int chull3d_test_func(struct chull3d_data *, simplex *, int, void *)
Definition: chull3d.cpp:618
HIDDEN long chull3d_site_numm(void *data, site p)
Definition: chull3d.cpp:1058
neighbor peak
Definition: chull3d.cpp:76
simplex ** simplex_block_table
Definition: chull3d.cpp:140
simplex * ns
Definition: chull3d.cpp:154
#define CHULL3D_VB(x)
Definition: chull3d.cpp:178
HIDDEN int chull3d_reduce(struct chull3d_data *cdata, basis_s **v, point p, simplex *s, int k)
Definition: chull3d.cpp:440
HIDDEN void chull3d_free_hull_storage(struct chull3d_data *cdata)
Definition: chull3d.cpp:558
Coord * normalp
Definition: chull3d.cpp:54
#define CHULL3D_NEWL(cdata, list, X, p)
Definition: chull3d.cpp:185
HIDDEN void chull3d_connect(struct chull3d_data *cdata, simplex *s)
Definition: chull3d.cpp:840
short mark
Definition: chull3d.cpp:94
#define CHULL3D_DELIFT
Definition: chull3d.cpp:180
double vol
Definition: chull3d.cpp:92
simplex ** st
Definition: chull3d.cpp:155
Tree ** Tree_block_table
Definition: chull3d.cpp:144
double dist
Definition: chull3d.cpp:92
#define BLOCKSIZE
Definition: chull3d.cpp:35
HIDDEN int chull3d_reduce_inner(struct chull3d_data *cdata, basis_s *v, simplex *s, int k)
Definition: chull3d.cpp:399
site key
Definition: chull3d.cpp:84
int input_vert_cnt
Definition: chull3d.cpp:163
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
HIDDEN void chull3d_get_normal_sede(struct chull3d_data *cdata, simplex *s)
Definition: chull3d.cpp:488
short mark
Definition: chull3d.cpp:74
double max_scale
Definition: chull3d.cpp:147
HIDDEN int chull3d_intermediate_set(point_t **vertices, int *num_vertices, const point_t *input_points_3d, int num_input_pnts)
Definition: chull3d.cpp:1213
size_t basis_s_size
Definition: chull3d.cpp:106
basis_s * infinity_basis
Definition: chull3d.cpp:120
Coord sqb
Definition: chull3d.cpp:61
basis_s * basis_s_list
Definition: chull3d.cpp:107
simplex * ch_root
Definition: chull3d.cpp:113
int free_point_lookup
Definition: chull3d.cpp:159
point * site_blocks
Definition: chull3d.cpp:111
basis_s * b
Definition: chull3d.cpp:151
int ref_count
Definition: chull3d.cpp:59
HIDDEN void chull3d_visit_fg_i(struct chull3d_data *cdata, void(*v_fg)(struct chull3d_data *, Tree *, int, int), Tree *t, int depth, int vn, int boundary)
Definition: chull3d.cpp:566
HIDDEN site chull3d_new_site(struct chull3d_data *cdata, site p, long j)
Definition: chull3d.cpp:1071
Coord * mins
Definition: chull3d.cpp:127
Tree * left
Definition: chull3d.cpp:83
int exact_bits
Definition: chull3d.cpp:116
struct simplex simplex
int * vert_cnt
Definition: chull3d.cpp:171
HIDDEN void chull3d_free_basis_s_storage(struct chull3d_data *cdata)
Definition: chull3d.cpp:270
basis_s ** basis_s_block_table
Definition: chull3d.cpp:142
HIDDEN void chull3d_buildhull(struct chull3d_data *cdata, simplex *root)
Definition: chull3d.cpp:988
site_n * site_num
Definition: chull3d.cpp:124
struct bu_ptbl * output_pnts
Definition: chull3d.cpp:164
basis_s * basis
Definition: chull3d.cpp:68
struct neighbor neighbor
int lscale
Definition: chull3d.cpp:60
HIDDEN void * chull3d_visit_hull(struct chull3d_data *cdata, simplex *root, chull3d_visit_func *visit)
Definition: chull3d.cpp:663
int ref_count
Definition: chull3d.cpp:95
double bn_randmt(void)
Mersenne Twister random number generation as defined by MT19937.
Definition: randmt.c:149
struct simplex * simp
Definition: chull3d.cpp:67
#define CHULL3D_pop(x)
Definition: chull3d.cpp:183