 # Fast exact digital differential analyzer for circle generation

In the first part of the paper we present a short review of applications of digital differential analyzers (DDA) to generation of circles showing that they can be treated as one-step numerical schemes. In the second part we present and discuss a novel fast algorithm based on a two-step numerical scheme (explicit midpoint rule). Although our algorithm is as cheap as the simplest one-step DDA algoritm (and can be represented in terms of shifts and additions), it generates circles with maximal accuracy, i.e., it is exact up to round-off errors.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

In spite of apparent simplicity fast generation of circles and other curves has always been a subject of numerous theoretical and practical studies, finding applications, among others, in digital plotting, graphical display and numerical machine tool control, compare [12, 15, 21, 25]. Leaving aside the fundamental algoritm of Bresenham  and techniques based on spline functions , in this paper we focus on digital differential analyzers (DDA), see for instance [11, 16, 19]. In section 2 we present brief but exhaustive discussion of digital differential analyzers from a unified point of view, treating them as special one-step difference numerical schemes. Then, in section 3, we introduce and discuss a new DDA scheme, cheap and accurate, based on a two-step numerical scheme.

In terms of the natural arc parameter the equation of the circle of radius is , . The corresponding differential equation is

 dxdϑ=−y ,dydϑ=x . (1)

Considering various discretizations of (1) we can obtain any DDA algorithm.

## 2 One-step numerical approximations

In this section we present a unified approach to digital differential analysers representing them by one-step numerical schemes. Therefore they can be represented by a general matrix difference equation of the first order:

 [xn+1yn+1]=[abcd][xnyn],[abcd]≡A , (2)

where depend on the time step . In order to approximate the system (1) these coefficients have to satisfy

 limh→0[abcd]=≡I,limh→0A−Ih=[0−110]≡J. (3)

In other words,  , i.e.,

 a=1+O(h2),  d=1+O(h2),  b=−h+O(h2),  c=h+O(h2).

The qualitative asymptotic behaviour of the discrete solution is determined by the characteristic equation

and its roots (eigenvalues of the matrix

) are given by

where .

### 2.1 Logarithmic spirals

The most popular DDA schemes have the form

 [xn+1yn+1]=[a−cca][xnyn]. (6)

It is convenient to denote

 ρ=√a2+c2 ,θ=arctanca , (7)

i.e.,  ,  . Then,

 [xnyn]=[ρcosθ−ρsinθρsinθρcosθ]n[x0y0]=ρn[cosnθ−sinnθsinnθcosnθ][x0y0].

Denoting   ,  ,  we finally get

 xn=r0ρncos(nθ+φ0) ,yn=r0ρnsin(nθ+φ0) . (8)
###### Corollary 2.1.

All points generated by (6) lie on the logarithmic spiral defined by:

 r=r0ek(φ−φ0) ,k:=θ−1lnρ , (9)

where are polar coordinates on the plane , compare .

From computational point of view the polynomial form of

is preferred. DDA algorithms can be classified by the order of these polynomials

. Below we present several schemes of this kind with the corresponding values of and (the exact algorithm should have and ).

#### First order simultaneous DDA algorithm

 a=1 ,c=h ,ρ2=1+h2 ,k=h2−h312+… . (10)

This is the conventional DDA method .

#### Second order simultaneous DDA algorithm

 a=1−12h2 ,c=h , ρ2=1+14h4 ,k=h38−h548+… , (11)

see [6, 15].

#### Third order simultaneous DDA algorithm

 a=1−12h2,c=h−16h3,ρ2=1−112h4+136h6,k=−h324+h572+…, (12)

compare .

#### Matsushiro’s analyzer

 a=1−12h2,c=h−14h3,ρ2=1−14h4+116h6,k=−h38+h548+…. (13)

An important point is that this scheme (with replaced by ) was implemented in terms of shifts and additions (without using multiplications) and patented, see [17, 18].

#### The best circular interpolator of the third order

Interestingly enough, the Matsushiro algoritm can be easily improved by changing the last term, see [8, 20]:

 a=1−12h2,c=h−18h3,ρ2=1+164h6,k=h5128+…. (14)

We are going to show that this is the best spiral-like circular interpolator of the third order in

.

###### Lemma 2.2.

The most accurate circular interpolator of the form (6), where are polynomials of the third order in , is given by (14)

Proof:  We assume

 a=1+a1h+a2h2+a3h3 ,c=h+c2h2+c3h3 (15)

Then we compute :

 ρ2=1+2a1h+(2a2+a21+1)h2+2(a3+a1a2+c2)h3+…+(a23+c23)h6.

