A circle or (its more general form) an ellipse is often regarded as a simpler shape than a cubic polynomial and so may be preferable to cubic splines.  Elliptical arcs are called conic splines and are popular in applications such as the generation of font outlines, where speed is critical.

Here we’ll investigate an algorithm for generating points along an elliptical

arc. The points are separated by a fixed angular increment specified in

radians of elliptical arc. The algorithm is based on a parametric representation of the ellipse. It is particularly inexpensive in terms of the amount of computation required. Only a few multiplications (shifts, in binary form), additions, and subtractions are needed to generate each pointwithout approximation.

A Conic Spline Drawing Algorithm

The idea is that we know two points, say P and Q, on a curve and the slopes of the tangent lines at those two points, then there is a unique ellipse determined by this information.

Exercise: Justify that statement and find the equation for the ellipse.

SOLN: A circle will do the job.  Suppose we have points P(x₁,y₁) and Q(x₂,y₂) and tangent lines with slopes m₁ and m₂ at P and Q, respectively. Construct perpendiculars

and compute the point of intersection as the center of the circle:

and the radius is the distance from this point to either one of the points of tangency:

And then plug these into .

Note that an ellipse has more degree of freedom than a circle, so there are likely many different ellipses that satisfy the initial conditions.   

The plan is to create the elliptical arc as a series of  pixels.  

Suppose the ellipse is centered at coordinates (h,k) and that ux and vx  are the displacements in x of each point from the center; uy and vy the corresponding y displacements.

The parameterization of any ellipse centered at (0,0) can be written as

                                                                                 (1.1)

The shape of the ellipse is determined by the phase and amplitude relationships of the two sinusoids, where A and B are the amplitudes, and x and y  are the phase angles. The full ellipse is formed as independent parameter t increases from 0 to 2π radians.

.

Exercise:  Show that x and y represent coordinates on an ellipse.

SOLN: Taking a right ellipse of the form .  Applying a rotation angle θ yields

Equating this with the trig expansion of parameterization (1.1) yields

Equating coefficients, we have

Equating ratios of these we have

The product of these is then

Also, equating sums of squares yields

Also, we can use identities to express various sums of products such as

…and similarly,

so that

What to do about this lingering CD product?  How about substituting ?

...which is an equation we could solve for θ.

Exercise: Find the angle of rotation for

and the lengths of the major and minor axes.

Even for a simple circle, explicitly calculating the sine and cosine terms can be quite time-consuming. A sequence of points on a circle centered at the origin can be computed by applying a (small) rotation angle Δt  to a point (x_n-1, y_n-1) to get the next point (x_n, y_n) through the transformation represented by the pair of equations:

Now if Δt is very small, then  and , so we have the approximations:

                                                                                     (1.2)

Given the coordinates (x_n-1, y_n-1) of a point on a circle, these expressions calculate a new point (x_n, y_n) that is also on the (approximate) circle. Beginning at the initial (x₀, y₀), these equations are evaluated iteratively to generate a series of points that lie roughly along a circular arc. Successive values of x_n and y_n approximate two sinusoids that have a relative phase difference of π/2. 

In terms of initial coordinates (x₀, y₀), the effect of n such iterations can be represented as

Now, suppose that we tweak our approximate iterator to use the freshly computed x_n to compute the new y_n, instead of the old one:

                                                                                     (1.3)

That means that   so that we can reformulate the matrix iteration like so:

The gemstone in this study is the idea that that iteration is equivalent to this one:

Exercise:  Show that the determinant of the transformation matrix is 1.

SOLN:   

Exercise:  Simplify a formula for the nth power of the transformation matrix.

SOLN:

So look at the first component, it works out (with little help from trigonometric identities):  

The off-diagonal element,  are straightforward and the other diagonal element works just like the first.  The higher powers work the same way, so that

 

Now, taking  we have

 Plugging these values back into the basic we have
  

