- Email: [email protected]
The generation of polygons representing circles, ellipses and hyperbolas
C’OMI’I.TlCR
GRAPHICS
.\ND
IM.SGI:
PROCESSING
10,
84-93 11!17!1)
NOTE The
Generation
of Polygons Representing Ellipses and Hyperbolas
P. C. Department
of Computer
MAXWELL
Science,
AND
University
P. W.
Circles,
BAKER
of 1Vew South
Wales,
Sydney,
Australia
Received July 1, 1977; revised December 22, 1977 Algorithms are presented for the computation of vertices of polygons representing circles, ellipses, and hyperbolas. Integers only are used, together with simple operations such as addition and shifting. The curves are represented by parametric differential equations which are solved using techniques baaed on digital differential analyzer (DDA) principles, but using only the registers of a general purpose computer. By employing modifications to increase the integration stepsize and hence the speed of curve generation, polygon points are generated which may be joined by the display vector generator to produce more pleasing results. The deviation of the generated curve from the true curve may be estimated from a detailed examination of all potential sources of error. Finally, by employing the DDA residue retention principle, these algorithms allow effective, fast curve generation with computer word lengths as small as 12 bits. 1.
INTRODUCTION
The digital generation of curves on a graphic display may be accomplished by algorithms based on two different approaches. The first is to consider incremental steps along fixed coordinate axes, where movement is restricted to one of eight neighboring grid points, in a similar fashion to a digital plotter. The second is to utilize the capability of the vector generator in many displays to connect arbitrary points with straight lines, resulting in a polygon approximation to the desired curve. In either approach, the algorithms should operate on integers since many graphics processors do not have floating point hardware and the accuracy of floating representation is not needed. Furthermore, it is desirable that the operations involved in the algorithms be simple ones of addition, subtract,ion, and shifting, since repeated multiplication and division can be excessively time consuming. Several algorithms have been reported which are based on the first approach, using a nonparametric equation of the curve [l-5]. These suffer various disadvantages such as loss of symmetry [Z], need for multiplications [4], special modifications for changing octants or quadrants [l, 4, 51, and problems in specifying and terminating arcs. Denert [6] describes an algorithm based on the 84 0146-664X/79/050084-10$02.00/0 Copyright All rights
0 1979 by Academic Press. Inc. of reproduction in any form reserved.
(:I’,NF:I:ATION
OF CIIICLIl:S,
I~LLIPSES,
HYPI~:I~ROL.4S
8.7
second approach and involving only simple opclrations on int,clgcxrs,although only circles may bc gcneratcd and a large number of operations are rcquircd. Morcx recently, Horn [a] describes a LLconnect the dots” schrtmc which produces elliptical approximations to circles, with satisfactory approximations being determined by appropriate selection of a parameter (this parameter also determines the number of polygon sides). The algorithms presented here arc similar to that drscribed by Horn [a] but possess the following advantages : (i) true circles (which can bc guarant’eed to be within one grid unit accuracy) are produced, rather than elliptical approximations, over a much wider range of the number of polygon sides; (ii) either contiguous set’s of points or polygon points may bc produced; (iii) ellipses and hyperbolas may also be generated; (iv) errors due to truncation using integer arithmet)ic are suppressed; (v) significantly fewer iterations are needed, rrsulting in increascad speed. The basic principle is the same, namely, to consider difference equation approximations to the differential equations describing a circle. However, the view taken htlrc is that t,he difference equations represent a numerical inktgration method for the solution of the differential equations. Since each step of the difference rquations involvcxs computing an increment, WV can USC’tc>chniqucs based on digital differential analyzer (DDA) principles. A DDA, as a special-purpose computer to carry out integration, is normally constructed in hardware, but its operations can be efftlrtivdy carric>d out using the available registers of a gcncral purpose computer. 2. GENERATION
The parametric
OF CIRCLES
equations represent,ing a circle are IE:= -y,
s(0) = ?“,
g = 2,
y(0) = 0.
(1)
A simple method of solving these is to use Euler integration, replacing Eqs. (1) with the difference equations zj+l
= zi + Ax;
yi+1
=
yi
+
which involves
AYi
= xi +
i:;At
= yi +
?j,At
= xi -
yiAt,
= yi +
ziAt.
(2)
Equations (2) form, in fact, the basis for Horn’s algorithm except that his “programming bug” results in using xi+1 rather than pi in the expression for ~i+~. In practice Euler integration cannot be used: due to the high discretization error, and a second-order formula is necessary. However, the simplicity of Eqs. (2) makes it convenient to discuss their solution first, and to indicate later what changes are necessary to produce a practical sy&cm. For a screen resolution of n bits it is convenient to make the step size At = 2-k, where k < 12,since multiplication by At merely constitutes a right shift of k places.
If we allow h: hits t#o spwify the radius r, thca cvc~rj~point on t8hc:circle byill b(h plott~cd if At = 2-“. The total number of itclrations wquired is thw hT = %r/‘At
= 2n X 2,
and I? is given by k = Llog, ?,J + 1. If Eqs. (2) are to be implemented using fixed point, finite-precision arithmetic, special provision must be made to avoid excessive errors due to truncation. This may be done by residue retention, whereby we retain and accumulate that portion of xi (or uI) which would have been lost as a consequence of the multiplicat,ion by At. Since At << 1, the value ciAt contains bits which are less significant, than bhe least significant bit of pi. These bits therefore cannot be added direct’ly to yi and must be stored and accumulated as residues until their contribution is large enough t.o modify the solution. Thrl residues are denoted by R, (for accumulating bits of xiAt) and R, (accumulating bits of $iAt). Mathematically, the retention of residues suppresses the accumulation of wcond-order truncation wror [7]. The reader is referred to [7] for an explanation of DDA operation. A diagrammatic view of the product’ion of Agi (a t’runcatcd version of x,At) is shown in Fig. 1: while Fig. 2a shows the operation in flowchart form, both figures applying to the ternary transfer case whew At = 2- li and the values of Al/i are 0, 1, or -1. The arithmetic in Eqs. (2) and Fig. 2 would be performed in the computer rcgistcrs. If there are at least six rc>gistcrs available (one each for R,, R,, xi, yi, AL~, A!!,) the entire operation can be performed within these registers, otherwise intermcdiat#c results would be wturncd to memory. Since all quantities (&her in the registers or memory) are right-justified, the upper pa& of each word (dotted sections in Fig. 1) merely carry sign information. The arrangement for Axi is exactly analogous. In practice thrrc is no need to test the value of R, as suggested in Fig. 2a, since the correct output value may be obtained by an arithmetic right shift. of R by k places, t’he sign being automatically taken care of. It is imperative to note that in order to preserve residue ret,ention, the value of R must be copied into another register before shifting, otherwise t,hc residue would be lost. The equivalent flow chart, is shown in Fig. 2b, \\-here it is assumc>d that 2’s complement arithmet,ic is
most significant bit of Rx (hence bit of byi)
I--:----
-------
-sign _-__--
7
bits
tnagnitudc lnagnitude I' I I I /
__--_
these bits in R catmut be added toXyi h
-------
I bit
I Rx
0
4
1
----------rI sign bits L--.----------------
FIG.
1. Schematic
--'i
of the production
I
of Ayi.
(ll~?jP:ItATION
OF
CIRCLl
lI:LLIPSES,
IIYI’ISl~I~OT,AS
1
Rx := Rx + xi
si
I arithmetic shift right k places obtain Ayi
. twncate R, to form Ayi = -1
truncate Rx to form Ayi = 0
trunca1e Rx to fan Ayi = +l
I
I
I
R' toX
retain cnly k least significarlt bits of Rx I
I
,
&I i :=
,
it1
+,
a Fro.
2. (a) Flow
chart
for
the
production
of Ayi.
(b) Practical
implementation
of Fig.
2a
used. To minimize the per step error, the correct initial value for both R, and R, is 2k-1, as shown in [7]. From Eqs. (1): z and y are initialized to r and O! respectively. The only problem with using Eys. (2) is t,hat the discrctization error is of a significance such that a spiral, rather than a circle, is generated, as has been shown in [5, 81. Rather than adopt the scheme [S] which results in an elliptical approximation, a more exact alternative is to use a second-order integration method. Adams (trapezoidal) integration may be implemented by replacing Eqs. (2) with xi+1 = xi + &At + +A&-lAt = xi -
(2~~+ +A2/;&~t,
ylibl = yi -I- giAt + sA?ji-lAt = yi + (xi + +AT~-~)A~.
(4
Equations (4) are very easily implemented since they merely involve the addition of $ Ay to R, (and $Ax to R,) before the normal Euler iteration. No truncation occurs if all quantities are multiplied by 2, which leaves the least significant bit of R, and R, untouched except for the addition of $Ax or $Ay. For
values of 1~corresponding to all practical screen resolutions (say in the range 7 to 12), the resultant, circle closes exactly and is accurate to within half a screen increment, this being the maximum DDA deviation. The generation of arcs is readily accomplished by specifying the angle of arc in radians, together with the starting point, the x and y coordinates of which serve as initial values. The number of iterations required to complete the arc is directly related to the angle of arc, since number of iterations
= angle of arc X 2k.
The above assumes the circle center is at (0,O). For other centers a relative initial displacement is all that is required since the outputs AX and Ay are incremental. 3. GENERATION
OF
POLYGON
POINTS
The above algorithms may be modified in order to generate polygon approximations to a circle. Rather than generate every point, the problem is to generate, for example, only 100 points, these being connected by straight lines using the vector generator of the display. For a system with lo-bit resolution, from Eq. (3) we need 27r X 2”’ = 6433 iterations for a full screen circle. A simple method would be to accumulate the Ax and Ay and only output a point every 64 iterations, giving a polygon of 6433/64 ‘v 100 sides. The accumulated value of At between output points is then 64 X 2- IO = 2-4. Clearly, the outputting of points every 64 iterations would be highly inefficient, and the process may be speeded up 64 times if a At of 2-4, rather than 2-lo, is used. A step size as large as this is possible only if the combined errors are such that the generated curve does not excessively diverge from the true circle. This is examined in detail in the Appendix, where it is shown that a step size of 2-5 or less guarantees that the maximum divergence is no more than 1.02 screen increments. However, visually more pleasing curves are generated for At = 2-4, even though the maximum divergence is 2.03 screen increments. In practice, a step size of this order may be achieved by allowing the outputs Axi and Ayi to be several bits in length. The only changes are that these outputs are obtained by a shift of 4 places (rather than 10) and that only the 4 least significant bits of R, and R, are kept. The multibit outputs are fed directly to the display vector generator in order to generate one side of the polygon. The complete circle is generated by sending enough point’s such that the last point is no more than two, but no less than one polygon side length away from the starting point. Since one revolution is 27r radians, the number of points t,o be generated, n, for a given At = 2-k, may be obtained from n = L27r
x
akJ
-
1.
Although the last polygon side is generally longer than the others, experience shows that the results are more visually satisfactory than the alternative of having a shorter final segment. Note that this technique dramatically increases the generation speed since we now need only 99 iterations, rather than 6433, an improvement of 64 times.
(;II:NI~:RATION
OF
il. EXTENSION
CIRCLES, TO
ELLIPSES,
ELLIPSES
AND
HYP~RI3OL.4S
S!)
HYPERBOLAS
-4 more general solution to the harmonic equation is 2 = a cos t, y = b cos (t + S), which gives the ellipse x2 y2 2 cos 6 2 + b2 - abxy = sin2 6.
((9
Such a solution results in Eqs. (4) being no longer valid since Ai and Ag can no longer be easily obtained from Ay and Ax. This means that the previous system, where both Ax and Ay are produced, must be duplicated, utilizing the following relations for y and 2: 21 = --x2, $1 = -yz, i2
=
Xl,
?jz
=
(7)
Yl,
where we define x = x1 and y = yl. Thus there are two exactly identical processes (apart, from initial values) going on, one producing Ax and the other Ay, and each process is the same as that described in the previous section. As before, eit’her discrete dots or polygon points may be generated. The initial value are Xl = a, yl = b cos 6, x2 = 0,
y2 = b sin 6,
(8)
which requires determination of a, b, and 6. Simple coordinate this since for an ellipse given in the form
geometry allows
Ax2 + By2 + 2Hxy = 1, then $=--
B
b2 =
AB -Hz’
A
HZ co52 6 = __ * AR
AB -Hz’
(9)
Ellipses are more likely to be specified by the lengths of the semiaxes, say k1 and k2, and the angle (Ywhich the major axis makes with the x axis. In this case we can utilize the relationships tan ICY = 2H/(A l/k12 + l/W
- B),
= A + B,
l/k~~kz~ = AB - Hz.
The necessary computation of starting values needs to be carried out once for each curve and may therefore be done by the host machine, rather the display processor. Note t’hat a circle will be generated when ICI = k2. A hyperbola may also be generated using the same algorithm by omitting sign change in Eqs. (7). This constitutes a trivial change in the algorithm. the curve is asymptotic, a suitable termination criterion must be employed, as the curve going off the screen.
WC) only than the Since such
00
FIG.
3. Cirrle
and
ellipses
with
100 sides.
The performance for circles and cllipscls on a 1024 X 1024 screen is shown in Fig. 3, where each curve is a polygon of 100 sides. The accuracy is clearly demonstrated since any pronounced deviations would show as some ellipses either intersecting or not touching the circk, which was drawn with a radius equal to the major semiaxis length of the ellipses. Figure 4 shows three examples of hyperbolas which were generated using the same value of At as used to produce Fig. 3. Note that although At is constant during gc>neration, the velocit,y of drawing is such that regions of high(>r curvaturcl arc’ displayed using shorter vectors, resulting in no loss of d&ail. Other advantages are : (1) The circles arc exact and not elliptical approximations. (2) The use of residue retention permits a wide range of At without dcgradation, as well as automatically rounding the circle point coordinates. (3) No truncation means that At can be sufficiently small to generate every point, even with word lengths as small as 12 bits. (4) The same algorithm may bc uschd to gencratc points or polygons. (5) There are no tests to determine o&ants or quadrants. (6) Only fixed point addition, subtraction, and shifting are employed. (7) It is possible to guarantee> that t,hc generated curve is within 1.02 screen increments of the true curve, using At = 2”.
Frc.
4. Examples
of hypeybolas.
(;~:~l~:RATION
OF CIRCLIM,
I:,LI,IPSKS,
HYPE:ItROI,AS
!I1
Considering the last point, n-e may compare thr algorithm here with that of Horn [S]. His elliptical approximation has a half-length of r/(1 - $1c)” \\-here At = l/k. Taking quantization and rounding errors int’o account’, the Appendix shows that with our scheme, a value of At = 2-” guarantees the maximum divergence to be ~1 scrwn incromc~nt for r = 512. Substituting li = 2j in the halflongth clxprrssion results in an error of owr four screen incrcmcbnts. To ensure comparable performance requires that h: hc no smaller than 28. This means Horn’s algorithm requires 8 times as many iterations, and thrwfore polygon sides, as the algorit)hm how, and would hrw.~~ run more slowly. Howwr, such a large value of I; falls in the rangt’, stated by Horn, which products wv(w curve degradation due to truncation in the arit8hmctic. This undesirable effect is easily eliminated, in computers with word lcngt,hs grt3at.w than 20 hi&, by left-justifying all quantit’ics. hltcrrnativcly, with shorter word lengths, rwiduc retcnt8ion may be cxmployed as is done in thr algorithms prcwntctd here. Finally, the algorithm prcsentcd hcrcin will give optimum performance when coded in assembly language, although a high-level language: could bc used. APPENI)IX:
ACCURACY
CONSII)EIIATIONS
Thcrc arc four factors which contribute from the true curve. These arc: 1. st’raight 2. 3. 4.
FOR CIRCLE3
towards the gencratcd curve deviat’ing
Gcomctriral--the polygon sidw diwrgt: from the true arr since they are lirws ; Discretization error in the intclgration formula; Quantization error due to tho screen resolution; Truncation error which accumulates in the DDA.
It is desirable t,hat the polygon sides do not diverge by more than say one screen increment from the true curve. This divergence is now discussed. Results are given for lo-bit screen resolution; similar arguments apply for other resolutions. h worst-case geometrical condition is shown in Fig. 5, where t’wo consecutive points A and R arc a distance s in from the true circle, and t’hc maximum diver-
FIG. 5. Worst-case geometrical condition.
92
MAXWELL
AND
BAKER
gence of the polygon side from this circle is cl. From Fig. 5, 68 [(r - s)2 - (r - d)‘-j sin - = _______--2 r-s - C27-(d r-s
s)l”
>
for
s, d << T,
For a radius with n significant bits, the worst case is r = 2” - 1 + 2”. If we take an approximat.ion of Eq. (11) for small 60 and using d = 1 as a basic criterion, a value of 5 as large as 3 gives At = &J)& 2-(+1)/z
(12)
If we limit At to be 2-4 or smaller, Eq. (12) implies that up to n = 9 (corresponding to a full screen circle) may be used. In practice the situation is better than this for large radii since discretization error causes the generated radius to increase, thereby decreasing the value of s. Exhaustive computer simulation has shown that values of s larger than 0.75 occurred at values of At and radii such that the inward divergence of the polygon side was again less than one screen increment. Hence we may safely ignore any geometrical errors, even for the last polygon segment. From considerations of discretization error alone [9, lo] the maximum step size is 2-5 in order to guarantee that the accumulated error over one cycle is less than one screen increment. For a radius of 512, the maximum discretization error <0.3 increments (in the radius), while for At = 2-4, the error < 1.2 increments. This, however, assumes infinite register lengths and screen resolution. The algorithm as presented effectively uses only lo-bit registers, so we must take quantization and DDA truncation error into account. This latter effect is essentially random, although simulation studies have shown that it produces very similar results to those obtained using high precision followed by quantization and rounding. The final error in the generated points is a combination of discretization, DDA truncation accumulation, and finally quantization and rounding (which is automat,ically effected by proper residue initialization). Exhaustive computer simulation reveals that the maximum point deviations for At = 2-5 and At = 2-4 are 1.02 and 2.03 increments, respectively (r 5 512). Closure is guaranteed, as explained in the body of the paper, by stopping point generation one iteration early and drawing a vector to the start point. So for At = 2-5 we can guarantee that, the circle is always within 1.02 screen units from the true curve. In fact, for almost all radii, the divergence for At = 2-5 is less than 1.00 screen increments. Although the above gives a reliable basis for predicting and limiting curve divergence, in many cases visual appearance is more important than mathematical accuracy. Also, many graphics units cannot display 1024 X 1024 points, a
(:ENEltATION
OF CIRCLES,
ELLIPSES,
HYPEHBOLAS
!)X
common figure being 1024 X 768 (as, for example, in the DEC GT40 and Tektronix 4015). Hence a full circle has a radius of 384, rat,her than 512, thereby reducing the discretization error, which is proportional to the true radius. In all cases tried by the authors, a step size of 2-4 gave visually more pleasing results than 2-j, since with only half as many polygon sides in the former, undesirable kinks tended to be smoothed out. The results in Fig. 3 support this. 6. ACKNOWLEDGMENTS
The authors wish to thank Dr. I’. G. McCrea for his many helpful discussions. This work was supported by the Australian Research Grants Committee. REFERENCES 1. M. L. V. Pitteway, Algorithm for drawing ellipses or hyperbolae with a digital plotter, Cornput. J. 10, 1968, 282-289. 2. P. E. Danielsson, Incremental curve generation, ZEEE Trans. Computers C-19,1970,783-793. 3. B. W. Jordan, W. J. Lennon, and B. D. Holm, An improved algorithm for the generation of nonparametric curves, IEEE Trans. Computers C-22, 1973, 1052-1060. 4. R. A. Metzger, Computer generated graphic segments in a raster display, in Proceedings, AFZPS 1969 Spring Joint Computer Conjerence, Vol. 34, pp. 161-172, AFIPS Press, Montvale, N.J., 1969. 5. B. K. P. Horn, Circle generators for display devices, Computer Graphics and Image Processing 5, 1976, 280-288. 6. E. Denert, A method for computing points of a circle using only integers, Computer Graphics and Zmage Processing 2, 1973, 83-91. 7. P. W. Baker, Suppression of 3rd-order truncation error in linear digital differential analysers, Electron. Lett. 9, 1973, 582-583. X. P. G. McCrea and P. W. Baker, On DDA circle generation for computer graphics, IEEE Trans. Computers C-24, 1975, 1109-1110. 9. P. Henrici, Error Propagation for Difference Methods, Wiley, New York, 1963. 10. E. G. Gilbert, Dynamic-error analysis of digital and combined analog-digital computer systems, Simulation 6, 1966, 241-257.
GRAPHICS
.\ND
IM.SGI:
PROCESSING
10,
84-93 11!17!1)
NOTE The
Generation
of Polygons Representing Ellipses and Hyperbolas
P. C. Department
of Computer
MAXWELL
Science,
AND
University
P. W.
Circles,
BAKER
of 1Vew South
Wales,
Sydney,
Australia
Received July 1, 1977; revised December 22, 1977 Algorithms are presented for the computation of vertices of polygons representing circles, ellipses, and hyperbolas. Integers only are used, together with simple operations such as addition and shifting. The curves are represented by parametric differential equations which are solved using techniques baaed on digital differential analyzer (DDA) principles, but using only the registers of a general purpose computer. By employing modifications to increase the integration stepsize and hence the speed of curve generation, polygon points are generated which may be joined by the display vector generator to produce more pleasing results. The deviation of the generated curve from the true curve may be estimated from a detailed examination of all potential sources of error. Finally, by employing the DDA residue retention principle, these algorithms allow effective, fast curve generation with computer word lengths as small as 12 bits. 1.
INTRODUCTION
The digital generation of curves on a graphic display may be accomplished by algorithms based on two different approaches. The first is to consider incremental steps along fixed coordinate axes, where movement is restricted to one of eight neighboring grid points, in a similar fashion to a digital plotter. The second is to utilize the capability of the vector generator in many displays to connect arbitrary points with straight lines, resulting in a polygon approximation to the desired curve. In either approach, the algorithms should operate on integers since many graphics processors do not have floating point hardware and the accuracy of floating representation is not needed. Furthermore, it is desirable that the operations involved in the algorithms be simple ones of addition, subtract,ion, and shifting, since repeated multiplication and division can be excessively time consuming. Several algorithms have been reported which are based on the first approach, using a nonparametric equation of the curve [l-5]. These suffer various disadvantages such as loss of symmetry [Z], need for multiplications [4], special modifications for changing octants or quadrants [l, 4, 51, and problems in specifying and terminating arcs. Denert [6] describes an algorithm based on the 84 0146-664X/79/050084-10$02.00/0 Copyright All rights
0 1979 by Academic Press. Inc. of reproduction in any form reserved.
(:I’,NF:I:ATION
OF CIIICLIl:S,
I~LLIPSES,
HYPI~:I~ROL.4S
8.7
second approach and involving only simple opclrations on int,clgcxrs,although only circles may bc gcneratcd and a large number of operations are rcquircd. Morcx recently, Horn [a] describes a LLconnect the dots” schrtmc which produces elliptical approximations to circles, with satisfactory approximations being determined by appropriate selection of a parameter (this parameter also determines the number of polygon sides). The algorithms presented here arc similar to that drscribed by Horn [a] but possess the following advantages : (i) true circles (which can bc guarant’eed to be within one grid unit accuracy) are produced, rather than elliptical approximations, over a much wider range of the number of polygon sides; (ii) either contiguous set’s of points or polygon points may bc produced; (iii) ellipses and hyperbolas may also be generated; (iv) errors due to truncation using integer arithmet)ic are suppressed; (v) significantly fewer iterations are needed, rrsulting in increascad speed. The basic principle is the same, namely, to consider difference equation approximations to the differential equations describing a circle. However, the view taken htlrc is that t,he difference equations represent a numerical inktgration method for the solution of the differential equations. Since each step of the difference rquations involvcxs computing an increment, WV can USC’tc>chniqucs based on digital differential analyzer (DDA) principles. A DDA, as a special-purpose computer to carry out integration, is normally constructed in hardware, but its operations can be efftlrtivdy carric>d out using the available registers of a gcncral purpose computer. 2. GENERATION
The parametric
OF CIRCLES
equations represent,ing a circle are IE:= -y,
s(0) = ?“,
g = 2,
y(0) = 0.
(1)
A simple method of solving these is to use Euler integration, replacing Eqs. (1) with the difference equations zj+l
= zi + Ax;
yi+1
=
yi
+
which involves
AYi
= xi +
i:;At
= yi +
?j,At
= xi -
yiAt,
= yi +
ziAt.
(2)
Equations (2) form, in fact, the basis for Horn’s algorithm except that his “programming bug” results in using xi+1 rather than pi in the expression for ~i+~. In practice Euler integration cannot be used: due to the high discretization error, and a second-order formula is necessary. However, the simplicity of Eqs. (2) makes it convenient to discuss their solution first, and to indicate later what changes are necessary to produce a practical sy&cm. For a screen resolution of n bits it is convenient to make the step size At = 2-k, where k < 12,since multiplication by At merely constitutes a right shift of k places.
If we allow h: hits t#o spwify the radius r, thca cvc~rj~point on t8hc:circle byill b(h plott~cd if At = 2-“. The total number of itclrations wquired is thw hT = %r/‘At
= 2n X 2,
and I? is given by k = Llog, ?,J + 1. If Eqs. (2) are to be implemented using fixed point, finite-precision arithmetic, special provision must be made to avoid excessive errors due to truncation. This may be done by residue retention, whereby we retain and accumulate that portion of xi (or uI) which would have been lost as a consequence of the multiplicat,ion by At. Since At << 1, the value ciAt contains bits which are less significant, than bhe least significant bit of pi. These bits therefore cannot be added direct’ly to yi and must be stored and accumulated as residues until their contribution is large enough t.o modify the solution. Thrl residues are denoted by R, (for accumulating bits of xiAt) and R, (accumulating bits of $iAt). Mathematically, the retention of residues suppresses the accumulation of wcond-order truncation wror [7]. The reader is referred to [7] for an explanation of DDA operation. A diagrammatic view of the product’ion of Agi (a t’runcatcd version of x,At) is shown in Fig. 1: while Fig. 2a shows the operation in flowchart form, both figures applying to the ternary transfer case whew At = 2- li and the values of Al/i are 0, 1, or -1. The arithmetic in Eqs. (2) and Fig. 2 would be performed in the computer rcgistcrs. If there are at least six rc>gistcrs available (one each for R,, R,, xi, yi, AL~, A!!,) the entire operation can be performed within these registers, otherwise intermcdiat#c results would be wturncd to memory. Since all quantities (&her in the registers or memory) are right-justified, the upper pa& of each word (dotted sections in Fig. 1) merely carry sign information. The arrangement for Axi is exactly analogous. In practice thrrc is no need to test the value of R, as suggested in Fig. 2a, since the correct output value may be obtained by an arithmetic right shift. of R by k places, t’he sign being automatically taken care of. It is imperative to note that in order to preserve residue ret,ention, the value of R must be copied into another register before shifting, otherwise t,hc residue would be lost. The equivalent flow chart, is shown in Fig. 2b, \\-here it is assumc>d that 2’s complement arithmet,ic is
most significant bit of Rx (hence bit of byi)
I--:----
-------
-sign _-__--
7
bits
tnagnitudc lnagnitude I' I I I /
__--_
these bits in R catmut be added toXyi h
-------
I bit
I Rx
0
4
1
----------rI sign bits L--.----------------
FIG.
1. Schematic
--'i
of the production
I
of Ayi.
(ll~?jP:ItATION
OF
CIRCLl
lI:LLIPSES,
IIYI’ISl~I~OT,AS
1
Rx := Rx + xi
si
I arithmetic shift right k places obtain Ayi
. twncate R, to form Ayi = -1
truncate Rx to form Ayi = 0
trunca1e Rx to fan Ayi = +l
I
I
I
R' toX
retain cnly k least significarlt bits of Rx I
I
,
&I i :=
,
it1
+,
a Fro.
2. (a) Flow
chart
for
the
production
of Ayi.
(b) Practical
implementation
of Fig.
2a
used. To minimize the per step error, the correct initial value for both R, and R, is 2k-1, as shown in [7]. From Eqs. (1): z and y are initialized to r and O! respectively. The only problem with using Eys. (2) is t,hat the discrctization error is of a significance such that a spiral, rather than a circle, is generated, as has been shown in [5, 81. Rather than adopt the scheme [S] which results in an elliptical approximation, a more exact alternative is to use a second-order integration method. Adams (trapezoidal) integration may be implemented by replacing Eqs. (2) with xi+1 = xi + &At + +A&-lAt = xi -
(2~~+ +A2/;&~t,
ylibl = yi -I- giAt + sA?ji-lAt = yi + (xi + +AT~-~)A~.
(4
Equations (4) are very easily implemented since they merely involve the addition of $ Ay to R, (and $Ax to R,) before the normal Euler iteration. No truncation occurs if all quantities are multiplied by 2, which leaves the least significant bit of R, and R, untouched except for the addition of $Ax or $Ay. For
values of 1~corresponding to all practical screen resolutions (say in the range 7 to 12), the resultant, circle closes exactly and is accurate to within half a screen increment, this being the maximum DDA deviation. The generation of arcs is readily accomplished by specifying the angle of arc in radians, together with the starting point, the x and y coordinates of which serve as initial values. The number of iterations required to complete the arc is directly related to the angle of arc, since number of iterations
= angle of arc X 2k.
The above assumes the circle center is at (0,O). For other centers a relative initial displacement is all that is required since the outputs AX and Ay are incremental. 3. GENERATION
OF
POLYGON
POINTS
The above algorithms may be modified in order to generate polygon approximations to a circle. Rather than generate every point, the problem is to generate, for example, only 100 points, these being connected by straight lines using the vector generator of the display. For a system with lo-bit resolution, from Eq. (3) we need 27r X 2”’ = 6433 iterations for a full screen circle. A simple method would be to accumulate the Ax and Ay and only output a point every 64 iterations, giving a polygon of 6433/64 ‘v 100 sides. The accumulated value of At between output points is then 64 X 2- IO = 2-4. Clearly, the outputting of points every 64 iterations would be highly inefficient, and the process may be speeded up 64 times if a At of 2-4, rather than 2-lo, is used. A step size as large as this is possible only if the combined errors are such that the generated curve does not excessively diverge from the true circle. This is examined in detail in the Appendix, where it is shown that a step size of 2-5 or less guarantees that the maximum divergence is no more than 1.02 screen increments. However, visually more pleasing curves are generated for At = 2-4, even though the maximum divergence is 2.03 screen increments. In practice, a step size of this order may be achieved by allowing the outputs Axi and Ayi to be several bits in length. The only changes are that these outputs are obtained by a shift of 4 places (rather than 10) and that only the 4 least significant bits of R, and R, are kept. The multibit outputs are fed directly to the display vector generator in order to generate one side of the polygon. The complete circle is generated by sending enough point’s such that the last point is no more than two, but no less than one polygon side length away from the starting point. Since one revolution is 27r radians, the number of points t,o be generated, n, for a given At = 2-k, may be obtained from n = L27r
x
akJ
-
1.
Although the last polygon side is generally longer than the others, experience shows that the results are more visually satisfactory than the alternative of having a shorter final segment. Note that this technique dramatically increases the generation speed since we now need only 99 iterations, rather than 6433, an improvement of 64 times.
(;II:NI~:RATION
OF
il. EXTENSION
CIRCLES, TO
ELLIPSES,
ELLIPSES
AND
HYP~RI3OL.4S
S!)
HYPERBOLAS
-4 more general solution to the harmonic equation is 2 = a cos t, y = b cos (t + S), which gives the ellipse x2 y2 2 cos 6 2 + b2 - abxy = sin2 6.
((9
Such a solution results in Eqs. (4) being no longer valid since Ai and Ag can no longer be easily obtained from Ay and Ax. This means that the previous system, where both Ax and Ay are produced, must be duplicated, utilizing the following relations for y and 2: 21 = --x2, $1 = -yz, i2
=
Xl,
?jz
=
(7)
Yl,
where we define x = x1 and y = yl. Thus there are two exactly identical processes (apart, from initial values) going on, one producing Ax and the other Ay, and each process is the same as that described in the previous section. As before, eit’her discrete dots or polygon points may be generated. The initial value are Xl = a, yl = b cos 6, x2 = 0,
y2 = b sin 6,
(8)
which requires determination of a, b, and 6. Simple coordinate this since for an ellipse given in the form
geometry allows
Ax2 + By2 + 2Hxy = 1, then $=--
B
b2 =
AB -Hz’
A
HZ co52 6 = __ * AR
AB -Hz’
(9)
Ellipses are more likely to be specified by the lengths of the semiaxes, say k1 and k2, and the angle (Ywhich the major axis makes with the x axis. In this case we can utilize the relationships tan ICY = 2H/(A l/k12 + l/W
- B),
= A + B,
l/k~~kz~ = AB - Hz.
The necessary computation of starting values needs to be carried out once for each curve and may therefore be done by the host machine, rather the display processor. Note t’hat a circle will be generated when ICI = k2. A hyperbola may also be generated using the same algorithm by omitting sign change in Eqs. (7). This constitutes a trivial change in the algorithm. the curve is asymptotic, a suitable termination criterion must be employed, as the curve going off the screen.
WC) only than the Since such
00
FIG.
3. Cirrle
and
ellipses
with
100 sides.
The performance for circles and cllipscls on a 1024 X 1024 screen is shown in Fig. 3, where each curve is a polygon of 100 sides. The accuracy is clearly demonstrated since any pronounced deviations would show as some ellipses either intersecting or not touching the circk, which was drawn with a radius equal to the major semiaxis length of the ellipses. Figure 4 shows three examples of hyperbolas which were generated using the same value of At as used to produce Fig. 3. Note that although At is constant during gc>neration, the velocit,y of drawing is such that regions of high(>r curvaturcl arc’ displayed using shorter vectors, resulting in no loss of d&ail. Other advantages are : (1) The circles arc exact and not elliptical approximations. (2) The use of residue retention permits a wide range of At without dcgradation, as well as automatically rounding the circle point coordinates. (3) No truncation means that At can be sufficiently small to generate every point, even with word lengths as small as 12 bits. (4) The same algorithm may bc uschd to gencratc points or polygons. (5) There are no tests to determine o&ants or quadrants. (6) Only fixed point addition, subtraction, and shifting are employed. (7) It is possible to guarantee> that t,hc generated curve is within 1.02 screen increments of the true curve, using At = 2”.
Frc.
4. Examples
of hypeybolas.
(;~:~l~:RATION
OF CIRCLIM,
I:,LI,IPSKS,
HYPE:ItROI,AS
!I1
Considering the last point, n-e may compare thr algorithm here with that of Horn [S]. His elliptical approximation has a half-length of r/(1 - $1c)” \\-here At = l/k. Taking quantization and rounding errors int’o account’, the Appendix shows that with our scheme, a value of At = 2-” guarantees the maximum divergence to be ~1 scrwn incromc~nt for r = 512. Substituting li = 2j in the halflongth clxprrssion results in an error of owr four screen incrcmcbnts. To ensure comparable performance requires that h: hc no smaller than 28. This means Horn’s algorithm requires 8 times as many iterations, and thrwfore polygon sides, as the algorit)hm how, and would hrw.~~ run more slowly. Howwr, such a large value of I; falls in the rangt’, stated by Horn, which products wv(w curve degradation due to truncation in the arit8hmctic. This undesirable effect is easily eliminated, in computers with word lcngt,hs grt3at.w than 20 hi&, by left-justifying all quantit’ics. hltcrrnativcly, with shorter word lengths, rwiduc retcnt8ion may be cxmployed as is done in thr algorithms prcwntctd here. Finally, the algorithm prcsentcd hcrcin will give optimum performance when coded in assembly language, although a high-level language: could bc used. APPENI)IX:
ACCURACY
CONSII)EIIATIONS
Thcrc arc four factors which contribute from the true curve. These arc: 1. st’raight 2. 3. 4.
FOR CIRCLE3
towards the gencratcd curve deviat’ing
Gcomctriral--the polygon sidw diwrgt: from the true arr since they are lirws ; Discretization error in the intclgration formula; Quantization error due to tho screen resolution; Truncation error which accumulates in the DDA.
It is desirable t,hat the polygon sides do not diverge by more than say one screen increment from the true curve. This divergence is now discussed. Results are given for lo-bit screen resolution; similar arguments apply for other resolutions. h worst-case geometrical condition is shown in Fig. 5, where t’wo consecutive points A and R arc a distance s in from the true circle, and t’hc maximum diver-
FIG. 5. Worst-case geometrical condition.
92
MAXWELL
AND
BAKER
gence of the polygon side from this circle is cl. From Fig. 5, 68 [(r - s)2 - (r - d)‘-j sin - = _______--2 r-s - C27-(d r-s
s)l”
>
for
s, d << T,
For a radius with n significant bits, the worst case is r = 2” - 1 + 2”. If we take an approximat.ion of Eq. (11) for small 60 and using d = 1 as a basic criterion, a value of 5 as large as 3 gives At = &J)& 2-(+1)/z
(12)
If we limit At to be 2-4 or smaller, Eq. (12) implies that up to n = 9 (corresponding to a full screen circle) may be used. In practice the situation is better than this for large radii since discretization error causes the generated radius to increase, thereby decreasing the value of s. Exhaustive computer simulation has shown that values of s larger than 0.75 occurred at values of At and radii such that the inward divergence of the polygon side was again less than one screen increment. Hence we may safely ignore any geometrical errors, even for the last polygon segment. From considerations of discretization error alone [9, lo] the maximum step size is 2-5 in order to guarantee that the accumulated error over one cycle is less than one screen increment. For a radius of 512, the maximum discretization error <0.3 increments (in the radius), while for At = 2-4, the error < 1.2 increments. This, however, assumes infinite register lengths and screen resolution. The algorithm as presented effectively uses only lo-bit registers, so we must take quantization and DDA truncation error into account. This latter effect is essentially random, although simulation studies have shown that it produces very similar results to those obtained using high precision followed by quantization and rounding. The final error in the generated points is a combination of discretization, DDA truncation accumulation, and finally quantization and rounding (which is automat,ically effected by proper residue initialization). Exhaustive computer simulation reveals that the maximum point deviations for At = 2-5 and At = 2-4 are 1.02 and 2.03 increments, respectively (r 5 512). Closure is guaranteed, as explained in the body of the paper, by stopping point generation one iteration early and drawing a vector to the start point. So for At = 2-5 we can guarantee that, the circle is always within 1.02 screen units from the true curve. In fact, for almost all radii, the divergence for At = 2-5 is less than 1.00 screen increments. Although the above gives a reliable basis for predicting and limiting curve divergence, in many cases visual appearance is more important than mathematical accuracy. Also, many graphics units cannot display 1024 X 1024 points, a
(:ENEltATION
OF CIRCLES,
ELLIPSES,
HYPEHBOLAS
!)X
common figure being 1024 X 768 (as, for example, in the DEC GT40 and Tektronix 4015). Hence a full circle has a radius of 384, rat,her than 512, thereby reducing the discretization error, which is proportional to the true radius. In all cases tried by the authors, a step size of 2-4 gave visually more pleasing results than 2-j, since with only half as many polygon sides in the former, undesirable kinks tended to be smoothed out. The results in Fig. 3 support this. 6. ACKNOWLEDGMENTS
The authors wish to thank Dr. I’. G. McCrea for his many helpful discussions. This work was supported by the Australian Research Grants Committee. REFERENCES 1. M. L. V. Pitteway, Algorithm for drawing ellipses or hyperbolae with a digital plotter, Cornput. J. 10, 1968, 282-289. 2. P. E. Danielsson, Incremental curve generation, ZEEE Trans. Computers C-19,1970,783-793. 3. B. W. Jordan, W. J. Lennon, and B. D. Holm, An improved algorithm for the generation of nonparametric curves, IEEE Trans. Computers C-22, 1973, 1052-1060. 4. R. A. Metzger, Computer generated graphic segments in a raster display, in Proceedings, AFZPS 1969 Spring Joint Computer Conjerence, Vol. 34, pp. 161-172, AFIPS Press, Montvale, N.J., 1969. 5. B. K. P. Horn, Circle generators for display devices, Computer Graphics and Image Processing 5, 1976, 280-288. 6. E. Denert, A method for computing points of a circle using only integers, Computer Graphics and Zmage Processing 2, 1973, 83-91. 7. P. W. Baker, Suppression of 3rd-order truncation error in linear digital differential analysers, Electron. Lett. 9, 1973, 582-583. X. P. G. McCrea and P. W. Baker, On DDA circle generation for computer graphics, IEEE Trans. Computers C-24, 1975, 1109-1110. 9. P. Henrici, Error Propagation for Difference Methods, Wiley, New York, 1963. 10. E. G. Gilbert, Dynamic-error analysis of digital and combined analog-digital computer systems, Simulation 6, 1966, 241-257.