diff options
| -rw-r--r-- | libdimension/torus.c | 54 |
1 files changed, 51 insertions, 3 deletions
diff --git a/libdimension/torus.c b/libdimension/torus.c index 8e53f99..df78212 100644 --- a/libdimension/torus.c +++ b/libdimension/torus.c @@ -25,7 +25,6 @@ */ /* Torus object callbacks */ - static bool dmnsn_torus_intersection_fn(const dmnsn_object *torus, dmnsn_line line, dmnsn_intersection *intersection); @@ -59,6 +58,51 @@ dmnsn_new_torus(double major, double minor) return torus; } +/* Bound the torus in a cylindrical shell */ +static inline bool +dmnsn_torus_bound_intersection(const dmnsn_torus_payload *payload, dmnsn_line l) +{ + double R = payload->major, r = payload->minor; + double rmax = R + r, rmin = R - r; + double rmax2 = rmax*rmax, rmin2 = rmin*rmin; + + /* Try the caps first */ + double tlower = (-r - l.x0.y)/l.n.y; + double tupper = (+r - l.x0.y)/l.n.y; + dmnsn_vector lower = dmnsn_line_point(l, tlower); + dmnsn_vector upper = dmnsn_line_point(l, tupper); + double ldist2 = lower.x*lower.x + lower.z*lower.z; + double udist2 = upper.x*upper.x + upper.z*upper.z; + if ((ldist2 < rmin2 || ldist2 > rmax2) && (udist2 < rmin2 || udist2 > rmax2)) + { + /* No valid intersection with the caps, try the cylinder walls */ + double dist2 = l.x0.x*l.x0.x + l.x0.z*l.x0.z; + double bigcyl[3], smallcyl[3]; + bigcyl[2] = smallcyl[2] = l.n.x*l.n.x + l.n.z*l.n.z; + bigcyl[1] = smallcyl[1] = 2.0*(l.n.x*l.x0.x + l.n.z*l.x0.z); + bigcyl[0] = dist2 - rmax2; + smallcyl[0] = dist2 - rmin2; + + double x[4]; + size_t n = dmnsn_solve_polynomial(bigcyl, 2, x); + n += dmnsn_solve_polynomial(smallcyl, 2, x + n); + + size_t i; + for (i = 0; i < n; ++i) { + dmnsn_vector p = dmnsn_line_point(l, x[i]); + if (p.y >= -r && p.y <= r) + break; + } + + if (i == n) { + /* No valid intersection found */ + return false; + } + } + + return true; +} + /* Returns the closest intersection of `line' with `torus' */ static bool dmnsn_torus_intersection_fn(const dmnsn_object *torus, dmnsn_line line, @@ -66,9 +110,13 @@ dmnsn_torus_intersection_fn(const dmnsn_object *torus, dmnsn_line line, { dmnsn_line l = dmnsn_transform_line(torus->trans_inv, line); const dmnsn_torus_payload *payload = torus->ptr; - double R2 = payload->major*payload->major, r2 = payload->minor*payload->minor; + double R = payload->major, r = payload->minor; + double R2 = R*R, r2 = r*r; + + if (!dmnsn_torus_bound_intersection(payload, l)) + return false; - /* This bit of algebra here is correct, and can be verified with a CAS */ + /* This bit of algebra here is correct */ dmnsn_vector x0mod = dmnsn_new_vector(l.x0.x, -l.x0.y, l.x0.z); dmnsn_vector nmod = dmnsn_new_vector(l.n.x, -l.n.y, l.n.z); double nn = dmnsn_vector_dot(l.n, l.n); |