BRL-CAD
ell.c
Go to the documentation of this file.
1 /* E L L . 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/ell/ell.c
23  *
24  * Intersect a ray with a Generalized Ellipsoid.
25  *
26  */
27 
28 #include "common.h"
29 
30 #include <stddef.h>
31 #include <string.h>
32 #include <math.h>
33 #include "bio.h"
34 
35 #include "bu/cv.h"
36 #include "vmath.h"
37 #include "db.h"
38 #include "nmg.h"
39 #include "rtgeom.h"
40 #include "raytrace.h"
41 #include "nurb.h"
42 
43 #include "../../librt_private.h"
44 
45 
46 extern int rt_sph_prep(struct soltab *stp, struct rt_db_internal *ip,
47  struct rt_i *rtip);
48 
49 const struct bu_structparse rt_ell_parse[] = {
50  { "%f", 3, "V", bu_offsetofarray(struct rt_ell_internal, v, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
51  { "%f", 3, "A", bu_offsetofarray(struct rt_ell_internal, a, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
52  { "%f", 3, "B", bu_offsetofarray(struct rt_ell_internal, b, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
53  { "%f", 3, "C", bu_offsetofarray(struct rt_ell_internal, c, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
54  { {'\0', '\0', '\0', '\0'}, 0, (char *)NULL, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL }
55 };
56 
57 
58 HIDDEN void nmg_sphere_face_snurb(struct faceuse *fu, const matp_t m);
59 
60 /*
61  * Algorithm:
62  *
63  * Given V, A, B, and C, there is a set of points on this ellipsoid
64  *
65  * { (x, y, z) | (x, y, z) is on ellipsoid defined by V, A, B, C }
66  *
67  * Through a series of Affine Transformations, this set will be
68  * transformed into a set of points on a unit sphere at the origin
69  *
70  * { (x', y', z') | (x', y', z') is on Sphere at origin }
71  *
72  * The transformation from X to X' is accomplished by:
73  *
74  * X' = S(R(X - V))
75  *
76  * where R(X) = (A/(|A|))
77  * (B/(|B|)) . X
78  * (C/(|C|))
79  *
80  * and S(X) = (1/|A| 0 0)
81  * (0 1/|B| 0) . X
82  * (0 0 1/|C|)
83  *
84  * To find the intersection of a line with the ellipsoid, consider the
85  * parametric line L:
86  *
87  * L : { P(n) | P + t(n) . D }
88  *
89  * Call W the actual point of intersection between L and the
90  * ellipsoid. Let W' be the point of intersection between L' and the
91  * unit sphere.
92  *
93  * L' : { P'(n) | P' + t(n) . D' }
94  *
95  * W = invR(invS(W')) + V
96  *
97  * Where W' = k D' + P'.
98  *
99  * Let dp = D' dot P'
100  * Let dd = D' dot D'
101  * Let pp = P' dot P'
102  *
103  * and k = [ -dp +/- sqrt(dp*dp - dd * (pp - 1)) ] / dd
104  * which is constant.
105  *
106  * Now, D' = S(R(D))
107  * and P' = S(R(P - V))
108  *
109  * Substituting,
110  *
111  * W = V + invR(invS[ k *(S(R(D))) + S(R(P - V)) ])
112  * = V + invR((k * R(D)) + R(P - V))
113  * = V + k * D + P - V
114  * = k * D + P
115  *
116  * Note that ``k'' is constant, and is the same in the formulations
117  * for both W and W'.
118  *
119  * NORMALS. Given the point W on the ellipsoid, what is the vector
120  * normal to the tangent plane at that point?
121  *
122  * Map W onto the unit sphere, i.e.: W' = S(R(W - V)).
123  *
124  * Plane on unit sphere at W' has a normal vector of the same
125  * value(!). N' = W'
126  *
127  * The plane transforms back to the tangent plane at W, and this new
128  * plane (on the ellipsoid) has a normal vector of N, viz:
129  *
130  * N = inverse[ transpose(inverse[ S o R ]) ] (N')
131  *
132  * because if H is perpendicular to plane Q, and matrix M maps from Q
133  * to Q', then inverse[ transpose(M) ] (H) is perpendicular to Q'.
134  * Here, H and Q are in "prime space" with the unit sphere. [Somehow,
135  * the notation here is backwards]. So, the mapping matrix M =
136  * inverse(S o R), because S o R maps from normal space to the unit
137  * sphere.
138  *
139  * N = inverse[ transpose(inverse[ S o R ]) ] (N')
140  * = inverse[ transpose(invR o invS) ] (N')
141  * = inverse[ transpose(invS) o transpose(invR) ] (N')
142  * = inverse[ inverse(S) o R ] (N')
143  * = invR o S (N')
144  *
145  * = invR o S (W')
146  * = invR(S(S(R(W - V))))
147  *
148  * because inverse(R) = transpose(R), so R = transpose(invR),
149  * and S = transpose(S).
150  *
151  * Note that the normal vector N produced above will not have unit
152  * length.
153  */
154 
155 struct ell_specific {
156  vect_t ell_V; /* Vector to center of ellipsoid */
157  vect_t ell_Au; /* unit-length A vector */
158  vect_t ell_Bu;
159  vect_t ell_Cu;
160  vect_t ell_invsq; /* [ 1/(|A|**2), 1/(|B|**2), 1/(|C|**2) ] */
161  mat_t ell_SoR; /* Scale(Rot(vect)) */
162  mat_t ell_invRSSR; /* invRot(Scale(Scale(Rot(vect)))) */
163 };
164 
165 
166 #define ELL_NULL ((struct ell_specific *)0)
167 
168 /**
169  * Compute the bounding RPP for an ellipsoid
170  */
171 int
172 rt_ell_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol)) {
173  vect_t w1, w2, P;
174  vect_t Au, Bu, Cu; /* A, B, C with unit length */
175  fastf_t magsq_a, magsq_b, magsq_c, f;
176  mat_t R;
177  struct rt_ell_internal *eip = (struct rt_ell_internal *)ip->idb_ptr;
178  RT_ELL_CK_MAGIC(eip);
179 
180  magsq_a = MAGSQ(eip->a);
181  magsq_b = MAGSQ(eip->b);
182  magsq_c = MAGSQ(eip->c);
183 
184  /* Try a shortcut - if this is a sphere, the calculation simplifies */
185  /* Check whether |A|, |B|, and |C| are nearly equal */
186  if (EQUAL(magsq_a, magsq_b) && EQUAL(magsq_a, magsq_c)) {
187  fastf_t sph_rad = sqrt(magsq_a);
188  (*min)[X] = eip->v[X] - sph_rad;
189  (*max)[X] = eip->v[X] + sph_rad;
190  (*min)[Y] = eip->v[Y] - sph_rad;
191  (*max)[Y] = eip->v[Y] + sph_rad;
192  (*min)[Z] = eip->v[Z] - sph_rad;
193  (*max)[Z] = eip->v[Z] + sph_rad;
194  return 0;
195  }
196 
197  /* Create unit length versions of A, B, C */
198  f = 1.0/sqrt(magsq_a);
199  VSCALE(Au, eip->a, f);
200  f = 1.0/sqrt(magsq_b);
201  VSCALE(Bu, eip->b, f);
202  f = 1.0/sqrt(magsq_c);
203  VSCALE(Cu, eip->c, f);
204 
205  MAT_IDN(R);
206  VMOVE(&R[0], Au);
207  VMOVE(&R[4], Bu);
208  VMOVE(&R[8], Cu);
209 
210  /* Compute bounding RPP */
211  VSET(w1, magsq_a, magsq_b, magsq_c);
212 
213  /* X */
214  VSET(P, 1.0, 0, 0); /* bounding plane normal */
215  MAT3X3VEC(w2, R, P); /* map plane to local coord syst */
216  VELMUL(w2, w2, w2); /* square each term */
217  f = VDOT(w1, w2);
218  f = sqrt(f);
219  (*min)[X] = eip->v[X] - f; /* V.P +/- f */
220  (*max)[X] = eip->v[X] + f;
221 
222  /* Y */
223  VSET(P, 0, 1.0, 0); /* bounding plane normal */
224  MAT3X3VEC(w2, R, P); /* map plane to local coord syst */
225  VELMUL(w2, w2, w2); /* square each term */
226  f = VDOT(w1, w2);
227  f = sqrt(f);
228  (*min)[Y] = eip->v[Y] - f; /* V.P +/- f */
229  (*max)[Y] = eip->v[Y] + f;
230 
231  /* Z */
232  VSET(P, 0, 0, 1.0); /* bounding plane normal */
233  MAT3X3VEC(w2, R, P); /* map plane to local coord syst */
234  VELMUL(w2, w2, w2); /* square each term */
235  f = VDOT(w1, w2);
236  f = sqrt(f);
237  (*min)[Z] = eip->v[Z] - f; /* V.P +/- f */
238  (*max)[Z] = eip->v[Z] + f;
239  return 0;
240 }
241 
242 
243 /**
244  * Given a pointer to a GED database record, and a transformation
245  * matrix, determine if this is a valid ellipsoid, and if so,
246  * precompute various terms of the formula.
247  *
248  * Returns -
249  * 0 ELL is OK
250  * !0 Error in description
251  *
252  * Implicit return -
253  * A struct ell_specific is created, and its address is stored in
254  * stp->st_specific for use by rt_ell_shot().
255  */
256 int
257 rt_ell_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
258 {
259  register struct ell_specific *ell;
260  struct rt_ell_internal *eip;
261  fastf_t magsq_a, magsq_b, magsq_c;
262  mat_t R;
263  mat_t Rinv;
264  mat_t SS;
265  mat_t mtemp;
266  vect_t Au, Bu, Cu; /* A, B, C with unit length */
267  fastf_t f;
268 
269  eip = (struct rt_ell_internal *)ip->idb_ptr;
270  RT_ELL_CK_MAGIC(eip);
271 
272  /*
273  * For a fast way out, hand this solid off to the SPH routine. If
274  * it takes it, then there is nothing to do, otherwise the solid
275  * is an ELL.
276  */
277  if (rt_sph_prep(stp, ip, rtip) == 0) {
278  return 0; /* OK */
279  }
280 
281  /* Validate that |A| > 0, |B| > 0, |C| > 0 */
282  magsq_a = MAGSQ(eip->a);
283  magsq_b = MAGSQ(eip->b);
284  magsq_c = MAGSQ(eip->c);
285 
286  if (magsq_a < rtip->rti_tol.dist_sq || magsq_b < rtip->rti_tol.dist_sq || magsq_c < rtip->rti_tol.dist_sq) {
287  bu_log("rt_ell_prep(): ell(%s) zero length A(%g), B(%g), or C(%g) vector\n",
288  stp->st_name, magsq_a, magsq_b, magsq_c);
289  return 1; /* BAD */
290  }
291 
292  /* Create unit length versions of A, B, C */
293  f = 1.0/sqrt(magsq_a);
294  VSCALE(Au, eip->a, f);
295  f = 1.0/sqrt(magsq_b);
296  VSCALE(Bu, eip->b, f);
297  f = 1.0/sqrt(magsq_c);
298  VSCALE(Cu, eip->c, f);
299 
300  /* Validate that A.B == 0, B.C == 0, A.C == 0 (check dir only) */
301  f = VDOT(Au, Bu);
302  if (! NEAR_ZERO(f, rtip->rti_tol.dist)) {
303  bu_log("rt_ell_prep(): ell(%s) A not perpendicular to B, f=%f\n", stp->st_name, f);
304  return 1; /* BAD */
305  }
306  f = VDOT(Bu, Cu);
307  if (! NEAR_ZERO(f, rtip->rti_tol.dist)) {
308  bu_log("rt_ell_prep(): ell(%s) B not perpendicular to C, f=%f\n", stp->st_name, f);
309  return 1; /* BAD */
310  }
311  f = VDOT(Au, Cu);
312  if (! NEAR_ZERO(f, rtip->rti_tol.dist)) {
313  bu_log("rt_ell_prep(): ell(%s) A not perpendicular to C, f=%f\n", stp->st_name, f);
314  return 1; /* BAD */
315  }
316 
317  /* Solid is OK, compute constant terms now */
318  BU_GET(ell, struct ell_specific);
319  stp->st_specific = (void *)ell;
320 
321  VMOVE(ell->ell_V, eip->v);
322 
323  VSET(ell->ell_invsq, 1.0/magsq_a, 1.0/magsq_b, 1.0/magsq_c);
324  VMOVE(ell->ell_Au, Au);
325  VMOVE(ell->ell_Bu, Bu);
326  VMOVE(ell->ell_Cu, Cu);
327 
328  MAT_IDN(ell->ell_SoR);
329  MAT_IDN(R);
330 
331  /* Compute R and Rinv matrices */
332  VMOVE(&R[0], Au);
333  VMOVE(&R[4], Bu);
334  VMOVE(&R[8], Cu);
335  bn_mat_trn(Rinv, R); /* inv of rot mat is trn */
336 
337  /* Compute SoS (Affine transformation) */
338  MAT_IDN(SS);
339  SS[ 0] = ell->ell_invsq[0];
340  SS[ 5] = ell->ell_invsq[1];
341  SS[10] = ell->ell_invsq[2];
342 
343  /* Compute invRSSR */
344  bn_mat_mul(mtemp, SS, R);
345  bn_mat_mul(ell->ell_invRSSR, Rinv, mtemp);
346 
347  /* Compute SoR */
348  VSCALE(&ell->ell_SoR[0], eip->a, ell->ell_invsq[0]);
349  VSCALE(&ell->ell_SoR[4], eip->b, ell->ell_invsq[1]);
350  VSCALE(&ell->ell_SoR[8], eip->c, ell->ell_invsq[2]);
351 
352  /* Compute bounding sphere */
353  VMOVE(stp->st_center, eip->v);
354  f = magsq_a;
355  if (magsq_b > f)
356  f = magsq_b;
357  if (magsq_c > f)
358  f = magsq_c;
359  stp->st_aradius = stp->st_bradius = sqrt(f);
360 
361  /* Compute bounding RPP */
362  if (stp->st_meth->ft_bbox(ip, &(stp->st_min), &(stp->st_max), &(rtip->rti_tol))) return 1;
363  return 0; /* OK */
364 }
365 
366 
367 void
368 rt_ell_print(register const struct soltab *stp)
369 {
370  register struct ell_specific *ell =
371  (struct ell_specific *)stp->st_specific;
372 
373  VPRINT("V", ell->ell_V);
374  bn_mat_print("S o R", ell->ell_SoR);
375  bn_mat_print("invRSSR", ell->ell_invRSSR);
376 }
377 
378 
379 /**
380  * Intersect a ray with an ellipsoid, where all constant terms have
381  * been precomputed by rt_ell_prep(). If an intersection occurs, a
382  * struct seg will be acquired and filled in.
383  *
384  * Returns -
385  * 0 MISS
386  * >0 HIT
387  */
388 int
389 rt_ell_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
390 {
391  register struct ell_specific *ell =
392  (struct ell_specific *)stp->st_specific;
393  register struct seg *segp;
394 
395  vect_t dprime; /* D' */
396  vect_t pprime; /* P' */
397  vect_t xlated; /* translated vector */
398  fastf_t dp, dd; /* D' dot P', D' dot D' */
399  fastf_t k1, k2; /* distance constants of solution */
400  fastf_t root; /* root of radical */
401 
402  /* out, Mat, vect */
403  MAT4X3VEC(dprime, ell->ell_SoR, rp->r_dir);
404  VSUB2(xlated, rp->r_pt, ell->ell_V);
405  MAT4X3VEC(pprime, ell->ell_SoR, xlated);
406 
407  dp = VDOT(dprime, pprime);
408  dd = VDOT(dprime, dprime);
409 
410  if ((root = dp*dp - dd * (VDOT(pprime, pprime)-1.0)) < 0)
411  return 0; /* No hit */
412  root = sqrt(root);
413 
414  RT_GET_SEG(segp, ap->a_resource);
415  segp->seg_stp = stp;
416  if ((k1=(-dp+root)/dd) <= (k2=(-dp-root)/dd)) {
417  /* k1 is entry, k2 is exit */
418  segp->seg_in.hit_dist = k1;
419  segp->seg_out.hit_dist = k2;
420  } else {
421  /* k2 is entry, k1 is exit */
422  segp->seg_in.hit_dist = k2;
423  segp->seg_out.hit_dist = k1;
424  }
425  segp->seg_in.hit_surfno = 0;
426  segp->seg_out.hit_surfno = 0;
427  BU_LIST_INSERT(&(seghead->l), &(segp->l));
428  return 2; /* HIT */
429 }
430 
431 
432 #define RT_ELL_SEG_MISS(SEG) (SEG).seg_stp=RT_SOLTAB_NULL
433 /**
434  * This is the Becker vector version.
435  */
436 void
437 rt_ell_vshot(struct soltab **stp, struct xray **rp, struct seg *segp, int n, struct application *ap)
438 /* An array of solid pointers */
439 /* An array of ray pointers */
440 /* array of segs (results returned) */
441 /* Number of ray/object pairs */
442 
443 {
444  register int i;
445  register struct ell_specific *ell;
446 
447  vect_t dprime; /* D' */
448  vect_t pprime; /* P' */
449  vect_t xlated; /* translated vector */
450  fastf_t dp, dd; /* D' dot P', D' dot D' */
451  fastf_t k1, k2; /* distance constants of solution */
452  fastf_t root; /* root of radical */
453 
454  if (ap) RT_CK_APPLICATION(ap);
455 
456  /* for each ray/ellipse pair */
457  for (i = 0; i < n; i++) {
458  if (stp[i] == 0) continue; /* stp[i] == 0 signals skip ray */
459 
460  ell = (struct ell_specific *)stp[i]->st_specific;
461 
462  MAT4X3VEC(dprime, ell->ell_SoR, rp[i]->r_dir);
463  VSUB2(xlated, rp[i]->r_pt, ell->ell_V);
464  MAT4X3VEC(pprime, ell->ell_SoR, xlated);
465 
466  dp = VDOT(dprime, pprime);
467  dd = VDOT(dprime, dprime);
468 
469  if ((root = dp*dp - dd * (VDOT(pprime, pprime)-1.0)) < 0) {
470  RT_ELL_SEG_MISS(segp[i]); /* No hit */
471  } else {
472  root = sqrt(root);
473 
474  segp[i].seg_stp = stp[i];
475 
476  if ((k1=(-dp+root)/dd) <= (k2=(-dp-root)/dd)) {
477  /* k1 is entry, k2 is exit */
478  segp[i].seg_in.hit_dist = k1;
479  segp[i].seg_out.hit_dist = k2;
480  } else {
481  /* k2 is entry, k1 is exit */
482  segp[i].seg_in.hit_dist = k2;
483  segp[i].seg_out.hit_dist = k1;
484  }
485  segp[i].seg_in.hit_surfno = 0;
486  segp[i].seg_out.hit_surfno = 0;
487  }
488  }
489 }
490 
491 
492 /**
493  * Given ONE ray distance, return the normal and entry/exit point.
494  */
495 void
496 rt_ell_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
497 {
498  register struct ell_specific *ell =
499  (struct ell_specific *)stp->st_specific;
500  vect_t xlated;
501  fastf_t scale;
502 
503  VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
504  VSUB2(xlated, hitp->hit_point, ell->ell_V);
505  MAT4X3VEC(hitp->hit_normal, ell->ell_invRSSR, xlated);
506  scale = 1.0 / MAGNITUDE(hitp->hit_normal);
507  VSCALE(hitp->hit_normal, hitp->hit_normal, scale);
508 
509  /* tuck away this scale for the curvature routine */
510  hitp->hit_vpriv[X] = scale;
511 }
512 
513 
514 /**
515  * Return the curvature of the ellipsoid.
516  */
517 void
518 rt_ell_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
519 {
520  register struct ell_specific *ell =
521  (struct ell_specific *)stp->st_specific;
522  vect_t u, v; /* basis vectors (with normal) */
523  vect_t vec1, vec2; /* eigen vectors */
524  vect_t tmp;
525  fastf_t a, b, c, scale;
526 
527  /*
528  * choose a tangent plane coordinate system (u, v, normal) form a
529  * right-handed triple
530  */
531  bn_vec_ortho(u, hitp->hit_normal);
532  VCROSS(v, hitp->hit_normal, u);
533 
534  /* get the saved away scale factor */
535  scale = - hitp->hit_vpriv[X];
536 
537  /* find the second fundamental form */
538  MAT4X3VEC(tmp, ell->ell_invRSSR, u);
539  a = VDOT(u, tmp) * scale;
540  b = VDOT(v, tmp) * scale;
541  MAT4X3VEC(tmp, ell->ell_invRSSR, v);
542  c = VDOT(v, tmp) * scale;
543 
544  bn_eigen2x2(&cvp->crv_c1, &cvp->crv_c2, vec1, vec2, a, b, c);
545  VCOMB2(cvp->crv_pdir, vec1[X], u, vec1[Y], v);
546  VUNITIZE(cvp->crv_pdir);
547 }
548 
549 
550 /**
551  * For a hit on the surface of an ELL, return the (u, v) coordinates
552  * of the hit point, 0 <= u, v <= 1.
553  *
554  * u = azimuth
555  * v = elevation
556  */
557 void
558 rt_ell_uv(struct application *ap, struct soltab *stp, register struct hit *hitp, register struct uvcoord *uvp)
559 {
560  register struct ell_specific *ell =
561  (struct ell_specific *)stp->st_specific;
562  vect_t work;
563  vect_t pprime;
564  fastf_t r;
565 
566  /* hit_point is on surface; project back to unit sphere, creating
567  * a vector from vertex to hit point which always has length=1.0
568  */
569  VSUB2(work, hitp->hit_point, ell->ell_V);
570  MAT4X3VEC(pprime, ell->ell_SoR, work);
571  /* Assert that pprime has unit length */
572 
573  /* U is azimuth, atan() range: -pi to +pi */
574  uvp->uv_u = bn_atan2(pprime[Y], pprime[X]) * M_1_2PI;
575  if (uvp->uv_u < 0)
576  uvp->uv_u += 1.0;
577  /* V is elevation, atan() range: -pi/2 to +pi/2, because sqrt()
578  * ensures that X parameter is always >0
579  */
580  uvp->uv_v = bn_atan2(pprime[Z],
581  sqrt(pprime[X] * pprime[X] + pprime[Y] * pprime[Y])) *
582  M_1_PI + 0.5;
583 
584  /* approximation: r / (circumference, 2 * pi * aradius) */
585  r = ap->a_rbeam + ap->a_diverge * hitp->hit_dist;
586  uvp->uv_du = uvp->uv_dv =
587  M_1_2PI * r / stp->st_aradius;
588 }
589 
590 
591 void
592 rt_ell_free(register struct soltab *stp)
593 {
594  register struct ell_specific *ell =
595  (struct ell_specific *)stp->st_specific;
596 
597  BU_PUT(ell, struct ell_specific);
598 }
599 
600 
601 /**
602  * Also used by the TGC code
603  */
604 #define ELLOUT(n) ov+(n-1)*3
605 void
607  fastf_t *V,
608  fastf_t *A,
609  fastf_t *B)
610 {
611  static fastf_t c, d, e, f, g, h;
612 
613  e = h = .92388; /* cos(22.5 degrees) */
614  c = d = M_SQRT1_2; /* cos(45 degrees) */
615  g = f = .382683; /* cos(67.5 degrees) */
616 
617  /*
618  * For angle theta, compute surface point as
619  *
620  * V + cos(theta) * A + sin(theta) * B
621  *
622  * note that sin(theta) is cos(90-theta), with arguments in degrees.
623  */
624 
625  VADD2(ELLOUT(1), V, A);
626  VJOIN2(ELLOUT(2), V, e, A, f, B);
627  VJOIN2(ELLOUT(3), V, c, A, d, B);
628  VJOIN2(ELLOUT(4), V, g, A, h, B);
629  VADD2(ELLOUT(5), V, B);
630  VJOIN2(ELLOUT(6), V, -g, A, h, B);
631  VJOIN2(ELLOUT(7), V, -c, A, d, B);
632  VJOIN2(ELLOUT(8), V, -e, A, f, B);
633  VSUB2(ELLOUT(9), V, A);
634  VJOIN2(ELLOUT(10), V, -e, A, -f, B);
635  VJOIN2(ELLOUT(11), V, -c, A, -d, B);
636  VJOIN2(ELLOUT(12), V, -g, A, -h, B);
637  VSUB2(ELLOUT(13), V, B);
638  VJOIN2(ELLOUT(14), V, g, A, -h, B);
639  VJOIN2(ELLOUT(15), V, c, A, -d, B);
640  VJOIN2(ELLOUT(16), V, e, A, -f, B);
641 }
642 
644  struct bu_list *vhead;
645  vect_t ell_center;
651 };
652 
654  vect_t a;
655  vect_t b;
656  vect_t translation;
659 };
660 
661 static double
662 cross_section_axis_magnitude(
663  struct ellipse_cross_section cross_section,
664  double target_axis_mag)
665 {
666  double travel_mag, travel_pos;
667  double fixed_term;
668 
669  travel_mag = cross_section.ellipsoid_travel_axis_mag;
670  travel_pos = cross_section.ellipsoid_travel_axis_position;
671 
672  fixed_term = (travel_pos * travel_pos) / (travel_mag * travel_mag);
673 
674  return sqrt(target_axis_mag * target_axis_mag * (1.0 - fixed_term));
675 }
676 
677 /* Draws elliptical cross-sections along an ell vector (the "travel vector").
678  *
679  * algorithm overview:
680  * 1 Pretend the ellipsoid is at the origin with the ell vectors axis-aligned.
681  * 2 Step along the travel vector t from -mag(t) to mag(t).
682  * 3 Calculate the axis vectors for the cross-section ellipse at that step.
683  * 4 Draw multiple points along the ellipse to define the cross-section,
684  * using a parametric vector equation to calculate ellipse points for
685  * given radian values.
686  */
687 static void
688 draw_cross_sections_along_ell_vector(struct ell_draw_configuration config)
689 {
690  int i;
691  vect_t ell_a, ell_b, ell_t;
692  double ellipse_a_mag, ellipse_b_mag;
693  double ell_a_mag, ell_b_mag, ell_t_mag;
694  double num_cross_sections, points_per_section;
695  double position_step;
696  struct ellipse_cross_section cross_section = {VINIT_ZERO, VINIT_ZERO, VINIT_ZERO, 0.0, 0.0};
697 
698  VMOVE(ell_a, config.ell_axis_vector_a);
699  VMOVE(ell_b, config.ell_axis_vector_b);
700  VMOVE(ell_t, config.ell_travel_vector);
701 
702  ell_a_mag = MAGNITUDE(ell_a);
703  ell_b_mag = MAGNITUDE(ell_b);
704  ell_t_mag = MAGNITUDE(ell_t);
705 
706  points_per_section = config.points_per_section;
707  num_cross_sections = config.num_cross_sections;
708 
709  position_step = 2.0 * ell_t_mag / (num_cross_sections + 1);
710 
711  cross_section.ellipsoid_travel_axis_mag = ell_t_mag;
712  cross_section.ellipsoid_travel_axis_position = -ell_t_mag;
713  for (i = 0; i < num_cross_sections; ++i) {
714  cross_section.ellipsoid_travel_axis_position += position_step;
715 
716  ellipse_a_mag = cross_section_axis_magnitude(cross_section, ell_a_mag);
717  ellipse_b_mag = cross_section_axis_magnitude(cross_section, ell_b_mag);
718 
719  VSCALE(cross_section.a, ell_a, ellipse_a_mag / ell_a_mag);
720  VSCALE(cross_section.b, ell_b, ellipse_b_mag / ell_b_mag);
721  VSCALE(cross_section.translation, ell_t,
722  cross_section.ellipsoid_travel_axis_position / ell_t_mag);
723  VADD2(cross_section.translation, cross_section.translation, config.ell_center);
724 
725  plot_ellipse(config.vhead, cross_section.translation, cross_section.a,
726  cross_section.b, points_per_section);
727  }
728 }
729 
730 /* decide how many ellipse points are needed to satisfy point spacing */
731 static int
732 ell_ellipse_points(
733  const struct rt_ell_internal *ell,
734  const struct rt_view_info *info)
735 {
736  fastf_t avg_radius, avg_circumference;
737  fastf_t ell_mag_a, ell_mag_b, ell_mag_c;
738 
739  ell_mag_a = MAGNITUDE(ell->a);
740  ell_mag_b = MAGNITUDE(ell->b);
741  ell_mag_c = MAGNITUDE(ell->c);
742 
743  avg_radius = (ell_mag_a + ell_mag_b + ell_mag_c) / 3.0;
744  avg_circumference = M_2PI * avg_radius;
745 
746  return avg_circumference / info->point_spacing;
747 }
748 
749 int
750 rt_ell_adaptive_plot(struct rt_db_internal *ip, const struct rt_view_info *info)
751 {
752  struct ell_draw_configuration config;
753  struct rt_ell_internal *eip;
754 
755  BU_CK_LIST_HEAD(info->vhead);
756  RT_CK_DB_INTERNAL(ip);
757  eip = (struct rt_ell_internal *)ip->idb_ptr;
758  RT_ELL_CK_MAGIC(eip);
759 
760  config.vhead = info->vhead;
761  VMOVE(config.ell_center, eip->v);
762 
763  config.points_per_section = ell_ellipse_points(eip, info);
764 
765  if (config.points_per_section < 4) {
766  RT_ADD_VLIST(info->vhead, eip->v, BN_VLIST_POINT_DRAW);
767  return 0;
768  }
769 
770  config.num_cross_sections = primitive_curve_count(ip, info);
771 
772  VMOVE(config.ell_travel_vector, eip->a);
773  VMOVE(config.ell_axis_vector_a, eip->b);
774  VMOVE(config.ell_axis_vector_b, eip->c);
775  draw_cross_sections_along_ell_vector(config);
776 
777  config.num_cross_sections = 1;
778 
779  VMOVE(config.ell_travel_vector, eip->b);
780  VMOVE(config.ell_axis_vector_a, eip->a);
781  VMOVE(config.ell_axis_vector_b, eip->c);
782  draw_cross_sections_along_ell_vector(config);
783 
784  VMOVE(config.ell_travel_vector, eip->c);
785  VMOVE(config.ell_axis_vector_a, eip->a);
786  VMOVE(config.ell_axis_vector_b, eip->b);
787  draw_cross_sections_along_ell_vector(config);
788 
789  return 0;
790 }
791 
792 int
793 rt_ell_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))
794 {
795  register int i;
796  struct rt_ell_internal *eip;
797  fastf_t top[16*3];
798  fastf_t middle[16*3];
799  fastf_t bottom[16*3];
800 
801  BU_CK_LIST_HEAD(vhead);
802  RT_CK_DB_INTERNAL(ip);
803  eip = (struct rt_ell_internal *)ip->idb_ptr;
804  RT_ELL_CK_MAGIC(eip);
805 
806  rt_ell_16pts(top, eip->v, eip->a, eip->b);
807  rt_ell_16pts(bottom, eip->v, eip->b, eip->c);
808  rt_ell_16pts(middle, eip->v, eip->a, eip->c);
809 
810  RT_ADD_VLIST(vhead, &top[15*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
811  for (i = 0; i < 16; i++) {
812  RT_ADD_VLIST(vhead, &top[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
813  }
814 
815  RT_ADD_VLIST(vhead, &bottom[15*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
816  for (i = 0; i < 16; i++) {
817  RT_ADD_VLIST(vhead, &bottom[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
818  }
819 
820  RT_ADD_VLIST(vhead, &middle[15*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
821  for (i = 0; i < 16; i++) {
822  RT_ADD_VLIST(vhead, &middle[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
823  }
824 
825  return 0;
826 }
827 
828 
829 /* Vertices of a unit octahedron. These need to be organized properly
830  * to give reasonable normals.
831  *
832  * static struct usvert {
833  * int a;
834  * int b;
835  * int c;
836  * } octahedron[8] = {
837  * { XPLUS, ZPLUS, YPLUS },
838  * { YPLUS, ZPLUS, XMIN },
839  * { XMIN , ZPLUS, YMIN },
840  * { YMIN , ZPLUS, XPLUS },
841  * { XPLUS, YPLUS, ZMIN },
842  * { YPLUS, XMIN , ZMIN },
843  * { XMIN , YMIN , ZMIN },
844  * { YMIN , XPLUS, ZMIN }
845  * };
846  */
847 struct ell_state {
848  struct shell *s;
849  struct rt_ell_internal *eip;
850  mat_t invRinvS;
851  mat_t invRoS;
853 };
854 
855 
858  int nverts;
859  struct vertex **vp;
860  vect_t *norms;
861  int nfaces;
862  struct faceuse **fu;
863 };
864 
865 
866 /**
867  * Tessellate an ellipsoid.
868  *
869  * The strategy is based upon the approach of Jon Leech 3/24/89, from
870  * program "sphere", which generates a polygon mesh approximating a
871  * sphere by recursive subdivision. First approximation is an
872  * octahedron; each level of refinement increases the number of
873  * polygons by a factor of 4. Level 3 (128 polygons) is a good
874  * tradeoff if gouraud shading is used to render the database.
875  *
876  * At the start, points ABC lie on surface of the unit sphere. Pick
877  * DEF as the midpoints of the three edges of ABC. Normalize the new
878  * points to lie on surface of the unit sphere.
879  *
880  * 1
881  * B
882  * /\
883  * 3 / \ 4
884  * D /____\ E
885  * /\ /\
886  * / \ / \
887  * /____\/____\
888  * A F C
889  * 0 5 2
890  *
891  * Returns -
892  * -1 failure
893  * 0 OK. *r points to nmgregion that holds this tessellation.
894  */
895 int
896 rt_ell_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
897 {
898  mat_t R;
899  mat_t S;
900  mat_t invR;
901  mat_t invS;
902  vect_t Au, Bu, Cu; /* A, B, C with unit length */
903  fastf_t Alen, Blen, Clen;
904  fastf_t invAlen, invBlen, invClen;
905  fastf_t magsq_a, magsq_b, magsq_c;
906  fastf_t f;
907  struct ell_state state;
908  register int i;
909  fastf_t radius;
910  int nsegs;
911  int nstrips;
912  struct ell_vert_strip *strips;
913  int j;
914  struct vertex **vertp[4];
915  int faceno;
916  int stripno;
917  int boff; /* base offset */
918  int toff; /* top offset */
919  int blim; /* base subscript limit */
920  int tlim; /* top subscript limit */
921  fastf_t dtol; /* Absolutized relative tolerance */
922 
923  RT_CK_DB_INTERNAL(ip);
924  state.eip = (struct rt_ell_internal *)ip->idb_ptr;
925  RT_ELL_CK_MAGIC(state.eip);
926 
927  /* Validate that |A| > 0, |B| > 0, |C| > 0 */
928  magsq_a = MAGSQ(state.eip->a);
929  magsq_b = MAGSQ(state.eip->b);
930  magsq_c = MAGSQ(state.eip->c);
931  if (magsq_a < tol->dist_sq || magsq_b < tol->dist_sq || magsq_c < tol->dist_sq) {
932  bu_log("rt_ell_tess(): zero length A, B, or C vector\n");
933  return -2; /* BAD */
934  }
935 
936  /* Create unit length versions of A, B, C */
937  invAlen = 1.0/(Alen = sqrt(magsq_a));
938  VSCALE(Au, state.eip->a, invAlen);
939  invBlen = 1.0/(Blen = sqrt(magsq_b));
940  VSCALE(Bu, state.eip->b, invBlen);
941  invClen = 1.0/(Clen = sqrt(magsq_c));
942  VSCALE(Cu, state.eip->c, invClen);
943 
944  /* Validate that A.B == 0, B.C == 0, A.C == 0 (check dir only) */
945  f = VDOT(Au, Bu);
946  if (! NEAR_ZERO(f, tol->dist)) {
947  bu_log("rt_ell_tess(): A not perpendicular to B, f=%f\n", f);
948  return -3; /* BAD */
949  }
950  f = VDOT(Bu, Cu);
951  if (! NEAR_ZERO(f, tol->dist)) {
952  bu_log("rt_ell_tess(): B not perpendicular to C, f=%f\n", f);
953  return -3; /* BAD */
954  }
955  f = VDOT(Au, Cu);
956  if (! NEAR_ZERO(f, tol->dist)) {
957  bu_log("rt_ell_tess(): A not perpendicular to C, f=%f\n", f);
958  return -3; /* BAD */
959  }
960 
961  {
962  vect_t axb;
963  VCROSS(axb, Au, Bu);
964  f = VDOT(axb, Cu);
965  if (f < 0) {
966  VREVERSE(Cu, Cu);
967  VREVERSE(state.eip->c, state.eip->c);
968  }
969  }
970 
971  /* Compute R and Rinv matrices */
972  MAT_IDN(R);
973  VMOVE(&R[0], Au);
974  VMOVE(&R[4], Bu);
975  VMOVE(&R[8], Cu);
976  bn_mat_trn(invR, R); /* inv of rot mat is trn */
977 
978  /* Compute S and invS matrices */
979  /* invS is just 1/diagonal elements */
980  MAT_IDN(S);
981  S[ 0] = invAlen;
982  S[ 5] = invBlen;
983  S[10] = invClen;
984  bn_mat_inv(invS, S);
985 
986  /* invRinvS, for converting points from unit sphere to model */
987  bn_mat_mul(state.invRinvS, invR, invS);
988 
989  /* invRoS, for converting normals from unit sphere to model */
990  bn_mat_mul(state.invRoS, invR, S);
991 
992  /* Compute radius of bounding sphere */
993  radius = Alen;
994  if (Blen > radius)
995  radius = Blen;
996  if (Clen > radius)
997  radius = Clen;
998 
999  dtol = primitive_get_absolute_tolerance(ttol, radius);
1000 
1001  if (dtol > radius) {
1002  dtol = radius;
1003  }
1004 
1005  /* Convert distance tolerance into a maximum permissible angle
1006  * tolerance. 'radius' is largest radius.
1007  */
1008  state.theta_tol = 2 * acos(1.0 - dtol / radius);
1009 
1010  /* To ensure normal tolerance, remain below this angle */
1011  if (ttol->norm > 0.0 && ttol->norm < state.theta_tol) {
1012  state.theta_tol = ttol->norm;
1013  }
1014 
1015  *r = nmg_mrsv(m); /* Make region, empty shell, vertex */
1016  state.s = BU_LIST_FIRST(shell, &(*r)->s_hd);
1017 
1018  /* Find the number of segments to divide 90 degrees worth into */
1019  nsegs = (int)(M_PI_2 / state.theta_tol + 0.999);
1020  if (nsegs < 2) nsegs = 2;
1021 
1022  /* Find total number of strips of vertices that will be needed.
1023  * nsegs for each hemisphere, plus the equator. Note that faces
1024  * are listed in the stripe ABOVE, i.e., toward the poles.
1025  * Thus, strips[0] will have 4 faces.
1026  */
1027  nstrips = 2 * nsegs + 1;
1028  strips = (struct ell_vert_strip *)bu_calloc(nstrips,
1029  sizeof(struct ell_vert_strip), "strips[]");
1030 
1031  /* North pole */
1032  strips[0].nverts = 1;
1033  strips[0].nverts_per_strip = 0;
1034  strips[0].nfaces = 4;
1035  /* South pole */
1036  strips[nstrips-1].nverts = 1;
1037  strips[nstrips-1].nverts_per_strip = 0;
1038  strips[nstrips-1].nfaces = 4;
1039  /* equator */
1040  strips[nsegs].nverts = nsegs * 4;
1041  strips[nsegs].nverts_per_strip = nsegs;
1042  strips[nsegs].nfaces = 0;
1043 
1044  for (i = 1; i < nsegs; i++) {
1045  strips[i].nverts_per_strip =
1046  strips[nstrips-1-i].nverts_per_strip = i;
1047  strips[i].nverts =
1048  strips[nstrips-1-i].nverts = i * 4;
1049  strips[i].nfaces =
1050  strips[nstrips-1-i].nfaces = (2 * i + 1)*4;
1051  }
1052  /* All strips have vertices and normals */
1053  for (i = 0; i < nstrips; i++) {
1054  strips[i].vp = (struct vertex **)bu_calloc(strips[i].nverts,
1055  sizeof(struct vertex *), "strip vertex[]");
1056  strips[i].norms = (vect_t *)bu_calloc(strips[i].nverts,
1057  sizeof(vect_t), "strip normals[]");
1058  }
1059  /* All strips have faces, except for the equator */
1060  for (i = 0; i < nstrips; i++) {
1061  if (strips[i].nfaces <= 0) continue;
1062  strips[i].fu = (struct faceuse **)bu_calloc(strips[i].nfaces,
1063  sizeof(struct faceuse *), "strip faceuse[]");
1064  }
1065 
1066  /* First, build the triangular mesh topology */
1067  /* Do the top. "toff" in i-1 is UP, towards +B */
1068  for (i = 1; i <= nsegs; i++) {
1069  faceno = 0;
1070  tlim = strips[i-1].nverts;
1071  blim = strips[i].nverts;
1072  for (stripno = 0; stripno < 4; stripno++) {
1073  toff = stripno * strips[i-1].nverts_per_strip;
1074  boff = stripno * strips[i].nverts_per_strip;
1075 
1076  /* Connect this quarter strip */
1077  for (j = 0; j < strips[i].nverts_per_strip; j++) {
1078 
1079  /* "Right-side-up" triangle */
1080  vertp[0] = &(strips[i].vp[j+boff]);
1081  vertp[1] = &(strips[i-1].vp[(j+toff)%tlim]);
1082  vertp[2] = &(strips[i].vp[(j+1+boff)%blim]);
1083  if ((strips[i-1].fu[faceno++] = nmg_cmface(state.s, vertp, 3)) == 0) {
1084  bu_log("rt_ell_tess() nmg_cmface failure\n");
1085  goto fail;
1086  }
1087  if (j+1 >= strips[i].nverts_per_strip) break;
1088 
1089  /* Follow with interior "Up-side-down" triangle */
1090  vertp[0] = &(strips[i].vp[(j+1+boff)%blim]);
1091  vertp[1] = &(strips[i-1].vp[(j+toff)%tlim]);
1092  vertp[2] = &(strips[i-1].vp[(j+1+toff)%tlim]);
1093  if ((strips[i-1].fu[faceno++] = nmg_cmface(state.s, vertp, 3)) == 0) {
1094  bu_log("rt_ell_tess() nmg_cmface failure\n");
1095  goto fail;
1096  }
1097  }
1098  }
1099  }
1100  /* Do the bottom. Everything is upside down. "toff" in i+1 is DOWN */
1101  for (i = nsegs; i < nstrips-1; i++) {
1102  faceno = 0;
1103  tlim = strips[i+1].nverts;
1104  blim = strips[i].nverts;
1105  for (stripno = 0; stripno < 4; stripno++) {
1106  toff = stripno * strips[i+1].nverts_per_strip;
1107  boff = stripno * strips[i].nverts_per_strip;
1108 
1109  /* Connect this quarter strip */
1110  for (j = 0; j < strips[i].nverts_per_strip; j++) {
1111 
1112  /* "Right-side-up" triangle */
1113  vertp[0] = &(strips[i].vp[j+boff]);
1114  vertp[1] = &(strips[i].vp[(j+1+boff)%blim]);
1115  vertp[2] = &(strips[i+1].vp[(j+toff)%tlim]);
1116  if ((strips[i+1].fu[faceno++] = nmg_cmface(state.s, vertp, 3)) == 0) {
1117  bu_log("rt_ell_tess() nmg_cmface failure\n");
1118  goto fail;
1119  }
1120  if (j+1 >= strips[i].nverts_per_strip) break;
1121 
1122  /* Follow with interior "Up-side-down" triangle */
1123  vertp[0] = &(strips[i].vp[(j+1+boff)%blim]);
1124  vertp[1] = &(strips[i+1].vp[(j+1+toff)%tlim]);
1125  vertp[2] = &(strips[i+1].vp[(j+toff)%tlim]);
1126  if ((strips[i+1].fu[faceno++] = nmg_cmface(state.s, vertp, 3)) == 0) {
1127  bu_log("rt_ell_tess() nmg_cmface failure\n");
1128  goto fail;
1129  }
1130  }
1131  }
1132  }
1133 
1134  /* Compute the geometry of each vertex. Start with the location
1135  * in the unit sphere, and project back. i=0 is "straight up"
1136  * along +B.
1137  */
1138  for (i = 0; i < nstrips; i++) {
1139  double alpha; /* decline down from B to A */
1140  double beta; /* angle around equator (azimuth) */
1141  fastf_t cos_alpha, sin_alpha;
1142  fastf_t cos_beta, sin_beta;
1143  point_t sphere_pt;
1144  point_t model_pt;
1145 
1146  alpha = (((double)i) / (nstrips-1));
1147  cos_alpha = cos(alpha*M_PI);
1148  sin_alpha = sin(alpha*M_PI);
1149  for (j = 0; j < strips[i].nverts; j++) {
1150 
1151  beta = ((double)j) / strips[i].nverts;
1152  cos_beta = cos(beta*M_2PI);
1153  sin_beta = sin(beta*M_2PI);
1154  VSET(sphere_pt,
1155  cos_beta * sin_alpha,
1156  cos_alpha,
1157  sin_beta * sin_alpha);
1158  /* Convert from ideal sphere coordinates */
1159  MAT4X3PNT(model_pt, state.invRinvS, sphere_pt);
1160  VADD2(model_pt, model_pt, state.eip->v);
1161  /* Associate vertex geometry */
1162  nmg_vertex_gv(strips[i].vp[j], model_pt);
1163 
1164  /* Convert sphere normal to ellipsoid normal */
1165  MAT4X3VEC(strips[i].norms[j], state.invRoS, sphere_pt);
1166  /* May not be unit length anymore */
1167  VUNITIZE(strips[i].norms[j]);
1168  }
1169  }
1170 
1171  /* Associate face geometry. Equator has no faces */
1172  for (i = 0; i < nstrips; i++) {
1173  for (j = 0; j < strips[i].nfaces; j++) {
1174  if (nmg_fu_planeeqn(strips[i].fu[j], tol) < 0)
1175  goto fail;
1176  }
1177  }
1178 
1179  /* Associate normals with vertexuses */
1180  for (i = 0; i < nstrips; i++) {
1181  for (j = 0; j < strips[i].nverts; j++) {
1182  struct faceuse *fu;
1183  struct vertexuse *vu;
1184  vect_t norm_opp;
1185 
1186  NMG_CK_VERTEX(strips[i].vp[j]);
1187  VREVERSE(norm_opp, strips[i].norms[j]);
1188 
1189  for (BU_LIST_FOR(vu, vertexuse, &strips[i].vp[j]->vu_hd)) {
1190  fu = nmg_find_fu_of_vu(vu);
1191  NMG_CK_FACEUSE(fu);
1192  /* get correct direction of normals depending on
1193  * faceuse orientation
1194  */
1195  if (fu->orientation == OT_SAME)
1196  nmg_vertexuse_nv(vu, strips[i].norms[j]);
1197  else if (fu->orientation == OT_OPPOSITE)
1198  nmg_vertexuse_nv(vu, norm_opp);
1199  }
1200  }
1201  }
1202 
1203  /* Compute "geometry" for region and shell */
1204  nmg_region_a(*r, tol);
1205 
1206  /* Release memory */
1207  /* All strips have vertices and normals */
1208  for (i = 0; i < nstrips; i++) {
1209  bu_free((char *)strips[i].vp, "strip vertex[]");
1210  bu_free((char *)strips[i].norms, "strip norms[]");
1211  }
1212  /* All strips have faces, except for equator */
1213  for (i = 0; i < nstrips; i++) {
1214  if (strips[i].fu == (struct faceuse **)0) continue;
1215  bu_free((char *)strips[i].fu, "strip faceuse[]");
1216  }
1217  bu_free((char *)strips, "strips[]");
1218  return 0;
1219 fail:
1220  /* Release memory */
1221  /* All strips have vertices and normals */
1222  for (i = 0; i < nstrips; i++) {
1223  bu_free((char *)strips[i].vp, "strip vertex[]");
1224  bu_free((char *)strips[i].norms, "strip norms[]");
1225  }
1226  /* All strips have faces, except for equator */
1227  for (i = 0; i < nstrips; i++) {
1228  if (strips[i].fu == (struct faceuse **)0) continue;
1229  bu_free((char *)strips[i].fu, "strip faceuse[]");
1230  }
1231  bu_free((char *)strips, "strips[]");
1232  return -1;
1233 }
1234 
1235 
1236 /**
1237  * Import an ellipsoid/sphere from the database format to the internal
1238  * structure. Apply modeling transformations as well.
1239  */
1240 int
1241 rt_ell_import4(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
1242 {
1243  struct rt_ell_internal *eip;
1244  union record *rp;
1245  fastf_t vec[3*4];
1246 
1247  if (dbip) RT_CK_DBI(dbip);
1248 
1249  BU_CK_EXTERNAL(ep);
1250  rp = (union record *)ep->ext_buf;
1251  /* Check record type */
1252  if (rp->u_id != ID_SOLID) {
1253  bu_log("rt_ell_import4: defective record\n");
1254  return -1;
1255  }
1256 
1257  RT_CK_DB_INTERNAL(ip);
1258  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1259  ip->idb_type = ID_ELL;
1260  ip->idb_meth = &OBJ[ID_ELL];
1261  BU_ALLOC(ip->idb_ptr, struct rt_ell_internal);
1262 
1263  eip = (struct rt_ell_internal *)ip->idb_ptr;
1264  eip->magic = RT_ELL_INTERNAL_MAGIC;
1265 
1266  /* Convert from database to internal format */
1267  flip_fastf_float(vec, rp->s.s_values, 4, dbip->dbi_version < 0 ? 1 : 0);
1268 
1269  /* Apply modeling transformations */
1270  if (mat == NULL) mat = bn_mat_identity;
1271  MAT4X3PNT(eip->v, mat, &vec[0*3]);
1272  MAT4X3VEC(eip->a, mat, &vec[1*3]);
1273  MAT4X3VEC(eip->b, mat, &vec[2*3]);
1274  MAT4X3VEC(eip->c, mat, &vec[3*3]);
1275 
1276  return 0; /* OK */
1277 }
1278 
1279 
1280 int
1281 rt_ell_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1282 {
1283  struct rt_ell_internal *tip;
1284  union record *rec;
1285 
1286  if (dbip) RT_CK_DBI(dbip);
1287 
1288  RT_CK_DB_INTERNAL(ip);
1289  if (ip->idb_type != ID_ELL && ip->idb_type != ID_SPH) return -1;
1290  tip = (struct rt_ell_internal *)ip->idb_ptr;
1291  RT_ELL_CK_MAGIC(tip);
1292 
1293  BU_CK_EXTERNAL(ep);
1294  ep->ext_nbytes = sizeof(union record);
1295  ep->ext_buf = (uint8_t *)bu_calloc(1, ep->ext_nbytes, "ell external");
1296  rec = (union record *)ep->ext_buf;
1297 
1298  rec->s.s_id = ID_SOLID;
1299  rec->s.s_type = GENELL;
1300 
1301  /* NOTE: This also converts to dbfloat_t */
1302  VSCALE(&rec->s.s_values[0], tip->v, local2mm);
1303  VSCALE(&rec->s.s_values[3], tip->a, local2mm);
1304  VSCALE(&rec->s.s_values[6], tip->b, local2mm);
1305  VSCALE(&rec->s.s_values[9], tip->c, local2mm);
1306 
1307  return 0;
1308 }
1309 
1310 
1311 /**
1312  * Import an ellipsoid/sphere from the database format to the internal
1313  * structure. Apply modeling transformations as well.
1314  */
1315 int
1316 rt_ell_import5(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
1317 {
1318  struct rt_ell_internal *eip;
1319 
1320  /* must be double for import and export */
1321  double vec[ELEMENTS_PER_VECT*4];
1322 
1323  if (dbip) RT_CK_DBI(dbip);
1324  RT_CK_DB_INTERNAL(ip);
1325  BU_CK_EXTERNAL(ep);
1326 
1327  BU_ASSERT_LONG(ep->ext_nbytes, ==, SIZEOF_NETWORK_DOUBLE * ELEMENTS_PER_VECT*4);
1328 
1329  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1330  ip->idb_type = ID_ELL;
1331  ip->idb_meth = &OBJ[ID_ELL];
1332  BU_ALLOC(ip->idb_ptr, struct rt_ell_internal);
1333 
1334  eip = (struct rt_ell_internal *)ip->idb_ptr;
1335  eip->magic = RT_ELL_INTERNAL_MAGIC;
1336 
1337  /* Convert from database (network) to internal (host) format */
1338  bu_cv_ntohd((unsigned char *)vec, ep->ext_buf, ELEMENTS_PER_VECT*4);
1339 
1340  /* Apply modeling transformations */
1341  if (mat == NULL) mat = bn_mat_identity;
1342  MAT4X3PNT(eip->v, mat, &vec[0*ELEMENTS_PER_VECT]);
1343  MAT4X3VEC(eip->a, mat, &vec[1*ELEMENTS_PER_VECT]);
1344  MAT4X3VEC(eip->b, mat, &vec[2*ELEMENTS_PER_VECT]);
1345  MAT4X3VEC(eip->c, mat, &vec[3*ELEMENTS_PER_VECT]);
1346 
1347  return 0; /* OK */
1348 }
1349 
1350 
1351 /**
1352  * The external format is:
1353  * V point
1354  * A vector
1355  * B vector
1356  * C vector
1357  */
1358 int
1359 rt_ell_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1360 {
1361  struct rt_ell_internal *eip;
1362 
1363  /* must be double for import and export */
1364  double vec[ELEMENTS_PER_VECT*4];
1365 
1366  if (dbip) RT_CK_DBI(dbip);
1367 
1368  RT_CK_DB_INTERNAL(ip);
1369  if (ip->idb_type != ID_ELL && ip->idb_type != ID_SPH) return -1;
1370  eip = (struct rt_ell_internal *)ip->idb_ptr;
1371  RT_ELL_CK_MAGIC(eip);
1372 
1373  BU_CK_EXTERNAL(ep);
1374  ep->ext_nbytes = SIZEOF_NETWORK_DOUBLE * ELEMENTS_PER_VECT*4;
1375  ep->ext_buf = (uint8_t *)bu_malloc(ep->ext_nbytes, "ell external");
1376 
1377  /* scale 'em into local buffer */
1378  VSCALE(&vec[0*ELEMENTS_PER_VECT], eip->v, local2mm);
1379  VSCALE(&vec[1*ELEMENTS_PER_VECT], eip->a, local2mm);
1380  VSCALE(&vec[2*ELEMENTS_PER_VECT], eip->b, local2mm);
1381  VSCALE(&vec[3*ELEMENTS_PER_VECT], eip->c, local2mm);
1382 
1383  /* Convert from internal (host) to database (network) format */
1384  bu_cv_htond(ep->ext_buf, (unsigned char *)vec, ELEMENTS_PER_VECT*4);
1385 
1386  return 0;
1387 }
1388 
1389 
1390 /**
1391  * Make human-readable formatted presentation of this solid. First
1392  * line describes type of solid. Additional lines are indented one
1393  * tab, and give parameter values.
1394  */
1395 int
1396 rt_ell_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
1397 {
1398  register struct rt_ell_internal *tip =
1399  (struct rt_ell_internal *)ip->idb_ptr;
1400  fastf_t mag_a, mag_b, mag_c;
1401  char buf[256];
1402  double angles[5];
1403  vect_t unitv;
1404 
1405  RT_ELL_CK_MAGIC(tip);
1406  bu_vls_strcat(str, "ellipsoid (ELL)\n");
1407 
1408  sprintf(buf, "\tV (%g, %g, %g)\n",
1409  INTCLAMP(tip->v[X] * mm2local),
1410  INTCLAMP(tip->v[Y] * mm2local),
1411  INTCLAMP(tip->v[Z] * mm2local));
1412  bu_vls_strcat(str, buf);
1413 
1414  mag_a = MAGNITUDE(tip->a);
1415  mag_b = MAGNITUDE(tip->b);
1416  mag_c = MAGNITUDE(tip->c);
1417 
1418  sprintf(buf, "\tA (%g, %g, %g) mag=%g\n",
1419  INTCLAMP(tip->a[X] * mm2local),
1420  INTCLAMP(tip->a[Y] * mm2local),
1421  INTCLAMP(tip->a[Z] * mm2local),
1422  INTCLAMP(mag_a * mm2local));
1423  bu_vls_strcat(str, buf);
1424 
1425  sprintf(buf, "\tB (%g, %g, %g) mag=%g\n",
1426  INTCLAMP(tip->b[X] * mm2local),
1427  INTCLAMP(tip->b[Y] * mm2local),
1428  INTCLAMP(tip->b[Z] * mm2local),
1429  INTCLAMP(mag_b * mm2local));
1430  bu_vls_strcat(str, buf);
1431 
1432  sprintf(buf, "\tC (%g, %g, %g) mag=%g\n",
1433  INTCLAMP(tip->c[X] * mm2local),
1434  INTCLAMP(tip->c[Y] * mm2local),
1435  INTCLAMP(tip->c[Z] * mm2local),
1436  INTCLAMP(mag_c * mm2local));
1437  bu_vls_strcat(str, buf);
1438 
1439  if (!verbose) return 0;
1440 
1441  VSCALE(unitv, tip->a, 1/mag_a);
1442  rt_find_fallback_angle(angles, unitv);
1443  rt_pr_fallback_angle(str, "\tA", angles);
1444 
1445  VSCALE(unitv, tip->b, 1/mag_b);
1446  rt_find_fallback_angle(angles, unitv);
1447  rt_pr_fallback_angle(str, "\tB", angles);
1448 
1449  VSCALE(unitv, tip->c, 1/mag_c);
1450  rt_find_fallback_angle(angles, unitv);
1451  rt_pr_fallback_angle(str, "\tC", angles);
1452 
1453  return 0;
1454 }
1455 
1456 
1457 /**
1458  * Free the storage associated with the rt_db_internal version of this
1459  * solid.
1460  */
1461 void
1463 {
1464  RT_CK_DB_INTERNAL(ip);
1465 
1466  bu_free(ip->idb_ptr, "ell ifree");
1467  ip->idb_ptr = ((void *)0);
1468 }
1469 
1470 
1471 /* The U parameter runs south to north. In order to orient loop CCW,
1472  * need to start with 0, 1-->0, 0 transition at the south pole.
1473  */
1474 static const fastf_t rt_ell_uvw[5*ELEMENTS_PER_VECT] = {
1475  0, 1, 0,
1476  0, 0, 0,
1477  1, 0, 0,
1478  1, 1, 0,
1479  0, 1, 0
1480 };
1481 
1482 
1483 int
1484 rt_ell_tnurb(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct bn_tol *tol)
1485 {
1486  mat_t R;
1487  mat_t S;
1488  mat_t invR;
1489  mat_t invS;
1490  mat_t invRinvS;
1491  mat_t invRoS;
1492  mat_t unit2model;
1493  mat_t xlate;
1494  vect_t Au, Bu, Cu; /* A, B, C with unit length */
1495  fastf_t Alen, Blen, Clen;
1496  fastf_t invAlen, invBlen, invClen;
1497  fastf_t magsq_a, magsq_b, magsq_c;
1498  fastf_t f;
1499  register int i;
1500  fastf_t radius;
1501  struct rt_ell_internal *eip;
1502  struct vertex *verts[8];
1503  struct vertex **vertp[4];
1504  struct faceuse *fu;
1505  struct shell *s;
1506  struct loopuse *lu;
1507  struct edgeuse *eu;
1508  point_t pole;
1509 
1510  RT_CK_DB_INTERNAL(ip);
1511  eip = (struct rt_ell_internal *)ip->idb_ptr;
1512  RT_ELL_CK_MAGIC(eip);
1513 
1514  /* Validate that |A| > 0, |B| > 0, |C| > 0 */
1515  magsq_a = MAGSQ(eip->a);
1516  magsq_b = MAGSQ(eip->b);
1517  magsq_c = MAGSQ(eip->c);
1518  if (magsq_a < tol->dist_sq || magsq_b < tol->dist_sq || magsq_c < tol->dist_sq) {
1519  bu_log("rt_ell_tnurb(): zero length A, B, or C vector\n");
1520  return -2; /* BAD */
1521  }
1522 
1523  /* Create unit length versions of A, B, C */
1524  invAlen = 1.0/(Alen = sqrt(magsq_a));
1525  VSCALE(Au, eip->a, invAlen);
1526  invBlen = 1.0/(Blen = sqrt(magsq_b));
1527  VSCALE(Bu, eip->b, invBlen);
1528  invClen = 1.0/(Clen = sqrt(magsq_c));
1529  VSCALE(Cu, eip->c, invClen);
1530 
1531  /* Validate that A.B == 0, B.C == 0, A.C == 0 (check dir only) */
1532  f = VDOT(Au, Bu);
1533  if (! NEAR_ZERO(f, tol->dist)) {
1534  bu_log("rt_ell_tnurb(): A not perpendicular to B, f=%f\n", f);
1535  return -3; /* BAD */
1536  }
1537  f = VDOT(Bu, Cu);
1538  if (! NEAR_ZERO(f, tol->dist)) {
1539  bu_log("rt_ell_tnurb(): B not perpendicular to C, f=%f\n", f);
1540  return -3; /* BAD */
1541  }
1542  f = VDOT(Au, Cu);
1543  if (! NEAR_ZERO(f, tol->dist)) {
1544  bu_log("rt_ell_tnurb(): A not perpendicular to C, f=%f\n", f);
1545  return -3; /* BAD */
1546  }
1547 
1548  {
1549  vect_t axb;
1550  VCROSS(axb, Au, Bu);
1551  f = VDOT(axb, Cu);
1552  if (f < 0) {
1553  VREVERSE(Cu, Cu);
1554  VREVERSE(eip->c, eip->c);
1555  }
1556  }
1557 
1558  /* Compute R and Rinv matrices */
1559  MAT_IDN(R);
1560  VMOVE(&R[0], Au);
1561  VMOVE(&R[4], Bu);
1562  VMOVE(&R[8], Cu);
1563  bn_mat_trn(invR, R); /* inv of rot mat is trn */
1564 
1565  /* Compute S and invS matrices */
1566  /* invS is just 1/diagonal elements */
1567  MAT_IDN(S);
1568  S[ 0] = invAlen;
1569  S[ 5] = invBlen;
1570  S[10] = invClen;
1571  bn_mat_inv(invS, S);
1572 
1573  /* invRinvS, for converting points from unit sphere to model */
1574  bn_mat_mul(invRinvS, invR, invS);
1575 
1576  /* invRoS, for converting normals from unit sphere to model */
1577  bn_mat_mul(invRoS, invR, S);
1578 
1579  /* Compute radius of bounding sphere */
1580  radius = Alen;
1581  if (Blen > radius)
1582  radius = Blen;
1583  if (Clen > radius)
1584  radius = Clen;
1585 
1586  MAT_IDN(xlate);
1587  MAT_DELTAS_VEC(xlate, eip->v);
1588  bn_mat_mul(unit2model, xlate, invRinvS);
1589 
1590  /*
1591  * --- Build Topology ---
1592  *
1593  * There is a vertex at either pole, and a single longitude line.
1594  * There is a single face, an snurb with singularities.
1595  *
1596  * vert[0] is the south pole, and is the first row of the ctl_points.
1597  * vert[1] is the north pole, and is the last row of the ctl_points.
1598  *
1599  * Somewhat surprisingly, the U parameter runs from south to north.
1600  */
1601  for (i = 0; i < 8; i++) verts[i] = (struct vertex *)0;
1602 
1603  *r = nmg_mrsv(m); /* Make region, empty shell, vertex */
1604  s = BU_LIST_FIRST(shell, &(*r)->s_hd);
1605 
1606  vertp[0] = &verts[0];
1607  vertp[1] = &verts[0];
1608  vertp[2] = &verts[1];
1609  vertp[3] = &verts[1];
1610 
1611  if ((fu = nmg_cmface(s, vertp, 4)) == 0) {
1612  bu_log("rt_ell_tnurb(): nmg_cmface() fail on face\n");
1613  return -1;
1614  }
1615 
1616  /* March around the fu's loop assigning uv parameter values */
1617  lu = BU_LIST_FIRST(loopuse, &fu->lu_hd);
1618  NMG_CK_LOOPUSE(lu);
1619  eu = BU_LIST_FIRST(edgeuse, &lu->down_hd);
1620  NMG_CK_EDGEUSE(eu);
1621 
1622  /* Loop always has Counter-Clockwise orientation (CCW) */
1623  for (i = 0; i < 4; i++) {
1624  nmg_vertexuse_a_cnurb(eu->vu_p, &rt_ell_uvw[i*ELEMENTS_PER_VECT]);
1625  nmg_vertexuse_a_cnurb(eu->eumate_p->vu_p, &rt_ell_uvw[(i+1)*ELEMENTS_PER_VECT]);
1626  eu = BU_LIST_NEXT(edgeuse, &eu->l);
1627  }
1628 
1629  /* Associate vertex geometry */
1630  VSUB2(pole, eip->v, eip->c); /* south pole */
1631  nmg_vertex_gv(verts[0], pole);
1632  VADD2(pole, eip->v, eip->c);
1633  nmg_vertex_gv(verts[1], pole); /* north pole */
1634 
1635  /* Build snurb, transformed into final position */
1636  nmg_sphere_face_snurb(fu, unit2model);
1637 
1638  /* Associate edge geometry (trimming curve) -- linear in param space */
1639  eu = BU_LIST_FIRST(edgeuse, &lu->down_hd);
1640  NMG_CK_EDGEUSE(eu);
1641  for (i = 0; i < 4; i++) {
1643  eu = BU_LIST_NEXT(edgeuse, &eu->l);
1644  }
1645 
1646  /* Compute "geometry" for region and shell */
1647  nmg_region_a(*r, tol);
1648 
1649  return 0;
1650 }
1651 
1652 
1653 /* u, v=(0, 0) is supposed to be the south pole, at Z=-1.0
1654  *
1655  * The V direction runs from the south to the north pole.
1656  */
1657 HIDDEN void
1658 nmg_sphere_face_snurb(struct faceuse *fu, const matp_t m)
1659 {
1660  struct face_g_snurb *fg;
1661  fastf_t root2_2;
1662  register fastf_t *op;
1663 
1664  NMG_CK_FACEUSE(fu);
1665  root2_2 = sqrt(2.0)*0.5;
1666 
1667  /* Let the library allocate all the storage */
1668  /* The V direction runs from south to north pole */
1669  nmg_face_g_snurb(fu,
1670  3, 3, /* u, v order */
1671  8, 12, /* Number of knots, u, v */
1672  NULL, NULL, /* initial u, v knot vectors */
1673  9, 5, /* n_rows, n_cols */
1674  RT_NURB_MAKE_PT_TYPE(4, RT_NURB_PT_XYZ, RT_NURB_PT_RATIONAL),
1675  NULL); /* initial mesh */
1676 
1677  fg = fu->f_p->g.snurb_p;
1678  NMG_CK_FACE_G_SNURB(fg);
1679 
1680  fg->v.knots[ 0] = 0;
1681  fg->v.knots[ 1] = 0;
1682  fg->v.knots[ 2] = 0;
1683  fg->v.knots[ 3] = 0.25;
1684  fg->v.knots[ 4] = 0.25;
1685  fg->v.knots[ 5] = 0.5;
1686  fg->v.knots[ 6] = 0.5;
1687  fg->v.knots[ 7] = 0.75;
1688  fg->v.knots[ 8] = 0.75;
1689  fg->v.knots[ 9] = 1;
1690  fg->v.knots[10] = 1;
1691  fg->v.knots[11] = 1;
1692 
1693  fg->u.knots[0] = 0;
1694  fg->u.knots[1] = 0;
1695  fg->u.knots[2] = 0;
1696  fg->u.knots[3] = 0.5;
1697  fg->u.knots[4] = 0.5;
1698  fg->u.knots[5] = 1;
1699  fg->u.knots[6] = 1;
1700  fg->u.knots[7] = 1;
1701 
1702  op = fg->ctl_points;
1703 
1704 /* Inspired by MAT4X4PNT */
1705 #define M(x, y, z, w) { \
1706  *op++ = m[ 0]*(x) + m[ 1]*(y) + m[ 2]*(z) + m[ 3]*(w);\
1707  *op++ = m[ 4]*(x) + m[ 5]*(y) + m[ 6]*(z) + m[ 7]*(w);\
1708  *op++ = m[ 8]*(x) + m[ 9]*(y) + m[10]*(z) + m[11]*(w);\
1709  *op++ = m[12]*(x) + m[13]*(y) + m[14]*(z) + m[15]*(w); }
1710 
1711 
1712  M(0 , 0 , -1.0 , 1.0);
1713  M(root2_2 , 0 , -root2_2 , root2_2);
1714  M(1.0 , 0 , 0 , 1.0);
1715  M(root2_2 , 0 , root2_2 , root2_2);
1716  M(0 , 0 , 1.0 , 1.0);
1717 
1718  M(0 , 0 , -root2_2 , root2_2);
1719  M(0.5 , -0.5 , -0.5 , 0.5);
1720  M(root2_2 , -root2_2 , 0 , root2_2);
1721  M(0.5 , -0.5 , 0.5 , 0.5);
1722  M(0 , 0 , root2_2 , root2_2);
1723 
1724  M(0 , 0 , -1.0 , 1.0);
1725  M(0 , -root2_2 , -root2_2 , root2_2);
1726  M(0 , -1.0 , 0 , 1.0);
1727  M(0 , -root2_2 , root2_2 , root2_2);
1728  M(0 , 0 , 1.0 , 1.0);
1729 
1730  M(0 , 0 , -root2_2 , root2_2);
1731  M(-0.5 , -0.5 , -0.5 , 0.5);
1732  M(-root2_2 , -root2_2 , 0 , root2_2);
1733  M(-0.5 , -0.5 , 0.5 , 0.5);
1734  M(0 , 0 , root2_2 , root2_2);
1735 
1736  M(0 , 0 , -1.0 , 1.0);
1737  M(-root2_2 , 0 , -root2_2 , root2_2);
1738  M(-1.0 , 0 , 0 , 1.0);
1739  M(-root2_2 , 0 , root2_2 , root2_2);
1740  M(0 , 0 , 1.0 , 1.0);
1741 
1742  M(0 , 0 , -root2_2 , root2_2);
1743  M(-0.5 , 0.5 , -0.5 , 0.5);
1744  M(-root2_2 , root2_2 , 0 , root2_2);
1745  M(-0.5 , 0.5 , 0.5 , 0.5);
1746  M(0 , 0 , root2_2 , root2_2);
1747 
1748  M(0 , 0 , -1.0 , 1.0);
1749  M(0 , root2_2 , -root2_2 , root2_2);
1750  M(0 , 1.0 , 0 , 1.0);
1751  M(0 , root2_2 , root2_2 , root2_2);
1752  M(0 , 0 , 1.0 , 1.0);
1753 
1754  M(0 , 0 , -root2_2 , root2_2);
1755  M(0.5 , 0.5 , -0.5 , 0.5);
1756  M(root2_2 , root2_2 , 0 , root2_2);
1757  M(0.5 , 0.5 , 0.5 , 0.5);
1758  M(0 , 0 , root2_2 , root2_2);
1759 
1760  M(0 , 0 , -1.0 , 1.0);
1761  M(root2_2 , 0 , -root2_2 , root2_2);
1762  M(1.0 , 0 , 0 , 1.0);
1763  M(root2_2 , 0 , root2_2 , root2_2);
1764  M(0 , 0 , 1.0 , 1.0);
1765 }
1766 
1767 
1768 /**
1769  * @brief Method for declaration of parameters so that it can be used
1770  * by Parametrics and Constraints library.
1771  *
1772  * allocates space for PC set , User is expect to free after use
1773  *
1774  * @return 0 on success
1775  * @return -1 on failure
1776  */
1777 int
1778 rt_ell_params(struct pc_pc_set *UNUSED(pcs), const struct rt_db_internal *UNUSED(ip))
1779 {
1780  return -1; /* FAIL */
1781 }
1782 
1783 
1784 /*
1785  * Used by EHY, EPA, HYP. See librt_private.h for details.
1786  */
1787 fastf_t
1789 {
1790  fastf_t dist, intr, m, theta0, theta1;
1791  point_t mpt, p0;
1792  vect_t norm_line, norm_ell;
1793 
1794  VSET(p0, a, 0., 0.);
1795  /* slope and intercept of segment */
1796  m = (p1[Y] - p0[Y]) / (p1[X] - p0[X]);
1797  intr = p0[Y] - m * p0[X];
1798  /* point on ellipse with max dist between ellipse and line */
1799  mpt[X] = a / sqrt(b*b / (m*m*a*a) + 1);
1800  mpt[Y] = b * sqrt(1 - mpt[X] * mpt[X] / (a*a));
1801  mpt[Z] = 0;
1802  /* max distance between that point and line */
1803  dist = fabs(m * mpt[X] - mpt[Y] + intr) / sqrt(m * m + 1);
1804  /* angles between normal of line and of ellipse at line endpoints */
1805  VSET(norm_line, m, -1., 0.);
1806  VSET(norm_ell, b * b * p0[X], a * a * p0[Y], 0.);
1807  VUNITIZE(norm_line);
1808  VUNITIZE(norm_ell);
1809  theta0 = fabs(acos(VDOT(norm_line, norm_ell)));
1810  VSET(norm_ell, b * b * p1[X], a * a * p1[Y], 0.);
1811  VUNITIZE(norm_ell);
1812  theta1 = fabs(acos(VDOT(norm_line, norm_ell)));
1813  /* split segment at widest point if not within error tolerances */
1814  if (dist > dtol || theta0 > ntol || theta1 > ntol) {
1815  /* split segment */
1816  return ell_angle(mpt, a, b, dtol, ntol);
1817  } else
1818  return (acos(VDOT(p0, p1)
1819  / (MAGNITUDE(p0) * MAGNITUDE(p1))));
1820 }
1821 
1822 
1823 /**
1824  * Computes volume of a ellipsoid.
1825  */
1826 void
1827 rt_ell_volume(fastf_t *volume, const struct rt_db_internal *ip)
1828 {
1829  fastf_t mag_a, mag_b, mag_c;
1830  struct rt_ell_internal *eip = (struct rt_ell_internal *)ip->idb_ptr;
1831  RT_ELL_CK_MAGIC(eip);
1832 
1833  mag_a = MAGNITUDE(eip->a);
1834  mag_b = MAGNITUDE(eip->b);
1835  mag_c = MAGNITUDE(eip->c);
1836  *volume = 4.0/3.0 * M_PI * mag_a * mag_b * mag_c;
1837 }
1838 
1839 
1840 /**
1841  * Computes centroid of an ellipsoid
1842  */
1843 void
1844 rt_ell_centroid(point_t *cent, const struct rt_db_internal *ip)
1845 {
1846  struct rt_ell_internal *eip = (struct rt_ell_internal *)ip->idb_ptr;
1847  RT_ELL_CK_MAGIC(eip);
1848  VMOVE(*cent,eip->v);
1849 }
1850 
1851 
1852 /**
1853  * Computes the surface area of an ellipsoid.
1854  * logs error if triaxial ellipsoid (cannot find surface area).
1855  */
1856 #define PROLATE 1
1857 #define OBLATE 2
1858 void
1859 rt_ell_surf_area(fastf_t *area, const struct rt_db_internal *ip)
1860 {
1861  fastf_t mag_a, mag_b, mag_c;
1862 #ifdef major /* Some systems have these defined as macros!!! */
1863 #undef major
1864 #endif
1865 #ifdef minor
1866 #undef minor
1867 #endif
1868  fastf_t major = 0, minor = 0;
1869  fastf_t major2, minor2;
1870  fastf_t ecc;
1871  int ell_type = 0;
1872 
1873  struct rt_ell_internal *eip = (struct rt_ell_internal *)ip->idb_ptr;
1874  RT_ELL_CK_MAGIC(eip);
1875 
1876  mag_a = MAGNITUDE(eip->a);
1877  mag_b = MAGNITUDE(eip->b);
1878  mag_c = MAGNITUDE(eip->c);
1879 
1880  if (EQUAL(mag_a, mag_b) && EQUAL(mag_b, mag_c)) {
1881  /* case: sphere */
1882  *area = 4.0 * M_PI * mag_a * mag_a;
1883  return;
1884  }
1885 
1886  if (EQUAL(mag_a, mag_b)) {
1887  if (mag_a > mag_c) {
1888  /* case: prolate spheroid */
1889  ell_type = PROLATE;
1890  major = mag_a;
1891  minor = mag_c;
1892  } else {
1893  /* case: oblate spheroid */
1894  ell_type = OBLATE;
1895  major = mag_c;
1896  minor = mag_a;
1897  }
1898  } else if (EQUAL(mag_a, mag_c)) {
1899  if (mag_a > mag_b) {
1900  /* case: prolate spheroid */
1901  ell_type = PROLATE;
1902  major = mag_a;
1903  minor = mag_b;
1904  } else {
1905  /* case: oblate spheroid */
1906  ell_type = OBLATE;
1907  major = mag_b;
1908  minor = mag_a;
1909  }
1910  } else if (EQUAL(mag_b, mag_c)) {
1911  if (mag_a > mag_c) {
1912  /* case: prolate spheroid */
1913  ell_type = PROLATE;
1914  major = mag_a;
1915  minor = mag_c;
1916  } else {
1917  /* case: oblate spheroid */
1918  ell_type = OBLATE;
1919  major = mag_c;
1920  minor = mag_a;
1921  }
1922  }
1923 
1924  major2 = major * major;
1925  minor2 = minor * minor;
1926  ecc = sqrt(1.0 - (minor2 / major2));
1927 
1928  switch (ell_type) {
1929  case PROLATE:
1930  *area = (M_2PI * minor2) + (M_2PI * major * minor / ecc) * asin(ecc);
1931  break;
1932  case OBLATE:
1933  *area = (M_2PI * major2) + (M_PI * minor2 / ecc) * log((1.0 + ecc) / (1.0 - ecc));
1934  break;
1935  default:
1936  bu_log("rt_ell_surf_area(): triaxial ellipsoid, cannot find surface area");
1937  }
1938 }
1939 
1940 
1941 /** @} */
1942 /*
1943  * Local Variables:
1944  * mode: C
1945  * tab-width: 8
1946  * indent-tabs-mode: t
1947  * c-file-style: "stroustrup"
1948  * End:
1949  * ex: shiftwidth=4 tabstop=8
1950  */
#define OBLATE
Definition: ell.c:1857
void rt_ell_ifree(struct rt_db_internal *ip)
Definition: ell.c:1462
#define BU_LIST_FOR(p, structure, hp)
Definition: list.h:365
vect_t ell_axis_vector_b
Definition: ell.c:647
Definition: raytrace.h:800
double ellipsoid_travel_axis_position
Definition: ell.c:658
int rt_ell_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
Definition: ell.c:389
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define BU_LIST_INSERT(old, new)
Definition: list.h:183
void rt_ell_centroid(point_t *cent, const struct rt_db_internal *ip)
Definition: ell.c:1844
#define SIZEOF_NETWORK_DOUBLE
Definition: cv.h:48
#define ELLOUT(n)
Definition: ell.c:604
struct hit seg_in
IN information.
Definition: raytrace.h:370
int(* ft_bbox)(struct rt_db_internal *, point_t *, point_t *, const struct bn_tol *)
Definition: raytrace.h:2196
struct vertex ** vp
Definition: ell.c:859
Definition: list.h:118
int nmg_fu_planeeqn(struct faceuse *fu, const struct bn_tol *tol)
Definition: nmg_mod.c:1311
vect_t ell_Cu
Definition: ell.c:159
#define RT_CK_APPLICATION(_p)
Definition: raytrace.h:1675
double dist
>= 0
Definition: tol.h:73
const struct bu_structparse rt_ell_parse[]
Definition: ell.c:49
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 points_per_section
Definition: ell.c:650
struct soltab * seg_stp
pointer back to soltab
Definition: raytrace.h:372
if lu s
Definition: nmg_mod.c:3860
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
Definition: raytrace.h:215
mat_t invRoS
Definition: ell.c:851
#define M_PI
Definition: fft.h:35
double ellipsoid_travel_axis_mag
Definition: ell.c:657
Definition: pc.h:108
Definition: raytrace.h:368
vect_t translation
Definition: ell.c:656
#define BU_ASSERT_LONG(_lhs, _relation, _rhs)
Definition: defines.h:240
Definition: raytrace.h:248
int rt_sph_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: sph.c:300
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
fastf_t ell_angle(fastf_t *p1, fastf_t a, fastf_t b, fastf_t dtol, fastf_t ntol)
Definition: ell.c:1788
int rt_ell_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: ell.c:1281
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)
#define M(x, y, z, w)
void flip_fastf_float(fastf_t *ff, const dbfloat_t *fp, int n, int flip)
Definition: db_flip.c:74
int rt_ell_params(struct pc_pc_set *pcs, const struct rt_db_internal *ip)
Method for declaration of parameters so that it can be used by Parametrics and Constraints library...
Definition: ell.c:1778
struct resource * a_resource
dynamic memory resources
Definition: raytrace.h:1591
int rt_ell_import5(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
Definition: ell.c:1316
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)
void rt_ell_vshot(struct soltab **stp, struct xray **rp, struct seg *segp, int n, struct application *ap)
Definition: ell.c:437
#define HIDDEN
Definition: common.h:86
NMG_CK_LOOPUSE(lu)
int rt_ell_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: ell.c:793
int rt_ell_import4(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
Definition: ell.c:1241
struct bu_list l
Definition: raytrace.h:369
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
struct shell * s
Definition: ell.c:848
if(share_geom)
Definition: nmg_mod.c:3829
int idb_major_type
Definition: raytrace.h:192
void rt_ell_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
Definition: ell.c:496
Definition: color.c:49
int rt_ell_tnurb(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct bn_tol *tol)
Definition: ell.c:1484
int nverts
Definition: ell.c:858
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
int nverts_per_strip
Definition: ell.c:857
int nfaces
Definition: ell.c:861
#define RT_CK_DB_INTERNAL(_p)
Definition: raytrace.h:207
#define RT_ELL_INTERNAL_MAGIC
Definition: magic.h:91
#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
vect_t ell_invsq
Definition: ell.c:160
#define BN_VLIST_LINE_MOVE
Definition: vlist.h:82
fastf_t a_diverge
slope of beam divergence/mm
Definition: raytrace.h:1600
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
void bn_mat_inv(mat_t output, const mat_t input)
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()
void rt_ell_print(register const struct soltab *stp)
Definition: ell.c:368
#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)
struct hit seg_out
OUT information.
Definition: raytrace.h:371
#define BN_VLIST_LINE_DRAW
Definition: vlist.h:83
vect_t ell_V
Definition: ell.c:156
void bn_mat_trn(mat_t om, const mat_t im)
#define UNUSED(parameter)
Definition: common.h:239
#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
#define ID_SPH
Sphere.
Definition: raytrace.h:468
mat_t invRinvS
Definition: ell.c:850
void rt_ell_16pts(fastf_t *ov, fastf_t *V, fastf_t *A, fastf_t *B)
Definition: ell.c:606
int num_cross_sections
Definition: ell.c:649
double bn_atan2(double x, double y)
Definition: mat.c:104
vect_t ell_travel_vector
Definition: ell.c:648
void rt_ell_surf_area(fastf_t *area, const struct rt_db_internal *ip)
Definition: ell.c:1859
#define bu_offsetofarray(_t, _a, _d, _i)
Definition: parse.h:65
void rt_ell_uv(struct application *ap, struct soltab *stp, register struct hit *hitp, register struct uvcoord *uvp)
Definition: ell.c:558
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
ustring alpha
vect_t * norms
Definition: ell.c:860
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
void nmg_edge_g_cnurb_plinear(struct edgeuse *eu)
Definition: nmg_mk.c:2019
int rt_ell_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: ell.c:257
vect_t ell_Au
Definition: ell.c:157
fastf_t primitive_curve_count(struct rt_db_internal *ip, const struct rt_view_info *info)
#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
int rt_ell_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
Definition: ell.c:1396
mat_t ell_SoR
Definition: ell.c:161
void * idb_ptr
Definition: raytrace.h:195
void rt_ell_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
Definition: ell.c:518
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
vect_t ell_Bu
Definition: ell.c:158
int rt_ell_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: ell.c:1359
void bn_vec_ortho(vect_t out, const vect_t in)
#define R
Definition: msr.c:55
fastf_t theta_tol
Definition: ell.c:852
int rt_ell_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
Definition: ell.c:896
void * st_specific
-> ID-specific (private) struct
Definition: raytrace.h:435
fastf_t uv_du
delta in u
Definition: raytrace.h:343
vect_t ell_axis_vector_a
Definition: ell.c:646
vect_t ell_center
Definition: ell.c:645
fastf_t crv_c1
curvature in principle dir
Definition: raytrace.h:308
#define A
Definition: msr.c:51
void rt_ell_free(register struct soltab *stp)
Definition: ell.c:592
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
HIDDEN void nmg_sphere_face_snurb(struct faceuse *fu, const matp_t m)
Definition: ell.c:1658
void rt_ell_volume(fastf_t *volume, const struct rt_db_internal *ip)
Definition: ell.c:1827
mat_t ell_invRSSR
Definition: ell.c:162
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
fastf_t a_rbeam
initial beam radius (mm)
Definition: raytrace.h:1599
#define BU_CK_LIST_HEAD(_p)
Definition: list.h:142
int rt_ell_adaptive_plot(struct rt_db_internal *ip, const struct rt_view_info *info)
Definition: ell.c:750
#define BU_CK_EXTERNAL(_p)
Definition: parse.h:224
int rt_ell_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *tol)
Definition: ell.c:172
struct rt_ell_internal * eip
Definition: ell.c:849
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
#define RT_ELL_SEG_MISS(SEG)
Definition: ell.c:432
const struct rt_functab * st_meth
pointer to per-solid methods
Definition: raytrace.h:428
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
Definition: ell.c:847
#define M_SQRT1_2
Definition: fft.h:38
double fastf_t
Definition: defines.h:300
#define VPRINT(a, b)
Definition: raytrace.h:1881
struct faceuse ** fu
Definition: ell.c:862
#define ID_ELL
Ellipsoid.
Definition: raytrace.h:461
#define BU_LIST_NEXT(structure, hp)
Definition: list.h:316
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 BN_VLIST_POINT_DRAW
Draw a single point.
Definition: vlist.h:94
struct bu_list * vhead
Definition: ell.c:644
#define PROLATE
Definition: ell.c:1856
void nmg_region_a(struct nmgregion *r, const struct bn_tol *tol)
Definition: nmg_mk.c:2557