AA Fast Parametric Ellipse Algorithm
Jerry R. Van AkenOctober 21, 2020
Abstract
This paper describes a 2-D graphics algorithm that uses shifts andadds to precisely plot a series of points on an ellipse of any shape and ori-entation. The algorithm can also plot an elliptic arc that starts and endsat arbitrary angles. The ellipse algorithm described here is largely basedon earlier papers by Van Aken and Simar [1, 2]. A new flatness test is pre-sented for automatically controlling the spacing between points plotted onan ellipse or elliptic arc. Most of the calculations performed by the ellipsealgorithm and flatness test use fixed-point addition and shift operations,and thus are well-suited to run on less-powerful processors. C++ sourcecode listings are included.
Keywords: parametric ellipse algorithm, rotated ellipse, Minsky circle al-gorithm, flatness, elliptic arc, conjugate diameters, affine invariance
This paper describes a 2-D graphics algorithm that uses fixed-point additionand shift operations to precisely plot a series of points on an ellipse of anyshape and orientation. The ellipse algorithm is largely based on earlier papersby Van Aken and Simar [1, 2], which extend Marvin Minsky’s well-known circlealgorithm [3, 4, 5] to ellipses, and show how to cancel out the sources of errorin Minsky’s original algorithm. The ShapeGen graphics library [6] uses theresulting ellipse algorithm to generate circles, ellipses, elliptic arcs, and ellipticsplines.For the convenience of the graphics library user, an ellipse can be specifiedin terms of the parallelogram in which it is inscribed. The center of the ellipsecoincides with the center of the enclosing parallelogram. The ellipse touches(and is tangent to) the parallelogram at the midpoint of each of its four sides.In fact, the ellipse algorithm described here requires only three points tocompletely specify the ellipse: the center point, and the end points of a pairof conjugate diameters of the ellipse [1, 2, 7]. As shown in Figure 1, the twoconjugate diameter end points, labeled P and Q , are simply the midpoints oftwo adjacent sides of the enclosing parallelogram.In some applications, the direction of rotation in which the ellipse or arc isdrawn is important. The direction in which the ellipse algorithm plots points is a r X i v : . [ c s . G R ] O c t igure 1: The end points P and Q of a pair of conjugate diameters of an ellipseare simply the midpoints of two adjacent sides of the enclosing square, rectangle,or parallelogram.from the first conjugate diameter end point, P , toward the second, Q . In thecase of an elliptic arc, the arc’s starting and ending angles are measured relativeto P , and are positive in the direction of Q .The ShapeGen graphics library [6] approximates a curved shape by plottinga series of points at more-or-less regular intervals along the curve, and thenconnecting each pair of adjacent points on the curve with a straight line segment.The library user controls the spacing between successive points on an ellipse orother curve by specifying a flatness parameter. The process of approximatinga curve with a series of chords or polygonal edges is called flattening [8]. Theflatness parameter is the maximum error tolerance of a flattened representation,as specified by the user, and is defined to be the largest gap, in pixels, that ispermitted between the chords and the curve segments that they represent.The smoothness of a flattened curve can be improved by decreasing the spac-ing between adjacent points, but doing so creates more points, which increasesprocessing time; it also increases memory usage if the points are stored to awaitfurther processing.To improve the tradeoff between smoothness and speed, the ellipse algorithmis designed to minimize the processing time needed to calculate each point. Thealgorithm’s inner loop requires only four additions and four right-shift operationsto calculate each new point on an ellipse centered at the x - y coordinate origin.Two more addition operations are required to translate the ellipse center to aposition other than the origin. Fixed-point arithmetic is used. The points plottedby the ellipse algorithm fall precisely on the ellipse, to within the precisionafforded by the fixed-point calculations. In this section, we will derive the parametric equations of an ellipse, given theend points, P and Q , of a pair of conjugate diameters of the ellipse. To simplifythe calculations, the ellipse is centered at the origin.Figure 2 shows that an ellipse and its bounding parallelogram are the affine-transformed images of a unit circle and its bounding square. We use u - v coor-dinates for the unit circle, and x - y coordinates for the ellipse. Points (1 ,
0) and(0 ,
1) on the unit circle are transformed to points P = ( x P , y P ) and Q = ( x Q , y Q )on the ellipse. Any pair of perpendicular diameters of the circle are conjugatediameters of the circle. An affine transformation of the circle to an ellipse mapsa pair of conjugate diameters of the circle to conjugate diameters of the ellipse.2igure 2: Affine transformation of a unit circle and its bounding square to anellipse and its bounding parallelogram.Thus, P and Q in Figure 2 are end points of conjugate diameters of the ellipse.The affine transformation of a point ( u, v ) on the unit circle to a point ( x, y )on the ellipse can be expressed as (cid:20) xy (cid:21) = M (cid:20) uv (cid:21) where M = (cid:20) m m m m (cid:21) To find the coefficients m ij of transformation matrix M , consider that conjugatediameter end points P and Q on the ellipse are related to points (1 ,
0) and (0 , P = M (cid:104) (cid:105) T and Q = M (cid:104) (cid:105) T . Thesetwo expressions expand to the following: (cid:20) x P y P (cid:21) = (cid:20) m m m m (cid:21) (cid:20) (cid:21) = (cid:20) m m (cid:21)(cid:20) x Q y Q (cid:21) = (cid:20) m m m m (cid:21) (cid:20) (cid:21) = (cid:20) m m (cid:21) From inspection, we see that M = (cid:20) x P x Q y P y Q (cid:21) The affine transformation of a point ( u, v ) on the unit circle to a point ( x, y ) onthe ellipse can now be expressed as (cid:20) xy (cid:21) = (cid:20) x P x Q y P y Q (cid:21) (cid:20) uv (cid:21) (1)where the matrix coefficients are the x and y coordinates of end points P and Q of a pair of conjugate diameters of the ellipse.3he parametric equations for the unit circle are u ( θ ) = cos θ (2) v ( θ ) = sin θ for 0 ≤ θ ≤ π . The corresponding parametric equations for the ellipse (Foley et al [9]) can be obtained by substituting equations (2) into equation (1), whichyields x ( θ ) = x P cos θ + x Q sin θ (3) y ( θ ) = y P cos θ + y Q sin θ for 0 ≤ θ ≤ π .These equations can also be expressed in the following sequential form, whichcan be used to plot a series of N points on the ellipse by increasing angle θ from α to 2 π in incremental steps of size α = 2 π/N radians: x n = x P cos( nα ) + x Q sin( nα ) (4) y n = y P cos( nα ) + y Q sin( nα ) (5)for n = 1 , , ..., N , where P = ( x P , y P ) and Q = ( x Q , y Q ) are the end pointsof a pair of conjugate diameters of the ellipse. However, directly using theseequations requires evaluating the sine and cosine terms in the inner loop of theellipse algorithm, which could be slow. To speed up the inner loop of the ellipse algorithm, we will adapt a well-knowncircle-generating algorithm that was invented by Marvin Minsky in the early1960s (Beeler et al [3], Paeth [4], Ziegler Hunts et al [5]). This circle generatorcan rotate a point about the origin using only shifts and adds. Here’s Minsky’sinner loop: u n = u n − − (cid:15) v n − (6) v n = (cid:15) u n + v n − (7)where 0 < (cid:15) ≤
1, and 1 ≤ n ≤ (cid:98) π/(cid:15) (cid:99) . We use u - v coordinates here to avoidconfusion with the x - y coordinates used for the ellipse in equations (4) and (5).The circle generator initially sets the u and v values to a starting point( u , v ) on a circle centered at the origin. Each iteration of equations (6) and(7) plots the next point on the (approximate) circle by rotating the previouspoint about the center by an (approximate) angular increment of (cid:15) radians.The circle generator draws an “almost circle”—a nearly circular ellipse—andthe ellipse becomes a better approximation to a circle as the size of angularincrement (cid:15) is decreased.The angular increment can be set to a negative power of two ( (cid:15) = 1 / k for k = 0 , , , ... ), for which the multiplications in equations (6) and (7) can be Given conjugate diameter end points P = ( x P , y P ) and Q = ( x Q , y Q ) on an ellipse, equa-tions (3) can be used to generate additional conjugate diameter end points P (cid:48) = ( x ( θ ) , y ( θ ))and Q (cid:48) = ( x ( θ ± π ) , y ( θ ± π )) on the same ellipse. Also, see equations (14) and (15). u and v values in these equations can berepresented as 16.16 fixed-point values; these values are stored as signed 32-bitintegers, but the 16 least-significant bits are assumed to lie to the right of thebinary point, and represent fractional values.The circle generator looks similar to a standard rotation that uses the ap-proximations cos (cid:15) ≈ (cid:15) ≈ (cid:15) , which are useful approximations if theangular increment (cid:15) is small. However, the calculation of v n in equation (7) uses u n rather than u n − . Although this subscript value might look like a typingerror, it is actually key to the algorithm’s behavior (Newman and Sproull [10],Blinn [11]). The result is that the determinant of equations (6) and (7) is unity.This fact can be verified by substituting the right-hand side of equation (6) inplace of u n in equation (7) and then expressing the equations in matrix form,as follows: (cid:20) u n v n (cid:21) = (cid:20) − (cid:15)(cid:15) − (cid:15) (cid:21) (cid:20) u n − v n − (cid:21) Because the determinant is unity, the curve plotted by Minsky’s circle generatorcloses as it completes a full rotation. In fact, the circle generator is remarkablystable—if the inner loop is allowed to free-run for many rotations, it will con-tinue to plot points on the same approximate circle without spiraling inward oroutward.The major drawback to using Minsky’s circle generator is that it draws onlyapproximate circles, whereas we want an algorithm to draw ellipses precisely.To fix this problem, we will identify the sources of error in the circle generatorand cancel them out.Analyses by Van Aken and Simar [1, 2] and Ziegler Hunts et al [5] found thatthe result of n iterations of equations (6) and (7) is as follows: u n = u cos( nα ) − v − (cid:15) u (cid:113) − (cid:15) sin( nα ) (8) v n = u − (cid:15) v (cid:113) − (cid:15) sin( nα ) + v cos( nα ) (9)where ( u , v ) is the starting point on the circle, and α is the precise angular in-crement, to which the value (cid:15) in equations (6) and (7) is only an approximation.To draw a perfect circle, the messy-looking terms in brackets in equations (8)and (9) should equal v and u , respectively. As expected, however, the termsin brackets do approach the desired values, v and u , as (cid:15) decreases in size.The angular increment α in equations (8) and (9) is related to the approxi-mate angular increment (cid:15) as follows [1, 2, 5]:sin α (cid:15) α (cid:113) − (cid:15) v or u . For example, the error inequation (9) can be eliminated by replacing u with a value U specified so that u = U − (cid:15) v (cid:113) − (cid:15) Solving for U , we have U = u (cid:113) − (cid:15) + (cid:15) v (10)Substituting U for u in equations (8) and (9) yields u n = (cid:20) u (cid:113) − (cid:15) + (cid:15) v (cid:21) cos( nα ) − (cid:20) v (cid:113) − (cid:15) − (cid:15) u (cid:21) sin( nα ) (11) v n = v cos( nα ) + u sin( nα ) (12)Equation (12) shows that we can turn Minsky’s approximate circle generatorinto a precise sine-wave generator simply by replacing the initial value u withthe modified value U . By comparing equation (12) with equations (4) and (5), we see that the el-lipse algorithm will need to incorporate two copies of the circle generator. Thefirst copy will produce x n in equation (4), and the second will produce y n inequation (5).The following C++ function implements the circle generator in equations(6) and (7), and is called twice in the ellipse algorithm’s inner loop: inline void CircleGen(FIXED& u, FIXED& v, int k){ u -= v >> k;v += u >> k;} The parameters of type
FIXED in this listing are 16.16 fixed-point values; aspreviously discussed, these values are stored as signed 32-bit integers, but the16 least-significant bits are assumed to lie to the right of the binary point, andrepresent fractional values. Parameters u and v represent the values u and v inequations (6) and (7). The angular increment (cid:15) in these equations is set to anegative power of two, (cid:15) = 1 / k for k = 0 , , , ... , where exponent k is specifiedby parameter k . The inline qualifier means that the CircleGen function incursno function-call performance penalty.For equation (10), U can be calculated by using floating-point arithmeticand calling the sqrt function in standard C header file math.h . However, a faster This implementation of the
CircleGen function assumes that the C++ compiler supportsarithmetic right-shift, as is required by the C++20 standard. U can be expressed as U = u (cid:18) − (cid:15) − (cid:15) − (cid:15) − (cid:15) − ... (cid:19) + (cid:15) v (13)where (cid:15) = 1 / k for k = 0 , , , . . . In this form, U can be calculated using fixed-point arithmetic. The middle three terms inside the parentheses are negativepowers of two, and multiplications of these terms by u can be performed asright-shift operations. For many applications, sufficient accuracy is obtained bytruncating the Taylor series after the 6th-order term.The following C++ function uses equation (13) to calculate U , but omitsthe 8th-order and higher terms: FIXED InitialValue(FIXED u0, FIXED v0, int k){ int shift = 2*k + 3;FIXED w = u0 >> shift;FIXED U0 = u0 - w + (v0 >> (k + 1));w >>= shift + 1;U0 -= w;w >>= shift;U0 -= w;return U0;}
Parameters u0 and v0 represent the initial values, u and v , in equation (13).As before, angular increment (cid:15) = 1 / k , where exponent k is represented byparameter k .The ellipse algorithm calls the InitialValue function twice—the first calladjusts the initial value of the x Q coordinate, and the second adjusts the initialvalue of y Q . These adjustments ensure the accuracy of the x - y coordinates thatare calculated in the algorithm’s inner loop and used to plot points on the ellipse.The following C++ function implements the core ellipse algorithm, and canplot both ellipses and elliptic arcs: void EllipseCore(FIXED xC, FIXED yC, FIXED xP, FIXED yP,FIXED xQ, FIXED yQ, FIXED sweep, int k){ int count = sweep >> (16 - k);PlotPoint(xP+xC, yP+yC);xQ = InitialValue(xQ, xP, k);yQ = InitialValue(yQ, yP, k);for (int i = 0; i < count; ++i){ CircleGen(xQ, xP, k);CircleGen(yQ, yP, k);PlotPoint(xP+xC, yP+yC);}} CircleGen and
InitialValue , that werepreviously discussed. It also calls an inline function,
PlotPoint , that we canassume plots a point that’s specified by its x - y coordinates (in 16.16 fixed-point format). Parameters xC and yC are the x - y coordinates at the centerof the ellipse. Parameters xP and yP specify the first conjugate diameter endpoint, P = ( x P , y P ), and xQ and yQ specify the second conjugate diameter endpoint, Q = ( x Q , y Q ). These two end points are specified with center-relative coordinates; that is, as x - y offsets from the ellipse center. Parameter sweep specifies the sweep angle, which is the angle traversed by the elliptic arc. The sweep parameter is a positive, fixed-point value, and is expressed in radians.The approximate angular increment between points plotted by the EllipseCore function is (cid:15) = 1 / k , where k is specified by parameter k .An arc plotted by the EllipseCore function always starts at point P andextends in the direction of point Q . However, a calling function can set anarbitrary arc starting point on the ellipse, and can support both positive andnegative sweep angles. To do so, the caller modifies the conjugate diameter endpoints before passing them to EllipseCore , as will be described shortly.To determine how many points to plot, the
EllipseCore function dividesthe sweep angle by the angular increment (cid:15) and truncates the fixed-point resultto an integer value. This calculation could be performed by shifting the sweep parameter value left by k bits, and then shifting the result right by 16 bits, butbecause k will always be less than 16, it’s safe to combine these two operationsinto a single right-shift by 16 − k bits. A later section will discuss limits on thesize of the k parameter in more detail.The graphics library user doesn’t directly call the EllipseCore function toplot ellipses and elliptic arcs. Instead, the user calls two public interface func-tions,
PlotEllipse and
PlotEllipticArc , which will be described next, andthese functions call
EllipseCore . In calls to these two public functions, the userspecifies window-relative or viewport-relative coordinates for the ellipse centerand the two conjugate diameter end points, P and Q . However, the PlotEllipse and
PlotEllipticArc functions convert P and Q to center-relative coordinatesbefore passing them to the EllipseCore function .The PlotEllipse function is relatively simple. Here’s the C++ implemen-tation: void PlotEllipse(FIXED xC, FIXED yC, FIXED xP, FIXED yP,FIXED xQ, FIXED yQ, int k){ EllipseCore(xC, yC, xP-xC, yP-yC, xQ-xC, yQ-yC, FIX_2PI, k);}
The sweep angle passed to the
EllipseCore function is set to the constant
FIX_2PI , which is the quantity 2 π in 16.16 fixed-point format.The PlotEllipticArc function is more complex. It plots an arc of an ellipsegiven the arc’s starting angle and its sweep angle. If P and Q are the caller-specified end points of a pair of conjugate diameters on the ellipse, and the Users frequently transform lists of points to window- or viewport-relative coordinatesbefore passing them to user-callable functions like
PlotEllipse and
PlotEllipticArc . Forconvenience, these functions translate conjugate diameter end points to center-relative coor-dinates rather than requiring the user to do so. ϕ is nonzero, the function generates a new pair ofconjugate diameter end points, P (cid:48) and Q (cid:48) , where P (cid:48) is the arc starting point,and P (cid:48) and Q (cid:48) define the same ellipse as the original end points, P and Q . The PlotEllipticArc function then passes P (cid:48) and Q (cid:48) to the EllipseCore function,which uses P (cid:48) as the arc starting point.To obtain expressions for points P (cid:48) and Q (cid:48) , consider that affine transforma-tions preserve the conjugate diameters property. For example, points P and Q on the ellipse in Figure 2 must lie on conjugate diameters of the ellipse becausethey are the transformed images of points (1 ,
0) and (0 ,
1) on the unit circle,which lie on conjugate diameters of the circle. To aid in our discussion, we’lllabel these points A = (1 ,
0) and B = (0 , A and B on the unit circle aresimultaneously rotated by starting angle ϕ , the rotated points, A (cid:48) and B (cid:48) , lieon diameters that are perpendicular, and thus are conjugate diameters of thecircle. If these two rotated points on the circle are then transformed to points P (cid:48) and Q (cid:48) on the ellipse, P (cid:48) and Q (cid:48) must lie on conjugate diameters of the ellipse.First, to obtain the u - v coordinates of points A (cid:48) and B (cid:48) on the unit circle,we rotate points A and B by arc starting angle ϕ : A (cid:48) = (cid:20) cos ϕ − sin ϕ sin ϕ cos ϕ (cid:21) A = (cid:20) cos ϕ sin ϕ (cid:21) B (cid:48) = (cid:20) cos ϕ − sin ϕ sin ϕ cos ϕ (cid:21) B = (cid:20) − sin ϕ cos ϕ (cid:21) Next, we use equation (1) to transform points A (cid:48) and B (cid:48) on the unit circle topoints P (cid:48) = ( x (cid:48) P , y (cid:48) P ) and Q (cid:48) = ( x (cid:48) Q , y (cid:48) Q ) on the ellipse: P (cid:48) = (cid:20) x P x Q y P y Q (cid:21) A (cid:48) Q (cid:48) = (cid:20) x P x Q y P y Q (cid:21) B (cid:48) where the matrix coefficients are the x and y coordinates of the original conju-gate diameter end points, P = ( x P , y P ) and Q = ( x Q , y Q ). Finally, we expandthese expressions to obtain the x - y coordinates of P (cid:48) and Q (cid:48) in terms of P , Q ,and the arc starting angle, ϕ : x (cid:48) P = x P cos ϕ + x Q sin ϕ (14) y (cid:48) P = y P cos ϕ + y Q sin ϕx (cid:48) Q = x Q cos ϕ − x P sin ϕ (15) y (cid:48) Q = y Q cos ϕ − y P sin ϕ Equations (14) and (15) are incorporated into the following C++ implemen-tation of the
PlotEllipticArc function:9 oid PlotEllipticArc(FIXED xC, FIXED yC,FIXED xP, FIXED yP, FIXED xQ, FIXED yQ,float astart, float asweep, int k){ float cosb, sinb;FIXED swangle;xP -= xC;yP -= yC;xQ -= xC;yQ -= yC;if (astart != 0){ // Set new conjugate diameter end points P’ and Q’float cosa = cos(astart);float sina = sin(astart);FIXED x = xP*cosa + xQ*sina;FIXED y = yP*cosa + yQ*sina;xQ = xQ*cosa - xP*sina;yQ = yQ*cosa - yP*sina;xP = x;yP = y;}// If sweep angle is negative, switch directionif (asweep < 0){ xQ = -xQ;yQ = -yQ;asweep = -asweep;}swangle = 65536*asweep;EllipseCore(xC, yC, xP, yP, xQ, yQ, swangle, k);// Plot arc end pointcosb = cos(asweep);sinb = sin(asweep);xP = xP*cosb + xQ*sinb;yP = yP*cosb + yQ*sinb;PlotPoint(xP+xC, yP+yC);}
This function uses floating-point arithmetic to calculate the arc end point—andalso the arc starting point if the specified arc starting angle is nonzero. All thepoints in between are calculated by the
EllipseCore function using fixed-pointarithmetic. The cos and sin functions called in this listing are declared in thestandard C header file math.h .On entry, the
PlotEllipticArc function converts the conjugate diameterend points, P and Q , to center-relative coordinates. Function parameter astart PlotEllipticArc that share iden-tical arc start and sweep angles.is the starting angle of the arc, as measured from the first conjugate diameterend point, P , and is positive in the direction of the second conjugate diame-ter end point, Q . Parameter asweep is the sweep angle, and is positive in thesame direction as astart . Both angles are in radians. If astart is nonzero, thefunction uses equations (14) and (15) to calculate the new conjugate diameterend points, P (cid:48) and Q (cid:48) , which are then passed to the EllipseCore function inplace of P and Q . A sweep angle that is negative and thus extends away fromend point Q = ( x Q , y Q ) is converted to a positive angle that extends toward the point ( − x Q , − y Q ) on the opposite end of the same conjugate diameter as Q .The final point plotted by EllipseCore typically falls short of the arc end pointby a fraction of an angular increment. Before returning, the
PlotEllipticArc function plots the arc end point.Strictly speaking, the
PlotEllipticArc function’s astart and asweep pa-rameters specify the arc’s start and sweep angles on the unit circle, not on theellipse . It’s fair to ask whether this is a convenient and intuitive way to specifyan arc of an ellipse.To explore this issue, Figure 3 shows a collection of pie charts—and their en-closing squares, rectangles, and parallelograms. Each arc in the figure is drawnby a PlotEllipticArc function call. If the circular pie chart at top center is Recall that equations (2) describe the rotation of a point on the unit circle in Figure 2,and that equations (3) describe the corresponding movement of the affine-transformed imageof this point on the ellipse.
PlotEllipticArc function’s interface is in accord with the user’s intuition.Arcs plotted by the
PlotEllipticArc function share an important propertywith B´ezier curves: both are affine-invariant . For a B´ezier curve, applying anyaffine transformation to the vertexes of the control polygon produces the sametransformed image as does directly transforming the points on the curve (Rogersand Adams [12], Watt and Watt [13]). Similarly, for an arc constructed bythe
PlotEllipticArc function, application of any affine transformation to theellipse center point and the two conjugate diameter end points has the sameeffect as directly transforming the points on the arc.
As previously described, the ShapeGen graphics library [6] approximates acurved shape by plotting a series of points along the curve, and then connectingeach pair of adjacent points with a straight line segment—a chord. To controlthe spacing between points, the library user specifies a flatness parameter, whichis the maximum error tolerance, in pixels, between a chord and the curve seg-ment that it represents. Thereafter, as the user draws various curves, the libraryautomatically adjusts the spacing between plotted points on each curve to meetthe flatness requirement.The user can specify a smaller flatness value to improve the smoothness ofthe plotted curve, or can specify a larger flatness value to decrease the processingtime and memory required to plot the curve.In the previous section, the
DrawEllipse and
DrawEllipticArc functionsrequire the caller to supply a parameter value k , which specifies the approximateangular increment ( (cid:15) = 1 / k for k = 0 , , ... ) between successive points plottedon the ellipse or elliptic arc. Thus, the caller must either calculate the correctangular increment in advance of each function call, or make a series of manualadjustments based on trial-and-error inspections of the results.Clearly, a more convenient option is to allow the caller to specify a globalflatness parameter in advance of any calls to these functions.To make this improvement, we remove the last call parameter, k , from the DrawEllipse and
DrawEllipticArc functions, and also from the
EllipseCore function. The
EllipseCore function is further modified to call a new function,which we will name
AngularInc , that calculates and returns an appropriatenonnegative, integer value for the k parameter. The determination of this valueis based on the global flatness parameter and on the size and shape of the ellipse The name
AngularInc is a bit misleading because this function doesn’t return angularincrement (cid:15) . Instead, it returns k = − log (cid:15) . However, multiplication of a fixed-point value x by (cid:15) is conveniently implemented in C++ as x >> k . δ for a circle and an ellipse.specified by conjugate diameter end points P and Q .Here’s the C++ declaration for the AngularInc function: int AngularInc(FIXED xP, FIXED yP, FIXED xQ, FIXED yQ);
The input parameters to this function are the center-relative x - y coordinatesof end points P and Q , and the return value is the k parameter. The functiondetermines the smallest value of k (and, thus, the largest angular increment)for which the chord-to-arc error does not violate the flatness requirement.With this modification, the EllipseCore function plots points at regularangular increments that are specified by the k value returned by the AngularInc function.In Figure 4, regular angular increments result in uniformly spaced pointsaround the circle at the top; but when the same angular increments are appliedto the ellipse below, the spacing varies around the ellipse. To account for thevariance in the chord-to-arc error around the ellipse, the angular increment mustbe made small enough to accommodate the worst-case error.For the ellipse in Figure 4, the largest chord-to-arc error, δ , occurs at theends of the major axis. The chord-to-arc error in the circle above is also δ ; thiscircle’s diameter is equal to the length of the ellipse’s major axis. We concludethat the largest chord-to-arc error for an ellipse occurs at the ends of its majoraxis, and that this error matches the chord-to-arc error of a circle whose diameterequals the length of the ellipse’s major axis. In the lower part of Figure 4, the angles labeled α are obviously not equal if the ellipseis viewed as a plane figure. They are equal, however, if the ellipse is viewed as a copy of thecircle above that has been rotated in three dimensions about a horizontal line through itscenter. δ on the unit circle.The auxiliary circle of an ellipse is the circle whose diameter matches thelength of the major axis of the ellipse, and whose center coincides with the centerof the ellipse [14]. We can refine the description of the AngularInc function tobe that it determines the smallest k parameter (and, thus, the largest angularincrement (cid:15) = 1 / k ) that satisfies the flatness requirement for the auxiliarycircle of the ellipse .How does the chord-to-arc error for a circle depend on the angular increment?Figure 5 shows the error δ for a chord of length (cid:15) on the unit circle. This figuredistinguishes the precise angular increment α from the approximate angularincrement (cid:15) = 1 / k . The distance from the midpoint of the chord to the circle’scenter is (cid:113) − (cid:15) . Thus, the chord-to-arc error is δ = 1 − (cid:113) − (cid:15) . This resultfor the unit circle is easily generalized to a circle of radius r , for which the errorscales to δ = r (cid:18) − (cid:113) − (cid:15) (cid:19) , where (cid:15) = 1 / k for k = 0 , , , ... Given a radius r , the AngularInc function might need to evaluate this expressionseveral times while it searches for the smallest k value for which δ does notexceed the user-specified flatness. To avoid the use of floating-point arithmetic,a truncated Taylor series can be used as an approximation to the square root, aswas done previously in equation (13). With this change, the chord-to-arc erroris expressed as δ = r (cid:18) − ( − (cid:15) − (cid:15) − (cid:15) − ... ) (cid:19) (16)= r (cid:18) (cid:15) + 1128 (cid:15) + 11024 (cid:15) + ... (cid:19) In this form, δ can be calculated using fixed-point arithmetic. All terms inside14he parentheses are negative powers of two, and multiplications of these termsby r can be performed as right-shift operations.The following C++ implementation of the AngularInc function uses equa-tion (16) to calculate δ , but omits the 6th-order and higher terms in the Taylorseries: int AngularInc(FIXED xP, FIXED yP, FIXED xQ, FIXED yQ){ FIXED r = AuxRadius(xP, yP, xQ, yQ);FIXED err2 = r >> 3; // 2nd-order termFIXED err4 = r >> 7; // 4th-order termfor (int k = 0; k < KMAX; ++k){ if (_flatness >= err2 + err4)return k;err2 >>= 2;err4 >>= 4;}return KMAX;} On entry, this function obtains the radius r of the auxiliary circle by calling afunction named AuxRadius , which will be discussed shortly. The global flatnessparameter, _flatness , is a 16.16 fixed-point value.
KMAX is an integer constantthat limits how large k can grow, and, thus, limits how small the angular incre-ment can be.How big should KMAX be? Well, that largely depends on the resolution of thegraphics display. Let r max be the radius of the auxiliary circle for the largestellipse the user might reasonably try to draw on the display, and let δ min be thesmallest flatness setting that the graphics library supports. By truncating afterthe quadratic term in the Taylor series in equation (16), we have the followingapproximation: δ min ≥ r max ( − (cid:113) − (cid:15) ) ≈ (cid:15) r max where (cid:15) = 1 / k max , and k max is the value of the KMAX constant in the previouslisting. Solving for k max , we have k max (cid:38)
12 log ( r max δ min ) For example, with a maximum radius of 5,000 pixels and a minimum flatnessof a quarter of a pixel, we have k max (cid:38) log ( × . ) = 5 . KMAX = 6.In the preceding
AngularInc listing, the arguments passed to the
AuxRadius function are the center-relative coordinates of two conjugate diameter end points, P and Q , that define an ellipse; AuxRadius returns the radius of the ellipse’s aux-iliary circle. Precise, closed-form solutions to the problem of finding this radius15igure 6: A simple strategy for estimating the radius of an ellipse’s auxiliarycircle.exist (Said [15], McCartin [16], Weisstein [17]), but they require floating-pointarithmetic and might run slowly on less-powerful processors. A solution thatuses only fixed-point addition and shift operations might therefore be prefer-able, especially if it quickly generates a good approximation of the radius.To implement an
AuxRadius function that has the desired characteristics,we will use the simple strategy shown in Figure 6 for estimating the radius ofan ellipse’s auxiliary circle. The basic idea is to calculate vector lengths | OP | , | OQ | , | OA (cid:48) | , and | OB (cid:48) | . We then choose the greatest of these lengths as ourestimate for the radius of the ellipse’s auxiliary circle.In Figure 6, the circle and two ellipses are centered at origin O . For each ofthese three figures, we are given end points P = ( x P , y P ) and Q = ( x Q , y Q ) ofa pair of conjugate diameters. We use vector addition to find the corner points A = ( x P + x Q , y P + y Q ) and B = ( x P − x Q , y P − y Q ) of the enclosing squareor parallelogram. For the circle on the left, each of the following four lengthsequals radius r : | OP | , | OQ | , | OA (cid:48) | = | OA | / √
2, and | OB (cid:48) | = | OB | / √
2. For anellipse that’s nearly circular, any one of these four lengths is still a reasonablygood approximation to the radius of the auxiliary circle.For the two eccentric ellipses on the right side of Figure 6, however, someof the corresponding lengths serve as better approximations to the auxiliarycircle radius than others. On the top right, point P comes closest to touchingthe auxiliary circle, and thus r ≈ | OP | gives the best approximation. On thebottom right, point A (cid:48) comes closest to touching the auxiliary circle, and so r ≈ | OA (cid:48) | gives the best approximation.To determine the points A (cid:48) and B (cid:48) for the ellipses on the right side ofFigure 6, we exploit a useful property of affine transformations: they preserveratios between directed segments that lie on the same line. If we view these twoellipses as affine-transformed images of the circle on the left, we see that the16atios A (cid:48) : A and B (cid:48) : B must be preserved. Thus, we have | OA (cid:48) | = | OA | / √ | OB (cid:48) | = | OB | / √ AuxRadius function that uses the strategy just described requires a meansof estimating the lengths of the vectors shown in Figure 6. The following C++function takes as input parameters the x and y components of a vector, andreturns the approximate length of the vector: FIXED VLen(FIXED x, FIXED y){ x = abs(x);y = abs(y);if (x > y)return x + max(y/8, y/2 - x/8);return y + max(x/8, x/2 - y/8);}
This function uses fixed-point addition and shift operations to estimate thevector length. The abs function call in this listing returns the absolute valueof the input argument. The max function call returns the greater of the twoargument values. Assume that abs and max are both inline functions.The VLen function is largely based on a method set forth in a very old Tamilpoem [18] for estimating the length of a hypotenuse. Given a right triangle withsides of length x and y , where 0 ≤ y ≤ x , this method estimates the length of thehypotenuse as x + y/ − x/
8. For y < x/
3, however, a much better approximationis x + y/
8. With this enhancement, the error in the
VLen function’s return valuefalls within the range − . .
78 percent.The following C++ implementation of the
AuxRadius function follows thestrategy presented in the discussion of Figure 6 to estimate the radius of anellipse’s auxiliary circle:
FIXED AuxRadius(FIXED xP, FIXED yP, FIXED xQ, FIXED yQ){ FIXED dP = VLen(xP, yP);FIXED dQ = VLen(xQ, yQ);FIXED dA = VLen(xP + xQ, yP + yQ);FIXED dB = VLen(xP - xQ, yP - yQ);FIXED r1 = max(dP, dQ);FIXED r2 = max(dA, dB);return max(r1 + r1/16, r2 - r2/4);}
The input parameters to this function are the center-relative coordinates of twoconjugate diameter end points that define the ellipse. The
VLen function is calledto obtain the approximate lengths of the four vectors OP , OQ , OA , and OB shown in Figure 6. The longest of the vectors OP , OQ , OA (cid:48) , and OB (cid:48) typically The assumption here is that the compiler converts divisions by constants that are powersof two into right-shift operations.
VLen tends to slightly under-estimate the vector lengths it calculates. The
AuxRadius function compensatesfor these shortfalls by slightly inflating the radius lengths that it calculates. Thelarger of the two lengths | OP | and | OQ | is inflated by 100 × (1 /
16) percent, andthe larger of the two lengths | OA (cid:48) | and | OB (cid:48) | is inflated by 100 × (3 / − / √ AuxRadius function’s return value falls within the range − . .
1. Van Aken, J., Simar, R. (1988). “A Conic Spline Algorithm,”
TMS34010Application Guide , Texas Instruments Incorporated, 255-278.2. Van Aken, J., Simar, R. (1992). “A Parametric Elliptical Arc Algorithm,”
Graphics Gems III , Academic Press, Inc., 164-172.3. Beeler, M., Gosper, R.W., Schroppel, R. (1972). “Item 149 (Minsky),”
HAKMEM , Massachusetts Institute of Technology Artificial IntelligenceLaboratory AIM-239. https://dspace.mit.edu/handle/1721.1/6086
4. Paeth, A.W. (1990). “A Fast Algorithm for General Raster Rotation,”
Graphics Gems , Academic Press, Inc., 179-195.5. Ziegler Hunts, C., Ziegler Hunts, J., Gosper, R.W., Holloway, J. (2011).
Minskys & Trinskys: Exploring an Early Computer Algorithm , 3rd ed.,ZH2G&H (self-published through Blurb, Inc.).
6. Van Aken, J. (2020). “ShapeGen: A lightweight, open-source 2-D graphicslibrary written in C++,” ResearchGate.
7. Van Aken, J. (2018). “A rotated ellipse from three points,” ResearchGate. PostScript Language Reference Manual , 3rd ed. (1999). Adobe SystemsIncorporated, 181.9. Foley, J.D., van Dam, A., Feiner, S.K., Hughes, J.F. (1990).
ComputerGraphics: Principles and Practice , 2nd ed., Addison-Wesley PublishingCompany, 952-953.10. Newman, W.M., Sproull, R.F. (1979).
Principles of Interactive ComputerGraphics , McGraw-Hill, Inc., 27-28.11. Blinn, J. (1996). “How Many Ways Can You Draw a Circle?,”
A TripDown the Graphics Pipeline , Morgan Kaufmann Publishers, Inc., 4-5.12. Rogers, D.F., Hart, J.A (1990).
Mathematical Elements for ComputerGraphics , 2nd ed., McGraw-Hill, Inc., 291-307.183. Watt, A., Watt, M. (1992).
Advanced Animation and Rendering Tech-niques , Addison-Wesley, 69.14. Weisstein, E.W. “Auxiliary Circle.” From MathWorld—A Wolfram WebResource. https://mathworld.wolfram.com/AuxiliaryCircle.html
15. Said, M.A. (2003). “Calibration of an Ellipse’s Algebraic Equation andDirect Determination of its Parameters,”
Acta Mathematica AcademiaePaedagogicae Ny´ıregyh´aziensis , , 221-225.16. McCartin, B.J. (2013). “A Matrix Analytic Approach to Conjugate Di-ameters of an Ellipse,” Applied Mathematical Sciences , (36), 1797-1810.17. Weisstein, E.W. “Ellipse.” From MathWorld—A Wolfram Web Resource. https://mathworld.wolfram.com/Ellipse.html
18. Hattangadi, A.A. (2002).