BRL-CAD
hrt.c
Go to the documentation of this file.
1 /* H R T . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2013-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/hrt/hrt.c
23  *
24  * Intersect a ray with a Heart
25  *
26  * Algorithm:
27  *
28  * Given V, X, Y, Z and d, there is a set of points on this heart
29  *
30  * { (x, y, z) | (x, y, z) is on heart defined by V, X, Y, Z and d }
31  *
32  * Through a series of Transformations, this set will be transformed
33  * into a set of points on a heart centered at the origin
34  * which lies on the X-Y plane (i.e., Z is on the Z axis).
35  *
36  * { (x', y', z') | (x', y', z') is on heart at origin }
37  *
38  * The transformation from X to X' is accomplished by:
39  *
40  * X' = S(R(X - V))
41  *
42  * where R(X) = (Z/(|X|))
43  * (Y/(|Y|)) . X
44  * (Z/(|Z|))
45  *
46  * and S(X) = (1/|A| 0 0)
47  * (0 1/|A| 0) . X
48  * (0 0 1/|A|)
49  * where |A| = d
50  *
51  * To find the intersection of a line with the heart, consider the
52  * parametric line L:
53  *
54  * L : { P(n) | P + t(n) . D }
55  *
56  * Call W the actual point of intersection between L and the heart.
57  * Let W' be the point of intersection between L' and the heart.
58  *
59  * L' : { P'(n) | P' + t(n) . D' }
60  *
61  * W = invR(invS(W')) + V
62  *
63  * Where W' = k D' + P'.
64  *
65  * Given a line, find the equation of the
66  * heart in terms of the variable 't'.
67  *
68  * The equation for the heart is:
69  *
70  * [ X**2 + 9/4*Y**2 + Z**2 - 1 ]**3 - Z**3 * (X**2 + 9/80 * Y**2) = 0
71  *
72  * First, find X, Y, and Z in terms of 't' for this line, then
73  * substitute them into the equation above.
74  *
75  * Wx = Dx*t + Px
76  *
77  * Wx**3 = Dx**3 * t**3 + 3 * Dx**2 * Px * t**2 + 3 * Dx * Px**2 * t + Px**3
78  *
79  * The real roots of the equation in 't' are the intersect points
80  * along the parametric line.
81  *
82  * NORMALS. Given the point W on the heart, what is the vector normal
83  * to the tangent plane at that point?
84  *
85  * Map W onto the heart, i.e.: W' = S(R(W - V)). In this case,
86  * we find W' by solving the parametric line given k.
87  *
88  * The gradient of the heart at W' is in fact the normal vector.
89  *
90  * Given that the sextic equation for the heart is:
91  *
92  * [ X**2 + 9/4 * Y**2 + Z**2 - 1 ]**3 - Z**3 * (X**2 + 9/80 * Y**2) = 0.
93  *
94  * let w = X**2 + 9/4*Y**2 + Z**2 - 1 , then the sextic equation becomes:
95  *
96  * w**3 - Y**3 * (X**2 + 9/80 * Y**2) = 0.
97  *
98  * For f(x, y, z) = 0, the gradient of f() is (df/dx, df/dy, df/dz).
99  *
100  * df/dx = 6 * x * (w**2 - y**3)
101  * df/dy = 6 * (12/27 * w**2 * y**2 - 1/2 * y**2 * (x**2 + 9 / 80*z**3))
102  * df/dz = 6 * (w**2 * z - 160 / 9 * y**3 * z**2)
103  *
104  * Note that the normal vector produced above will not have
105  * length. Also, to make this useful for the original heart, it will
106  * have to be rotated back to the orientation of the original heart.
107  *
108  */
109 /** @} */
110 
111 #include "common.h"
112 
113 #include <stddef.h>
114 #include <string.h>
115 #include <math.h>
116 #include "bio.h"
117 
118 #include "bu/cv.h"
119 #include "vmath.h"
120 #include "db.h"
121 #include "nmg.h"
122 #include "rtgeom.h"
123 #include "raytrace.h"
124 
125 #include "../../librt_private.h"
126 
127 
128 const struct bu_structparse rt_hrt_parse[] = {
129  { "%f", 3, "V", bu_offsetofarray(struct rt_hrt_internal, v, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
130  { "%f", 3, "X", bu_offsetofarray(struct rt_hrt_internal, xdir, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
131  { "%f", 3, "Y", bu_offsetofarray(struct rt_hrt_internal, ydir, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
132  { "%f", 3, "Z", bu_offsetofarray(struct rt_hrt_internal, zdir, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
133  { "%g", 1, "d", bu_offsetof(struct rt_hrt_internal, d), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
134  { {'\0', '\0', '\0', '\0'}, 0, (char *)NULL, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL }
135 };
136 
137 
138 struct hrt_specific {
139  vect_t hrt_V; /* Vector to center of heart */
140  vect_t hrt_X; /* unit-length X vector */
141  vect_t hrt_Y; /* unit-length Y vector */
142  vect_t hrt_Z; /* unit-length Z vector */
143  fastf_t hrt_d; /* for distance to upper and lower cusps */
144  vect_t hrt_invsq; /* [ 1.0 / |X|**2, 1.0 / |Y|**2, 1.0 / |Z|**2 ] */
145  mat_t hrt_SoR; /* Scale(Rot(vect)) */
146  mat_t hrt_invRSSR; /* invRot(Scale(Scale(vect))) */
147  mat_t hrt_invR; /* transposed rotation matrix */
148 };
149 
150 
151 /**
152  * Compute the bounding RPP for a heart.
153  */
154 int
155 rt_hrt_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol))
156 {
157  struct rt_hrt_internal *hip;
158  fastf_t magsq_x, magsq_y, magsq_z;
159  vect_t Xu, Yu, Zu;
160  mat_t R;
161  fastf_t f;
162 
163  hip = (struct rt_hrt_internal *)ip->idb_ptr;
164  RT_HRT_CK_MAGIC(hip);
165 
166  magsq_x = MAGSQ(hip->xdir);
167  magsq_y = MAGSQ(hip->ydir);
168  magsq_z = MAGSQ(hip->zdir);
169 
170  /* Create unit length versions of X, Y, Z */
171  f = 1.0/sqrt(magsq_x);
172  VSCALE(Xu, hip->xdir, f);
173  f = 1.0/sqrt(magsq_y);
174  VSCALE(Yu, hip->ydir, f);
175  f = 1.0/sqrt(magsq_z);
176  VSCALE(Zu, hip->zdir, f);
177 
178  MAT_IDN(R);
179  VMOVE(&R[0], Xu);
180  VMOVE(&R[4], Yu);
181  VMOVE(&R[8], Zu);
182 
183  /**
184  * View-points from directly above and besides the heart indicate that it is
185  * elliptical, Y-axis is the major axis, X-axis the minor axis and the
186  * ratio of the radius of the minor axis to the radius of the major axis
187  * is 2/3
188  */
189  /* X */
190  f = hip->xdir[X];
191  (*min)[X] = hip->v[X] - f * 2/3;
192  (*max)[X] = hip->v[X] + f * 2/3;
193 
194  /* Y */
195  f = hip->ydir[Y];
196  (*min)[Y] = hip->v[Y] - f;
197  (*max)[Y] = hip->v[Y] + f;
198 
199  /* Z */
200  f = hip->zdir[Z];
201  /**
202  * 1.0 stands for the length of zdir vector and 0.25 closely approximates
203  * some value which encloses the displacement from upper cusp to highest
204  * point of either lobe in the Z direction
205  */
206  (*min)[Z] = hip->v[Z] - f;
207  (*max)[Z] = hip->v[Z] + f * 1.25;
208 
209  return 0;
210 }
211 
212 
213 /**
214  * Given a pointer to a GED database record, and a transformation
215  * matrix, determine if this is a valid heart, and if so, precompute
216  * various terms of the formula.
217  *
218  * Returns -
219  * 0 HRT is OK
220  * !0 Error in description
221  *
222  * Implicit return -
223  * A struct hrt_specific is created, and its address is stored in
224  * stp->st_specific for use by rt_hrt_shot().
225  */
226 int
227 rt_hrt_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
228 {
229  register struct hrt_specific *hrt;
230  struct rt_hrt_internal *hip;
231  fastf_t magsq_x, magsq_y, magsq_z;
232  mat_t R, TEMP;
233  vect_t Xu, Yu, Zu;
234  fastf_t f;
235 
236  hip = (struct rt_hrt_internal *)ip->idb_ptr;
237  RT_HRT_CK_MAGIC(hip);
238 
239  /* Validate that |X| > 0, |Y| > 0 and |Z| > 0 */
240  magsq_x = MAGSQ(hip->xdir);
241  magsq_y = MAGSQ(hip->ydir);
242  magsq_z = MAGSQ(hip->zdir);
243 
244  if (magsq_x < rtip->rti_tol.dist_sq || magsq_y < rtip->rti_tol.dist_sq || magsq_z < rtip->rti_tol.dist_sq) {
245  bu_log("rt_hrt_prep(): hrt(%s) near-zero length X(%g), Y(%g), or Z(%g) vector\n",
246  stp->st_name, magsq_x, magsq_y, magsq_z);
247  return 1;
248  }
249  if (hip->d < rtip->rti_tol.dist) {
250  bu_log("rt_hrt_prep(): hrt(%s) near-zero length <d> distance to cusps (%g) causes problems\n", stp->st_name, hip->d);
251  return 1;
252  }
253  if (hip->d > 10000.0) {
254  bu_log("rt_hrt_prep(): hrt(%s) very large <d> distance to cusps (%g) causes problems\n", stp->st_name, hip->d);
255  /* BAD */
256  }
257 
258  /* Create unit length versions of X, Y, Z */
259  f = 1.0/sqrt(magsq_x);
260  VSCALE(Xu, hip->xdir, f);
261  f = 1.0/sqrt(magsq_y);
262  VSCALE(Yu, hip->ydir, f);
263  f = 1.0/sqrt(magsq_z);
264  VSCALE(Zu, hip->zdir, f);
265 
266  /* Check that X, Y, Z are perpendicular to each other */
267  f = VDOT(Xu, Yu);
268  if (! NEAR_ZERO(f, rtip->rti_tol.dist)) {
269  bu_log("rt_hrt_prep(): hrt(%s) X not perpendicular to Y, f=%f\n", stp->st_name, f);
270  return 1;
271  }
272  f = VDOT(Xu, Zu);
273  if (! NEAR_ZERO(f, rtip->rti_tol.dist)) {
274  bu_log("rt_hrt_prep(): hrt(%s) X not perpendicular to Z, f=%f\n", stp->st_name, f);
275  return 1;
276  }
277  f = VDOT(Zu, Yu);
278  if (! NEAR_ZERO(f, rtip->rti_tol.dist)) {
279  bu_log("rt_hrt_prep(): hrt(%s) Z not perpendicular to Y, f=%f\n", stp->st_name, f);
280  return 1;
281  }
282 
283  /* Scale X and Y */
284  f = 0.8/sqrt(magsq_x);
285  VSCALE(Xu, hip->xdir, f);
286  f = 0.8/sqrt(magsq_y);
287  VSCALE(Yu, hip->ydir, f);
288 
289  /* Solid is OK, compute constant terms now */
290 
291  BU_GET(hrt, struct hrt_specific);
292  stp->st_specific = (void *)hrt;
293 
294  hrt->hrt_d = hip->d;
295 
296  VMOVE(hrt->hrt_V, hip->v);
297 
298  VSET(hrt->hrt_invsq, 1.0/magsq_x, 1.0/magsq_y, 1.0/magsq_z);
299  VMOVE(hrt->hrt_X, Xu);
300  VMOVE(hrt->hrt_Y, Yu);
301  VMOVE(hrt->hrt_Z, Zu);
302 
303  /* compute the rotation matrix */
304  MAT_IDN(R);
305  VMOVE(&R[0], Xu);
306  VMOVE(&R[4], Yu);
307  VMOVE(&R[8], Zu);
308  bn_mat_trn(hrt->hrt_invR, R);
309 
310  /* compute invRSSR */
311  MAT_IDN(hrt->hrt_invRSSR);
312  MAT_IDN(TEMP);
313  TEMP[0] = hrt->hrt_invsq[0];
314  TEMP[5] = hrt->hrt_invsq[1];
315  TEMP[10] = hrt->hrt_invsq[2];
316  bn_mat_mul(TEMP, TEMP, R);
317  bn_mat_mul(hrt->hrt_invRSSR, hrt->hrt_invR, TEMP);
318 
319  /* compute Scale(Rotate(vect)) */
320  MAT_IDN(hrt->hrt_SoR);
321  VSCALE(&hrt->hrt_SoR[0], hip->xdir, hrt->hrt_invsq[0]);
322  VSCALE(&hrt->hrt_SoR[4], hip->ydir, hrt->hrt_invsq[1]);
323  VSCALE(&hrt->hrt_SoR[8], hip->zdir, hrt->hrt_invsq[2]);
324 
325  /* compute bounding sphere */
326  VMOVE(stp->st_center, hrt->hrt_V);
327 
328  /**
329  * 1.0 stands for the length of zdir vector and 0.25 closely approximates
330  * some value which encloses the displacement from upper cusp to highest
331  * point of either lobe in the Z direction
332  */
333  f = hip->zdir[Z] * 1.25;
334  stp->st_aradius = stp->st_bradius = sqrt(f);
335 
336  /* compute bounding RPP */
337  if (rt_hrt_bbox(ip, &(stp->st_min), &(stp->st_max), &rtip->rti_tol)) return 1;
338 
339  return 0; /* OK */
340 }
341 
342 
343 void
344 rt_hrt_print(register const struct soltab *stp)
345 {
346  register struct hrt_specific *hrt =
347  (struct hrt_specific *)stp->st_specific;
348 
349  VPRINT("V", hrt->hrt_V);
350  bn_mat_print("S o R", hrt->hrt_SoR);
351  bn_mat_print("invRSSR", hrt->hrt_invRSSR);
352 }
353 
354 
355 /**
356  * Intersect a ray with a heart, where all constant terms have been
357  * precomputed by rt_hrt_prep(). If an intersection occurs, one or
358  * two struct seg(s) will be acquired and filled in.
359  *
360  * NOTE: All lines in this function are represented parametrically by
361  * a point, P(x0, y0, z0) and a direction normal, D = ax + by + cz.
362  * Any point on a line can be expressed by one variable 't', where
363  *
364  * X = a*t + x0, e.g., X = Dx*t + Px
365  * Y = b*t + y0,
366  * Z = c*t + z0.
367  *
368  * ^ ^ ^
369  * | | |
370  *
371  * W = D*t + P
372  *
373  * First, convert the line to the coordinate system of a "standard"
374  * heart. This is a heart which lies in the X-Y plane, circles the
375  * origin, and whose distance to cusps is one.
376  *
377  * Then find the equation of that line and the heart, which
378  * turns out (by substituting X, Y, Z above into the sextic equation above)
379  * to be a sextic equation S in 't' given below.
380  *
381  * S(t)=C6*t**6 + C5*t**5 + C4*t**4 + C3*t**3 + C2*t**2 + C1*t + C0 = 0.
382  *
383  * where C0, C1, C2, C3, C4, C5, C6 are coefficients of the equation.
384  *
385  * Solve the equation using a general polynomial root finder.
386  * Use those values of 't' to compute the points of intersection
387  * in the original coordinate system.
388  *
389  * Returns -
390  * 0 MISS
391  * >0 HIT
392  */
393 int
394 rt_hrt_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
395 {
396  register struct hrt_specific *hrt =
397  (struct hrt_specific *)stp->st_specific;
398  register struct seg *segp;
399  vect_t dprime; /* D' : The new shot direction */
400  vect_t pprime; /* P' : The new shot point */
401  vect_t trans; /* Translated shot vector */
402  vect_t norm_pprime; /* P' with normalized distance from heart */
403  bn_poly_t Xsqr, Ysqr, Zsqr; /* X^2, 9/4*Y^2, Z^2 - 1 */
404  bn_poly_t Acube, S; /* The sextic equation (of power 6) and A^3 */
405  bn_poly_t Zcube, A; /* Z^3 and X^2 + 9/4 * Y^2 + Z^2 - 1 */
406  bn_poly_t X2_Y2, Z3_X2_Y2; /* X^2 + 9/80*Y^2 and Z^3*(X^2 + 9/80*Y^2) */
407  bn_complex_t complex[6]; /* The complex roots */
408  double real[6]; /* The real roots */
409  int i;
410  int j;
411 
412  /* Translate the ray point */
413  (trans)[X] = (rp->r_pt)[X] - (hrt->hrt_V)[X];
414  (trans)[Y] = (rp->r_pt)[Y] - (hrt->hrt_V)[Y];
415  (trans)[Z] = (rp->r_pt)[Z] - (hrt->hrt_V)[Z];
416 
417  /* Scale and Rotate point to get P' */
418  pprime[X] = (hrt->hrt_SoR[0]*trans[X] + hrt->hrt_SoR[1]*trans[Y] + hrt->hrt_SoR[2]*trans[Z]) * 1.0/(hrt->hrt_SoR[15]);
419  pprime[Y] = (hrt->hrt_SoR[4]*trans[X] + hrt->hrt_SoR[5]*trans[Y] + hrt->hrt_SoR[6]*trans[Z]) * 1.0/(hrt->hrt_SoR[15]);
420  pprime[Z] = (hrt->hrt_SoR[8]*trans[X] + hrt->hrt_SoR[9]*trans[Y] + hrt->hrt_SoR[10]*trans[Z]) * 1.0/(hrt->hrt_SoR[15]);
421  /* Translate ray direction vector */
422  MAT4X3VEC(dprime, hrt->hrt_SoR, rp->r_dir);
423  VUNITIZE(dprime);
424 
425  /* Normalize distance from the heart. Substitutes a corrected ray
426  * point, which contains a translation along the ray direction to
427  * the closest approach to vertex of the heart.Translating the ray
428  * along the direction of the ray to the closest point near the
429  * heart's center vertex. Thus, the New ray origin is normalized.
430  */
431  VSCALE(norm_pprime, dprime, VDOT(pprime, dprime));
432  VSUB2(norm_pprime, pprime, norm_pprime);
433 
434  /**
435  * Generate the sextic equation S(t) = 0 to be passed through the root finder.
436  */
437  /* X**2 */
438  Xsqr.dgr = 2;
439  Xsqr.cf[0] = dprime[X] * dprime[X];
440  Xsqr.cf[1] = 2 * dprime[X] * pprime[X];
441  Xsqr.cf[2] = pprime[X] * pprime[X];
442 
443  /* 9/4 * Y**2*/
444  Ysqr.dgr = 2;
445  Ysqr.cf[0] = 9/4 * dprime[Y] * dprime[Y];
446  Ysqr.cf[1] = 9/2 * dprime[Y] * pprime[Y];
447  Ysqr.cf[2] = 9/4 * (pprime[Y] * pprime[Y]);
448 
449  /* Z**2 - 1 */
450  Zsqr.dgr = 2;
451  Zsqr.cf[0] = dprime[Z] * dprime[Z];
452  Zsqr.cf[1] = 2 * dprime[Z] * pprime[Z];
453  Zsqr.cf[2] = pprime[Z] * pprime[Z] - 1.0 ;
454 
455  /* A = X^2 + 9/4 * Y^2 + Z^2 - 1 */
456  A.dgr = 2;
457  A.cf[0] = Xsqr.cf[0] + Ysqr.cf[0] + Zsqr.cf[0];
458  A.cf[1] = Xsqr.cf[1] + Ysqr.cf[1] + Zsqr.cf[1];
459  A.cf[2] = Xsqr.cf[2] + Ysqr.cf[2] + Zsqr.cf[2];
460 
461  /* Z**3 */
462  Zcube.dgr = 3;
463  Zcube.cf[0] = dprime[Z] * Zsqr.cf[0];
464  Zcube.cf[1] = 1.5 * dprime[Z] * Zsqr.cf[1];
465  Zcube.cf[2] = 1.5 * pprime[Z] * Zsqr.cf[1];
466  Zcube.cf[3] = pprime[Z] * ( Zsqr.cf[2] + 1.0 );
467 
468  Acube.dgr = 6;
469  Acube.cf[0] = A.cf[0] * A.cf[0] * A.cf[0];
470  Acube.cf[1] = 3.0 * A.cf[0] * A.cf[0] * A.cf[1];
471  Acube.cf[2] = 3.0 * (A.cf[0] * A.cf[0] * A.cf[2] + A.cf[0] * A.cf[1] * A.cf[1]);
472  Acube.cf[3] = 6.0 * A.cf[0] * A.cf[1] * A.cf[2] + A.cf[1] * A.cf[1] * A.cf[1];
473  Acube.cf[4] = 3.0 * (A.cf[0] * A.cf[2] * A.cf[2] + A.cf[1] * A.cf[1] * A.cf[2]);
474  Acube.cf[5] = 3.0 * A.cf[1] * A.cf[2] * A.cf[2];
475  Acube.cf[6] = A.cf[2] * A.cf[2] * A.cf[2];
476 
477  /* X**2 + 9/80 Y**2 */
478  X2_Y2.dgr = 2;
479  X2_Y2.cf[0] = Xsqr.cf[0] + Ysqr.cf[0] / 20 ;
480  X2_Y2.cf[1] = Xsqr.cf[1] + Ysqr.cf[1] / 20 ;
481  X2_Y2.cf[2] = Xsqr.cf[2] + Ysqr.cf[2] / 20 ;
482 
483  /* Z**3 * (X**2 + 9/80 * Y**2) */
484  Z3_X2_Y2.dgr = 5;
485  Z3_X2_Y2.cf[0] = Zcube.cf[0] * X2_Y2.cf[0];
486  Z3_X2_Y2.cf[1] = X2_Y2.cf[0] * Zcube.cf[1];
487  Z3_X2_Y2.cf[2] = X2_Y2.cf[0] * Zcube.cf[2] + X2_Y2.cf[1] * Zcube.cf[0] + X2_Y2.cf[1] * Zcube.cf[1] + X2_Y2.cf[2] * Zcube.cf[0];
488  Z3_X2_Y2.cf[3] = X2_Y2.cf[0] * Zcube.cf[3] + X2_Y2.cf[1] * Zcube.cf[2] + X2_Y2.cf[2] * Zcube.cf[1];
489  Z3_X2_Y2.cf[4] = X2_Y2.cf[1] * Zcube.cf[3] + X2_Y2.cf[2] * Zcube.cf[2];
490  Z3_X2_Y2.cf[5] = X2_Y2.cf[2] * Zcube.cf[3];
491 
492  S.dgr = 6;
493  S.cf[0] = Acube.cf[0];
494  S.cf[1] = Acube.cf[1] - Z3_X2_Y2.cf[0];
495  S.cf[2] = Acube.cf[2] - Z3_X2_Y2.cf[1];
496  S.cf[3] = Acube.cf[3] - Z3_X2_Y2.cf[2];
497  S.cf[4] = Acube.cf[4] - Z3_X2_Y2.cf[3];
498  S.cf[5] = Acube.cf[5] - Z3_X2_Y2.cf[4];
499  S.cf[6] = Acube.cf[6] - Z3_X2_Y2.cf[5];
500 
501  /* It is known that the equation is sextic (of order 6). Therefore, if the
502  * root finder returns other than 6 roots, return an error.
503  */
504  if ((i = rt_poly_roots(&S, complex, stp->st_dp->d_namep)) != 6) {
505  if (i > 0) {
506  bu_log("hrt: rt_poly_roots() 6!=%d\n", i);
507  bn_pr_roots(stp->st_name, complex, i);
508  } else if (i < 0) {
509  static int reported = 0;
510  bu_log("The root solver failed to converge on a solution for %s\n", stp->st_dp->d_namep);
511  if (!reported) {
512  VPRINT("while shooting from:\t", rp->r_pt);
513  VPRINT("while shooting at:\t", rp->r_dir);
514  bu_log("Additional heart convergence failure details will be suppressed.\n");
515  reported = 1;
516  }
517  }
518  return 0; /* MISS */
519  }
520 
521  /* Only real roots indicate an intersection in real space.
522  *
523  * Look at each root returned; if the imaginary part is zero or
524  * sufficiently close, then use the real part as one value of 't'
525  * for the intersections
526  */
527  for (j = 0, i = 0; j < 6; j++) {
528  if (NEAR_ZERO(complex[j].im, ap->a_rt_i->rti_tol.dist))
529  real[i++] = complex[j].re;
530  }
531  /* Here, 'i' is number of points found */
532  switch (i) {
533  case 0:
534  return 0; /* No hit */
535 
536  default:
537  bu_log("rt_hrt_shot: reduced 6 to %d roots\n", i);
538  bn_pr_roots(stp->st_name, complex, 6);
539  return 0; /* No hit */
540 
541  case 2:
542  {
543  /* Sort most distant to least distant. */
544  fastf_t u;
545  if ((u=real[0]) < real[1]) {
546  /* bubble larger towards [0] */
547  real[0] = real[1];
548  real[1] = u;
549  }
550  }
551  break;
552  case 4:
553  {
554  short n;
555  short lim;
556 
557  /* Inline rt_pt_sort(). Sorts real[] into descending order. */
558  for (lim = i-1; lim > 0; lim--) {
559  for (n = 0; n < lim; n++) {
560  fastf_t u;
561  if ((u=real[n]) < real[n+1]) {
562  /* bubble larger towards [0] */
563  real[n] = real[n+1];
564  real[n+1] = u;
565  }
566  }
567  }
568  }
569  break;
570  case 6:
571  {
572  short num;
573  short limit;
574 
575  /* Inline rt_pt_sort(). Sorts real[] into descending order. */
576  for (limit = i-1; limit > 0; limit--) {
577  for (num = 0; num < limit; num++) {
578  fastf_t u;
579  if ((u=real[num]) < real[num+1]) {
580  /* bubble larger towards [0] */
581  real[num] = real[num+1];
582  real[num+1] = u;
583  }
584  }
585  }
586  }
587  break;
588  }
589 
590  /* Now, t[0] > t[npts-1] */
591  /* real[1] is entry point, and real[0] is farthest exit point */
592  RT_GET_SEG(segp, ap->a_resource);
593  segp->seg_stp = stp;
594  segp->seg_in.hit_dist = real[1];
595  segp->seg_out.hit_dist = real[0];
596  segp->seg_in.hit_surfno = segp->seg_out.hit_surfno = 0;
597  /* Set aside vector for rt_hrt_norm() later */
598  VJOIN1(segp->seg_in.hit_vpriv, pprime, real[1], dprime);
599  VJOIN1(segp->seg_out.hit_vpriv, pprime, real[0], dprime);
600  BU_LIST_INSERT(&(seghead->l), &(segp->l));
601 
602  if (i == 2)
603  return 2; /* HIT */
604 
605  /* 4 points */
606  /* real[3] is entry point, and real[2] is exit point */
607  RT_GET_SEG(segp, ap->a_resource);
608  segp->seg_stp = stp;
609  segp->seg_in.hit_dist = real[3];
610  segp->seg_out.hit_dist = real[2];
611  segp->seg_in.hit_surfno = segp->seg_out.hit_surfno = 0;
612  /* Set aside vector for rt_hrt_norm() later */
613  VJOIN1(segp->seg_in.hit_vpriv, pprime, real[3], dprime);
614  VJOIN1(segp->seg_out.hit_vpriv, pprime, real[2], dprime);
615  BU_LIST_INSERT(&(seghead->l), &(segp->l));
616 
617  if (i == 4)
618  return 4; /* HIT */
619 
620  /* 6 points */
621  /* real[5] is entry point, and real[4] is exit point */
622  RT_GET_SEG(segp, ap->a_resource);
623  segp->seg_stp = stp;
624  segp->seg_in.hit_dist = real[5];
625  segp->seg_out.hit_dist = real[4];
626  segp->seg_in.hit_surfno = segp->seg_out.hit_surfno = 0;
627  /* Set aside vector for rt_hrt_norm() later */
628  VJOIN1(segp->seg_in.hit_vpriv, pprime, real[5], dprime);
629  VJOIN1(segp->seg_out.hit_vpriv, pprime, real[4], dprime);
630  BU_LIST_INSERT(&(seghead->l), &(segp->l));
631  return 6;
632 }
633 
634 
635 #define RT_HRT_SEG_MISS(SEG) (SEG).seg_stp=(struct soltab *) 0;
636 /**
637  * This is the Becker vector version
638  */
639 void
640 rt_hrt_vshot(struct soltab **stp, struct xray **rp, struct seg *segp, int n, struct application *ap)
641 {
642  register struct hrt_specific *hrt;
643  vect_t dprime;
644  vect_t pprime;
645  vect_t work;
646  vect_t norm_pprime;
647  bn_poly_t Xsqr, Ysqr, Zsqr;
648  bn_poly_t A, Acube, Zcube;
649  bn_poly_t X2_Y2, Z3_X2_Y2;
650  bn_poly_t *S;
651  bn_complex_t (*complex)[6];
652  int num_roots, num_zero;
653  register int i;
654  fastf_t *cor_proj;
655 
656  /* Allocate space for polynomials and roots */
657  S = (bn_poly_t *)bu_malloc(n * sizeof(bn_poly_t) * 6, "hrt bn_complex_t");
658  complex = (bn_complex_t (*)[6])bu_malloc(n * sizeof(bn_complex_t) * 6, "hrt bn_complex_t");
659  cor_proj = (fastf_t *)bu_malloc(n * sizeof(fastf_t), "hrt_proj");
660 
661  /* Initialize seg_stp to assume hit (zero will then flag a miss) */
662  for (i = 0; i < n; i++ )
663  segp[i].seg_stp = stp[i];
664 
665  /* for each ray/heart pair */
666  for (i = 0; i < n; i++) {
667  if (segp[i].seg_stp == 0)
668  continue; /* Skip this iteration */
669 
670  hrt = (struct hrt_specific *)stp[i]->st_specific;
671 
672  /* Translate ray direction vector */
673  MAT4X3VEC(dprime, hrt->hrt_SoR, rp[i]->r_dir);
674  VUNITIZE(dprime);
675 
676  /* Use segp[i].seg_in.hit_normal as tmp to hold dprime */
677  VMOVE(segp[i].seg_in.hit_normal, dprime);
678  VSUB2(work, rp[i]->r_pt, hrt->hrt_V);
679  MAT4X3VEC(pprime, hrt->hrt_SoR, work);
680 
681  /* use segp[i].seg_out.hit_normal as tmp to hold pprime */
682  VMOVE(segp[i].seg_out.hit_normal, pprime);
683 
684  /* Normalize distance from the heart. Substitutes a corrected ray
685  * point, which contains a translation along the ray direction to
686  * the closest approach to vertex of the heart.Translating the ray
687  * along the direction of the ray to the closest point near the
688  * heart's center vertex. Thus, the New ray origin is normalized.
689  */
690  cor_proj[i] = VDOT(pprime, dprime);
691  VSCALE(norm_pprime, dprime, cor_proj[i]);
692  VSUB2(norm_pprime, pprime, norm_pprime);
693 
694  /*
695  * Generate sextic equation S(t) = 0 to be passed through the root finder
696  */
697  /* X**2 */
698  Xsqr.dgr = 2;
699  Xsqr.cf[0] = dprime[X] * dprime[X];
700  Xsqr.cf[1] = 2 * dprime[X] * pprime[X];
701  Xsqr.cf[2] = pprime[X] * pprime[X];
702 
703  /* 9/4 * Y**2*/
704  Ysqr.dgr = 2;
705  Ysqr.cf[0] = 9/4 * dprime[Y] * dprime[Y];
706  Ysqr.cf[1] = 9/2 * dprime[Y] * pprime[Y];
707  Ysqr.cf[2] = 9/4 * (pprime[Y] * pprime[Y]);
708 
709  /* Z**2 - 1 */
710  Zsqr.dgr = 2;
711  Zsqr.cf[0] = dprime[Z] * dprime[Z];
712  Zsqr.cf[1] = 2 * dprime[Z] * pprime[Z];
713  Zsqr.cf[2] = pprime[Z] * pprime[Z] - 1.0 ;
714 
715  /* A = X^2 + 9/4 * Y^2 + Z^2 - 1 */
716  A.dgr = 2;
717  A.cf[0] = Xsqr.cf[0] + Ysqr.cf[0] + Zsqr.cf[0];
718  A.cf[1] = Xsqr.cf[1] + Ysqr.cf[1] + Zsqr.cf[1];
719  A.cf[2] = Xsqr.cf[2] + Ysqr.cf[2] + Zsqr.cf[2];
720 
721  /* Z**3 */
722  Zcube.dgr = 3;
723  Zcube.cf[0] = dprime[Z] * Zsqr.cf[0];
724  Zcube.cf[1] = 1.5 * dprime[Z] * Zsqr.cf[1];
725  Zcube.cf[2] = 1.5 * pprime[Z] * Zsqr.cf[1];
726  Zcube.cf[3] = pprime[Z] * ( Zsqr.cf[2] + 1.0 );
727 
728  /* A**3 */
729  Acube.dgr = 6;
730  Acube.cf[0] = A.cf[0] * A.cf[0] * A.cf[0];
731  Acube.cf[1] = 3.0 * A.cf[0] * A.cf[0] * A.cf[1];
732  Acube.cf[2] = 3.0 * (A.cf[0] * A.cf[0] * A.cf[2] + A.cf[0] * A.cf[1] * A.cf[1]);
733  Acube.cf[3] = 6.0 * A.cf[0] * A.cf[1] * A.cf[2] + A.cf[1] * A.cf[1] * A.cf[1];
734  Acube.cf[4] = 3.0 * (A.cf[0] * A.cf[2] * A.cf[2] + A.cf[1] * A.cf[1] * A.cf[2]);
735  Acube.cf[5] = 3.0 * A.cf[1] * A.cf[2] * A.cf[2];
736  Acube.cf[6] = A.cf[2] * A.cf[2] * A.cf[2];
737 
738  /* X**2 + 9/80 Y**2 */
739  X2_Y2.dgr = 2;
740  X2_Y2.cf[0] = Xsqr.cf[0] + Ysqr.cf[0] / 20 ;
741  X2_Y2.cf[1] = Xsqr.cf[1] + Ysqr.cf[1] / 20 ;
742  X2_Y2.cf[2] = Xsqr.cf[2] + Ysqr.cf[2] / 20 ;
743 
744  /* Z**3 * (X**2 + 9/80 * Y**2) */
745  Z3_X2_Y2.dgr = 5;
746  Z3_X2_Y2.cf[0] = Zcube.cf[0] * X2_Y2.cf[0];
747  Z3_X2_Y2.cf[1] = X2_Y2.cf[0] * Zcube.cf[1];
748  Z3_X2_Y2.cf[2] = X2_Y2.cf[0] * Zcube.cf[2] + X2_Y2.cf[1] * Zcube.cf[0] + X2_Y2.cf[1] * Zcube.cf[1] + X2_Y2.cf[2] * Zcube.cf[0];
749  Z3_X2_Y2.cf[3] = X2_Y2.cf[0] * Zcube.cf[3] + X2_Y2.cf[1] * Zcube.cf[2] + X2_Y2.cf[2] * Zcube.cf[1];
750  Z3_X2_Y2.cf[4] = X2_Y2.cf[1] * Zcube.cf[3] + X2_Y2.cf[2] * Zcube.cf[2];
751  Z3_X2_Y2.cf[5] = X2_Y2.cf[2] * Zcube.cf[3];
752 
753  /* S(t) = 0 */
754  S[i].dgr = 6;
755  S[i].cf[0] = Acube.cf[0];
756  S[i].cf[1] = Acube.cf[1] - Z3_X2_Y2.cf[0];
757  S[i].cf[2] = Acube.cf[2] - Z3_X2_Y2.cf[1];
758  S[i].cf[3] = Acube.cf[3] - Z3_X2_Y2.cf[2];
759  S[i].cf[4] = Acube.cf[4] - Z3_X2_Y2.cf[3];
760  S[i].cf[5] = Acube.cf[5] - Z3_X2_Y2.cf[4];
761  S[i].cf[6] = Acube.cf[6] - Z3_X2_Y2.cf[5];
762 
763  }
764  for (i = 0; i < n; i++) {
765  continue; /* Skip this iteration */
766  /*
767  * It is known that the equation is sextic (of order 6).
768  * Therefore, if the root finder returns other than six roots, Error
769  */
770  if ((num_roots = rt_poly_roots(&(S[i]), &(complex[i][0]), (*stp)->st_dp->d_namep)) != 6) {
771  if (num_roots > 0) {
772  bu_log("hrt: rt_poly_roots() 6 != %d\n", num_roots);
773  bn_pr_roots("hrt", complex[i], num_roots);
774  } else if (num_roots < 0) {
775  static int reported = 0;
776  bu_log("The root solver failed to converge on a solution for %s\n", stp[i]->st_dp->d_namep);
777  if (!reported) {
778  VPRINT("while shooting from: \t", rp[i]->r_pt);
779  VPRINT("while shooting at:\t", rp[i]->r_dir);
780  bu_log("Additional heart convergence failure details will be suppressed.\n");
781  reported = 1;
782  }
783  }
784  RT_HRT_SEG_MISS(segp[i]);
785  }
786 
787  /* for each ray/heart pair */
788  for (i = 0; i < n; i++) {
789  if (segp[i].seg_stp == 0)
790  continue; /* Skip current iteration */
791 
792  /* Only real roots indicate an intersection in real space.
793  *
794  * Look at each root returned; if the imaginary part is zero or
795  * sufficiently close, then use the real part as one value of 't'
796  * for the intersections
797  *
798  * Reverse translation by adding distance to all 'real' values
799  * Reuse S to hold 'real' values
800  */
801  num_zero = 0;
802  if (NEAR_ZERO(complex[i][0].im, ap->a_rt_i->rti_tol.dist)) {
803  S[i].cf[num_zero++] = complex[i][0].re - cor_proj[i];
804  }
805  if (NEAR_ZERO(complex[i][1].im, ap->a_rt_i->rti_tol.dist)) {
806  S[i].cf[num_zero++] = complex[i][1].re - cor_proj[i];
807  }
808  if (NEAR_ZERO(complex[i][2].im, ap->a_rt_i->rti_tol.dist)) {
809  S[i].cf[num_zero++] = complex[i][2].re - cor_proj[i];
810  }
811  if (NEAR_ZERO(complex[i][3].im, ap->a_rt_i->rti_tol.dist)) {
812  S[i].cf[num_zero++] = complex[i][3].re - cor_proj[i];
813  }
814  if (NEAR_ZERO(complex[i][4].im, ap->a_rt_i->rti_tol.dist)) {
815  S[i].cf[num_zero++] = complex[i][4].re - cor_proj[i];
816  }
817  if (NEAR_ZERO(complex[i][5].im, ap->a_rt_i->rti_tol.dist)) {
818  S[i].cf[num_zero++] = complex[i][5].re - cor_proj[i];
819  }
820  S[i].dgr = num_zero;
821 
822  /* Here 'i' is the number of points found */
823  if (num_zero == 0) {
824  RT_HRT_SEG_MISS(segp[i]); /* MISS */
825  } else if ((num_zero != 2 && num_zero != 4 ) && (num_zero != 6)) {
826  RT_HRT_SEG_MISS(segp[i]); /* MISS */
827  }
828  }
829 
830  /* Process each one segment hit */
831  for (i = 0; i < n; i++) {
832  if (segp[i].seg_stp == 0)
833  continue; /* Skip This Iteration */
834  if (S[i].dgr != 2)
835  continue; /* Not one segment */
836 
837  hrt = (struct hrt_specific *)stp[i]->st_specific;
838 
839  /* segp[i].seg_in.hit_normal holds dprime */
840  /* segp[i].seg_out.hit_normal holds pprime */
841  if (S[i].cf[1] < S[i].cf[0]) {
842  /* S[i].cf[1] is entry point */
843  segp[i].seg_in.hit_dist = S[i].cf[1];
844  segp[i].seg_in.hit_dist = S[i].cf[0];
845  /* Set aside vector for rt_hrt_norm() later */
846  VJOIN1(segp[i].seg_in.hit_vpriv, segp[i].seg_out.hit_normal, S[i].cf[1], segp[i].seg_in.hit_normal);
847  VJOIN1(segp[i].seg_in.hit_vpriv, segp[i].seg_out.hit_normal, S[i].cf[0], segp[i].seg_in.hit_normal);
848  } else {
849  /* S[i].cf[0] is entry point */
850  segp[i].seg_in.hit_dist = S[i].cf[0];
851  segp[i].seg_in.hit_dist = S[i].cf[1];
852  /* Set aside vector for rt_hrt_norm() later */
853  VJOIN1(segp[i].seg_in.hit_vpriv, segp[i].seg_out.hit_normal, S[i].cf[0], segp[i].seg_in.hit_normal);
854  VJOIN1(segp[i].seg_in.hit_vpriv, segp[i].seg_out.hit_normal, S[i].cf[1], segp[i].seg_in.hit_normal);
855  }
856  }
857 
858  /* Process each two segment hit */
859  for (i = 0; i < n; i++) {
860  if (segp[i].seg_stp == 0)
861  continue; /* Skip This Iteration */
862  if (S[i].dgr != 4)
863  continue; /* Not one segment */
864 
865  hrt = (struct hrt_specific *)stp[i]->st_specific;
866 
867  /* Sort most distant to least distant */
868  rt_pt_sort(S[i].cf, 4);
869  /* Now, t[0] > t[npts - 1] */
870  /* segp[i].seg_in.hit_normal holds dprime */
871  VMOVE(dprime, segp[i].seg_out.hit_normal);
872  /* segp[i].seg_out.hit_normal holds pprime */
873  VMOVE(pprime, segp[i].seg_out.hit_normal);
874 
875  /* S[i].cf[1] is entry point */
876  segp[i].seg_in.hit_dist = S[i].cf[3];
877  segp[i].seg_in.hit_dist = S[i].cf[2];
878  /* Set aside vector for rt_hrt_norm() later */
879  VJOIN1(segp[i].seg_in.hit_vpriv, segp[i].seg_out.hit_normal, S[i].cf[1], segp[i].seg_in.hit_normal);
880  VJOIN1(segp[i].seg_in.hit_vpriv, segp[i].seg_out.hit_normal, S[i].cf[0], segp[i].seg_in.hit_normal);
881  }
882 
883  /* Process each three segment hit */
884  for (i = 0; i < n; i++) {
885  if (segp[i].seg_stp == 0)
886  continue; /* Skip This Iteration */
887  if (S[i].dgr != 6)
888  continue; /* Not one segment */
889 
890  hrt = (struct hrt_specific *)stp[i]->st_specific;
891 
892  /* Sort most distant to least distant */
893  rt_pt_sort(S[i].cf, 6);
894  /* Now, t[0] > t[npts - 1] */
895  /* segp[i].seg_in.hit_normal holds dprime */
896  VMOVE(dprime, segp[i].seg_out.hit_normal);
897  /* segp[i].seg_out.hit_normal holds pprime */
898  VMOVE(pprime, segp[i].seg_out.hit_normal);
899 
900  /* S[i].cf[1] is entry point */
901  segp[i].seg_in.hit_dist = S[i].cf[5];
902  segp[i].seg_in.hit_dist = S[i].cf[4];
903  /* Set aside vector for rt_hrt_norm() later */
904  VJOIN1(segp[i].seg_in.hit_vpriv, pprime, S[i].cf[5], dprime);
905  VJOIN1(segp[i].seg_in.hit_vpriv, pprime, S[i].cf[4], dprime);
906  }
907 
908  }
909 
910  /* Free tmp space used */
911  bu_free((char *)S, "hrt S");
912  bu_free((char *)complex, "hrt complex");
913  bu_free((char *)cor_proj, "hrt_cor_proj");
914 
915 }
916 
917 
918 /**
919  * Compute the normal to the heart, given a point on the heart
920  * centered at the origin on the X-Y plane. The gradient of the heart
921  * at that point is in fact the normal vector, which will have to be
922  * given unit length. To make this useful for the original heart, it
923  * will have to be rotated back to the orientation of the original
924  * heart.
925  *
926  * Given that the equation for the heart is:
927  *
928  * [ X**2 + 9/4 * Y**2 + Z**2 - 1 ]**3 - Z**3 * (X**2 + 9/80*Y**2) = 0.
929  *
930  * let w = X**2 + 9/4 * Y**2 + Z**2 - 1, then the equation becomes:
931  *
932  * w**3 - Z**3 * (X**2 + 9/80 * Y**2) = 0.
933  *
934  * For f(x, y, z) = 0, the gradient of f() is (df/dx, df/dy, df/dz).
935  *
936  * df/dx = 6 * X * (w**2 - Z**3/3)
937  * df/dy = 6 * Y * (12/27 * w**2 - 80/3 * Z**3)
938  * df/dz = 6 * Z * ( w**2 - 1/2 * Z * (X**2 + 9/80 * Y**2))
939  *
940  * Since we rescale the gradient (normal) to unity, we divide the
941  * above equations by six here.
942  */
943 void
944 rt_hrt_norm(register struct hit *hitp, register struct xray *rp)
945 {
946 
947  fastf_t w, fx, fy, fz;
948  vect_t work;
949 
950  VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
951  w = hitp->hit_vpriv[X] * hitp->hit_vpriv[X]
952  + 9.0/4.0 * hitp->hit_vpriv[Y] * hitp->hit_vpriv[Y]
953  + hitp->hit_vpriv[Z] * hitp->hit_vpriv[Z] - 1.0;
954  fx = (w * w - 1/3 * hitp->hit_vpriv[Z] * hitp->hit_vpriv[Z] * hitp->hit_vpriv[Z]) * hitp->hit_vpriv[X];
955  fy = hitp->hit_vpriv[Y] * (12/27 * w * w - 80/3 * hitp->hit_vpriv[Z] * hitp->hit_vpriv[Z] * hitp->hit_vpriv[Z]);
956  fz = (w * w - 0.5 * hitp->hit_vpriv[Z] * (hitp->hit_vpriv[X] * hitp->hit_vpriv[X] + 9/80 * hitp->hit_vpriv[Y] * hitp->hit_vpriv[Y])) * hitp->hit_vpriv[Z];
957  VSET(work, fx, fy, fz);
958  VMOVE(hitp->hit_normal, work);
959 }
960 
961 
962 /**
963  * Return the curvature of the heart.
964  */
965 void
967 {
968  bu_log("rt_hrt_curve: Not implemented yet!\n");
969 }
970 
971 
972 void
974 {
975  bu_log("rt_hrt_uv: Not implemented yet!\n");
976 }
977 
978 
979 void
980 rt_hrt_free(struct soltab *stp)
981 {
982  struct hrt_specific *hrt =
983  (struct hrt_specific *)stp->st_specific;
984 
985  BU_PUT(hrt, struct hrt_specific);
986 }
987 
988 
989 /**
990  * Similar to code used by the ELL
991  */
992 
993 #define HRTOUT(n) ov+(n-1)*3
994 void
996 {
997  static fastf_t c, d, e, f, g, h, i, j, k, l;
998 
999  c = k = 0.95105; /* cos(18.0 degrees) */
1000  d = j = 0.80901; /* cos(36.0 degrees) */
1001  e = i = 0.58778; /* cos(54.0 degrees) */
1002  f = h = 0.30901; /* cos(72.0 degrees) */
1003  g = 0.00000; /* cos(90.0 degrees) */
1004  l = 1.00000; /* sin(90.0 degrees) */
1005 
1006  /*
1007  * For angle theta, compute surface point as
1008  *
1009  * V + cos(theta) * A + sin(theta) * B
1010  *
1011  * note that sin(theta) is cos(90-theta); arguments in degrees.
1012  */
1013 
1014  VADD2(HRTOUT(1), V, A);
1015  VJOIN2(HRTOUT(2), V, c, A, h, B);
1016  VJOIN2(HRTOUT(3), V, d, A, i, B);
1017  VJOIN2(HRTOUT(4), V, e, A, j, B);
1018  VJOIN2(HRTOUT(5), V, f, A, k, B);
1019  VJOIN2(HRTOUT(6), V, g, A, l, B);
1020  VADD2(HRTOUT(7), V, B);
1021  VJOIN2(HRTOUT(8), V, -g, A, l, B);
1022  VJOIN2(HRTOUT(9), V, -f, A, k, B);
1023  VJOIN2(HRTOUT(10), V, -e, A, j, B);
1024  VJOIN2(HRTOUT(11), V, -d, A, i, B);
1025  VJOIN2(HRTOUT(12), V, -c, A, h, B);
1026  VSUB2(HRTOUT(13), V, A);
1027  VJOIN2(HRTOUT(14), V, -c, A, -h, B);
1028  VJOIN2(HRTOUT(15), V, -d, A, -i, B);
1029  VJOIN2(HRTOUT(16), V, -e, A, -j, B);
1030  VJOIN2(HRTOUT(17), V, -f, A, -k, B);
1031  VJOIN2(HRTOUT(18), V, -g, A, -l, B);
1032  VSUB2(HRTOUT(19), V, B);
1033  VJOIN2(HRTOUT(20), V, g, A, -l, B);
1034  VJOIN2(HRTOUT(21), V, f, A, -k, B);
1035  VJOIN2(HRTOUT(22), V, e, A, -j, B);
1036  VJOIN2(HRTOUT(23), V, d, A, -i, B);
1037  VJOIN2(HRTOUT(24), V, c, A, -h, B);
1038 }
1039 
1040 
1041 int
1043 {
1044  bu_log("rt_adaptive_plot: Not implemented yet!\n");
1045  return 0;
1046 }
1047 
1048 
1049 int
1050 rt_hrt_plot(struct bu_list *vhead, struct rt_db_internal *ip,const struct rt_tess_tol *ttol, const struct bn_tol *UNUSED(tol), const struct rt_view_info *UNUSED(info))
1051 {
1052  fastf_t c, dtol, mag_h, ntol = M_PI, r1, r2, **ellipses, theta_prev, theta_new;
1053  int *pts_dbl;
1054  int nseg; /* The number of line segments in a particular ellipse */
1055  int j, k, jj, na, nb;
1056  int nell; /* The number of ellipses in the parabola */
1057  int recalc_b, i, ellipse_below, ellipse_above;
1058  mat_t R;
1059  mat_t invR;
1060  point_t p1;
1061  struct rt_pt_node *pos_a, *pos_b, *pts_a, *pts_b;
1062  vect_t A, Au, B, Bu, Cu;
1063  vect_t V, Work;
1064  vect_t lower_cusp, upper_cusp, upper_cusp_xdir;
1065  vect_t highest_point_left, highest_point_right, top_xdir, top_ydir;
1066  vect_t top01_center, top01_xdir, top01_ydir;
1067  vect_t top1_center, top1_xdir, top1_ydir;
1068  vect_t top02_center, top02_xdir, top02_ydir;
1069  vect_t v1_left, xdir1_left, ydir1_left, top1_center_left, top1_center_right, v1_right, xdir1_right, ydir1_right;
1070  vect_t v2_left, xdir2_left, ydir2_left, top2_center_left, top2_center_right, v2_right, xdir2_right, ydir2_right;
1071  vect_t v3_left, xdir3_left, ydir3_left, top3_center_left, top3_center_right, v3_right, xdir3_right, ydir3_right;
1072  vect_t v4_left, xdir4_left, ydir4_left, top4_center_left, top4_center_right, v4_right, xdir4_right, ydir4_right;
1073 
1074  struct rt_hrt_internal *hip;
1075 
1076  fastf_t top[24*3];
1077  fastf_t top01[24*3];
1078  fastf_t top02[24*3];
1079  fastf_t top1[24*3];
1080  fastf_t top1_left_lobe[24*3];
1081  fastf_t top1_right_lobe[24*3];
1082  fastf_t top2_left_lobe[24*3];
1083  fastf_t top2_right_lobe[24*3];
1084  fastf_t top3_left_lobe[24*3];
1085  fastf_t top3_right_lobe[24*3];
1086  fastf_t top4_left_lobe[24*3];
1087  fastf_t top4_right_lobe[24*3];
1088  fastf_t top5_upper_cusp[24*3];
1089 
1090  fastf_t angle_c = 0.64656; /* sin(40.2) */
1091  fastf_t angle_d = 0.80901; /* sin(54.0) */
1092  fastf_t angle_e = 0.90350; /* sin(64.5) */
1093  fastf_t angle_f = 0.95000; /* sin(71.8) */
1094  fastf_t angle_g = 0.33990; /* sin(19.8) */
1095  fastf_t angle_h = 0.48357; /* sin(29.0) */
1096 
1097  BU_CK_LIST_HEAD(vhead);
1098  RT_CK_DB_INTERNAL(ip);
1099  hip = (struct rt_hrt_internal *)ip->idb_ptr;
1100  RT_HRT_CK_MAGIC(hip);
1101 
1102  VSET(top01_center, 0 , 0 , hip->zdir[Z] * angle_g );
1103  VSET(top02_center, 0 , 0 , hip->zdir[Z] * angle_h );
1104  VSET(top1_center_left, -0.77 * hip->xdir[X], 0, hip->zdir[Z]);
1105  VSET(top1_center_right, hip->xdir[X] * 0.77, 0, hip->zdir[Z]);
1106  VSET(top1_center, 0, 0, hip->zdir[Z] * 0.64);
1107  VSET(top2_center_left, -1.00 * hip->xdir[X] * 0.63, 0, hip->zdir[Z]);
1108  VSET(top2_center_right, hip->xdir[X] * 0.63, 0, hip->zdir[Z]);
1109  VSET(top3_center_left, -1.00 * hip->xdir[X] * 0.55, 0, hip->zdir[Z]);
1110  VSET(top3_center_right, hip->xdir[X] * 0.55, 0, hip->zdir[Z]);
1111  VSET(top4_center_left, -1.00 * hip->xdir[X] * 0.38, 0, hip->zdir[Z] * 0.77 );
1112  VSET(top4_center_right, hip->xdir[X] * 0.38, 0, hip->zdir[Z] * 0.77 );
1113 
1114  VSCALE(v1_left, top1_center_left, angle_c);
1115  VSCALE(v1_right, top1_center_right, angle_c);
1116  VSCALE(v2_left, top2_center_left, angle_d );
1117  VSCALE(v2_right, top2_center_right, angle_d );
1118  VSCALE(v3_left, top3_center_left, angle_e );
1119  VSCALE(v3_right, top3_center_right, angle_e );
1120  VSCALE(v4_left, top4_center_left, angle_f );
1121  VSCALE(v4_right, top4_center_right, angle_f);
1122 
1123  VSET(upper_cusp, (top3_center_left[X] + top3_center_right[X]) * 0.70,
1124  (top3_center_left[Y] + top3_center_right[Y]) * 0.50,
1125  hip->zdir[Z] * 0.64);
1126  VSET(highest_point_left, hip->xdir[X] * -0.50 , 0 , hip->zdir[Z] );
1127  VSET(highest_point_right, hip->xdir[X] * 0.50 , 0 , hip->zdir[Z] );
1128 
1129  VSCALE(top_xdir, hip->xdir, 0.80);
1130  VSET(top01_xdir, hip->xdir[X] * 0.99, 0, 0);
1131  VSET(top02_xdir, hip->xdir[X] * 1.05, 0, 0);
1132  VSET(top1_xdir, hip->xdir[X] * 1.01, 0, 0);
1133  VSET(xdir1_left, v4_left[Z] * 0.68, 0, 0 );
1134  VSET(xdir1_right, v4_right[Z] * 0.68, 0, 0 );
1135  VSET(xdir2_left, v3_left[Z] * 0.45, 0, 0 );
1136  VSET(xdir2_right, v3_right[Z] * 0.45, 0, 0 );
1137  VSET(xdir3_left, v2_left[Z] * 0.40, 0, 0 );
1138  VSET(xdir3_right, v2_right[Z] * 0.40, 0, 0 );
1139  VSET(xdir4_left, v1_left[Z] * 0.01, 0, 0 );
1140  VSET(xdir4_right, v1_right[Z] * 0.01, 0, 0 );
1141 
1142  VSCALE(top_ydir, hip->ydir, 0.50);
1143  VSET(top01_ydir, 0, hip->ydir[Y] * 0.61, 0 );
1144  VSET(top02_ydir, 0, hip->ydir[Y] * 0.61, 0 );
1145  VSET(top1_ydir, 0, hip->ydir[Y] * 0.53, 0 );
1146  VSET(ydir1_left, 0, v4_left[Z] * 0.71, 0 );
1147  VSET(ydir1_right, 0, v4_right[Z] * 0.71, 0 );
1148  VSET(ydir2_left, 0, v3_left[Z] * 0.45, 0 );
1149  VSET(ydir2_right, 0, v3_right[Z] * 0.45, 0 );
1150  VSET(ydir3_left, 0, v2_left[Z] * 0.27, 0 );
1151  VSET(ydir3_right, 0, v2_right[Z] * 0.27, 0 );
1152  VSET(ydir4_left, 0, v1_left[Z] * 0.01, 0 );
1153  VSET(ydir4_right, 0, v1_right[Z] * 0.01, 0 );
1154  VSET(upper_cusp_xdir, 0 , v3_left[Z] * 0.01 , 0 );
1155 
1156  rt_hrt_24pts(top, hip->v, top_xdir, top_ydir);
1157  rt_hrt_24pts(top01, top01_center, top01_xdir, top01_ydir );
1158  rt_hrt_24pts(top02, top02_center, top02_xdir, top02_ydir );
1159  rt_hrt_24pts(top1, top1_center, top1_xdir, top1_ydir );
1160  rt_hrt_24pts(top1_left_lobe, v1_left, xdir1_left, ydir1_left);
1161  rt_hrt_24pts(top1_right_lobe, v1_right, xdir1_right, ydir1_right);
1162  rt_hrt_24pts(top2_left_lobe, v2_left, xdir2_left, ydir2_left);
1163  rt_hrt_24pts(top2_right_lobe, v2_right, xdir2_right, ydir2_right);
1164  rt_hrt_24pts(top3_left_lobe, v3_left, xdir3_left, ydir3_left);
1165  rt_hrt_24pts(top3_right_lobe, v3_right, xdir3_right, ydir3_right);
1166  rt_hrt_24pts(top4_left_lobe, highest_point_left, xdir4_left, ydir4_left);
1167  rt_hrt_24pts(top4_right_lobe, highest_point_right, xdir4_right, ydir4_right);
1168  rt_hrt_24pts(top5_upper_cusp, upper_cusp, upper_cusp_xdir, upper_cusp_xdir);
1169 
1170  RT_ADD_VLIST(vhead, &top[23*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1171  for (i = 0; i < 24; i++) {
1172  RT_ADD_VLIST(vhead, &top[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1173  }
1174 
1175  RT_ADD_VLIST(vhead, &top01[23*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1176  for (i = 0; i < 24; i++) {
1177  RT_ADD_VLIST(vhead, &top01[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1178  }
1179 
1180  RT_ADD_VLIST(vhead, &top02[23*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1181  for (i = 0; i < 24; i++) {
1182  RT_ADD_VLIST(vhead, &top02[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1183  }
1184 
1185  RT_ADD_VLIST(vhead, &top1[23*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1186  for (i = 0; i < 24; i++) {
1187  RT_ADD_VLIST(vhead, &top1[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1188  }
1189 
1190  RT_ADD_VLIST(vhead, &top1_left_lobe[23*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1191  for (i = 0; i < 24; i++) {
1192  RT_ADD_VLIST(vhead, &top1_left_lobe[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1193  }
1194 
1195  RT_ADD_VLIST(vhead, &top1_right_lobe[23*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1196  for (i = 0; i < 24; i++) {
1197  RT_ADD_VLIST(vhead, &top1_right_lobe[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1198  }
1199 
1200  RT_ADD_VLIST(vhead, &top2_left_lobe[23*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1201  for (i = 0; i < 24; i++) {
1202  RT_ADD_VLIST(vhead, &top2_left_lobe[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1203  }
1204 
1205  RT_ADD_VLIST(vhead, &top2_right_lobe[23*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1206  for (i = 0; i < 24; i++) {
1207  RT_ADD_VLIST(vhead, &top2_right_lobe[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1208  }
1209 
1210  RT_ADD_VLIST(vhead, &top3_left_lobe[23*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1211  for (i = 0; i < 24; i++) {
1212  RT_ADD_VLIST(vhead, &top3_left_lobe[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1213  }
1214 
1215  RT_ADD_VLIST(vhead, &top3_right_lobe[23*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1216  for (i = 0; i < 24; i++) {
1217  RT_ADD_VLIST(vhead, &top3_right_lobe[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1218  }
1219 
1220  RT_ADD_VLIST(vhead, &top4_left_lobe[23*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1221  for (i = 0; i < 24; i++) {
1222  RT_ADD_VLIST(vhead, &top4_left_lobe[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1223  }
1224 
1225  RT_ADD_VLIST(vhead, &top4_right_lobe[23*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1226  for (i = 0; i < 24; i++) {
1227  RT_ADD_VLIST(vhead, &top4_right_lobe[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1228  }
1229 
1230  RT_ADD_VLIST(vhead, &top5_upper_cusp[23*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1231  for (i = 0; i < 24; i++) {
1232  RT_ADD_VLIST(vhead, &top5_upper_cusp[i*ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1233  }
1234 
1235  /*
1236  Make connections between ellipses
1237  This for repetition statement draws connectors between
1238  the nell-1 ellipses drawn above
1239  */
1240 
1241  for (k = 1; k < 24; k++) {
1242  RT_ADD_VLIST(vhead,
1243  &top01[k*ELEMENTS_PER_VECT],
1245  RT_ADD_VLIST(vhead,
1246  &top02[k*ELEMENTS_PER_VECT],
1248  RT_ADD_VLIST(vhead,
1249  &top02[k*ELEMENTS_PER_VECT],
1251  RT_ADD_VLIST(vhead,
1252  &top1[k*ELEMENTS_PER_VECT],
1254  RT_ADD_VLIST(vhead,
1255  &top1_left_lobe[k*ELEMENTS_PER_VECT],
1257  RT_ADD_VLIST(vhead,
1258  &top2_left_lobe[k*ELEMENTS_PER_VECT],
1260  RT_ADD_VLIST(vhead,
1261  &top1_right_lobe[k*ELEMENTS_PER_VECT],
1263  RT_ADD_VLIST(vhead,
1264  &top2_right_lobe[k*ELEMENTS_PER_VECT],
1266  RT_ADD_VLIST(vhead,
1267  &top2_left_lobe[k*ELEMENTS_PER_VECT],
1269  RT_ADD_VLIST(vhead,
1270  &top3_left_lobe[k*ELEMENTS_PER_VECT],
1272  RT_ADD_VLIST(vhead,
1273  &top2_right_lobe[k*ELEMENTS_PER_VECT],
1275  RT_ADD_VLIST(vhead,
1276  &top3_right_lobe[k*ELEMENTS_PER_VECT],
1278  RT_ADD_VLIST(vhead,
1279  &top3_left_lobe[k*ELEMENTS_PER_VECT],
1281  RT_ADD_VLIST(vhead,
1282  &top4_left_lobe[k*ELEMENTS_PER_VECT],
1284  RT_ADD_VLIST(vhead,
1285  &top3_right_lobe[k*ELEMENTS_PER_VECT],
1287  RT_ADD_VLIST(vhead,
1288  &top4_right_lobe[k*ELEMENTS_PER_VECT],
1290  }
1291 
1292  mag_h = MAGNITUDE(hip->zdir);
1293  r1 = MAGNITUDE(hip->xdir) * 0.80;
1294  r2 = MAGNITUDE(hip->ydir) * 0.61;
1295  c = hip->d;
1296 
1297  /* make unit vectors in A, B, and AxB directions */
1298  VMOVE(Cu, hip->zdir);
1299  VUNITIZE(Cu);
1300  VMOVE(Au, hip->xdir);
1301  VUNITIZE(Au);
1302  VMOVE(Bu, hip->ydir);
1303  VUNITIZE(Bu);
1304 
1305  /* Compute R and Rinv matrices */
1306  MAT_IDN(R);
1307  VREVERSE(&R[0], Bu);
1308  VMOVE(&R[4], Au);
1309  VREVERSE(&R[8], Cu);
1310  bn_mat_trn(invR, R); /* inv of rot mat is trn */
1311 
1312  dtol = primitive_get_absolute_tolerance(ttol, r2 * 2.00);
1313 
1314  /*
1315  * build ehy from 2 hyperbolas
1316  */
1317 
1318  /* approximate positive half of hyperbola along semi-minor axis */
1319  BU_ALLOC(pts_b, struct rt_pt_node);
1320  BU_ALLOC(pts_b->next, struct rt_pt_node);
1321 
1322  pts_b->next->next = NULL;
1323  VSET(pts_b->p, 0, 0, 0);
1324  VSET(pts_b->next->p, 0, 0, mag_h/3);
1325  /* 2 endpoints in 1st approximation */
1326  nb = 2;
1327  /* recursively break segment 'til within error tolerances */
1328  nb += rt_mk_hyperbola(pts_b, mag_h/3, mag_h, c, dtol, ntol);
1329  nell = nb - 1; /* Number of ellipses needed */
1330 
1331  /*
1332  * construct positive half of hyperbola along semi-major axis of
1333  * ehy using same z coords as hyperbola along semi-minor axis
1334  */
1335 
1336  BU_ALLOC(pts_a, struct rt_pt_node);
1337  VMOVE(pts_a->p, pts_b->p); /* 1st pt is the apex */
1338  pts_a->next = NULL;
1339  pos_b = pts_b->next;
1340  pos_a = pts_a;
1341  while (pos_b) {
1342  /* copy node from b_hyperbola to a_hyperbola */
1343  BU_ALLOC(pos_a->next, struct rt_pt_node);
1344  pos_a = pos_a->next;
1345  pos_a->p[Z] = pos_b->p[Z];
1346  /* at given z, find y on hyperbola */
1347  pos_a->p[Y] = r1
1348  * sqrt(((pos_b->p[Z] + mag_h + c)
1349  * (pos_b->p[Z] + mag_h + c) - c*c)
1350  / (mag_h*(mag_h + 2*c)));
1351  pos_a->p[X] = 0;
1352  pos_b = pos_b->next;
1353  }
1354  pos_a->next = NULL;
1355 
1356  /* see if any segments need further breaking up */
1357  recalc_b = 0;
1358  pos_a = pts_a;
1359  while (pos_a->next) {
1360  na = rt_mk_hyperbola(pos_a, r1, mag_h, c, dtol, ntol);
1361  if (na != 0) {
1362  recalc_b = 1;
1363  nell += na;
1364  }
1365  pos_a = pos_a->next;
1366  }
1367  /* if any segments were broken, recalculate hyperbola 'a' */
1368  if (recalc_b) {
1369  /* free mem for old approximation of hyperbola */
1370  pos_b = pts_b;
1371  while (pos_b) {
1372  struct rt_pt_node *next;
1373 
1374  /* get next node before freeing */
1375  next = pos_b->next;
1376  bu_free((char *)pos_b, "rt_pt_node");
1377  pos_b = next;
1378  }
1379 
1380  /* construct hyperbola along semi-major axis of ehy using same
1381  * z coords as parabola along semi-minor axis
1382  */
1383  BU_ALLOC(pts_b, struct rt_pt_node);
1384  pts_b->p[Z] = pts_a->p[Z];
1385  pts_b->next = NULL;
1386  pos_a = pts_a->next;
1387  pos_b = pts_b;
1388  while (pos_a) {
1389  /* copy node from a_hyperbola to b_hyperbola */
1390  BU_ALLOC(pos_b->next, struct rt_pt_node);
1391  pos_b = pos_b->next;
1392  pos_b->p[Z] = pos_a->p[Z];
1393  /* at given z, find y on hyperbola */
1394  pos_b->p[Y] = r2 * sqrt(((pos_a->p[Z] + mag_h + c)
1395  * (pos_a->p[Z] + mag_h + c) - c*c)
1396  / (mag_h*(mag_h + 2*c)));
1397  pos_b->p[X] = 0;
1398  pos_a = pos_a->next;
1399  }
1400  pos_b->next = NULL;
1401  }
1402 
1403  /* make array of ptrs to ehy ellipses */
1404  ellipses = (fastf_t **)bu_malloc(nell * sizeof(fastf_t *), "fastf_t ell[]");
1405  /* keep track of whether pts in each ellipse are doubled or not */
1406  pts_dbl = (int *)bu_malloc(nell * sizeof(int), "dbl ints");
1407 
1408  /* make ellipses at different levels in the +Z direction */
1409  i = 0;
1410  nseg = 0;
1411  theta_prev = M_2PI;
1412  pos_a = pts_a->next; /* skip over lower cusp of heart ( at pts_a ) */
1413  pos_b = pts_b->next;
1414  while (pos_a) {
1415  VSCALE(A, Au, pos_a->p[Y]); /* semimajor axis */
1416  VSCALE(B, Bu, pos_b->p[Y] * 0.80); /* semiminor axis */
1417  VJOIN1(V, hip->v, pos_a->p[Z], Cu);
1418 
1419  VSET(p1, 0.00, pos_b->p[Y], 0.00);
1420  theta_new = ell_angle(p1, pos_a->p[Y], pos_b->p[Y], dtol, ntol);
1421  if (nseg == 0) {
1422  nseg = (int)(M_2PI / theta_new) + 1;
1423  pts_dbl[i] = 0;
1424  } else if (theta_new < theta_prev) {
1425  nseg *= 2;
1426  pts_dbl[i] = 1;
1427  } else
1428  pts_dbl[i] = 0;
1429  theta_prev = theta_new;
1430 
1431  ellipses[i] = (fastf_t *)bu_malloc(3*(nseg+1)*sizeof(fastf_t),"pts ell");
1432  rt_ell(ellipses[i], V, A, B, nseg);
1433 
1434  i++;
1435  pos_a = pos_a->next;
1436  pos_b = pos_b->next;
1437  }
1438 
1439  /* Draw the top (largest) ellipse in the XY plane */
1440  RT_ADD_VLIST(vhead,
1441  &ellipses[nell-1][(nseg-1)*ELEMENTS_PER_VECT],
1443  for (i = 0; i < nseg; i++) {
1444  RT_ADD_VLIST(vhead,
1445  &ellipses[nell-1][i*ELEMENTS_PER_VECT],
1447  }
1448 
1449  /* Connect ellipses skipping the top (largest) ellipse which is at index nell - 1 */
1450  for (i = nell-2; i >= 0; i--) {
1451  ellipse_above = i + 1;
1452  ellipse_below = i;
1453  if (pts_dbl[ellipse_above])
1454  nseg /= 2; /* Number of line segments in 'ellipse_below' ellipse is halved if points double */
1455 
1456  /* Draw the current ellipse */
1457  RT_ADD_VLIST(vhead,
1458  &ellipses[ellipse_below][(nseg-1)*ELEMENTS_PER_VECT],
1460  for (j = 0; j < nseg; j++) {
1461  RT_ADD_VLIST(vhead,
1462  &ellipses[ellipse_below][j*ELEMENTS_PER_VECT],
1464  }
1465 
1466  /*
1467  * Make connections between ellipses
1468  * This for repetition statement draws connectors between
1469  * the nell-1 ellipses drawn above
1470  */
1471 
1472  for (j = 1; j < nseg; j++) {
1473  if (pts_dbl[ellipse_above])
1474  jj = j + j;
1475  else
1476  jj = j;
1477  RT_ADD_VLIST(vhead,
1478  &ellipses[ellipse_below][j*ELEMENTS_PER_VECT],
1480  RT_ADD_VLIST(vhead,
1481  &ellipses[ellipse_above][jj*ELEMENTS_PER_VECT],
1483  }
1484  }
1485 
1486  VSET(lower_cusp, 0, 0, hip->zdir[Z] * -1.00);
1487  VADD2(Work, hip->v, lower_cusp);
1488  for (i = 0; i < nseg; i++) {
1489  /* Draw connector */
1490  RT_ADD_VLIST(vhead, Work, BN_VLIST_LINE_MOVE);
1491  RT_ADD_VLIST(vhead,
1492  &ellipses[0][i*ELEMENTS_PER_VECT],
1494  }
1495 
1496  /* free memory */
1497  for (i = 0; i < nell; i++) {
1498  bu_free((char *)ellipses[i], "pts ell");
1499  }
1500  bu_free((char *)ellipses, "fastf_t ell[]");
1501  bu_free((char *)pts_dbl, "dbl ints");
1502 
1503  return 0;
1504 }
1505 
1506 
1507 int
1509 {
1510  bu_log("rt_hrt_tess: Not implemented yet!\n");
1511  return 0;
1512 }
1513 
1514 
1515 /**
1516  * The external form is:
1517  * V point
1518  * Xdir vector
1519  * Ydir vector
1520  * Zdir vector
1521  */
1522 int
1523 rt_hrt_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1524 {
1525  struct rt_hrt_internal *hip;
1526 
1527  /* must be double for import and export */
1528  double hec[ELEMENTS_PER_VECT*4 + 1];
1529 
1530  if (dbip) RT_CK_DBI(dbip);
1531 
1532  RT_CK_DB_INTERNAL(ip);
1533  if (ip->idb_type != ID_HRT) return -1;
1534  hip = (struct rt_hrt_internal *)ip->idb_ptr;
1535  RT_HRT_CK_MAGIC(hip);
1536 
1537  BU_CK_EXTERNAL(ep);
1538  ep->ext_nbytes = SIZEOF_NETWORK_DOUBLE * (ELEMENTS_PER_VECT*4 + 1);
1539  ep->ext_buf = (uint8_t *)bu_malloc(ep->ext_nbytes, "hrt external");
1540 
1541  /* scale values to local buffer */
1542  VSCALE(&hec[0*ELEMENTS_PER_VECT], hip->v, local2mm);
1543  VSCALE(&hec[1*ELEMENTS_PER_VECT], hip->xdir, local2mm);
1544  VSCALE(&hec[2*ELEMENTS_PER_VECT], hip->ydir, local2mm);
1545  VSCALE(&hec[3*ELEMENTS_PER_VECT], hip->zdir, local2mm);
1546  hec[4*ELEMENTS_PER_VECT] = hip->d;
1547 
1548  /* Convert from internal (host) to database (network) format */
1549  bu_cv_htond(ep->ext_buf, (unsigned char *)hec, ELEMENTS_PER_VECT*4 + 1);
1550 
1551  return 0;
1552 }
1553 
1554 
1555 /**
1556  * Import a heart from the database format to the internal format.
1557  *
1558  */
1559 int
1560 rt_hrt_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
1561 {
1562  struct rt_hrt_internal *hip;
1563 
1564  /* must be double for import and export */
1565  double hec[ELEMENTS_PER_VECT*4 + 1];
1566 
1567  if (dbip) RT_CK_DBI(dbip);
1568 
1569  RT_CK_DB_INTERNAL(ip);
1570  BU_CK_EXTERNAL(ep);
1571 
1572  BU_ASSERT_LONG(ep->ext_nbytes, ==, SIZEOF_NETWORK_DOUBLE * (ELEMENTS_PER_VECT*4 + 1));
1573 
1574  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1575  ip->idb_type = ID_HRT;
1576  ip->idb_meth = &OBJ[ID_HRT];
1577  BU_ALLOC(ip->idb_ptr, struct rt_hrt_internal);
1578 
1579  hip = (struct rt_hrt_internal *)ip->idb_ptr;
1580  hip->hrt_magic = RT_HRT_INTERNAL_MAGIC;
1581 
1582  /* Convert from database(network) to internal (host) format */
1583  bu_cv_ntohd((unsigned char *)hec, ep->ext_buf, ELEMENTS_PER_VECT*4 + 1);
1584 
1585  /* Apply modelling transformations */
1586  if (mat == NULL) mat = bn_mat_identity;
1587  MAT4X3PNT(hip->v, mat, &hec[0*ELEMENTS_PER_VECT]);
1588  MAT4X3PNT(hip->xdir, mat, &hec[1*ELEMENTS_PER_VECT]);
1589  MAT4X3PNT(hip->ydir, mat, &hec[2*ELEMENTS_PER_VECT]);
1590  MAT4X3PNT(hip->zdir, mat, &hec[3*ELEMENTS_PER_VECT]);
1591  hip->d = hec[4*ELEMENTS_PER_VECT];
1592 
1593  return 0; /* OK */
1594 }
1595 
1596 
1597 /**
1598  * Make human-readable formatted presentation of this solid. First
1599  * line describes type of solid. Additional lines are indented one
1600  * tab, and give parameter values.
1601  */
1602 int
1603 rt_hrt_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
1604 {
1605  struct rt_hrt_internal *hip =
1606  (struct rt_hrt_internal *)ip->idb_ptr;
1607  fastf_t mag_x, mag_y, mag_z;
1608  char buf[256];
1609  double angles[5];
1610  vect_t unit_vector;
1611 
1612  RT_HRT_CK_MAGIC(hip);
1613  bu_vls_strcat(str, "Heart (HRT)\n");
1614 
1615  sprintf(buf, "\tV (%g, %g, %g)\n",
1616  INTCLAMP(hip->v[X] * mm2local),
1617  INTCLAMP(hip->v[Y] * mm2local),
1618  INTCLAMP(hip->v[Z] * mm2local));
1619  bu_vls_strcat(str, buf);
1620 
1621  mag_x = MAGNITUDE(hip->xdir);
1622  mag_y = MAGNITUDE(hip->ydir);
1623  mag_z = MAGNITUDE(hip->zdir);
1624 
1625  sprintf(buf, "\tXdir (%g, %g, %g) mag=%g\n",
1626  INTCLAMP(hip->xdir[X] * mm2local),
1627  INTCLAMP(hip->xdir[Y] * mm2local),
1628  INTCLAMP(hip->xdir[Z] * mm2local),
1629  INTCLAMP(mag_x * mm2local));
1630  bu_vls_strcat(str, buf);
1631 
1632  sprintf(buf, "\tYdir (%g, %g, %g) mag=%g\n",
1633  INTCLAMP(hip->ydir[X] * mm2local),
1634  INTCLAMP(hip->ydir[Y] * mm2local),
1635  INTCLAMP(hip->ydir[Z] * mm2local),
1636  INTCLAMP(mag_y * mm2local));
1637  bu_vls_strcat(str, buf);
1638 
1639  sprintf(buf, "\tZdir (%g, %g, %g) mag=%g\n",
1640  INTCLAMP(hip->zdir[X] * mm2local),
1641  INTCLAMP(hip->zdir[Y] * mm2local),
1642  INTCLAMP(hip->zdir[Z] * mm2local),
1643  INTCLAMP(mag_z * mm2local));
1644  bu_vls_strcat(str, buf);
1645 
1646  sprintf(buf, "\td=%g\n", INTCLAMP(hip->d));
1647  bu_vls_strcat(str, buf);
1648 
1649  if (!verbose) return 0;
1650 
1651  VSCALE(unit_vector, hip->xdir, 1/mag_x);
1652  rt_find_fallback_angle(angles, unit_vector);
1653  rt_pr_fallback_angle(str, "\tXdir", angles);
1654 
1655  VSCALE(unit_vector, hip->ydir, 1/mag_y);
1656  rt_find_fallback_angle(angles, unit_vector);
1657  rt_pr_fallback_angle(str, "\tYdir", angles);
1658 
1659  VSCALE(unit_vector, hip->zdir, 1/mag_z);
1660  rt_find_fallback_angle(angles, unit_vector);
1661  rt_pr_fallback_angle(str, "\tZdir", angles);
1662 
1663  return 0;
1664 }
1665 
1666 
1667 /**
1668  * Free the storage associated with the rt_db_internal version of this
1669  * solid.
1670  */
1671 void
1673 {
1674  register struct rt_hrt_internal *hip;
1675 
1676  RT_CK_DB_INTERNAL(ip);
1677 
1678  hip = (struct rt_hrt_internal *)ip->idb_ptr;
1679  RT_HRT_CK_MAGIC(hip);
1680 
1681  bu_free((char *)hip, "rt_hrt_");
1682  ip->idb_ptr = ((void *)0);
1683 }
1684 
1685 
1686 int
1687 rt_hrt_params(struct pc_pc_set *UNUSED(ps), const struct rt_db_internal *ip)
1688 {
1689  if (ip) RT_CK_DB_INTERNAL(ip);
1690 
1691  return 0; /* OK */
1692 }
1693 
1694 
1695 void
1696 rt_hrt_surf_area(fastf_t *area, const struct rt_db_internal *ip)
1697 {
1698  fastf_t area_hrt_YZ_plane;
1699  struct rt_hrt_internal *hip = (struct rt_hrt_internal *)ip->idb_ptr;
1700  RT_HRT_CK_MAGIC(hip);
1701 
1702  /* Area of ellipse in YZ-plane is PI * ydir[Y] * ( zdir[Z] * 1.125 ) */
1703  area_hrt_YZ_plane = M_PI * hip->ydir[Y] * hip->zdir[Z] * 1.125;
1704 
1705  /* Area of heart = 180 * M_PI * Area of ellipse in YZ-plane
1706  * The 180 * M_PI scalar comes from http://mathworld.wolfram.com/HeartCurve.html
1707  */
1708  *area = 180 * M_PI * area_hrt_YZ_plane;
1709 }
1710 
1711 
1712 void
1714 {
1715  bu_log("rt_hrt_volume: Not implemented yet!\n");
1716 }
1717 
1718 /**
1719  * Computes centroid of a heart
1720  */
1721 void
1722 rt_hrt_centroid(point_t *cent, const struct rt_db_internal *ip)
1723 {
1724  struct rt_hrt_internal *hip = (struct rt_hrt_internal *)ip->idb_ptr;
1725  RT_HRT_CK_MAGIC(hip);
1726  VSET(*cent, hip->xdir[X], hip->ydir[Y], hip->zdir[Z] * 0.125);
1727 }
1728 
1729 
1730 /*
1731  * Local Variables:
1732  * mode: C
1733  * tab-width: 8
1734  * indent-tabs-mode: t
1735  * c-file-style: "stroustrup"
1736  * End:
1737  * ex: shiftwidth=4 tabstop=8
1738  */
vect_t hrt_V
Definition: hrt.c:139
void rt_hrt_ifree(struct rt_db_internal *ip)
Definition: hrt.c:1672
char * d_namep
pointer to name string
Definition: raytrace.h:859
vect_t hrt_invsq
Definition: hrt.c:144
Definition: raytrace.h:800
fastf_t cf[BN_MAX_POLY_DEGREE+1]
Definition: poly.h:50
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define BU_LIST_INSERT(old, new)
Definition: list.h:183
#define SIZEOF_NETWORK_DOUBLE
Definition: cv.h:48
struct hit seg_in
IN information.
Definition: raytrace.h:370
Definition: list.h:118
const struct directory * st_dp
Directory entry of solid.
Definition: raytrace.h:436
double dist
>= 0
Definition: tol.h:73
void rt_hrt_volume(void)
Definition: hrt.c:1713
const mat_t bn_mat_identity
Matrix and vector functionality.
Definition: mat.c:46
void rt_hrt_print(register const struct soltab *stp)
Definition: hrt.c:344
int rt_hrt_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: hrt.c:1523
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)
#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
void rt_hrt_curve(void)
Definition: hrt.c:966
struct bn_complex bn_complex_t
Complex numbers.
fastf_t st_aradius
Radius of APPROXIMATING sphere.
Definition: raytrace.h:433
mat_t hrt_invRSSR
Definition: hrt.c:146
fastf_t ell_angle(fastf_t *p1, fastf_t a, fastf_t b, fastf_t dtol, fastf_t ntol)
Definition: ell.c:1788
Header file for the BRL-CAD common definitions.
Definition: poly.h:47
struct resource * a_resource
dynamic memory resources
Definition: raytrace.h:1591
void bu_cv_htond(unsigned char *out, const unsigned char *in, size_t count)
void rt_pt_sort(fastf_t t[], int npts)
Definition: tgc.c:1394
struct bu_list l
Definition: raytrace.h:369
int rt_mk_hyperbola(struct rt_pt_node *pts, fastf_t r, fastf_t b, fastf_t c, fastf_t dtol, fastf_t ntol)
Definition: rhc.c:1101
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
void rt_hrt_centroid(point_t *cent, const struct rt_db_internal *ip)
Definition: hrt.c:1722
vect_t hrt_Z
Definition: hrt.c:142
int idb_major_type
Definition: raytrace.h:192
size_t dgr
Definition: poly.h:49
vect_t hrt_Y
Definition: hrt.c:141
Definition: color.c:49
void rt_hrt_free(struct soltab *stp)
Definition: hrt.c:980
struct rt_i * a_rt_i
this librt instance
Definition: raytrace.h:1588
int rt_hrt_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *tol)
Definition: hrt.c:155
vect_t hit_vpriv
PRIVATE vector for xxx_*()
Definition: raytrace.h:253
#define RT_ADD_VLIST(hd, pnt, draw)
Definition: raytrace.h:1865
void bn_mat_print(const char *title, const mat_t m)
Definition: mat.c:81
#define RT_CK_DB_INTERNAL(_p)
Definition: raytrace.h:207
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
fastf_t st_bradius
Radius of BOUNDING sphere.
Definition: raytrace.h:434
#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 primitive_get_absolute_tolerance(const struct rt_tess_tol *ttol, fastf_t rel_to_abs)
int rt_hrt_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: hrt.c:227
#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_hrt_norm(register struct hit *hitp, register struct xray *rp)
Definition: hrt.c:944
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
int rt_hrt_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
Definition: hrt.c:394
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)
Support for uniform tolerances.
Definition: tol.h:71
#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
int rt_hrt_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
Definition: hrt.c:1560
#define RT_HRT_INTERNAL_MAGIC
Definition: magic.h:116
struct bn_tol rti_tol
Math tolerances for this model.
Definition: raytrace.h:1765
struct rt_pt_node * next
ptr to next pt
Definition: raytrace.h:1912
#define bu_offsetof(_t, _m)
Definition: parse.h:64
#define RT_CK_DBI(_p)
Definition: raytrace.h:829
const struct bu_structparse rt_hrt_parse[]
Definition: hrt.c:128
#define RT_GET_SEG(p, res)
Definition: raytrace.h:379
void * idb_ptr
Definition: raytrace.h:195
point_t r_pt
Point at which ray starts.
Definition: raytrace.h:218
point_t st_min
min X, Y, Z of bounding RPP
Definition: raytrace.h:437
void bu_cv_ntohd(unsigned char *out, const unsigned char *in, size_t count)
const struct rt_functab OBJ[]
Definition: table.c:159
int rt_poly_roots(bn_poly_t *eqn, bn_complex_t roots[], const char *name)
mat_t hrt_SoR
Definition: hrt.c:145
void rt_hrt_vshot(struct soltab **stp, struct xray **rp, struct seg *segp, int n, struct application *ap)
Definition: hrt.c:640
int rt_hrt_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
Definition: hrt.c:1603
void rt_hrt_surf_area(fastf_t *area, const struct rt_db_internal *ip)
Definition: hrt.c:1696
mat_t hrt_invR
Definition: hrt.c:147
#define R
Definition: msr.c:55
void rt_ell(fastf_t *ov, const fastf_t *V, const fastf_t *A, const fastf_t *B, int sides)
Definition: epa.c:1141
void * st_specific
-> ID-specific (private) struct
Definition: raytrace.h:435
point_t p
a point
Definition: raytrace.h:1911
#define A
Definition: msr.c:51
Definition: color.c:51
void rt_hrt_24pts(fastf_t *ov, fastf_t *V, fastf_t *A, fastf_t *B)
Definition: hrt.c:995
#define RT_HRT_SEG_MISS(SEG)
Definition: hrt.c:635
int rt_hrt_adaptive_plot(void)
Definition: hrt.c:1042
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
vect_t hit_normal
DEPRECATED: Surface Normal at hit_point, use RT_HIT_NORMAL.
Definition: raytrace.h:252
#define BU_CK_LIST_HEAD(_p)
Definition: list.h:142
#define BU_CK_EXTERNAL(_p)
Definition: parse.h:224
fastf_t hrt_d
Definition: hrt.c:143
Complex numbers.
Definition: complex.h:39
size_t ext_nbytes
Definition: parse.h:210
#define ID_HRT
Heart.
Definition: raytrace.h:514
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
int rt_hrt_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: hrt.c:1050
double fastf_t
Definition: defines.h:300
#define VPRINT(a, b)
Definition: raytrace.h:1881
int rt_hrt_params(struct pc_pc_set *ps, const struct rt_db_internal *ip)
Definition: hrt.c:1687
void rt_hrt_uv(void)
Definition: hrt.c:973
vect_t hrt_X
Definition: hrt.c:140
Definition: color.c:50
#define HRTOUT(n)
Definition: hrt.c:993
point_t st_center
Centroid of solid.
Definition: raytrace.h:432
int rt_hrt_tess(void)
Definition: hrt.c:1508