BRL-CAD
epa.c
Go to the documentation of this file.
1 /* E P A . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1990-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/epa/epa.c
23  *
24  * Intersect a ray with an Elliptical Paraboloid.
25  *
26  * Algorithm -
27  *
28  * Given V, H, R, and B, there is a set of points on this epa
29  *
30  * { (x, y, z) | (x, y, z) is on epa }
31  *
32  * Through a series of Affine Transformations, this set of points will
33  * be transformed into a set of points on an epa located at the origin
34  * with a semi-major axis R1 along the +Y axis, a semi-minor axis R2
35  * along the -X axis, a height H along the -Z axis, and a vertex V at
36  * the origin.
37  *
38  *
39  * { (x', y', z') | (x', y', z') is on epa at origin }
40  *
41  * The transformation from X to X' is accomplished by:
42  *
43  * X' = S(R(X - V))
44  *
45  * where R(X) = (R1/(-|R1|))
46  * (R2/(|R2|)) . X
47  * (H /(-|H |))
48  *
49  * and S(X) = (1/|R1| 0 0)
50  * (0 1/|R2| 0) . X
51  * (0 0 1/|H |)
52  *
53  * To find the intersection of a line with the surface of the epa,
54  * consider the parametric line L:
55  *
56  * L : { P(n) | P + t(n) . D }
57  *
58  * Call W the actual point of intersection between L and the epa. Let
59  * W' be the point of intersection between L' and the unit epa.
60  *
61  * L' : { P'(n) | P' + t(n) . D' }
62  *
63  * W = invR(invS(W')) + V
64  *
65  * Where W' = k D' + P'.
66  *
67  * If Dy' and Dz' are both 0, then there is no hit on the epa; but the
68  * end plates need checking. If there is now only 1 hit point, the
69  * top plate needs to be checked as well.
70  *
71  * Line L' hits the infinitely long epa at W' when
72  *
73  * A * k**2 + B * k + C = 0
74  *
75  * where
76  *
77  * A = Dx'**2 + Dy'**2
78  * B = 2 * (Dx' * Px' + Dy' * Py') - Dz'
79  * C = Px'**2 + Py'**2 - Pz' - 1
80  * b = |Breadth| = 1.0
81  * h = |Height| = 1.0
82  * r = 1.0
83  *
84  * The quadratic formula yields k (which is constant):
85  *
86  * k = [ -B +/- sqrt(B**2 - 4 * A * C)] / (2.0 * A)
87  *
88  * Now, D' = S(R(D))
89  * and P' = S(R(P - V))
90  *
91  * Substituting,
92  *
93  * W = V + invR(invS[ k *(S(R(D))) + S(R(P - V)) ])
94  * = V + invR((k * R(D)) + R(P - V))
95  * = V + k * D + P - V
96  * = k * D + P
97  *
98  * Note that ``k'' is constant, and is the same in the formulations
99  * for both W and W'.
100  *
101  * The hit at ``k'' is a hit on the canonical epa IFF Wz' <= 0.
102  *
103  * NORMALS. Given the point W on the surface of the epa, what is the
104  * vector normal to the tangent plane at that point?
105  *
106  * Map W onto the unit epa, i.e.: W' = S(R(W - V)).
107  *
108  * Plane on unit epa at W' has a normal vector N' where
109  *
110  * N' = <Wx', Wy', -.5>.
111  *
112  * The plane transforms back to the tangent plane at W, and this new
113  * plane (on the original epa) has a normal vector of N, viz:
114  *
115  * N = inverse[ transpose(inverse[ S o R ]) ] (N')
116  *
117  * because if H is perpendicular to plane Q, and matrix M maps from Q
118  * to Q', then inverse[ transpose(M) ] (H) is perpendicular to Q'.
119  * Here, H and Q are in "prime space" with the unit sphere. [Somehow,
120  * the notation here is backwards]. So, the mapping matrix M =
121  * inverse(S o R), because S o R maps from normal space to the unit
122  * sphere.
123  *
124  * N = inverse[ transpose(inverse[ S o R ]) ] (N')
125  * = inverse[ transpose(invR o invS) ] (N')
126  * = inverse[ transpose(invS) o transpose(invR) ] (N')
127  * = inverse[ inverse(S) o R ] (N')
128  * = invR o S (N')
129  *
130  * because inverse(R) = transpose(R), so R = transpose(invR),
131  * and S = transpose(S).
132  *
133  * Note that the normal vector produced above will not have unit
134  * length.
135  *
136  * THE TOP PLATE.
137  *
138  * If Dz' == 0, line L' is parallel to the top plate, so there is no
139  * hit on the top plate. Otherwise, rays intersect the top plate with
140  * k = (0 - Pz')/Dz'. The solution is within the top plate IFF Wx'**2
141  * + Wy'**2 <= 1.
142  *
143  * The normal for a hit on the top plate is -Hunit.
144  *
145  */
146 
147 #include "common.h"
148 
149 #include <stddef.h>
150 #include <string.h>
151 #include <math.h>
152 #include "bio.h"
153 
154 #include "bu/cv.h"
155 #include "vmath.h"
156 #include "db.h"
157 #include "nmg.h"
158 #include "rtgeom.h"
159 #include "raytrace.h"
160 
161 #include "../../librt_private.h"
162 
163 static int epa_is_valid(struct rt_epa_internal *epa);
164 
165 struct epa_specific {
166  point_t epa_V; /* vector to epa origin */
167  vect_t epa_Hunit; /* unit H vector */
168  vect_t epa_Aunit; /* unit vector along semi-major axis */
169  vect_t epa_Bunit; /* unit vector, A x H */
170  fastf_t epa_h; /* |H| */
171  fastf_t epa_inv_r1sq; /* 1/(r1 * r1) */
172  fastf_t epa_inv_r2sq; /* 1/(r2 * r2) */
173  mat_t epa_SoR; /* Scale(Rot(vect)) */
174  mat_t epa_invRoS; /* invRot(Scale(vect)) */
175 };
176 
177 
178 const struct bu_structparse rt_epa_parse[] = {
179  { "%f", 3, "V", bu_offsetofarray(struct rt_epa_internal, epa_V, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
180  { "%f", 3, "H", bu_offsetofarray(struct rt_epa_internal, epa_H, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
181  { "%f", 3, "A", bu_offsetofarray(struct rt_epa_internal, epa_Au, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
182  { "%f", 1, "r_1", bu_offsetof(struct rt_epa_internal, epa_r1), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
183  { "%f", 1, "r_2", bu_offsetof(struct rt_epa_internal, epa_r2), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
184  { {'\0', '\0', '\0', '\0'}, 0, (char *)NULL, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL }
185 };
186 
187 /**
188  * Create a bounding RPP for an epa
189  */
190 int
191 rt_epa_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol)) {
192  struct rt_epa_internal *xip;
193  vect_t epa_A, epa_B, epa_An, epa_Bn, epa_H;
194  vect_t pt1, pt2, pt3, pt4, pt5, pt6, pt7, pt8;
195  RT_CK_DB_INTERNAL(ip);
196  xip = (struct rt_epa_internal *)ip->idb_ptr;
197  RT_EPA_CK_MAGIC(xip);
198 
199  VMOVE(epa_H, xip->epa_H);
200  VUNITIZE(epa_H);
201  VMOVE(epa_A, xip->epa_Au);
202  VCROSS(epa_B, epa_A, epa_H);
203 
204  VSETALL((*min), INFINITY);
205  VSETALL((*max), -INFINITY);
206 
207  VSCALE(epa_A, epa_A, xip->epa_r1);
208  VSCALE(epa_B, epa_B, xip->epa_r2);
209  VREVERSE(epa_An, epa_A);
210  VREVERSE(epa_Bn, epa_B);
211 
212  VADD3(pt1, xip->epa_V, epa_A, epa_B);
213  VADD3(pt2, xip->epa_V, epa_A, epa_Bn);
214  VADD3(pt3, xip->epa_V, epa_An, epa_B);
215  VADD3(pt4, xip->epa_V, epa_An, epa_Bn);
216  VADD4(pt5, xip->epa_V, epa_A, epa_B, xip->epa_H);
217  VADD4(pt6, xip->epa_V, epa_A, epa_Bn, xip->epa_H);
218  VADD4(pt7, xip->epa_V, epa_An, epa_B, xip->epa_H);
219  VADD4(pt8, xip->epa_V, epa_An, epa_Bn, xip->epa_H);
220 
221  /* Find the RPP of the rotated axis-aligned epa bbox - that is,
222  * the bounding box the given epa would have if its height
223  * vector were in the positive Z direction. This does not give
224  * us an optimal bbox except in the case where the epa is
225  * actually axis aligned to start with, but it's usually
226  * at least a bit better than the bounding sphere RPP. */
227  VMINMAX((*min), (*max), pt1);
228  VMINMAX((*min), (*max), pt2);
229  VMINMAX((*min), (*max), pt3);
230  VMINMAX((*min), (*max), pt4);
231  VMINMAX((*min), (*max), pt5);
232  VMINMAX((*min), (*max), pt6);
233  VMINMAX((*min), (*max), pt7);
234  VMINMAX((*min), (*max), pt8);
235  return 0;
236 }
237 
238 
239 /**
240  * Given a pointer to a GED database record, and a transformation
241  * matrix, determine if this is a valid EPA, and if so, precompute
242  * various terms of the formula.
243  *
244  * Returns -
245  * 0 EPA is OK
246  * !0 Error in description
247  *
248  * Implicit return -
249  * A struct epa_specific is created, and its address is stored in
250  * stp->st_specific for use by epa_shot().
251  */
252 int
253 rt_epa_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
254 {
255  struct rt_epa_internal *xip;
256  struct epa_specific *epa;
257 
258  fastf_t magsq_h;
259  fastf_t /* mag_a, */ mag_h;
260  fastf_t r1, r2;
261  mat_t R;
262  mat_t Rinv;
263  mat_t S;
264 
265  RT_CK_DB_INTERNAL(ip);
266 
267  xip = (struct rt_epa_internal *)ip->idb_ptr;
268  if (!epa_is_valid(xip)) {
269  return 1;
270  }
271 
272  /* compute |A| |H| */
273  /* mag_a = sqrt(MAGSQ(xip->epa_Au)); */
274  mag_h = sqrt(magsq_h = MAGSQ(xip->epa_H));
275  r1 = xip->epa_r1;
276  r2 = xip->epa_r2;
277 
278 
279  /*
280  * EPA is ok
281  */
282  stp->st_id = ID_EPA; /* set soltab ID */
283  stp->st_meth = &OBJ[ID_EPA];
284 
285  BU_GET(epa, struct epa_specific);
286  stp->st_specific = (void *)epa;
287 
288  epa->epa_h = mag_h;
289  epa->epa_inv_r1sq = 1 / (r1 * r1);
290  epa->epa_inv_r2sq = 1 / (r2 * r2);
291 
292  /* make unit vectors in A, H, and BxH directions */
293  VMOVE(epa->epa_Hunit, xip->epa_H);
294  VUNITIZE(epa->epa_Hunit);
295  VMOVE(epa->epa_Aunit, xip->epa_Au);
296  VCROSS(epa->epa_Bunit, epa->epa_Aunit, epa->epa_Hunit);
297 
298  VMOVE(epa->epa_V, xip->epa_V);
299 
300  /* Compute R and Rinv matrices */
301  MAT_IDN(R);
302  VREVERSE(&R[0], epa->epa_Bunit);
303  VMOVE(&R[4], epa->epa_Aunit);
304  VREVERSE(&R[8], epa->epa_Hunit);
305  bn_mat_trn(Rinv, R); /* inv of rot mat is trn */
306 
307  /* Compute S */
308  MAT_IDN(S);
309  S[ 0] = 1.0/r2;
310  S[ 5] = 1.0/r1;
311  S[10] = 1.0/mag_h;
312 
313  /* Compute SoR and invRoS */
314  bn_mat_mul(epa->epa_SoR, S, R);
315  bn_mat_mul(epa->epa_invRoS, Rinv, S);
316 
317  /* Compute bounding sphere */
318  /* bounding sphere center */
319  VJOIN1(stp->st_center, epa->epa_V, mag_h / 2.0, epa->epa_Hunit);
320  /* bounding radius */
321  stp->st_bradius = sqrt(0.25*magsq_h + r2*r2 + r1*r1);
322  /* approximate bounding radius */
323  stp->st_aradius = stp->st_bradius;
324 
325  /* Calculate bounding box (RPP) */
326  if (rt_epa_bbox(ip, &(stp->st_min), &(stp->st_max), &rtip->rti_tol)) return 1;
327 
328  return 0; /* OK */
329 }
330 
331 
332 void
333 rt_epa_print(const struct soltab *stp)
334 {
335  const struct epa_specific *epa =
336  (struct epa_specific *)stp->st_specific;
337 
338  VPRINT("V", epa->epa_V);
339  VPRINT("Hunit", epa->epa_Hunit);
340  VPRINT("Aunit", epa->epa_Aunit);
341  VPRINT("Bunit", epa->epa_Bunit);
342  bn_mat_print("S o R", epa->epa_SoR);
343  bn_mat_print("invR o S", epa->epa_invRoS);
344 }
345 
346 
347 /* hit_surfno is set to one of these */
348 #define EPA_NORM_BODY (1) /* compute normal */
349 #define EPA_NORM_TOP (2) /* copy epa_N */
350 
351 
352 /**
353  * Intersect a ray with a epa. If an intersection occurs, a struct
354  * seg will be acquired and filled in.
355  *
356  * Returns -
357  * 0 MISS
358  * >0 HIT
359  */
360 int
361 rt_epa_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
362 {
363  struct epa_specific *epa =
364  (struct epa_specific *)stp->st_specific;
365  vect_t dprime; /* D' */
366  vect_t pprime; /* P' */
367  fastf_t k1, k2; /* distance constants of solution */
368  vect_t xlated; /* translated vector */
369  struct hit hits[3]; /* 2 potential hit points */
370  struct hit *hitp; /* pointer to hit point */
371 
372  hitp = &hits[0];
373 
374  /* out, Mat, vect */
375  MAT4X3VEC(dprime, epa->epa_SoR, rp->r_dir);
376  VSUB2(xlated, rp->r_pt, epa->epa_V);
377  MAT4X3VEC(pprime, epa->epa_SoR, xlated);
378 
379  /* Find roots of the equation, using formula for quadratic */
380  {
381  fastf_t a, b, c; /* coeffs of polynomial */
382  fastf_t disc; /* disc of radical */
383 
384  a = dprime[X] * dprime[X] + dprime[Y] * dprime[Y];
385  b = 2*(dprime[X] * pprime[X] + dprime[Y] * pprime[Y])
386  - dprime[Z];
387  c = pprime[X] * pprime[X]
388  + pprime[Y] * pprime[Y] - pprime[Z] - 1.0;
389  if (!NEAR_ZERO(a, RT_PCOEF_TOL)) {
390  disc = b*b - 4 * a * c;
391  if (disc <= 0)
392  goto check_plates;
393  disc = sqrt(disc);
394 
395  k1 = (-b + disc) / (2.0 * a);
396  k2 = (-b - disc) / (2.0 * a);
397 
398  /* k1 and k2 are potential solutions to intersection with
399  * side. See if they fall in range.
400  */
401  VJOIN1(hitp->hit_vpriv, pprime, k1, dprime); /* hit' */
402  if (hitp->hit_vpriv[Z] <= 0.0) {
403  hitp->hit_magic = RT_HIT_MAGIC;
404  hitp->hit_dist = k1;
405  hitp->hit_surfno = EPA_NORM_BODY; /* compute N */
406  hitp++;
407  }
408 
409  VJOIN1(hitp->hit_vpriv, pprime, k2, dprime); /* hit' */
410  if (hitp->hit_vpriv[Z] <= 0.0) {
411  hitp->hit_magic = RT_HIT_MAGIC;
412  hitp->hit_dist = k2;
413  hitp->hit_surfno = EPA_NORM_BODY; /* compute N */
414  hitp++;
415  }
416  } else if (!NEAR_ZERO(b, RT_PCOEF_TOL)) {
417  k1 = -c/b;
418  VJOIN1(hitp->hit_vpriv, pprime, k1, dprime); /* hit' */
419  if (hitp->hit_vpriv[Z] <= 0.0) {
420  hitp->hit_magic = RT_HIT_MAGIC;
421  hitp->hit_dist = k1;
422  hitp->hit_surfno = EPA_NORM_BODY; /* compute N */
423  hitp++;
424  }
425  }
426  }
427 
428 
429  /*
430  * Check for hitting the top plate.
431  */
432  check_plates:
433  /* check top plate */
434  if (hitp == &hits[1] && !ZERO(dprime[Z])) {
435  /* 1 hit so far, this is worthwhile */
436  k1 = -pprime[Z] / dprime[Z]; /* top plate */
437 
438  VJOIN1(hitp->hit_vpriv, pprime, k1, dprime); /* hit' */
439  if (hitp->hit_vpriv[X] * hitp->hit_vpriv[X] +
440  hitp->hit_vpriv[Y] * hitp->hit_vpriv[Y] <= 1.0) {
441  hitp->hit_magic = RT_HIT_MAGIC;
442  hitp->hit_dist = k1;
443  hitp->hit_surfno = EPA_NORM_TOP; /* -H */
444  hitp++;
445  }
446  }
447 
448  if (hitp != &hits[2])
449  return 0; /* MISS */
450 
451  if (hits[0].hit_dist < hits[1].hit_dist) {
452  /* entry is [0], exit is [1] */
453  struct seg *segp;
454 
455  RT_GET_SEG(segp, ap->a_resource);
456  segp->seg_stp = stp;
457  segp->seg_in = hits[0]; /* struct copy */
458  segp->seg_out = hits[1]; /* struct copy */
459  BU_LIST_INSERT(&(seghead->l), &(segp->l));
460  } else {
461  /* entry is [1], exit is [0] */
462  struct seg *segp;
463 
464  RT_GET_SEG(segp, ap->a_resource);
465  segp->seg_stp = stp;
466  segp->seg_in = hits[1]; /* struct copy */
467  segp->seg_out = hits[0]; /* struct copy */
468  BU_LIST_INSERT(&(seghead->l), &(segp->l));
469  }
470  return 2; /* HIT */
471 }
472 
473 
474 /**
475  * Given ONE ray distance, return the normal and entry/exit point.
476  */
477 void
478 rt_epa_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
479 {
480  fastf_t scale;
481  vect_t can_normal; /* normal to canonical epa */
482  struct epa_specific *epa =
483  (struct epa_specific *)stp->st_specific;
484 
485  VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
486  switch (hitp->hit_surfno) {
487  case EPA_NORM_BODY:
488  VSET(can_normal,
489  hitp->hit_vpriv[X],
490  hitp->hit_vpriv[Y],
491  -0.5);
492  MAT4X3VEC(hitp->hit_normal, epa->epa_invRoS, can_normal);
493  scale = 1.0 / MAGNITUDE(hitp->hit_normal);
494  VSCALE(hitp->hit_normal, hitp->hit_normal, scale);
495 
496  /* tuck away this scale for the curvature routine */
497  hitp->hit_vpriv[X] = scale;
498  break;
499  case EPA_NORM_TOP:
500  VREVERSE(hitp->hit_normal, epa->epa_Hunit);
501  break;
502  default:
503  bu_log("rt_epa_norm: surfno=%d bad\n", hitp->hit_surfno);
504  break;
505  }
506 }
507 
508 
509 /**
510  * Return the curvature of the epa.
511  */
512 void
513 rt_epa_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
514 {
515  fastf_t a, b, c, scale;
516  mat_t M1, M2;
517  struct epa_specific *epa =
518  (struct epa_specific *)stp->st_specific;
519  vect_t u, v; /* basis vectors (with normal) */
520  vect_t vec1, vec2; /* eigen vectors */
521  vect_t tmp;
522 
523  switch (hitp->hit_surfno) {
524  case EPA_NORM_BODY:
525  /* choose a tangent plane coordinate system (u, v, normal)
526  * form a right-handed triple
527  */
528  bn_vec_ortho(u, hitp->hit_normal);
529  VCROSS(v, hitp->hit_normal, u);
530 
531  /* get the saved away scale factor */
532  scale = - hitp->hit_vpriv[X];
533 
534  MAT_IDN(M1);
535  M1[10] = 0; /* M1[3, 3] = 0 */
536  /* M1 = invR * S * M1 * S * R */
537  bn_mat_mul(M2, epa->epa_invRoS, M1);
538  bn_mat_mul(M1, M2, epa->epa_SoR);
539 
540  /* find the second fundamental form */
541  MAT4X3VEC(tmp, M1, u);
542  a = VDOT(u, tmp) * scale;
543  b = VDOT(v, tmp) * scale;
544  MAT4X3VEC(tmp, M1, v);
545  c = VDOT(v, tmp) * scale;
546 
547  bn_eigen2x2(&cvp->crv_c1, &cvp->crv_c2, vec1, vec2, a, b, c);
548  VCOMB2(cvp->crv_pdir, vec1[X], u, vec1[Y], v);
549  VUNITIZE(cvp->crv_pdir);
550  break;
551  case EPA_NORM_TOP:
552  cvp->crv_c1 = cvp->crv_c2 = 0;
553  /* any tangent direction */
554  bn_vec_ortho(cvp->crv_pdir, hitp->hit_normal);
555  break;
556  }
557 }
558 
559 
560 /**
561  * For a hit on the surface of an epa, return the (u, v) coordinates
562  * of the hit point, 0 <= u, v <= 1.
563  *
564  * u = azimuth
565  * v = elevation
566  */
567 void
568 rt_epa_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
569 {
570  struct epa_specific *epa =
571  (struct epa_specific *)stp->st_specific;
572 
573  vect_t work;
574  vect_t pprime;
575  fastf_t len;
576 
577  if (ap) RT_CK_APPLICATION(ap);
578 
579  /* hit_point is on surface; project back to unit epa, creating a
580  * vector from vertex to hit point.
581  */
582  VSUB2(work, hitp->hit_point, epa->epa_V);
583  MAT4X3VEC(pprime, epa->epa_SoR, work);
584 
585  switch (hitp->hit_surfno) {
586  case EPA_NORM_BODY:
587  /* top plate, polar coords */
588  if (ZERO(pprime[Z] + 1.0)) { /* i.e., == -1.0 */
589  /* bottom pt of body */
590  uvp->uv_u = 0;
591  } else {
592  len = sqrt(pprime[X]*pprime[X] + pprime[Y]*pprime[Y]);
593  uvp->uv_u = acos(pprime[X]/len) * M_1_2PI;
594  }
595  uvp->uv_v = -pprime[Z];
596  break;
597  case EPA_NORM_TOP:
598  /* top plate, polar coords */
599  len = sqrt(pprime[X]*pprime[X] + pprime[Y]*pprime[Y]);
600  uvp->uv_u = acos(pprime[X]/len) * M_1_2PI;
601  uvp->uv_v = 1.0 - len;
602  break;
603  }
604  /* Handle other half of acos() domain */
605  if (pprime[Y] < 0)
606  uvp->uv_u = 1.0 - uvp->uv_u;
607 
608  /* uv_du should be relative to rotation, uv_dv relative to height */
609  uvp->uv_du = uvp->uv_dv = 0;
610 }
611 
612 
613 void
614 rt_epa_free(struct soltab *stp)
615 {
616  struct epa_specific *epa =
617  (struct epa_specific *)stp->st_specific;
618 
619  BU_PUT(epa, struct epa_specific);
620 }
621 
622 
623 /* A canonical parabola in the Y-Z plane has equation z = y^2 / 4p, and opens
624  * toward positive z with vertex at the origin.
625  *
626  * The contour of an epa in the plane H-R (where R is one of the epa axes A or
627  * B) is a parabola with vertex at H, opening toward -H. We can transform this
628  * parabola to get an equivalent canonical parabola in the Y-Z plane, opening
629  * toward positive Z (-H) with vertex at the origin (H).
630  *
631  * This parabola passes through the point (r, |H|) (where r = |A| or |B|). If
632  * we plug the point (r, |H|) into our canonical equation, we see how p relates
633  * to r and |H|:
634  *
635  * |H| = r^2 / 4p
636  * p = (r^2) / (4|H|)
637  */
638 static fastf_t
639 epa_parabola_p(fastf_t r, fastf_t mag_h)
640 {
641  return (r * r) / (4 * mag_h);
642 }
643 
644 /* The contour of an epa in the plane H-R (where R is one of the epa axes A or
645  * B) is a parabola with vertex at H, opening toward -H. We can transform this
646  * parabola to get an equivalent parabola in the Y-Z plane, opening toward
647  * positive Z (-H) with vertex at (0, -|H|).
648  *
649  * The part of this parabola that passes between (0, -|H|) and (r, 0) is
650  * approximated by num_points points (including (0, -|H|) and (r, 0)).
651  *
652  * The constructed point list is returned (NULL returned on error). Because the
653  * above transformation puts the epa vertex at the origin and the parabola
654  * vertex at (0, -|H|), multiplying the z values by -1 gives corresponding
655  * distances along the epa height vector H.
656  */
657 static struct rt_pt_node *
658 epa_parabolic_curve(fastf_t mag_h, fastf_t r, int num_points)
659 {
660  int count;
661  struct rt_pt_node *curve;
662 
663  BU_ALLOC(curve, struct rt_pt_node);
664  BU_ALLOC(curve->next, struct rt_pt_node);
665 
666  curve->next->next = NULL;
667  VSET(curve->p, 0, 0, -mag_h);
668  VSET(curve->next->p, 0, r, 0);
669 
670  count = approximate_parabolic_curve(curve, epa_parabola_p(r, mag_h), num_points - 2);
671 
672  if (count != (num_points - 2)) {
673  return NULL;
674  }
675 
676  return curve;
677 }
678 
679 /* The contour of an epa in the plane H-R (where R is one of the epa axes A or
680  * B) is a parabola with vertex at H, opening toward -H. We can transform this
681  * parabola into an equivalent one in the Y-Z plane which has vertext at (0, |H|)
682  * and opens toward -Z.
683  *
684  * The equation for this parabola is a variant of the equation for a canonical
685  * parabola in the Y-Z plane (z = y^2 / 4p):
686  * z = |H| - y^2 / 4p
687  *
688  * Solving this equation for y yields:
689  * y = sqrt(4p * (|H| - z))
690  *
691  * Substituting p = r^2 / 4|H| (see above comment):
692  * y = sqrt(r^2 * (|H| - z) / |H|)
693  * = r * sqrt((|H| - z) / |H|)
694  * = r * sqrt(1 - z / |H|)
695  */
696 static fastf_t
697 epa_parabola_y(fastf_t r, fastf_t mag_H, fastf_t z)
698 {
699  return r * sqrt(1 - z / mag_H);
700 }
701 
702 /* Plot the elliptical cross section of the given epa at distance h along the
703  * epa height vector (h >= 0, h <= |H|) consisting of num_points points.
704  */
705 static void
706 epa_plot_ellipse(
707  struct bu_list *vhead,
708  struct rt_epa_internal *epa,
709  fastf_t h,
710  fastf_t num_points)
711 {
712  fastf_t mag_H;
713  vect_t V, Hu, Au, Bu, A, B, cross_section_plane;
714 
715  VMOVE(V, epa->epa_V);
716 
717  mag_H = MAGNITUDE(epa->epa_H);
718  VSCALE(Hu, epa->epa_H, 1.0 / mag_H);
719 
720  VMOVE(Au, epa->epa_Au);
721  VCROSS(Bu, Au, Hu);
722 
723  /* calculate semi-major and semi-minor axis for the elliptical
724  * cross-section at distance h along H
725  */
726  VSCALE(A, Au, epa_parabola_y(epa->epa_r1, mag_H, h));
727  VSCALE(B, Bu, epa_parabola_y(epa->epa_r2, mag_H, h));
728  VJOIN1(cross_section_plane, V, h, Hu);
729 
730  plot_ellipse(vhead, cross_section_plane, A, B, num_points);
731 }
732 
733 static void
734 epa_plot_parabola(
735  struct bu_list *vhead,
736  struct rt_epa_internal *epa,
737  struct rt_pt_node *pts,
738  vect_t Ru,
739  fastf_t r)
740 {
741  point_t p;
742  vect_t epa_V, Hu;
743  fastf_t mag_H, z;
744  struct rt_pt_node *node;
745 
746  VMOVE(epa_V, epa->epa_V);
747  mag_H = MAGNITUDE(epa->epa_H);
748  VSCALE(Hu, epa->epa_H, 1.0 / mag_H);
749 
750  z = pts->p[Z];
751  VJOIN2(p, epa_V, epa_parabola_y(r, mag_H, -z), Ru, -z, Hu);
752  RT_ADD_VLIST(vhead, p, BN_VLIST_LINE_MOVE);
753 
754  node = pts->next;
755  while (node != NULL) {
756  z = node->p[Z];
757  VJOIN2(p, epa_V, epa_parabola_y(r, mag_H, -z), Ru, -z, Hu);
758 
759  RT_ADD_VLIST(vhead, p, BN_VLIST_LINE_DRAW);
760 
761  node = node->next;
762  }
763 }
764 
765 static int
766 epa_curve_points(
767  const struct rt_epa_internal *epa,
768  const struct rt_view_info *info)
769 {
770  fastf_t avg_r, approx_curve_len;
771  point_t p1, p2;
772 
773  avg_r = (epa->epa_r1 + epa->epa_r2) / 2.0;
774 
775  VADD2(p1, epa->epa_V, epa->epa_H);
776  VJOIN1(p2, epa->epa_V, avg_r, epa->epa_Au);
777 
778  approx_curve_len = 2.0 * DIST_PT_PT(p1, p2);
779 
780  return approx_curve_len / info->point_spacing;
781 }
782 
783 static int
784 epa_ellipse_points(
785  struct rt_epa_internal *epa,
786  const struct rt_view_info *info)
787 {
788  fastf_t avg_radius, avg_circumference;
789 
790  avg_radius = (epa->epa_r1 + epa->epa_r2) / 2.0;
791  avg_circumference = M_2PI * avg_radius;
792 
793  return avg_circumference / info->point_spacing;
794 }
795 
796 int
797 rt_epa_adaptive_plot(struct rt_db_internal *ip, const struct rt_view_info *info)
798 {
799  vect_t epa_H, Hu, Au, Bu;
800  fastf_t mag_H, z, z_step, r1, r2;
801  int i, num_curve_points, num_ellipse_points, num_curves;
802  struct rt_epa_internal *epa;
803  struct rt_pt_node *pts_r1, *pts_r2, *node, *node1, *node2;
804 
805  BU_CK_LIST_HEAD(info->vhead);
806  RT_CK_DB_INTERNAL(ip);
807 
808  epa = (struct rt_epa_internal *)ip->idb_ptr;
809  if (!epa_is_valid(epa)) {
810  return -2;
811  }
812 
813  num_curve_points = epa_curve_points(epa, info);
814 
815  if (num_curve_points < 3) {
816  num_curve_points = 3;
817  }
818 
819  num_ellipse_points = epa_ellipse_points(epa, info);
820 
821  if (num_ellipse_points < 6) {
822  num_ellipse_points = 6;
823  }
824 
825 
826  VMOVE(epa_H, epa->epa_H);
827 
828  mag_H = MAGNITUDE(epa_H);
829  VSCALE(Hu, epa->epa_H, 1.0 / mag_H);
830 
831  VMOVE(Au, epa->epa_Au);
832  VCROSS(Bu, Au, Hu);
833 
834  r1 = epa->epa_r1;
835  r2 = epa->epa_r2;
836 
837  pts_r1 = epa_parabolic_curve(mag_H, r1, num_curve_points);
838  pts_r2 = epa_parabolic_curve(mag_H, r2, num_curve_points);
839 
840  if (pts_r1 == NULL || pts_r2 == NULL) {
841  return -1;
842  }
843 
844  num_curves = mag_H / info->curve_spacing;
845  if (num_curves < 2) {
846  num_curves = 2;
847  }
848 
849  z_step = mag_H / num_curves;
850  z = 0.0;
851  for (i = 0; i < num_curves; ++i) {
852  epa_plot_ellipse(info->vhead, epa, z, num_ellipse_points);
853 
854  z += z_step;
855  }
856 
857  epa_plot_parabola(info->vhead, epa, pts_r1, Au, r1);
858  epa_plot_parabola(info->vhead, epa, pts_r1, Au, -r1);
859  epa_plot_parabola(info->vhead, epa, pts_r1, Bu, r2);
860  epa_plot_parabola(info->vhead, epa, pts_r1, Bu, -r2);
861 
862  node1 = pts_r1;
863  node2 = pts_r2;
864  for (i = 0; i < num_curve_points; ++i) {
865  node = node1;
866  node1 = node1->next;
867  bu_free(node, "rt_pt_node");
868 
869  node = node2;
870  node2 = node2->next;
871  bu_free(node, "rt_pt_node");
872  }
873 
874  return 0;
875 }
876 
877 int
878 rt_epa_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *UNUSED(tol), const struct rt_view_info *UNUSED(info))
879 {
880  fastf_t dtol, mag_h, ntol, r1, r2;
881  fastf_t **ellipses, theta_new, theta_prev;
882  int *pts_dbl, i, j, nseg;
883  int jj, na, nb, nell, recalc_b;
884  mat_t R;
885  mat_t invR;
886  struct rt_epa_internal *xip;
887  point_t p1;
888  struct rt_pt_node *pos_a, *pos_b, *pts_a, *pts_b;
889  vect_t A, Au, B, Bu, Hu, V, Work;
890 
891  BU_CK_LIST_HEAD(vhead);
892  RT_CK_DB_INTERNAL(ip);
893 
894  xip = (struct rt_epa_internal *)ip->idb_ptr;
895  if (!epa_is_valid(xip)) {
896  return -2;
897  }
898 
899  /* compute |H| */
900  mag_h = MAGNITUDE(xip->epa_H);
901  r1 = xip->epa_r1;
902  r2 = xip->epa_r2;
903 
904  /* make unit vectors in A, H, and BxH directions */
905  VMOVE(Hu, xip->epa_H);
906  VUNITIZE(Hu);
907  VMOVE(Au, xip->epa_Au);
908  VCROSS(Bu, Au, Hu);
909 
910  /* Compute R and Rinv matrices */
911  MAT_IDN(R);
912  VREVERSE(&R[0], Bu);
913  VMOVE(&R[4], Au);
914  VREVERSE(&R[8], Hu);
915  bn_mat_trn(invR, R); /* inv of rot mat is trn */
916 
917  dtol = primitive_get_absolute_tolerance(ttol, 2.0 * r2);
918 
919  /* To ensure normal tolerance, remain below this angle */
920  if (ttol->norm > 0.0)
921  ntol = ttol->norm;
922  else
923  /* tolerate everything */
924  ntol = M_PI;
925 
926  /*
927  * build epa from 2 parabolas
928  */
929 
930  /* approximate positive half of parabola along semi-minor axis */
931  BU_ALLOC(pts_b, struct rt_pt_node);
932  BU_ALLOC(pts_b->next, struct rt_pt_node);
933 
934  pts_b->next->next = NULL;
935  VSET(pts_b->p, 0, 0, -mag_h);
936  VSET(pts_b->next->p, 0, r2, 0);
937  /* 2 endpoints in 1st approximation */
938  nb = 2;
939 
940  /* recursively break segment 'til within error tolerances */
941  nb += rt_mk_parabola(pts_b, r2, mag_h, dtol, ntol);
942  nell = nb - 1; /* # of ellipses needed */
943 
944  /* construct positive half of parabola along semi-major axis of
945  * epa using same z coords as parab along semi-minor axis
946  */
947  BU_ALLOC(pts_a, struct rt_pt_node);
948  VMOVE(pts_a->p, pts_b->p); /* 1st pt is the apex */
949  pts_a->next = NULL;
950  pos_b = pts_b->next;
951  pos_a = pts_a;
952  while (pos_b) {
953  /* copy node from b_parabola to a_parabola */
954  BU_ALLOC(pos_a->next, struct rt_pt_node);
955  pos_a = pos_a->next;
956  pos_a->p[Z] = pos_b->p[Z];
957  /* at given z, find y on parabola */
958  pos_a->p[Y] = r1*sqrt(pos_a->p[Z] / mag_h + 1);
959  pos_a->p[X] = 0;
960  pos_b = pos_b->next;
961  }
962  pos_a->next = NULL;
963 
964  /* see if any segments need further breaking up */
965  recalc_b = 0;
966  na = 0;
967  pos_a = pts_a;
968  while (pos_a->next) {
969  na = rt_mk_parabola(pos_a, r1, mag_h, dtol, ntol);
970  if (na != 0) {
971  recalc_b = 1;
972  nell += na;
973  }
974  pos_a = pos_a->next;
975  }
976  /* if any were broken, recalculate parabola 'a' */
977  if (recalc_b) {
978  /* free mem for old approximation of parabola */
979  pos_b = pts_b;
980  while (pos_b) {
981  struct rt_pt_node *next;
982 
983  /* get next node before freeing */
984  next = pos_b->next;
985  bu_free((char *)pos_b, "rt_pt_node");
986  pos_b = next;
987  }
988  /* construct parabola along semi-major axis of epa using same
989  * z coords as parab along semi-minor axis
990  */
991  BU_ALLOC(pts_b, struct rt_pt_node);
992  pts_b->p[Z] = pts_a->p[Z];
993  pts_b->next = NULL;
994  pos_a = pts_a->next;
995  pos_b = pts_b;
996  while (pos_a) {
997  /* copy node from a_parabola to b_parabola */
998  BU_ALLOC(pos_b->next, struct rt_pt_node);
999  pos_b = pos_b->next;
1000  pos_b->p[Z] = pos_a->p[Z];
1001  /* at given z, find y on parabola */
1002  pos_b->p[Y] = r2*sqrt(pos_b->p[Z] / mag_h + 1);
1003  pos_b->p[X] = 0;
1004  pos_a = pos_a->next;
1005  }
1006  pos_b->next = NULL;
1007  }
1008 
1009  /* make array of ptrs to epa ellipses */
1010  ellipses = (fastf_t **)bu_malloc(nell * sizeof(fastf_t *), "fastf_t ell[]");
1011  /* keep track of whether pts in each ellipse are doubled or not */
1012  pts_dbl = (int *)bu_malloc(nell * sizeof(int), "dbl ints");
1013 
1014  /* make ellipses at each z level */
1015  i = 0;
1016  nseg = 0;
1017  theta_prev = M_2PI;
1018  pos_a = pts_a->next; /* skip over apex of epa */
1019  pos_b = pts_b->next;
1020  while (pos_a) {
1021  VSCALE(A, Au, pos_a->p[Y]); /* semimajor axis */
1022  VSCALE(B, Bu, pos_b->p[Y]); /* semiminor axis */
1023  VJOIN1(V, xip->epa_V, -pos_a->p[Z], Hu);
1024 
1025  VSET(p1, 0., pos_b->p[Y], 0.);
1026  theta_new = ell_angle(p1, pos_a->p[Y], pos_b->p[Y], dtol, ntol);
1027  if (nseg == 0) {
1028  nseg = (int)(M_2PI / theta_new) + 1;
1029  pts_dbl[i] = 0;
1030  } else if (theta_new < theta_prev) {
1031  nseg *= 2;
1032  pts_dbl[i] = 1;
1033  } else
1034  pts_dbl[i] = 0;
1035  theta_prev = theta_new;
1036 
1037  ellipses[i] = (fastf_t *)bu_malloc(3*(nseg+1)*sizeof(fastf_t),
1038  "pts ell");
1039  rt_ell(ellipses[i], V, A, B, nseg);
1040 
1041  i++;
1042  pos_a = pos_a->next;
1043  pos_b = pos_b->next;
1044  }
1045  /* Draw the top ellipse */
1046  RT_ADD_VLIST(vhead,
1047  &ellipses[nell-1][(nseg-1)*ELEMENTS_PER_VECT],
1049  for (i = 0; i < nseg; i++) {
1050  RT_ADD_VLIST(vhead,
1051  &ellipses[nell-1][i*ELEMENTS_PER_VECT],
1053  }
1054 
1055  /* connect ellipses */
1056  for (i = nell-2; i >= 0; i--) {
1057  /* skip top ellipse */
1058  int bottom, top;
1059 
1060  top = i + 1;
1061  bottom = i;
1062  if (pts_dbl[top])
1063  nseg /= 2; /* # segs in 'bottom' ellipse */
1064 
1065  /* Draw the current ellipse */
1066  RT_ADD_VLIST(vhead,
1067  &ellipses[bottom][(nseg-1)*ELEMENTS_PER_VECT],
1069  for (j = 0; j < nseg; j++) {
1070  RT_ADD_VLIST(vhead,
1071  &ellipses[bottom][j*ELEMENTS_PER_VECT],
1073  }
1074 
1075  /* make connections between ellipses */
1076  for (j = 0; j < nseg; j++) {
1077  if (pts_dbl[top])
1078  jj = j + j; /* top ellipse index */
1079  else
1080  jj = j;
1081  RT_ADD_VLIST(vhead,
1082  &ellipses[bottom][j*ELEMENTS_PER_VECT],
1084  RT_ADD_VLIST(vhead,
1085  &ellipses[top][jj*ELEMENTS_PER_VECT],
1087  }
1088  }
1089 
1090  VADD2(Work, xip->epa_V, xip->epa_H);
1091  for (i = 0; i < nseg; i++) {
1092  /* Draw connector */
1093  RT_ADD_VLIST(vhead, Work, BN_VLIST_LINE_MOVE);
1094  RT_ADD_VLIST(vhead,
1095  &ellipses[0][i*ELEMENTS_PER_VECT],
1097  }
1098 
1099  /* free mem */
1100  for (i = 0; i < nell; i++) {
1101  bu_free((char *)ellipses[i], "pts ell");
1102  }
1103  bu_free((char *)ellipses, "fastf_t ell[]");
1104  bu_free((char *)pts_dbl, "dbl ints");
1105 
1106  return 0;
1107 }
1108 
1109 
1110 #define ELLOUT(n) ov+(n-1)*3
1111 
1112 void
1113 rt_ell_norms(fastf_t *ov, fastf_t *A, fastf_t *B, fastf_t *h_vec, fastf_t t, int sides)
1114 {
1115  fastf_t ang, theta, x, y, sqrt_1mt;
1116  int n;
1117  vect_t partial_t, partial_ang;
1118 
1119  sqrt_1mt = sqrt(1.0 - t);
1120  if (sqrt_1mt <= SMALL_FASTF)
1121  bu_bomb("rt_epa_tess: rt_ell_norms: sqrt(1.0 -t) is zero\n");
1122  theta = M_2PI / sides;
1123  ang = 0.;
1124 
1125  for (n = 1; n <= sides; n++, ang += theta) {
1126  x = cos(ang);
1127  y = sin(ang);
1128  VJOIN2(partial_t, h_vec, -x/(2.0*sqrt_1mt), A, -y/(2.0*sqrt_1mt), B);
1129  VBLEND2(partial_ang, x*sqrt_1mt, B, -y*sqrt_1mt, A);
1130  VCROSS(ELLOUT(n), partial_t, partial_ang);
1131  VUNITIZE(ELLOUT(n));
1132  }
1133  VMOVE(ELLOUT(n), ELLOUT(1));
1134 }
1135 
1136 
1137 /**
1138  * Generate an ellipsoid with the specified number of sides approximating it.
1139  */
1140 void
1141 rt_ell(fastf_t *ov, const fastf_t *V, const fastf_t *A, const fastf_t *B, int sides)
1142 {
1143  fastf_t ang, theta, x, y;
1144  int n;
1145 
1146  theta = M_2PI / sides;
1147  ang = 0.;
1148  /* make ellipse regardless of whether it meets req's */
1149  for (n = 1; n <= sides; n++, ang += theta) {
1150  x = cos(ang);
1151  y = sin(ang);
1152  VJOIN2(ELLOUT(n), V, x, A, y, B);
1153  }
1154  VMOVE(ELLOUT(n), ELLOUT(1));
1155 }
1156 
1157 
1158 /**
1159  * Returns -
1160  * -1 failure
1161  * 0 OK. *r points to nmgregion that holds this tessellation.
1162  */
1163 int
1164 rt_epa_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
1165 {
1166  fastf_t dtol, mag_h, ntol, r1, r2;
1167  fastf_t **ellipses, **normals, theta_new, theta_prev;
1168  int *pts_dbl, face, i, j, nseg;
1169  int *segs_per_ell;
1170  int jj, na, nb, nell, recalc_b;
1171  mat_t R;
1172  mat_t invR;
1173  struct rt_epa_internal *xip;
1174  point_t p1;
1175  struct rt_pt_node *pos_a, *pos_b, *pts_a, *pts_b;
1176  struct shell *s;
1177  struct faceuse **outfaceuses = NULL;
1178  struct vertex *vertp[3];
1179  struct vertex ***vells = (struct vertex ***)NULL;
1180  vect_t A, Au, B, Bu, Hu, V;
1181  vect_t apex_norm, rev_norm;
1182  vect_t A_orig, B_orig;
1183  struct vertex *apex_v;
1184  struct vertexuse *vu;
1185  struct faceuse *fu;
1186 
1187  RT_CK_DB_INTERNAL(ip);
1188 
1189  xip = (struct rt_epa_internal *)ip->idb_ptr;
1190  if (!epa_is_valid(xip)) {
1191  return -2;
1192  }
1193 
1194  /* compute |H| */
1195  mag_h = MAGNITUDE(xip->epa_H);
1196  r1 = xip->epa_r1;
1197  r2 = xip->epa_r2;
1198 
1199 
1200  /* make unit vectors in A, H, and BxH directions */
1201  VMOVE(Hu, xip->epa_H);
1202  VUNITIZE(Hu);
1203  VMOVE(Au, xip->epa_Au);
1204  VCROSS(Bu, Au, Hu);
1205 
1206  VSCALE(A_orig, Au, xip->epa_r1);
1207  VSCALE(B_orig, Bu, xip->epa_r2);
1208 
1209  /* Compute R and Rinv matrices */
1210  MAT_IDN(R);
1211  VREVERSE(&R[0], Bu);
1212  VMOVE(&R[4], Au);
1213  VREVERSE(&R[8], Hu);
1214  bn_mat_trn(invR, R); /* inv of rot mat is trn */
1215 
1216  dtol = primitive_get_absolute_tolerance(ttol, 2.0 * r2);
1217 
1218  /* To ensure normal tolerance, remain below this angle */
1219  if (ttol->norm > 0.0)
1220  ntol = ttol->norm;
1221  else
1222  /* tolerate everything */
1223  ntol = M_PI;
1224 
1225  /*
1226  * build epa from 2 parabolas
1227  */
1228 
1229  /* approximate positive half of parabola along semi-minor axis */
1230  BU_ALLOC(pts_b, struct rt_pt_node);
1231  BU_ALLOC(pts_b->next, struct rt_pt_node);
1232 
1233  pts_b->next->next = NULL;
1234  VSET(pts_b->p, 0, 0, -mag_h);
1235  VSET(pts_b->next->p, 0, r2, 0);
1236  /* 2 endpoints in 1st approximation */
1237  nb = 2;
1238  /* recursively break segment 'til within error tolerances */
1239  nb += rt_mk_parabola(pts_b, r2, mag_h, dtol, ntol);
1240  nell = nb - 1; /* # of ellipses needed */
1241 
1242  /* construct positive half of parabola along semi-major axis of
1243  * epa using same z coords as parab along semi-minor axis
1244  */
1245  BU_ALLOC(pts_a, struct rt_pt_node);
1246  VMOVE(pts_a->p, pts_b->p); /* 1st pt is the apex */
1247  pts_a->next = NULL;
1248  pos_b = pts_b->next;
1249  pos_a = pts_a;
1250  while (pos_b) {
1251  /* copy node from b_parabola to a_parabola */
1252  BU_ALLOC(pos_a->next, struct rt_pt_node);
1253  pos_a = pos_a->next;
1254  pos_a->p[Z] = pos_b->p[Z];
1255  /* at given z, find y on parabola */
1256  pos_a->p[Y] = r1*sqrt(pos_a->p[Z] / mag_h + 1);
1257  pos_a->p[X] = 0;
1258  pos_b = pos_b->next;
1259  }
1260  pos_a->next = NULL;
1261 
1262  /* see if any segments need further breaking up */
1263  recalc_b = 0;
1264  na = 0;
1265  pos_a = pts_a;
1266  while (pos_a->next) {
1267  na = rt_mk_parabola(pos_a, r1, mag_h, dtol, ntol);
1268  if (na != 0) {
1269  recalc_b = 1;
1270  nell += na;
1271  }
1272  pos_a = pos_a->next;
1273  }
1274  /* if any were broken, recalculate parabola 'a' */
1275  if (recalc_b) {
1276  /* free mem for old approximation of parabola */
1277  pos_b = pts_b;
1278  while (pos_b) {
1279  struct rt_pt_node *tmp;
1280 
1281  tmp = pos_b->next;
1282  bu_free((char *)pos_b, "rt_pt_node");
1283  pos_b = tmp;
1284  }
1285  /* construct parabola along semi-major axis of epa
1286  * using same z coords as parab along semi-minor axis
1287  */
1288  BU_ALLOC(pts_b, struct rt_pt_node);
1289  pts_b->p[Z] = pts_a->p[Z];
1290  pts_b->next = NULL;
1291  pos_a = pts_a->next;
1292  pos_b = pts_b;
1293  while (pos_a) {
1294  /* copy node from a_parabola to b_parabola */
1295  BU_ALLOC(pos_b->next, struct rt_pt_node);
1296  pos_b = pos_b->next;
1297  pos_b->p[Z] = pos_a->p[Z];
1298  /* at given z, find y on parabola */
1299  pos_b->p[Y] = r2*sqrt(pos_b->p[Z] / mag_h + 1);
1300  pos_b->p[X] = 0;
1301  pos_a = pos_a->next;
1302  }
1303  pos_b->next = NULL;
1304  }
1305 
1306  /* make array of ptrs to epa ellipses */
1307  ellipses = (fastf_t **)bu_malloc(nell * sizeof(fastf_t *), "fastf_t ell[]");
1308  /* keep track of whether pts in each ellipse are doubled or not */
1309  pts_dbl = (int *)bu_malloc(nell * sizeof(int), "dbl ints");
1310  /* I don't understand this pts_dbl, so here is an array containing
1311  * the length of each ellipses array
1312  */
1313  segs_per_ell = (int *)bu_calloc(nell, sizeof(int), "rt_epa_tess: segs_per_ell");
1314 
1315  /* and an array of normals */
1316  normals = (fastf_t **)bu_malloc(nell * sizeof(fastf_t *), "fastf_t normals[]");
1317 
1318  /* make ellipses at each z level */
1319  i = 0;
1320  nseg = 0;
1321  theta_prev = M_2PI;
1322  pos_a = pts_a->next; /* skip over apex of epa */
1323  pos_b = pts_b->next;
1324  while (pos_a) {
1325  fastf_t t;
1326 
1327  t = (-pos_a->p[Z] / mag_h);
1328  VSCALE(A, Au, pos_a->p[Y]); /* semimajor axis */
1329  VSCALE(B, Bu, pos_b->p[Y]); /* semiminor axis */
1330  VJOIN1(V, xip->epa_V, -pos_a->p[Z], Hu);
1331 
1332  VSET(p1, 0., pos_b->p[Y], 0.);
1333  theta_new = ell_angle(p1, pos_a->p[Y], pos_b->p[Y], dtol, ntol);
1334  if (nseg == 0) {
1335  nseg = (int)(M_2PI / theta_new) + 1;
1336  pts_dbl[i] = 0;
1337  /* maximum number of faces needed for epa */
1338  face = nseg*(1 + 3*((1 << (nell-1)) - 1));
1339  /* array for each triangular face */
1340  outfaceuses = (struct faceuse **)
1341  bu_malloc((face+1) * sizeof(struct faceuse *), "faceuse []");
1342  } else if (theta_new < theta_prev) {
1343  nseg *= 2;
1344  pts_dbl[i] = 1;
1345  } else {
1346  pts_dbl[i] = 0;
1347  }
1348  theta_prev = theta_new;
1349 
1350  ellipses[i] = (fastf_t *)bu_malloc(3*(nseg+1)*sizeof(fastf_t),
1351  "pts ell");
1352  segs_per_ell[i] = nseg;
1353  normals[i] = (fastf_t *)bu_malloc(3*(nseg+1)*sizeof(fastf_t), "rt_epa_tess_ normals");
1354  rt_ell(ellipses[i], V, A, B, nseg);
1355  rt_ell_norms(normals[i], A_orig, B_orig, xip->epa_H, t, nseg);
1356 
1357  i++;
1358  pos_a = pos_a->next;
1359  pos_b = pos_b->next;
1360  }
1361 
1362  /*
1363  * put epa geometry into nmg data structures
1364  */
1365 
1366  *r = nmg_mrsv(m); /* Make region, empty shell, vertex */
1367  s = BU_LIST_FIRST(shell, &(*r)->s_hd);
1368 
1369  /* vertices of ellipses of epa */
1370  vells = (struct vertex ***)
1371  bu_malloc(nell*sizeof(struct vertex **), "vertex [][]");
1372  j = nseg;
1373  for (i = nell-1; i >= 0; i--) {
1374  vells[i] = (struct vertex **)
1375  bu_malloc(j*sizeof(struct vertex *), "vertex []");
1376  if (i && pts_dbl[i])
1377  j /= 2;
1378  }
1379 
1380  /* top face of epa */
1381  for (i = 0; i < nseg; i++)
1382  vells[nell-1][i] = (struct vertex *)0;
1383  face = 0;
1384  if ((outfaceuses[face++] = nmg_cface(s, vells[nell-1], nseg)) == 0) {
1385  bu_log("rt_epa_tess() failure, top face\n");
1386  goto fail;
1387  }
1388  for (i = 0; i < nseg; i++) {
1389  NMG_CK_VERTEX(vells[nell-1][i]);
1390  nmg_vertex_gv(vells[nell-1][i], &ellipses[nell-1][3*i]);
1391  }
1392 
1393  /* Mark the edges of this face as real, this is the only real edge */
1394  (void)nmg_mark_edges_real(&outfaceuses[0]->l.magic);
1395 
1396  /* connect ellipses with triangles */
1397  for (i = nell-2; i >= 0; i--) {
1398  /* skip top ellipse */
1399  int bottom, top;
1400 
1401  top = i + 1;
1402  bottom = i;
1403  if (pts_dbl[top])
1404  nseg /= 2; /* # segs in 'bottom' ellipse */
1405  vertp[0] = (struct vertex *)0;
1406 
1407  /* make triangular faces */
1408  for (j = 0; j < nseg; j++) {
1409  jj = j + j; /* top ellipse index */
1410 
1411  if (pts_dbl[top]) {
1412  /* first triangle */
1413  vertp[1] = vells[top][jj+1];
1414  vertp[2] = vells[top][jj];
1415  if ((outfaceuses[face++] = nmg_cface(s, vertp, 3)) == 0) {
1416  bu_log("rt_epa_tess() failure\n");
1417  goto fail;
1418  }
1419  if (j == 0)
1420  vells[bottom][j] = vertp[0];
1421 
1422  /* second triangle */
1423  vertp[2] = vertp[1];
1424  if (j == nseg-1)
1425  vertp[1] = vells[bottom][0];
1426  else
1427  vertp[1] = (struct vertex *)0;
1428  if ((outfaceuses[face++] = nmg_cface(s, vertp, 3)) == 0) {
1429  bu_log("rt_epa_tess() failure\n");
1430  goto fail;
1431  }
1432  if (j != nseg-1)
1433  vells[bottom][j+1] = vertp[1];
1434 
1435  /* third triangle */
1436  vertp[0] = vertp[1];
1437  if (j == nseg-1)
1438  vertp[1] = vells[top][0];
1439  else
1440  vertp[1] = vells[top][jj+2];
1441  if ((outfaceuses[face++] = nmg_cface(s, vertp, 3)) == 0) {
1442  bu_log("rt_epa_tess() failure\n");
1443  goto fail;
1444  }
1445  } else {
1446  /* first triangle */
1447  if (j == nseg-1)
1448  vertp[1] = vells[top][0];
1449  else
1450  vertp[1] = vells[top][j+1];
1451  vertp[2] = vells[top][j];
1452  if ((outfaceuses[face++] = nmg_cface(s, vertp, 3)) == 0) {
1453  bu_log("rt_epa_tess() failure\n");
1454  goto fail;
1455  }
1456  if (j == 0)
1457  vells[bottom][j] = vertp[0];
1458 
1459  /* second triangle */
1460  vertp[2] = vertp[0];
1461  if (j == nseg-1)
1462  vertp[0] = vells[bottom][0];
1463  else
1464  vertp[0] = (struct vertex *)0;
1465  if ((outfaceuses[face++] = nmg_cface(s, vertp, 3)) == 0) {
1466  bu_log("rt_epa_tess() failure\n");
1467  goto fail;
1468  }
1469  if (j != nseg-1)
1470  vells[bottom][j+1] = vertp[0];
1471  }
1472  }
1473 
1474  /* associate geometry with each vertex */
1475  for (j = 0; j < nseg; j++) {
1476  NMG_CK_VERTEX(vells[bottom][j]);
1477  nmg_vertex_gv(vells[bottom][j],
1478  &ellipses[bottom][3*j]);
1479  }
1480  }
1481 
1482  /* connect bottom of ellipse to apex of epa */
1483  VADD2(V, xip->epa_V, xip->epa_H);
1484  vertp[0] = (struct vertex *)0;
1485  vertp[1] = vells[0][1];
1486  vertp[2] = vells[0][0];
1487  if ((outfaceuses[face++] = nmg_cface(s, vertp, 3)) == 0) {
1488  bu_log("rt_epa_tess() failure\n");
1489  goto fail;
1490  }
1491  /* associate geometry with topology */
1492  NMG_CK_VERTEX(vertp[0]);
1493  nmg_vertex_gv(vertp[0], (fastf_t *)V);
1494  apex_v = vertp[0];
1495  /* create rest of faces around apex */
1496  for (i = 1; i < nseg; i++) {
1497  vertp[2] = vertp[1];
1498  if (i == nseg-1)
1499  vertp[1] = vells[0][0];
1500  else
1501  vertp[1] = vells[0][i+1];
1502  if ((outfaceuses[face++] = nmg_cface(s, vertp, 3)) == 0) {
1503  bu_log("rt_epa_tess() failure\n");
1504  goto fail;
1505  }
1506  }
1507 
1508  /* Associate the face geometry */
1509  for (i = 0; i < face; i++) {
1510  if (nmg_fu_planeeqn(outfaceuses[i], tol) < 0)
1511  goto fail;
1512  }
1513 
1514  /* Associate vertexuse normals */
1515  for (i = 0; i < nell; i++) {
1516  for (j = 0; j < segs_per_ell[i]; j++) {
1517  VREVERSE(rev_norm, &normals[i][j*3]);
1518  for (BU_LIST_FOR(vu, vertexuse, &vells[i][j]->vu_hd)) {
1519 
1520  fu = nmg_find_fu_of_vu(vu);
1521  NMG_CK_FACEUSE(fu);
1522 
1523  if (fu == outfaceuses[0] || fu->fumate_p == outfaceuses[0])
1524  continue; /* don't assign normals to top faceuse (flat) */
1525 
1526  if (fu->orientation == OT_SAME)
1527  nmg_vertexuse_nv(vu, &normals[i][j*3]);
1528  else if (fu->orientation == OT_OPPOSITE)
1529  nmg_vertexuse_nv(vu, rev_norm);
1530  }
1531  }
1532  }
1533  /* and don't forget the apex */
1534  VMOVE(apex_norm, xip->epa_H);
1535  VUNITIZE(apex_norm);
1536  VREVERSE(rev_norm, apex_norm);
1537  for (BU_LIST_FOR(vu, vertexuse, &apex_v->vu_hd)) {
1538  NMG_CK_VERTEXUSE(vu);
1539  fu = nmg_find_fu_of_vu(vu);
1540  NMG_CK_FACEUSE(fu);
1541  if (fu->orientation == OT_SAME)
1542  nmg_vertexuse_nv(vu, apex_norm);
1543  else if (fu->orientation == OT_OPPOSITE)
1544  nmg_vertexuse_nv(vu, rev_norm);
1545  }
1546 
1547  /* Glue the edges of different outward pointing face uses together */
1548  nmg_gluefaces(outfaceuses, face, tol);
1549 
1550  /* Compute "geometry" for region and shell */
1551  nmg_region_a(*r, tol);
1552 
1553  /* XXX just for testing, to make up for loads of triangles ... */
1554  nmg_shell_coplanar_face_merge(s, tol, 1);
1555 
1556  /* free mem */
1557  bu_free((char *)outfaceuses, "faceuse []");
1558  for (i = 0; i < nell; i++) {
1559  bu_free((char *)ellipses[i], "pts ell");
1560  bu_free((char *)vells[i], "vertex []");
1561  }
1562  bu_free((char *)ellipses, "fastf_t ell[]");
1563  bu_free((char *)pts_dbl, "dbl ints");
1564  bu_free((char *)vells, "vertex [][]");
1565 
1566  return 0;
1567 
1568  fail:
1569  /* free mem */
1570  bu_free((char *)outfaceuses, "faceuse []");
1571  for (i = 0; i < nell; i++) {
1572  bu_free((char *)ellipses[i], "pts ell");
1573  bu_free((char *)vells[i], "vertex []");
1574  }
1575  bu_free((char *)ellipses, "fastf_t ell[]");
1576  bu_free((char *)pts_dbl, "dbl ints");
1577  bu_free((char *)vells, "vertex [][]");
1578 
1579  return -1;
1580 }
1581 
1582 
1583 /**
1584  * Import an EPA from the database format to the internal format.
1585  * Apply modeling transformations as well.
1586  */
1587 int
1588 rt_epa_import4(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
1589 {
1590  struct rt_epa_internal *xip;
1591  union record *rp;
1592  vect_t v1, v2, v3;
1593 
1594  if (dbip) RT_CK_DBI(dbip);
1595 
1596  BU_CK_EXTERNAL(ep);
1597  rp = (union record *)ep->ext_buf;
1598  /* Check record type */
1599  if (rp->u_id != ID_SOLID) {
1600  bu_log("rt_epa_import4: defective record\n");
1601  return -1;
1602  }
1603 
1604  RT_CK_DB_INTERNAL(ip);
1605  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1606  ip->idb_type = ID_EPA;
1607  ip->idb_meth = &OBJ[ID_EPA];
1608  BU_ALLOC(ip->idb_ptr, struct rt_epa_internal);
1609 
1610  xip = (struct rt_epa_internal *)ip->idb_ptr;
1611  xip->epa_magic = RT_EPA_INTERNAL_MAGIC;
1612 
1613  /* Warning: type conversion */
1614  if (mat == NULL) mat = bn_mat_identity;
1615 
1616  if (dbip->dbi_version < 0) {
1617  flip_fastf_float(v1, &rp->s.s_values[0*3], 1, 1);
1618  flip_fastf_float(v2, &rp->s.s_values[1*3], 1, 1);
1619  flip_fastf_float(v3, &rp->s.s_values[2*3], 1, 1);
1620  } else {
1621  VMOVE(v1, &rp->s.s_values[0*3]);
1622  VMOVE(v2, &rp->s.s_values[1*3]);
1623  VMOVE(v3, &rp->s.s_values[2*3]);
1624  }
1625 
1626  MAT4X3PNT(xip->epa_V, mat, v1);
1627  MAT4X3VEC(xip->epa_H, mat, v2);
1628  MAT4X3VEC(xip->epa_Au, mat, v3);
1629 
1630  VUNITIZE(xip->epa_Au);
1631 
1632  if (dbip->dbi_version < 0) {
1633  v1[X] = flip_dbfloat(rp->s.s_values[3*3+0]);
1634  v1[Y] = flip_dbfloat(rp->s.s_values[3*3+1]);
1635  } else {
1636  v1[X] = rp->s.s_values[3*3+0];
1637  v1[Y] = rp->s.s_values[3*3+1];
1638  }
1639 
1640  xip->epa_r1 = v1[X] / mat[15];
1641  xip->epa_r2 = v1[Y] / mat[15];
1642 
1643  if (xip->epa_r1 <= SMALL_FASTF || xip->epa_r2 <= SMALL_FASTF) {
1644  bu_log("rt_epa_import4: r1 or r2 are zero\n");
1645  bu_free((char *)ip->idb_ptr, "rt_epa_import4: ip->idb_ptr");
1646  return -1;
1647  }
1648 
1649  return 0; /* OK */
1650 }
1651 
1652 
1653 /**
1654  * The name is added by the caller, in the usual place.
1655  */
1656 int
1657 rt_epa_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1658 {
1659  struct rt_epa_internal *xip;
1660  union record *epa;
1661  fastf_t mag_h;
1662 
1663  if (dbip) RT_CK_DBI(dbip);
1664 
1665  RT_CK_DB_INTERNAL(ip);
1666  if (ip->idb_type != ID_EPA) return -1;
1667  xip = (struct rt_epa_internal *)ip->idb_ptr;
1668  RT_EPA_CK_MAGIC(xip);
1669 
1670  BU_CK_EXTERNAL(ep);
1671  ep->ext_nbytes = sizeof(union record);
1672  ep->ext_buf = (uint8_t *)bu_calloc(1, ep->ext_nbytes, "epa external");
1673  epa = (union record *)ep->ext_buf;
1674 
1675  epa->s.s_id = ID_SOLID;
1676  epa->s.s_type = EPA;
1677 
1678  if (!NEAR_EQUAL(MAGNITUDE(xip->epa_Au), 1.0, RT_LEN_TOL)) {
1679  bu_log("rt_epa_export4: Au not a unit vector!\n");
1680  return -1;
1681  }
1682 
1683  mag_h = MAGNITUDE(xip->epa_H);
1684 
1685  if (mag_h < RT_LEN_TOL
1686  || xip->epa_r1 < RT_LEN_TOL
1687  || xip->epa_r2 < RT_LEN_TOL) {
1688  bu_log("rt_epa_export4: not all dimensions positive!\n");
1689  return -1;
1690  }
1691 
1692  if (!NEAR_ZERO(VDOT(xip->epa_Au, xip->epa_H)/mag_h, RT_DOT_TOL)) {
1693  bu_log("rt_epa_export4: Au and H are not perpendicular!\n");
1694  return -1;
1695  }
1696 
1697  if (xip->epa_r2 > xip->epa_r1) {
1698  bu_log("rt_epa_export4: semi-minor axis cannot be longer than semi-major axis!\n");
1699  return -1;
1700  }
1701 
1702  /* Warning: type conversion */
1703  VSCALE(&epa->s.s_values[0*3], xip->epa_V, local2mm);
1704  VSCALE(&epa->s.s_values[1*3], xip->epa_H, local2mm);
1705  VMOVE(&epa->s.s_values[2*3], xip->epa_Au); /* don't scale a unit vector */
1706  epa->s.s_values[3*3] = xip->epa_r1 * local2mm;
1707  epa->s.s_values[3*3+1] = xip->epa_r2 * local2mm;
1708 
1709  return 0;
1710 }
1711 
1712 
1713 /**
1714  * Import an EPA from the database format to the internal format.
1715  * Apply modeling transformations as well.
1716  */
1717 int
1718 rt_epa_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
1719 {
1720  struct rt_epa_internal *xip;
1721 
1722  /* must be double for import and export */
1723  double vec[11];
1724 
1725  if (dbip) RT_CK_DBI(dbip);
1726  BU_CK_EXTERNAL(ep);
1727 
1729 
1730  RT_CK_DB_INTERNAL(ip);
1731  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1732  ip->idb_type = ID_EPA;
1733  ip->idb_meth = &OBJ[ID_EPA];
1734  BU_ALLOC(ip->idb_ptr, struct rt_epa_internal);
1735 
1736  xip = (struct rt_epa_internal *)ip->idb_ptr;
1737  xip->epa_magic = RT_EPA_INTERNAL_MAGIC;
1738 
1739  /* Convert from database (network) to internal (host) format */
1740  bu_cv_ntohd((unsigned char *)vec, ep->ext_buf, 11);
1741 
1742  /* Apply modeling transformations */
1743  if (mat == NULL) mat = bn_mat_identity;
1744  MAT4X3PNT(xip->epa_V, mat, &vec[0*3]);
1745  MAT4X3VEC(xip->epa_H, mat, &vec[1*3]);
1746  MAT4X3VEC(xip->epa_Au, mat, &vec[2*3]);
1747  VUNITIZE(xip->epa_Au);
1748  xip->epa_r1 = vec[3*3] / mat[15];
1749  xip->epa_r2 = vec[3*3+1] / mat[15];
1750 
1751  if (xip->epa_r1 <= SMALL_FASTF || xip->epa_r2 <= SMALL_FASTF) {
1752  bu_log("rt_epa_import4: r1 or r2 are zero\n");
1753  bu_free((char *)ip->idb_ptr, "rt_epa_import4: ip->idb_ptr");
1754  return -1;
1755  }
1756 
1757  return 0; /* OK */
1758 }
1759 
1760 
1761 /**
1762  * The name is added by the caller, in the usual place.
1763  */
1764 int
1765 rt_epa_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1766 {
1767  struct rt_epa_internal *xip;
1768  fastf_t mag_h;
1769 
1770  /* must be double for import and export */
1771  double vec[11];
1772 
1773  if (dbip) RT_CK_DBI(dbip);
1774 
1775  RT_CK_DB_INTERNAL(ip);
1776  if (ip->idb_type != ID_EPA) return -1;
1777  xip = (struct rt_epa_internal *)ip->idb_ptr;
1778  RT_EPA_CK_MAGIC(xip);
1779 
1780  BU_CK_EXTERNAL(ep);
1781  ep->ext_nbytes = SIZEOF_NETWORK_DOUBLE * 11;
1782  ep->ext_buf = (uint8_t *)bu_malloc(ep->ext_nbytes, "epa external");
1783 
1784  if (!NEAR_EQUAL(MAGNITUDE(xip->epa_Au), 1.0, RT_LEN_TOL)) {
1785  bu_log("rt_epa_export4: Au not a unit vector!\n");
1786  return -1;
1787  }
1788 
1789  mag_h = MAGNITUDE(xip->epa_H);
1790 
1791  if (mag_h < RT_LEN_TOL
1792  || xip->epa_r1 < RT_LEN_TOL
1793  || xip->epa_r2 < RT_LEN_TOL) {
1794  bu_log("rt_epa_export4: not all dimensions positive!\n");
1795  return -1;
1796  }
1797 
1798  if (!NEAR_ZERO(VDOT(xip->epa_Au, xip->epa_H)/mag_h, RT_DOT_TOL)) {
1799  bu_log("rt_epa_export4: Au and H are not perpendicular!\n");
1800  return -1;
1801  }
1802 
1803  if (xip->epa_r2 > xip->epa_r1) {
1804  bu_log("rt_epa_export4: semi-minor axis cannot be longer than semi-major axis!\n");
1805  return -1;
1806  }
1807 
1808  /* scale 'em into local buffer */
1809  VSCALE(&vec[0*3], xip->epa_V, local2mm);
1810  VSCALE(&vec[1*3], xip->epa_H, local2mm);
1811  VMOVE(&vec[2*3], xip->epa_Au); /* don't scale a unit vector */
1812  vec[3*3] = xip->epa_r1 * local2mm;
1813  vec[3*3+1] = xip->epa_r2 * local2mm;
1814 
1815  /* Convert from internal (host) to database (network) format */
1816  bu_cv_htond(ep->ext_buf, (unsigned char *)vec, 11);
1817 
1818  return 0;
1819 }
1820 
1821 
1822 /**
1823  * Make human-readable formatted presentation of this solid. First
1824  * line describes type of solid. Additional lines are indented one
1825  * tab, and give parameter values.
1826  */
1827 int
1828 rt_epa_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
1829 {
1830  struct rt_epa_internal *xip = (struct rt_epa_internal *)ip->idb_ptr;
1831  char buf[256];
1832 
1833  RT_EPA_CK_MAGIC(xip);
1834  bu_vls_strcat(str, "Elliptical Paraboloid (EPA)\n");
1835 
1836  sprintf(buf, "\tV (%g, %g, %g)\n",
1837  INTCLAMP(xip->epa_V[X] * mm2local),
1838  INTCLAMP(xip->epa_V[Y] * mm2local),
1839  INTCLAMP(xip->epa_V[Z] * mm2local));
1840  bu_vls_strcat(str, buf);
1841 
1842  sprintf(buf, "\tH (%g, %g, %g) mag=%g\n",
1843  INTCLAMP(xip->epa_H[X] * mm2local),
1844  INTCLAMP(xip->epa_H[Y] * mm2local),
1845  INTCLAMP(xip->epa_H[Z] * mm2local),
1846  INTCLAMP(MAGNITUDE(xip->epa_H) * mm2local));
1847  bu_vls_strcat(str, buf);
1848 
1849  if (!verbose)
1850  return 0;
1851 
1852  sprintf(buf, "\tA=%g\n", INTCLAMP(xip->epa_r1 * mm2local));
1853  bu_vls_strcat(str, buf);
1854 
1855  sprintf(buf, "\tB=%g\n", INTCLAMP(xip->epa_r2 * mm2local));
1856  bu_vls_strcat(str, buf);
1857 
1858  return 0;
1859 }
1860 
1861 
1862 /**
1863  * Free the storage associated with the rt_db_internal version of this
1864  * solid.
1865  */
1866 void
1868 {
1869  struct rt_epa_internal *xip;
1870 
1871  RT_CK_DB_INTERNAL(ip);
1872 
1873  xip = (struct rt_epa_internal *)ip->idb_ptr;
1874  RT_EPA_CK_MAGIC(xip);
1875  xip->epa_magic = 0; /* sanity */
1876 
1877  bu_free((char *)xip, "epa ifree");
1878  ip->idb_ptr = ((void *)0); /* sanity */
1879 }
1880 
1881 
1882 int
1883 rt_epa_params(struct pc_pc_set *ps, const struct rt_db_internal *ip)
1884 {
1885  if (!ps) return 0;
1886  if (ip) RT_CK_DB_INTERNAL(ip);
1887 
1888  return 0; /* OK */
1889 }
1890 
1891 
1892 void
1893 rt_epa_volume(fastf_t *vol, const struct rt_db_internal *ip)
1894 {
1895  fastf_t mag_h;
1896  struct rt_epa_internal *xip = (struct rt_epa_internal *)ip->idb_ptr;
1897  RT_EPA_CK_MAGIC(xip);
1898 
1899  mag_h = MAGNITUDE(xip->epa_H);
1900  *vol = M_PI_2 * xip->epa_r1 * xip->epa_r2 * mag_h;
1901 }
1902 
1903 
1904 void
1905 rt_epa_centroid(point_t *cent, const struct rt_db_internal *ip)
1906 {
1907  struct rt_epa_internal *xip = (struct rt_epa_internal *)ip->idb_ptr;
1908  RT_EPA_CK_MAGIC(xip);
1909  VJOIN1(*cent, xip->epa_V, 1.0/3.0, xip->epa_H);
1910 }
1911 
1912 
1913 void
1914 rt_epa_surf_area(fastf_t *area, const struct rt_db_internal *ip)
1915 {
1916  fastf_t magsq_h, m;
1917  struct rt_epa_internal *xip = (struct rt_epa_internal *)ip->idb_ptr;
1918  RT_EPA_CK_MAGIC(xip);
1919 
1920  magsq_h = MAGSQ(xip->epa_H);
1921  m = sqrt(1.0 + (4.0 * magsq_h) / (xip->epa_r1 * xip->epa_r2));
1922  *area = 2.0/3.0 * M_PI * xip->epa_r1 * xip->epa_r2 * (m + (1.0 / (m + 1.0)));
1923 }
1924 
1925 static int
1926 epa_is_valid(struct rt_epa_internal *epa)
1927 {
1928  fastf_t mag_h, cos_angle_ah;
1929  vect_t epa_H, epa_Au;
1930 
1931  RT_EPA_CK_MAGIC(epa);
1932 
1933  if (!(epa->epa_r1 > 0.0 && epa->epa_r2 > 0.0)) {
1934  return 0;
1935  }
1936 
1937  VMOVE(epa_H, epa->epa_H);
1938  VMOVE(epa_Au, epa->epa_Au);
1939 
1940  /* Check that Au is a unit vector. If it is, then it should be true that
1941  * |Au| == |Au|^2 == 1.0.
1942  */
1943  if (!NEAR_EQUAL(MAGSQ(epa_Au), 1.0, RT_LEN_TOL)) {
1944  return 0;
1945  }
1946 
1947  /* check that |H| > 0.0 */
1948  mag_h = MAGNITUDE(epa_H);
1949  if (NEAR_ZERO(mag_h, RT_LEN_TOL)) {
1950  return 0;
1951  }
1952 
1953  /* check that A and H are orthogonal */
1954  cos_angle_ah = VDOT(epa_Au, epa_H) / mag_h;
1955  if (!NEAR_ZERO(cos_angle_ah, RT_DOT_TOL)) {
1956  return 0;
1957  }
1958 
1959  return 1;
1960 }
1961 
1962 /** @} */
1963 /*
1964  * Local Variables:
1965  * mode: C
1966  * tab-width: 8
1967  * indent-tabs-mode: t
1968  * c-file-style: "stroustrup"
1969  * End:
1970  * ex: shiftwidth=4 tabstop=8
1971  */
#define BU_LIST_FOR(p, structure, hp)
Definition: list.h:365
Definition: raytrace.h:800
#define RT_LEN_TOL
Definition: raytrace.h:169
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define BU_LIST_INSERT(old, new)
Definition: list.h:183
#define SIZEOF_NETWORK_DOUBLE
Definition: cv.h:48
struct faceuse * nmg_cface(struct shell *s, struct vertex **verts, int n)
Definition: nmg_mod.c:1130
#define ELLOUT(n)
Definition: epa.c:1110
struct hit seg_in
IN information.
Definition: raytrace.h:370
Definition: list.h:118
int nmg_fu_planeeqn(struct faceuse *fu, const struct bn_tol *tol)
Definition: nmg_mod.c:1311
#define RT_DOT_TOL
Definition: raytrace.h:170
#define RT_CK_APPLICATION(_p)
Definition: raytrace.h:1675
int rt_epa_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
Definition: epa.c:1828
vect_t crv_pdir
Principle direction.
Definition: raytrace.h:307
const mat_t bn_mat_identity
Matrix and vector functionality.
Definition: mat.c:46
fastf_t uv_u
Range 0..1.
Definition: raytrace.h:341
struct soltab * seg_stp
pointer back to soltab
Definition: raytrace.h:372
void rt_ell_norms(fastf_t *ov, fastf_t *A, fastf_t *B, fastf_t *h_vec, fastf_t t, int sides)
Definition: epa.c:1113
if lu s
Definition: nmg_mod.c:3860
void bu_vls_strcat(struct bu_vls *vp, const char *s)
Definition: vls.c:368
#define VSET(a, b, c, d)
Definition: color.c:53
#define VSETALL(a, s)
Definition: color.c:54
Definition: raytrace.h:215
const struct bu_structparse rt_epa_parse[]
Definition: epa.c:178
#define M_PI
Definition: fft.h:35
Definition: pc.h:108
Definition: raytrace.h:368
#define BU_ASSERT_LONG(_lhs, _relation, _rhs)
Definition: defines.h:240
void rt_epa_volume(fastf_t *vol, const struct rt_db_internal *ip)
Definition: epa.c:1893
int rt_epa_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
Definition: epa.c:361
int rt_epa_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: epa.c:253
Definition: raytrace.h:248
void nmg_vertex_gv(struct vertex *v, const fastf_t *pt)
Definition: nmg_mk.c:1668
#define SMALL_FASTF
Definition: defines.h:342
fastf_t st_aradius
Radius of APPROXIMATING sphere.
Definition: raytrace.h:433
fastf_t ell_angle(fastf_t *p1, fastf_t a, fastf_t b, fastf_t dtol, fastf_t ntol)
Definition: ell.c:1788
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)
void flip_fastf_float(fastf_t *ff, const dbfloat_t *fp, int n, int flip)
Definition: db_flip.c:74
struct resource * a_resource
dynamic memory resources
Definition: raytrace.h:1591
#define EPA_NORM_BODY
Definition: epa.c:348
void bu_cv_htond(unsigned char *out, const unsigned char *in, size_t count)
#define RT_EPA_INTERNAL_MAGIC
Definition: magic.h:92
void rt_epa_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
Definition: epa.c:513
struct bu_list l
Definition: raytrace.h:369
fastf_t curve_spacing
Definition: raytrace.h:1940
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
void rt_epa_surf_area(fastf_t *area, const struct rt_db_internal *ip)
Definition: epa.c:1914
int rt_epa_import4(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
Definition: epa.c:1588
if(share_geom)
Definition: nmg_mod.c:3829
int idb_major_type
Definition: raytrace.h:192
Definition: color.c:49
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
mat_t epa_SoR
Definition: epa.c:173
mat_t epa_invRoS
Definition: epa.c:174
void nmg_shell_coplanar_face_merge(struct shell *s, const struct bn_tol *tol, const int simplify)
Definition: nmg_mod.c:93
#define RT_CK_DB_INTERNAL(_p)
Definition: raytrace.h:207
void rt_epa_centroid(point_t *cent, const struct rt_db_internal *ip)
Definition: epa.c:1905
#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
fastf_t crv_c2
curvature in other direction
Definition: raytrace.h:309
#define BN_VLIST_LINE_MOVE
Definition: vlist.h:82
const struct rt_functab * idb_meth
for ft_ifree(), etc.
Definition: raytrace.h:194
int rt_epa_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: epa.c:1657
fastf_t uv_dv
delta in v
Definition: raytrace.h:344
fastf_t primitive_get_absolute_tolerance(const struct rt_tess_tol *ttol, fastf_t rel_to_abs)
#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
void rt_epa_free(struct soltab *stp)
Definition: epa.c:614
int rt_mk_parabola(struct rt_pt_node *pts, fastf_t r, fastf_t b, fastf_t dtol, fastf_t ntol)
Definition: rpc.c:1009
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
struct hit seg_out
OUT information.
Definition: raytrace.h:371
#define BN_VLIST_LINE_DRAW
Definition: vlist.h:83
fastf_t epa_inv_r2sq
Definition: epa.c:172
fastf_t flip_dbfloat(dbfloat_t d)
Definition: db_flip.c:58
fastf_t epa_inv_r1sq
Definition: epa.c:171
void bn_mat_trn(mat_t om, const mat_t im)
void rt_epa_ifree(struct rt_db_internal *ip)
Definition: epa.c:1867
#define UNUSED(parameter)
Definition: common.h:239
int nmg_mark_edges_real(const uint32_t *magic_p)
Definition: nmg_misc.c:846
#define BU_PUT(_ptr, _type)
Definition: malloc.h:215
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
void rt_ell(fastf_t *ov, const fastf_t *V, const fastf_t *A, const fastf_t *B, int sides)
Definition: epa.c:1141
#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
vect_t r_dir
Direction of ray (UNIT Length)
Definition: raytrace.h:219
struct bn_tol rti_tol
Math tolerances for this model.
Definition: raytrace.h:1765
struct rt_pt_node * next
ptr to next pt
Definition: raytrace.h:1912
int rt_epa_params(struct pc_pc_set *ps, const struct rt_db_internal *ip)
Definition: epa.c:1883
#define bu_offsetof(_t, _m)
Definition: parse.h:64
#define ID_EPA
Elliptical Paraboloid.
Definition: raytrace.h:477
#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
vect_t epa_Hunit
Definition: epa.c:167
#define ZERO(val)
Definition: units.c:38
int rt_epa_adaptive_plot(struct rt_db_internal *ip, const struct rt_view_info *info)
Definition: epa.c:797
void * idb_ptr
Definition: raytrace.h:195
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
void bu_cv_ntohd(unsigned char *out, const unsigned char *in, size_t count)
const struct rt_functab OBJ[]
Definition: table.c:159
fastf_t epa_h
Definition: epa.c:170
void bn_vec_ortho(vect_t out, const vect_t in)
#define R
Definition: msr.c:55
int rt_epa_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
Definition: epa.c:1718
vect_t epa_Bunit
Definition: epa.c:169
void * st_specific
-> ID-specific (private) struct
Definition: raytrace.h:435
fastf_t uv_du
delta in u
Definition: raytrace.h:343
int rt_epa_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *tol)
Definition: epa.c:191
point_t p
a point
Definition: raytrace.h:1911
fastf_t crv_c1
curvature in principle dir
Definition: raytrace.h:308
#define A
Definition: msr.c:51
int rt_epa_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
Definition: epa.c:1164
Definition: color.c:51
int dbi_version
PRIVATE: use db_version()
Definition: raytrace.h:824
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
int hit_surfno
solid-specific surface indicator
Definition: raytrace.h:255
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
void rt_epa_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
Definition: epa.c:568
#define BU_CK_LIST_HEAD(_p)
Definition: list.h:142
#define BU_CK_EXTERNAL(_p)
Definition: parse.h:224
void nmg_gluefaces(struct faceuse **fulist, int n, const struct bn_tol *tol)
Definition: nmg_mod.c:1408
int rt_epa_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: epa.c:878
vect_t epa_Aunit
Definition: epa.c:168
const struct rt_functab * st_meth
pointer to per-solid methods
Definition: raytrace.h:428
#define RT_PCOEF_TOL
Definition: raytrace.h:171
int st_id
Solid ident.
Definition: raytrace.h:431
void rt_epa_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
Definition: epa.c:478
size_t ext_nbytes
Definition: parse.h:210
fastf_t hit_dist
dist from r_pt to hit_point
Definition: raytrace.h:250
HIDDEN void verbose(struct human_data_t *dude)
Definition: human.c:2008
Definition: vls.h:56
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
double fastf_t
Definition: defines.h:300
#define VPRINT(a, b)
Definition: raytrace.h:1881
int rt_epa_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: epa.c:1765
#define EPA_NORM_TOP
Definition: epa.c:349
point_t epa_V
Definition: epa.c:166
void rt_epa_print(const struct soltab *stp)
Definition: epa.c:333
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
int approximate_parabolic_curve(struct rt_pt_node *pts, fastf_t p, int num_new_points)
point_t st_center
Centroid of solid.
Definition: raytrace.h:432
#define RT_HIT_MAGIC
Definition: magic.h:161
void nmg_region_a(struct nmgregion *r, const struct bn_tol *tol)
Definition: nmg_mk.c:2557