Equating to zero five first coefficients we obtain a system of 5 equations for 5 unknowns, which can be solved easily:

 2a1=0 ⟹a1=0 ,2a2+a21+1=0 ⟹a2=−12 ,a3+a1a2+c2=0 ⟹a3=−c2 ,a22+2a1a3+c22+2c3=0 ⟹c3=−18−12c22 ,a2a3+c2c3=0 ⟹12c2(34−c22)=0 . (16)

Then, the coefficient by is given by  . The last equation of (16) yields two possibilities: either or . Corresponding coefficients by are given by and , respectively. Therefore the best accuracy is attained in the first case, which leads to (14) (the second case has another disadvantage: coefficients and are irrational).

Interestingly enough, considering -expansion of the function we get the same final results (although the starting system of equations is more complicated). The case yields while in the second case, , we get .

### 2.2 Elliptical deformations

Discrete points generated by (2) lie on an ellipse if and only if   and , which is equivalent to

The conditions are known as Barkhausen criteria .

#### First order sequential DDA algorithm

 [xn+1yn+1]=[1−hh1−h2][xnyn]. (18)

This is the generator of the so called “magic circle” [17, 21, 23].

### 2.3 Elliptical spirals

Two complex eigenvalues () are obtained when

 (a−d)2+4bc<0 , (19)

and, for the generated circle deforms into an elliptical spiral. The following example can be found in .

#### Second order sequential DDA algorithm

 [xn+1yn+1]=[1−12h2−hh1−32h2][xnyn] (20)

One can easily verify, that . Therefore (19) is satisfied for .

Sequentials algorithms are derived from simulataneous algorithms by replacing by in the equation for , i.e., simultaneous algorithm (6) of order is transformed into

 xn+1=axn−cyn ,yn+1=cxn+1+ayn=acxn+(a−c2)yn , (21)

and then truncated up to the order , see . In the case we get the scheme (20).

### 2.4 Exactly circular interpolation

Digital differential analyzer of the form (2) yields an exact circle (up to round-off errors) if and only if , and . Then, the matrix can be parameterized by a single parameter

 A=[cosϕ−sinϕsinϕcosϕ], (22)

where depends on time step . In other words,

 [xn+1yn+1]=[C(h)−S(h)S(h)C(h)][xnyn] (23)

where . Taking we obtain the exact numerical scheme (compare [21, 24])

 [xn+1yn+1]=[cosh−sinhsinhcosh][xnyn]. (24)

Exact numerical schemes can be constructed for any matrix differential equation of the first order, see [9, 10].

The Taylor expansion of (24) up to the third order yields the scheme (12), compare . We recall that this is not the best scheme of the third order, see Lemma 2.2.

#### Implicit midpoint rule

In the case of linear equations the implicit midpoint rule can be represented in an explicit form. Indeed,

 xn+1−xnh=−yn+yn+12 ,yn+1−ynh=−xn+xn+12 , (25)

can be rewritten in a matrix form as follows

 [112h−12h1][xn+1yn+1]=[1−12h12h1][xnyn]. (26)

Thus we obtain another algortihm of the form (6), where

 a=4−h24+h2 ,c=4h4+h2 . (27)

see . Note that and , so the implicit midpoint rule generates the circle exactly. However, this algorithm has a relatively high cost (multiplications in every step). Actually the scheme (27) was a motivation for the derivation of the Matsushiro algorithm (13), see . Indeed, (13) is the Taylor truncation of (27) up to the third order.

## 3 New fast circular interpolator

The main result of our paper consists in applying a two step method to the circular interpolation. The explicit midpoint rule (see ), known also as 2-step Nyström method 

, applied to a general ordinary differential equation

(where and is a given function) reads:

 zn+2=zn+2hf(tn+1,zn+1) . (28)

In the case of the circular interpolation   and . Therefore we obtain:

 xn+2=xn−2hyn+1 ,yn+2=yn+2hxn+1 , (29)

compare . This is a two step method and we have to prescribe not only but also  . As usual, we can rewrite this scheme in a one-step form increasing the dimension of the matrix problem, namely:

 ⎡⎢ ⎢ ⎢ ⎢⎣xn+2yn+2xn+1yn+1⎤⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢⎣0−2h102h00110000100⎤⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢⎣xn+1yn+1xnyn⎤⎥ ⎥ ⎥ ⎥⎦. (30)

### 3.1 First integrals

In this section we are going to show that (29) preserves circular trajectories exactly provided that initial conditions are appropriately chosen.

###### Lemma 3.1.

If    is a first integral of the system (29), then    is also a first integral. What is more, in this case both first integrals are linearly dependent, namely

 xnyn+1−xn+1yn=h(x2n+y2n) (31)

