BRL-CAD
rec.c
Go to the documentation of this file.
1 /* R E C . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1985-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 /** @addtogroup primitives */
21 /** @{ */
22 /** @file primitives/rec/rec.c
23  *
24  * Intersect a ray with a Right Elliptical Cylinder. This is a special
25  * (but common) case of the TGC, which is handled separately.
26  *
27  * Algorithm -
28  *
29  * Given V, H, A, and B, there is a set of points on this cylinder
30  *
31  * { (x, y, z) | (x, y, z) is on cylinder }
32  *
33  * Through a series of Affine Transformations, this set of points will
34  * be transformed into a set of points on a unit cylinder located at
35  * the origin with a radius of 1, and a height of +1 along the +Z
36  * axis.
37  *
38  * { (x', y', z') | (x', y', z') is on cylinder at origin }
39  *
40  * The transformation from X to X' is accomplished by:
41  *
42  * X' = S(R(X - V))
43  *
44  * where R(X) = (A/(|A|))
45  * (B/(|B|)) . X
46  * (H/(|H|))
47  *
48  * and S(X) = (1/|A| 0 0)
49  * (0 1/|B| 0) . X
50  * (0 0 1/|H|)
51  *
52  * To find the intersection of a line with the surface of the
53  * cylinder, consider the parametric line L:
54  *
55  * L : { P(n) | P + t(n) . D }
56  *
57  * Call W the actual point of intersection between L and the cylinder.
58  * Let W' be the point of intersection between L' and the unit
59  * cylinder.
60  *
61  * L' : { P'(n) | P' + t(n) . D' }
62  *
63  * W = invR(invS(W')) + V
64  *
65  * Where W' = k D' + P'.
66  *
67  * If Dx' and Dy' are both 0, then there is no hit on the cylinder;
68  * but the end plates need checking.
69  *
70  * Line L' hits the infinitely tall unit cylinder at W' when
71  *
72  * k**2 + b * k + c = 0
73  *
74  * where
75  *
76  * b = 2 * (Dx' * Px' + Dy' * Py') / (Dx'**2 + Dy'**2)
77  * c = ((Px'**2 + Py'**2) - r**2) / (Dx'**2 + Dy'**2)
78  * r = 1.0
79  *
80  * The quadratic formula yields k (which is constant):
81  *
82  * k = [ -b +/- sqrt(b**2 - 4 * c ] / 2.0
83  *
84  * Now, D' = S(R(D))
85  * and P' = S(R(P - V))
86  *
87  * Substituting,
88  *
89  * W = V + invR(invS[ k *(S(R(D))) + S(R(P - V)) ])
90  * = V + invR((k * R(D)) + R(P - V))
91  * = V + k * D + P - V
92  * = k * D + P
93  *
94  * Note that ``k'' is constant, and is the same in the formulations
95  * for both W and W'.
96  *
97  * The hit at ``k'' is a hit on the height=1 unit cylinder IFF
98  * 0 <= Wz' <= 1.
99  *
100  * NORMALS. Given the point W on the surface of the cylinder, what is
101  * the vector normal to the tangent plane at that point?
102  *
103  * Map W onto the unit cylinder, i.e.: W' = S(R(W - V)).
104  *
105  * Plane on unit cylinder at W' has a normal vector N' of the same
106  * value as W' in x and y, with z set to zero, i.e., (Wx', Wy', 0)
107  *
108  * The plane transforms back to the tangent plane at W, and this new
109  * plane (on the original cylinder) has a normal vector of N, viz:
110  *
111  * N = inverse[ transpose(invR o invS) ] (N')
112  * = inverse[ transpose(invS) o transpose(invR) ] (N')
113  * = inverse[ inverse(S) o R ] (N')
114  * = invR o S (N')
115  *
116  * Note that the normal vector produced above will not have unit
117  * length.
118  *
119  * THE END PLATES.
120  *
121  * If Dz' == 0, line L' is parallel to the end plates, so there is no
122  * hit.
123  *
124  * Otherwise, the line L' hits the bottom plate with k = (0 - Pz') / Dz',
125  * and hits the top plate with k = (1 - Pz') / Dz'.
126  *
127  * The solution W' is within the end plate IFF
128  *
129  * Wx'**2 + Wy'**2 <= 1.0
130  *
131  * The normal for a hit on the bottom plate is -Hunit, and the normal
132  * for a hit on the top plate is +Hunit.
133  *
134  */
135 
136 #include "common.h"
137 
138 #include <string.h>
139 #include <math.h>
140 #include "bio.h"
141 
142 #include "vmath.h"
143 #include "rtgeom.h"
144 #include "raytrace.h"
145 
146 
147 struct rec_specific {
148  vect_t rec_V; /* Vector to center of base of cylinder */
149  vect_t rec_Hunit; /* Unit H vector */
150  mat_t rec_SoR; /* Scale(Rot(vect)) */
151  mat_t rec_invRoS; /* invRot(Scale(vect)) */
152  vect_t rec_A; /* One axis of ellipse */
153  vect_t rec_B; /* Other axis of ellipse */
154  fastf_t rec_iAsq; /* 1/MAGSQ(A) */
155  fastf_t rec_iBsq; /* 1/MAGSQ(B) */
156 };
157 
158 /**
159  * Calculate the RPP for an REC
160  */
161 int
162 rt_rec_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol)) {
163  mat_t R;
164  vect_t P, w1;
165  fastf_t f, tmp, z;
166  double magsq_h, magsq_a, magsq_b, magsq_c, magsq_d;
167  double mag_h, mag_a, mag_b;
168  struct rt_tgc_internal *tip;
169 
170  RT_CK_DB_INTERNAL(ip);
171 
172  tip = (struct rt_tgc_internal *)ip->idb_ptr;
173  RT_TGC_CK_MAGIC(tip);
174 
175  mag_h = sqrt(magsq_h = MAGSQ(tip->h));
176  mag_a = sqrt(magsq_a = MAGSQ(tip->a));
177  mag_b = sqrt(magsq_b = MAGSQ(tip->b));
178  magsq_c = MAGSQ(tip->c);
179  magsq_d = MAGSQ(tip->d);
180 
181  MAT_IDN(R);
182  f = 1.0/mag_a;
183  VSCALE(&R[0], tip->a, f);
184  f = 1.0/mag_b;
185  VSCALE(&R[4], tip->b, f);
186  f = 1.0/mag_h;
187  VSCALE(&R[8], tip->h, f);
188 
189  /* X */
190  VSET(P, 1.0, 0, 0); /* bounding plane normal */
191  MAT3X3VEC(w1, R, P); /* map plane into local coord syst */
192  /* 1st end ellipse (no Z part) */
193  tmp = magsq_a * w1[X] * w1[X] + magsq_b * w1[Y] * w1[Y];
194  if (tmp > SMALL)
195  f = sqrt(tmp); /* XY part */
196  else
197  f = 0;
198  (*min)[X] = tip->v[X] - f; /* V.P +/- f */
199  (*max)[X] = tip->v[X] + f;
200  /* 2nd end ellipse */
201  z = w1[Z] * mag_h; /* Z part */
202  tmp = magsq_c * w1[X] * w1[X] + magsq_d * w1[Y] * w1[Y];
203  if (tmp > SMALL)
204  f = sqrt(tmp); /* XY part */
205  else
206  f = 0;
207  if (tip->v[X] - f + z < (*min)[X])
208  (*min)[X] = tip->v[X] - f + z;
209  if (tip->v[X] + f + z > (*max)[X])
210  (*max)[X] = tip->v[X] + f + z;
211 
212  /* Y */
213  VSET(P, 0, 1.0, 0); /* bounding plane normal */
214  MAT3X3VEC(w1, R, P); /* map plane into local coord syst */
215  /* 1st end ellipse (no Z part) */
216  tmp = magsq_a * w1[X] * w1[X] + magsq_b * w1[Y] * w1[Y];
217  if (tmp > SMALL)
218  f = sqrt(tmp); /* XY part */
219  else
220  f = 0;
221  (*min)[Y] = tip->v[Y] - f; /* V.P +/- f */
222  (*max)[Y] = tip->v[Y] + f;
223  /* 2nd end ellipse */
224  z = w1[Z] * mag_h; /* Z part */
225  tmp = magsq_c * w1[X] * w1[X] + magsq_d * w1[Y] * w1[Y];
226  if (tmp > SMALL)
227  f = sqrt(tmp); /* XY part */
228  else
229  f = 0;
230  if (tip->v[Y] - f + z < (*min)[Y])
231  (*min)[Y] = tip->v[Y] - f + z;
232  if (tip->v[Y] + f + z > (*max)[Y])
233  (*max)[Y] = tip->v[Y] + f + z;
234 
235  /* Z */
236  VSET(P, 0, 0, 1.0); /* bounding plane normal */
237  MAT3X3VEC(w1, R, P); /* map plane into local coord syst */
238  /* 1st end ellipse (no Z part) */
239  tmp = magsq_a * w1[X] * w1[X] + magsq_b * w1[Y] * w1[Y];
240  if (tmp > SMALL)
241  f = sqrt(tmp); /* XY part */
242  else
243  f = 0;
244  (*min)[Z] = tip->v[Z] - f; /* V.P +/- f */
245  (*max)[Z] = tip->v[Z] + f;
246  /* 2nd end ellipse */
247  z = w1[Z] * mag_h; /* Z part */
248  tmp = magsq_c * w1[X] * w1[X] + magsq_d * w1[Y] * w1[Y];
249  if (tmp > SMALL)
250  f = sqrt(tmp); /* XY part */
251  else
252  f = 0;
253  if (tip->v[Z] - f + z < (*min)[Z])
254  (*min)[Z] = tip->v[Z] - f + z;
255  if (tip->v[Z] + f + z > (*max)[Z])
256  (*max)[Z] = tip->v[Z] + f + z;
257 
258  return 0;
259 }
260 
261 /**
262  * Given a pointer to a GED database record, and a transformation matrix,
263  * determine if this is a valid REC,
264  * and if so, precompute various terms of the formulas.
265  *
266  * Returns -
267  * 0 REC is OK
268  * !0 Error in description
269  *
270  * Implicit return - A struct rec_specific is created, and its
271  * address is stored in stp->st_specific for use by rt_rec_shot(). If
272  * the TGC is really an REC, stp->st_id is modified to ID_REC.
273  */
274 int
275 rt_rec_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
276 {
277  struct rt_tgc_internal *tip;
278  struct rec_specific *rec;
279  double magsq_h, magsq_a, magsq_b;
280  double mag_h, mag_a, mag_b;
281  mat_t R;
282  mat_t Rinv;
283  mat_t S;
284  vect_t invsq; /* [ 1/(|A|**2), 1/(|B|**2), 1/(|Hv|**2) ] */
285  vect_t work;
286  fastf_t f;
287 
288  if (!stp || !ip)
289  return -1;
290  RT_CK_SOLTAB(stp);
291  RT_CK_DB_INTERNAL(ip);
292  if (rtip) RT_CK_RTI(rtip);
293 
294  tip = (struct rt_tgc_internal *)ip->idb_ptr;
295  RT_TGC_CK_MAGIC(tip);
296 
297  /* Validate that |H| > 0, compute |A| |B| |C| |D| */
298  mag_h = sqrt(magsq_h = MAGSQ(tip->h));
299  mag_a = sqrt(magsq_a = MAGSQ(tip->a));
300  mag_b = sqrt(magsq_b = MAGSQ(tip->b));
301 
302  /* Check for |H| > 0, |A| > 0, |B| > 0 */
303  if (NEAR_ZERO(mag_h, RT_LEN_TOL) || NEAR_ZERO(mag_a, RT_LEN_TOL)
304  || NEAR_ZERO(mag_b, RT_LEN_TOL)) {
305  return 1; /* BAD, too small */
306  }
307 
308  /* Make sure that A == C, B == D */
309  VSUB2(work, tip->a, tip->c);
310  f = MAGNITUDE(work);
311  if (! NEAR_ZERO(f, RT_LEN_TOL)) {
312  return 1; /* BAD, !cylinder */
313  }
314  VSUB2(work, tip->b, tip->d);
315  f = MAGNITUDE(work);
316  if (! NEAR_ZERO(f, RT_LEN_TOL)) {
317  return 1; /* BAD, !cylinder */
318  }
319 
320  /* Check for A.B == 0, H.A == 0 and H.B == 0 */
321  f = VDOT(tip->a, tip->b) / (mag_a * mag_b);
322  if (! NEAR_ZERO(f, RT_DOT_TOL)) {
323  return 1; /* BAD */
324  }
325  f = VDOT(tip->h, tip->a) / (mag_h * mag_a);
326  if (! NEAR_ZERO(f, RT_DOT_TOL)) {
327  return 1; /* BAD */
328  }
329  f = VDOT(tip->h, tip->b) / (mag_h * mag_b);
330  if (! NEAR_ZERO(f, RT_DOT_TOL)) {
331  return 1; /* BAD */
332  }
333 
334  /*
335  * This TGC is really an REC
336  */
337  stp->st_id = ID_REC; /* "fix" soltab ID */
338  stp->st_meth = &OBJ[ID_REC];
339 
340  BU_GET(rec, struct rec_specific);
341  stp->st_specific = (void *)rec;
342 
343  VMOVE(rec->rec_Hunit, tip->h);
344  VUNITIZE(rec->rec_Hunit);
345 
346  VMOVE(rec->rec_V, tip->v);
347  VMOVE(rec->rec_A, tip->a);
348  VMOVE(rec->rec_B, tip->b);
349  rec->rec_iAsq = 1.0/magsq_a;
350  rec->rec_iBsq = 1.0/magsq_b;
351 
352  VSET(invsq, 1.0/magsq_a, 1.0/magsq_b, 1.0/magsq_h);
353 
354  /* Compute R and Rinv matrices */
355  MAT_IDN(R);
356  f = 1.0/mag_a;
357  VSCALE(&R[0], tip->a, f);
358  f = 1.0/mag_b;
359  VSCALE(&R[4], tip->b, f);
360  f = 1.0/mag_h;
361  VSCALE(&R[8], tip->h, f);
362  bn_mat_trn(Rinv, R); /* inv of rot mat is trn */
363 
364  /* Compute S */
365  MAT_IDN(S);
366  S[ 0] = sqrt(invsq[0]);
367  S[ 5] = sqrt(invsq[1]);
368  S[10] = sqrt(invsq[2]);
369 
370  /* Compute SoR and invRoS */
371  bn_mat_mul(rec->rec_SoR, S, R);
372  bn_mat_mul(rec->rec_invRoS, Rinv, S);
373 
374  /* Compute bounding sphere and RPP */
375  {
376  fastf_t dx, dy, dz; /* For bounding sphere */
377 
378  if (stp->st_meth->ft_bbox(ip, &(stp->st_min), &(stp->st_max), &(rtip->rti_tol))) return 1;
379 
380  VSET(stp->st_center,
381  (stp->st_max[X] + stp->st_min[X])/2,
382  (stp->st_max[Y] + stp->st_min[Y])/2,
383  (stp->st_max[Z] + stp->st_min[Z])/2);
384 
385  dx = (stp->st_max[X] - stp->st_min[X])/2;
386  f = dx;
387  dy = (stp->st_max[Y] - stp->st_min[Y])/2;
388  if (dy > f) f = dy;
389  dz = (stp->st_max[Z] - stp->st_min[Z])/2;
390  if (dz > f) f = dz;
391  stp->st_aradius = f;
392  stp->st_bradius = sqrt(dx*dx + dy*dy + dz*dz);
393  }
394  return 0; /* OK */
395 }
396 
397 
398 void
399 rt_rec_print(const struct soltab *stp)
400 {
401  const struct rec_specific *rec =
402  (struct rec_specific *)stp->st_specific;
403 
404  VPRINT("V", rec->rec_V);
405  VPRINT("Hunit", rec->rec_Hunit);
406  bn_mat_print("S o R", rec->rec_SoR);
407  bn_mat_print("invR o S", rec->rec_invRoS);
408 }
409 
410 
411 /* hit_surfno is set to one of these */
412 #define REC_NORM_BODY (1) /* compute normal */
413 #define REC_NORM_TOP (2) /* copy tgc_N */
414 #define REC_NORM_BOT (3) /* copy reverse tgc_N */
415 
416 
417 /**
418  * Intersect a ray with a right elliptical cylinder,
419  * where all constant terms have
420  * been precomputed by rt_rec_prep(). If an intersection occurs,
421  * a struct seg will be acquired and filled in.
422  *
423  * Returns -
424  * 0 MISS
425  * >0 HIT
426  */
427 int
428 rt_rec_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
429 {
430  struct rec_specific *rec;
431  vect_t dprime; /* D' */
432  vect_t pprime; /* P' */
433  fastf_t k1, k2; /* distance constants of solution */
434  vect_t xlated; /* translated vector */
435  struct hit hits[4] = /* 4 potential hit points */
436  {RT_HIT_INIT_ZERO, RT_HIT_INIT_ZERO, RT_HIT_INIT_ZERO, RT_HIT_INIT_ZERO};
437  struct hit *hitp = NULL; /* pointer to hit point */
438  int nhits = 0; /* Number of hit points */
439  fastf_t tol_dist;
440 
441  if (UNLIKELY(!stp || !rp || !ap || !seghead))
442  return 0;
443 
444  rec = (struct rec_specific *)stp->st_specific;
445 
446  hitp = &hits[0];
447 
448  /* out, Mat, vect */
449  MAT4X3VEC(dprime, rec->rec_SoR, rp->r_dir);
450  VSUB2(xlated, rp->r_pt, rec->rec_V);
451  MAT4X3VEC(pprime, rec->rec_SoR, xlated);
452 
453  /*
454  * Check for hitting the end plates.
455  */
456  if (!ZERO(dprime[Z])) {
457  k1 = -pprime[Z] / dprime[Z]; /* bottom plate */
458  k2 = (1.0 - pprime[Z]) / dprime[Z]; /* top plate */
459 
460  VJOIN1(hitp->hit_vpriv, pprime, k1, dprime); /* hit' */
461  if (hitp->hit_vpriv[X] * hitp->hit_vpriv[X] +
462  hitp->hit_vpriv[Y] * hitp->hit_vpriv[Y] < (1.0 + SMALL_FASTF)) {
463  hitp->hit_magic = RT_HIT_MAGIC;
464  hitp->hit_dist = k1;
465  hitp->hit_surfno = REC_NORM_BOT; /* -H */
466  hitp++; nhits++;
467  }
468 
469  VJOIN1(hitp->hit_vpriv, pprime, k2, dprime); /* hit' */
470  if (hitp->hit_vpriv[X] * hitp->hit_vpriv[X] +
471  hitp->hit_vpriv[Y] * hitp->hit_vpriv[Y] < (1.0 + SMALL_FASTF)) {
472  hitp->hit_magic = RT_HIT_MAGIC;
473  hitp->hit_dist = k2;
474  hitp->hit_surfno = REC_NORM_TOP; /* +H */
475  hitp++; nhits++;
476  }
477  }
478 
479  /* Check for hitting the cylinder. Find roots of the equation,
480  * using formula for quadratic w/ a=1
481  */
482  if (nhits != 2) {
483  fastf_t b; /* coeff of polynomial */
484  fastf_t discriminant; /* root of radical, the discriminant */
485  fastf_t dx2dy2;
486 
487  dx2dy2 = 1 / (dprime[X]*dprime[X] + dprime[Y]*dprime[Y]);
488  b = 2 * (dprime[X]*pprime[X] + dprime[Y]*pprime[Y]) * dx2dy2;
489  discriminant = b*b - 4 * dx2dy2 * (pprime[X]*pprime[X] + pprime[Y]*pprime[Y] - 1);
490 
491  /* might want compare against tol_dist here? */
492 
493  if (NEAR_ZERO(discriminant, SMALL_FASTF)) {
494  /* if the discriminant is zero, it's a double-root grazer */
495  k1 = -b * 0.5;
496  VJOIN1(hitp->hit_vpriv, pprime, k1, dprime); /* hit' */
497  if (hitp->hit_vpriv[Z] > -SMALL_FASTF && hitp->hit_vpriv[Z] < (1.0 + SMALL_FASTF)) {
498  hitp->hit_magic = RT_HIT_MAGIC;
499  hitp->hit_dist = k1;
500  hitp->hit_surfno = REC_NORM_BODY; /* compute N */
501  hitp++; nhits++;
502  }
503 
504  } else if (discriminant > SMALL_FASTF) {
505  /* if the discriminant is positive, there are two roots */
506 
507  discriminant = sqrt(discriminant);
508  k1 = (-b+discriminant) * 0.5;
509  k2 = (-b-discriminant) * 0.5;
510 
511  /*
512  * k1 and k2 are potential solutions to intersection with side.
513  * See if they fall in range.
514  */
515  VJOIN1(hitp->hit_vpriv, pprime, k1, dprime); /* hit' */
516  if (hitp->hit_vpriv[Z] > -SMALL_FASTF && hitp->hit_vpriv[Z] < (1.0 + SMALL_FASTF)) {
517  hitp->hit_magic = RT_HIT_MAGIC;
518  hitp->hit_dist = k1;
519  hitp->hit_surfno = REC_NORM_BODY; /* compute N */
520  hitp++; nhits++;
521  }
522 
523  VJOIN1(hitp->hit_vpriv, pprime, k2, dprime); /* hit' */
524  if (hitp->hit_vpriv[Z] > -SMALL_FASTF && hitp->hit_vpriv[Z] < (1.0 + SMALL_FASTF)) {
525  hitp->hit_magic = RT_HIT_MAGIC;
526  hitp->hit_dist = k2;
527  hitp->hit_surfno = REC_NORM_BODY; /* compute N */
528  hitp++; nhits++;
529  }
530  }
531  }
532 
533  /* missed both ends and side? */
534  if (nhits == 0)
535  return 0;
536 
537  /* Prepare to collapse duplicate points. Check for case where two
538  * or more of the hits have the same distance, e.g. hitting at the
539  * rim or down an edge.
540  */
541 
542  tol_dist = ap->a_rt_i->rti_tol.dist;
543 
544  if (nhits > 3) {
545  /* collapse just one duplicate (4->3) */
546  if (NEAR_EQUAL(hits[0].hit_dist, hits[3].hit_dist, tol_dist)) {
548  bu_log("rt_rec_shot(%s): repeat hit, collapsing 0&3\n", stp->st_name);
549  nhits--; /* discard [3] */
550  } else if (NEAR_EQUAL(hits[1].hit_dist, hits[3].hit_dist, tol_dist)) {
552  bu_log("rt_rec_shot(%s): repeat hit, collapsing 1&3\n", stp->st_name);
553  nhits--; /* discard [3] */
554  } else if (NEAR_EQUAL(hits[2].hit_dist, hits[3].hit_dist, tol_dist)) {
556  bu_log("rt_rec_shot(%s): repeat hit, collapsing 2&3\n", stp->st_name);
557  nhits--; /* discard [3] */
558  }
559  }
560  if (nhits > 2) {
561  /* collapse any other duplicate (3->2) */
562  if (NEAR_EQUAL(hits[0].hit_dist, hits[2].hit_dist, tol_dist)) {
564  bu_log("rt_rec_shot(%s): repeat hit, collapsing 0&2\n", stp->st_name);
565  nhits--; /* discard [2] */
566  } else if (NEAR_EQUAL(hits[1].hit_dist, hits[2].hit_dist, tol_dist)) {
568  bu_log("rt_rec_shot(%s): repeat hit, collapsing 1&2\n", stp->st_name);
569  nhits--; /* discard [2] */
570  } else if (NEAR_EQUAL(hits[0].hit_dist, hits[1].hit_dist, tol_dist)) {
572  bu_log("rt_rec_shot(%s): repeat hit, collapsing 0&1\n", stp->st_name);
573  hits[1] = hits[2]; /* struct copy */
574  nhits--; /* moved [2] to [1], discarded [2] */
575  }
576  }
577 
578  /* sanity check that we don't end up with too many hits */
579  if (nhits > 3) {
580  bu_log("rt_rec_shot(%s): %d unique hits?!? %g, %g, %g, %g\n",
581  stp->st_name, nhits,
582  hits[0].hit_dist,
583  hits[1].hit_dist,
584  hits[2].hit_dist,
585  hits[3].hit_dist);
586  /* count just the first two hits, to have something */
587  nhits-=2;
588  } else if (nhits > 2) {
589  bu_log("rt_rec_shot(%s): %d unique hits?!? %g, %g, %g\n",
590  stp->st_name, nhits,
591  hits[0].hit_dist,
592  hits[1].hit_dist,
593  hits[2].hit_dist);
594  /* count just the first two hits, to have something */
595  nhits--;
596  } else if (nhits == 1) {
597  /* Ray is probably tangent to body of cylinder or a single hit
598  * on only an end plate. This could be considered a MISS, but
599  * to signal the condition, return 0-thickness hit.
600  */
601  hits[1] = hits[0]; /* struct copy */
602  nhits++; /* replicate [0] to [1] */
603  }
604 
605  if (hits[0].hit_dist < hits[1].hit_dist) {
606  /* entry is [0], exit is [1] */
607  struct seg *segp;
608 
609  RT_GET_SEG(segp, ap->a_resource);
610  segp->seg_stp = stp;
611  segp->seg_in = hits[0]; /* struct copy */
612  segp->seg_out = hits[1]; /* struct copy */
613  BU_LIST_INSERT(&(seghead->l), &(segp->l));
614  } else {
615  /* entry is [1], exit is [0] */
616  struct seg *segp;
617 
618  RT_GET_SEG(segp, ap->a_resource);
619  segp->seg_stp = stp;
620  segp->seg_in = hits[1]; /* struct copy */
621  segp->seg_out = hits[0]; /* struct copy */
622  BU_LIST_INSERT(&(seghead->l), &(segp->l));
623  }
624  return 2; /* HIT */
625 }
626 
627 
628 #define RT_REC_SEG_MISS(SEG) (SEG).seg_stp=(struct soltab *) 0;
629 /**
630  * This is the Becker vector version
631  */
632 void
633 rt_rec_vshot(struct soltab **stp, struct xray **rp, struct seg *segp, int n, struct application *ap)
634  /* An array of solid pointers */
635  /* An array of ray pointers */
636  /* array of segs (results returned) */
637  /* Number of ray/object pairs */
638 
639 {
640  int i;
641  struct rec_specific *rec;
642  vect_t dprime; /* D' */
643  vect_t pprime; /* P' */
644  fastf_t k1, k2; /* distance constants of solution */
645  vect_t xlated; /* translated vector */
646  struct hit hits[3]; /* 4 potential hit points */
647  struct hit *hitp; /* pointer to hit point */
648  fastf_t b; /* coeff of polynomial */
649  fastf_t root; /* root of radical */
650  fastf_t dx2dy2;
651 
652  if (ap) RT_CK_APPLICATION(ap);
653 
654  /* for each ray/right_elliptical_cylinder pair */
655  for (i = 0; i < n; i++) {
656  if (stp[i] == 0) continue; /* stp[i] == 0 signals skip ray */
657 
658  rec = (struct rec_specific *)stp[i]->st_specific;
659  hitp = &hits[0];
660 
661  /* out, Mat, vect */
662  MAT4X3VEC(dprime, rec->rec_SoR, rp[i]->r_dir);
663  VSUB2(xlated, rp[i]->r_pt, rec->rec_V);
664  MAT4X3VEC(pprime, rec->rec_SoR, xlated);
665 
666  if (ZERO(dprime[X]) && ZERO(dprime[Y]))
667  goto check_plates;
668 
669  /* Find roots of eqn, using formula for quadratic w/ a=1 */
670  b = 2 * (dprime[X]*pprime[X] + dprime[Y]*pprime[Y]) *
671  (dx2dy2 = 1 / (dprime[X]*dprime[X] + dprime[Y]*dprime[Y]));
672  if ((root = b*b - 4 * dx2dy2 *
673  (pprime[X]*pprime[X] + pprime[Y]*pprime[Y] - 1)) <= 0)
674  goto check_plates;
675 
676  root = sqrt(root);
677  k1 = (root-b) * 0.5;
678  k2 = (root+b) * (-0.5);
679 
680  /*
681  * k1 and k2 are potential solutions to intersection with side.
682  * See if they fall in range.
683  */
684  VJOIN1(hitp->hit_vpriv, pprime, k1, dprime); /* hit' */
685  if (hitp->hit_vpriv[Z] >= 0.0 && hitp->hit_vpriv[Z] <= 1.0) {
686  hitp->hit_dist = k1;
687  hitp->hit_surfno = REC_NORM_BODY; /* compute N */
688  hitp++;
689  }
690 
691  VJOIN1(hitp->hit_vpriv, pprime, k2, dprime); /* hit' */
692  if (hitp->hit_vpriv[Z] >= 0.0 && hitp->hit_vpriv[Z] <= 1.0) {
693  hitp->hit_dist = k2;
694  hitp->hit_surfno = REC_NORM_BODY; /* compute N */
695  hitp++;
696  }
697 
698  /*
699  * Check for hitting the end plates.
700  */
701  check_plates:
702  if (hitp < &hits[2] && !ZERO(dprime[Z])) {
703  /* 0 or 1 hits so far, this is worthwhile */
704  k1 = -pprime[Z] / dprime[Z]; /* bottom plate */
705  k2 = (1.0 - pprime[Z]) / dprime[Z]; /* top plate */
706 
707  VJOIN1(hitp->hit_vpriv, pprime, k1, dprime);/* hit' */
708  if (hitp->hit_vpriv[X] * hitp->hit_vpriv[X] +
709  hitp->hit_vpriv[Y] * hitp->hit_vpriv[Y] <= 1.0) {
710  hitp->hit_dist = k1;
711  hitp->hit_surfno = REC_NORM_BOT; /* -H */
712  hitp++;
713  }
714 
715  VJOIN1(hitp->hit_vpriv, pprime, k2, dprime);/* hit' */
716  if (hitp->hit_vpriv[X] * hitp->hit_vpriv[X] +
717  hitp->hit_vpriv[Y] * hitp->hit_vpriv[Y] <= 1.0) {
718  hitp->hit_dist = k2;
719  hitp->hit_surfno = REC_NORM_TOP; /* +H */
720  hitp++;
721  }
722  }
723 
724  if (hitp != &hits[2]) {
725  RT_REC_SEG_MISS(segp[i]); /* MISS */
726  } else {
727  segp[i].seg_stp = stp[i];
728 
729  if (hits[0].hit_dist < hits[1].hit_dist) {
730  /* entry is [0], exit is [1] */
731  VMOVE(segp[i].seg_in.hit_vpriv, hits[0].hit_vpriv);
732  segp[i].seg_in.hit_dist = hits[0].hit_dist;
733  segp[i].seg_in.hit_surfno = hits[0].hit_surfno;
734  VMOVE(segp[i].seg_out.hit_vpriv, hits[1].hit_vpriv);
735  segp[i].seg_out.hit_dist = hits[1].hit_dist;
736  segp[i].seg_out.hit_surfno = hits[1].hit_surfno;
737  } else {
738  /* entry is [1], exit is [0] */
739  VMOVE(segp[i].seg_in.hit_vpriv, hits[1].hit_vpriv);
740  segp[i].seg_in.hit_dist = hits[1].hit_dist;
741  segp[i].seg_in.hit_surfno = hits[1].hit_surfno;
742  VMOVE(segp[i].seg_out.hit_vpriv, hits[0].hit_vpriv);
743  segp[i].seg_out.hit_dist = hits[0].hit_dist;
744  segp[i].seg_out.hit_surfno = hits[0].hit_surfno;
745  }
746  }
747  }
748 }
749 
750 
751 /**
752  * Given ONE ray distance, return the normal and entry/exit point.
753  * hit_surfno is a flag indicating if normal needs to be computed or not.
754  */
755 void
756 rt_rec_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
757 {
758  struct rec_specific *rec =
759  (struct rec_specific *)stp->st_specific;
760 
761  VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
762  switch (hitp->hit_surfno) {
763  case REC_NORM_BODY:
764  /* compute it */
765  hitp->hit_vpriv[Z] = 0.0;
766  MAT4X3VEC(hitp->hit_normal, rec->rec_invRoS,
767  hitp->hit_vpriv);
768  VUNITIZE(hitp->hit_normal);
769  break;
770  case REC_NORM_TOP:
771  VMOVE(hitp->hit_normal, rec->rec_Hunit);
772  break;
773  case REC_NORM_BOT:
774  VREVERSE(hitp->hit_normal, rec->rec_Hunit);
775  break;
776  default:
777  bu_log("rt_rec_norm: surfno=%d bad\n", hitp->hit_surfno);
778  break;
779  }
780 }
781 
782 
783 /**
784  * Return the "curvature" of the cylinder. If an endplate,
785  * pick a principle direction orthogonal to the normal, and
786  * indicate no curvature. Otherwise, compute curvature.
787  * Normal must have been computed before calling this routine.
788  */
789 void
790 rt_rec_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
791 {
792  struct rec_specific *rec =
793  (struct rec_specific *)stp->st_specific;
794  vect_t uu;
795  fastf_t ax, bx, q;
796 
797  switch (hitp->hit_surfno) {
798  case REC_NORM_BODY:
799  /* This could almost certainly be simpler if we used
800  * inverse A rather than inverse A squared, right Ed?
801  */
802  VMOVE(cvp->crv_pdir, rec->rec_Hunit);
803  VSUB2(uu, hitp->hit_point, rec->rec_V);
804  cvp->crv_c1 = 0;
805  ax = VDOT(uu, rec->rec_A) * rec->rec_iAsq;
806  bx = VDOT(uu, rec->rec_B) * rec->rec_iBsq;
807  q = sqrt(ax * ax * rec->rec_iAsq + bx * bx * rec->rec_iBsq);
808  cvp->crv_c2 = - rec->rec_iAsq * rec->rec_iBsq / (q*q*q);
809  break;
810  case REC_NORM_TOP:
811  case REC_NORM_BOT:
812  bn_vec_ortho(cvp->crv_pdir, hitp->hit_normal);
813  cvp->crv_c1 = cvp->crv_c2 = 0;
814  break;
815  default:
816  bu_log("rt_rec_curve: bad surfno %d\n", hitp->hit_surfno);
817  break;
818  }
819 }
820 
821 
822 /**
823  * For a hit on the surface of an REC, return the (u, v) coordinates
824  * of the hit point, 0 <= u, v <= 1.
825  *
826  * u is the rotation around the cylinder, and
827  * v is the displacement along H.
828  */
829 void
830 rt_rec_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
831 {
832  struct rec_specific *rec =
833  (struct rec_specific *)stp->st_specific;
834 
835  vect_t work;
836  vect_t pprime;
837  fastf_t len;
838  fastf_t ratio;
839 
840  if (ap) RT_CK_APPLICATION(ap);
841 
842  /* hit_point is on surface; project back to unit cylinder,
843  * creating a vector from vertex to hit point.
844  */
845  VSUB2(work, hitp->hit_point, rec->rec_V);
846  MAT4X3VEC(pprime, rec->rec_SoR, work);
847 
848  switch (hitp->hit_surfno) {
849  case REC_NORM_BODY:
850  /* Skin. x, y coordinates define rotation. radius = 1 */
851  ratio = pprime[Y];
852  if (ratio > 1.0)
853  ratio = 1.0;
854  if (ratio < -1.0)
855  ratio = -1.0;
856  uvp->uv_u = acos(ratio) * M_1_2PI;
857  uvp->uv_v = pprime[Z]; /* height */
858  break;
859  case REC_NORM_TOP:
860  /* top plate */
861  len = sqrt(pprime[X]*pprime[X]+pprime[Y]*pprime[Y]);
862  ratio = pprime[Y]/len;
863  if (ratio > 1.0)
864  ratio = 1.0;
865  if (ratio < -1.0)
866  ratio = -1.0;
867  uvp->uv_u = acos(ratio) * M_1_2PI;
868  uvp->uv_v = len; /* rim v = 1 */
869  break;
870  case REC_NORM_BOT:
871  /* bottom plate */
872  len = sqrt(pprime[X]*pprime[X]+pprime[Y]*pprime[Y]);
873  ratio = pprime[Y]/len;
874  if (ratio > 1.0)
875  ratio = 1.0;
876  if (ratio < -1.0)
877  ratio = -1.0;
878  uvp->uv_u = acos(ratio) * M_1_2PI;
879  uvp->uv_v = 1 - len; /* rim v = 0 */
880  break;
881  }
882  /* Handle other half of acos() domain */
883  if (pprime[X] < 0)
884  uvp->uv_u = 1.0 - uvp->uv_u;
885 
886  if (uvp->uv_u < 0) uvp->uv_u = 0;
887  else if (uvp->uv_u > 1) uvp->uv_u = 1;
888  if (uvp->uv_v < 0) uvp->uv_v = 0;
889  else if (uvp->uv_v > 1) uvp->uv_v = 1;
890 
891  /* XXX uv_du should be relative to rotation, uv_dv relative to height */
892  uvp->uv_du = uvp->uv_dv = 0;
893 }
894 
895 
896 int
897 rt_rec_params(struct pc_pc_set *UNUSED(ps), const struct rt_db_internal *ip)
898 {
899  if (ip) RT_CK_DB_INTERNAL(ip);
900 
901  return 0; /* OK */
902 }
903 
904 
905 void
906 rt_rec_free(struct soltab *stp)
907 {
908  struct rec_specific *rec =
909  (struct rec_specific *)stp->st_specific;
910 
911  BU_PUT(rec, struct rec_specific);
912 }
913 
914 
915 /* plot and tess are handled by g_tgc.c */
916 /* import, export, ifree, and describe are also handled by g_tgc.c */
917 
918 /** @} */
919 /*
920  * Local Variables:
921  * mode: C
922  * tab-width: 8
923  * indent-tabs-mode: t
924  * c-file-style: "stroustrup"
925  * End:
926  * ex: shiftwidth=4 tabstop=8
927  */
uint32_t hit_magic
Definition: raytrace.h:249
#define RT_LEN_TOL
Definition: raytrace.h:169
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define BU_LIST_INSERT(old, new)
Definition: list.h:183
vect_t rec_A
Definition: rec.c:152
struct hit seg_in
IN information.
Definition: raytrace.h:370
int(* ft_bbox)(struct rt_db_internal *, point_t *, point_t *, const struct bn_tol *)
Definition: raytrace.h:2196
fastf_t rec_iBsq
Definition: rec.c:155
void rt_rec_free(struct soltab *stp)
Definition: rec.c:906
#define RT_DOT_TOL
Definition: raytrace.h:170
#define SMALL
Definition: defines.h:351
vect_t rec_V
Definition: rec.c:148
#define RT_CK_APPLICATION(_p)
Definition: raytrace.h:1675
#define RT_CK_RTI(_p)
Definition: raytrace.h:1833
vect_t rec_Hunit
Definition: rec.c:149
int rt_rec_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: rec.c:275
double dist
>= 0
Definition: tol.h:73
vect_t crv_pdir
Principle direction.
Definition: raytrace.h:307
fastf_t uv_u
Range 0..1.
Definition: raytrace.h:341
void rt_rec_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
Definition: rec.c:830
struct soltab * seg_stp
pointer back to soltab
Definition: raytrace.h:372
#define VSET(a, b, c, d)
Definition: color.c:53
Definition: raytrace.h:215
mat_t rec_SoR
Definition: rec.c:150
Definition: pc.h:108
Definition: raytrace.h:368
Definition: raytrace.h:248
void rt_rec_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
Definition: rec.c:790
#define SMALL_FASTF
Definition: defines.h:342
fastf_t st_aradius
Radius of APPROXIMATING sphere.
Definition: raytrace.h:433
Header file for the BRL-CAD common definitions.
struct resource * a_resource
dynamic memory resources
Definition: raytrace.h:1591
fastf_t rec_iAsq
Definition: rec.c:154
struct bu_list l
Definition: raytrace.h:369
#define REC_NORM_BODY
Definition: rec.c:412
if(share_geom)
Definition: nmg_mod.c:3829
void rt_rec_print(const struct soltab *stp)
Definition: rec.c:399
Definition: color.c:49
struct rt_i * a_rt_i
this librt instance
Definition: raytrace.h:1588
#define RT_G_DEBUG
Definition: raytrace.h:1718
vect_t hit_vpriv
PRIVATE vector for xxx_*()
Definition: raytrace.h:253
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
mat_t rec_invRoS
Definition: rec.c:151
fastf_t st_bradius
Radius of BOUNDING sphere.
Definition: raytrace.h:434
fastf_t crv_c2
curvature in other direction
Definition: raytrace.h:309
fastf_t uv_dv
delta in v
Definition: raytrace.h:344
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
#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
struct hit seg_out
OUT information.
Definition: raytrace.h:371
#define REC_NORM_BOT
Definition: rec.c:414
void bn_mat_trn(mat_t om, const mat_t im)
#define UNUSED(parameter)
Definition: common.h:239
void rt_rec_vshot(struct soltab **stp, struct xray **rp, struct seg *segp, int n, struct application *ap)
Definition: rec.c:633
#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
vect_t rec_B
Definition: rec.c:153
vect_t r_dir
Direction of ray (UNIT Length)
Definition: raytrace.h:219
#define ID_REC
Right Elliptical Cylinder [TGC special].
Definition: raytrace.h:465
struct bn_tol rti_tol
Math tolerances for this model.
Definition: raytrace.h:1765
#define RT_GET_SEG(p, res)
Definition: raytrace.h:379
#define ZERO(val)
Definition: units.c:38
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
int rt_rec_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *tol)
Definition: rec.c:162
const struct rt_functab OBJ[]
Definition: table.c:159
void bn_vec_ortho(vect_t out, const vect_t in)
#define R
Definition: msr.c:55
#define RT_CK_SOLTAB(_p)
Definition: raytrace.h:453
void * st_specific
-> ID-specific (private) struct
Definition: raytrace.h:435
fastf_t uv_du
delta in u
Definition: raytrace.h:343
fastf_t crv_c1
curvature in principle dir
Definition: raytrace.h:308
void rt_rec_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
Definition: rec.c:756
Definition: color.c:51
int hit_surfno
solid-specific surface indicator
Definition: raytrace.h:255
#define REC_NORM_TOP
Definition: rec.c:413
vect_t hit_normal
DEPRECATED: Surface Normal at hit_point, use RT_HIT_NORMAL.
Definition: raytrace.h:252
#define RT_HIT_INIT_ZERO
Definition: raytrace.h:260
const struct rt_functab * st_meth
pointer to per-solid methods
Definition: raytrace.h:428
int st_id
Solid ident.
Definition: raytrace.h:431
#define RT_REC_SEG_MISS(SEG)
Definition: rec.c:628
int rt_rec_params(struct pc_pc_set *ps, const struct rt_db_internal *ip)
Definition: rec.c:897
fastf_t hit_dist
dist from r_pt to hit_point
Definition: raytrace.h:250
int rt_rec_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
Definition: rec.c:428
double fastf_t
Definition: defines.h:300
#define VPRINT(a, b)
Definition: raytrace.h:1881
Definition: color.c:50
fastf_t uv_v
Range 0..1.
Definition: raytrace.h:342
point_t st_center
Centroid of solid.
Definition: raytrace.h:432
#define RT_HIT_MAGIC
Definition: magic.h:161
#define UNLIKELY(expression)
Definition: common.h:282
#define DEBUG_ARB8
8 Print voluminous ARB8 details
Definition: raytrace.h:91