summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--libdimension/torus.c54
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);