## The offset ellipse

Thursday, December 7th, 2006 at 2:29 pm I have made this lovely drawing of an offset ellipse in SVG which I do declare is in the public domain. Click to download.

The geometry of the offset ellipse is fundamental to many CAM operations that involve the bull-nosed/toroidal cutter. I’ll explain the applications at some later time. When your CAM system is generating cutter paths, it may be spending up to 20% of its time working on the calculation illustrated here. It would be a great resource if someone were to make a public domain version of this calculation which we could all rely on, in the same way that we rely on other sophisticated calculations, such as the cosine or sqrt.

In the diagram, the ellipse is aligned along the xy axes, with a major radius a and minor radius b. Points on the ellipse have coordinates (au, bv) where u2 + v2 = 1.

The normal to the curve at this point is parallel to (bu, av). You can prove this by setting v = (1 – u2)0.5, differentiating with respect to u to get a tangent vector that is perpendicular to this normal direction.

While the direction of the normal vector is correct, its magnitude changes for different points along the ellipse. The curve we want, which is an offset by a fixed distance t, is defined by normalizing (changing its length to become a unit) this vector.

Let the length of our normal vector be

n = |(bu, av)| = ((bu)2 + (av)2)0.5

Then a point on the offset ellipse curve has coordinates

(au + but/n, bv + avt/n)

It’s easy to plot this curve, as I have done, by looping through successive values of (u, v), say by using (cos z, sin z), find n and apply the transformation. But what we really need to have written down in some complex subroutine is a function called E(k) which, for any x = k, provides for the point (k, E(k)) to be on this offset ellipse curve.

To solve for this, we need to find u such that k =au + but/n. Then the value of E(k) = bv + avt/n is easy to determine.

The equation simplifies to (k – au)n = but. We know that:

n2 = b2u2 + a2v2 = b2u2 + a2(1 – u2) = (b2 – a2)u2 + a2

so

(k – au)2((b2 – a2)u2 + a2) = b2u2t2

is the equation we must solve to get the answer.

In the past I used the approach of multiplying everything out into the great big quartic equation as follows:

a2(b2 – a2)u4 – 2ak(b2 – a2)u3 + (a4 – k2(b2 – a2) – b2t2)u2 – 2a3ku + a2k2 = 0

and solve it by applying the standard methods. I’m not so convinced with this way of doing things these days, because the value of v = (1 – u2)0.5 becomes unstable near the edge when u=1. I favour numeric methods these days which respect the geometry, rather than projecting everything into one dimension that happens to have an analytic solution.

You’ll notice that the u4 and u3 terms vanish when a=b and the ellipse is a perfect circle, and more amenable to calculation, because the offset curve of a circle is also a circle. As you can see, the offset of an ellipse is not an ellipse.

I’ll update any corrections to these equations, because I am sure I’ve typed in a mistake somewhere.

• 1. Cheng-Chang Wu replies at 9th December 2006, 4:44 pm :

Hi,

what do you mean by “numerical methods which respect the geometry”? Could you give an example?

• 2. Julian replies at 9th December 2006, 9:47 pm :

Sure. You can use the geometry and the layout of the curves in 2D to make a good first starting point from which to apply the Newton-Rhapson method. The importance of this is easy to underestimate.

For an example, imagine you are finding the intersection between the unit circle and a line. You can express the circle in the form (u, sqrt(1 – u2)). The N-R method takes an approximation at u0, and returns a value d0, so that the next and closer approximation is u1 = (u0 + d0). This will get very sensitive when u is close to 1 where there is a singularity in this expression.

But there are no singularities on the circle, from a geometric point of view. You can express your circle in the form (u, v) where u2+v2=1, and get the N-R method to work with actual tangent vectors, giving a vector (d0, e0), so that your next approximation will be (u0 + d0, v0 + e0). This point won’t be on the unit circle, so you’ll have to normalize it, but the calculation won’t have any problems at u=1.

• 3. Anders replies at 14th January 2007, 9:43 pm :

Julian, Martin: is this still a ‘current’ problem that you are interested in ?

I solved it fairly easily using the parametrization
x=a*cos(t)
y=b*sin(t)

my E(k) looks for which normal ends exactly on x=k. I din’t even look seriously for an analytic solution to the resulting equation, since I understood from your text that numerical methods are ‘allowed’!

see http://imagebin.org/6913 (link valid for 7 days)

let me know if you are interested in the code.

Anders

• 4. Julian replies at 15th January 2007, 6:58 pm :

If you’ve found it easy, then it can’t be right. Are you sure your normal has length exactly t so you’re not just working with a bigger ellipse?

• 5. Anders replies at 16th January 2007, 9:51 am :

Hi again Julian, define “easy” 🙂

math markup is not very nice using only ascii, but here goes:

points on the ellipse lie on
e_x = a*cos(t)
e_y = b*sin(t))
where t is in [0,pi/2] (since it’s symmetrical, we don’t have to worry about the other quadrants)

for each t value, the normal has a direction:

(b*cos(t), a*sin(t))

but as you noted earlier, this needs to be normalized, and then multiplied by T to get the correct offset:
so the coordinates of the sought normal vector are:

n_x = (b*cos(t)/l)*T
n_y = (a*sin(t)/l)*T

where l is the length of the normal

l = sqtr( ( b*cost(t) )^2 + ( a*sin(t) )^2 )

now we can exress where the endpoint of each normal lies:

N_x(t) = e_x + n_x = (a+T*b/l)*cos(t)
N_y(t) = e_y + n_y = (b+T*a/l)*sin(t)

and our E(k) needs to find which of these normals end at x=k, i.e. for what t value, does N_x = k hold, or

(a+T*b/l)*cos(t) = k

now this would be really easy if l was some simple expression, but it’s not, so we need to insert it :
(I’ve also moved over k)

(a+T*b/sqtr( ( b*cost(t) )^2 + ( a*sin(t) )^2 ))*cos(t) -k = 0

using numerical methods, the solution tk to this equation is quite easily found. I used MATLAB’s function fzero (google tells me that it “uses a combination of bisection, secant, and inverse quadratic interpolation methods.”)

then ofcourse the sought point on the offset ellipse is

N_x(tk)
N_y(tk)

hope this helps and regards,

Anders

• 6. Freesteel » Blog Ar&hellip replies at 14th February 2007, 11:36 am :

[…] t the raw SVG file, which you are now previewing. I did this to make the diagram for my previous posting on the offset ellipse, but then lost the program. This time I&#82 […]

• 7. Anders replies at 16th February 2007, 6:24 pm :

A diagram with perhaps a better explanation of my idea for finding the offset ellipse is now at
http://www.anderswallin.net/2007/02/offset-ellipse/

• 8. Freesteel&hellip replies at 3rd October 2008, 1:05 pm :

[…] for the last thirty years.” Let’s not be too hasty and get beyond into stuff involving tori, and triangulations, and hundreds of other things it would be useful for someone to look at. You […]

• 9. anderswallin.net ›&hellip replies at 17th December 2011, 11:41 am :

[…] Again the edge-test is the most difficult, and I never really understood where the geometry of the offset-ellipse shows up… anyone care to explain? This was written by admin. Posted on Sunday, March 21, 2010, […]