Proof:  Using (29) we compute

 x2n+2+y2n+2=x2n+y2n−4h(xnyn+1−xn+1yn)+4h2(x2n+1+y2n+1) .

If  , then, obviously,  .

Motivated by above results, we define:

 ξn=⎡⎢ ⎢⎣x2n+y2nx2n+1+y2n+1xnyn+1−xn+1yn⎤⎥ ⎥⎦,B=⎡⎢⎣01014h2−4h02h−1⎤⎥⎦. (32)
###### Lemma 3.2.

For any satisfying (29) we have

 ξn+1=Bξn, (33)

Proof:  A simple straightforward calculation.

###### Theorem 3.3.

All points generated by (29) lie exactly (up to round-off errors) on the circle provided that initial conditions satisfy

 x20+y20=r2,x21+y21=r2,x0y1−x1y0=hr2. (34)

Proof:  We easily verify that

 Bf1=f1 ,f1:=⎡⎢⎣11h⎤⎥⎦. (35)

Initial conditions (46) can be shortly written as . By (33) we have for any . It means that    for any .

###### Remark 3.4.

The simplest initial conditions satisfying (46) read

 x0=r ,y0=0 ,x1=r√1−h2 ,y1=hr , (36)

compare .

### 3.2 Low computational cost

The scheme (29) has a very low computation cost, to be compared only with the simplest one-step scheme (10). Indeed, both (29) and (10) can be implemented in terms of two additions and two shifts. However, (29) produces an exact circle while the accuracy of (10) is rather limited.

There is only one more complicated computation, at the start, when initial conditions are are determined. However, taking into account the Taylor expansion

 √1−h2=1−12h2−18h4−116h6−5128h8−7256h10+O(h12) (37)

we see that even in that case one can use few shifts to produce . Indeed, let and . Then

 √1−h2=1−2−2m−1−2−4m−3−2−6m−4−2−8m−5−2−8m−7+… (38)

and

 x1=2N−2N−2m−1−2N−4m−3−2N−6m−4−2N−8m−5−2N−8m−7+… (39)

where only integer coefficients are left.

The case (i.e., ) corresponds to the circle approximated by regular dodecagon (because then , , compare (36)). Assuming we may confine ourselves to the first three terms.

### 3.3 Stability

 (1−λ)(λ2+2(1−2h2)λ+1)=0, (40)

hence eigenvalues of are given by:

 λ1=1,  λ2=2h2−1+2h√h2−1,  λ3=2h2−1−2h√h2−1. (41)

All the eigenvalues lie on the unit circle: for

. Corresponding eigenvectors are denoted by

, , namely

 Bf1=f1 ,Bf2=λ2f2 ,Bf3=λ3f3, (42)

( is given by (35) and the form of other eigenvectors is not important). The initial condition can always be represented as a linear combination of eigenvectors: , where we assume (compare Theorem 3.3) and , are small perturbations. Using (42) and (33) we obtain

 ξn=f1+bλn2f2+cλn3f3. (43)

Taking into account  , we see that the initial perturbations remain unchanged during the evolution.

### 3.4 Exact circular interpolator

The scheme (29) produces exact circular trajectory but the period is modified. For instance, in the case of regular dodecagon () the total change of after 12 steps is while the exact value is, of course, . In general, the period of the circular interpolator (29) is given by

 T=2πharcsinh. (44)

In order to preserve exactly not only the trajectory but also the period of the circular motion, we can use a nonstandard modification of (29), following the approach of [9, 10]:

 xn+2=xn−2δyn+1 ,yn+2=yn+2δxn+1 , (45)

where is a given function of .

###### Theorem 3.5.

For any all points generated by (45) lie exactly (up to round-off errors) on the circle provided that initial conditions satisfy

 x20+y20=r2,x21+y21=r2,x0y1−x1y0=δr2. (46)

Proof: It is enough to repeat the proof of Theorem 3.3 (with replaced by in appropriate places).

The period of the circular interpolator (45) is

 T=2πharcsinδ(h). (47)
###### Theorem 3.6.

If  , then scheme

 xn+2=xn−2yn+1sinh ,yn+2=yn+2xn+1sinh ,x0y1−x1y0=r2sinh , (48)

generates the circle exactly (up to round-off errors) preserving the period.

Proof: The scheme (45) will preserve the period exactly if   and . Taking into account well known trigonometric identities

 cos(hn+2h)−cos(hn)=−2sinhsin(hn+h) ,sin(hn+2h)−sin(hn)=2sinhcos(hn+h) , (49)

we conclude that scheme (45) is exact for .

The exact interpolator (48) uses two multiplication by at every step, hence it is relatively expensive (but two times cheaper than its one-step conterpart (24)).

