BRL-CAD
arbn.c
Go to the documentation of this file.
1 /* A R B N . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1989-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 /** @addtogroup primitives */
21 /** @{ */
22 /** @file primitives/arbn/arbn.c
23  *
24  * Intersect a ray with an Arbitrary Regular Polyhedron with an
25  * arbitrary number of faces.
26  *
27  */
28 
29 #include "common.h"
30 
31 #include <stdlib.h>
32 #include <stdio.h>
33 #include <string.h>
34 #include <math.h>
35 #include <ctype.h>
36 #include "bnetwork.h"
37 
38 #include "tcl.h"
39 #include "bu/cv.h"
40 #include "vmath.h"
41 #include "nmg.h"
42 #include "db.h"
43 #include "rtgeom.h"
44 #include "raytrace.h"
45 
46 /**
47  * Calculate a bounding RPP for an ARBN
48  */
49 int
50 rt_arbn_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol)) {
51  size_t i, j, k;
52  struct rt_arbn_internal *aip;
54  aip = (struct rt_arbn_internal *)ip->idb_ptr;
55  RT_ARBN_CK_MAGIC(aip);
56 
57  VSETALL((*min), INFINITY);
58  VSETALL((*max), -INFINITY);
59 
60  /* Discover all vertices, use to calculate RPP */
61  for (i = 0; i < aip->neqn-2; i++) {
62  for (j = i+1; j < aip->neqn-1; j++) {
63  double dot;
64 
65  /* If normals are parallel, no intersection */
66  dot = VDOT(aip->eqn[i], aip->eqn[j]);
67  if (((dot) <= -SMALL_FASTF) ? (NEAR_EQUAL((dot), -1.0, RT_DOT_TOL)) : (NEAR_EQUAL((dot), 1.0, RT_DOT_TOL))) continue;
68 
69  /* Have an edge line, isect with higher numbered planes */
70  for (k = j + 1; k < aip->neqn; k++) {
71  size_t m;
72  size_t next_k;
73  point_t pt;
74 
75  next_k = 0;
76 
77  if (bn_mkpoint_3planes(pt, aip->eqn[i], aip->eqn[j], aip->eqn[k]) < 0) continue;
78 
79  /* See if point is outside arb */
80  for (m = 0; m < aip->neqn; m++) {
81  if (i == m || j == m || k == m)
82  continue;
83  if (VDOT(pt, aip->eqn[m])-aip->eqn[m][3] > RT_LEN_TOL) {
84  next_k = 1;
85  break;
86  }
87  }
88  if (next_k != 0) continue;
89 
90  VMINMAX((*min), (*max), pt);
91  }
92  }
93  }
94 
95  return 0;
96 }
97 
98 
99 /**
100  * Returns -
101  * 0 OK
102  * !0 failure
103  */
104 int
105 rt_arbn_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
106 {
107  struct rt_arbn_internal *aip;
108  vect_t work;
109  fastf_t f;
110  size_t i;
111  size_t j;
112  size_t k;
113  int *used = (int *)0; /* plane eqn use count */
114  const struct bn_tol *tol = &rtip->rti_tol;
115 
116  RT_CK_DB_INTERNAL(ip);
117  aip = (struct rt_arbn_internal *)ip->idb_ptr;
118  RT_ARBN_CK_MAGIC(aip);
119 
120  used = (int *)bu_malloc(aip->neqn*sizeof(int), "arbn used[]");
121 
122  /*
123  * ARBN must be convex. Test for concavity.
124  * Byproduct is an enumeration of all the vertices,
125  * which are used to make the bounding RPP. No need
126  * to call the bbox routine, as the work must be duplicated
127  * here to count faces.
128  */
129 
130  /* Zero face use counts
131  * and make sure normal vectors are unit vectors
132  */
133  for (i = 0; i < aip->neqn; i++) {
134  double normalLen = MAGNITUDE(aip->eqn[i]);
135  double scale;
136  if (ZERO(normalLen)) {
137  bu_log("arbn has zero length normal vector\n");
138  return 1;
139  }
140  scale = 1.0 / normalLen;
141  HSCALE(aip->eqn[i], aip->eqn[i], scale);
142  used[i] = 0;
143  }
144  for (i = 0; i < aip->neqn-2; i++) {
145  for (j=i+1; j<aip->neqn-1; j++) {
146  double dot;
147 
148  /* If normals are parallel, no intersection */
149  dot = VDOT(aip->eqn[i], aip->eqn[j]);
150  if (BN_VECT_ARE_PARALLEL(dot, tol)) continue;
151 
152  /* Have an edge line, isect with higher numbered planes */
153  for (k=j+1; k<aip->neqn; k++) {
154  size_t m;
155  size_t next_k;
156  point_t pt;
157 
158  next_k = 0;
159 
160  if (bn_mkpoint_3planes(pt, aip->eqn[i], aip->eqn[j], aip->eqn[k]) < 0) continue;
161 
162  /* See if point is outside arb */
163  for (m = 0; m < aip->neqn; m++) {
164  if (i == m || j == m || k == m)
165  continue;
166  if (VDOT(pt, aip->eqn[m])-aip->eqn[m][3] > tol->dist) {
167  next_k = 1;
168  break;
169  }
170  }
171  if (next_k != 0) continue;
172 
173  VMINMAX(stp->st_min, stp->st_max, pt);
174 
175  /* Increment "face used" counts */
176  used[i]++;
177  used[j]++;
178  used[k]++;
179  }
180  }
181  }
182 
183  /* If any planes were not used, then arbn is not convex */
184  for (i = 0; i < aip->neqn; i++) {
185  if (used[i] != 0) continue; /* face was used */
186  bu_log("arbn(%s) face %zu unused, solid is not convex\n",
187  stp->st_name, i);
188  bu_free((char *)used, "arbn used[]");
189  return -1; /* BAD */
190  }
191  bu_free((char *)used, "arbn used[]");
192 
193  stp->st_specific = (void *)aip;
194  ip->idb_ptr = ((void *)0); /* indicate we stole it */
195 
196  VADD2SCALE(stp->st_center, stp->st_min, stp->st_max, 0.5);
197  VSUB2SCALE(work, stp->st_max, stp->st_min, 0.5);
198 
199  f = work[X];
200  if (work[Y] > f) f = work[Y];
201  if (work[Z] > f) f = work[Z];
202  stp->st_aradius = f;
203  stp->st_bradius = MAGNITUDE(work);
204  return 0; /* OK */
205 }
206 
207 
208 void
209 rt_arbn_print(const struct soltab *stp)
210 {
211  size_t i;
212  struct rt_arbn_internal *arbp = (struct rt_arbn_internal *)stp->st_specific;
213 
214  RT_ARBN_CK_MAGIC(arbp);
215  bu_log("arbn bounded by %zu planes\n", arbp->neqn);
216 
217  for (i = 0; i < arbp->neqn; i++) {
218  bu_log("\t%zu: (%g, %g, %g) %g\n",
219  i,
220  INTCLAMP(arbp->eqn[i][X]), /* should have unit length */
221  INTCLAMP(arbp->eqn[i][Y]),
222  INTCLAMP(arbp->eqn[i][Z]),
223  INTCLAMP(arbp->eqn[i][W]));
224  }
225 }
226 
227 
228 /**
229  * Intersect a ray with an ARBN.
230  * Find the largest "in" distance and the smallest "out" distance.
231  * Cyrus & Beck algorithm for convex polyhedra.
232  *
233  * Returns -
234  * 0 MISS
235  * >0 HIT
236  */
237 int
238 rt_arbn_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
239 {
240  struct rt_arbn_internal *aip =
241  (struct rt_arbn_internal *)stp->st_specific;
242  int i;
243  int iplane, oplane;
244  fastf_t in, out; /* ray in/out distances */
245 
246  in = -INFINITY;
247  out = INFINITY;
248  iplane = oplane = -1;
249 
250  for (i = aip->neqn-1; i >= 0; i--) {
251  fastf_t slant_factor; /* Direction dot Normal */
252  fastf_t norm_dist;
253  fastf_t s;
254 
255  norm_dist = VDOT(aip->eqn[i], rp->r_pt) - aip->eqn[i][3];
256  if ((slant_factor = -VDOT(aip->eqn[i], rp->r_dir)) < -1.0e-10) {
257  /* exit point, when dir.N < 0. out = min(out, s) */
258  if (out > (s = norm_dist/slant_factor)) {
259  out = s;
260  oplane = i;
261  }
262  } else if (slant_factor > 1.0e-10) {
263  /* entry point, when dir.N > 0. in = max(in, s) */
264  if (in < (s = norm_dist/slant_factor)) {
265  in = s;
266  iplane = i;
267  }
268  } else {
269  /* ray is parallel to plane when dir.N == 0.
270  * If it is outside the solid, stop now
271  * Allow very small amount of slop, to catch
272  * rays that lie very nearly in the plane of a face.
273  */
274  if (norm_dist > SQRT_SMALL_FASTF)
275  return 0; /* MISS */
276  }
277  if (in > out)
278  return 0; /* MISS */
279  }
280 
281  /* Validate */
282  if (iplane == -1 || oplane == -1) {
283  bu_log("rt_arbn_shoot(%s): 1 hit => MISS\n",
284  stp->st_name);
285  return 0; /* MISS */
286  }
287  if (in >= out || out >= INFINITY)
288  return 0; /* MISS */
289 
290  {
291  struct seg *segp;
292 
293  RT_GET_SEG(segp, ap->a_resource);
294  segp->seg_stp = stp;
295  segp->seg_in.hit_dist = in;
296  segp->seg_in.hit_surfno = iplane;
297 
298  segp->seg_out.hit_dist = out;
299  segp->seg_out.hit_surfno = oplane;
300  BU_LIST_INSERT(&(seghead->l), &(segp->l));
301  }
302  return 2; /* HIT */
303 }
304 
305 
306 /**
307  * Given ONE ray distance, return the normal and entry/exit point.
308  */
309 void
310 rt_arbn_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
311 {
312  struct rt_arbn_internal *aip =
313  (struct rt_arbn_internal *)stp->st_specific;
314  size_t h;
315 
316  VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
317  h = hitp->hit_surfno;
318  if (h > aip->neqn) {
319  bu_log("rt_arbn_norm(): hit_surfno=%zu?\n", h);
320  VSETALL(hitp->hit_normal, 0);
321  return;
322  }
323  VMOVE(hitp->hit_normal, aip->eqn[h]);
324 }
325 
326 
327 /**
328  * Return the "curvature" of the ARB face.
329  * Pick a principle direction orthogonal to normal, and
330  * indicate no curvature.
331  */
332 void
333 rt_arbn_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
334 {
335  struct rt_arbn_internal *arbn = (struct rt_arbn_internal *)stp->st_specific;
336 
337  RT_ARBN_CK_MAGIC(arbn);
338 
339  bn_vec_ortho(cvp->crv_pdir, hitp->hit_normal);
340  cvp->crv_c1 = cvp->crv_c2 = 0;
341 }
342 
343 
344 /**
345  * For a hit on a face of an ARB, return the (u, v) coordinates
346  * of the hit point. 0 <= u, v <= 1.
347  * u extends along the arb_U direction defined by B-A,
348  * v extends along the arb_V direction defined by Nx(B-A).
349  */
350 void
351 rt_arbn_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
352 {
353  struct rt_arbn_internal *arbn = (struct rt_arbn_internal *)stp->st_specific;
354 
355  if (ap) RT_CK_APPLICATION(ap);
356  RT_ARBN_CK_MAGIC(arbn);
357  if (hitp) RT_CK_HIT(hitp);
358 
359  if (uvp) {
360  uvp->uv_u = uvp->uv_v = 0;
361  uvp->uv_du = uvp->uv_dv = 0;
362  }
363 }
364 
365 
366 void
367 rt_arbn_free(struct soltab *stp)
368 {
369  struct rt_arbn_internal *aip =
370  (struct rt_arbn_internal *)stp->st_specific;
371 
372  bu_free((char *)aip->eqn, "rt_arbn_internal eqn[]");
373  bu_free((char *)aip, "rt_arbn_internal");
374 }
375 
376 
377 /**
378  * Brute force through all possible plane intersections.
379  * Generate all edge lines, then intersect the line with all
380  * the other faces to find the vertices on that line.
381  * If the geometry is correct, there will be no more than two.
382  * While not the fastest strategy, this will produce an accurate
383  * plot without requiring extra bookkeeping.
384  * Note that the vectors will be drawn in no special order.
385  */
386 int
387 rt_arbn_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct rt_tess_tol *UNUSED(ttol), const struct bn_tol *tol, const struct rt_view_info *UNUSED(info))
388 {
389  struct rt_arbn_internal *aip;
390  size_t i;
391  size_t j;
392  size_t k;
393 
394  BU_CK_LIST_HEAD(vhead);
395  RT_CK_DB_INTERNAL(ip);
396  aip = (struct rt_arbn_internal *)ip->idb_ptr;
397  RT_ARBN_CK_MAGIC(aip);
398 
399  for (i = 0; i < aip->neqn - 1; i++) {
400  for (j = i + 1; j < aip->neqn; j++) {
401  double dot;
402  int point_count; /* # points on this line */
403  point_t a, b; /* start and end points */
404  vect_t dist;
405 
406  VSETALL(a, 0);
407  VSETALL(b, 0);
408 
409  /* If normals are parallel, no intersection */
410  dot = VDOT(aip->eqn[i], aip->eqn[j]);
411  if (BN_VECT_ARE_PARALLEL(dot, tol)) continue;
412 
413  /* Have an edge line, isect with all other planes */
414  point_count = 0;
415  for (k = 0; k < aip->neqn; k++) {
416  size_t m;
417  point_t pt;
418  size_t next_k;
419 
420  next_k = 0;
421 
422  if (k == i || k == j) continue;
423  if (bn_mkpoint_3planes(pt, aip->eqn[i], aip->eqn[j], aip->eqn[k]) < 0) continue;
424 
425  /* See if point is outside arb */
426  for (m = 0; m < aip->neqn; m++) {
427  if (i == m || j == m || k == m) continue;
428  if (VDOT(pt, aip->eqn[m])-aip->eqn[m][3] > tol->dist) {
429  next_k = 1;
430  break;
431  }
432  }
433 
434  if (next_k != 0) continue;
435 
436  if (point_count <= 0) {
437  RT_ADD_VLIST(vhead, pt, BN_VLIST_LINE_MOVE);
438  VMOVE(a, pt);
439  } else if (point_count == 1) {
440  VSUB2(dist, pt, a);
441  if (MAGSQ(dist) < tol->dist_sq) continue;
442  RT_ADD_VLIST(vhead, pt, BN_VLIST_LINE_DRAW);
443  VMOVE(b, pt);
444  } else {
445  VSUB2(dist, pt, a);
446  if (MAGSQ(dist) < tol->dist_sq) continue;
447  VSUB2(dist, pt, b);
448  if (MAGSQ(dist) < tol->dist_sq) continue;
449  bu_log("rt_arbn_plot() error, point_count=%d (>2) on edge %zu/%zu, non-convex\n",
450  point_count+1,
451  i, j);
452  VPRINT(" a", a);
453  VPRINT(" b", b);
454  VPRINT("pt", pt);
455  RT_ADD_VLIST(vhead, pt, BN_VLIST_LINE_DRAW); /* draw it */
456  }
457  point_count++;
458  }
459  /* Point counts of 1 are (generally) not harmful,
460  * occurring on pyramid peaks and the like.
461  */
462  }
463  }
464  return 0;
465 }
466 
467 
468 /* structures used by arbn tessellator */
469 struct arbn_pts
470 {
471  point_t pt; /* coordinates for vertex */
472  int plane_no[3]; /* which planes intersect here */
473  struct vertex **vp; /* pointer to vertex struct pointer for NMG's */
474 };
475 
476 
478 {
479  int v1_no, v2_no; /* index into arbn_pts for endpoints of edge */
480 };
481 
482 
483 #define LOC(i, j) i*(aip->neqn)+j
484 
485 static void
486 Sort_edges(struct arbn_edges *edges, size_t *edge_count, const struct rt_arbn_internal *aip)
487 {
488  size_t face;
489 
490  for (face = 0; face < aip->neqn; face++) {
491  int done = 0;
492  size_t edge1, edge2;
493 
494  if (edge_count[face] < 3)
495  continue; /* nothing to sort */
496 
497  edge1 = 0;
498  edge2 = 0;
499  while (!done) {
500  size_t edge3;
501  size_t tmp;
502 
503  /* Look for out of order edge (edge2) */
504  while (++edge2 < edge_count[face] &&
505  edges[LOC(face, edge1)].v2_no == edges[LOC(face, edge2)].v1_no)
506  edge1++;
507  if (edge2 == edge_count[face]) {
508  /* all edges are in order */
509  done = 1;
510  continue;
511  }
512 
513  /* look for edge (edge3) that belongs where edge2 is */
514  edge3 = edge2 - 1;
515  while (++edge3 < edge_count[face] &&
516  edges[LOC(face, edge1)].v2_no != edges[LOC(face, edge3)].v1_no &&
517  edges[LOC(face, edge1)].v2_no != edges[LOC(face, edge3)].v2_no);
518 
519  if (edge3 == edge_count[face])
520  bu_bomb("rt_arbn_tess: Sort_edges: Cannot find next edge in loop\n");
521 
522  if (edge2 != edge3) {
523  /* swap edge2 and edge3 */
524  tmp = edges[LOC(face, edge2)].v1_no;
525  edges[LOC(face, edge2)].v1_no = edges[LOC(face, edge3)].v1_no;
526  edges[LOC(face, edge3)].v1_no = tmp;
527  tmp = edges[LOC(face, edge2)].v2_no;
528  edges[LOC(face, edge2)].v2_no = edges[LOC(face, edge3)].v2_no;
529  edges[LOC(face, edge3)].v2_no = tmp;
530  }
531  if (edges[LOC(face, edge1)].v2_no == edges[LOC(face, edge2)].v2_no) {
532  /* reverse order of edge */
533  tmp = edges[LOC(face, edge2)].v1_no;
534  edges[LOC(face, edge2)].v1_no = edges[LOC(face, edge2)].v2_no;
535  edges[LOC(face, edge2)].v2_no = tmp;
536  }
537 
538  edge1 = edge2;
539  }
540  }
541 }
542 
543 
544 /**
545  * "Tessellate" an ARB into an NMG data structure.
546  * Purely a mechanical transformation of one faceted object
547  * into another.
548  *
549  * Returns -
550  * -1 failure
551  * 0 OK. *r points to nmgregion that holds this tessellation.
552  */
553 int
554 rt_arbn_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *UNUSED(ttol), const struct bn_tol *tol)
555 {
556  struct rt_arbn_internal *aip;
557  struct shell *s;
558  struct faceuse **fu; /* array of faceuses */
559  size_t nverts; /* maximum possible number of vertices = neqn!/(3!(neqn-3)! */
560  size_t point_count = 0; /* actual number of vertices */
561  size_t face_count = 0; /* actual number of faces built */
562  size_t i, j, k, l, n;
563  struct arbn_pts *pts;
564  struct arbn_edges *edges; /* A list of edges for each plane eqn (each face) */
565  size_t *edge_count; /* number of edges for each face */
566  size_t max_edge_count; /* maximum number of edges for any face */
567  struct vertex **verts; /* Array of pointers to vertex structs */
568  struct vertex ***loop_verts; /* Array of pointers to vertex structs to pass to nmg_cmface */
569 
570  RT_CK_DB_INTERNAL(ip);
571  aip = (struct rt_arbn_internal *)ip->idb_ptr;
572  RT_ARBN_CK_MAGIC(aip);
573 
574  /* Allocate memory for the vertices */
575  nverts = aip->neqn * (aip->neqn-1) * (aip->neqn-2) / 6;
576  pts = (struct arbn_pts *)bu_calloc(nverts, sizeof(struct arbn_pts), "rt_arbn_tess: pts");
577 
578  /* Allocate memory for arbn_edges */
579  edges = (struct arbn_edges *)bu_calloc(aip->neqn*aip->neqn, sizeof(struct arbn_edges) ,
580  "rt_arbn_tess: edges");
581  edge_count = (size_t *)bu_calloc(aip->neqn, sizeof(size_t), "rt_arbn_tess: edge_count");
582 
583  /* Allocate memory for faceuses */
584  fu = (struct faceuse **)bu_calloc(aip->neqn, sizeof(struct faceuse *), "rt_arbn_tess: fu");
585 
586  /* Calculate all vertices */
587  for (i = 0; i < aip->neqn; i++) {
588  for (j = i + 1; j < aip->neqn; j++) {
589  for (k = j + 1; k < aip->neqn; k++) {
590  int keep_point = 1;
591 
592  if (bn_mkpoint_3planes(pts[point_count].pt, aip->eqn[i], aip->eqn[j], aip->eqn[k]))
593  continue;
594 
595  for (l = 0; l < aip->neqn; l++) {
596  if (l == i || l == j || l == k)
597  continue;
598  if (DIST_PT_PLANE(pts[point_count].pt, aip->eqn[l]) > tol->dist) {
599  keep_point = 0;
600  break;
601  }
602  }
603  if (keep_point) {
604  pts[point_count].plane_no[0] = i;
605  pts[point_count].plane_no[1] = j;
606  pts[point_count].plane_no[2] = k;
607  point_count++;
608  }
609  }
610  }
611  }
612 
613  /* Allocate memory for the NMG vertex pointers */
614  verts = (struct vertex **)bu_calloc(point_count, sizeof(struct vertex *) ,
615  "rt_arbn_tess: verts");
616 
617  /* Associate points with vertices */
618  for (i = 0; i < point_count; i++)
619  pts[i].vp = &verts[i];
620 
621  /* Check for duplicate points */
622  for (i = 0; i < point_count; i++) {
623  for (j = i + 1; j < point_count; j++) {
624  if (DIST_PT_PT_SQ(pts[i].pt, pts[j].pt) < tol->dist_sq) {
625  /* These two points should point to the same vertex */
626  pts[j].vp = pts[i].vp;
627  }
628  }
629  }
630 
631  /* Make list of edges for each face */
632  for (i = 0; i < aip->neqn; i++) {
633  /* look for a point that lies in this face */
634  for (j = 0; j < point_count; j++) {
635  if (pts[j].plane_no[0] != (int)i
636  && pts[j].plane_no[1] != (int)i
637  && pts[j].plane_no[2] != (int)i)
638  continue;
639 
640  /* look for another point that shares plane "i" and another with this one */
641  for (k = j + 1; k < point_count; k++) {
642  size_t match = (size_t)-1;
643  size_t pt1, pt2;
644  int duplicate = 0;
645 
646  /* skip points not on plane "i" */
647  if (pts[k].plane_no[0] != (int)i
648  && pts[k].plane_no[1] != (int)i
649  && pts[k].plane_no[2] != (int)i)
650  continue;
651 
652  for (l = 0; l < 3; l++) {
653  for (n = 0; n < 3; n++) {
654  if (pts[j].plane_no[l] == pts[k].plane_no[n]
655  && pts[j].plane_no[l] != (int)i) {
656  match = pts[j].plane_no[l];
657  break;
658  }
659  }
660  if (match != (size_t)-1)
661  break;
662  }
663 
664  if (match == (size_t)-1)
665  continue;
666 
667  /* convert equivalent points to lowest point number */
668  pt1 = j;
669  pt2 = k;
670  for (l = 0; l < pt1; l++) {
671  if (pts[pt1].vp == pts[l].vp) {
672  pt1 = l;
673  break;
674  }
675  }
676  for (l = 0; l < pt2; l++) {
677  if (pts[pt2].vp == pts[l].vp) {
678  pt2 = l;
679  break;
680  }
681  }
682 
683  /* skip null edges */
684  if (pt1 == pt2)
685  continue;
686 
687  /* check for duplicate edge */
688  for (l = 0; l < edge_count[i]; l++) {
689  if ((edges[LOC(i, l)].v1_no == (int)pt1
690  && edges[LOC(i, l)].v2_no == (int)pt2)
691  || (edges[LOC(i, l)].v2_no == (int)pt1
692  && edges[LOC(i, l)].v1_no == (int)pt2))
693  {
694  duplicate = 1;
695  break;
696  }
697  }
698  if (duplicate)
699  continue;
700 
701  /* found an edge belonging to faces "i" and "match" */
702  if (edge_count[i] == aip->neqn) {
703  bu_log("Too many edges found for one face\n");
704  goto fail;
705  }
706  edges[LOC(i, edge_count[i])].v1_no = pt1;
707  edges[LOC(i, edge_count[i])].v2_no = pt2;
708  edge_count[i]++;
709  }
710  }
711  }
712 
713  /* for each face, sort the list of edges into a loop */
714  Sort_edges(edges, edge_count, aip);
715 
716  /* Get max number of edges for any face */
717  max_edge_count = 0;
718  for (i = 0; i < aip->neqn; i++)
719  if (edge_count[i] > max_edge_count)
720  max_edge_count = edge_count[i];
721 
722  /* Allocate memory for array to pass to nmg_cmface */
723  loop_verts = (struct vertex ***) bu_calloc(max_edge_count, sizeof(struct vertex **) ,
724  "rt_arbn_tess: loop_verts");
725 
726  *r = nmg_mrsv(m); /* Make region, empty shell, vertex */
727  s = BU_LIST_FIRST(shell, &(*r)->s_hd);
728 
729  /* Make the faces */
730  for (i = 0; i < aip->neqn; i++) {
731  int loop_length = 0;
732 
733  for (j = 0; j < edge_count[i]; j++) {
734  /* skip zero length edges */
735  if (pts[edges[LOC(i, j)].v1_no].vp == pts[edges[LOC(i, j)].v2_no].vp)
736  continue;
737 
738  /* put vertex pointers into loop_verts array */
739  loop_verts[loop_length] = pts[edges[LOC(i, j)].v2_no].vp;
740  loop_length++;
741  }
742 
743  /* Make the face if there is are least 3 vertices */
744  if (loop_length > 2)
745  fu[face_count++] = nmg_cmface(s, loop_verts, loop_length);
746  }
747 
748  /* Associate vertex geometry */
749  for (i = 0; i < point_count; i++) {
750  if (!(*pts[i].vp))
751  continue;
752 
753  if ((*pts[i].vp)->vg_p)
754  continue;
755 
756  nmg_vertex_gv(*pts[i].vp, pts[i].pt);
757  }
758 
759  bu_free((char *)pts, "rt_arbn_tess: pts");
760  bu_free((char *)edges, "rt_arbn_tess: edges");
761  bu_free((char *)edge_count, "rt_arbn_tess: edge_count");
762  bu_free((char *)verts, "rt_arbn_tess: verts");
763  bu_free((char *)loop_verts, "rt_arbn_tess: loop_verts");
764 
765  /* Associate face geometry */
766  for (i = 0; i < face_count; i++) {
767  if (nmg_fu_planeeqn(fu[i], tol)) {
768  bu_log("Failed to calculate face plane equation\n");
769  bu_free((char *)fu, "rt_arbn_tess: fu");
770  nmg_kr(*r);
771  *r = (struct nmgregion *)NULL;
772  return -1;
773  }
774  }
775 
776  bu_free((char *)fu, "rt_arbn_tess: fu");
777 
778  nmg_fix_normals(s, tol);
779 
780  (void)nmg_mark_edges_real(&s->l.magic);
781 
782  /* Compute "geometry" for region and shell */
783  nmg_region_a(*r, tol);
784 
785  return 0;
786 
787 fail:
788  bu_free((char *)pts, "rt_arbn_tess: pts");
789  bu_free((char *)edges, "rt_arbn_tess: edges");
790  bu_free((char *)edge_count, "rt_arbn_tess: edge_count");
791  bu_free((char *)verts, "rt_arbn_tess: verts");
792  return -1;
793 }
794 
795 
796 /**
797  * Convert from "network" doubles to machine specific.
798  * Transform
799  */
800 int
801 rt_arbn_import4(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
802 {
803  union record *rp;
804  struct rt_arbn_internal *aip;
805  size_t i;
806 
807  /* must be double for import and export */
808  double *scan;
809 
810  if (dbip) RT_CK_DBI(dbip);
811 
812  BU_CK_EXTERNAL(ep);
813  rp = (union record *)ep->ext_buf;
814  if (rp->u_id != DBID_ARBN) {
815  bu_log("rt_arbn_import4: defective record, id=x%x\n", rp->u_id);
816  return -1;
817  }
818 
819  RT_CK_DB_INTERNAL(ip);
820  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
821  ip->idb_type = ID_ARBN;
822  ip->idb_meth = &OBJ[ID_ARBN];
823  BU_ALLOC(ip->idb_ptr, struct rt_arbn_internal);
824 
825  aip = (struct rt_arbn_internal *)ip->idb_ptr;
826  aip->magic = RT_ARBN_INTERNAL_MAGIC;
827  aip->neqn = ntohl(*(uint32_t *)rp->n.n_neqn);
828  if (aip->neqn <= 0) return -1;
829 
830  aip->eqn = (plane_t *)bu_malloc(aip->neqn*sizeof(plane_t), "arbn plane eqn[]");
831  scan = (double *)bu_malloc(aip->neqn*sizeof(double)*ELEMENTS_PER_PLANE, "scan array");
832 
833  bu_cv_ntohd((unsigned char *)scan, (unsigned char *)(&rp[1]), aip->neqn*ELEMENTS_PER_PLANE);
834  for (i = 0; i < aip->neqn; i++) {
835  aip->eqn[i][X] = scan[(i*ELEMENTS_PER_PLANE)+0]; /* convert double to fastf_t */
836  aip->eqn[i][Y] = scan[(i*ELEMENTS_PER_PLANE)+1]; /* convert double to fastf_t */
837  aip->eqn[i][Z] = scan[(i*ELEMENTS_PER_PLANE)+2]; /* convert double to fastf_t */
838  aip->eqn[i][W] = scan[(i*ELEMENTS_PER_PLANE)+3]; /* convert double to fastf_t */
839  }
840  bu_free(scan, "scan array");
841 
842  /* Transform by the matrix */
843  if (mat == NULL) mat = bn_mat_identity;
844  for (i = 0; i < aip->neqn; i++) {
845  point_t orig_pt;
846  point_t pt;
847  vect_t norm;
848  fastf_t factor;
849 
850  /* unitize the plane equation first */
851  factor = 1.0 / MAGNITUDE(aip->eqn[i]);
852  VSCALE(aip->eqn[i], aip->eqn[i], factor);
853  aip->eqn[i][W] = aip->eqn[i][W] * factor;
854 
855 
856  /* Pick a point on the original halfspace */
857  VSCALE(orig_pt, aip->eqn[i], aip->eqn[i][W]);
858 
859  /* Transform the point, and the normal */
860  MAT4X3VEC(norm, mat, aip->eqn[i]);
861  MAT4X3PNT(pt, mat, orig_pt);
862 
863  /* Measure new distance from origin to new point */
864  VUNITIZE(norm);
865  VMOVE(aip->eqn[i], norm);
866  aip->eqn[i][W] = VDOT(pt, norm);
867  }
868 
869  return 0;
870 }
871 
872 
873 int
874 rt_arbn_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
875 {
876  struct rt_arbn_internal *aip;
877  union record *rec;
878  size_t ngrans;
879  size_t i;
880 
881  /* scaling buffer must be double, not fastf_t */
882  double *sbuf;
883  double *sp;
884 
885  if (dbip) RT_CK_DBI(dbip);
886 
887  RT_CK_DB_INTERNAL(ip);
888  if (ip->idb_type != ID_ARBN) return -1;
889  aip = (struct rt_arbn_internal *)ip->idb_ptr;
890  RT_ARBN_CK_MAGIC(aip);
891 
892  if (aip->neqn <= 0) return -1;
893 
894  /*
895  * The network format for a double is 8 bytes and there are 4
896  * doubles per plane equation.
897  */
898  ngrans = (aip->neqn * 8 * ELEMENTS_PER_PLANE + sizeof(union record)-1) /
899  sizeof(union record);
900 
901  BU_CK_EXTERNAL(ep);
902  ep->ext_nbytes = (ngrans + 1) * sizeof(union record);
903  ep->ext_buf = (uint8_t *)bu_calloc(1, ep->ext_nbytes, "arbn external");
904  rec = (union record *)ep->ext_buf;
905 
906  rec[0].n.n_id = DBID_ARBN;
907  *(uint32_t *)rec[0].n.n_neqn = htonl(aip->neqn);
908  *(uint32_t *)rec[0].n.n_grans = htonl(ngrans);
909 
910  /* Take the data from the caller, and scale it, into sbuf */
911  sp = sbuf = (double *)bu_malloc(
912  aip->neqn * sizeof(double) * ELEMENTS_PER_PLANE, "arbn temp");
913  for (i = 0; i < aip->neqn; i++) {
914  /* Normal is unscaled, should have unit length; d is scaled */
915  *sp++ = aip->eqn[i][X];
916  *sp++ = aip->eqn[i][Y];
917  *sp++ = aip->eqn[i][Z];
918  *sp++ = aip->eqn[i][W] * local2mm;
919  }
920 
921  bu_cv_htond((unsigned char *)&rec[1], (unsigned char *)sbuf, aip->neqn * ELEMENTS_PER_PLANE);
922 
923  bu_free((char *)sbuf, "arbn temp");
924  return 0; /* OK */
925 }
926 
927 
928 /**
929  * Convert from "network" doubles to machine specific.
930  * Transform
931  */
932 int
933 rt_arbn_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
934 {
935  struct rt_arbn_internal *aip;
936  size_t i;
937  unsigned long neqn;
938  int double_count;
939  size_t byte_count;
940 
941  /* must be double for import and export */
942  double *eqn;
943 
944  RT_CK_DB_INTERNAL(ip);
945  BU_CK_EXTERNAL(ep);
946  if (dbip) RT_CK_DBI(dbip);
947 
948  neqn = ntohl(*(uint32_t *)ep->ext_buf);
949  double_count = neqn * ELEMENTS_PER_PLANE;
950  byte_count = double_count * SIZEOF_NETWORK_DOUBLE;
951 
952  BU_ASSERT_LONG(ep->ext_nbytes, ==, 4+ byte_count);
953 
954  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
955  ip->idb_type = ID_ARBN;
956  ip->idb_meth = &OBJ[ID_ARBN];
957  BU_ALLOC(ip->idb_ptr, struct rt_arbn_internal);
958 
959  aip = (struct rt_arbn_internal *)ip->idb_ptr;
960  aip->magic = RT_ARBN_INTERNAL_MAGIC;
961  aip->neqn = neqn;
962  if (aip->neqn <= 0) return -1;
963 
964  eqn = (double *)bu_malloc(byte_count, "arbn plane eqn[] temp buf");
965  bu_cv_ntohd((unsigned char *)eqn, (unsigned char *)ep->ext_buf + ELEMENTS_PER_PLANE, double_count);
966  aip->eqn = (plane_t *)bu_malloc(double_count * sizeof(fastf_t), "arbn plane eqn[]");
967  for (i = 0; i < aip->neqn; i++) {
968  HMOVE(aip->eqn[i], &eqn[i*ELEMENTS_PER_PLANE]);
969  }
970  bu_free(eqn, "arbn plane eqn[] temp buf");
971 
972  /* Transform by the matrix, if we have one that is not the identity */
973  if (mat && !bn_mat_is_identity(mat)) {
974  for (i = 0; i < aip->neqn; i++) {
975  point_t orig_pt;
976  point_t pt;
977  vect_t norm;
978  fastf_t factor;
979 
980  /* unitize the plane equation first */
981  factor = 1.0 / MAGNITUDE(aip->eqn[i]);
982  VSCALE(aip->eqn[i], aip->eqn[i], factor);
983  aip->eqn[i][W] = aip->eqn[i][W] * factor;
984 
985  /* Pick a point on the original halfspace */
986  VSCALE(orig_pt, aip->eqn[i], aip->eqn[i][W]);
987 
988  /* Transform the point, and the normal */
989  MAT4X3VEC(norm, mat, aip->eqn[i]);
990  MAT4X3PNT(pt, mat, orig_pt);
991 
992  /* Measure new distance from origin to new point */
993  VUNITIZE(norm);
994  VMOVE(aip->eqn[i], norm);
995  aip->eqn[i][W] = VDOT(pt, norm);
996  }
997  }
998 
999  return 0;
1000 }
1001 
1002 
1003 int
1004 rt_arbn_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1005 {
1006  struct rt_arbn_internal *aip;
1007  size_t i;
1008  int double_count;
1009  int byte_count;
1010 
1011  /* must be double for export */
1012  double *vec;
1013  double *sp;
1014 
1015  RT_CK_DB_INTERNAL(ip);
1016  if (dbip) RT_CK_DBI(dbip);
1017 
1018  if (ip->idb_type != ID_ARBN) return -1;
1019  aip = (struct rt_arbn_internal *)ip->idb_ptr;
1020  RT_ARBN_CK_MAGIC(aip);
1021 
1022  if (aip->neqn <= 0) return -1;
1023 
1024  double_count = aip->neqn * ELEMENTS_PER_PLANE;
1025  byte_count = double_count * SIZEOF_NETWORK_DOUBLE;
1026 
1027  BU_CK_EXTERNAL(ep);
1028  ep->ext_nbytes = SIZEOF_NETWORK_LONG + byte_count;
1029  ep->ext_buf = (uint8_t *)bu_malloc(ep->ext_nbytes, "arbn external");
1030 
1031  *(uint32_t *)ep->ext_buf = htonl(aip->neqn);
1032 
1033  /* Take the data from the caller, and scale it, into vec */
1034  sp = vec = (double *)bu_malloc(byte_count, "arbn temp");
1035  for (i = 0; i < aip->neqn; i++) {
1036  /* Normal is unscaled, should have unit length; d is scaled */
1037  *sp++ = aip->eqn[i][X];
1038  *sp++ = aip->eqn[i][Y];
1039  *sp++ = aip->eqn[i][Z];
1040  *sp++ = aip->eqn[i][W] * local2mm;
1041  }
1042 
1043  /* Convert from internal (host) to database (network) format */
1044  bu_cv_htond((unsigned char *)ep->ext_buf + SIZEOF_NETWORK_LONG, (unsigned char *)vec, double_count);
1045 
1046  bu_free((char *)vec, "arbn temp");
1047  return 0; /* OK */
1048 }
1049 
1050 
1051 /**
1052  * Make human-readable formatted presentation of this solid.
1053  * First line describes type of solid.
1054  * Additional lines are indented one tab, and give parameter values.
1055  */
1056 int
1057 rt_arbn_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
1058 {
1059  struct rt_arbn_internal *aip =
1060  (struct rt_arbn_internal *)ip->idb_ptr;
1061  char buf[256];
1062  size_t i;
1063 
1064  RT_ARBN_CK_MAGIC(aip);
1065  sprintf(buf, "arbn bounded by %lu planes\n", (long unsigned)aip->neqn);
1066  bu_vls_strcat(str, buf);
1067 
1068  if (!verbose) return 0;
1069 
1070  for (i = 0; i < aip->neqn; i++) {
1071  sprintf(buf, "\t%lu: (%g, %g, %g) %g\n",
1072  (long unsigned)i,
1073  INTCLAMP(aip->eqn[i][X]), /* should have unit length */
1074  INTCLAMP(aip->eqn[i][Y]),
1075  INTCLAMP(aip->eqn[i][Z]),
1076  INTCLAMP(aip->eqn[i][W] * mm2local));
1077  bu_vls_strcat(str, buf);
1078  }
1079  return 0;
1080 }
1081 
1082 
1083 /**
1084  * Free the storage associated with the rt_db_internal version of this solid.
1085  */
1086 void
1088 {
1089  struct rt_arbn_internal *aip;
1090 
1091  RT_CK_DB_INTERNAL(ip);
1092  aip = (struct rt_arbn_internal *)ip->idb_ptr;
1093  RT_ARBN_CK_MAGIC(aip);
1094 
1095  if (aip->neqn > 0)
1096  bu_free((char *)aip->eqn, "rt_arbn_internal eqn[]");
1097  bu_free((char *)aip, "rt_arbn_internal");
1098 
1099  ip->idb_ptr = (void *)0; /* sanity */
1100 }
1101 
1102 
1103 /**
1104  * Routine to format the parameters of an ARBN primitive for "db get"
1105  *
1106  * Legal requested parameters include:
1107  * "N" - number of equations
1108  * "P" - list of all the planes
1109  * "P#" - the specified plane number (0 based)
1110  * no arguments returns everything
1111  */
1112 int
1113 rt_arbn_get(struct bu_vls *logstr, const struct rt_db_internal *intern, const char *attr)
1114 {
1115  struct rt_arbn_internal *arbn=(struct rt_arbn_internal *)intern->idb_ptr;
1116  size_t i;
1117  long val;
1118 
1119  RT_ARBN_CK_MAGIC(arbn);
1120 
1121  if (attr == (char *)NULL) {
1122  bu_vls_strcpy(logstr, "arbn");
1123  bu_vls_printf(logstr, " N %zu", arbn->neqn);
1124  for (i = 0; i < arbn->neqn; i++) {
1125  bu_vls_printf(logstr, " P%zu {%.25g %.25g %.25g %.25g}", i,
1126  V4ARGS(arbn->eqn[i]));
1127  }
1128  } else if (BU_STR_EQUAL(attr, "N")) {
1129  bu_vls_printf(logstr, "%zu", arbn->neqn);
1130  } else if (BU_STR_EQUAL(attr, "P")) {
1131  for (i = 0; i < arbn->neqn; i++) {
1132  bu_vls_printf(logstr, " P%zu {%.25g %.25g %.25g %.25g}", i,
1133  V4ARGS(arbn->eqn[i]));
1134  }
1135  } else if (attr[0] == 'P') {
1136  if (isdigit((int)attr[1]) == 0) {
1137  bu_vls_printf(logstr, "ERROR: Illegal plane number\n");
1138  return BRLCAD_ERROR;
1139  }
1140 
1141  val = atol(&attr[1]);
1142  if (val < 0 || (size_t)val >= arbn->neqn) {
1143  bu_vls_printf(logstr, "ERROR: Illegal plane number [%ld]\n", val);
1144  return BRLCAD_ERROR;
1145  }
1146  i = (size_t)val;
1147 
1148  bu_vls_printf(logstr, "%.25g %.25g %.25g %.25g", V4ARGS(arbn->eqn[i]));
1149  } else {
1150  bu_vls_printf(logstr, "ERROR: Unknown attribute, choices are N, P, or P#\n");
1151  return BRLCAD_ERROR;
1152  }
1153 
1154  return BRLCAD_OK;
1155 }
1156 
1157 
1158 /**
1159  * Routine to modify an arbn via the "db adjust" command
1160  *
1161  * Legal parameters are:
1162  * "N" - adjust the number of planes (new ones will be zeroed)
1163  * "P" - adjust the entire list of planes
1164  * "P#" - adjust a specific plane (0 based)
1165  * "P+" - add a new plane to the list of planes
1166  */
1167 int
1168 rt_arbn_adjust(struct bu_vls *logstr, struct rt_db_internal *intern, int argc, const char **argv)
1169 {
1170  struct rt_arbn_internal *arbn;
1171  unsigned char *c;
1172  int len;
1173  size_t i, j;
1174  long val;
1175  fastf_t *new_planes;
1176  fastf_t *array;
1177 
1178  RT_CK_DB_INTERNAL(intern);
1179 
1180  arbn = (struct rt_arbn_internal *)intern->idb_ptr;
1181  RT_ARBN_CK_MAGIC(arbn);
1182 
1183  while (argc >= 2) {
1184  if (BU_STR_EQUAL(argv[0], "N")) {
1185  val = atol(argv[1]);
1186  if (val < 0) {
1187  bu_vls_printf(logstr, "ERROR: number of planes [%ld] must be greater than 0\n", val);
1188  val = 1;
1189  }
1190 
1191  i = (size_t)val;
1192  if (i == arbn->neqn)
1193  goto cont;
1194 
1195  arbn->eqn = (plane_t *)bu_realloc(arbn->eqn,
1196  i * sizeof(plane_t),
1197  "arbn->eqn");
1198  for (j=arbn->neqn; j<i; j++) {
1199  HSETALL(arbn->eqn[j], 0.0);
1200  }
1201  arbn->neqn = i;
1202 
1203  } else if (BU_STR_EQUAL(argv[0], "P")) {
1204  /* eliminate all the '{' and '}' chars */
1205  c = (unsigned char *)argv[1];
1206  while (*c != '\0') {
1207  if (*c == '{' || *c == '}')
1208  *c = ' ';
1209  c++;
1210  }
1211  len = 0;
1212  (void)tcl_list_to_fastf_array(brlcad_interp, argv[1], &new_planes, &len);
1213 
1214  if (len%ELEMENTS_PER_PLANE) {
1215  bu_vls_printf(logstr,
1216  "ERROR: Incorrect number of plane coefficients\n");
1217  if (len)
1218  bu_free((char *)new_planes, "new_planes");
1219  return BRLCAD_ERROR;
1220  }
1221  if (arbn->eqn)
1222  bu_free((char *)arbn->eqn, "arbn->eqn");
1223  arbn->eqn = (plane_t *)new_planes;
1224  arbn->neqn = (size_t)len / ELEMENTS_PER_PLANE;
1225  for (i = 0; i < arbn->neqn; i++)
1226  VUNITIZE(arbn->eqn[i]);
1227  } else if (argv[0][0] == 'P') {
1228  if (argv[0][1] == '+') {
1229  i = arbn->neqn;
1230  arbn->neqn++;
1231  arbn->eqn = (plane_t *)bu_realloc(arbn->eqn,
1232  (arbn->neqn) * sizeof(plane_t),
1233  "arbn->eqn");
1234  } else if (isdigit((int)argv[0][1])) {
1235  i = atoi(&argv[0][1]);
1236  } else {
1237  bu_vls_printf(logstr,
1238  "ERROR: illegal argument, choices are P, P#, P+, or N\n");
1239  return TCL_ERROR;
1240  }
1241  if (i >= arbn->neqn) {
1242  bu_vls_printf(logstr, "ERROR: plane number out of range\n");
1243  return BRLCAD_ERROR;
1244  }
1245  len = ELEMENTS_PER_PLANE;
1246  array = (fastf_t *)&arbn->eqn[i];
1248  &array, &len) != ELEMENTS_PER_PLANE) {
1249  bu_vls_printf(logstr,
1250  "ERROR: incorrect number of coefficients for a plane\n");
1251  return BRLCAD_ERROR;
1252  }
1253  VUNITIZE(arbn->eqn[i]);
1254  } else {
1255  bu_vls_printf(logstr,
1256  "ERROR: illegal argument, choices are P, P#, P+, or N\n");
1257  return BRLCAD_ERROR;
1258  }
1259  cont:
1260  argc -= 2;
1261  argv += 2;
1262  }
1263  return BRLCAD_OK;
1264 }
1265 
1266 
1267 int
1268 rt_arbn_params(struct pc_pc_set *UNUSED(ps), const struct rt_db_internal *ip)
1269 {
1270  struct rt_arbn_internal *aip;
1271 
1272  RT_CK_DB_INTERNAL(ip);
1273  aip = (struct rt_arbn_internal *)ip->idb_ptr;
1274  RT_ARBN_CK_MAGIC(aip);
1275 
1276  return 0; /* OK */
1277 }
1278 
1279 
1280 /* contains information used to analyze a polygonal face */
1282 {
1283  char label[5];
1284  size_t npts;
1285  point_t *pts;
1286  plane_t plane_eqn;
1289  point_t cent_pyramid;
1290  point_t cent;
1291 };
1292 
1293 
1294 static void
1295 rt_arbn_faces_area(struct poly_face* faces, struct rt_arbn_internal* aip)
1296 {
1297  size_t i;
1298  size_t *npts = (size_t *)bu_calloc(aip->neqn, sizeof(size_t), "rt_arbn_faces_area: npts");
1299  point_t **tmp_pts = (point_t **)bu_calloc(aip->neqn, sizeof(point_t *), "rt_arbn_faces_area: tmp_pts");
1300  plane_t *eqs = (plane_t *)bu_calloc(aip->neqn, sizeof(plane_t), "rt_arbn_faces_area: eqs");
1301 
1302  for (i = 0; i < aip->neqn; i++) {
1303  HMOVE(faces[i].plane_eqn, aip->eqn[i]);
1304  VUNITIZE(faces[i].plane_eqn);
1305  tmp_pts[i] = faces[i].pts;
1306  HMOVE(eqs[i], faces[i].plane_eqn);
1307  }
1308  bn_polygon_mk_pts_planes(npts, tmp_pts, aip->neqn, (const plane_t *)eqs);
1309  for (i = 0; i < aip->neqn; i++) {
1310  faces[i].npts = npts[i];
1311  bn_polygon_sort_ccw(faces[i].npts, faces[i].pts, faces[i].plane_eqn);
1312  bn_polygon_area(&faces[i].area, faces[i].npts, (const point_t *)faces[i].pts);
1313  }
1314  bu_free((char *)tmp_pts, "rt_arbn_faces_area: tmp_pts");
1315  bu_free((char *)npts, "rt_arbn_faces_area: npts");
1316  bu_free((char *)eqs, "rt_arbn_faces_area: eqs");
1317 }
1318 
1319 
1320 void
1321 rt_arbn_surf_area(fastf_t *area, const struct rt_db_internal *ip)
1322 {
1323  struct poly_face *faces;
1324  struct rt_arbn_internal *aip = (struct rt_arbn_internal *)ip->idb_ptr;
1325  size_t i;
1326 
1327  /* allocate array of face structs */
1328  faces = (struct poly_face *)bu_calloc(aip->neqn, sizeof(struct poly_face), "arbn_surf_area: faces");
1329  for (i = 0; i < aip->neqn; i++) {
1330  /* allocate array of pt structs, max number of verts per faces = (# of faces) - 1 */
1331  faces[i].pts = (point_t *)bu_calloc(aip->neqn - 1, sizeof(point_t), "arbn_surf_area: pts");
1332  }
1333  rt_arbn_faces_area(faces, aip);
1334  for (i = 0; i < aip->neqn; i++) {
1335  *area += faces[i].area;
1336  bu_free((char *)faces[i].pts, "rt_arbn_surf_area: pts");
1337  }
1338  bu_free((char *)faces, "rt_arbn_surf_area: faces");
1339 }
1340 
1341 
1342 void
1343 rt_arbn_centroid(point_t *cent, const struct rt_db_internal *ip)
1344 {
1345  struct poly_face *faces;
1346  struct rt_arbn_internal *aip = (struct rt_arbn_internal *)ip->idb_ptr;
1347  size_t i;
1348  point_t arbit_point = VINIT_ZERO;
1349  fastf_t volume = 0.0;
1350 
1351  *cent[0] = 0.0;
1352  *cent[1] = 0.0;
1353  *cent[2] = 0.0;
1354 
1355  if (cent == NULL)
1356  return;
1357 
1358  /* allocate array of face structs */
1359  faces = (struct poly_face *)bu_calloc(aip->neqn, sizeof(struct poly_face), "rt_arbn_centroid: faces");
1360  for (i = 0; i < aip->neqn; i++) {
1361  /* allocate array of pt structs, max number of verts per faces = (# of faces) - 1 */
1362  faces[i].pts = (point_t *)bu_calloc(aip->neqn - 1, sizeof(point_t), "rt_arbn_centroid: pts");
1363  }
1364  rt_arbn_faces_area(faces, aip);
1365  for (i = 0; i < aip->neqn; i++) {
1366  bn_polygon_centroid(&faces[i].cent, faces[i].npts, (const point_t *) faces[i].pts);
1367  VADD2(arbit_point, arbit_point, faces[i].cent);
1368 
1369  }
1370  VSCALE(arbit_point, arbit_point, (1/aip->neqn));
1371 
1372  for (i = 0; i < aip->neqn; i++) {
1373  vect_t tmp = VINIT_ZERO;
1374  /* calculate volume */
1375  VSCALE(tmp, faces[i].plane_eqn, faces[i].area);
1376  faces[i].vol_pyramid = (VDOT(faces[i].pts[0], tmp)/3);
1377  volume += faces[i].vol_pyramid;
1378  /*Vector from arbit_point to centroid of face, results in h of pyramid */
1379  VSUB2(faces[i].cent_pyramid, faces[i].cent, arbit_point);
1380  /*centroid of pyramid is 1/4 up from the bottom */
1381  VSCALE(faces[i].cent_pyramid, faces[i].cent_pyramid, 0.75f);
1382  /* now cent_pyramid is back in the polyhedron */
1383  VADD2(faces[i].cent_pyramid, faces[i].cent_pyramid, arbit_point);
1384  /* weight centroid of pyramid by pyramid's volume */
1385  VSCALE(faces[i].cent_pyramid, faces[i].cent_pyramid, faces[i].vol_pyramid);
1386  /* add cent_pyramid to the centroid of the polyhedron */
1387  VADD2(*cent, *cent, faces[i].cent_pyramid);
1388  }
1389  /* reverse the weighting */
1390  VSCALE(*cent, *cent, (1/volume));
1391  for (i = 0; i < aip->neqn; i++) {
1392  bu_free((char *)faces[i].pts, "rt_arbn_centroid: pts");
1393  }
1394  bu_free((char *)faces, "rt_arbn_centroid: faces");
1395 }
1396 
1397 
1398 void
1399 rt_arbn_volume(fastf_t *volume, const struct rt_db_internal *ip)
1400 {
1401  struct poly_face *faces;
1402  struct rt_arbn_internal *aip = (struct rt_arbn_internal *)ip->idb_ptr;
1403  size_t i;
1404 
1405  *volume = 0.0;
1406  /* allocate array of face structs */
1407  faces = (struct poly_face *)bu_calloc(aip->neqn, sizeof(struct poly_face), "rt_arbn_volume: faces");
1408  for (i = 0; i < aip->neqn; i++) {
1409  /* allocate array of pt structs, max number of verts per faces = (# of faces) - 1 */
1410  faces[i].pts = (point_t *)bu_calloc(aip->neqn - 1, sizeof(point_t), "rt_arbn_volume: pts");
1411  }
1412  rt_arbn_faces_area(faces, aip);
1413  for (i = 0; i < aip->neqn; i++) {
1414  vect_t tmp;
1415 
1416  /* calculate volume of pyramid */
1417  VSCALE(tmp, faces[i].plane_eqn, faces[i].area);
1418  *volume += VDOT(faces[i].pts[0], tmp)/3;
1419  }
1420  for (i = 0; i < aip->neqn; i++) {
1421  bu_free((char *)faces[i].pts, "rt_arbn_volume: pts");
1422  }
1423  bu_free((char *)faces, "rt_arbn_volume: faces");
1424 }
1425 
1426 
1427 /** @} */
1428 /*
1429  * Local Variables:
1430  * mode: C
1431  * tab-width: 8
1432  * indent-tabs-mode: t
1433  * c-file-style: "stroustrup"
1434  * End:
1435  * ex: shiftwidth=4 tabstop=8
1436  */
void nmg_fix_normals(struct shell *s_orig, const struct bn_tol *tol)
Definition: nmg_misc.c:3505
int rt_arbn_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
Definition: arbn.c:1057
int rt_arbn_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: arbn.c:1004
Definition: db_flip.c:35
Definition: raytrace.h:800
#define RT_LEN_TOL
Definition: raytrace.h:169
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define BU_LIST_INSERT(old, new)
Definition: list.h:183
int rt_arbn_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: arbn.c:874
#define SIZEOF_NETWORK_DOUBLE
Definition: cv.h:48
struct hit seg_in
IN information.
Definition: raytrace.h:370
int bn_mat_is_identity(const mat_t m)
Definition: mat.c:980
Definition: list.h:118
int nmg_fu_planeeqn(struct faceuse *fu, const struct bn_tol *tol)
Definition: nmg_mod.c:1311
#define RT_DOT_TOL
Definition: raytrace.h:170
#define RT_CK_APPLICATION(_p)
Definition: raytrace.h:1675
int rt_arbn_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
Definition: arbn.c:554
double dist
>= 0
Definition: tol.h:73
vect_t crv_pdir
Principle direction.
Definition: raytrace.h:307
const mat_t bn_mat_identity
Matrix and vector functionality.
Definition: mat.c:46
fastf_t uv_u
Range 0..1.
Definition: raytrace.h:341
struct soltab * seg_stp
pointer back to soltab
Definition: raytrace.h:372
int rt_arbn_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
Definition: arbn.c:238
if lu s
Definition: nmg_mod.c:3860
int plane_no[3]
Definition: arbn.c:472
point_t cent
Definition: arbn.c:1290
void bu_vls_strcat(struct bu_vls *vp, const char *s)
Definition: vls.c:368
#define VSETALL(a, s)
Definition: color.c:54
Definition: raytrace.h:215
int nmg_kr(struct nmgregion *r)
Definition: nmg_mk.c:1595
#define SIZEOF_NETWORK_LONG
Definition: cv.h:46
Definition: pc.h:108
double dist_sq
dist * dist
Definition: tol.h:74
Definition: raytrace.h:368
#define BU_ASSERT_LONG(_lhs, _relation, _rhs)
Definition: defines.h:240
Definition: raytrace.h:248
void nmg_vertex_gv(struct vertex *v, const fastf_t *pt)
Definition: nmg_mk.c:1668
#define SMALL_FASTF
Definition: defines.h:342
fastf_t st_aradius
Radius of APPROXIMATING sphere.
Definition: raytrace.h:433
Header file for the BRL-CAD common definitions.
void rt_arbn_surf_area(fastf_t *area, const struct rt_db_internal *ip)
Definition: arbn.c:1321
void rt_arbn_centroid(point_t *cent, const struct rt_db_internal *ip)
Definition: arbn.c:1343
int bn_mkpoint_3planes(point_t pt, const plane_t a, const plane_t b, const plane_t c)
Given the description of three planes, compute the point of intersection, if any. The direction vecto...
struct resource * a_resource
dynamic memory resources
Definition: raytrace.h:1591
void bu_cv_htond(unsigned char *out, const unsigned char *in, size_t count)
point_t pt
Definition: arbn.c:471
struct bu_list l
Definition: raytrace.h:369
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
#define BN_VECT_ARE_PARALLEL(_dot, _tol)
Definition: tol.h:115
if(share_geom)
Definition: nmg_mod.c:3829
int idb_major_type
Definition: raytrace.h:192
Definition: color.c:49
int bn_polygon_mk_pts_planes(size_t *npts, point_t **pts, size_t neqs, const plane_t *eqs)
Calculate for an array of plane_eqs, which build a polyhedron, the point_t's for each face...
Definition: polygon.c:137
#define RT_ADD_VLIST(hd, pnt, draw)
Definition: raytrace.h:1865
int rt_arbn_adjust(struct bu_vls *logstr, struct rt_db_internal *intern, int argc, const char **argv)
Definition: arbn.c:1168
#define RT_CK_DB_INTERNAL(_p)
Definition: raytrace.h:207
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
fastf_t st_bradius
Radius of BOUNDING sphere.
Definition: raytrace.h:434
Definition: arbn.c:469
fastf_t crv_c2
curvature in other direction
Definition: raytrace.h:309
#define RT_CK_HIT(_p)
Definition: raytrace.h:259
point_t cent_pyramid
Definition: arbn.c:1289
int v1_no
Definition: arbn.c:479
#define BN_VLIST_LINE_MOVE
Definition: vlist.h:82
const struct rt_functab * idb_meth
for ft_ifree(), etc.
Definition: raytrace.h:194
#define LOC(i, j)
Definition: arbn.c:483
fastf_t uv_dv
delta in v
Definition: raytrace.h:344
void rt_arbn_free(struct soltab *stp)
Definition: arbn.c:367
#define BRLCAD_OK
Definition: defines.h:71
uint8_t * ext_buf
Definition: parse.h:216
void rt_arbn_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
Definition: arbn.c:310
point_t hit_point
DEPRECATED: Intersection point, use VJOIN1 hit_dist.
Definition: raytrace.h:251
point_t st_max
max X, Y, Z of bounding RPP
Definition: raytrace.h:438
#define SQRT_SMALL_FASTF
Definition: defines.h:346
struct hit seg_out
OUT information.
Definition: raytrace.h:371
#define BN_VLIST_LINE_DRAW
Definition: vlist.h:83
void * bu_realloc(void *ptr, size_t siz, const char *str)
point_t * pts
Definition: arbn.c:1285
#define UNUSED(parameter)
Definition: common.h:239
int bn_polygon_area(fastf_t *area, size_t npts, const point_t *pts)
Functions for working with polygons.
Definition: polygon.c:28
int nmg_mark_edges_real(const uint32_t *magic_p)
Definition: nmg_misc.c:846
goto out
Definition: nmg_mod.c:3846
Support for uniform tolerances.
Definition: tol.h:71
void rt_arbn_volume(fastf_t *volume, const struct rt_db_internal *ip)
Definition: arbn.c:1399
int bn_polygon_centroid(point_t *cent, size_t npts, const point_t *pts)
Calculate the centroid of a non self-intersecting polygon.
Definition: polygon.c:72
size_t npts
Definition: arbn.c:1284
struct nmgregion * nmg_mrsv(struct model *m)
Definition: nmg_mk.c:306
int rt_arbn_params(struct pc_pc_set *ps, const struct rt_db_internal *ip)
Definition: arbn.c:1268
vect_t r_dir
Direction of ray (UNIT Length)
Definition: raytrace.h:219
int rt_arbn_get(struct bu_vls *logstr, const struct rt_db_internal *intern, const char *attr)
Definition: arbn.c:1113
#define RT_ARBN_INTERNAL_MAGIC
Definition: magic.h:81
struct bn_tol rti_tol
Math tolerances for this model.
Definition: raytrace.h:1765
#define ID_ARBN
ARB with N faces.
Definition: raytrace.h:472
#define RT_CK_DBI(_p)
Definition: raytrace.h:829
#define RT_GET_SEG(p, res)
Definition: raytrace.h:379
int rt_arbn_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
Definition: arbn.c:933
char label[5]
Definition: arbn.c:1283
#define ZERO(val)
Definition: units.c:38
void * idb_ptr
Definition: raytrace.h:195
void rt_arbn_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
Definition: arbn.c:333
point_t r_pt
Point at which ray starts.
Definition: raytrace.h:218
point_t st_min
min X, Y, Z of bounding RPP
Definition: raytrace.h:437
plane_t plane_eqn
Definition: arbn.c:1286
void bu_cv_ntohd(unsigned char *out, const unsigned char *in, size_t count)
int rt_arbn_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *tol)
Definition: arbn.c:50
const struct rt_functab OBJ[]
Definition: table.c:159
void bn_vec_ortho(vect_t out, const vect_t in)
void bu_vls_printf(struct bu_vls *vls, const char *fmt,...) _BU_ATTR_PRINTF23
Definition: vls.c:694
void * st_specific
-> ID-specific (private) struct
Definition: raytrace.h:435
fastf_t uv_du
delta in u
Definition: raytrace.h:343
int tcl_list_to_fastf_array(Tcl_Interp *interp, const char *char_list, fastf_t **array, int *array_len)
Definition: tcl.c:829
fastf_t crv_c1
curvature in principle dir
Definition: raytrace.h:308
Definition: color.c:51
void bu_vls_strcpy(struct bu_vls *vp, const char *s)
Definition: vls.c:310
int hit_surfno
solid-specific surface indicator
Definition: raytrace.h:255
int rt_arbn_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: arbn.c:105
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
vect_t hit_normal
DEPRECATED: Surface Normal at hit_point, use RT_HIT_NORMAL.
Definition: raytrace.h:252
#define BU_CK_LIST_HEAD(_p)
Definition: list.h:142
int rt_arbn_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol, const struct rt_view_info *info)
Definition: arbn.c:387
#define BU_CK_EXTERNAL(_p)
Definition: parse.h:224
void rt_arbn_ifree(struct rt_db_internal *ip)
Definition: arbn.c:1087
struct faceuse * nmg_cmface(struct shell *s, struct vertex ***verts, int n)
Definition: nmg_mod.c:979
int rt_arbn_import4(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
Definition: arbn.c:801
size_t ext_nbytes
Definition: parse.h:210
void rt_arbn_print(const struct soltab *stp)
Definition: arbn.c:209
fastf_t hit_dist
dist from r_pt to hit_point
Definition: raytrace.h:250
HIDDEN void verbose(struct human_data_t *dude)
Definition: human.c:2008
Definition: vls.h:56
#define BRLCAD_ERROR
Definition: defines.h:72
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
void rt_arbn_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
Definition: arbn.c:351
fastf_t area
Definition: arbn.c:1287
double fastf_t
Definition: defines.h:300
#define VPRINT(a, b)
Definition: raytrace.h:1881
Tcl_Interp * brlcad_interp
Definition: tcl.c:41
struct vertex ** vp
Definition: arbn.c:473
Definition: color.c:50
#define BU_LIST_FIRST(structure, hp)
Definition: list.h:312
fastf_t uv_v
Range 0..1.
Definition: raytrace.h:342
fastf_t vol_pyramid
Definition: arbn.c:1288
point_t st_center
Centroid of solid.
Definition: raytrace.h:432
int bn_polygon_sort_ccw(size_t npts, point_t *pts, plane_t cmp)
Sort an array of point_ts, building a convex polygon, counter-clockwise.
Definition: polygon.c:182
#define BU_STR_EQUAL(s1, s2)
Definition: str.h:126
int v2_no
Definition: arbn.c:479
void nmg_region_a(struct nmgregion *r, const struct bn_tol *tol)
Definition: nmg_mk.c:2557