This is precisely what we got for swapping x_n-1 for x_n in the simple rotational iterative scheme.

So we can iterate this well-behaved rotational transformation (elliptical?) without computing a sine or cosine or even needing a look-up table.

The transformation we have can be written as a pair of iterative equations:

Whereas Δt is the approximate angle increment in this iteration, the precise increment is

Here Δt is the chord on the unit circle corresponding to arc α.

Now, if we reposition the intial value by setting  and replacing all values of x₀ in the iterations with this we get

or, simplifying

                       (1.4)

Each iterative formula is then a sinusoid (the sum or difference of two sinusoids with the same frequency is another sinusoid with the same frequence, but a different amplitude.)  Together they generate an ellipse that can be calculated precisely by means of the very simple (no sines to evaluate) iteration represented by Eq. (1.3).  A separate copy of the iteration formulas is needed for each of the two dimensions; were the arc three dimensional, then three copies would be needed. The error in the

approximation can be entirely eliminated by modifying only the initial values. The inner loop calculation itself remains unmodified and the determinant of the corresponding rotation matrix is unity.

Computing χ₀ involves a nasty square root, so it may be better to use a binomial series approximation:

Conjugate Diameters

The iterative formulas of equations (1.1) can be split up into sine and cosine components:

                                                                            (1.5)

Exercise: Verify this claim.

SOLN:  so  and .  Similarly for y(t).  Also, .  

For an ellipse centered at the origin, the geometric interpretation of the parameterization in Eq. 1.5 is that the point on the ellipse corresponding to t = 0 is at P = (AP, BP); the point corresponding to t = p/2 is at Q = (AQ, BQ).  As t increments from 0 to p/2, the arc is a blend of P and Q, with the respective contributions at each point determined by the cosine and sine functions. The two diameters formed by extending lines from points P and Q through the center to the other side of the ellipse are referred to as conjugate diameters of the ellipse.  This is equivalent to defining the ellipse in terms of an enclosing parallelogram, as shown below

The core idea here is that an ellipse inscribed in the parallelogram is an affine transformation of a unit circle inscribed in a square. An affine transformation (a transformation that preserves lines and parallelism) transforms the circle into an ellipse and the enclosing square into a parallelogram.  An affine transformation is a combination of shears and rotations.

The ellipse's two diameters (shown as dotted lines with endpoints including P and Q) are not necessarily perpendicular (the affine transformation doesn’t preserve angles), but still connect the midpoints of opposing sides and divide the larger parallelogram into four congruent smaller ones, each containing one-quarter of the ellipse. The conic spline can be described in terms triangle PQK.    These points also define a parallelogram whose fourth vertex J is is the center of the ellipse. The coordinates of the fourth vertex, J, can be calculated as

This is needed to translate the origin of the ellipse to wherever we want to put it.

The Code

Here’s pseudocode to carry out this algorithm.

procedure ConSpline (x_p ,y_p ,x_Q ,y_Q ,x_k ,y_k ,Dt : real)

x_J ,y_J ,u_x ,v_x ,u_y ,v_y : real ;

i, n, x, y : integer :

begin

vx ← x_K - x_Q :

ux ← x_K - x_P :

vy ← y_K - y_Q :

uy ← y_K - y_P :

xJ ← x_P - v_x :

yJ ← y_P - v_y :

ux ← u_x*sqrt(10.25*Dt^2)  0.5*v_xDt ;

uy ← uy*sqrt(10.25*Dt^2)  0.5*v_yDt ;

n ← floor(π/Dt) :

for i←0 to n do

DrawPoint (x←round(vx + xJ ), y←round(vy + yJ )) ;

ux ←ux - vxDt ;

vx ←vx - uxDt ;

uy ←uy - vyDt ;

vy ←vy - uyDt ;

endloop

end

Exercise:  Write code to employ this algorithm.

SOLN: Here.

Here’s some sample output from one person’s attempt to implement this algorithm.