The scheme (45) can serve as a starting point for deriving cheaper algorithms by taking -polynomials instead of . All these schemes preserve exactly the circular trajectory. The best approximation of the period are obtained for Taylor truncations. In particular, taking

 δ=h−16h3 (50)

we get the best scheme (45) of order 3. The multiplication by can be replaced by shifts because

 16=18(1−14)=18∞∑k=04−k=∞∑k=02−3−2k, (51)

and we can approximate by given by

 δN=h−h3N∑k=02−3−2k, (52)

where is chosen to assure the required accuracy.

## 4 Conclusions

We discussed DDA circular interpolators based on one-step numerical schemes. We have shown that the best cheap interpolator of the third order is given by (14), see Lemma 2.2. This algorithm can be expressed in terms of shifts and additions only, and it has quite good accuracy. Exact interpolators based on one-step methods are more expensive because multiplications in every step are required.

In section 3 we presented a family of circular interpolators, based on a two-step method, namely the explicit midpoint rule. All of them interpolate the circle exactly (up to round-off errors). The simplest interpolator, given by (29), is very cheap and is very good for graphical applications (when the period of the circular motion is irrelevant). The scheme (48) preserves exactly also the period but is more expensive. Taking an appropriate polynomial we can control the accuracy and computational cost of the scheme (45), obtaining interpolators suitable for different purposes.

Acknowledgement. The first author (J.L.C.) is partly supported by the National Science Centre (NCN) grant no. 2011/01/B/ST1/05137.

## References

• 
• 
• 
• 
•  C.Baumgarten, G.Farin: Approximation of logarithmic spirals, Computer Aided Geometric Design 14 (6) (1997) 515-532.
•  C.A.Bergren: A simple algorithm for circular interpolation, Control Engineering 18 (9) (1971) 57-59.
•  J.E.Bresenham: A linear algorithm for incremental digital display of circular arcs, Commun. ACM 20 (2) (1977) 100-106.
•  I.Chami: A high precision digital differential analyzer for circular interpolation, Tishreen Univ. J. Stud. Sci. Res., Eng. Sci. Ser. 29 (1) (2007).
•  J.L.Cieśliński, B.Ratkiewicz: On simulations of the classical harmonic oscillator equation by difference equations, Adv. Difference Equ.  2006 (2006) 40171 (17pp).
•  J.L.Cieśliński: On the exact discretization of the classical harmonic oscillator equation, J. Difference Equ. Appl. 17 (2011) 1673-1694.
•  P.E.Danielsson: Incremental curve generation, IEEE Trans. Comput. C-19 (9) (1970) 783-793.
•  X.Y.Fan, Y.H.Guo, S.C.Li: New digital differential analyzer interpolation algorithm, Advanced Science Letters 6 (2012) 692-695.
•  E.Hairer, S.P.Nørsett, G.Wanner: Solving ordinary differential equations I. Nonstiff problems, Springer, Berlin 1987.
•  A.Iserles: A first course in the numerical analysis of differential equations, second edition, Cambridge Univ. Press 2009.
•  Y.Koren: Interpolator for a computer numerical control system, IEEE Trans. Comput. C-25 (1) (1976) 32-37.
•  F.S.Lim, Y.S.Wong, M.Rahman: Circular interpolators for numerical control: A comparison of the modified DDA techniques and an LSI interpolator, Computers in Industry 18 (1992) 41-52.
•  N.Matsushiro: A new digital differential analyzer for circle generation, IEICE Trans. Inf. Syst. E81-D (2) (1998) 239-242.
•  N.Matsushiro, I.Oyake: Method of and device for circle generation, US Patent No. 4,999,797 (March 1991).
•  P.G.McCrea, P.W.Baker: On digital differential analyzer (DDA) circle generation for computer graphics, IEEE Trans. Comput. C-24 (11) (1975) 1109-1110.
•  L.V.Moroz: A method for studying characteristics of digital sine-cosine generators, Information Extraction and Processing 36 (112) (2012) 84-90 [in Ukrainian].
•  W.M.Newman, R.F.Sproull: Principles of interactive computer graphics, McGraw-Hill 1979.
•  L.A.Piegl, W.Tiller: Circle approximation using integral B-splines, Computer-Aided Design 35 (2003) 601-607.
•  C.S.Turner: Recursive discrete-time sinusoidal oscillators, IEEE Signal Processing Magazine, May 2003.
•  W.P.Wang, C.Y.Wang: Difference method for generation of circular arcs and ellipses, Computer-Aided Design 21 (1) (1989) 33-37.
•  L.Yong-Kui: Algorithm for circle approximation and generation, Computer-Aided Design 25 (3) (1993) 169-171.