BRL-CAD
superell.c
Go to the documentation of this file.
1 /* S U P E R 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/superell/superell.c
23  *
24  * Intersect a ray with a Superquadratic Ellipsoid.
25  *
26  * NOTICE: this primitive is incomplete and should be considered
27  * experimental. This primitive will exhibit several
28  * instabilities in the existing root solver method.
29  *
30  */
31 /** @} */
32 
33 #include "common.h"
34 
35 #include <stddef.h>
36 #include <string.h>
37 #include <math.h>
38 #include "bio.h"
39 
40 #include "bu/cv.h"
41 #include "vmath.h"
42 #include "db.h"
43 #include "nmg.h"
44 #include "rtgeom.h"
45 #include "raytrace.h"
46 
47 #include "../../librt_private.h"
48 
49 
51  { "%f", 3, "V", bu_offsetofarray(struct rt_superell_internal, v, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
52  { "%f", 3, "A", bu_offsetofarray(struct rt_superell_internal, a, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
53  { "%f", 3, "B", bu_offsetofarray(struct rt_superell_internal, b, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
54  { "%f", 3, "C", bu_offsetofarray(struct rt_superell_internal, c, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
55  { "%g", 1, "n", bu_offsetof(struct rt_superell_internal, n), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
56  { "%g", 1, "e", bu_offsetof(struct rt_superell_internal, e), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
57  { {'\0', '\0', '\0', '\0'}, 0, (char *)NULL, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL }
58 };
59 
60 
61 /*
62  * Algorithm:
63  *
64  * Given V, A, B, and C, there is a set of points on this superellipsoid
65  *
66  * { (x, y, z) | (x, y, z) is on superellipsoid defined by V, A, B, C }
67  *
68  * Through a series of Affine Transformations, this set will be
69  * transformed into a set of points on a unit sphere at the origin
70  *
71  * { (x', y', z') | (x', y', z') is on Sphere at origin }
72  *
73  * The transformation from X to X' is accomplished by:
74  *
75  * X' = S(R(X - V))
76  *
77  * where R(X) = (A/(|A|))
78  * (B/(|B|)) . X
79  * (C/(|C|))
80  *
81  * and S(X) = (1/|A| 0 0)
82  * (0 1/|B| 0) . X
83  * (0 0 1/|C|)
84  *
85  * To find the intersection of a line with the superellipsoid,
86  * consider the parametric line L:
87  *
88  * L : { P(n) | P + t(n) . D }
89  *
90  * Call W the actual point of intersection between L and the
91  * superellipsoid. Let W' be the point of intersection between L' and
92  * the unit sphere.
93  *
94  * L' : { P'(n) | P' + t(n) . D' }
95  *
96  * W = invR(invS(W')) + V
97  *
98  * Where W' = k D' + P'.
99  *
100  * Let dp = D' dot P'
101  * Let dd = D' dot D'
102  * Let pp = P' dot P'
103  *
104  * and k = [ -dp +/- sqrt(dp*dp - dd * (pp - 1)) ] / dd
105  * which is constant.
106  *
107  * Now, D' = S(R(D))
108  * and P' = S(R(P - V))
109  *
110  * Substituting,
111  *
112  * W = V + invR(invS[ k *(S(R(D))) + S(R(P - V)) ])
113  * = V + invR((k * R(D)) + R(P - V))
114  * = V + k * D + P - V
115  * = k * D + P
116  *
117  * Note that ``k'' is constant, and is the same in the formulations
118  * for both W and W'.
119  *
120  * NORMALS. Given the point W on the superellipsoid, what is the
121  * vector normal to the tangent plane at that point?
122  *
123  * Map W onto the unit sphere, i.e.: W' = S(R(W - V)).
124  *
125  * Plane on unit sphere at W' has a normal vector of the same value(!).
126  * N' = W'
127  *
128  * The plane transforms back to the tangent plane at W, and this new
129  * plane (on the superellipsoid) has a normal vector of N, viz:
130  *
131  * N = inverse[ transpose(inverse[ S o R ]) ] (N')
132  *
133  * because if H is perpendicular to plane Q, and matrix M maps from Q
134  * to Q', then inverse[ transpose(M) ] (H) is perpendicular to Q'.
135  * Here, H and Q are in "prime space" with the unit sphere. [Somehow,
136  * the notation here is backwards]. So, the mapping matrix M =
137  * inverse(S o R), because S o R maps from normal space to the unit
138  * sphere.
139  *
140  * N = inverse[ transpose(inverse[ S o R ]) ] (N')
141  * = inverse[ transpose(invR o invS) ] (N')
142  * = inverse[ transpose(invS) o transpose(invR) ] (N')
143  * = inverse[ inverse(S) o R ] (N')
144  * = invR o S (N')
145  *
146  * = invR o S (W')
147  * = invR(S(S(R(W - V))))
148  *
149  * because inverse(R) = transpose(R), so R = transpose(invR),
150  * and S = transpose(S).
151  *
152  * Note that the normal vector N produced above will not have unit
153  * length.
154  */
155 
157  vect_t superell_V; /* Vector to center of superellipsoid */
158  vect_t superell_Au; /* unit-length A vector */
159  vect_t superell_Bu;
160  vect_t superell_Cu;
161  double superell_n; /* north-south curvature power */
162  double superell_e; /* east-west curvature power */
163  double superell_invmsAu; /* 1.0 / |Au|^2 */
164  double superell_invmsBu; /* 1.0 / |Bu|^2 */
165  double superell_invmsCu; /* 1.0 / |Cu|^2 */
167  mat_t superell_SoR; /* matrix for local coordinate system, Scale(Rotate(V))*/
168  mat_t superell_invRSSR; /* invR(Scale(Scale(Rot(V)))) */
169  mat_t superell_invR; /* transposed rotation matrix */
170 };
171 #define SUPERELL_NULL ((struct superell_specific *)0)
172 
173 /**
174  * Calculate a bounding RPP for a superell
175  */
176 int
177 rt_superell_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol)) {
178 
179  struct rt_superell_internal *eip;
180  fastf_t magsq_a, magsq_b, magsq_c;
181  vect_t Au, Bu, Cu;
182  mat_t R;
183  vect_t w1, w2, P; /* used for bounding RPP */
184  fastf_t f;
185 
186  eip = (struct rt_superell_internal *)ip->idb_ptr;
187  RT_SUPERELL_CK_MAGIC(eip);
188 
189  magsq_a = MAGSQ(eip->a);
190  magsq_b = MAGSQ(eip->b);
191  magsq_c = MAGSQ(eip->c);
192 
193 
194  /* Create unit length versions of A, B, C */
195  f = 1.0/sqrt(magsq_a);
196  VSCALE(Au, eip->a, f);
197  f = 1.0/sqrt(magsq_b);
198  VSCALE(Bu, eip->b, f);
199  f = 1.0/sqrt(magsq_c);
200  VSCALE(Cu, eip->c, f);
201 
202  MAT_IDN(R);
203  VMOVE(&R[0], Au);
204  VMOVE(&R[4], Bu);
205  VMOVE(&R[8], Cu);
206 
207  /* Compute bounding RPP */
208  VSET(w1, magsq_a, magsq_b, magsq_c);
209 
210  /* X */
211  VSET(P, 1.0, 0, 0); /* bounding plane normal */
212  MAT3X3VEC(w2, R, P); /* map plane to local coord syst */
213  VELMUL(w2, w2, w2); /* square each term */
214  f = VDOT(w1, w2);
215  f = sqrt(f);
216  (*min)[X] = eip->v[X] - f; /* V.P +/- f */
217  (*max)[X] = eip->v[X] + f;
218 
219  /* Y */
220  VSET(P, 0, 1.0, 0); /* bounding plane normal */
221  MAT3X3VEC(w2, R, P); /* map plane to local coord syst */
222  VELMUL(w2, w2, w2); /* square each term */
223  f = VDOT(w1, w2);
224  f = sqrt(f);
225  (*min)[Y] = eip->v[Y] - f; /* V.P +/- f */
226  (*max)[Y] = eip->v[Y] + f;
227 
228  /* Z */
229  VSET(P, 0, 0, 1.0); /* bounding plane normal */
230  MAT3X3VEC(w2, R, P); /* map plane to local coord syst */
231  VELMUL(w2, w2, w2); /* square each term */
232  f = VDOT(w1, w2);
233  f = sqrt(f);
234  (*min)[Z] = eip->v[Z] - f; /* V.P +/- f */
235  (*max)[Z] = eip->v[Z] + f;
236 
237  return 0;
238 }
239 
240 
241 /**
242  * Given a pointer to a GED database record, and a transformation
243  * matrix, determine if this is a valid superellipsoid, and if so,
244  * precompute various terms of the formula.
245  *
246  * Returns -
247  * 0 SUPERELL is OK
248  * !0 Error in description
249  *
250  * Implicit return - A struct superell_specific is created, and its
251  * address is stored in stp->st_specific for use by rt_superell_shot()
252  */
253 int
254 rt_superell_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
255 {
256 
257  struct superell_specific *superell;
258  struct rt_superell_internal *eip;
259  fastf_t magsq_a, magsq_b, magsq_c;
260  mat_t R, TEMP;
261  vect_t Au, Bu, Cu; /* A, B, C with unit length */
262  fastf_t f;
263 
264  eip = (struct rt_superell_internal *)ip->idb_ptr;
265  RT_SUPERELL_CK_MAGIC(eip);
266 
267  /* Validate that |A| > 0, |B| > 0, |C| > 0 */
268  magsq_a = MAGSQ(eip->a);
269  magsq_b = MAGSQ(eip->b);
270  magsq_c = MAGSQ(eip->c);
271 
272  if (magsq_a < rtip->rti_tol.dist_sq || magsq_b < rtip->rti_tol.dist_sq || magsq_c < rtip->rti_tol.dist_sq) {
273  bu_log("rt_superell_prep(): superell(%s) near-zero length A(%g), B(%g), or C(%g) vector\n",
274  stp->st_name, magsq_a, magsq_b, magsq_c);
275  return 1; /* BAD */
276  }
277  if (eip->n < rtip->rti_tol.dist || eip->e < rtip->rti_tol.dist) {
278  bu_log("rt_superell_prep(): superell(%s) near-zero length <n, e> curvature (%g, %g) causes problems\n",
279  stp->st_name, eip->n, eip->e);
280  /* BAD */
281  }
282  if (eip->n > 10000.0 || eip->e > 10000.0) {
283  bu_log("rt_superell_prep(): superell(%s) very large <n, e> curvature (%g, %g) causes problems\n",
284  stp->st_name, eip->n, eip->e);
285  /* BAD */
286  }
287 
288  /* Create unit length versions of A, B, C */
289  f = 1.0/sqrt(magsq_a);
290  VSCALE(Au, eip->a, f);
291  f = 1.0/sqrt(magsq_b);
292  VSCALE(Bu, eip->b, f);
293  f = 1.0/sqrt(magsq_c);
294  VSCALE(Cu, eip->c, f);
295 
296  /* Validate that A.B == 0, B.C == 0, A.C == 0 (check dir only) */
297  f = VDOT(Au, Bu);
298  if (! NEAR_ZERO(f, rtip->rti_tol.dist)) {
299  bu_log("rt_superell_prep(): superell(%s) A not perpendicular to B, f=%f\n", stp->st_name, f);
300  return 1; /* BAD */
301  }
302  f = VDOT(Bu, Cu);
303  if (! NEAR_ZERO(f, rtip->rti_tol.dist)) {
304  bu_log("rt_superell_prep(): superell(%s) B not perpendicular to C, f=%f\n", stp->st_name, f);
305  return 1; /* BAD */
306  }
307  f = VDOT(Au, Cu);
308  if (! NEAR_ZERO(f, rtip->rti_tol.dist)) {
309  bu_log("rt_superell_prep(): superell(%s) A not perpendicular to C, f=%f\n", stp->st_name, f);
310  return 1; /* BAD */
311  }
312 
313  /* Solid is OK, compute constant terms now */
314 
315  BU_GET(superell, struct superell_specific);
316  stp->st_specific = (void *)superell;
317 
318  superell->superell_n = eip->n;
319  superell->superell_e = eip->e;
320 
321  VMOVE(superell->superell_V, eip->v);
322 
323  VSET(superell->superell_invsq, 1.0/magsq_a, 1.0/magsq_b, 1.0/magsq_c);
324  VMOVE(superell->superell_Au, Au);
325  VMOVE(superell->superell_Bu, Bu);
326  VMOVE(superell->superell_Cu, Cu);
327 
328  /* compute the inverse magnitude square for equations during shot */
329  superell->superell_invmsAu = 1.0 / magsq_a;
330  superell->superell_invmsBu = 1.0 / magsq_b;
331  superell->superell_invmsCu = 1.0 / magsq_c;
332 
333  /* compute the rotation matrix */
334  MAT_IDN(R);
335  VMOVE(&R[0], Au);
336  VMOVE(&R[4], Bu);
337  VMOVE(&R[8], Cu);
338  bn_mat_trn(superell->superell_invR, R);
339 
340  /* computer invRSSR */
341  MAT_IDN(superell->superell_invRSSR);
342  MAT_IDN(TEMP);
343  TEMP[0] = superell->superell_invsq[0];
344  TEMP[5] = superell->superell_invsq[1];
345  TEMP[10] = superell->superell_invsq[2];
346  bn_mat_mul(TEMP, TEMP, R);
347  bn_mat_mul(superell->superell_invRSSR, superell->superell_invR, TEMP);
348 
349  /* compute Scale(Rotate(vect)) */
350  MAT_IDN(superell->superell_SoR);
351  VSCALE(&superell->superell_SoR[0], eip->a, superell->superell_invsq[0]);
352  VSCALE(&superell->superell_SoR[4], eip->b, superell->superell_invsq[1]);
353  VSCALE(&superell->superell_SoR[8], eip->c, superell->superell_invsq[2]);
354 
355  /* Compute bounding sphere */
356  VMOVE(stp->st_center, eip->v);
357  f = magsq_a;
358  if (magsq_b > f)
359  f = magsq_b;
360  if (magsq_c > f)
361  f = magsq_c;
362  stp->st_aradius = stp->st_bradius = sqrt(f);
363 
364  /* Compute bounding RPP */
365  if (rt_superell_bbox(ip, &(stp->st_min), &(stp->st_max), &rtip->rti_tol)) return 1;
366 
367  return 0; /* OK */
368 }
369 
370 
371 void
372 rt_superell_print(const struct soltab *stp)
373 {
374  struct superell_specific *superell =
375  (struct superell_specific *)stp->st_specific;
376 
377  VPRINT("V", superell->superell_V);
378 }
379 
380 
381 /* Equation of a superellipsoid:
382  *
383  * f(x) = [ (x^(2/e2) + y^(2/e2))^(e2/e1) + z^(2/e1) ]^(e1/2) - 1
384  */
385 
386 /**
387  * Intersect a ray with an superellipsoid, where all constant terms
388  * have been precomputed by rt_superell_prep(). If an intersection
389  * occurs, a struct seg will be acquired and filled in.
390  *
391  * Returns -
392  * 0 MISS
393  * >0 HIT
394  */
395 int
396 rt_superell_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
397 {
398  static int counter=10;
399 
400  struct superell_specific *superell = (struct superell_specific *)stp->st_specific;
401  bn_poly_t equation; /* equation of superell to be solved */
402  vect_t translated; /* translated shot vector */
403  vect_t newShotPoint; /* P' */
404  vect_t newShotDir; /* D' */
405  vect_t normalizedShotPoint; /* P' with normalized dist from superell */
406  bn_complex_t complexRoot[4]; /* roots returned from poly solver */
407  double realRoot[4]; /* real ray distance values */
408  int i, j;
409  struct seg *segp;
410 
411  /* translate ray point */
412  /* VSUB2(translated, rp->r_pt, superell->superell_V); */
413  (translated)[X] = (rp->r_pt)[X] - (superell->superell_V)[X];
414  (translated)[Y] = (rp->r_pt)[Y] - (superell->superell_V)[Y];
415  (translated)[Z] = (rp->r_pt)[Z] - (superell->superell_V)[Z];
416 
417  /* scale and rotate point to get P' */
418 
419  /* MAT4X3VEC(newShotPoint, superell->superell_SoR, translated); */
420  newShotPoint[X] = (superell->superell_SoR[0]*translated[X] + superell->superell_SoR[1]*translated[Y] + superell->superell_SoR[ 2]*translated[Z]) * 1.0/(superell->superell_SoR[15]);
421  newShotPoint[Y] = (superell->superell_SoR[4]*translated[X] + superell->superell_SoR[5]*translated[Y] + superell->superell_SoR[ 6]*translated[Z]) * 1.0/(superell->superell_SoR[15]);
422  newShotPoint[Z] = (superell->superell_SoR[8]*translated[X] + superell->superell_SoR[9]*translated[Y] + superell->superell_SoR[10]*translated[Z]) * 1.0/(superell->superell_SoR[15]);
423 
424  /* translate ray direction vector */
425  MAT4X3VEC(newShotDir, superell->superell_SoR, rp->r_dir);
426  VUNITIZE(newShotDir);
427 
428  /* normalize distance from the superell. substitutes a corrected ray
429  * point, which contains a translation along the ray direction to the
430  * closest approach to vertex of the superell. Translating the ray
431  * along the direction of the ray to the closest point near the
432  * primitive's center vertex. New ray origin is hence, normalized.
433  */
434  VSCALE(normalizedShotPoint, newShotDir,
435  VDOT(newShotPoint, newShotDir));
436  VSUB2(normalizedShotPoint, newShotPoint, normalizedShotPoint);
437 
438  /* Now generate the polynomial equation for passing to the root finder */
439 
440  equation.dgr = 2;
441 
442  /* (x^2 / A) + (y^2 / B) + (z^2 / C) - 1 */
443  equation.cf[0] = newShotPoint[X] * newShotPoint[X] * superell->superell_invmsAu + newShotPoint[Y] * newShotPoint[Y] * superell->superell_invmsBu + newShotPoint[Z] * newShotPoint[Z] * superell->superell_invmsCu - 1;
444  /* (2xX / A) + (2yY / B) + (2zZ / C) */
445  equation.cf[1] = 2 * newShotDir[X] * newShotPoint[X] * superell->superell_invmsAu + 2 * newShotDir[Y] * newShotPoint[Y] * superell->superell_invmsBu + 2 * newShotDir[Z] * newShotPoint[Z] * superell->superell_invmsCu;
446  /* (X^2 / A) + (Y^2 / B) + (Z^2 / C) */
447  equation.cf[2] = newShotDir[X] * newShotDir[X] * superell->superell_invmsAu + newShotDir[Y] * newShotDir[Y] * superell->superell_invmsBu + newShotDir[Z] * newShotDir[Z] * superell->superell_invmsCu;
448 
449  if ((i = rt_poly_roots(&equation, complexRoot, stp->st_dp->d_namep)) != 2) {
450  if (i > 0) {
451  bu_log("rt_superell_shot(): poly roots %d != 2\n", i);
452  bn_pr_roots(stp->st_name, complexRoot, i);
453  } else if (i < 0) {
454  static int reported=0;
455  bu_log("rt_superell_shot(): The root solver failed to converge on a solution for %s\n", stp->st_dp->d_namep);
456  if (!reported) {
457  VPRINT("while shooting from:\t", rp->r_pt);
458  VPRINT("while shooting at:\t", rp->r_dir);
459  bu_log("rt_superell_shot(): Additional superellipsoid convergence failure details will be suppressed.\n");
460  reported=1;
461  }
462  }
463  return 0; /* MISS */
464  }
465 
466  /* XXX BEGIN CUT */
467  /* Only real roots indicate an intersection in real space.
468  *
469  * Look at each root returned; if the imaginary part is zero
470  * or sufficiently close, then use the real part as one value
471  * of 't' for the intersections
472  */
473  for (j=0, i=0; j < 2; j++) {
474  if (NEAR_ZERO(complexRoot[j].im, 0.001))
475  realRoot[i++] = complexRoot[j].re;
476  }
477 
478  /* reverse above translation by adding distance to all 'k' values. */
479  /* for (j = 0; j < i; ++j)
480  realRoot[j] -= VDOT(newShotPoint, newShotDir);
481  */
482 
483  /* Here, 'i' is number of points found */
484  switch (i) {
485  case 0:
486  return 0; /* No hit */
487 
488  default:
489  bu_log("rt_superell_shot(): reduced 4 to %d roots\n", i);
490  bn_pr_roots(stp->st_name, complexRoot, 4);
491  return 0; /* No hit */
492 
493  case 2:
494  {
495  /* Sort most distant to least distant. */
496  fastf_t u;
497  if ((u=realRoot[0]) < realRoot[1]) {
498  /* bubble larger towards [0] */
499  realRoot[0] = realRoot[1];
500  realRoot[1] = u;
501  }
502  }
503  break;
504  case 4:
505  {
506  short n;
507  short lim;
508 
509  /* Inline rt_pt_sort(). Sorts realRoot[] into descending order. */
510  for (lim = i-1; lim > 0; lim--) {
511  for (n = 0; n < lim; n++) {
512  fastf_t u;
513  if ((u=realRoot[n]) < realRoot[n+1]) {
514  /* bubble larger towards [0] */
515  realRoot[n] = realRoot[n+1];
516  realRoot[n+1] = u;
517  }
518  }
519  }
520  }
521  break;
522  }
523 
524  if (counter > 0) {
525  bu_log("rt_superell_shot(): realroot in %f out %f\n", realRoot[1], realRoot[0]);
526  counter--;
527  }
528 
529 
530  /* Now, t[0] > t[npts-1] */
531  /* realRoot[1] is entry point, and realRoot[0] is farthest exit point */
532  RT_GET_SEG(segp, ap->a_resource);
533  segp->seg_stp = stp;
534  segp->seg_in.hit_dist = realRoot[1];
535  segp->seg_out.hit_dist = realRoot[0];
536  /* segp->seg_in.hit_surfno = segp->seg_out.hit_surfno = 0; */
537  /* Set aside vector for rt_superell_norm() later */
538  /* VJOIN1(segp->seg_in.hit_vpriv, newShotPoint, realRoot[1], newShotDir); */
539  /* VJOIN1(segp->seg_out.hit_vpriv, newShotPoint, realRoot[0], newShotDir); */
540  BU_LIST_INSERT(&(seghead->l), &(segp->l));
541 
542  if (i == 2) {
543  return 2; /* HIT */
544  }
545 
546  /* 4 points */
547  /* realRoot[3] is entry point, and realRoot[2] is exit point */
548  RT_GET_SEG(segp, ap->a_resource);
549  segp->seg_stp = stp;
550  segp->seg_in.hit_dist = realRoot[3]*superell->superell_e;
551  segp->seg_out.hit_dist = realRoot[2]*superell->superell_e;
552  segp->seg_in.hit_surfno = segp->seg_out.hit_surfno = 1;
553  VJOIN1(segp->seg_in.hit_vpriv, newShotPoint, realRoot[3], newShotDir);
554  VJOIN1(segp->seg_out.hit_vpriv, newShotPoint, realRoot[2], newShotDir);
555  BU_LIST_INSERT(&(seghead->l), &(segp->l));
556  return 4; /* HIT */
557  /* XXX END CUT */
558 
559 }
560 
561 
562 /**
563  * Given ONE ray distance, return the normal and entry/exit point.
564  */
565 void
566 rt_superell_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
567 {
568  struct superell_specific *superell =
569  (struct superell_specific *)stp->st_specific;
570 
571  vect_t xlated;
572  fastf_t scale;
573 
574  VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
575  VSUB2(xlated, hitp->hit_point, superell->superell_V);
576  MAT4X3VEC(hitp->hit_normal, superell->superell_invRSSR, xlated);
577  scale = 1.0 / MAGNITUDE(hitp->hit_normal);
578  VSCALE(hitp->hit_normal, hitp->hit_normal, scale);
579 
580  return;
581 }
582 
583 
584 /**
585  * Return the curvature of the superellipsoid.
586  */
587 void
588 rt_superell_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
589 {
590  if (!cvp || !hitp || !stp)
591  return;
592  RT_CK_HIT(hitp);
593  RT_CK_SOLTAB(stp);
594 
595  bu_log("called rt_superell_curve()\n");
596  return;
597 }
598 
599 
600 /**
601  * For a hit on the surface of an SUPERELL, return the (u, v) coordinates
602  * of the hit point, 0 <= u, v <= 1.
603  * u = azimuth
604  * v = elevation
605  */
606 void
607 rt_superell_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
608 {
609  if (ap) RT_CK_APPLICATION(ap);
610  if (!stp || !hitp || !uvp)
611  return;
612  RT_CK_SOLTAB(stp);
613  RT_CK_HIT(hitp);
614 
615  bu_log("called rt_superell_uv()\n");
616  return;
617 }
618 
619 
620 void
622 {
623  struct superell_specific *superell =
624  (struct superell_specific *)stp->st_specific;
625 
626  BU_PUT(superell, struct superell_specific);
627 }
628 
629 
630 /**
631  * Also used by the TGC code
632  */
633 #define SUPERELLOUT(n) ov+(n-1)*3
634 void
636  fastf_t *V,
637  fastf_t *A,
638  fastf_t *B)
639 {
640  static fastf_t c, d, e, f, g, h;
641 
642  e = h = .92388; /* cos(22.5 degrees) */
643  c = d = M_SQRT1_2; /* cos(45 degrees) */
644  g = f = .382683; /* cos(67.5 degrees) */
645 
646  /*
647  * For angle theta, compute surface point as
648  *
649  * V + cos(theta) * A + sin(theta) * B
650  *
651  * note that sin(theta) is cos(90-theta); arguments in degrees.
652  */
653 
654  VADD2(SUPERELLOUT(1), V, A);
655  VJOIN2(SUPERELLOUT(2), V, e, A, f, B);
656  VJOIN2(SUPERELLOUT(3), V, c, A, d, B);
657  VJOIN2(SUPERELLOUT(4), V, g, A, h, B);
658  VADD2(SUPERELLOUT(5), V, B);
659  VJOIN2(SUPERELLOUT(6), V, -g, A, h, B);
660  VJOIN2(SUPERELLOUT(7), V, -c, A, d, B);
661  VJOIN2(SUPERELLOUT(8), V, -e, A, f, B);
662  VSUB2(SUPERELLOUT(9), V, A);
663  VJOIN2(SUPERELLOUT(10), V, -e, A, -f, B);
664  VJOIN2(SUPERELLOUT(11), V, -c, A, -d, B);
665  VJOIN2(SUPERELLOUT(12), V, -g, A, -h, B);
666  VSUB2(SUPERELLOUT(13), V, B);
667  VJOIN2(SUPERELLOUT(14), V, g, A, -h, B);
668  VJOIN2(SUPERELLOUT(15), V, c, A, -d, B);
669  VJOIN2(SUPERELLOUT(16), V, e, A, -f, B);
670 }
671 
672 
673 int
674 rt_superell_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))
675 {
676  int i;
677  struct rt_superell_internal *eip;
678  fastf_t top[16*3];
679  fastf_t middle[16*3];
680  fastf_t bottom[16*3];
681 
682  BU_CK_LIST_HEAD(vhead);
683  RT_CK_DB_INTERNAL(ip);
684  eip = (struct rt_superell_internal *)ip->idb_ptr;
685  RT_SUPERELL_CK_MAGIC(eip);
686 
687  rt_superell_16pts(top, eip->v, eip->a, eip->b);
688  rt_superell_16pts(bottom, eip->v, eip->b, eip->c);
689  rt_superell_16pts(middle, eip->v, eip->a, eip->c);
690 
691  RT_ADD_VLIST(vhead, &top[15*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
692  for (i=0; i<16; i++) {
693  RT_ADD_VLIST(vhead, &top[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
694  }
695 
696  RT_ADD_VLIST(vhead, &bottom[15*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
697  for (i=0; i<16; i++) {
698  RT_ADD_VLIST(vhead, &bottom[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
699  }
700 
701  RT_ADD_VLIST(vhead, &middle[15*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
702  for (i=0; i<16; i++) {
703  RT_ADD_VLIST(vhead, &middle[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
704  }
705  return 0;
706 }
707 
708 
710  struct shell *s;
711  struct rt_superell_internal *eip;
712  mat_t invRinvS;
713  mat_t invRoS;
715 };
716 
717 
720  int nverts;
721  struct vertex **vp;
722  vect_t *norms;
723  int nfaces;
724  struct faceuse **fu;
725 };
726 
727 
728 /**
729  * Tesssuperellate an superellipsoid.
730  *
731  * The strategy is based upon the approach of Jon Leech 3/24/89, from
732  * program "sphere", which generates a polygon mesh approximating a
733  * sphere by recursive subdivision. First approximation is an
734  * octahedron; each level of refinement increases the number of
735  * polygons by a factor of 4. Level 3 (128 polygons) is a good
736  * tradeoff if gouraud shading is used to render the database.
737  *
738  * At the start, points ABC lie on surface of the unit sphere. Pick
739  * DEF as the midpoints of the three edges of ABC. Normalize the new
740  * points to lie on surface of the unit sphere.
741  *
742  * 1
743  * B
744  * /\
745  * 3 / \ 4
746  * D /____\ E
747  * /\ /\
748  * / \ / \
749  * /____\/____\
750  * A F C
751  * 0 5 2
752  *
753  * Returns -
754  * -1 failure
755  * 0 OK. *r points to nmgregion that holds this tesssuperellation.
756  */
757 int
758 rt_superell_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *UNUSED(ttol), const struct bn_tol *UNUSED(tol))
759 {
760  if (r) NMG_CK_REGION(*r);
761  if (m) NMG_CK_MODEL(m);
762  if (ip) RT_CK_DB_INTERNAL(ip);
763 
764  bu_log("called rt_superell_tess()\n");
765  return -1;
766 }
767 
768 
769 /**
770  * Import an superellipsoid/sphere from the database format to the
771  * internal structure. Apply modeling transformations as wsuperell.
772  */
773 int
774 rt_superell_import4(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
775 {
776  struct rt_superell_internal *eip;
777  union record *rp;
778  fastf_t vec[3*4 + 2];
779 
780  if (dbip) RT_CK_DBI(dbip);
781 
782  BU_CK_EXTERNAL(ep);
783  rp = (union record *)ep->ext_buf;
784  /* Check record type */
785  if (rp->u_id != ID_SOLID) {
786  bu_log("rt_superell_import4(): defective record\n");
787  return -1;
788  }
789 
790  RT_CK_DB_INTERNAL(ip);
791  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
792  ip->idb_type = ID_SUPERELL;
793  ip->idb_meth = &OBJ[ID_SUPERELL];
794  BU_ALLOC(ip->idb_ptr, struct rt_superell_internal);
795 
796  eip = (struct rt_superell_internal *)ip->idb_ptr;
797  eip->magic = RT_SUPERELL_INTERNAL_MAGIC;
798 
799  /* Convert from database to internal format */
800  flip_fastf_float(vec, rp->s.s_values, 4, dbip->dbi_version < 0 ? 1 : 0);
801 
802  /* Apply modeling transformations */
803  if (mat == NULL) mat = bn_mat_identity;
804  MAT4X3PNT(eip->v, mat, &vec[0*3]);
805  MAT4X3VEC(eip->a, mat, &vec[1*3]);
806  MAT4X3VEC(eip->b, mat, &vec[2*3]);
807  MAT4X3VEC(eip->c, mat, &vec[3*3]);
808 
809  if (dbip->dbi_version < 0) {
810  eip->n = flip_dbfloat(rp->s.s_values[12]);
811  eip->e = flip_dbfloat(rp->s.s_values[13]);
812  } else {
813  eip->n = rp->s.s_values[12];
814  eip->e = rp->s.s_values[13];
815  }
816 
817  return 0; /* OK */
818 }
819 
820 
821 int
822 rt_superell_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
823 {
824  struct rt_superell_internal *tip;
825  union record *rec;
826 
827  if (dbip) RT_CK_DBI(dbip);
828 
829  RT_CK_DB_INTERNAL(ip);
830  if (ip->idb_type != ID_SUPERELL) return -1;
831  tip = (struct rt_superell_internal *)ip->idb_ptr;
832  RT_SUPERELL_CK_MAGIC(tip);
833 
834  BU_CK_EXTERNAL(ep);
835  ep->ext_nbytes = sizeof(union record);
836  ep->ext_buf = (uint8_t *)bu_calloc(1, ep->ext_nbytes, "superell external");
837  rec = (union record *)ep->ext_buf;
838 
839  rec->s.s_id = ID_SOLID;
840  rec->s.s_type = SUPERELL;
841 
842  /* NOTE: This also converts to dbfloat_t */
843  VSCALE(&rec->s.s_values[0], tip->v, local2mm);
844  VSCALE(&rec->s.s_values[3], tip->a, local2mm);
845  VSCALE(&rec->s.s_values[6], tip->b, local2mm);
846  VSCALE(&rec->s.s_values[9], tip->c, local2mm);
847 
848  printf("SUPERELL: %g %g\n", tip->n, tip->e);
849 
850  rec->s.s_values[12] = tip->n;
851  rec->s.s_values[13] = tip->e;
852 
853  return 0;
854 }
855 
856 
857 /**
858  * Import an superellipsoid/sphere from the database format to the
859  * internal structure. Apply modeling transformations as wsuperell.
860  */
861 int
862 rt_superell_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
863 {
864  struct rt_superell_internal *eip;
865 
866  /* must be double for import and export */
867  double vec[ELEMENTS_PER_VECT*4 + 2];
868 
869  if (dbip) RT_CK_DBI(dbip);
870 
871  RT_CK_DB_INTERNAL(ip);
872  BU_CK_EXTERNAL(ep);
873 
874  BU_ASSERT_LONG(ep->ext_nbytes, ==, SIZEOF_NETWORK_DOUBLE * (ELEMENTS_PER_VECT*4 + 2));
875 
876  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
877  ip->idb_type = ID_SUPERELL;
878  ip->idb_meth = &OBJ[ID_SUPERELL];
879  BU_ALLOC(ip->idb_ptr, struct rt_superell_internal);
880 
881  eip = (struct rt_superell_internal *)ip->idb_ptr;
882  eip->magic = RT_SUPERELL_INTERNAL_MAGIC;
883 
884  /* Convert from database (network) to internal (host) format */
885  bu_cv_ntohd((unsigned char *)vec, ep->ext_buf, ELEMENTS_PER_VECT*4 + 2);
886 
887  /* Apply modeling transformations */
888  if (mat == NULL) mat = bn_mat_identity;
889  MAT4X3PNT(eip->v, mat, &vec[0*ELEMENTS_PER_VECT]);
890  MAT4X3VEC(eip->a, mat, &vec[1*ELEMENTS_PER_VECT]);
891  MAT4X3VEC(eip->b, mat, &vec[2*ELEMENTS_PER_VECT]);
892  MAT4X3VEC(eip->c, mat, &vec[3*ELEMENTS_PER_VECT]);
893  eip->n = vec[4*ELEMENTS_PER_VECT];
894  eip->e = vec[4*ELEMENTS_PER_VECT + 1];
895 
896  return 0; /* OK */
897 }
898 
899 
900 /**
901  * The external format is:
902  * V point
903  * A vector
904  * B vector
905  * C vector
906  */
907 int
908 rt_superell_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
909 {
910  struct rt_superell_internal *eip;
911 
912  /* must be double for import and export */
913  double vec[ELEMENTS_PER_VECT*4 + 2];
914 
915  if (dbip) RT_CK_DBI(dbip);
916 
917  RT_CK_DB_INTERNAL(ip);
918  if (ip->idb_type != ID_SUPERELL) return -1;
919  eip = (struct rt_superell_internal *)ip->idb_ptr;
920  RT_SUPERELL_CK_MAGIC(eip);
921 
922  BU_CK_EXTERNAL(ep);
923  ep->ext_nbytes = SIZEOF_NETWORK_DOUBLE * (ELEMENTS_PER_VECT*4 + 2);
924  ep->ext_buf = (uint8_t *)bu_malloc(ep->ext_nbytes, "superell external");
925 
926  /* scale 'em into local buffer */
927  VSCALE(&vec[0*ELEMENTS_PER_VECT], eip->v, local2mm);
928  VSCALE(&vec[1*ELEMENTS_PER_VECT], eip->a, local2mm);
929  VSCALE(&vec[2*ELEMENTS_PER_VECT], eip->b, local2mm);
930  VSCALE(&vec[3*ELEMENTS_PER_VECT], eip->c, local2mm);
931 
932  vec[4*ELEMENTS_PER_VECT] = eip->n;
933  vec[4*ELEMENTS_PER_VECT + 1] = eip->e;
934 
935  /* Convert from internal (host) to database (network) format */
936  bu_cv_htond(ep->ext_buf, (unsigned char *)vec, ELEMENTS_PER_VECT*4 + 2);
937 
938  return 0;
939 }
940 
941 
942 /**
943  * Make human-readable formatted presentation of this solid. First
944  * line describes type of solid. Additional lines are indented one
945  * tab, and give parameter values.
946  */
947 int
948 rt_superell_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
949 {
950  struct rt_superell_internal *tip =
951  (struct rt_superell_internal *)ip->idb_ptr;
952  fastf_t mag_a, mag_b, mag_c;
953  char buf[256];
954  double angles[5];
955  vect_t unitv;
956 
957  RT_SUPERELL_CK_MAGIC(tip);
958  bu_vls_strcat(str, "superellipsoid (SUPERELL)\n");
959 
960  sprintf(buf, "\tV (%g, %g, %g)\n",
961  INTCLAMP(tip->v[X] * mm2local),
962  INTCLAMP(tip->v[Y] * mm2local),
963  INTCLAMP(tip->v[Z] * mm2local));
964  bu_vls_strcat(str, buf);
965 
966  mag_a = MAGNITUDE(tip->a);
967  mag_b = MAGNITUDE(tip->b);
968  mag_c = MAGNITUDE(tip->c);
969 
970  sprintf(buf, "\tA (%g, %g, %g) mag=%g\n",
971  INTCLAMP(tip->a[X] * mm2local),
972  INTCLAMP(tip->a[Y] * mm2local),
973  INTCLAMP(tip->a[Z] * mm2local),
974  INTCLAMP(mag_a * mm2local));
975  bu_vls_strcat(str, buf);
976 
977  sprintf(buf, "\tB (%g, %g, %g) mag=%g\n",
978  INTCLAMP(tip->b[X] * mm2local),
979  INTCLAMP(tip->b[Y] * mm2local),
980  INTCLAMP(tip->b[Z] * mm2local),
981  INTCLAMP(mag_b * mm2local));
982  bu_vls_strcat(str, buf);
983 
984  sprintf(buf, "\tC (%g, %g, %g) mag=%g\n",
985  INTCLAMP(tip->c[X] * mm2local),
986  INTCLAMP(tip->c[Y] * mm2local),
987  INTCLAMP(tip->c[Z] * mm2local),
988  INTCLAMP(mag_c * mm2local));
989  bu_vls_strcat(str, buf);
990 
991  sprintf(buf, "\t<n, e> (%g, %g)\n", INTCLAMP(tip->n), INTCLAMP(tip->e));
992  bu_vls_strcat(str, buf);
993 
994  if (!verbose) return 0;
995 
996  VSCALE(unitv, tip->a, 1/mag_a);
997  rt_find_fallback_angle(angles, unitv);
998  rt_pr_fallback_angle(str, "\tA", angles);
999 
1000  VSCALE(unitv, tip->b, 1/mag_b);
1001  rt_find_fallback_angle(angles, unitv);
1002  rt_pr_fallback_angle(str, "\tB", angles);
1003 
1004  VSCALE(unitv, tip->c, 1/mag_c);
1005  rt_find_fallback_angle(angles, unitv);
1006  rt_pr_fallback_angle(str, "\tC", angles);
1007 
1008  return 0;
1009 }
1010 
1011 
1012 /**
1013  * Free the storage associated with the rt_db_internal version of this
1014  * solid.
1015  */
1016 void
1018 {
1019  RT_CK_DB_INTERNAL(ip);
1020 
1021  bu_free(ip->idb_ptr, "superell ifree");
1022  ip->idb_ptr = ((void *)0);
1023 }
1024 
1025 
1026 /* The U parameter runs south to north. In order to orient loop CCW,
1027  * need to start with 0, 1-->0, 0 transition at the south pole.
1028  */
1029 /* unused var:
1030 static const fastf_t rt_superell_uvw[5*ELEMENTS_PER_VECT] = {
1031  0, 1, 0,
1032  0, 0, 0,
1033  1, 0, 0,
1034  1, 1, 0,
1035  0, 1, 0
1036 };
1037 */
1038 
1039 int
1040 rt_superell_params(struct pc_pc_set *UNUSED(ps), const struct rt_db_internal *ip)
1041 {
1042  if (ip) RT_CK_DB_INTERNAL(ip);
1043 
1044  return 0; /* OK */
1045 }
1046 
1047 
1048 /**
1049  * Computes the volume of a superellipsoid
1050  *
1051  * Volume equation from http://lrv.fri.uni-lj.si/~franc/SRSbook/geometry.pdf
1052  * which also includes a derivation on page 32.
1053  */
1054 void
1055 rt_superell_volume(fastf_t *volume, const struct rt_db_internal *ip)
1056 {
1057 #ifdef HAVE_TGAMMA
1058  struct rt_superell_internal *sip;
1059  double mag_a, mag_b, mag_c;
1060 #endif
1061 
1062  if (volume == NULL || ip == NULL) {
1063  return;
1064  }
1065 
1066 #ifdef HAVE_TGAMMA
1067  RT_CK_DB_INTERNAL(ip);
1068  sip = (struct rt_superell_internal *)ip->idb_ptr;
1069  RT_SUPERELL_CK_MAGIC(sip);
1070 
1071  mag_a = MAGNITUDE(sip->a);
1072  mag_b = MAGNITUDE(sip->b);
1073  mag_c = MAGNITUDE(sip->c);
1074 
1075  *volume = 2.0 * mag_a * mag_b * mag_c * sip->e * sip->n * (tgamma(sip->n/2.0 + 1.0) * tgamma(sip->n) / tgamma(3.0 * sip->n/2.0 + 1.0)) * (tgamma(sip->e / 2.0) * tgamma(sip->e / 2.0) / tgamma(sip->e));
1076 #endif
1077 }
1078 
1079 
1080 static fastf_t
1081 superell_surf_area_box(vect_t mags)
1082 {
1083  return 8 * (mags[0] * mags[1] + mags[0] * mags[2] + mags[1] * mags[2]);
1084 }
1085 
1086 
1087 static fastf_t
1088 sign(fastf_t x)
1089 {
1090  if (x > 0) {
1091  return 1.0;
1092  } else if (x < 0) {
1093  return -1.0;
1094  } else {
1095  return 0.0;
1096  }
1097 }
1098 
1099 
1100 /* This is the c auxiliary function for superell_xyz_from_uv.
1101  * See http://en.wikipedia.org/wiki/Superellipsoid for more
1102  * information.
1103  */
1104 static fastf_t
1105 superell_c(fastf_t w, fastf_t m)
1106 {
1107  return sign(cos(w)) * pow(fabs(cos(w)), m);
1108 }
1109 
1110 
1111 /* This is the s auxiliary function for superell_xyz_from_uv.
1112  * See http://en.wikipedia.org/wiki/Superellipsoid for more
1113  * information.
1114  */
1115 static fastf_t
1116 superell_s(fastf_t w, fastf_t m)
1117 {
1118  return sign(sin(w)) * pow(fabs(sin(w)), m);
1119 }
1120 
1121 
1122 /* This calculates the xyz coordinates of a set of uv coordinates on
1123  * a superellipsoid, using the algorithm detailed in these places:
1124  *
1125  * http://en.wikipedia.org/wiki/Superellipsoid
1126  * http://paulbourke.net/geometry/superellipse/
1127  *
1128  * Since this code is called inside a loop through v-values inside a
1129  * loop through u-values, superell_c(u, sip->e) and superell_s(u,
1130  * sip->e) can be precomputed, which results in a significant
1131  * performance gain; for this reason, the function takes there values
1132  * and not he original u-value.
1133  */
1134 static void
1135 superell_xyz_from_uv(point_t pt, fastf_t u, fastf_t v, vect_t mags, const struct rt_superell_internal *sip)
1136 {
1137  pt[X] = mags[0] * superell_c(v, sip->n) * superell_c(u, sip->e);
1138  pt[Y] = mags[1] * superell_c(v, sip->n) * superell_s(u, sip->e);
1139  pt[Z] = mags[2] * superell_s(v, sip->n);
1140 }
1141 
1142 
1143 /* This attempts to find the surface area by subdividing the uv plane
1144  * into squares, and finding the approximate surface area of each
1145  * square.
1146  */
1147 static fastf_t
1148 superell_surf_area_general(const struct rt_superell_internal *sip, vect_t mags, fastf_t side_length)
1149 {
1150  fastf_t area = 0;
1151  fastf_t u, v;
1152 
1153  /* There are M_PI / side_length + 1 values because both ends have
1154  to be stored */
1155  int row_length = sizeof(point_t) * (M_PI / side_length + 1);
1156 
1157  point_t *row1 = (point_t *)bu_malloc(row_length, "superell_surf_area_general");
1158  point_t *row2 = (point_t *)bu_malloc(row_length, "superell_surf_area_general");
1159 
1160  int idx = 0;
1161 
1162 
1163  /* This function keeps a moving window of two rows at any time,
1164  * and calculates the area of space between those two rows. It
1165  * calculates the first row outside of the loop so that at each
1166  * iteration it only has to calculate the second. Following
1167  * standard definitions of u and v, u ranges from -pi to pi, and v
1168  * ranges from -pi/2 to pi/2. Using an extra index variable allows
1169  * the code to compute the index into the array very efficiently.
1170  */
1171  for (v = -M_PI_2; v < M_PI_2; v += side_length, idx++) {
1172  superell_xyz_from_uv(row1[idx], -M_PI, v, mags, sip);
1173  }
1174 
1175 
1176  /* This starts at -M_PI + side_length because the first row is
1177  * computed outside the loop, which allows the loop to always
1178  * calculate the second row of the pair.
1179  */
1180  for (u = -M_PI + side_length; u < M_PI; u += side_length) {
1181  idx = 0;
1182  for (v = -M_PI_2; v < M_PI_2; v += side_length, idx++) {
1183  superell_xyz_from_uv(row2[idx], u + side_length, v, mags, sip);
1184  }
1185 
1186  idx = 0;
1187 
1188  /* This ends at -M_PI / 2 - side_length because if it kept
1189  * going it would overflow the array, since it always looks at
1190  * the square to the right of its current index.
1191  */
1192  for (v = - M_PI_2; v < M_PI_2 - side_length; v += side_length, idx++) {
1193  area +=
1194  bn_dist_pt3_pt3(row1[idx], row1[idx + 1]) *
1195  bn_dist_pt3_pt3(row1[idx], row2[idx]);
1196  }
1197 
1198  memcpy(row1, row2, row_length);
1199  }
1200 
1201  bu_free(row1, "superell_surf_area_general");
1202  bu_free(row2, "superell_surf_area_general");
1203 
1204  return area;
1205 }
1206 
1207 
1208 void
1210 {
1211  struct rt_superell_internal *sip;
1212  vect_t mags;
1213 
1214  RT_CK_DB_INTERNAL(ip);
1215  sip = (struct rt_superell_internal *)ip->idb_ptr;
1216 
1217  mags[0] = MAGNITUDE(sip->a);
1218  mags[1] = MAGNITUDE(sip->b);
1219  mags[2] = MAGNITUDE(sip->c);
1220 
1221  /* The parametric representation does not work correctly at n = e
1222  * = 0, so this uses a special calculation for boxes.
1223  */
1224  if ((NEAR_EQUAL(sip->e, 0, BN_TOL_DIST)) && NEAR_EQUAL(sip->n, 0, BN_TOL_DIST)) {
1225  *area = superell_surf_area_box(mags);
1226  } else {
1227  /* This number specifies the initial precision used. The
1228  * precision is roughly the number of chunks that the uv-space
1229  * is divided into during the approximation; the larger the
1230  * number the higher the accuracy and the lower the
1231  * performance.
1232  */
1233  int precision = 1024;
1234 
1235  fastf_t previous_area = 0;
1236  fastf_t current_area = superell_surf_area_general(sip, mags, M_PI / precision);
1237  while (!(NEAR_EQUAL(current_area, previous_area, BN_TOL_DIST))) {
1238  /* The precision is multiplied by a constant on each round
1239  * through the loop to make sure that the precision
1240  * continues to increase until a satisfactory value is
1241  * found. If this value is very small, this approximation
1242  * will likely converge fairly quickly and with lower
1243  * accuracy: the smaller the distance between the inputs
1244  * the more likely that the outputs will be within
1245  * BN_TOL_DIST. A large value will result in slow
1246  * convergence, but in high accuracy: when the values are
1247  * finally within BN_TOL_DIST of each other the
1248  * approximation must be very good, because a large
1249  * increase in input precision resulted in a small
1250  * increase in accuracy.
1251  */
1252  precision *= 2;
1253  previous_area = current_area;
1254  current_area = superell_surf_area_general(sip, mags, M_PI / precision);
1255  }
1256  *area = current_area;
1257  }
1258 }
1259 
1260 
1261 /*
1262  * Local Variables:
1263  * mode: C
1264  * tab-width: 8
1265  * indent-tabs-mode: t
1266  * c-file-style: "stroustrup"
1267  * End:
1268  * ex: shiftwidth=4 tabstop=8
1269  */
HIDDEN int sign(double val)
Definition: brep.cpp:1145
char * d_namep
pointer to name string
Definition: raytrace.h:859
Definition: raytrace.h:800
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
Definition: list.h:118
const struct directory * st_dp
Directory entry of solid.
Definition: raytrace.h:436
void rt_superell_volume(fastf_t *volume, const struct rt_db_internal *ip)
Definition: superell.c:1055
int rt_superell_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
Definition: superell.c:396
#define RT_CK_APPLICATION(_p)
Definition: raytrace.h:1675
vect_t superell_Bu
Definition: superell.c:159
int rt_superell_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: superell.c:254
double dist
>= 0
Definition: tol.h:73
const mat_t bn_mat_identity
Matrix and vector functionality.
Definition: mat.c:46
int rt_superell_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: superell.c:674
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
Definition: raytrace.h:215
void bn_pr_roots(const char *title, const struct bn_complex roots[], int n)
void rt_superell_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
Definition: superell.c:607
#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
Definition: raytrace.h:248
double superell_invmsCu
Definition: superell.c:165
fastf_t st_aradius
Radius of APPROXIMATING sphere.
Definition: raytrace.h:433
vect_t superell_Au
Definition: superell.c:158
Header file for the BRL-CAD common definitions.
Definition: poly.h:47
struct shell * s
Definition: superell.c:710
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
struct faceuse ** fu
Definition: superell.c:724
void bu_cv_htond(unsigned char *out, const unsigned char *in, size_t count)
const struct bu_structparse rt_superell_parse[]
Definition: superell.c:50
struct bu_list l
Definition: raytrace.h:369
double superell_invmsBu
Definition: superell.c:164
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
if(share_geom)
Definition: nmg_mod.c:3829
int idb_major_type
Definition: raytrace.h:192
Definition: color.c:49
void rt_superell_print(const struct soltab *stp)
Definition: superell.c:372
double superell_e
Definition: superell.c:162
int rt_superell_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
Definition: superell.c:862
vect_t superell_Cu
Definition: superell.c:160
#define RT_ADD_VLIST(hd, pnt, draw)
Definition: raytrace.h:1865
#define RT_CK_DB_INTERNAL(_p)
Definition: raytrace.h:207
double superell_n
Definition: superell.c:161
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
fastf_t st_bradius
Radius of BOUNDING sphere.
Definition: raytrace.h:434
#define RT_CK_HIT(_p)
Definition: raytrace.h:259
void rt_superell_16pts(fastf_t *ov, fastf_t *V, fastf_t *A, fastf_t *B)
Definition: superell.c:635
mat_t superell_SoR
Definition: superell.c:167
#define BN_VLIST_LINE_MOVE
Definition: vlist.h:82
const struct rt_functab * idb_meth
for ft_ifree(), etc.
Definition: raytrace.h:194
fastf_t theta_tol
Definition: superell.c:714
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
uint8_t * ext_buf
Definition: parse.h:216
static void top()
#define BU_GET(_ptr, _type)
Definition: malloc.h:201
point_t hit_point
DEPRECATED: Intersection point, use VJOIN1 hit_dist.
Definition: raytrace.h:251
point_t st_max
max X, Y, Z of bounding RPP
Definition: raytrace.h:438
void rt_pr_fallback_angle(struct bu_vls *str, const char *prefix, const double angles[5])
void rt_find_fallback_angle(double angles[5], const vect_t vec)
#define BN_VLIST_LINE_DRAW
Definition: vlist.h:83
#define BN_TOL_DIST
Definition: tol.h:109
fastf_t flip_dbfloat(dbfloat_t d)
Definition: db_flip.c:58
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
mat_t superell_invRSSR
Definition: superell.c:168
void bn_mat_mul(mat_t o, const mat_t a, const mat_t b)
void rt_superell_surf_area(fastf_t *area, const struct rt_db_internal *ip)
Definition: superell.c:1209
Support for uniform tolerances.
Definition: tol.h:71
vect_t superell_V
Definition: superell.c:157
#define bu_offsetofarray(_t, _a, _d, _i)
Definition: parse.h:65
#define BU_STRUCTPARSE_FUNC_NULL
Definition: parse.h:153
vect_t r_dir
Direction of ray (UNIT Length)
Definition: raytrace.h:219
mat_t superell_invR
Definition: superell.c:169
void rt_superell_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
Definition: superell.c:588
struct bn_tol rti_tol
Math tolerances for this model.
Definition: raytrace.h:1765
int rt_superell_params(struct pc_pc_set *ps, const struct rt_db_internal *ip)
Definition: superell.c:1040
#define bu_offsetof(_t, _m)
Definition: parse.h:64
#define RT_CK_DBI(_p)
Definition: raytrace.h:829
#define RT_GET_SEG(p, res)
Definition: raytrace.h:379
int rt_superell_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *tol)
Definition: superell.c:177
void * idb_ptr
Definition: raytrace.h:195
double bn_dist_pt3_pt3(const point_t a, const point_t b)
Returns distance between two points.
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
int rt_poly_roots(bn_poly_t *eqn, bn_complex_t roots[], const char *name)
struct rt_superell_internal * eip
Definition: superell.c:711
struct vertex ** vp
Definition: superell.c:721
#define R
Definition: msr.c:55
#define ID_SUPERELL
Superquadratic ellipsoid.
Definition: raytrace.h:507
int rt_superell_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
Definition: superell.c:948
#define RT_CK_SOLTAB(_p)
Definition: raytrace.h:453
void * st_specific
-> ID-specific (private) struct
Definition: raytrace.h:435
int rt_superell_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: superell.c:822
#define A
Definition: msr.c:51
mat_t invRinvS
Definition: superell.c:712
Definition: color.c:51
int dbi_version
PRIVATE: use db_version()
Definition: raytrace.h:824
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
vect_t hit_normal
DEPRECATED: Surface Normal at hit_point, use RT_HIT_NORMAL.
Definition: raytrace.h:252
#define BU_CK_LIST_HEAD(_p)
Definition: list.h:142
mat_t invRoS
Definition: superell.c:713
#define BU_CK_EXTERNAL(_p)
Definition: parse.h:224
Complex numbers.
Definition: complex.h:39
double superell_invmsAu
Definition: superell.c:163
size_t ext_nbytes
Definition: parse.h:210
int rt_superell_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
Definition: superell.c:758
void rt_superell_free(struct soltab *stp)
Definition: superell.c:621
void rt_superell_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
Definition: superell.c:566
void rt_superell_ifree(struct rt_db_internal *ip)
Definition: superell.c:1017
fastf_t hit_dist
dist from r_pt to hit_point
Definition: raytrace.h:250
HIDDEN void verbose(struct human_data_t *dude)
Definition: human.c:2008
Definition: vls.h:56
#define M_SQRT1_2
Definition: fft.h:38
double fastf_t
Definition: defines.h:300
#define VPRINT(a, b)
Definition: raytrace.h:1881
int rt_superell_import4(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
Definition: superell.c:774
vect_t superell_invsq
Definition: superell.c:166
int rt_superell_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: superell.c:908
#define RT_SUPERELL_INTERNAL_MAGIC
Definition: magic.h:110
Definition: color.c:50
#define SUPERELLOUT(n)
Definition: superell.c:633
point_t st_center
Centroid of solid.
Definition: raytrace.h:432