Freesteel Blog » Line intersecting a finite cylinder

Line intersecting a finite cylinder

Wednesday, October 31st, 2012 at 12:07 am Written by:

This is part one of a two part posting that deals with one component of the 5-axis drop cutter function. (It’s also the programming problem given to me when I applied for a job at Parasolid in about 1995.)

There’s a big problem with the finite cylinder because you need to deal with 3 separate surfaces.
As usual with mathematics, there’s a lot of setting up for the equations.

Point in cylinder

Consider the cylinder C of radius r whose axis goes from a to b, and the point p.


  v = b - a
  vsq = v . v
  m = (p - a) . v / vsq
      ==> (a + v m - p) . v == 0  [the vectors are perpendicular]
  d = | a + v m - p |

Then p is inside the cylinder when

  d <= r  and
  0 <= m <= 1
Line intersect with cylinder

Now take the line L that goes from p to q. Which values of s will the interpolated point (p (1 – s) + q s) on L be inside the cylinder C?

Because the cylinder is convex, the intersection is either empty or defined by the interval
s0 <= s <= s1

Without loss of generality, assume that a = 0, (so that b=v) and the line L is parallel to the Z-axis, so that points on it are of the form

  p = (qx, qy, z) = q + (0, 0, z)

We want to calculate the values of entry and exit, z0 <= z1.

The cylinder considered as independent surfaces

This method requires that we find the set Z of all points of intersection across each of the three surfaces (cylindrical surface, and two flat end disks) and, if Z is not empty, then z0 = min(Z), and z1 = max(Z).

At an intersection with the cylindrical surface:

    r2 = (v m - p)2 
       = (v m - p) . v m - (v m - p) . p  [first term is zero]
    r2 + v . p m - p . p = 0
      v . p = v . q + vz z
      p . p = q . q + z2
      m = p . v / vsq = (v . q + vz z) / vsq

    As an equation in z:

    0 = r2 vsq + (v . q + vz z)^2 - p . p vsq
      = z2 (vz2 - vsq) + 2 z (v . q vz) + (v . q)2 + r2 vsq - q . q vsq 

    This is a quadratic equation with either 0 or 2 solutions. For each solution value of z, calculate m and,
    if 0 <= m <= 1, add it to the set of intersections Z.

At an intersection with the flat end disk at a [at the origin]:

       0 = v . p
         = v . q + vz z
       za = -v . q / vz

    If this point on the plane is within the circular boundary of the disk, then:

    r2 >= p . p
        = q . q + za2
        = q . q + (v . q)2 / vz2

    in which case add it to the set of intersections Z.

At an intersection with the flat end disk at b [same as v because a is the origin]:

       0 = v . (p - v)
         = v . q + vz z - vsq
       zb = vsq - v . q / vz

    If this point is on the disk, then:

      r2 >= (p - v)2

    in which case add it to the set of intersections Z.

As already stated: z0 = min(Z), and z1 = max(Z). We expect Z to contain either 0 or 2 values.

This algorithm works for convex volumes bounded by surfaces. However, it suffers from the obvious problem of cracks. What happens if the line intersects exactly through the edge between the cylindrical surface and one of the flat disks at the ends? Well, since floating point calculations on such a boundary can go either way, there is a 25% chance of accepting an intersection from both surfaces, and a 25% change of missing both of them.

A surplus intersection isn’t a problem because it is masked by the min() and max() functions, but missing the intersections completely means we return a measurement of the cylinder that is paper thin, missing its bulky volume and probably causing a collision to result from an algorithm that depends on this function down the line.

How likely is it that the line will hit this edge? Very likely, because geometrical assemblies tend to line these edges up so that the worst always happens.

The epsilon fix

If it’s okay to over-estimate the range of the line that intersects with the cylinder, then all we need to do is extend the surfaces slightly by small value, such as 0.00001, and rewrite our three inequalities as:

 -0.00001 <= m <= 1.00001
(r + 0.00001)2 >= p . p
(r + 0.00001)2 >= (p - v)2

thus papering over any possible gap, so the bug is fixed and we break for lunch satisfied.

Unfortunately, this leads to code-rot by littering the system with fixed values that just felt right at the time, and which now limits the scale at which this function will be reliable. And what happens if v is close enough to horizontal or vertical such that the crack opens up again? If the bug recurs, your only option will be to increase the epsilon value from 0.00001, or scale it in a specific case
(eg setting it to 0.001 when |vx|+|vy|<0.001).

Pretty soon, everything is kludged, and the results are wider than they should be, and subsequent code becomes dependent on this over-estimate to cover for its flaws.

There’s no way back once you let the epsilon fix in.

The finite cylinder expressed as an intersection of two volumes

This is a small digression. The finite cylinder can be defined as the intersection volume between the infinite cylinder and the slab volume between the two planes through a and b perpendicular to v.

Take the interval [z, z+] from the quadratic equation for the intersections with the infinite cylinder, and intersect it with the interval [za, zb] from intersecting with the two planes (assuming vz > 0), and that’s the value.

This is much better, but it doesn’t necessarily generalize (see part 2 when I write it), and does not lend to optimization when we need to find the answer for multiple parallel lines at the same time.

Treating the line as a point within a 2D shape

Let’s look at that diagram again.

By projecting the geometry down the Z-axis, the problem becomes one of locating a 2D point within a shape made up of two ellipses and a rectangle between them.

This dimensional reduction is what makes the calculation work. Next we need to reduce the problem to 1 dimension, which is the only level that computational geometry operates.

First, let f be the flattened unit vector of v (assuming vz >= 0):

   f = (vx, vy) / |(vx, vy)|
   e = vz / |v|  [the eccentricity of the ellipses] 
   d = f . v

   m = f . p
   s = fT . p   [the 2D perpendicular of f]

We know that we completely miss the cylinder if:

   |s| > r   or
   m < r e   or
   m > d + r e

So, subject to the precalculations of f, e and d, we can filter out most cases when the line misses the cylinder with two dot-products. That’s very efficient.

And, given the value of s, we can see that the red line intersects the two ellipses at four different points
{-k, k, d-k, d+k} where k = e sqrt(r2 – s2).

From the position of m in between these four values, we can determine exactly which components of the cylinder with which we must find an intersection with the line, and produce an optimal solution that will always be reliable and does not depend on any epsilon fudge factors.

In part two, the line gets replaced by a rectangular volume in a problem that took me four complete attempts (of at least 600 lines of wasted code each) to get it right. The code is not intended to be looked at again, because it is as fast as it can be and should not have any failure cases, until it becomes necessary to extend it by the cylinder becoming a truncated cone.

1 Comment

  • 1. Line intersecting a finit&hellip replies at 31st October 2012, 8:59 am :

    […] more here: Line intersecting a finite cylinder – Freesteel This entry was posted in Cylinder Scales and tagged a-big-problem-, because-you, equations, […]

Leave a comment

XHTML: You can use these tags: <a href="" title=""> <abbr title=""> <acronym title=""> <blockquote cite=""> <code> <em> <strong>