BRL-CAD
tgc.c
Go to the documentation of this file.
1 /* T G C . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1985-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/tgc/tgc.c
23  *
24  * Intersect a ray with a Truncated General Cone.
25  *
26  * Method -
27  * TGC: solve quartic equation of cone and line
28  *
29  */
30 /** @} */
31 
32 #include "common.h"
33 
34 #include <stddef.h>
35 #include <string.h>
36 #include <math.h>
37 #include "bio.h"
38 
39 #include "bu/cv.h"
40 #include "vmath.h"
41 #include "db.h"
42 #include "nmg.h"
43 #include "rtgeom.h"
44 #include "raytrace.h"
45 #include "nurb.h"
46 
47 #include "../../librt_private.h"
48 
49 
50 struct tgc_specific {
51  vect_t tgc_V; /* Vector to center of base of TGC */
52  fastf_t tgc_sH; /* magnitude of sheared H vector */
53  fastf_t tgc_A; /* magnitude of A vector */
54  fastf_t tgc_B; /* magnitude of B vector */
55  fastf_t tgc_C; /* magnitude of C vector */
56  fastf_t tgc_D; /* magnitude of D vector */
57  fastf_t tgc_CdAm1; /* (C/A - 1) */
58  fastf_t tgc_DdBm1; /* (D/B - 1) */
59  fastf_t tgc_AAdCC; /* (|A|**2)/(|C|**2) */
60  fastf_t tgc_BBdDD; /* (|B|**2)/(|D|**2) */
61  vect_t tgc_N; /* normal at 'top' of cone */
62  mat_t tgc_ScShR; /* Scale(Shear(Rot(vect))) */
63  mat_t tgc_invRtShSc; /* invRot(trnShear(Scale(vect))) */
64  char tgc_AD_CB; /* boolean: A*D == C*B */
65 };
66 
67 
68 extern int rt_rec_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip);
69 
70 static void rt_tgc_rotate(fastf_t *A, fastf_t *B, fastf_t *Hv, fastf_t *Rot, fastf_t *Inv, struct tgc_specific *tgc);
71 static void rt_tgc_shear(const fastf_t *vect, int axis, fastf_t *Shr, fastf_t *Trn, fastf_t *Inv);
72 static void rt_tgc_scale(fastf_t a, fastf_t b, fastf_t h, fastf_t *Scl, fastf_t *Inv);
73 static void nmg_tgc_disk(struct faceuse *fu, fastf_t *rmat, fastf_t height, int flip);
74 static void nmg_tgc_nurb_cyl(struct faceuse *fu, fastf_t *top_mat, fastf_t *bot_mat);
75 
76 void rt_pt_sort(register fastf_t *t, int npts);
77 
78 
79 #define VLARGE 1000000.0
80 #define MAX_RATIO 10.0 /* maximum allowed height-to-width ration for triangles */
81 #define RAT M_SQRT1_2
82 
83 #define OUT 0
84 #define IN 1
85 
86 /* hit_surfno is set to one of these */
87 #define TGC_NORM_BODY (1) /* compute normal */
88 #define TGC_NORM_TOP (2) /* copy tgc_N */
89 #define TGC_NORM_BOT (3) /* copy reverse tgc_N */
90 
91 #define RT_TGC_SEG_MISS(SEG) (SEG).seg_stp=RT_SOLTAB_NULL
92 #define VEXCHANGE(a, b, tmp) { VMOVE(tmp, a); VMOVE(a, b); VMOVE(b, tmp); }
93 #define ALPHA(x, y, c, d) ((x)*(x)*(c) + (y)*(y)*(d))
94 
95 /* determines the class of tgc given vector magnitudes a, b, c, d */
96 #define GET_TGC_TYPE(type, a, b, c, d) \
97 { \
98  if (EQUAL((a), (b)) && EQUAL((c), (d))) { \
99  if (EQUAL((a), (c))) { \
100  (type) = RCC; \
101  } else { \
102  (type) = TRC; \
103  } \
104  } else { \
105  if (EQUAL((a), (c)) && EQUAL((b), (d))) { \
106  (type) = REC; \
107  } else { \
108  (type) = TEC; \
109  } \
110  } \
111 }
112 
113 
114 const struct bu_structparse rt_tgc_parse[] = {
115  { "%f", 3, "V", bu_offsetofarray(struct rt_tgc_internal, v, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
116  { "%f", 3, "H", bu_offsetofarray(struct rt_tgc_internal, h, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
117  { "%f", 3, "A", bu_offsetofarray(struct rt_tgc_internal, a, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
118  { "%f", 3, "B", bu_offsetofarray(struct rt_tgc_internal, b, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
119  { "%f", 3, "C", bu_offsetofarray(struct rt_tgc_internal, c, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
120  { "%f", 3, "D", bu_offsetofarray(struct rt_tgc_internal, d, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
121  { {'\0', '\0', '\0', '\0'}, 0, (char *)NULL, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL }
122 };
123 
124 
125 static const fastf_t nmg_tgc_unitcircle[36] = {
126  1.0, 0.0, 0.0, 1.0,
127  RAT, -RAT, 0.0, RAT,
128  0.0, -1.0, 0.0, 1.0,
129  -RAT, -RAT, 0.0, RAT,
130  -1.0, 0.0, 0.0, 1.0,
131  -RAT, RAT, 0.0, RAT,
132  0.0, 1.0, 0.0, 1.0,
133  RAT, RAT, 0.0, RAT,
134  1.0, 0.0, 0.0, 1.0
135 };
136 
137 
138 static const fastf_t nmg_uv_unitcircle[27] = {
139  1.0, .5, 1.0,
140  RAT, RAT, RAT,
141  .5, 1.0, 1.0,
142  0.0, RAT, RAT,
143  0.0, .5, 1.0,
144  0.0, 0.0, RAT,
145  .5, 0.0, 1.0,
146  RAT, 0.0, RAT,
147  1.0, .5, 1.0
148 };
149 
150 
151 /**
152  * Compute the bounding RPP for a truncated general cone
153  */
154 int
155 rt_tgc_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol)) {
156  vect_t work, temp;
157 
158  struct rt_tgc_internal *tip = (struct rt_tgc_internal *)ip->idb_ptr;
159  RT_TGC_CK_MAGIC(tip);
160 
161  VSETALL((*min), INFINITY);
162  VSETALL((*max), -INFINITY);
163 
164  /* There are 8 corners to the bounding RPP */
165  /* This may not be minimal, but does fully contain the TGC */
166  VADD2(temp, tip->v, tip->a);
167  VADD2(work, temp, tip->b);
168  VMINMAX((*min), (*max), work); /* V + A + B */
169  VSUB2(work, temp, tip->b);
170  VMINMAX((*min), (*max), work); /* V + A - B */
171 
172  VSUB2(temp, tip->v, tip->a);
173  VADD2(work, temp, tip->b);
174  VMINMAX((*min), (*max), work); /* V - A + B */
175  VSUB2(work, temp, tip->b);
176  VMINMAX((*min), (*max), work); /* V - A - B */
177 
178  VADD3(temp, tip->v, tip->h, tip->c);
179  VADD2(work, temp, tip->d);
180  VMINMAX((*min), (*max), work); /* V + H + C + D */
181  VSUB2(work, temp, tip->d);
182  VMINMAX((*min), (*max), work); /* V + H + C - D */
183 
184  VADD2(temp, tip->v, tip->h);
185  VSUB2(temp, temp, tip->c);
186  VADD2(work, temp, tip->d);
187  VMINMAX((*min), (*max), work); /* V + H - C + D */
188  VSUB2(work, temp, tip->d);
189  VMINMAX((*min), (*max), work); /* V + H - C - D */
190  return 0;
191 }
192 
193 
194 /**
195  * Given the parameters (in vector form) of a truncated general cone,
196  * compute the constant terms and a transformation matrix needed for
197  * solving the intersection of a ray with the cone.
198  *
199  * Also compute the return transformation for normals in the
200  * transformed space to the original space. This NOT the inverse of
201  * the transformation matrix (if you really want to know why, talk to
202  * Ed Davisson).
203  */
204 int
205 rt_tgc_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
206 {
207  struct rt_tgc_internal *tip;
208  register struct tgc_specific *tgc;
209  register fastf_t f;
210  fastf_t prod_ab, prod_cd;
211  fastf_t magsq_a, magsq_b, magsq_c, magsq_d;
212  fastf_t mag_h, mag_a, mag_b, mag_c, mag_d;
213  mat_t Rot, Shr, Scl;
214  mat_t iRot, tShr, iShr, iScl;
215  mat_t tmp;
216  vect_t nH;
217  vect_t work;
218 
219  tip = (struct rt_tgc_internal *)ip->idb_ptr;
220  RT_TGC_CK_MAGIC(tip);
221 
222  /*
223  * For a fast way out, hand this solid off to the REC routine.
224  * If it takes it, then there is nothing to do, otherwise
225  * the solid is a TGC.
226  */
227  if (rt_rec_prep(stp, ip, rtip) == 0)
228  return 0; /* OK */
229 
230  /* Validate that |H| > 0, compute |A| |B| |C| |D| */
231  mag_h = sqrt(MAGSQ(tip->h));
232  mag_a = sqrt(magsq_a = MAGSQ(tip->a));
233  mag_b = sqrt(magsq_b = MAGSQ(tip->b));
234  mag_c = sqrt(magsq_c = MAGSQ(tip->c));
235  mag_d = sqrt(magsq_d = MAGSQ(tip->d));
236  prod_ab = mag_a * mag_b;
237  prod_cd = mag_c * mag_d;
238 
239  if (NEAR_ZERO(mag_h, RT_LEN_TOL)) {
240  bu_log("tgc(%s): zero length H vector\n", stp->st_name);
241  return 1; /* BAD */
242  }
243 
244  /* Validate that figure is not two-dimensional */
245  if (NEAR_ZERO(mag_a, RT_LEN_TOL) &&
246  NEAR_ZERO(mag_c, RT_LEN_TOL)) {
247  bu_log("tgc(%s): vectors A, C zero length\n", stp->st_name);
248  return 1;
249  }
250  if (NEAR_ZERO(mag_b, RT_LEN_TOL) &&
251  NEAR_ZERO(mag_d, RT_LEN_TOL)) {
252  bu_log("tgc(%s): vectors B, D zero length\n", stp->st_name);
253  return 1;
254  }
255 
256  /* Validate that both ends are not degenerate */
257  if (prod_ab <= SMALL) {
258  /* AB end is degenerate */
259  if (prod_cd <= SMALL) {
260  bu_log("tgc(%s): Both ends degenerate\n", stp->st_name);
261  return 1; /* BAD */
262  }
263  /* Exchange ends, so that in solids with one degenerate end,
264  * the CD end always is always the degenerate one
265  */
266  VADD2(tip->v, tip->v, tip->h);
267  VREVERSE(tip->h, tip->h);
268  VEXCHANGE(tip->a, tip->c, work);
269  VEXCHANGE(tip->b, tip->d, work);
270  bu_log("NOTE: tgc(%s): degenerate end exchanged\n", stp->st_name);
271  }
272 
273  /* Ascertain whether H lies in A-B plane */
274  VCROSS(work, tip->a, tip->b);
275  f = VDOT(tip->h, work) / (prod_ab*mag_h);
276  if (NEAR_ZERO(f, RT_DOT_TOL)) {
277  bu_log("tgc(%s): H lies in A-B plane\n", stp->st_name);
278  return 1; /* BAD */
279  }
280 
281  if (prod_ab > SMALL) {
282  /* Validate that A.B == 0 */
283  f = VDOT(tip->a, tip->b) / prod_ab;
284  if (! NEAR_ZERO(f, RT_DOT_TOL)) {
285  bu_log("tgc(%s): A not perpendicular to B, f=%g\n",
286  stp->st_name, f);
287  bu_log("tgc: dot=%g / a*b=%g\n",
288  VDOT(tip->a, tip->b), prod_ab);
289  return 1; /* BAD */
290  }
291  }
292  if (prod_cd > SMALL) {
293  /* Validate that C.D == 0 */
294  f = VDOT(tip->c, tip->d) / prod_cd;
295  if (! NEAR_ZERO(f, RT_DOT_TOL)) {
296  bu_log("tgc(%s): C not perpendicular to D, f=%g\n",
297  stp->st_name, f);
298  bu_log("tgc: dot=%g / c*d=%g\n",
299  VDOT(tip->c, tip->d), prod_cd);
300  return 1; /* BAD */
301  }
302  }
303 
304  if (mag_a * mag_c > SMALL) {
305  /* Validate that A || C */
306  f = 1.0 - VDOT(tip->a, tip->c) / (mag_a * mag_c);
307  if (! NEAR_ZERO(f, RT_DOT_TOL)) {
308  bu_log("tgc(%s): A not parallel to C, f=%g\n",
309  stp->st_name, f);
310  return 1; /* BAD */
311  }
312  }
313 
314  if (mag_b * mag_d > SMALL) {
315  /* Validate that B || D, for parallel planes */
316  f = 1.0 - VDOT(tip->b, tip->d) / (mag_b * mag_d);
317  if (! NEAR_ZERO(f, RT_DOT_TOL)) {
318  bu_log("tgc(%s): B not parallel to D, f=%g\n",
319  stp->st_name, f);
320  return 1; /* BAD */
321  }
322  }
323 
324  /* solid is OK, compute constant terms, etc. */
325  BU_GET(tgc, struct tgc_specific);
326  stp->st_specific = (void *)tgc;
327 
328  VMOVE(tgc->tgc_V, tip->v);
329  tgc->tgc_A = mag_a;
330  tgc->tgc_B = mag_b;
331  tgc->tgc_C = mag_c;
332  tgc->tgc_D = mag_d;
333 
335  bu_log("%s: a is %.20f, b is %.20f, c is %.20f, d is %.20f\n",
336  stp->st_name, magsq_a, magsq_b, magsq_c, magsq_d);
337 
338  /* Part of computing ALPHA() */
339  if (ZERO(magsq_c)) {
340  tgc->tgc_AAdCC = VLARGE;
341  } else {
342  tgc->tgc_AAdCC = magsq_a / magsq_c;
343  }
344  if (ZERO(magsq_d)) {
345  tgc->tgc_BBdDD = VLARGE;
346  } else {
347  tgc->tgc_BBdDD = magsq_b / magsq_d;
348  }
349 
350  /* If the eccentricities of the two ellipses are the same, then
351  * the cone equation reduces to a much simpler quadratic form.
352  * Otherwise it is a (gah!) quartic equation.
353  */
354  tgc->tgc_AD_CB = NEAR_EQUAL((tgc->tgc_A*tgc->tgc_D), (tgc->tgc_C*tgc->tgc_B), 0.0001); /* A*D == C*B */
355  rt_tgc_rotate(tip->a, tip->b, tip->h, Rot, iRot, tgc);
356  MAT4X3VEC(nH, Rot, tip->h);
357  tgc->tgc_sH = nH[Z];
358 
359  tgc->tgc_CdAm1 = tgc->tgc_C/tgc->tgc_A - 1.0;
360  tgc->tgc_DdBm1 = tgc->tgc_D/tgc->tgc_B - 1.0;
361  if (ZERO(tgc->tgc_CdAm1))
362  tgc->tgc_CdAm1 = 0.0;
363  if (ZERO(tgc->tgc_DdBm1))
364  tgc->tgc_DdBm1 = 0.0;
365 
366  /*
367  * Added iShr parameter to tgc_shear(). Changed inverse
368  * transformation of normal vectors of std. solid intersection to
369  * include shear inverse (tgc_invRtShSc). Fold in scaling
370  * transformation into the transformation to std. space from
371  * target space (tgc_ScShR).
372  */
373  rt_tgc_shear(nH, Z, Shr, tShr, iShr);
374  rt_tgc_scale(tgc->tgc_A, tgc->tgc_B, tgc->tgc_sH, Scl, iScl);
375 
376  bn_mat_mul(tmp, Shr, Rot);
377  bn_mat_mul(tgc->tgc_ScShR, Scl, tmp);
378 
379  bn_mat_mul(tmp, tShr, Scl);
380  bn_mat_mul(tgc->tgc_invRtShSc, iRot, tmp);
381 
382  /* Compute bounding sphere and RPP */
383  {
384  fastf_t dx, dy, dz; /* For bounding sphere */
385  if (stp->st_meth->ft_bbox(ip, &(stp->st_min), &(stp->st_max), &(rtip->rti_tol))) return 1;
386  VSET(stp->st_center,
387  (stp->st_max[X] + stp->st_min[X])/2,
388  (stp->st_max[Y] + stp->st_min[Y])/2,
389  (stp->st_max[Z] + stp->st_min[Z])/2);
390 
391  dx = (stp->st_max[X] - stp->st_min[X])/2;
392  f = dx;
393  dy = (stp->st_max[Y] - stp->st_min[Y])/2;
394  if (dy > f) f = dy;
395  dz = (stp->st_max[Z] - stp->st_min[Z])/2;
396  if (dz > f) f = dz;
397  stp->st_aradius = f;
398  stp->st_bradius = sqrt(dx*dx + dy*dy + dz*dz);
399  }
400  return 0;
401 }
402 
403 
404 /**
405  * To rotate vectors A and B (where A is perpendicular to B) to the X
406  * and Y axes respectively, create a rotation matrix
407  *
408  * | A' |
409  * R = | B' |
410  * | C' |
411  *
412  * where A', B' and C' are vectors such that
413  *
414  * A' = A/|A| B' = B/|B| C' = C/|C|
415  *
416  * where C = H - (H.A')A' - (H.B')B'
417  *
418  * The last operation (Gram Schmidt method) finds the component of the
419  * vector H perpendicular A and to B. This is, therefore the normal
420  * for the planar sections of the truncated cone.
421  */
422 static void
423 rt_tgc_rotate(fastf_t *A, fastf_t *B, fastf_t *Hv, fastf_t *Rot, fastf_t *Inv, struct tgc_specific *tgc)
424 {
425  vect_t uA, uB, uC; /* unit vectors */
426  fastf_t mag_ha; /* magnitude of H in the */
427  fastf_t mag_hb; /* A and B directions */
428 
429  /* copy A and B, then 'unitize' the results */
430  VMOVE(uA, A);
431  VUNITIZE(uA);
432  VMOVE(uB, B);
433  VUNITIZE(uB);
434 
435  /* Find component of H in the A direction */
436  mag_ha = VDOT(Hv, uA);
437  /* Find component of H in the B direction */
438  mag_hb = VDOT(Hv, uB);
439 
440  /* Subtract the A and B components of H to find the component
441  * perpendicular to both, then 'unitize' the result.
442  */
443  VJOIN2(uC, Hv, -mag_ha, uA, -mag_hb, uB);
444  VUNITIZE(uC);
445  VMOVE(tgc->tgc_N, uC);
446 
447  MAT_IDN(Rot);
448  MAT_IDN(Inv);
449 
450  Rot[0] = Inv[0] = uA[X];
451  Rot[1] = Inv[4] = uA[Y];
452  Rot[2] = Inv[8] = uA[Z];
453 
454  Rot[4] = Inv[1] = uB[X];
455  Rot[5] = Inv[5] = uB[Y];
456  Rot[6] = Inv[9] = uB[Z];
457 
458  Rot[8] = Inv[2] = uC[X];
459  Rot[9] = Inv[6] = uC[Y];
460  Rot[10] = Inv[10] = uC[Z];
461 }
462 
463 
464 /**
465  * To shear the H vector to the Z axis, every point must be shifted in
466  * the X direction by -(Hx/Hz)*z, and in the Y direction by -(Hy/Hz)*z
467  * This operation makes the equation for the standard cone much easier
468  * to work with.
469  *
470  * NOTE: This computes the TRANSPOSE of the shear matrix rather than
471  * the inverse.
472  *
473  * Begin changes GSM, EOD -- Added INVERSE (Inv) calculation.
474  */
475 static void
476 rt_tgc_shear(const fastf_t *vect, int axis, fastf_t *Shr, fastf_t *Trn, fastf_t *Inv)
477 {
478  MAT_IDN(Shr);
479  MAT_IDN(Trn);
480  MAT_IDN(Inv);
481 
482  if (ZERO(vect[axis]))
483  bu_bomb("rt_tgc_shear() divide by zero\n");
484 
485  if (axis == X) {
486  Inv[4] = -(Shr[4] = Trn[1] = -vect[Y]/vect[X]);
487  Inv[8] = -(Shr[8] = Trn[2] = -vect[Z]/vect[X]);
488  } else if (axis == Y) {
489  Inv[1] = -(Shr[1] = Trn[4] = -vect[X]/vect[Y]);
490  Inv[9] = -(Shr[9] = Trn[6] = -vect[Z]/vect[Y]);
491  } else if (axis == Z) {
492  Inv[2] = -(Shr[2] = Trn[8] = -vect[X]/vect[Z]);
493  Inv[6] = -(Shr[6] = Trn[9] = -vect[Y]/vect[Z]);
494  }
495 }
496 
497 
498 static void
499 rt_tgc_scale(fastf_t a, fastf_t b, fastf_t h, fastf_t *Scl, fastf_t *Inv)
500 {
501  MAT_IDN(Scl);
502  MAT_IDN(Inv);
503  Scl[0] /= a;
504  Scl[5] /= b;
505  Scl[10] /= h;
506  Inv[0] = a;
507  Inv[5] = b;
508  Inv[10] = h;
509  return;
510 }
511 
512 
513 void
514 rt_tgc_print(register const struct soltab *stp)
515 {
516  register const struct tgc_specific *tgc =
517  (struct tgc_specific *)stp->st_specific;
518 
519  VPRINT("V", tgc->tgc_V);
520  bu_log("mag sheared H = %f\n", tgc->tgc_sH);
521  bu_log("mag A = %f\n", tgc->tgc_A);
522  bu_log("mag B = %f\n", tgc->tgc_B);
523  bu_log("mag C = %f\n", tgc->tgc_C);
524  bu_log("mag D = %f\n", tgc->tgc_D);
525  VPRINT("Top normal", tgc->tgc_N);
526 
527  bn_mat_print("Sc o Sh o R", tgc->tgc_ScShR);
528  bn_mat_print("invR o trnSh o Sc", tgc->tgc_invRtShSc);
529 
530  if (tgc->tgc_AD_CB) {
531  bu_log("A*D == C*B. Equal eccentricities gives quadratic equation.\n");
532  } else {
533  bu_log("A*D != C*B. Quartic equation.\n");
534  }
535  bu_log("(C/A - 1) = %f\n", tgc->tgc_CdAm1);
536  bu_log("(D/B - 1) = %f\n", tgc->tgc_DdBm1);
537  bu_log("(|A|**2)/(|C|**2) = %f\n", tgc->tgc_AAdCC);
538  bu_log("(|B|**2)/(|D|**2) = %f\n", tgc->tgc_BBdDD);
539 }
540 
541 
542 /**
543  * Intersect a ray with a truncated general cone, where all constant
544  * terms have been computed by rt_tgc_prep().
545  *
546  * NOTE: All lines in this function are represented parametrically by
547  * a point, P(Px, Py, Pz) and a unit direction vector, D = iDx + jDy +
548  * kDz. Any point on a line can be expressed by one variable 't',
549  * where
550  *
551  * X = Dx*t + Px,
552  * Y = Dy*t + Py,
553  * Z = Dz*t + Pz.
554  *
555  * First, convert the line to the coordinate system of a "stan- dard"
556  * cone. This is a cone whose base lies in the X-Y plane, and whose H
557  * (now H') vector is lined up with the Z axis.
558  *
559  * Then find the equation of that line and the standard cone as an
560  * equation in 't'. Solve the equation using a general polynomial
561  * root finder. Use those values of 't' to compute the points of
562  * intersection in the original coordinate system.
563  */
564 int
565 rt_tgc_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
566 {
567  register const struct tgc_specific *tgc =
568  (struct tgc_specific *)stp->st_specific;
569  register struct seg *segp;
570  vect_t pprime;
571  vect_t dprime;
572  vect_t work;
573  fastf_t k[6];
574  int hit_type[6];
575  fastf_t t, b, zval, dir;
576  fastf_t t_scale;
577  fastf_t alf1, alf2;
578  int npts;
579  int intersect;
580  vect_t cor_pprime; /* corrected P prime */
581  fastf_t cor_proj = 0; /* corrected projected dist */
582  int i;
583  bn_poly_t C; /* final equation */
584  bn_poly_t Xsqr, Ysqr;
585  bn_poly_t R, Rsqr;
586 
587  /* find rotated point and direction */
588  MAT4X3VEC(dprime, tgc->tgc_ScShR, rp->r_dir);
589 
590  /* A vector of unit length in model space (r_dir) changes length
591  * in the special unit-tgc space. This scale factor will restore
592  * proper length after hit points are found.
593  */
594  t_scale = MAGNITUDE(dprime);
595  if (ZERO(t_scale)) {
596  bu_log("tgc(%s) dprime=(%g, %g, %g), t_scale=%e, miss.\n", stp->st_dp->d_namep,
597  V3ARGS(dprime), t_scale);
598  return 0;
599  }
600  t_scale = 1/t_scale;
601  VSCALE(dprime, dprime, t_scale); /* VUNITIZE(dprime); */
602 
603  if (NEAR_ZERO(dprime[Z], RT_PCOEF_TOL)) {
604  dprime[Z] = 0.0; /* prevent rootfinder heartburn */
605  }
606 
607  VSUB2(work, rp->r_pt, tgc->tgc_V);
608  MAT4X3VEC(pprime, tgc->tgc_ScShR, work);
609 
610  /* Translating ray origin along direction of ray to closest pt. to
611  * origin of solids coordinate system, new ray origin is
612  * 'cor_pprime'.
613  */
614  cor_proj = -VDOT(pprime, dprime);
615  VJOIN1(cor_pprime, pprime, cor_proj, dprime);
616 
617  /* The TGC is defined in "unit" space, so the parametric distance
618  * from one side of the TGC to the other is on the order of 2.
619  * Therefore, any vector/point coordinates that are very small
620  * here may be considered to be zero, since double precision only
621  * has 18 digits of significance. If these tiny values were left
622  * in, then as they get squared (below) they will cause
623  * difficulties.
624  */
625  for (i=0; i<3; i++) {
626  /* Direction cosines */
627  if (NEAR_ZERO(dprime[i], 1e-10)) {
628  dprime[i] = 0;
629  }
630  /* Position in -1..+1 coordinates */
631  if (NEAR_ZERO(cor_pprime[i], 1e-20)) {
632  cor_pprime[i] = 0;
633  }
634  }
635 
636  /* Given a line and the parameters for a standard cone, finds the
637  * roots of the equation for that cone and line. Returns the
638  * number of real roots found.
639  *
640  * Given a line and the cone parameters, finds the equation of the
641  * cone in terms of the variable 't'.
642  *
643  * The equation for the cone is:
644  *
645  * X**2 * Q**2 + Y**2 * R**2 - R**2 * Q**2 = 0
646  *
647  * where R = a + ((c - a)/|H'|)*Z
648  * Q = b + ((d - b)/|H'|)*Z
649  *
650  * First, find X, Y, and Z in terms of 't' for this line, then
651  * substitute them into the equation above.
652  *
653  * Express each variable (X, Y, and Z) as a linear equation in
654  * 'k', e.g., (dprime[X] * k) + cor_pprime[X], and substitute into
655  * the cone equation.
656  */
657  Xsqr.dgr = 2;
658  Xsqr.cf[0] = dprime[X] * dprime[X];
659  Xsqr.cf[1] = 2.0 * dprime[X] * cor_pprime[X];
660  Xsqr.cf[2] = cor_pprime[X] * cor_pprime[X];
661 
662  Ysqr.dgr = 2;
663  Ysqr.cf[0] = dprime[Y] * dprime[Y];
664  Ysqr.cf[1] = 2.0 * dprime[Y] * cor_pprime[Y];
665  Ysqr.cf[2] = cor_pprime[Y] * cor_pprime[Y];
666 
667  R.dgr = 1;
668  R.cf[0] = dprime[Z] * tgc->tgc_CdAm1;
669  /* A vector is unitized (tgc->tgc_A == 1.0) */
670  R.cf[1] = (cor_pprime[Z] * tgc->tgc_CdAm1) + 1.0;
671 
672  /* (void) rt_poly_mul(&Rsqr, &R, &R); */
673  Rsqr.dgr = 2;
674  Rsqr.cf[0] = R.cf[0] * R.cf[0];
675  Rsqr.cf[1] = R.cf[0] * R.cf[1] * 2;
676  Rsqr.cf[2] = R.cf[1] * R.cf[1];
677 
678  /* If the eccentricities of the two ellipses are the same, then
679  * the cone equation reduces to a much simpler quadratic form.
680  * Otherwise it is a (gah!) quartic equation.
681  *
682  * this can only be done when C.cf[0] is not too small! (JRA)
683  */
684  C.cf[0] = Xsqr.cf[0] + Ysqr.cf[0] - Rsqr.cf[0];
685  if (tgc->tgc_AD_CB && !NEAR_ZERO(C.cf[0], 1.0e-10)) {
686  fastf_t roots;
687 
688  /*
689  * (void) bn_poly_add(&sum, &Xsqr, &Ysqr);
690  * (void) bn_poly_sub(&C, &sum, &Rsqr);
691  */
692  C.dgr = 2;
693  C.cf[1] = Xsqr.cf[1] + Ysqr.cf[1] - Rsqr.cf[1];
694  C.cf[2] = Xsqr.cf[2] + Ysqr.cf[2] - Rsqr.cf[2];
695 
696  /* Find the real roots the easy way. C.dgr==2 */
697  if ((roots = C.cf[1]*C.cf[1] - 4 * C.cf[0] * C.cf[2]) < 0) {
698  npts = 0; /* no real roots */
699  } else {
700  register fastf_t f;
701  roots = sqrt(roots);
702  k[0] = (roots - C.cf[1]) * (f = 0.5 / C.cf[0]);
703  hit_type[0] = TGC_NORM_BODY;
704  k[1] = (roots + C.cf[1]) * -f;
705  hit_type[1] = TGC_NORM_BODY;
706  npts = 2;
707  }
708  } else {
709  bn_poly_t Q, Qsqr;
710  bn_complex_t val[4]; /* roots of final equation */
711  register int l;
712  register int nroots;
713 
714  Q.dgr = 1;
715  Q.cf[0] = dprime[Z] * tgc->tgc_DdBm1;
716  /* B vector is unitized (tgc->tgc_B == 1.0) */
717  Q.cf[1] = (cor_pprime[Z] * tgc->tgc_DdBm1) + 1.0;
718 
719  /* (void) bn_poly_mul(&Qsqr, &Q, &Q); */
720  Qsqr.dgr = 2;
721  Qsqr.cf[0] = Q.cf[0] * Q.cf[0];
722  Qsqr.cf[1] = Q.cf[0] * Q.cf[1] * 2;
723  Qsqr.cf[2] = Q.cf[1] * Q.cf[1];
724 
725  /*
726  * (void) bn_poly_mul(&T1, &Qsqr, &Xsqr);
727  * (void) bn_poly_mul(&T2 &Rsqr, &Ysqr);
728  * (void) bn_poly_mul(&T1, &Rsqr, &Qsqr);
729  * (void) bn_poly_add(&sum, &T1, &T2);
730  * (void) bn_poly_sub(&C, &sum, &T3);
731  */
732  C.dgr = 4;
733  C.cf[0] = Qsqr.cf[0] * Xsqr.cf[0] +
734  Rsqr.cf[0] * Ysqr.cf[0] -
735  (Rsqr.cf[0] * Qsqr.cf[0]);
736  C.cf[1] = Qsqr.cf[0] * Xsqr.cf[1] + Qsqr.cf[1] * Xsqr.cf[0] +
737  Rsqr.cf[0] * Ysqr.cf[1] + Rsqr.cf[1] * Ysqr.cf[0] -
738  (Rsqr.cf[0] * Qsqr.cf[1] + Rsqr.cf[1] * Qsqr.cf[0]);
739  C.cf[2] = Qsqr.cf[0] * Xsqr.cf[2] + Qsqr.cf[1] * Xsqr.cf[1] +
740  Qsqr.cf[2] * Xsqr.cf[0] +
741  Rsqr.cf[0] * Ysqr.cf[2] + Rsqr.cf[1] * Ysqr.cf[1] +
742  Rsqr.cf[2] * Ysqr.cf[0] -
743  (Rsqr.cf[0] * Qsqr.cf[2] + Rsqr.cf[1] * Qsqr.cf[1] +
744  Rsqr.cf[2] * Qsqr.cf[0]);
745  C.cf[3] = Qsqr.cf[1] * Xsqr.cf[2] + Qsqr.cf[2] * Xsqr.cf[1] +
746  Rsqr.cf[1] * Ysqr.cf[2] + Rsqr.cf[2] * Ysqr.cf[1] -
747  (Rsqr.cf[1] * Qsqr.cf[2] + Rsqr.cf[2] * Qsqr.cf[1]);
748  C.cf[4] = Qsqr.cf[2] * Xsqr.cf[2] +
749  Rsqr.cf[2] * Ysqr.cf[2] -
750  (Rsqr.cf[2] * Qsqr.cf[2]);
751 
752  /* The equation is 4th order, so we expect 0 to 4 roots */
753  nroots = rt_poly_roots(&C, val, stp->st_dp->d_namep);
754 
755  /* bn_pr_roots("roots", val, nroots); */
756 
757  /* Only real roots indicate an intersection in real space.
758  *
759  * Look at each root returned; if the imaginary part is zero
760  * or sufficiently close, then use the real part as one value
761  * of 't' for the intersections
762  */
763  for (l=0, npts=0; l < nroots; l++) {
764  if (NEAR_ZERO(val[l].im, 1e-2)) {
765  hit_type[npts] = TGC_NORM_BODY;
766  k[npts++] = val[l].re;
767  }
768  }
769  /* bu_log("npts rooted is %d; ", npts); */
770 
771  /* Here, 'npts' is number of points being returned */
772  if (npts != 0 && npts != 2 && npts != 4 && npts > 0) {
773  bu_log("tgc: reduced %d to %d roots\n", nroots, npts);
774  bn_pr_roots(stp->st_name, val, nroots);
775  } else if (nroots < 0) {
776  static int reported=0;
777  bu_log("The root solver failed to converge on a solution for %s\n", stp->st_dp->d_namep);
778  if (!reported) {
779  VPRINT("while shooting from:\t", rp->r_pt);
780  VPRINT("while shooting at:\t", rp->r_dir);
781  bu_log("Additional truncated general cone convergence failure details will be suppressed.\n");
782  reported=1;
783  }
784  }
785  }
786 
787  /*
788  * Reverse above translation by adding distance to all 'k' values.
789  */
790  for (i = 0; i < npts; ++i) {
791  k[i] += cor_proj;
792  }
793 
794  /* bu_log("npts before elimination is %d; ", npts); */
795  /*
796  * Eliminate hits beyond the end planes
797  */
798  i = 0;
799  while (i < npts) {
800  zval = k[i]*dprime[Z] + pprime[Z];
801  /* Height vector is unitized (tgc->tgc_sH == 1.0) */
802  if (zval >= 1.0 || zval <= 0.0) {
803  int j;
804  /* drop this hit */
805  npts--;
806  for (j=i; j<npts; j++) {
807  hit_type[j] = hit_type[j+1];
808  k[j] = k[j+1];
809  }
810  } else {
811  i++;
812  }
813  }
814 
815  /*
816  * Consider intersections with the end ellipses
817  */
818  /* bu_log("npts before base is %d; ", npts); */
819  dir = VDOT(tgc->tgc_N, rp->r_dir);
820  if (!ZERO(dprime[Z]) && !NEAR_ZERO(dir, RT_DOT_TOL)) {
821  b = (-pprime[Z])/dprime[Z];
822  /* Height vector is unitized (tgc->tgc_sH == 1.0) */
823  t = (1.0 - pprime[Z])/dprime[Z];
824 
825  /* the top end */
826  VJOIN1(work, pprime, b, dprime);
827  /* A and B vectors are unitized (tgc->tgc_A == _B == 1.0) */
828  /* alf1 = ALPHA(work[X], work[Y], 1.0, 1.0) */
829  alf1 = work[X]*work[X] + work[Y]*work[Y];
830 
831  /* the bottom end */
832  VJOIN1(work, pprime, t, dprime);
833  /* Must scale C and D vectors */
834  alf2 = ALPHA(work[X], work[Y], tgc->tgc_AAdCC, tgc->tgc_BBdDD);
835 
836  /*
837  bu_log("alf1 is %f, alf2 is %f\n", alf1, alf2);
838  bu_log("work[x]=%f, work[y]=%f, aadcc=%f, bbddd=%f\n", work[X], work[Y], tgc->tgc_AAdCC, tgc->tgc_BBdDD);
839  */
840  if (alf1 <= 1.0) {
841  hit_type[npts] = TGC_NORM_BOT;
842  k[npts++] = b;
843  }
844  if (alf2 <= 1.0) {
845  hit_type[npts] = TGC_NORM_TOP;
846  k[npts++] = t;
847  }
848  }
849 
850  /* bu_log("npts FINAL is %d\n", npts); */
851 
852  /* Sort Most distant to least distant: rt_pt_sort(k, npts) */
853  {
854  register fastf_t u;
855  register short lim, n;
856  register int type;
857 
858  for (lim = npts-1; lim > 0; lim--) {
859  for (n = 0; n < lim; n++) {
860  if ((u=k[n]) < k[n+1]) {
861  /* bubble larger towards [0] */
862  type = hit_type[n];
863  hit_type[n] = hit_type[n+1];
864  hit_type[n+1] = type;
865  k[n] = k[n+1];
866  k[n+1] = u;
867  }
868  }
869  }
870  }
871  /* Now, k[0] > k[npts-1] */
872 
873  if (npts%2) {
874  /* odd number of hits!!!
875  * perhaps we got two hits on an edge
876  * check for duplicate hit distances
877  */
878 
879  for (i=npts-1; i>0; i--) {
880  fastf_t diff;
881 
882  diff = k[i-1] - k[i]; /* non-negative due to sorting */
883  if (diff < ap->a_rt_i->rti_tol.dist) {
884  /* remove this duplicate hit */
885  int j;
886 
887  npts--;
888  for (j=i; j<npts; j++) {
889  hit_type[j] = hit_type[j+1];
890  k[j] = k[j+1];
891  }
892 
893  /* now have even number of hits */
894  break;
895  }
896  }
897  }
898 
899  if (npts != 0 && npts != 2 && npts != 4) {
900  static size_t tgc_msgs = 0;
901  /* these are printed in 'mm' regardless of local units */
902 
903  if (tgc_msgs++ < 100) {
904  bu_log("tgc(%s): %d intersects != {0, 2, 4}\n", stp->st_name, npts);
905  bu_log("\tray: pt = (%g %g %g), dir = (%g %g %g), units in mm\n", V3ARGS(ap->a_ray.r_pt), V3ARGS(ap->a_ray.r_dir));
906  for (i = 0; i < npts; i++) {
907  bu_log("\t%g", k[i]*t_scale);
908  }
909  bu_log("\n");
910  } else if (tgc_msgs == 100) {
911  bu_log("tgc(%s): too many grazing intersections encountered. further reporting suppressed.\n",
912  stp->st_name);
913  tgc_msgs++;
914  }
915 
916  return 0; /* No hit */
917  }
918 
919  intersect = 0;
920  for (i=npts-1; i>0; i -= 2) {
921  RT_GET_SEG(segp, ap->a_resource);
922  segp->seg_stp = stp;
923 
924  segp->seg_in.hit_dist = k[i] * t_scale;
925  segp->seg_in.hit_surfno = hit_type[i];
926  if (segp->seg_in.hit_surfno == TGC_NORM_BODY) {
927  VJOIN1(segp->seg_in.hit_vpriv, pprime, k[i], dprime);
928  } else {
929  if (dir > 0.0) {
930  segp->seg_in.hit_surfno = TGC_NORM_BOT;
931  } else {
932  segp->seg_in.hit_surfno = TGC_NORM_TOP;
933  }
934  }
935 
936  segp->seg_out.hit_dist = k[i-1] * t_scale;
937  segp->seg_out.hit_surfno = hit_type[i-1];
938  if (segp->seg_out.hit_surfno == TGC_NORM_BODY) {
939  VJOIN1(segp->seg_out.hit_vpriv, pprime, k[i-1], dprime);
940  } else {
941  if (dir > 0.0) {
942  segp->seg_out.hit_surfno = TGC_NORM_TOP;
943  } else {
944  segp->seg_out.hit_surfno = TGC_NORM_BOT;
945  }
946  }
947  intersect++;
948  BU_LIST_INSERT(&(seghead->l), &(segp->l));
949  }
950 
951  return intersect;
952 }
953 
954 
955 /**
956  * The Homer vectorized version.
957  */
958 void
959 rt_tgc_vshot(struct soltab **stp, register struct xray **rp, struct seg *segp, int n, struct application *ap)
960 
961 
962 /* array of segs (results returned) */
963 /* Number of ray/object pairs */
964 
965 {
966  register struct tgc_specific *tgc;
967  register int ix;
968  vect_t pprime;
969  vect_t dprime;
970  vect_t work;
971  fastf_t k[4], pt[2];
972  fastf_t t, b, zval, dir;
973  fastf_t t_scale = 0;
974  fastf_t alf1, alf2;
975  int npts;
976  int intersect;
977  vect_t cor_pprime; /* corrected P prime */
978  fastf_t cor_proj = 0; /* corrected projected dist */
979  int i;
980  bn_poly_t *C; /* final equation */
981  bn_poly_t Xsqr, Ysqr;
982  bn_poly_t R, Rsqr;
983 
984  VSETALLN(k, 0, 4);
985  VSETALLN(pt, 0, 2);
986 
987  if (ap) RT_CK_APPLICATION(ap);
988 
989  /* Allocate space for polys and roots */
990  C = (bn_poly_t *)bu_malloc(n * sizeof(bn_poly_t), "tor bn_poly_t");
991 
992  /* Initialize seg_stp to assume hit (zero will then flag miss) */
993  for (ix = 0; ix < n; ix++) segp[ix].seg_stp = stp[ix];
994 
995  /* for each ray/cone pair */
996  for (ix = 0; ix < n; ix++) {
997  if (segp[ix].seg_stp == 0) continue; /* == 0 signals skip ray */
998 
999  tgc = (struct tgc_specific *)stp[ix]->st_specific;
1000 
1001  /* find rotated point and direction */
1002  MAT4X3VEC(dprime, tgc->tgc_ScShR, rp[ix]->r_dir);
1003 
1004  /*
1005  * A vector of unit length in model space (r_dir) changes length in
1006  * the special unit-tgc space. This scale factor will restore
1007  * proper length after hit points are found.
1008  */
1009  t_scale = 1/MAGNITUDE(dprime);
1010  VSCALE(dprime, dprime, t_scale); /* VUNITIZE(dprime); */
1011 
1012  if (NEAR_ZERO(dprime[Z], RT_PCOEF_TOL))
1013  dprime[Z] = 0.0; /* prevent rootfinder heartburn */
1014 
1015  /* Use segp[ix].seg_in.hit_normal as tmp to hold dprime */
1016  VMOVE(segp[ix].seg_in.hit_normal, dprime);
1017 
1018  VSUB2(work, rp[ix]->r_pt, tgc->tgc_V);
1019  MAT4X3VEC(pprime, tgc->tgc_ScShR, work);
1020 
1021  /* Use segp[ix].seg_out.hit_normal as tmp to hold pprime */
1022  VMOVE(segp[ix].seg_out.hit_normal, pprime);
1023 
1024  /* Translating ray origin along direction of ray to closest
1025  * pt. to origin of solids coordinate system, new ray origin
1026  * is 'cor_pprime'.
1027  */
1028  cor_proj = VDOT(pprime, dprime);
1029  VSCALE(cor_pprime, dprime, cor_proj);
1030  VSUB2(cor_pprime, pprime, cor_pprime);
1031 
1032  /*
1033  * Given a line and the parameters for a standard cone, finds
1034  * the roots of the equation for that cone and line. Returns
1035  * the number of real roots found.
1036  *
1037  * Given a line and the cone parameters, finds the equation of
1038  * the cone in terms of the variable 't'.
1039  *
1040  * The equation for the cone is:
1041  *
1042  * X**2 * Q**2 + Y**2 * R**2 - R**2 * Q**2 = 0
1043  *
1044  * where:
1045  * R = a + ((c - a)/|H'|)*Z
1046  * Q = b + ((d - b)/|H'|)*Z
1047  *
1048  * First, find X, Y, and Z in terms of 't' for this line, then
1049  * substitute them into the equation above.
1050  *
1051  * Express each variable (X, Y, and Z) as a linear equation in
1052  * 'k', e.g., (dprime[X] * k) + cor_pprime[X], and substitute
1053  * into the cone equation.
1054  */
1055  Xsqr.dgr = 2;
1056  Xsqr.cf[0] = dprime[X] * dprime[X];
1057  Xsqr.cf[1] = 2.0 * dprime[X] * cor_pprime[X];
1058  Xsqr.cf[2] = cor_pprime[X] * cor_pprime[X];
1059 
1060  Ysqr.dgr = 2;
1061  Ysqr.cf[0] = dprime[Y] * dprime[Y];
1062  Ysqr.cf[1] = 2.0 * dprime[Y] * cor_pprime[Y];
1063  Ysqr.cf[2] = cor_pprime[Y] * cor_pprime[Y];
1064 
1065  R.dgr = 1;
1066  R.cf[0] = dprime[Z] * tgc->tgc_CdAm1;
1067  /* A vector is unitized (tgc->tgc_A == 1.0) */
1068  R.cf[1] = (cor_pprime[Z] * tgc->tgc_CdAm1) + 1.0;
1069 
1070  /* (void) bn_poly_mul(&Rsqr, &R, &R); manual expansion: */
1071  Rsqr.dgr = 2;
1072  Rsqr.cf[0] = R.cf[0] * R.cf[0];
1073  Rsqr.cf[1] = R.cf[0] * R.cf[1] * 2;
1074  Rsqr.cf[2] = R.cf[1] * R.cf[1];
1075 
1076  /*
1077  * If the eccentricities of the two ellipses are the same,
1078  * then the cone equation reduces to a much simpler quadratic
1079  * form. Otherwise it is a (gah!) quartic equation.
1080  */
1081  if (tgc->tgc_AD_CB) {
1082  /* (void) bn_poly_add(&sum, &Xsqr, &Ysqr); and */
1083  /* (void) bn_poly_sub(&C, &sum, &Rsqr); manual expansion: */
1084  C[ix].dgr = 2;
1085  C[ix].cf[0] = Xsqr.cf[0] + Ysqr.cf[0] - Rsqr.cf[0];
1086  C[ix].cf[1] = Xsqr.cf[1] + Ysqr.cf[1] - Rsqr.cf[1];
1087  C[ix].cf[2] = Xsqr.cf[2] + Ysqr.cf[2] - Rsqr.cf[2];
1088  } else {
1089  bn_poly_t Q, Qsqr;
1090 
1091  Q.dgr = 1;
1092  Q.cf[0] = dprime[Z] * tgc->tgc_DdBm1;
1093  /* B vector is unitized (tgc->tgc_B == 1.0) */
1094  Q.cf[1] = (cor_pprime[Z] * tgc->tgc_DdBm1) + 1.0;
1095 
1096  /* (void) bn_poly_mul(&Qsqr, &Q, &Q); manual expansion: */
1097  Qsqr.dgr = 2;
1098  Qsqr.cf[0] = Q.cf[0] * Q.cf[0];
1099  Qsqr.cf[1] = Q.cf[0] * Q.cf[1] * 2;
1100  Qsqr.cf[2] = Q.cf[1] * Q.cf[1];
1101 
1102  /* (void) bn_poly_mul(&T1, &Qsqr, &Xsqr); manual expansion: */
1103  C[ix].dgr = 4;
1104  C[ix].cf[0] = Qsqr.cf[0] * Xsqr.cf[0];
1105  C[ix].cf[1] = Qsqr.cf[0] * Xsqr.cf[1] +
1106  Qsqr.cf[1] * Xsqr.cf[0];
1107  C[ix].cf[2] = Qsqr.cf[0] * Xsqr.cf[2] +
1108  Qsqr.cf[1] * Xsqr.cf[1] +
1109  Qsqr.cf[2] * Xsqr.cf[0];
1110  C[ix].cf[3] = Qsqr.cf[1] * Xsqr.cf[2] +
1111  Qsqr.cf[2] * Xsqr.cf[1];
1112  C[ix].cf[4] = Qsqr.cf[2] * Xsqr.cf[2];
1113 
1114  /* (void) bn_poly_mul(&T2, &Rsqr, &Ysqr); and */
1115  /* (void) bn_poly_add(&sum, &T1, &T2); manual expansion: */
1116  C[ix].cf[0] += Rsqr.cf[0] * Ysqr.cf[0];
1117  C[ix].cf[1] += Rsqr.cf[0] * Ysqr.cf[1] +
1118  Rsqr.cf[1] * Ysqr.cf[0];
1119  C[ix].cf[2] += Rsqr.cf[0] * Ysqr.cf[2] +
1120  Rsqr.cf[1] * Ysqr.cf[1] +
1121  Rsqr.cf[2] * Ysqr.cf[0];
1122  C[ix].cf[3] += Rsqr.cf[1] * Ysqr.cf[2] +
1123  Rsqr.cf[2] * Ysqr.cf[1];
1124  C[ix].cf[4] += Rsqr.cf[2] * Ysqr.cf[2];
1125 
1126  /* (void) bn_poly_mul(&T3, &Rsqr, &Qsqr); and */
1127  /* (void) bn_poly_sub(&C, &sum, &T3); manual expansion: */
1128  C[ix].cf[0] -= Rsqr.cf[0] * Qsqr.cf[0];
1129  C[ix].cf[1] -= Rsqr.cf[0] * Qsqr.cf[1] +
1130  Rsqr.cf[1] * Qsqr.cf[0];
1131  C[ix].cf[2] -= Rsqr.cf[0] * Qsqr.cf[2] +
1132  Rsqr.cf[1] * Qsqr.cf[1] +
1133  Rsqr.cf[2] * Qsqr.cf[0];
1134  C[ix].cf[3] -= Rsqr.cf[1] * Qsqr.cf[2] +
1135  Rsqr.cf[2] * Qsqr.cf[1];
1136  C[ix].cf[4] -= Rsqr.cf[2] * Qsqr.cf[2];
1137  }
1138 
1139  }
1140 
1141  /* It seems impractical to try to vectorize finding and sorting roots. */
1142  for (ix = 0; ix < n; ix++) {
1143  if (segp[ix].seg_stp == 0) continue; /* == 0 signals skip ray */
1144 
1145  /* Again, check for the equal eccentricities case. */
1146  if (C[ix].dgr == 2) {
1147  fastf_t roots;
1148 
1149  /* Find the real roots the easy way. */
1150  if ((roots = C[ix].cf[1]*C[ix].cf[1]-4*C[ix].cf[0]*C[ix].cf[2]
1151  ) < 0) {
1152  npts = 0; /* no real roots */
1153  } else {
1154  roots = sqrt(roots);
1155  k[0] = (roots - C[ix].cf[1]) * 0.5 / C[ix].cf[0];
1156  k[1] = (roots + C[ix].cf[1]) * (-0.5) / C[ix].cf[0];
1157  npts = 2;
1158  }
1159  } else {
1160  bn_complex_t val[4]; /* roots of final equation */
1161  register int l;
1162  register int nroots;
1163 
1164  /* The equation is 4th order, so we expect 0 to 4 roots */
1165  nroots = rt_poly_roots(&C[ix], val, (*stp)->st_dp->d_namep);
1166 
1167  /* Only real roots indicate an intersection in real space.
1168  *
1169  * Look at each root returned; if the imaginary part is zero
1170  * or sufficiently close, then use the real part as one value
1171  * of 't' for the intersections
1172  */
1173  for (l=0, npts=0; l < nroots; l++) {
1174  if (NEAR_ZERO(val[l].im, 0.0001))
1175  k[npts++] = val[l].re;
1176  }
1177  /* Here, 'npts' is number of points being returned */
1178  if (npts != 0 && npts != 2 && npts != 4 && npts > 0) {
1179  bu_log("tgc: reduced %d to %d roots\n", nroots, npts);
1180  bn_pr_roots("tgc", val, nroots);
1181  } else if (nroots < 0) {
1182  static int reported=0;
1183  bu_log("The root solver failed to converge on a solution for %s\n", stp[ix]->st_dp->d_namep);
1184  if (!reported) {
1185  VPRINT("while shooting from:\t", rp[ix]->r_pt);
1186  VPRINT("while shooting at:\t", rp[ix]->r_dir);
1187  bu_log("Additional truncated general cone convergence failure details will be suppressed.\n");
1188  reported=1;
1189  }
1190  }
1191  }
1192 
1193  /*
1194  * Reverse above translation by adding distance to all 'k' values.
1195  */
1196  for (i = 0; i < npts; ++i)
1197  k[i] -= cor_proj;
1198 
1199  if (npts != 0 && npts != 2 && npts != 4) {
1200  bu_log("tgc(%s): %d intersects != {0, 2, 4}\n",
1201  stp[ix]->st_name, npts);
1202  RT_TGC_SEG_MISS(segp[ix]); /* No hit */
1203  continue;
1204  }
1205 
1206  /* Most distant to least distant */
1207  rt_pt_sort(k, npts);
1208 
1209  /* Now, k[0] > k[npts-1] */
1210 
1211  /* General Cone may have 4 intersections, but Truncated Cone
1212  * may only have 2.
1213  */
1214 
1215  /* Truncation Procedure
1216  *
1217  * Determine whether any of the intersections found are
1218  * between the planes truncating the cone.
1219  */
1220  intersect = 0;
1221  tgc = (struct tgc_specific *)stp[ix]->st_specific;
1222  for (i=0; i < npts; i++) {
1223  /* segp[ix].seg_in.hit_normal holds dprime */
1224  /* segp[ix].seg_out.hit_normal holds pprime */
1225  zval = k[i]*segp[ix].seg_in.hit_normal[Z] +
1226  segp[ix].seg_out.hit_normal[Z];
1227  /* Height vector is unitized (tgc->tgc_sH == 1.0) */
1228  if (zval < 1.0 && zval > 0.0) {
1229  if (++intersect == 2) {
1230  pt[IN] = k[i];
1231  } else
1232  pt[OUT] = k[i];
1233  }
1234  }
1235  /* Reuse C to hold values of intersect and k. */
1236  C[ix].dgr = intersect;
1237  C[ix].cf[OUT] = pt[OUT];
1238  C[ix].cf[IN] = pt[IN];
1239  }
1240 
1241  /* for each ray/cone pair */
1242  for (ix = 0; ix < n; ix++) {
1243  if (segp[ix].seg_stp == 0) continue; /* Skip */
1244 
1245  tgc = (struct tgc_specific *)stp[ix]->st_specific;
1246  intersect = C[ix].dgr;
1247  pt[OUT] = C[ix].cf[OUT];
1248  pt[IN] = C[ix].cf[IN];
1249  /* segp[ix].seg_out.hit_normal holds pprime */
1250  VMOVE(pprime, segp[ix].seg_out.hit_normal);
1251  /* segp[ix].seg_in.hit_normal holds dprime */
1252  VMOVE(dprime, segp[ix].seg_in.hit_normal);
1253 
1254  if (intersect == 2) {
1255  /* If two between-plane intersections exist, they are
1256  * the hit points for the ray.
1257  */
1258  segp[ix].seg_in.hit_dist = pt[IN] * t_scale;
1259  segp[ix].seg_in.hit_surfno = TGC_NORM_BODY; /* compute N */
1260  VJOIN1(segp[ix].seg_in.hit_vpriv, pprime, pt[IN], dprime);
1261 
1262  segp[ix].seg_out.hit_dist = pt[OUT] * t_scale;
1263  segp[ix].seg_out.hit_surfno = TGC_NORM_BODY; /* compute N */
1264  VJOIN1(segp[ix].seg_out.hit_vpriv, pprime, pt[OUT], dprime);
1265  } else if (intersect == 1) {
1266  int nflag;
1267  /*
1268  * If only one between-plane intersection exists (pt[OUT]),
1269  * then the other intersection must be on
1270  * one of the planar surfaces (pt[IN]).
1271  *
1272  * Find which surface it lies on by calculating the
1273  * X and Y values of the line as it intersects each
1274  * plane (in the standard coordinate system), and test
1275  * whether this lies within the governing ellipse.
1276  */
1277  if (ZERO(dprime[Z])) {
1278  RT_TGC_SEG_MISS(segp[ix]);
1279  continue;
1280  }
1281  b = (-pprime[Z])/dprime[Z];
1282  /* Height vector is unitized (tgc->tgc_sH == 1.0) */
1283  t = (1.0 - pprime[Z])/dprime[Z];
1284 
1285  VJOIN1(work, pprime, b, dprime);
1286  /* A and B vectors are unitized (tgc->tgc_A == _B == 1.0) */
1287  /* alf1 = ALPHA(work[X], work[Y], 1.0, 1.0) */
1288  alf1 = work[X]*work[X] + work[Y]*work[Y];
1289 
1290  VJOIN1(work, pprime, t, dprime);
1291  /* Must scale C and D vectors */
1292  alf2 = ALPHA(work[X], work[Y], tgc->tgc_AAdCC, tgc->tgc_BBdDD);
1293 
1294  if (alf1 <= 1.0) {
1295  pt[IN] = b;
1296  nflag = TGC_NORM_BOT; /* copy reverse normal */
1297  } else if (alf2 <= 1.0) {
1298  pt[IN] = t;
1299  nflag = TGC_NORM_TOP; /* copy normal */
1300  } else {
1301  /* intersection apparently invalid */
1302  RT_TGC_SEG_MISS(segp[ix]);
1303  continue;
1304  }
1305 
1306  /* pt[OUT] on skin, pt[IN] on end */
1307  if (pt[OUT] >= pt[IN]) {
1308  segp[ix].seg_in.hit_dist = pt[IN] * t_scale;
1309  segp[ix].seg_in.hit_surfno = nflag;
1310 
1311  segp[ix].seg_out.hit_dist = pt[OUT] * t_scale;
1312  segp[ix].seg_out.hit_surfno = TGC_NORM_BODY; /* compute N */
1313  /* transform-space vector needed for normal */
1314  VJOIN1(segp[ix].seg_out.hit_vpriv, pprime, pt[OUT], dprime);
1315  } else {
1316  segp[ix].seg_in.hit_dist = pt[OUT] * t_scale;
1317  /* transform-space vector needed for normal */
1318  segp[ix].seg_in.hit_surfno = TGC_NORM_BODY; /* compute N */
1319  VJOIN1(segp[ix].seg_in.hit_vpriv, pprime, pt[OUT], dprime);
1320 
1321  segp[ix].seg_out.hit_dist = pt[IN] * t_scale;
1322  segp[ix].seg_out.hit_surfno = nflag;
1323  }
1324  } else {
1325 
1326  /* If all conic intersections lie outside the plane, then
1327  * check to see whether there are two planar intersections
1328  * inside the governing ellipses.
1329  *
1330  * But first, if the direction is parallel (or nearly so)
1331  * to the planes, it (obviously) won't intersect either of
1332  * them.
1333  */
1334  if (ZERO(dprime[Z])) {
1335  RT_TGC_SEG_MISS(segp[ix]);
1336  continue;
1337  }
1338 
1339  dir = VDOT(tgc->tgc_N, rp[ix]->r_dir); /* direc */
1340  if (NEAR_ZERO(dir, RT_DOT_TOL)) {
1341  RT_TGC_SEG_MISS(segp[ix]);
1342  continue;
1343  }
1344 
1345  b = (-pprime[Z])/dprime[Z];
1346  /* Height vector is unitized (tgc->tgc_sH == 1.0) */
1347  t = (1.0 - pprime[Z])/dprime[Z];
1348 
1349  VJOIN1(work, pprime, b, dprime);
1350  /* A and B vectors are unitized (tgc->tgc_A == _B == 1.0) */
1351  /* alpf = ALPHA(work[0], work[1], 1.0, 1.0) */
1352  alf1 = work[X]*work[X] + work[Y]*work[Y];
1353 
1354  VJOIN1(work, pprime, t, dprime);
1355  /* Must scale C and D vectors. */
1356  alf2 = ALPHA(work[X], work[Y], tgc->tgc_AAdCC, tgc->tgc_BBdDD);
1357 
1358  /* It should not be possible for one planar intersection
1359  * to be outside its ellipse while the other is inside ...
1360  * but I wouldn't take any chances.
1361  */
1362  if (alf1 > 1.0 || alf2 > 1.0) {
1363  RT_TGC_SEG_MISS(segp[ix]);
1364  continue;
1365  }
1366 
1367  /* Use the dot product (found earlier) of the plane normal
1368  * with the direction vector to determine the orientation
1369  * of the intersections.
1370  */
1371  if (dir > 0.0) {
1372  segp[ix].seg_in.hit_dist = b * t_scale;
1373  segp[ix].seg_in.hit_surfno = TGC_NORM_BOT; /* reverse normal */
1374 
1375  segp[ix].seg_out.hit_dist = t * t_scale;
1376  segp[ix].seg_out.hit_surfno = TGC_NORM_TOP; /* normal */
1377  } else {
1378  segp[ix].seg_in.hit_dist = t * t_scale;
1379  segp[ix].seg_in.hit_surfno = TGC_NORM_TOP; /* normal */
1380 
1381  segp[ix].seg_out.hit_dist = b * t_scale;
1382  segp[ix].seg_out.hit_surfno = TGC_NORM_BOT; /* reverse normal */
1383  }
1384  }
1385  } /* end for each ray/cone pair */
1386  bu_free((char *)C, "tor bn_poly_t");
1387 }
1388 
1389 
1390 /**
1391  * Sorts the values in t[] in descending order.
1392  */
1393 void
1394 rt_pt_sort(fastf_t t[], int npts)
1395 {
1396  fastf_t u;
1397  short lim, n;
1398 
1399  for (lim = npts-1; lim > 0; lim--) {
1400  for (n = 0; n < lim; n++) {
1401  if ((u=t[n]) < t[n+1]) {
1402  /* bubble larger towards [0] */
1403  t[n] = t[n+1];
1404  t[n+1] = u;
1405  }
1406  }
1407  }
1408 }
1409 
1410 
1411 /**
1412  * Compute the normal to the cone, given a point on the STANDARD CONE
1413  * centered at the origin of the X-Y plane.
1414  *
1415  * The gradient of the cone at that point is the normal vector in the
1416  * standard space. This vector will need to be transformed back to
1417  * the coordinate system of the original cone in order to be useful.
1418  * Then the transformed vector must be 'unitized.'
1419  *
1420  * NOTE: The transformation required is NOT the inverse of the of the
1421  * rotation to the standard cone, due to the shear involved in the
1422  * mapping. The inverse maps points back to points, but it is the
1423  * transpose which maps normals back to normals. If you really want
1424  * to know why, talk to Ed Davisson or Peter Stiller.
1425  *
1426  * The equation for the standard cone *without* scaling is:
1427  * (rotated the sheared)
1428  *
1429  * f(X, Y, Z) = X**2 * Q**2 + Y**2 * R**2 - R**2 * Q**2 = 0
1430  *
1431  * where:
1432  * R = a + ((c - a)/|H'|)*Z
1433  * Q = b + ((d - b)/|H'|)*Z
1434  *
1435  * When the equation is scaled so the A, B, and the sheared H are unit
1436  * length, as is done here, the equation can be coerced back into this
1437  * same form with R and Q now being:
1438  *
1439  * R = 1 + (c/a - 1)*Z
1440  * Q = 1 + (d/b - 1)*Z
1441  *
1442  * The gradient of f(x, y, z) = 0 is:
1443  *
1444  * df/dx = 2 * x * Q**2
1445  * df/dy = 2 * y * R**2
1446  * df/dz = x**2 * 2 * Q * dQ/dz + y**2 * 2 * R * dR/dz
1447  * - R**2 * 2 * Q * dQ/dz - Q**2 * 2 * R * dR/dz
1448  * = 2 [(x**2 - R**2) * Q * dQ/dz + (y**2 - Q**2) * R * dR/dz]
1449  *
1450  * where:
1451  * dR/dz = (c/a - 1)
1452  * dQ/dz = (d/b - 1)
1453  *
1454  * [in the *unscaled* case these would be:
1455  * (c - a)/|H'| and (d - b)/|H'|]
1456  *
1457  * Since the gradient (normal) needs to be rescaled to unit length
1458  * after mapping back to absolute coordinates, we divide the 2 out of
1459  * the above expressions.
1460  */
1461 void
1462 rt_tgc_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
1463 {
1464  register struct tgc_specific *tgc =
1465  (struct tgc_specific *)stp->st_specific;
1466 
1467  fastf_t Q;
1468  fastf_t R;
1469  vect_t stdnorm;
1470 
1471  /* Hit point */
1472  VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
1473 
1474  /* Hits on the end plates are easy */
1475  switch (hitp->hit_surfno) {
1476  case TGC_NORM_TOP:
1477  VMOVE(hitp->hit_normal, tgc->tgc_N);
1478  break;
1479  case TGC_NORM_BOT:
1480  VREVERSE(hitp->hit_normal, tgc->tgc_N);
1481  break;
1482  case TGC_NORM_BODY:
1483  /* Compute normal, given hit point on standard (unit) cone */
1484  R = 1 + tgc->tgc_CdAm1 * hitp->hit_vpriv[Z];
1485  Q = 1 + tgc->tgc_DdBm1 * hitp->hit_vpriv[Z];
1486  stdnorm[X] = hitp->hit_vpriv[X] * Q * Q;
1487  stdnorm[Y] = hitp->hit_vpriv[Y] * R * R;
1488  stdnorm[Z] = (hitp->hit_vpriv[X]*hitp->hit_vpriv[X] - R*R)
1489  * Q * tgc->tgc_DdBm1
1490  + (hitp->hit_vpriv[Y]*hitp->hit_vpriv[Y] - Q*Q)
1491  * R * tgc->tgc_CdAm1;
1492  MAT4X3VEC(hitp->hit_normal, tgc->tgc_invRtShSc, stdnorm);
1493  /*XXX - save scale */
1494  VUNITIZE(hitp->hit_normal);
1495  break;
1496  default:
1497  bu_log("rt_tgc_norm: bad surfno=%d\n", hitp->hit_surfno);
1498  break;
1499  }
1500 }
1501 
1502 
1503 void
1504 rt_tgc_uv(struct application *ap, struct soltab *stp, register struct hit *hitp, register struct uvcoord *uvp)
1505 {
1506  struct tgc_specific *tgc = (struct tgc_specific *)stp->st_specific;
1507 
1508  vect_t work;
1509  vect_t pprime;
1510  fastf_t len;
1511 
1512  if (ap) RT_CK_APPLICATION(ap);
1513 
1514  /* hit_point is on surface; project back to unit cylinder,
1515  * creating a vector from vertex to hit point.
1516  */
1517  VSUB2(work, hitp->hit_point, tgc->tgc_V);
1518  MAT4X3VEC(pprime, tgc->tgc_ScShR, work);
1519 
1520  switch (hitp->hit_surfno) {
1521  case TGC_NORM_BODY:
1522  /* scale coords to unit circle (they are already scaled by bottom plate radii) */
1523  pprime[X] *= tgc->tgc_A / (tgc->tgc_A*(1.0 - pprime[Z]) + tgc->tgc_C*pprime[Z]);
1524  pprime[Y] *= tgc->tgc_B / (tgc->tgc_B*(1.0 - pprime[Z]) + tgc->tgc_D*pprime[Z]);
1525  uvp->uv_u = atan2(pprime[Y], pprime[X]) / M_2PI + 0.5;
1526  uvp->uv_v = pprime[Z]; /* height */
1527  break;
1528  case TGC_NORM_TOP:
1529  /* top plate */
1530  /* scale coords to unit circle (they are already scaled by bottom plate radii) */
1531  pprime[X] *= tgc->tgc_A / tgc->tgc_C;
1532  pprime[Y] *= tgc->tgc_B / tgc->tgc_D;
1533  uvp->uv_u = atan2(pprime[Y], pprime[X]) / M_2PI + 0.5;
1534  len = sqrt(pprime[X]*pprime[X]+pprime[Y]*pprime[Y]);
1535  uvp->uv_v = len; /* rim v = 1 */
1536  break;
1537  case TGC_NORM_BOT:
1538  /* bottom plate */
1539  len = sqrt(pprime[X]*pprime[X]+pprime[Y]*pprime[Y]);
1540  uvp->uv_u = atan2(pprime[Y], pprime[X]) / M_2PI + 0.5;
1541  uvp->uv_v = 1 - len; /* rim v = 0 */
1542  break;
1543  }
1544 
1545  if (uvp->uv_u < 0.0)
1546  uvp->uv_u = 0.0;
1547  else if (uvp->uv_u > 1.0)
1548  uvp->uv_u = 1.0;
1549  if (uvp->uv_v < 0.0)
1550  uvp->uv_v = 0.0;
1551  else if (uvp->uv_v > 1.0)
1552  uvp->uv_v = 1.0;
1553 
1554  /* uv_du should be relative to rotation, uv_dv relative to height */
1555  uvp->uv_du = uvp->uv_dv = 0;
1556 }
1557 
1558 
1559 void
1560 rt_tgc_free(struct soltab *stp)
1561 {
1562  register struct tgc_specific *tgc =
1563  (struct tgc_specific *)stp->st_specific;
1564 
1565  BU_PUT(tgc, struct tgc_specific);
1566 }
1567 
1568 
1569 /**
1570  * Import a TGC from the database format to the internal format.
1571  * Apply modeling transformations as well.
1572  */
1573 int
1574 rt_tgc_import4(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
1575 {
1576  struct rt_tgc_internal *tip;
1577  union record *rp;
1578  fastf_t vec[3*6];
1579 
1580  if (dbip) RT_CK_DBI(dbip);
1581 
1582  BU_CK_EXTERNAL(ep);
1583  rp = (union record *)ep->ext_buf;
1584  /* Check record type */
1585  if (rp->u_id != ID_SOLID) {
1586  bu_log("rt_tgc_import4: defective record\n");
1587  return -1;
1588  }
1589 
1590  RT_CK_DB_INTERNAL(ip);
1591  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1592  ip->idb_type = ID_TGC;
1593  ip->idb_meth = &OBJ[ID_TGC];
1594  BU_ALLOC(ip->idb_ptr, struct rt_tgc_internal);
1595 
1596  tip = (struct rt_tgc_internal *)ip->idb_ptr;
1597  tip->magic = RT_TGC_INTERNAL_MAGIC;
1598 
1599  /* Convert from database to internal format */
1600  flip_fastf_float(vec, rp->s.s_values, 6, dbip->dbi_version < 0 ? 1 : 0);
1601 
1602  /* Apply modeling transformations */
1603  if (mat == NULL) mat = bn_mat_identity;
1604  MAT4X3PNT(tip->v, mat, &vec[0*3]);
1605  MAT4X3VEC(tip->h, mat, &vec[1*3]);
1606  MAT4X3VEC(tip->a, mat, &vec[2*3]);
1607  MAT4X3VEC(tip->b, mat, &vec[3*3]);
1608  MAT4X3VEC(tip->c, mat, &vec[4*3]);
1609  MAT4X3VEC(tip->d, mat, &vec[5*3]);
1610 
1611  return 0; /* OK */
1612 }
1613 
1614 
1615 int
1616 rt_tgc_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1617 {
1618  struct rt_tgc_internal *tip;
1619  union record *rec;
1620 
1621  if (dbip) RT_CK_DBI(dbip);
1622 
1623  RT_CK_DB_INTERNAL(ip);
1624  if (ip->idb_type != ID_TGC && ip->idb_type != ID_REC) return -1;
1625  tip = (struct rt_tgc_internal *)ip->idb_ptr;
1626  RT_TGC_CK_MAGIC(tip);
1627 
1628  BU_CK_EXTERNAL(ep);
1629  ep->ext_nbytes = sizeof(union record);
1630  ep->ext_buf = (uint8_t *)bu_calloc(1, ep->ext_nbytes, "tgc external");
1631  rec = (union record *)ep->ext_buf;
1632 
1633  rec->s.s_id = ID_SOLID;
1634  rec->s.s_type = GENTGC;
1635 
1636  /* NOTE: This also converts to dbfloat_t */
1637  VSCALE(&rec->s.s_values[0], tip->v, local2mm);
1638  VSCALE(&rec->s.s_values[3], tip->h, local2mm);
1639  VSCALE(&rec->s.s_values[6], tip->a, local2mm);
1640  VSCALE(&rec->s.s_values[9], tip->b, local2mm);
1641  VSCALE(&rec->s.s_values[12], tip->c, local2mm);
1642  VSCALE(&rec->s.s_values[15], tip->d, local2mm);
1643 
1644  return 0;
1645 }
1646 
1647 
1648 /**
1649  * Import a TGC from the database format to the internal format.
1650  * Apply modeling transformations as well.
1651  */
1652 int
1653 rt_tgc_import5(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
1654 {
1655  struct rt_tgc_internal *tip;
1656 
1657  /* must be double for import and export */
1658  double vec[ELEMENTS_PER_VECT*6];
1659 
1660  if (dbip) RT_CK_DBI(dbip);
1661 
1662  BU_CK_EXTERNAL(ep);
1663  BU_ASSERT_LONG(ep->ext_nbytes, ==, SIZEOF_NETWORK_DOUBLE * ELEMENTS_PER_VECT*6);
1664 
1665  RT_CK_DB_INTERNAL(ip);
1666  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1667  ip->idb_type = ID_TGC;
1668  ip->idb_meth = &OBJ[ID_TGC];
1669  BU_ALLOC(ip->idb_ptr, struct rt_tgc_internal);
1670 
1671  tip = (struct rt_tgc_internal *)ip->idb_ptr;
1672  tip->magic = RT_TGC_INTERNAL_MAGIC;
1673 
1674  /* Convert from database (network) to internal (host) format */
1675  bu_cv_ntohd((unsigned char *)vec, ep->ext_buf, ELEMENTS_PER_VECT*6);
1676 
1677  /* Apply modeling transformations */
1678  if (mat == NULL) mat = bn_mat_identity;
1679  MAT4X3PNT(tip->v, mat, &vec[0*ELEMENTS_PER_VECT]);
1680  MAT4X3VEC(tip->h, mat, &vec[1*ELEMENTS_PER_VECT]);
1681  MAT4X3VEC(tip->a, mat, &vec[2*ELEMENTS_PER_VECT]);
1682  MAT4X3VEC(tip->b, mat, &vec[3*ELEMENTS_PER_VECT]);
1683  MAT4X3VEC(tip->c, mat, &vec[4*ELEMENTS_PER_VECT]);
1684  MAT4X3VEC(tip->d, mat, &vec[5*ELEMENTS_PER_VECT]);
1685 
1686  return 0; /* OK */
1687 }
1688 
1689 
1690 int
1691 rt_tgc_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1692 {
1693  struct rt_tgc_internal *tip;
1694 
1695  /* must be double for import and export */
1696  double vec[ELEMENTS_PER_VECT*6];
1697 
1698  if (dbip) RT_CK_DBI(dbip);
1699 
1700  RT_CK_DB_INTERNAL(ip);
1701  if (ip->idb_type != ID_TGC && ip->idb_type != ID_REC) return -1;
1702  tip = (struct rt_tgc_internal *)ip->idb_ptr;
1703  RT_TGC_CK_MAGIC(tip);
1704 
1705  BU_CK_EXTERNAL(ep);
1706  ep->ext_nbytes = SIZEOF_NETWORK_DOUBLE * ELEMENTS_PER_VECT*6;
1707  ep->ext_buf = (uint8_t *)bu_malloc(ep->ext_nbytes, "tgc external");
1708 
1709  /* scale 'em into local buffer */
1710  VSCALE(&vec[0*ELEMENTS_PER_VECT], tip->v, local2mm);
1711  VSCALE(&vec[1*ELEMENTS_PER_VECT], tip->h, local2mm);
1712  VSCALE(&vec[2*ELEMENTS_PER_VECT], tip->a, local2mm);
1713  VSCALE(&vec[3*ELEMENTS_PER_VECT], tip->b, local2mm);
1714  VSCALE(&vec[4*ELEMENTS_PER_VECT], tip->c, local2mm);
1715  VSCALE(&vec[5*ELEMENTS_PER_VECT], tip->d, local2mm);
1716 
1717  /* Convert from internal (host) to database (network) format */
1718  bu_cv_htond(ep->ext_buf, (unsigned char *)vec, ELEMENTS_PER_VECT*6);
1719 
1720  return 0;
1721 }
1722 
1723 
1724 /**
1725  * Make human-readable formatted presentation of this solid.
1726  * First line describes type of solid.
1727  * Additional lines are indented one tab, and give parameter values.
1728  */
1729 int
1730 rt_tgc_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
1731 {
1732  struct rt_tgc_internal *tip = (struct rt_tgc_internal *)ip->idb_ptr;
1733  char buf[256];
1734  double angles[5];
1735  vect_t unitv;
1736  fastf_t Hmag;
1737 
1738  RT_TGC_CK_MAGIC(tip);
1739  bu_vls_strcat(str, "truncated general cone (TGC)\n");
1740 
1741  if (!verbose)
1742  return 0;
1743 
1744  sprintf(buf, "\tV (%g, %g, %g)\n",
1745  INTCLAMP(tip->v[X] * mm2local),
1746  INTCLAMP(tip->v[Y] * mm2local),
1747  INTCLAMP(tip->v[Z] * mm2local));
1748  bu_vls_strcat(str, buf);
1749 
1750  sprintf(buf, "\tTop (%g, %g, %g)\n",
1751  INTCLAMP((tip->v[X] + tip->h[X]) * mm2local),
1752  INTCLAMP((tip->v[Y] + tip->h[Y]) * mm2local),
1753  INTCLAMP((tip->v[Z] + tip->h[Z]) * mm2local));
1754  bu_vls_strcat(str, buf);
1755 
1756  Hmag = MAGNITUDE(tip->h);
1757  sprintf(buf, "\tH (%g, %g, %g) mag=%g\n",
1758  INTCLAMP(tip->h[X] * mm2local),
1759  INTCLAMP(tip->h[Y] * mm2local),
1760  INTCLAMP(tip->h[Z] * mm2local),
1761  INTCLAMP(Hmag * mm2local));
1762  bu_vls_strcat(str, buf);
1763  if (Hmag < VDIVIDE_TOL) {
1764  bu_vls_strcat(str, "H vector is zero!\n");
1765  } else {
1766  register double f = 1/Hmag;
1767  VSCALE(unitv, tip->h, f);
1768  rt_find_fallback_angle(angles, unitv);
1769  rt_pr_fallback_angle(str, "\tH", angles);
1770  }
1771 
1772  sprintf(buf, "\tA (%g, %g, %g) mag=%g\n",
1773  INTCLAMP(tip->a[X] * mm2local),
1774  INTCLAMP(tip->a[Y] * mm2local),
1775  INTCLAMP(tip->a[Z] * mm2local),
1776  INTCLAMP(MAGNITUDE(tip->a) * mm2local));
1777  bu_vls_strcat(str, buf);
1778 
1779  sprintf(buf, "\tB (%g, %g, %g) mag=%g\n",
1780  INTCLAMP(tip->b[X] * mm2local),
1781  INTCLAMP(tip->b[Y] * mm2local),
1782  INTCLAMP(tip->b[Z] * mm2local),
1783  INTCLAMP(MAGNITUDE(tip->b) * mm2local));
1784  bu_vls_strcat(str, buf);
1785 
1786  sprintf(buf, "\tC (%g, %g, %g) mag=%g\n",
1787  INTCLAMP(tip->c[X] * mm2local),
1788  INTCLAMP(tip->c[Y] * mm2local),
1789  INTCLAMP(tip->c[Z] * mm2local),
1790  INTCLAMP(MAGNITUDE(tip->c) * mm2local));
1791  bu_vls_strcat(str, buf);
1792 
1793  sprintf(buf, "\tD (%g, %g, %g) mag=%g\n",
1794  INTCLAMP(tip->d[X] * mm2local),
1795  INTCLAMP(tip->d[Y] * mm2local),
1796  INTCLAMP(tip->d[Z] * mm2local),
1797  INTCLAMP(MAGNITUDE(tip->d) * mm2local));
1798  bu_vls_strcat(str, buf);
1799 
1800  VCROSS(unitv, tip->c, tip->d);
1801  VUNITIZE(unitv);
1802  rt_find_fallback_angle(angles, unitv);
1803  rt_pr_fallback_angle(str, "\tAxB", angles);
1804 
1805  return 0;
1806 }
1807 
1808 
1809 /**
1810  * Free the storage associated with the rt_db_internal version of this solid.
1811  */
1812 void
1814 {
1815  RT_CK_DB_INTERNAL(ip);
1816 
1817  bu_free(ip->idb_ptr, "tgc ifree");
1818  ip->idb_ptr = ((void *)0);
1819 }
1820 
1821 struct ellipse {
1822  point_t center;
1823  vect_t axis_a;
1824  vect_t axis_b;
1825 };
1826 
1827 
1828 static void
1829 draw_lines_between_rec_ellipses(
1830  struct bu_list *vhead,
1831  struct ellipse ellipse1,
1832  vect_t h,
1833  int num_lines)
1834 {
1835  int i;
1836  point_t ellipse1_point, ellipse2_point;
1837  fastf_t radian_step = M_2PI / num_lines;
1838 
1839  for (i = 0; i < num_lines; ++i) {
1840  ellipse_point_at_radian(ellipse1_point, ellipse1.center,
1841  ellipse1.axis_a, ellipse1.axis_b, i * radian_step);
1842  VADD2(ellipse2_point, ellipse1_point, h);
1843 
1844  RT_ADD_VLIST(vhead, ellipse1_point, BN_VLIST_LINE_MOVE);
1845  RT_ADD_VLIST(vhead, ellipse2_point, BN_VLIST_LINE_DRAW);
1846  }
1847 }
1848 
1849 static void
1850 draw_lines_between_ellipses(
1851  struct bu_list *vhead,
1852  struct ellipse ellipse1,
1853  struct ellipse ellipse2,
1854  int num_lines)
1855 {
1856  int i;
1857  point_t ellipse1_point, ellipse2_point;
1858  fastf_t radian_step = M_2PI / num_lines;
1859 
1860  for (i = 0; i < num_lines; ++i) {
1861  ellipse_point_at_radian(ellipse1_point, ellipse1.center,
1862  ellipse1.axis_a, ellipse1.axis_b, i * radian_step);
1863  ellipse_point_at_radian(ellipse2_point, ellipse2.center,
1864  ellipse2.axis_a, ellipse2.axis_b, i * radian_step);
1865 
1866  RT_ADD_VLIST(vhead, ellipse1_point, BN_VLIST_LINE_MOVE);
1867  RT_ADD_VLIST(vhead, ellipse2_point, BN_VLIST_LINE_DRAW);
1868  }
1869 }
1870 
1871 static int
1872 tgc_points_per_ellipse(const struct rt_db_internal *ip, const struct rt_view_info *info)
1873 {
1874  struct rt_tgc_internal *tgc;
1875  fastf_t avg_radius, avg_circumference;
1876  fastf_t tgc_mag_a, tgc_mag_b, tgc_mag_c, tgc_mag_d;
1877 
1878  RT_CK_DB_INTERNAL(ip);
1879  tgc = (struct rt_tgc_internal *)ip->idb_ptr;
1880  RT_TGC_CK_MAGIC(tgc);
1881 
1882  tgc_mag_a = MAGNITUDE(tgc->a);
1883  tgc_mag_b = MAGNITUDE(tgc->b);
1884  tgc_mag_c = MAGNITUDE(tgc->c);
1885  tgc_mag_d = MAGNITUDE(tgc->d);
1886 
1887  avg_radius = (tgc_mag_a + tgc_mag_b + tgc_mag_c + tgc_mag_d) / 4.0;
1888  avg_circumference = M_2PI * avg_radius;
1889 
1890  return avg_circumference / info->point_spacing;
1891 }
1892 
1893 static fastf_t
1894 ramanujan_approx_circumference(fastf_t major_len, fastf_t minor_len)
1895 {
1896  fastf_t a = major_len, b = minor_len;
1897 
1898  return M_PI * (3.0 * (a + b) - sqrt(10.0 * a * b + 3.0 * (a * a + b * b)));
1899 }
1900 
1901 static int
1902 tgc_connecting_lines(
1903  const struct rt_tgc_internal *tgc,
1904  const struct rt_view_info *info)
1905 {
1906  fastf_t mag_a, mag_b, mag_c, mag_d, avg_circumference;
1907 
1908  mag_a = MAGNITUDE(tgc->a);
1909  mag_b = MAGNITUDE(tgc->b);
1910  mag_c = MAGNITUDE(tgc->c);
1911  mag_d = MAGNITUDE(tgc->d);
1912 
1913  avg_circumference = ramanujan_approx_circumference(mag_a, mag_b);
1914  avg_circumference += ramanujan_approx_circumference(mag_c, mag_d);
1915  avg_circumference /= 2.0;
1916 
1917  return avg_circumference / info->curve_spacing;
1918 }
1919 
1920 int
1921 rt_tgc_adaptive_plot(struct rt_db_internal *ip, const struct rt_view_info *info)
1922 {
1923  int points_per_ellipse, connecting_lines;
1924  struct rt_tgc_internal *tip;
1925  struct ellipse ellipse1, ellipse2;
1926 
1927  BU_CK_LIST_HEAD(info->vhead);
1928  RT_CK_DB_INTERNAL(ip);
1929  tip = (struct rt_tgc_internal *)ip->idb_ptr;
1930  RT_TGC_CK_MAGIC(tip);
1931 
1932  points_per_ellipse = tgc_points_per_ellipse(ip, info);
1933 
1934  if (points_per_ellipse < 6) {
1935  point_t p;
1936 
1937  VADD2(p, tip->v, tip->h);
1938  RT_ADD_VLIST(info->vhead, tip->v, BN_VLIST_LINE_MOVE);
1940 
1941  return 0;
1942  }
1943 
1944  connecting_lines = tgc_connecting_lines(tip, info);
1945 
1946  if (connecting_lines < 4) {
1947  connecting_lines = 4;
1948  }
1949 
1950  VMOVE(ellipse1.center, tip->v);
1951  VMOVE(ellipse1.axis_a, tip->a);
1952  VMOVE(ellipse1.axis_b, tip->b);
1953 
1954  VADD2(ellipse2.center, tip->v, tip->h);
1955  VMOVE(ellipse2.axis_a, tip->c);
1956  VMOVE(ellipse2.axis_b, tip->d);
1957 
1958  /* looks like a right elliptical cylinder */
1959  if (VNEAR_EQUAL(tip->a, tip->c, info->tol->dist) &&
1960  VNEAR_EQUAL(tip->b, tip->d, info->tol->dist))
1961  {
1962  int i;
1963  point_t *pts;
1964  fastf_t radian, radian_step;
1965 
1966  pts = (point_t *)bu_malloc(sizeof(point_t) * points_per_ellipse,
1967  "tgc points");
1968 
1969  radian_step = M_2PI / points_per_ellipse;
1970 
1971  /* calculate and plot first ellipse */
1972  ellipse_point_at_radian(pts[0], tip->v, tip->a, tip->b,
1973  radian_step * (points_per_ellipse - 1));
1974  RT_ADD_VLIST(info->vhead, pts[0], BN_VLIST_LINE_MOVE);
1975 
1976  radian = 0;
1977  for (i = 0; i < points_per_ellipse; ++i) {
1978  ellipse_point_at_radian(pts[i], tip->v, tip->a, tip->b, radian);
1979  RT_ADD_VLIST(info->vhead, pts[i], BN_VLIST_LINE_DRAW);
1980 
1981  radian += radian_step;
1982  }
1983 
1984  /* calculate and plot second ellipse */
1985  for (i = 0; i < points_per_ellipse; ++i) {
1986  VADD2(pts[i], tip->h, pts[i]);
1987  }
1988 
1989  RT_ADD_VLIST(info->vhead, pts[points_per_ellipse - 1], BN_VLIST_LINE_MOVE);
1990  for (i = 0; i < points_per_ellipse; ++i) {
1991  RT_ADD_VLIST(info->vhead, pts[i], BN_VLIST_LINE_DRAW);
1992  }
1993 
1994  bu_free(pts, "tgc points");
1995 
1996  draw_lines_between_rec_ellipses(info->vhead, ellipse1, tip->h,
1997  connecting_lines);
1998  } else {
1999  plot_ellipse(info->vhead, ellipse1.center, ellipse1.axis_a, ellipse1.axis_b,
2000  points_per_ellipse);
2001 
2002  plot_ellipse(info->vhead, ellipse2.center, ellipse2.axis_a, ellipse2.axis_b,
2003  points_per_ellipse);
2004 
2005  draw_lines_between_ellipses(info->vhead, ellipse1, ellipse2,
2006  connecting_lines);
2007  }
2008 
2009  return 0;
2010 }
2011 
2012 int
2013 rt_tgc_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct rt_tess_tol *UNUSED(ttol), const struct bn_tol *UNUSED(tol), const struct rt_view_info *UNUSED(info))
2014 {
2015  struct rt_tgc_internal *tip;
2016  register int i;
2017  fastf_t top[16*ELEMENTS_PER_VECT];
2018  fastf_t bottom[16*ELEMENTS_PER_VECT];
2019  vect_t work; /* Vec addition work area */
2020 
2021  BU_CK_LIST_HEAD(vhead);
2022  RT_CK_DB_INTERNAL(ip);
2023  tip = (struct rt_tgc_internal *)ip->idb_ptr;
2024  RT_TGC_CK_MAGIC(tip);
2025 
2026  rt_ell_16pts(bottom, tip->v, tip->a, tip->b);
2027  VADD2(work, tip->v, tip->h);
2028  rt_ell_16pts(top, work, tip->c, tip->d);
2029 
2030  /* Draw the top */
2031  RT_ADD_VLIST(vhead, &top[15*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
2032  for (i=0; i<16; i++) {
2033  RT_ADD_VLIST(vhead, &top[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
2034  }
2035 
2036  /* Draw the bottom */
2037  RT_ADD_VLIST(vhead, &bottom[15*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
2038  for (i=0; i<16; i++) {
2039  RT_ADD_VLIST(vhead, &bottom[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
2040  }
2041 
2042  /* Draw connections */
2043  for (i=0; i<16; i += 4) {
2044  RT_ADD_VLIST(vhead, &top[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
2045  RT_ADD_VLIST(vhead, &bottom[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
2046  }
2047  return 0;
2048 }
2049 
2050 
2051 /**
2052  * Return the curvature of the TGC.
2053  */
2054 void
2055 rt_tgc_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
2056 {
2057  register struct tgc_specific *tgc =
2058  (struct tgc_specific *)stp->st_specific;
2059  fastf_t R, Q, R2, Q2;
2060  mat_t M, dN, mtmp;
2061  vect_t gradf, tmp, u, v;
2062  fastf_t a, b, c, scale;
2063  vect_t vec1, vec2;
2064 
2065  if (hitp->hit_surfno != TGC_NORM_BODY) {
2066  /* We hit an end plate. Choose any tangent vector. */
2067  bn_vec_ortho(cvp->crv_pdir, hitp->hit_normal);
2068  cvp->crv_c1 = cvp->crv_c2 = 0;
2069  return;
2070  }
2071 
2072  R = 1 + tgc->tgc_CdAm1 * hitp->hit_vpriv[Z];
2073  Q = 1 + tgc->tgc_DdBm1 * hitp->hit_vpriv[Z];
2074  R2 = R*R;
2075  Q2 = Q*Q;
2076 
2077  /*
2078  * Compute derivatives of the gradient (normal) field
2079  * in ideal coords. This is a symmetric matrix with
2080  * the columns (dNx, dNy, dNz).
2081  */
2082  MAT_IDN(dN);
2083  dN[0] = Q2;
2084  dN[2] = dN[8] = 2.0*Q*tgc->tgc_DdBm1 * hitp->hit_vpriv[X];
2085  dN[5] = R2;
2086  dN[6] = dN[9] = 2.0*R*tgc->tgc_CdAm1 * hitp->hit_vpriv[Y];
2087  dN[10] = tgc->tgc_DdBm1*tgc->tgc_DdBm1 * hitp->hit_vpriv[X]*hitp->hit_vpriv[X]
2088  + tgc->tgc_CdAm1*tgc->tgc_CdAm1 * hitp->hit_vpriv[Y]*hitp->hit_vpriv[Y]
2089  - tgc->tgc_DdBm1*tgc->tgc_DdBm1 * R2
2090  - tgc->tgc_CdAm1*tgc->tgc_CdAm1 * Q2
2091  - 4.0*tgc->tgc_CdAm1*tgc->tgc_DdBm1 * R*Q;
2092 
2093  /* M = At * dN * A */
2094  bn_mat_mul(mtmp, dN, tgc->tgc_ScShR);
2095  bn_mat_mul(M, tgc->tgc_invRtShSc, mtmp);
2096 
2097  /* XXX - determine the scaling */
2098  gradf[X] = Q2 * hitp->hit_vpriv[X];
2099  gradf[Y] = R2 * hitp->hit_vpriv[Y];
2100  gradf[Z] = (hitp->hit_vpriv[X]*hitp->hit_vpriv[X] - R2) * Q * tgc->tgc_DdBm1 +
2101  (hitp->hit_vpriv[Y]*hitp->hit_vpriv[Y] - Q2) * R * tgc->tgc_CdAm1;
2102  MAT4X3VEC(tmp, tgc->tgc_invRtShSc, gradf);
2103  scale = -1.0 / MAGNITUDE(tmp);
2104  /* XXX */
2105 
2106  /*
2107  * choose a tangent plane coordinate system
2108  * (u, v, normal) form a right-handed triple
2109  */
2110  bn_vec_ortho(u, hitp->hit_normal);
2111  VCROSS(v, hitp->hit_normal, u);
2112 
2113  /* find the second fundamental form */
2114  MAT4X3VEC(tmp, M, u);
2115  a = VDOT(u, tmp) * scale;
2116  b = VDOT(v, tmp) * scale;
2117  MAT4X3VEC(tmp, M, v);
2118  c = VDOT(v, tmp) * scale;
2119 
2120  bn_eigen2x2(&cvp->crv_c1, &cvp->crv_c2, vec1, vec2, a, b, c);
2121  VCOMB2(cvp->crv_pdir, vec1[X], u, vec1[Y], v);
2122  VUNITIZE(cvp->crv_pdir);
2123 }
2124 
2125 
2126 /**
2127  * Tessellation of the TGC.
2128  *
2129  * Returns -
2130  * -1 failure
2131  * 0 OK. *r points to nmgregion that holds this tessellation.
2132  */
2133 
2134 struct tgc_pts
2135 {
2136  point_t pt;
2137  vect_t tan_axb;
2138  struct vertex *v;
2139  char dont_use;
2140 };
2141 
2142 
2143 /* version using tolerances */
2144 int
2145 rt_tgc_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
2146 {
2147  struct shell *s; /* shell to hold facetted TGC */
2148  struct faceuse *fu, *fu_top, *fu_base;
2149  struct rt_tgc_internal *tip;
2150  fastf_t radius; /* bounding sphere radius */
2151  fastf_t max_radius, min_radius; /* max/min of a, b, c, d */
2152  fastf_t h, a, b, c, d; /* lengths of TGC vectors */
2153  fastf_t inv_length; /* 1.0/length of a vector */
2154  vect_t unit_a, unit_b, unit_c, unit_d; /* units vectors in a, b, c, d directions */
2155  fastf_t rel, absolute, norm; /* interpreted tolerances */
2156  fastf_t alpha_tol; /* final tolerance for ellipse parameter */
2157  fastf_t abs_tol; /* handle invalid ttol->abs */
2158  int nells; /* total number of ellipses */
2159  int nsegs; /* number of vertices/ellipse */
2160  vect_t *A; /* array of A vectors for ellipses */
2161  vect_t *B; /* array of B vectors for ellipses */
2162  fastf_t *factors; /* array of ellipse locations along height vector */
2163  vect_t vtmp;
2164  vect_t normal; /* normal vector */
2165  vect_t rev_norm; /* reverse normal */
2166  struct tgc_pts **pts; /* array of points (pts[ellipse#][seg#]) */
2167  struct bu_ptbl verts; /* table of vertices used for top and bottom faces */
2168  struct bu_ptbl faces; /* table of faceuses for nmg_gluefaces */
2169  struct vertex **v[3]; /* array for making triangular faces */
2170 
2171  int i;
2172 
2173  VSETALL(unit_a, 0);
2174  VSETALL(unit_b, 0);
2175  VSETALL(unit_c, 0);
2176  VSETALL(unit_d, 0);
2177 
2178  RT_CK_DB_INTERNAL(ip);
2179  tip = (struct rt_tgc_internal *)ip->idb_ptr;
2180  RT_TGC_CK_MAGIC(tip);
2181 
2182  if (ttol->abs > 0.0 && ttol->abs < tol->dist) {
2183  bu_log("WARNING: tessellation tolerance is %fmm while calculational tolerance is %fmm\n",
2184  ttol->abs, tol->dist);
2185  bu_log("Cannot tessellate a TGC to finer tolerance than the calculational tolerance\n");
2186  abs_tol = tol->dist;
2187  } else {
2188  abs_tol = ttol->abs;
2189  }
2190 
2191  h = MAGNITUDE(tip->h);
2192  a = MAGNITUDE(tip->a);
2193  if (2.0*a <= tol->dist)
2194  a = 0.0;
2195  b = MAGNITUDE(tip->b);
2196  if (2.0*b <= tol->dist)
2197  b = 0.0;
2198  c = MAGNITUDE(tip->c);
2199  if (2.0*c <= tol->dist)
2200  c = 0.0;
2201  d = MAGNITUDE(tip->d);
2202  if (2.0*d <= tol->dist)
2203  d = 0.0;
2204 
2205  if (ZERO(a) && ZERO(b) && (ZERO(c) || ZERO(d))) {
2206  bu_log("Illegal TGC a, b, and c or d less than tolerance\n");
2207  return -1;
2208  } else if (ZERO(c) && ZERO(d) && (ZERO(a) || ZERO(b))) {
2209  bu_log("Illegal TGC c, d, and a or b less than tolerance\n");
2210  return -1;
2211  }
2212 
2213  if (a > 0.0) {
2214  inv_length = 1.0/a;
2215  VSCALE(unit_a, tip->a, inv_length);
2216  }
2217  if (b > 0.0) {
2218  inv_length = 1.0/b;
2219  VSCALE(unit_b, tip->b, inv_length);
2220  }
2221  if (c > 0.0) {
2222  inv_length = 1.0/c;
2223  VSCALE(unit_c, tip->c, inv_length);
2224  }
2225  if (d > 0.0) {
2226  inv_length = 1.0/d;
2227  VSCALE(unit_d, tip->d, inv_length);
2228  }
2229 
2230  /* get bounding sphere radius for relative tolerance */
2231  radius = h/2.0;
2232  max_radius = 0.0;
2233  if (a > max_radius)
2234  max_radius = a;
2235  if (b > max_radius)
2236  max_radius = b;
2237  if (c > max_radius)
2238  max_radius = c;
2239  if (d > max_radius)
2240  max_radius = d;
2241 
2242  if (max_radius > radius)
2243  radius = max_radius;
2244 
2245  min_radius = MAX_FASTF;
2246  if (a < min_radius && a > 0.0)
2247  min_radius = a;
2248  if (b < min_radius && b > 0.0)
2249  min_radius = b;
2250  if (c < min_radius && c > 0.0)
2251  min_radius = c;
2252  if (d < min_radius && d > 0.0)
2253  min_radius = d;
2254 
2255  if (abs_tol <= 0.0 && ttol->rel <= 0.0 && ttol->norm <= 0.0) {
2256  /* no tolerances specified, use 10% relative tolerance */
2257  if ((radius * 0.2) < max_radius)
2258  alpha_tol = 2.0 * acos(1.0 - 2.0 * radius * 0.1 / max_radius);
2259  else
2260  alpha_tol = M_PI_2;
2261  } else {
2262  if (abs_tol > 0.0)
2263  absolute = 2.0 * acos(1.0 - abs_tol/max_radius);
2264  else
2265  absolute = M_PI_2;
2266 
2267  if (ttol->rel > 0.0) {
2268  if (ttol->rel * 2.0 * radius < max_radius)
2269  rel = 2.0 * acos(1.0 - ttol->rel * 2.0 * radius/max_radius);
2270  else
2271  rel = M_PI_2;
2272  } else
2273  rel = M_PI_2;
2274 
2275  if (ttol->norm > 0.0) {
2276  fastf_t norm_top, norm_bot;
2277 
2278  if (a<b)
2279  norm_bot = 2.0 * atan(tan(ttol->norm) * (a/b));
2280  else
2281  norm_bot = 2.0 * atan(tan(ttol->norm) * (b/a));
2282 
2283  if (c<d)
2284  norm_top = 2.0 * atan(tan(ttol->norm) * (c/d));
2285  else
2286  norm_top = 2.0 * atan(tan(ttol->norm) * (d/c));
2287 
2288  if (norm_bot < norm_top)
2289  norm = norm_bot;
2290  else
2291  norm = norm_top;
2292  } else
2293  norm = M_PI_2;
2294 
2295  if (absolute < rel)
2296  alpha_tol = absolute;
2297  else
2298  alpha_tol = rel;
2299  if (norm < alpha_tol)
2300  alpha_tol = norm;
2301  }
2302 
2303  /* get number of segments per quadrant */
2304  nsegs = (int)(M_PI_2 / alpha_tol + 0.9999);
2305  if (nsegs < 2)
2306  nsegs = 2;
2307 
2308  /* and for complete ellipse */
2309  nsegs *= 4;
2310 
2311  /* get number and placement of intermediate ellipses */
2312  {
2313  fastf_t ratios[4], max_ratio;
2314  fastf_t new_ratio = 0;
2315  int which_ratio;
2316  fastf_t len_ha, len_hb;
2317  vect_t ha, hb;
2318  fastf_t ang;
2319  fastf_t sin_ang, cos_ang, cos_m_1_sq, sin_sq;
2320  fastf_t len_A, len_B, len_C, len_D;
2321  int bot_ell=0, top_ell=1;
2322  int reversed=0;
2323 
2324  nells = 2;
2325 
2326  max_ratio = MAX_RATIO + 1.0;
2327 
2328  factors = (fastf_t *)bu_malloc(nells*sizeof(fastf_t), "factors");
2329  A = (vect_t *)bu_malloc(nells*sizeof(vect_t), "A vectors");
2330  B = (vect_t *)bu_malloc(nells*sizeof(vect_t), "B vectors");
2331 
2332  factors[bot_ell] = 0.0;
2333  factors[top_ell] = 1.0;
2334  VMOVE(A[bot_ell], tip->a);
2335  VMOVE(A[top_ell], tip->c);
2336  VMOVE(B[bot_ell], tip->b);
2337  VMOVE(B[top_ell], tip->d);
2338 
2339  /* make sure that AxB points in the general direction of H */
2340  VCROSS(vtmp, A[0], B[0]);
2341  if (VDOT(vtmp, tip->h) < 0.0) {
2342  VMOVE(A[bot_ell], tip->b);
2343  VMOVE(A[top_ell], tip->d);
2344  VMOVE(B[bot_ell], tip->a);
2345  VMOVE(B[top_ell], tip->c);
2346  reversed = 1;
2347  }
2348  ang = M_2PI/((double)nsegs);
2349  sin_ang = sin(ang);
2350  cos_ang = cos(ang);
2351  cos_m_1_sq = (cos_ang - 1.0)*(cos_ang - 1.0);
2352  sin_sq = sin_ang*sin_ang;
2353 
2354  VJOIN2(ha, tip->h, 1.0, tip->c, -1.0, tip->a);
2355  VJOIN2(hb, tip->h, 1.0, tip->d, -1.0, tip->b);
2356  len_ha = MAGNITUDE(ha);
2357  len_hb = MAGNITUDE(hb);
2358 
2359  while (max_ratio > MAX_RATIO) {
2360  fastf_t tri_width;
2361 
2362  len_A = MAGNITUDE(A[bot_ell]);
2363  if (2.0*len_A <= tol->dist)
2364  len_A = 0.0;
2365  len_B = MAGNITUDE(B[bot_ell]);
2366  if (2.0*len_B <= tol->dist)
2367  len_B = 0.0;
2368  len_C = MAGNITUDE(A[top_ell]);
2369  if (2.0*len_C <= tol->dist)
2370  len_C = 0.0;
2371  len_D = MAGNITUDE(B[top_ell]);
2372  if (2.0*len_D <= tol->dist)
2373  len_D = 0.0;
2374 
2375  if ((len_B > 0.0 && len_D > 0.0) ||
2376  (len_B > 0.0 && (ZERO(len_D) && ZERO(len_C))))
2377  {
2378  tri_width = sqrt(cos_m_1_sq*len_A*len_A + sin_sq*len_B*len_B);
2379  ratios[0] = (factors[top_ell] - factors[bot_ell])*len_ha
2380  /tri_width;
2381  } else
2382  ratios[0] = 0.0;
2383 
2384  if ((len_A > 0.0 && len_C > 0.0) ||
2385  (len_A > 0.0 && (ZERO(len_C) && ZERO(len_D))))
2386  {
2387  tri_width = sqrt(sin_sq*len_A*len_A + cos_m_1_sq*len_B*len_B);
2388  ratios[1] = (factors[top_ell] - factors[bot_ell])*len_hb
2389  /tri_width;
2390  } else
2391  ratios[1] = 0.0;
2392 
2393  if ((len_D > 0.0 && len_B > 0.0) ||
2394  (len_D > 0.0 && (ZERO(len_A) && ZERO(len_B))))
2395  {
2396  tri_width = sqrt(cos_m_1_sq*len_C*len_C + sin_sq*len_D*len_D);
2397  ratios[2] = (factors[top_ell] - factors[bot_ell])*len_ha
2398  /tri_width;
2399  } else
2400  ratios[2] = 0.0;
2401 
2402  if ((len_C > 0.0 && len_A > 0.0) ||
2403  (len_C > 0.0 && (ZERO(len_A) && ZERO(len_B))))
2404  {
2405  tri_width = sqrt(sin_sq*len_C*len_C + cos_m_1_sq*len_D*len_D);
2406  ratios[3] = (factors[top_ell] - factors[bot_ell])*len_hb
2407  /tri_width;
2408  } else
2409  ratios[3] = 0.0;
2410 
2411  which_ratio = -1;
2412  max_ratio = 0.0;
2413 
2414  for (i=0; i<4; i++) {
2415  if (ratios[i] > max_ratio) {
2416  max_ratio = ratios[i];
2417  which_ratio = i;
2418  }
2419  }
2420 
2421  if (ZERO(len_A) && ZERO(len_B) && ZERO(len_C) && ZERO(len_D)) {
2422  if (top_ell == nells - 1) {
2423  VMOVE(A[top_ell-1], A[top_ell]);
2424  VMOVE(B[top_ell-1], A[top_ell]);
2425  factors[top_ell-1] = factors[top_ell];
2426  } else if (bot_ell == 0) {
2427  for (i=0; i<nells-1; i++) {
2428  VMOVE(A[i], A[i+1]);
2429  VMOVE(B[i], B[i+1]);
2430  factors[i] = factors[i+1];
2431  }
2432  }
2433 
2434  nells -= 1;
2435  break;
2436  }
2437 
2438  if (max_ratio <= MAX_RATIO)
2439  break;
2440 
2441  if (which_ratio == 0 || which_ratio == 1) {
2442  new_ratio = MAX_RATIO/max_ratio;
2443  if (bot_ell == 0 && new_ratio > 0.5)
2444  new_ratio = 0.5;
2445  } else if (which_ratio == 2 || which_ratio == 3) {
2446  new_ratio = 1.0 - MAX_RATIO/max_ratio;
2447  if (top_ell == nells - 1 && new_ratio < 0.5)
2448  new_ratio = 0.5;
2449  } else {
2450  /* no MAX??? */
2451  bu_log("rt_tgc_tess: Should never get here!!\n");
2452  bu_bomb("rt_tgc_tess: Should never get here!!\n");
2453  }
2454 
2455  nells++;
2456  factors = (fastf_t *)bu_realloc(factors, nells*sizeof(fastf_t), "factors");
2457  A = (vect_t *)bu_realloc(A, nells*sizeof(vect_t), "A vectors");
2458  B = (vect_t *)bu_realloc(B, nells*sizeof(vect_t), "B vectors");
2459 
2460  for (i=nells-1; i>top_ell; i--) {
2461  factors[i] = factors[i-1];
2462  VMOVE(A[i], A[i-1]);
2463  VMOVE(B[i], B[i-1]);
2464  }
2465 
2466  factors[top_ell] = factors[bot_ell] +
2467  new_ratio*(factors[top_ell+1] - factors[bot_ell]);
2468 
2469  if (reversed) {
2470  VBLEND2(A[top_ell], (1.0-factors[top_ell]), tip->b, factors[top_ell], tip->d);
2471  VBLEND2(B[top_ell], (1.0-factors[top_ell]), tip->a, factors[top_ell], tip->c);
2472  } else {
2473  VBLEND2(A[top_ell], (1.0-factors[top_ell]), tip->a, factors[top_ell], tip->c);
2474  VBLEND2(B[top_ell], (1.0-factors[top_ell]), tip->b, factors[top_ell], tip->d);
2475  }
2476 
2477  if (which_ratio == 0 || which_ratio == 1) {
2478  top_ell++;
2479  bot_ell++;
2480  }
2481 
2482  }
2483 
2484  }
2485 
2486  /* get memory for points */
2487  pts = (struct tgc_pts **)bu_calloc(nells, sizeof(struct tgc_pts *), "rt_tgc_tess: pts");
2488  for (i=0; i<nells; i++)
2489  pts[i] = (struct tgc_pts *)bu_calloc(nsegs, sizeof(struct tgc_pts), "rt_tgc_tess: pts");
2490 
2491  /* calculate geometry for points */
2492  for (i=0; i<nells; i++) {
2493  fastf_t h_factor;
2494  int j;
2495 
2496  h_factor = factors[i];
2497  for (j=0; j<nsegs; j++) {
2498  double alpha;
2499  double sin_alpha, cos_alpha;
2500 
2501  alpha = M_2PI * (double)(2*j+1)/(double)(2*nsegs);
2502  sin_alpha = sin(alpha);
2503  cos_alpha = cos(alpha);
2504 
2505  /* vertex geometry */
2506  if (i == 0 && ZERO(a) && ZERO(b)) {
2507  VMOVE(pts[i][j].pt, tip->v);
2508  } else if (i == nells-1 && ZERO(c) && ZERO(d)) {
2509  VADD2(pts[i][j].pt, tip->v, tip->h);
2510  } else {
2511  VJOIN3(pts[i][j].pt, tip->v, h_factor, tip->h, cos_alpha, A[i], sin_alpha, B[i]);
2512  }
2513 
2514  /* Storing the tangent here while sines and cosines are available */
2515  if (i == 0 && ZERO(a) && ZERO(b)) {
2516  VCOMB2(pts[0][j].tan_axb, -sin_alpha, unit_c, cos_alpha, unit_d);
2517  } else if (i == nells-1 && ZERO(c) && ZERO(d)) {
2518  VCOMB2(pts[i][j].tan_axb, -sin_alpha, unit_a, cos_alpha, unit_b);
2519  } else {
2520  VCOMB2(pts[i][j].tan_axb, -sin_alpha, A[i], cos_alpha, B[i]);
2521  }
2522  }
2523  }
2524 
2525  /* make sure no edges will be created with length < tol->dist */
2526  for (i=0; i<nells; i++) {
2527  int j;
2528  point_t curr_pt;
2529  vect_t edge_vect;
2530 
2531  if (i == 0 && (ZERO(a) || ZERO(b)))
2532  continue;
2533  else if (i == nells-1 && (ZERO(c) || ZERO(d)))
2534  continue;
2535 
2536  VMOVE(curr_pt, pts[i][0].pt);
2537  for (j=1; j<nsegs; j++) {
2538  fastf_t edge_len_sq;
2539 
2540  VSUB2(edge_vect, curr_pt, pts[i][j].pt);
2541  edge_len_sq = MAGSQ(edge_vect);
2542  if (edge_len_sq > tol->dist_sq) {
2543  VMOVE(curr_pt, pts[i][j].pt);
2544  } else {
2545  /* don't use this point, it will create a too short edge */
2546  pts[i][j].dont_use = 'n';
2547  }
2548  }
2549  }
2550 
2551  /* make region, shell, vertex */
2552  *r = nmg_mrsv(m);
2553  s = BU_LIST_FIRST(shell, &(*r)->s_hd);
2554 
2555  bu_ptbl_init(&verts, 64, " &verts ");
2556  bu_ptbl_init(&faces, 64, " &faces ");
2557  /* Make bottom face */
2558  if (a > 0.0 && b > 0.0) {
2559  for (i=nsegs-1; i>=0; i--) {
2560  /* reverse order to get outward normal */
2561  if (!pts[0][i].dont_use)
2562  bu_ptbl_ins(&verts, (long *)&pts[0][i].v);
2563  }
2564 
2565  if (BU_PTBL_END(&verts) > 2) {
2566  fu_base = nmg_cmface(s, (struct vertex ***)BU_PTBL_BASEADDR(&verts), BU_PTBL_END(&verts));
2567  bu_ptbl_ins(&faces, (long *)fu_base);
2568  } else
2569  fu_base = (struct faceuse *)NULL;
2570  } else
2571  fu_base = (struct faceuse *)NULL;
2572 
2573 
2574  /* Make top face */
2575  if (c > 0.0 && d > 0.0) {
2576  bu_ptbl_reset(&verts);
2577  for (i=0; i<nsegs; i++) {
2578  if (!pts[nells-1][i].dont_use)
2579  bu_ptbl_ins(&verts, (long *)&pts[nells-1][i].v);
2580  }
2581 
2582  if (BU_PTBL_END(&verts) > 2) {
2583  fu_top = nmg_cmface(s, (struct vertex ***)BU_PTBL_BASEADDR(&verts), BU_PTBL_END(&verts));
2584  bu_ptbl_ins(&faces, (long *)fu_top);
2585  } else
2586  fu_top = (struct faceuse *)NULL;
2587  } else
2588  fu_top = (struct faceuse *)NULL;
2589 
2590  /* Free table of vertices */
2591  bu_ptbl_free(&verts);
2592 
2593  /* Make triangular faces */
2594  for (i=0; i<nells-1; i++) {
2595  int j;
2596  struct vertex **curr_top;
2597  struct vertex **curr_bot;
2598 
2599  curr_bot = &pts[i][0].v;
2600  curr_top = &pts[i+1][0].v;
2601  for (j=0; j<nsegs; j++) {
2602  int k;
2603 
2604  k = j+1;
2605  if (k == nsegs)
2606  k = 0;
2607  if (i != 0 || a > 0.0 || b > 0.0) {
2608  if (!pts[i][k].dont_use) {
2609  v[0] = curr_bot;
2610  v[1] = &pts[i][k].v;
2611  if (i+1 == nells-1 && ZERO(c) && ZERO(d))
2612  v[2] = &pts[i+1][0].v;
2613  else
2614  v[2] = curr_top;
2615  fu = nmg_cmface(s, v, 3);
2616  bu_ptbl_ins(&faces, (long *)fu);
2617  curr_bot = &pts[i][k].v;
2618  }
2619  }
2620 
2621  if (i != nells-2 || c > 0.0 || d > 0.0) {
2622  if (!pts[i+1][k].dont_use) {
2623  v[0] = &pts[i+1][k].v;
2624  v[1] = curr_top;
2625  if (i == 0 && ZERO(a) && ZERO(b))
2626  v[2] = &pts[i][0].v;
2627  else
2628  v[2] = curr_bot;
2629  fu = nmg_cmface(s, v, 3);
2630  bu_ptbl_ins(&faces, (long *)fu);
2631  curr_top = &pts[i+1][k].v;
2632  }
2633  }
2634  }
2635  }
2636 
2637  /* Assign geometry */
2638  for (i=0; i<nells; i++) {
2639  int j;
2640 
2641  for (j=0; j<nsegs; j++) {
2642  point_t pt_geom;
2643  double alpha;
2644  double sin_alpha, cos_alpha;
2645 
2646  alpha = M_2PI * (double)(2*j+1)/(double)(2*nsegs);
2647  sin_alpha = sin(alpha);
2648  cos_alpha = cos(alpha);
2649 
2650  /* vertex geometry */
2651  if (i == 0 && ZERO(a) && ZERO(b)) {
2652  if (j == 0)
2653  nmg_vertex_gv(pts[0][0].v, tip->v);
2654  } else if (i == nells-1 && ZERO(c) && ZERO(d)) {
2655  if (j == 0) {
2656  VADD2(pt_geom, tip->v, tip->h);
2657  nmg_vertex_gv(pts[i][0].v, pt_geom);
2658  }
2659  } else if (pts[i][j].v)
2660  nmg_vertex_gv(pts[i][j].v, pts[i][j].pt);
2661 
2662  /* Storing the tangent here while sines and cosines are available */
2663  if (i == 0 && ZERO(a) && ZERO(b)) {
2664  VCOMB2(pts[0][j].tan_axb, -sin_alpha, unit_c, cos_alpha, unit_d);
2665  } else if (i == nells-1 && ZERO(c) && ZERO(d)) {
2666  VCOMB2(pts[i][j].tan_axb, -sin_alpha, unit_a, cos_alpha, unit_b);
2667  } else {
2668  VCOMB2(pts[i][j].tan_axb, -sin_alpha, A[i], cos_alpha, B[i]);
2669  }
2670  }
2671  }
2672 
2673  /* Associate face plane equations */
2674  for (i=0; i<BU_PTBL_END(&faces); i++) {
2675  fu = (struct faceuse *)BU_PTBL_GET(&faces, i);
2676  NMG_CK_FACEUSE(fu);
2677 
2678  if (nmg_calc_face_g(fu)) {
2679  bu_log("rt_tess_tgc: failed to calculate plane equation\n");
2680  nmg_pr_fu_briefly(fu, "");
2681  return -1;
2682  }
2683  }
2684 
2685 
2686  /* Calculate vertexuse normals */
2687  for (i=0; i<nells; i++) {
2688  int j, k;
2689 
2690  k = i + 1;
2691  if (k == nells)
2692  k = i - 1;
2693 
2694  for (j=0; j<nsegs; j++) {
2695  vect_t tan_h; /* vector tangent from one ellipse to next */
2696  struct vertexuse *vu;
2697 
2698  /* normal at vertex */
2699  if (i == nells - 1) {
2700  if (ZERO(c) && ZERO(d)) {
2701  VSUB2(tan_h, pts[i][0].pt, pts[k][j].pt);
2702  } else if (k == 0 && ZERO(c) && ZERO(d)) {
2703  VSUB2(tan_h, pts[i][j].pt, pts[k][0].pt);
2704  } else {
2705  VSUB2(tan_h, pts[i][j].pt, pts[k][j].pt);
2706  }
2707  } else if (i == 0) {
2708  if (ZERO(a) && ZERO(b)) {
2709  VSUB2(tan_h, pts[k][j].pt, pts[i][0].pt);
2710  } else if (k == nells-1 && ZERO(c) && ZERO(d)) {
2711  VSUB2(tan_h, pts[k][0].pt, pts[i][j].pt);
2712  } else {
2713  VSUB2(tan_h, pts[k][j].pt, pts[i][j].pt);
2714  }
2715  } else if (k == 0 && ZERO(a) && ZERO(b)) {
2716  VSUB2(tan_h, pts[k][0].pt, pts[i][j].pt);
2717  } else if (k == nells-1 && ZERO(c) && ZERO(d)) {
2718  VSUB2(tan_h, pts[k][0].pt, pts[i][j].pt);
2719  } else {
2720  VSUB2(tan_h, pts[k][j].pt, pts[i][j].pt);
2721  }
2722 
2723  VCROSS(normal, pts[i][j].tan_axb, tan_h);
2724  VUNITIZE(normal);
2725  VREVERSE(rev_norm, normal);
2726 
2727  if (!(i == 0 && ZERO(a) && ZERO(b)) &&
2728  !(i == nells-1 && ZERO(c) && ZERO(d)) &&
2729  pts[i][j].v)
2730  {
2731  for (BU_LIST_FOR(vu, vertexuse, &pts[i][j].v->vu_hd)) {
2732  NMG_CK_VERTEXUSE(vu);
2733 
2734  fu = nmg_find_fu_of_vu(vu);
2735  NMG_CK_FACEUSE(fu);
2736 
2737  /* don't need vertexuse normals for faces that are really flat */
2738  if (fu == fu_base || fu->fumate_p == fu_base ||
2739  fu == fu_top || fu->fumate_p == fu_top)
2740  continue;
2741 
2742  if (fu->orientation == OT_SAME)
2743  nmg_vertexuse_nv(vu, normal);
2744  else if (fu->orientation == OT_OPPOSITE)
2745  nmg_vertexuse_nv(vu, rev_norm);
2746  }
2747  }
2748  }
2749  }
2750 
2751  /* Finished with storage, so free it */
2752  bu_free((char *)factors, "rt_tgc_tess: factors");
2753  bu_free((char *)A, "rt_tgc_tess: A");
2754  bu_free((char *)B, "rt_tgc_tess: B");
2755  for (i=0; i<nells; i++)
2756  bu_free((char *)pts[i], "rt_tgc_tess: pts[i]");
2757  bu_free((char *)pts, "rt_tgc_tess: pts");
2758 
2759  /* mark real edges for top and bottom faces */
2760  for (i=0; i<2; i++) {
2761  struct loopuse *lu;
2762 
2763  if (i == 0)
2764  fu = fu_base;
2765  else
2766  fu = fu_top;
2767 
2768  if (fu == (struct faceuse *)NULL)
2769  continue;
2770 
2771  NMG_CK_FACEUSE(fu);
2772 
2773  for (BU_LIST_FOR(lu, loopuse, &fu->lu_hd)) {
2774  struct edgeuse *eu;
2775 
2776  NMG_CK_LOOPUSE(lu);
2777 
2778  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
2779  continue;
2780 
2781  for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
2782  struct edge *e;
2783 
2784  NMG_CK_EDGEUSE(eu);
2785  e = eu->e_p;
2786  NMG_CK_EDGE(e);
2787  e->is_real = 1;
2788  }
2789  }
2790  }
2791 
2792  nmg_region_a(*r, tol);
2793 
2794  /* glue faces together */
2795  nmg_gluefaces((struct faceuse **)BU_PTBL_BASEADDR(&faces), BU_PTBL_END(&faces), tol);
2796  bu_ptbl_free(&faces);
2797 
2798  return 0;
2799 }
2800 
2801 
2802 /**
2803  * "Tessellate an TGC into a trimmed-NURB-NMG data structure.
2804  * Computing NURB surfaces and trimming curves to interpolate
2805  * the parameters of the TGC
2806  *
2807  * The process is to create the nmg topology of the TGC fill it
2808  * in with a unit cylinder geometry (i.e. unitcircle at the top (0, 0, 1)
2809  * unit cylinder of radius 1, and unitcircle at the bottom), and then
2810  * scale it with a perspective matrix derived from the parameters of the
2811  * tgc. The result is three trimmed nurb surfaces which interpolate the
2812  * parameters of the original TGC.
2813  *
2814  * Returns -
2815  * -1 failure
2816  * 0 OK. *r points to nmgregion that holds this tessellation
2817  */
2818 
2819 int
2820 rt_tgc_tnurb(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct bn_tol *tol)
2821 {
2822  struct rt_tgc_internal *tip;
2823 
2824  struct shell *s;
2825  struct vertex *verts[2];
2826  struct vertex **vertp[4];
2827  struct faceuse * top_fu;
2828  struct faceuse * cyl_fu;
2829  struct faceuse * bot_fu;
2830  vect_t uvw;
2831  struct loopuse *lu;
2832  struct edgeuse *eu;
2833  struct edgeuse *top_eu;
2834  struct edgeuse *bot_eu;
2835 
2836  mat_t mat;
2837  mat_t imat, omat, top_mat, bot_mat;
2838  vect_t anorm;
2839  vect_t bnorm;
2840  vect_t cnorm;
2841 
2842 
2843  /* Get the internal representation of the tgc */
2844 
2845  RT_CK_DB_INTERNAL(ip);
2846  tip = (struct rt_tgc_internal *) ip->idb_ptr;
2847  RT_TGC_CK_MAGIC(tip);
2848 
2849  /* Create the NMG Topology */
2850 
2851  *r = nmg_mrsv(m); /* Make region, empty shell, vertex */
2852  s = BU_LIST_FIRST(shell, &(*r)->s_hd);
2853 
2854 
2855  /* Create transformation matrix for the top cap surface*/
2856 
2857  MAT_IDN(omat);
2858  MAT_IDN(mat);
2859 
2860  omat[0] = MAGNITUDE(tip->c);
2861  omat[5] = MAGNITUDE(tip->d);
2862  omat[3] = tip->v[0] + tip->h[0];
2863  omat[7] = tip->v[1] + tip->h[1];
2864  omat[11] = tip->v[2] + tip->h[2];
2865 
2866  bn_mat_mul(imat, mat, omat);
2867 
2868  VMOVE(anorm, tip->c);
2869  VMOVE(bnorm, tip->d);
2870  VCROSS(cnorm, tip->c, tip->d);
2871  VUNITIZE(anorm);
2872  VUNITIZE(bnorm);
2873  VUNITIZE(cnorm);
2874 
2875  MAT_IDN(omat);
2876 
2877  VMOVE(&omat[0], anorm);
2878  VMOVE(&omat[4], bnorm);
2879  VMOVE(&omat[8], cnorm);
2880 
2881 
2882  bn_mat_mul(top_mat, omat, imat);
2883 
2884  /* Create topology for top cap surface */
2885 
2886  verts[0] = verts[1] = NULL;
2887  vertp[0] = &verts[0];
2888  top_fu = nmg_cmface(s, vertp, 1);
2889 
2890  lu = BU_LIST_FIRST(loopuse, &top_fu->lu_hd);
2891  NMG_CK_LOOPUSE(lu);
2892  eu= BU_LIST_FIRST(edgeuse, &lu->down_hd);
2893  NMG_CK_EDGEUSE(eu);
2894  top_eu = eu;
2895 
2896  VSET(uvw, 0, 0, 0);
2897  nmg_vertexuse_a_cnurb(eu->vu_p, uvw);
2898  VSET(uvw, 1, 0, 0);
2899  nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw);
2900  eu = BU_LIST_NEXT(edgeuse, &eu->l);
2901 
2902  /* Top cap surface */
2903  nmg_tgc_disk(top_fu, top_mat, 0.0, 0);
2904 
2905  /* Create transformation matrix for the bottom cap surface*/
2906 
2907  MAT_IDN(omat);
2908  MAT_IDN(mat);
2909 
2910  omat[0] = MAGNITUDE(tip->a);
2911  omat[5] = MAGNITUDE(tip->b);
2912  omat[3] = tip->v[0];
2913  omat[7] = tip->v[1];
2914  omat[11] = tip->v[2];
2915 
2916  bn_mat_mul(imat, mat, omat);
2917 
2918  VMOVE(anorm, tip->a);
2919  VMOVE(bnorm, tip->b);
2920  VCROSS(cnorm, tip->a, tip->b);
2921  VUNITIZE(anorm);
2922  VUNITIZE(bnorm);
2923  VUNITIZE(cnorm);
2924 
2925  MAT_IDN(omat);
2926 
2927  VMOVE(&omat[0], anorm);
2928  VMOVE(&omat[4], bnorm);
2929  VMOVE(&omat[8], cnorm);
2930 
2931  bn_mat_mul(bot_mat, omat, imat);
2932 
2933  /* Create topology for bottom cap surface */
2934 
2935  vertp[0] = &verts[1];
2936  bot_fu = nmg_cmface(s, vertp, 1);
2937 
2938  lu = BU_LIST_FIRST(loopuse, &bot_fu->lu_hd);
2939  NMG_CK_LOOPUSE(lu);
2940  eu= BU_LIST_FIRST(edgeuse, &lu->down_hd);
2941  NMG_CK_EDGEUSE(eu);
2942  bot_eu = eu;
2943 
2944  VSET(uvw, 0, 0, 0);
2945  nmg_vertexuse_a_cnurb(eu->vu_p, uvw);
2946  VSET(uvw, 1, 0, 0);
2947  nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw);
2948 
2949 
2950  nmg_tgc_disk(bot_fu, bot_mat, 0.0, 1);
2951 
2952  /* Create topology for cylinder surface */
2953 
2954  vertp[0] = &verts[0];
2955  vertp[1] = &verts[0];
2956  vertp[2] = &verts[1];
2957  vertp[3] = &verts[1];
2958  cyl_fu = nmg_cmface(s, vertp, 4);
2959 
2960  nmg_tgc_nurb_cyl(cyl_fu, top_mat, bot_mat);
2961 
2962  /* fuse top cylinder edge to matching edge on body of cylinder */
2963  lu = BU_LIST_FIRST(loopuse, &cyl_fu->lu_hd);
2964 
2965  eu= BU_LIST_FIRST(edgeuse, &lu->down_hd);
2966 
2967  nmg_je(top_eu, eu);
2968 
2969  eu= BU_LIST_LAST(edgeuse, &lu->down_hd);
2970  eu= BU_LIST_LAST(edgeuse, &eu->l);
2971 
2972  nmg_je(bot_eu, eu);
2973  nmg_region_a(*r, tol);
2974 
2975  return 0;
2976 }
2977 
2978 
2979 static void
2980 nmg_tgc_disk(struct faceuse *fu, fastf_t *rmat, fastf_t height, int flip)
2981 {
2982  struct face_g_snurb * fg;
2983  struct loopuse * lu;
2984  struct edgeuse * eu;
2985  struct edge_g_cnurb * eg;
2986  fastf_t *mptr;
2987  int i;
2988  vect_t vect;
2989  point_t point;
2990 
2991  nmg_face_g_snurb(fu,
2992  2, 2, /* u, v order */
2993  4, 4, /* number of knots */
2994  NULL, NULL, /* initial knot vectors */
2995  2, 2, /* n_rows, n_cols */
2996  RT_NURB_MAKE_PT_TYPE(3, 2, 0),
2997  NULL); /* Initial mesh */
2998 
2999  fg = fu->f_p->g.snurb_p;
3000 
3001 
3002  fg->u.knots[0] = 0.0;
3003  fg->u.knots[1] = 0.0;
3004  fg->u.knots[2] = 1.0;
3005  fg->u.knots[3] = 1.0;
3006 
3007  fg->v.knots[0] = 0.0;
3008  fg->v.knots[1] = 0.0;
3009  fg->v.knots[2] = 1.0;
3010  fg->v.knots[3] = 1.0;
3011 
3012  if (flip) {
3013  fg->ctl_points[0] = 1.;
3014  fg->ctl_points[1] = -1.;
3015  fg->ctl_points[2] = height;
3016 
3017  fg->ctl_points[3] = -1;
3018  fg->ctl_points[4] = -1.;
3019  fg->ctl_points[5] = height;
3020 
3021  fg->ctl_points[6] = 1.;
3022  fg->ctl_points[7] = 1.;
3023  fg->ctl_points[8] = height;
3024 
3025  fg->ctl_points[9] = -1.;
3026  fg->ctl_points[10] = 1.;
3027  fg->ctl_points[11] = height;
3028  } else {
3029 
3030  fg->ctl_points[0] = -1.;
3031  fg->ctl_points[1] = -1.;
3032  fg->ctl_points[2] = height;
3033 
3034  fg->ctl_points[3] = 1;
3035  fg->ctl_points[4] = -1.;
3036  fg->ctl_points[5] = height;
3037 
3038  fg->ctl_points[6] = -1.;
3039  fg->ctl_points[7] = 1.;
3040  fg->ctl_points[8] = height;
3041 
3042  fg->ctl_points[9] = 1.;
3043  fg->ctl_points[10] = 1.;
3044  fg->ctl_points[11] = height;
3045  }
3046 
3047  /* multiple the matrix to get orientation and scaling that we want */
3048  mptr = fg->ctl_points;
3049 
3050  i = fg->s_size[0] * fg->s_size[1];
3051 
3052  for (; i> 0; i--) {
3053  MAT4X3PNT(vect, rmat, mptr);
3054  mptr[0] = vect[0];
3055  mptr[1] = vect[1];
3056  mptr[2] = vect[2];
3057  mptr += 3;
3058  }
3059 
3060  lu = BU_LIST_FIRST(loopuse, &fu->lu_hd);
3061  NMG_CK_LOOPUSE(lu);
3062  eu= BU_LIST_FIRST(edgeuse, &lu->down_hd);
3063  NMG_CK_EDGEUSE(eu);
3064 
3065 
3066  if (!flip) {
3067  rt_nurb_s_eval(fu->f_p->g.snurb_p,
3068  nmg_uv_unitcircle[0], nmg_uv_unitcircle[1], point);
3069  nmg_vertex_gv(eu->vu_p->v_p, point);
3070  } else {
3071  rt_nurb_s_eval(fu->f_p->g.snurb_p,
3072  nmg_uv_unitcircle[12], nmg_uv_unitcircle[13], point);
3073  nmg_vertex_gv(eu->vu_p->v_p, point);
3074  }
3075 
3076  nmg_edge_g_cnurb(eu, 3, 12, NULL, 9, RT_NURB_MAKE_PT_TYPE(3, 3, 1),
3077  NULL);
3078 
3079  eg = eu->g.cnurb_p;
3080  eg->order = 3;
3081 
3082  eg->k.knots[0] = 0.0;
3083  eg->k.knots[1] = 0.0;
3084  eg->k.knots[2] = 0.0;
3085  eg->k.knots[3] = .25;
3086  eg->k.knots[4] = .25;
3087  eg->k.knots[5] = .5;
3088  eg->k.knots[6] = .5;
3089  eg->k.knots[7] = .75;
3090  eg->k.knots[8] = .75;
3091  eg->k.knots[9] = 1.0;
3092  eg->k.knots[10] = 1.0;
3093  eg->k.knots[11] = 1.0;
3094 
3095  if (!flip) {
3096  for (i = 0; i < 27; i++)
3097  eg->ctl_points[i] = nmg_uv_unitcircle[i];
3098  } else {
3099 
3100  VSET(&eg->ctl_points[0], 0.0, .5, 1.0);
3101  VSET(&eg->ctl_points[3], 0.0, 0.0, RAT);
3102  VSET(&eg->ctl_points[6], 0.5, 0.0, 1.0);
3103  VSET(&eg->ctl_points[9], RAT, 0.0, RAT);
3104  VSET(&eg->ctl_points[12], 1.0, .5, 1.0);
3105  VSET(&eg->ctl_points[15], RAT, RAT, RAT);
3106  VSET(&eg->ctl_points[18], .5, 1.0, 1.0);
3107  VSET(&eg->ctl_points[21], 0.0, RAT, RAT);
3108  VSET(&eg->ctl_points[24], 0.0, .5, 1.0);
3109  }
3110 }
3111 
3112 
3113 /* Create a cylinder with a top surface and a bottom surface
3114  * defined by the ellipsoids at the top and bottom of the
3115  * cylinder, the top_mat, and bot_mat are applied to a unit circle
3116  * for the top row of the surface and the bot row of the surface
3117  * respectively.
3118  */
3119 static void
3120 nmg_tgc_nurb_cyl(struct faceuse *fu, fastf_t *top_mat, fastf_t *bot_mat)
3121 {
3122  struct face_g_snurb * fg;
3123  struct loopuse * lu;
3124  struct edgeuse * eu;
3125  fastf_t * mptr;
3126  int i;
3127  hvect_t vect;
3128  point_t uvw;
3129  point_t point;
3130  hvect_t hvect;
3131 
3132  nmg_face_g_snurb(fu,
3133  3, 2,
3134  12, 4,
3135  NULL, NULL,
3136  2, 9,
3137  RT_NURB_MAKE_PT_TYPE(4, 3, 1),
3138  NULL);
3139 
3140  fg = fu->f_p->g.snurb_p;
3141 
3142  fg->v.knots[0] = 0.0;
3143  fg->v.knots[1] = 0.0;
3144  fg->v.knots[2] = 1.0;
3145  fg->v.knots[3] = 1.0;
3146 
3147  fg->u.knots[0] = 0.0;
3148  fg->u.knots[1] = 0.0;
3149  fg->u.knots[2] = 0.0;
3150  fg->u.knots[3] = .25;
3151  fg->u.knots[4] = .25;
3152  fg->u.knots[5] = .5;
3153  fg->u.knots[6] = .5;
3154  fg->u.knots[7] = .75;
3155  fg->u.knots[8] = .75;
3156  fg->u.knots[9] = 1.0;
3157  fg->u.knots[10] = 1.0;
3158  fg->u.knots[11] = 1.0;
3159 
3160  mptr = fg->ctl_points;
3161 
3162  for (i = 0; i < 9; i++) {
3163  MAT4X4PNT(vect, top_mat, &nmg_tgc_unitcircle[i*4]);
3164  mptr[0] = vect[0];
3165  mptr[1] = vect[1];
3166  mptr[2] = vect[2];
3167  mptr[3] = vect[3];
3168  mptr += 4;
3169  }
3170 
3171  for (i = 0; i < 9; i++) {
3172  MAT4X4PNT(vect, bot_mat, &nmg_tgc_unitcircle[i*4]);
3173  mptr[0] = vect[0];
3174  mptr[1] = vect[1];
3175  mptr[2] = vect[2];
3176  mptr[3] = vect[3];
3177  mptr += 4;
3178  }
3179 
3180  /* Assign edgeuse parameters & vertex geometry */
3181 
3182  lu = BU_LIST_FIRST(loopuse, &fu->lu_hd);
3183  NMG_CK_LOOPUSE(lu);
3184  eu = BU_LIST_FIRST(edgeuse, &lu->down_hd);
3185  NMG_CK_EDGEUSE(eu);
3186 
3187  /* March around the fu's loop assigning uv parameter values */
3188 
3189  rt_nurb_s_eval(fg, 0.0, 0.0, hvect);
3190  HDIVIDE(point, hvect);
3191  nmg_vertex_gv(eu->vu_p->v_p, point); /* 0, 0 vertex */
3192 
3193  VSET(uvw, 0, 0, 0);
3194  nmg_vertexuse_a_cnurb(eu->vu_p, uvw);
3195  VSET(uvw, 1, 0, 0);
3196  nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw);
3197  eu = BU_LIST_NEXT(edgeuse, &eu->l);
3198 
3199  VSET(uvw, 1, 0, 0);
3200  nmg_vertexuse_a_cnurb(eu->vu_p, uvw);
3201  VSET(uvw, 1, 1, 0);
3202  nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw);
3203  eu = BU_LIST_NEXT(edgeuse, &eu->l);
3204 
3205  VSET(uvw, 1, 1, 0);
3206  nmg_vertexuse_a_cnurb(eu->vu_p, uvw);
3207  VSET(uvw, 0, 1, 0);
3208  nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw);
3209 
3210  rt_nurb_s_eval(fg, 1., 1., hvect);
3211  HDIVIDE(point, hvect);
3212  nmg_vertex_gv(eu->vu_p->v_p, point); /* 4, 1 vertex */
3213 
3214  eu = BU_LIST_NEXT(edgeuse, &eu->l);
3215 
3216  VSET(uvw, 0, 1, 0);
3217  nmg_vertexuse_a_cnurb(eu->vu_p, uvw);
3218  VSET(uvw, 0, 0, 0);
3219  nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, uvw);
3220  eu = BU_LIST_NEXT(edgeuse, &eu->l);
3221 
3222  /* Create the edge loop geometry */
3223 
3224  lu = BU_LIST_FIRST(loopuse, &fu->lu_hd);
3225  NMG_CK_LOOPUSE(lu);
3226  eu = BU_LIST_FIRST(edgeuse, &lu->down_hd);
3227  NMG_CK_EDGEUSE(eu);
3228 
3229  for (i = 0; i < 4; i++) {
3231  eu = BU_LIST_NEXT(edgeuse, &eu->l);
3232  }
3233 }
3234 
3235 
3236 int
3237 rt_tgc_params(struct pc_pc_set *UNUSED(ps), const struct rt_db_internal *ip)
3238 {
3239  if (ip) RT_CK_DB_INTERNAL(ip);
3240 
3241  return 0; /* OK */
3242 }
3243 
3244 
3245 void
3246 rt_tgc_volume(fastf_t *vol, const struct rt_db_internal *ip)
3247 {
3248  int tgc_type = 0;
3249  fastf_t mag_a, mag_b, mag_c, mag_d, mag_h;
3250  struct rt_tgc_internal *tip = (struct rt_tgc_internal *)ip->idb_ptr;
3251  RT_TGC_CK_MAGIC(tip);
3252 
3253  mag_a = MAGNITUDE(tip->a);
3254  mag_b = MAGNITUDE(tip->b);
3255  mag_c = MAGNITUDE(tip->c);
3256  mag_d = MAGNITUDE(tip->d);
3257  mag_h = MAGNITUDE(tip->h);
3258 
3259  GET_TGC_TYPE(tgc_type, mag_a, mag_b, mag_c, mag_d);
3260 
3261  switch (tgc_type) {
3262  case RCC:
3263  case REC:
3264  *vol = M_PI * mag_h * mag_a * mag_b;
3265  break;
3266  case TRC:
3267  /* TRC could fall through, but this formula avoids a sqrt and
3268  * so will probably be more accurate */
3269  *vol = M_PI * mag_h * (mag_a * mag_a + mag_c * mag_c + mag_a * mag_c) / 3.0;
3270  break;
3271  case TEC:
3272  *vol = M_PI * mag_h * (mag_a * mag_b + mag_c * mag_d + sqrt(mag_a * mag_b * mag_c * mag_d)) / 3.0;
3273  break;
3274  default:
3275  /* never reached */
3276  bu_log("rt_tgc_volume(): cannot find volume\n");
3277  }
3278 }
3279 
3280 void
3281 rt_tgc_surf_area(fastf_t *area, const struct rt_db_internal *ip)
3282 {
3283  int tgc_type = 0;
3284  fastf_t mag_a, mag_b, mag_c, mag_d, mag_h;
3285  fastf_t magsq_a, magsq_c, magsq_h;
3286  fastf_t c;
3287  fastf_t area_base;
3288  struct rt_tgc_internal *tip = (struct rt_tgc_internal *)ip->idb_ptr;
3289  RT_TGC_CK_MAGIC(tip);
3290 
3291  magsq_a = MAGSQ(tip->a);
3292  magsq_c = MAGSQ(tip->c);
3293  magsq_h = MAGSQ(tip->h);
3294 
3295  mag_a = sqrt(magsq_a);
3296  mag_b = MAGNITUDE(tip->b);
3297  mag_c = sqrt(magsq_c);
3298  mag_d = MAGNITUDE(tip->d);
3299  mag_h = sqrt(magsq_h);
3300 
3301  GET_TGC_TYPE(tgc_type, mag_a, mag_b, mag_c, mag_d);
3302 
3303  switch (tgc_type) {
3304  case RCC:
3305  *area = M_2PI * mag_a * (mag_a + mag_h);
3306  break;
3307  case TRC:
3308  *area = M_PI * ((mag_a + mag_c) * sqrt((mag_a - mag_c) * (mag_a - mag_c) + magsq_h) + magsq_a + magsq_c);
3309  break;
3310  case REC:
3311  area_base = M_PI * mag_a * mag_b;
3312  /* approximation */
3313  c = ELL_CIRCUMFERENCE(mag_a, mag_b);
3314  *area = c * mag_h + 2.0 * area_base;
3315  break;
3316  case TEC:
3317  default:
3318  bu_log("rt_tgc_surf_area(): cannot find surface area\n");
3319  }
3320 }
3321 
3322 
3323 void
3324 rt_tgc_centroid(point_t *cent, const struct rt_db_internal *ip)
3325 {
3326  int tgc_type = 0;
3327  fastf_t mag_a, mag_b, mag_c, mag_d;
3328  fastf_t magsq_a, magsq_c;
3329  fastf_t scalar;
3330  struct rt_tgc_internal *tip = (struct rt_tgc_internal *)ip->idb_ptr;
3331  RT_TGC_CK_MAGIC(tip);
3332 
3333  magsq_a = MAGSQ(tip->a);
3334  magsq_c = MAGSQ(tip->c);
3335 
3336  mag_a = sqrt(magsq_a);
3337  mag_b = MAGNITUDE(tip->b);
3338  mag_c = sqrt(magsq_c);
3339  mag_d = MAGNITUDE(tip->d);
3340 
3341  GET_TGC_TYPE(tgc_type, mag_a, mag_b, mag_c, mag_d);
3342 
3343  switch (tgc_type) {
3344  case RCC:
3345  case REC:
3346  scalar = 0.5;
3347  VJOIN1(*cent, tip->v, scalar, tip->h);
3348  break;
3349  case TRC:
3350  scalar = 0.25 * (magsq_a + 2.0 * mag_a * mag_c + 3.0 * magsq_c) /
3351  (magsq_a + mag_a * mag_c + magsq_c);
3352  VJOIN1(*cent, tip->v, scalar, tip->h);
3353  break;
3354  case TEC:
3355  /* need to confirm formula */
3356  default:
3357  bu_log("rt_tgc_centroid(): cannot find centroid\n");
3358  }
3359 }
3360 
3361 
3362 /*
3363  * Local Variables:
3364  * mode: C
3365  * tab-width: 8
3366  * indent-tabs-mode: t
3367  * c-file-style: "stroustrup"
3368  * End:
3369  * ex: shiftwidth=4 tabstop=8
3370  */
Definition: db_flip.c:35
#define RT_TGC_INTERNAL_MAGIC
Definition: magic.h:111
char * d_namep
pointer to name string
Definition: raytrace.h:859
struct xray a_ray
Actual ray to be shot.
Definition: raytrace.h:1583
#define BU_LIST_FOR(p, structure, hp)
Definition: list.h:365
Definition: raytrace.h:800
fastf_t cf[BN_MAX_POLY_DEGREE+1]
Definition: poly.h:50
#define NMG_EDGEUSE_MAGIC
Definition: magic.h:120
int rt_tgc_import5(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
Definition: tgc.c:1653
void nmg_edge_g_cnurb(struct edgeuse *eu, int order, int n_knots, fastf_t *kv, int n_pts, int pt_type, fastf_t *points)
Definition: nmg_mk.c:1899
#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_tgc_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *tol)
Definition: tgc.c:155
fastf_t tgc_sH
Definition: tgc.c:52
#define SIZEOF_NETWORK_DOUBLE
Definition: cv.h:48
struct hit seg_in
IN information.
Definition: raytrace.h:370
fastf_t tgc_DdBm1
Definition: tgc.c:58
int(* ft_bbox)(struct rt_db_internal *, point_t *, point_t *, const struct bn_tol *)
Definition: raytrace.h:2196
fastf_t C[2 *MAX_CNT+1][2 *MAX_CNT+1]
Definition: dsp_brep.cpp:38
#define IN
Definition: tgc.c:84
Definition: list.h:118
const struct directory * st_dp
Directory entry of solid.
Definition: raytrace.h:436
#define BU_LIST_LAST(structure, hp)
Definition: list.h:306
#define RT_DOT_TOL
Definition: raytrace.h:170
const struct bn_tol * tol
Definition: raytrace.h:1927
fastf_t tgc_CdAm1
Definition: tgc.c:57
#define SMALL
Definition: defines.h:351
#define RT_CK_APPLICATION(_p)
Definition: raytrace.h:1675
int rt_tgc_params(struct pc_pc_set *ps, const struct rt_db_internal *ip)
Definition: tgc.c:3237
fastf_t tgc_C
Definition: tgc.c:55
int rt_rec_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: rec.c:275
Definition: tgc.c:2134
double dist
>= 0
Definition: tol.h:73
#define RT_TGC_SEG_MISS(SEG)
Definition: tgc.c:91
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
int rt_tgc_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: tgc.c:1616
void nmg_je(struct edgeuse *eudst, struct edgeuse *eusrc)
Definition: nmg_mk.c:2765
if lu s
Definition: nmg_mod.c:3860
void bu_ptbl_init(struct bu_ptbl *b, size_t len, const char *str)
Definition: ptbl.c:32
void bu_vls_strcat(struct bu_vls *vp, const char *s)
Definition: vls.c:368
lu
Definition: nmg_mod.c:3855
#define VSET(a, b, c, d)
Definition: color.c:53
#define ALPHA(x, y, c, d)
Definition: tgc.c:93
int rt_tgc_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
Definition: tgc.c:1730
#define VSETALL(a, s)
Definition: color.c:54
Definition: raytrace.h:215
void bn_pr_roots(const char *title, const struct bn_complex roots[], int n)
char dont_use
Definition: tgc.c:2139
#define GET_TGC_TYPE(type, a, b, c, d)
Definition: tgc.c:96
#define M_PI
Definition: fft.h:35
Definition: pc.h:108
double dist_sq
dist * dist
Definition: tol.h:74
Definition: raytrace.h:368
void rt_tgc_surf_area(fastf_t *area, const struct rt_db_internal *ip)
Definition: tgc.c:3281
#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
fastf_t st_aradius
Radius of APPROXIMATING sphere.
Definition: raytrace.h:433
vect_t tgc_N
Definition: tgc.c:61
void rt_tgc_ifree(struct rt_db_internal *ip)
Definition: tgc.c:1813
Header file for the BRL-CAD common definitions.
void nmg_vertexuse_nv(struct vertexuse *vu, const fastf_t *norm)
Definition: nmg_mk.c:1719
void bn_eigen2x2(fastf_t *val1, fastf_t *val2, vect_t vec1, vect_t vec2, fastf_t a, fastf_t b, fastf_t c)
Definition: poly.h:47
fastf_t tgc_A
Definition: tgc.c:53
int rt_tgc_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
Definition: tgc.c:2145
void flip_fastf_float(fastf_t *ff, const dbfloat_t *fp, int n, int flip)
Definition: db_flip.c:74
int bu_ptbl_ins(struct bu_ptbl *b, long *p)
struct resource * a_resource
dynamic memory resources
Definition: raytrace.h:1591
void nmg_face_g_snurb(struct faceuse *fu, int u_order, int v_order, int n_u_knots, int n_v_knots, fastf_t *ukv, fastf_t *vkv, int n_rows, int n_cols, int pt_type, fastf_t *mesh)
Definition: nmg_mk.c:2347
void bu_cv_htond(unsigned char *out, const unsigned char *in, size_t count)
#define MAX_FASTF
Definition: defines.h:340
void bu_ptbl_reset(struct bu_ptbl *b)
Definition: ptbl.c:49
int nmg_calc_face_g(struct faceuse *fu)
Definition: nmg_misc.c:1786
#define OUT
Definition: tgc.c:83
NMG_CK_LOOPUSE(lu)
Definition: tgc.c:1821
struct bu_list l
Definition: raytrace.h:369
fastf_t curve_spacing
Definition: raytrace.h:1940
double re
real part
Definition: complex.h:40
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
Definition: ptbl.h:62
double rel
rel dist tol
Definition: raytrace.h:181
point_t pt
Definition: tgc.c:2136
void rt_tgc_uv(struct application *ap, struct soltab *stp, register struct hit *hitp, register struct uvcoord *uvp)
Definition: tgc.c:1504
if(share_geom)
Definition: nmg_mod.c:3829
int idb_major_type
Definition: raytrace.h:192
point_t center
Definition: tgc.c:1822
size_t dgr
Definition: poly.h:49
int rt_tgc_tnurb(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct bn_tol *tol)
Definition: tgc.c:2820
Definition: color.c:49
fastf_t tgc_B
Definition: tgc.c:54
#define TGC_NORM_TOP
Definition: tgc.c:88
vect_t axis_a
Definition: tgc.c:1823
#define RT_G_DEBUG
Definition: raytrace.h:1718
vect_t hit_vpriv
PRIVATE vector for xxx_*()
Definition: raytrace.h:253
struct bu_list * vhead
Definition: raytrace.h:1926
#define RT_ADD_VLIST(hd, pnt, draw)
Definition: raytrace.h:1865
void bn_mat_print(const char *title, const mat_t m)
Definition: mat.c:81
#define RT_CK_DB_INTERNAL(_p)
Definition: raytrace.h:207
int rt_tgc_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: tgc.c:205
#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
#define VEXCHANGE(a, b, tmp)
Definition: tgc.c:92
void rt_tgc_volume(fastf_t *vol, const struct rt_db_internal *ip)
Definition: tgc.c:3246
fastf_t crv_c2
curvature in other direction
Definition: raytrace.h:309
#define BU_PTBL_GET(ptbl, i)
Definition: ptbl.h:108
mat_t tgc_invRtShSc
Definition: tgc.c:63
#define BN_VLIST_LINE_MOVE
Definition: vlist.h:82
const struct rt_functab * idb_meth
for ft_ifree(), etc.
Definition: raytrace.h:194
fastf_t uv_dv
delta in v
Definition: raytrace.h:344
#define V3ARGS(a)
Definition: color.c:56
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
uint8_t * ext_buf
Definition: parse.h:216
static void top()
#define BU_GET(_ptr, _type)
Definition: malloc.h:201
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
void rt_pr_fallback_angle(struct bu_vls *str, const char *prefix, const double angles[5])
void rt_find_fallback_angle(double angles[5], const vect_t vec)
#define TGC_NORM_BOT
Definition: tgc.c:89
struct hit seg_out
OUT information.
Definition: raytrace.h:371
mat_t tgc_ScShR
Definition: tgc.c:62
#define BN_VLIST_LINE_DRAW
Definition: vlist.h:83
Coord * point
Definition: chull3d.cpp:52
void * bu_realloc(void *ptr, size_t siz, const char *str)
#define UNUSED(parameter)
Definition: common.h:239
#define BU_PUT(_ptr, _type)
Definition: malloc.h:215
void rt_tgc_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
Definition: tgc.c:1462
void bn_mat_mul(mat_t o, const mat_t a, const mat_t b)
struct faceuse * nmg_find_fu_of_vu(const struct vertexuse *vu)
Definition: nmg_info.c:304
Support for uniform tolerances.
Definition: tol.h:71
#define BU_LIST_FIRST_MAGIC(hp)
Definition: list.h:416
#define bu_offsetofarray(_t, _a, _d, _i)
Definition: parse.h:65
struct nmgregion * nmg_mrsv(struct model *m)
Definition: nmg_mk.c:306
#define BU_STRUCTPARSE_FUNC_NULL
Definition: parse.h:153
struct fg_node fg
Definition: chull3d.cpp:80
int rt_tgc_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: tgc.c:2013
ustring alpha
vect_t r_dir
Direction of ray (UNIT Length)
Definition: raytrace.h:219
#define ID_REC
Right Elliptical Cylinder [TGC special].
Definition: raytrace.h:465
void bu_ptbl_free(struct bu_ptbl *b)
Definition: ptbl.c:226
struct bn_tol rti_tol
Math tolerances for this model.
Definition: raytrace.h:1765
void nmg_edge_g_cnurb_plinear(struct edgeuse *eu)
Definition: nmg_mk.c:2019
struct vertex * v
Definition: tgc.c:2138
#define RT_CK_DBI(_p)
Definition: raytrace.h:829
#define RT_GET_SEG(p, res)
Definition: raytrace.h:379
fastf_t point_spacing
Definition: raytrace.h:1934
#define ZERO(val)
Definition: units.c:38
void * idb_ptr
Definition: raytrace.h:195
void rt_tgc_print(register const struct soltab *stp)
Definition: tgc.c:514
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
fastf_t tgc_AAdCC
Definition: tgc.c:59
void bu_cv_ntohd(unsigned char *out, const unsigned char *in, size_t count)
const struct bu_structparse rt_tgc_parse[]
Definition: tgc.c:114
fastf_t tgc_BBdDD
Definition: tgc.c:60
const struct rt_functab OBJ[]
Definition: table.c:159
int rt_poly_roots(bn_poly_t *eqn, bn_complex_t roots[], const char *name)
vect_t tgc_V
Definition: tgc.c:51
#define BU_PTBL_END(ptbl)
Definition: ptbl.h:106
axis
Definition: color.c:48
#define TGC_NORM_BODY
Definition: tgc.c:87
void rt_tgc_centroid(point_t *cent, const struct rt_db_internal *ip)
Definition: tgc.c:3324
void bn_vec_ortho(vect_t out, const vect_t in)
#define R
Definition: msr.c:55
int rt_tgc_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
Definition: tgc.c:565
double abs
absolute dist tol
Definition: raytrace.h:180
void * st_specific
-> ID-specific (private) struct
Definition: raytrace.h:435
void rt_nurb_s_eval(const struct face_g_snurb *srf, fastf_t u, fastf_t v, fastf_t *final_value)
Definition: nurb_eval.c:49
void rt_tgc_vshot(struct soltab **stp, register struct xray **rp, struct seg *segp, int n, struct application *ap)
Definition: tgc.c:959
fastf_t uv_du
delta in u
Definition: raytrace.h:343
#define BU_PTBL_BASEADDR(ptbl)
Definition: ptbl.h:104
#define VLARGE
Definition: tgc.c:79
fastf_t crv_c1
curvature in principle dir
Definition: raytrace.h:308
#define A
Definition: msr.c:51
Definition: color.c:51
int dbi_version
PRIVATE: use db_version()
Definition: raytrace.h:824
void rt_pt_sort(register fastf_t *t, int npts)
void plot_ellipse(struct bu_list *vhead, const vect_t t, const vect_t a, const vect_t b, int num_points)
double norm
normal tol
Definition: raytrace.h:182
fastf_t tgc_D
Definition: tgc.c:56
int hit_surfno
solid-specific surface indicator
Definition: raytrace.h:255
vect_t axis_b
Definition: tgc.c:1824
#define Q
Definition: msr.c:54
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
#define ELL_CIRCUMFERENCE(a, b)
Definition: librt_private.h:40
vect_t tan_axb
Definition: tgc.c:2137
#define BU_CK_EXTERNAL(_p)
Definition: parse.h:224
void nmg_pr_fu_briefly(const struct faceuse *fu, char *h)
Definition: nmg_pr.c:359
void nmg_gluefaces(struct faceuse **fulist, int n, const struct bn_tol *tol)
Definition: nmg_mod.c:1408
struct faceuse * nmg_cmface(struct shell *s, struct vertex ***verts, int n)
Definition: nmg_mod.c:979
void nmg_vertexuse_a_cnurb(struct vertexuse *vu, const fastf_t *uvw)
Definition: nmg_mk.c:1757
Complex numbers.
Definition: complex.h:39
const struct rt_functab * st_meth
pointer to per-solid methods
Definition: raytrace.h:428
char tgc_AD_CB
Definition: tgc.c:64
#define RT_PCOEF_TOL
Definition: raytrace.h:171
void rt_ell_16pts(fastf_t *ov, fastf_t *V, fastf_t *A, fastf_t *B)
Definition: ell.c:606
size_t ext_nbytes
Definition: parse.h:210
void ellipse_point_at_radian(point_t result, const vect_t center, const vect_t axis_a, const vect_t axis_b, fastf_t radian)
int rt_tgc_import4(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
Definition: tgc.c:1574
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
#define ID_TGC
Generalized Truncated General Cone.
Definition: raytrace.h:460
#define MAX_RATIO
Definition: tgc.c:80
Definition: vls.h:56
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
#define st_name
Definition: raytrace.h:448
void rt_tgc_free(struct soltab *stp)
Definition: tgc.c:1560
double fastf_t
Definition: defines.h:300
#define RAT
Definition: tgc.c:81
#define VPRINT(a, b)
Definition: raytrace.h:1881
#define BU_LIST_NEXT(structure, hp)
Definition: list.h:316
#define DEBUG_SOLIDS
6 Print prep'ed solids
Definition: raytrace.h:89
void rt_tgc_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
Definition: tgc.c:2055
int rt_tgc_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: tgc.c:1691
int rt_tgc_adaptive_plot(struct rt_db_internal *ip, const struct rt_view_info *info)
Definition: tgc.c:1921
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
point_t st_center
Centroid of solid.
Definition: raytrace.h:432
#define M
Definition: msr.c:52
void nmg_region_a(struct nmgregion *r, const struct bn_tol *tol)
Definition: nmg_mk.c:2557