# A Comparison of Standard One-Step DDA Circular Interpolators with a New Cheap Two-Step Algorithm

Authors: Leonid Moroz; Jan L. Cieśliński; Marta Stakhiv; Volodymyr Maksymovych
Journal: Modelling and Simulation in Engineering
 (2014)

---

## Abstract

We present and study existing digital differential analyzer (DDA) algorithms for circle generation, including an improved two-step DDA algorithm which can be implemented solely in terms of elementary shifts, addition, and subtraction.


---

## Body

## 1. Introduction

Digital interpolation algorithms are widely used in machine tools with numerical control, graphics displays and plotters, and manipulation robots. Circles and circular arcs frequently appear in computer graphics, computer-controlled printing, and automated control; see [1]. One of most popular methods for generation of circles and arcs is known as digital differential analyzer (DDA). Digital circular interpolators (or angular sweep generators) based on DDA approach are widely used [2–6]. There are many papers where characteristics research results of such interpolators are systematized [7, 8].In this paper a new improved two-step algorithm for DDA circle generation is presented. In paper [9] the idea of applying a Nystrom two-step scheme to circle generation appeared for the first time. Our next paper [10] developed theoretical and geometric aspects of this method. In the present paper we focus on practical and experimental issues. Moreover, here we consider circles of arbitrary radius R, while in [9, 10] we assumed R
=
1.The accuracy of this method is higher than the accuracy of other known algorithms. Because of its simplicity (it uses only elementary shift, addition, and subtraction); this method can also be used in numerical control, planning mechanisms, and so forth.


## 2. DDA Algorithms

The general class of DDA algorithms for circles generation is based on obvious trigonometric transformations describing rotation of vectorR in the coordinate plane x-y (Figure 1);
(1)
x
n
+
1
=
R
sin
(
α
+
Δ
α
)
=
R
sin
α
cos
⁡
Δ
α
+
R
cos
⁡
α
sin
Δ
α
=
x
n
cos
⁡
Δ
α
+
y
n
sin
Δ
α
,
y
n
+
1
=
R
cos
⁡
(
α
+
Δ
α
)
=
R
cos
⁡
α
cos
⁡
Δ
α
+
R
sin
α
sin
Δ
α
=
y
n
cos
⁡
Δ
α
-
x
n
sin
Δ
α
.
Advantages of DDA algorithms include simplicity and high speed in generating circle point coordinate x
n
+
1, y
n
+
1.Figure 1
Rotation of vectorR in the coordinate plane x-y.SubstitutingΔ
α
=
ε
=
2
-
m, where m is an integer (usually m
≥
3), we rewrite (1) as
(2)
x
n
+
1
=
x
n
cos
⁡
ε
+
y
n
sin
ε
,
y
n
+
1
=
y
n
cos
⁡
ε
-
x
n
sin
ε
,
which, obviously, can be expressed by a rotation matrix A:
(3)
A
=
[
cos
⁡
ε
-
sin
ε
sin
ε
cos
⁡
ε
]
.To avoid expensive computation of trigonometric functions DDA algorithms use simpler (cheaper) expressions instead ofcos
⁡
ε and sin
ε. For example, by replacing trigonometric functions by truncated Taylor series one gets so-called simultaneous DDA algorithms [3]. The determinant of the rotation matrix is 1. Matrices of DDA algorithms have different determinants and det
⁡
A
≈
1 only approximately. The closer it equals 1 the more accurate is the corresponding circular interpolator [2]. Approximating cos
⁡
ε
≈
1 and sin
ε
≈
ε, we obtain the simplest (and least accurate) DDA algorithm [2, 11]:
(4)
x
n
+
1
=
x
n
+
ε
y
n
,
y
n
+
1
=
y
n
-
ε
x
n
,
det
⁡
A
=
1
+
ε
2
.
Much more accurate algorithm is obtained using cos
⁡
ε
≈
1
-
(
ε
2
/
2
) and sin
ε
≈
ε [7, 12, 13]:
(5)
x
n
+
1
=
(
1
-
ε
2
2
)
x
n
+
ε
y
n
,
y
n
+
1
=
(
1
-
ε
2
2
)
y
n
-
ε
x
n
,
det
⁡
A
=
1
+
ε
4
4
.
Approximating cos
⁡
ε and sin
ε with first two terms of Taylor series, we obtain the following algorithm [14]:
(6)
x
n
+
1
=
(
1
-
ε
2
2
)
x
n
+
(
ε
-
ε
3
6
)
y
n
,
y
n
+
1
=
(
1
-
ε
2
2
)
y
n
-
(
ε
-
ε
3
6
)
x
n
,
det
⁡
A
=
1
-
ε
4
4
+
ε
6
36
.
Another second order scheme, different form truncated Taylor series, is proposed in [2]:
(7)
x
n
+
1
=
(
1
-
ε
2
2
)
x
n
+
(
ε
-
ε
3
4
)
y
n
,
y
n
+
1
=
(
1
-
ε
2
2
)
y
n
-
(
ε
-
ε
3
4
)
x
n
,
det
⁡
A
=
1
-
ε
4
4
+
ε
6
16
.
All the above DDA algorithms have systematic radial error at each step, which is defined as
(8)
Δ
R
n
=
x
n
2
+
y
n
2
-
R
.
In [7], given the recurrence relations for algorithms that do not have radial errors
(9)
x
n
+
1
=
4
-
ε
2
4
+
ε
2
x
n
+
4
ε
4
+
ε
2
y
n
,
y
n
+
1
=
4
-
ε
2
4
+
ε
2
y
n
-
4
ε
4
+
ε
2
x
n
.
However, implementation of these algorithms requires floating point operations [7]. Replacing the coefficients of (9) with the first two terms of Taylor series [2], we obtain
(10)
4
-
ε
2
4
+
ε
2
≈
1
-
ε
2
2
,
4
ε
4
+
ε
2
≈
ε
-
ε
3
4
.
As a result we obtain, obviously, algorithm given by (7).A lower level of radial error with much simpler hardware implementation can be achieved using the following algorithm [3, 13]:
(11)
x
n
+
1
=
(
1
-
ε
2
2
)
x
n
+
(
ε
-
ε
3
8
)
y
n
,
y
n
+
1
=
(
1
-
ε
2
2
)
y
n
-
(
ε
-
ε
3
8
)
x
n
,
det
⁡
A
=
1
+
ε
6
64
.
Another scheme worthwhile to be mentioned is the so-called “magic circle” interpolator [5, 7, 15]:
(12)
x
n
+
1
=
x
n
+
ε
y
n
,
y
n
+
1
=
(
1
-
ε
2
)
y
n
-
ε
x
n
,
det
⁡
A
=
1
.
This interpolator is more accurate than the simplest scheme (4) and requires the use of only two multipliers. However, in practice it is not used for large ε, because then the generated “circle” becomes an ellipse [2, 15].There is also a group of sequential DDA algorithms [12]. Here are some of them:
(13)
x
n
+
1
=
x
n
+
ε
y
n
,
y
n
+
1
=
y
n
-
ε
x
n
+
1
,
(14)
x
n
+
1
=
x
n
+
ε
y
n
,
y
n
+
1
=
y
n
-
ε
x
n
+
1
,
y
n
+
2
=
y
n
+
1
-
ε
x
n
+
1
,
x
n
+
2
=
x
n
+
1
+
ε
y
n
+
2
,
(15)
x
n
+
1
=
(
1
-
ε
2
2
)
x
n
+
ε
y
n
,
y
n
+
1
=
(
1
-
ε
2
2
)
y
n
-
ε
x
n
+
1
.
Two-step algorithm:
(16)
x
n
+
2
=
x
n
+
1
cos
⁡
ε
-
x
n
,
y
n
+
2
=
y
n
+
1
cos
⁡
ε
-
y
n
with initial conditions: x
0
=
R, y
0
=
0, x
1
=
R
cos
⁡
ε, and y
1
=
R
sin
ε, is known as a “direct” or biquad form of the circular interpolator [14]. Note that its implementation needs only two multipliers.In the process of circle points generation radial errors accumulate, and the law of accumulation is specific to a given algorithm. Note that in most cases, the circle radius is chosen to satisfy condition2
m
-
1
≤
R
≤
2
m [2, 12, 15]. Obviously, the greatest error will occur when R
=
2
m. In Table 1 we give maximum values of the absolute radial systematic errors for R
=
2
m and n
=
2
π
(
1
/
ε
)
=
[
2
π
2
m
], that is, for the full circle generation. Initial conditions are chosen to be x
0
=
0, y
0
=
R.Table 1
Maximum values of the absolute radial error (R
=
2
m).
Number of algorithms
m
=
3
m
=
5
m
=
6
m
=
8
m
=
10
m
=
12
(4)
3.972
3.332
3.234
3.159
3.145
3.142
(5)
0.012
0.771
e
-
3
0.192
e
-
3
0.120
e
-
4
0.749
e
-
6
0.468
e
-
7
(6)
-
0.41
e
-
2
-
0.257
e
-
3
-
0.641
e
-
4
-
0.400
e
-
5
-
0.25
e
-
6
-
0.156
e
-
7
(7)
-
0.012
-
0.771
e
-
3
-
0.192
e
-
3
-
0.120
e
-
4
-
0.74
e
-
6
-
0.468
e
-
7
(11)
0.12
e
-
4
0.47
e
-
7
0.29
e
-
8
0.11
e
-
10
0.45
e
-
13
0.17
e
-
15
(12)
0.262
0.253
0.251
0.250
0.250
0.250
(13)
0.262
0.252
0.251
0.250
0.250
-
0.250
(14)
-
0.063
-
0.015
-
0.781
e
-
2
-
0.195
e
-
2
-
0.48
e
-
3
-
0.122
e
-
3
(15)
-
2.677
-
3.020
-
3.079
-
3.124
-
3.136
-
3.140
(16)
0.18
e
-
7
0.31
e
-
6
0.12
e
-
5
0.19
e
-
4
0.30
e
-
3
0.49
e
-
2Apart from radial error an important role in such devices (especially in computer numerical control (CNC) interpolators) is played by the chord errorE
ch [7], which should not exceed one basic length unit (BLU). Unlike radial error, the chord error does not accumulate. In [7] it is shown that for faster circle point generation for most accurate methods (5)–(7), (9) Δ
α can be chosen from condition Δ
α
=
8
/
R at E
ch
=
1 (we recall that Δ
α
≡
ε); that is
(17)
R
≤
2
2
m
+
3
E
ch
or
R
≤
2
2
m
+
3
.
Since the CNC systems radius can reach value R
max
⁡
=
1
0
8 BLU, the CNC interpolators angles range Δ
α
=
ε
=
2
-
m, where m
=
3,12
¯ overlapping all possible radius values range (the though adopted algorithms can work and at m
>
12).In Table2 we give maximum values of absolute radial systematic error for most accurate algorithms in case R
=
2
2
m
+
3.Table 2
Maximum values of the absolute radial error (R
=
2
2
m
+
3).
Number  of algorithms
m
=
3
m
=
5
m
=
6
m
=
8
m
=
10
m
=
12
(5)
0.7975
0.1973
0.984
e
-
1
0.246
e
-
1
0.614
e
-
2
0.153
e
-
2
(6)
-
0.2642
-
0.657
e
-
1
-
0.328
e
-
1
-
0.818
e
-
2
-
0.205
e
-
2
-
0.511
e
-
3
(7)
-
0.7932
-
0.1972
-
0.984
e
-
1
-
0.246
e
-
1
0.614
e
-
2
0.153
e
-
2
(11)
0.78
e
-
3
0.12
e
-
4
0.15
e
-
5
0.23
e
-
7
0.37
e
-
9
0.57
e
-
11From these results we see that algorithm (11) is the most accurate. Therefore, in next sections we choose algorithm (11) in order to compare it with two new circular interpolators introduced in [10].


## 3. New DDA Algorithms

In the proposed two-step DDA algorithm the valuesx
n
+
2 and y
n
+
2 are determined by the following formula:
(18)
x
n
+
2
=
x
n
-
2
·
sin
ε
·
y
n
+
1
,
y
n
+
2
=
y
n
+
2
·
sin
ε
·
x
n
+
1
.
This is an exact algorithm which can be derived as follows (for another approach, see [10]). We consider the rotation from the point x
n, y
n to x
n
+
2, y
n
+
2; see Figure 2. Hence, we obtain the following equations:
(19)
x
n
+
1
=
x
n
cos
⁡
ε
-
y
n
sin
ε
,
y
n
+
1
=
y
n
cos
⁡
ε
+
x
n
sin
ε
,
(20)
x
n
+
2
=
x
n
+
1
cos
⁡
ε
-
y
n
+
1
sin
ε
,
y
n
+
2
=
y
n
+
1
cos
⁡
ε
+
x
n
+
1
sin
ε
.
From (20) we obtain values x
n
+
1 and -
y
n
+
1:
(21)
x
n
+
1
=
y
n
+
2
-
y
n
+
1
cos
⁡
ε
sin
ε
,
-
y
n
+
1
=
x
n
+
2
-
x
n
+
1
cos
⁡
ε
sin
ε
.
Substituting values given by (21) into (20) we get
(22)
x
n
+
2
-
2
x
n
+
1
cos
⁡
ε
+
x
n
=
0
,
y
n
+
2
-
2
y
n
+
1
cos
⁡
ε
+
y
n
=
0
.
Here we recognize the exact discretization of the classical harmonic oscillator equation [16, 17]. From (22) it follows that
(23)
x
n
+
1
=
x
n
+
2
+
x
n
2
cos
⁡
ε
,
y
n
+
1
=
y
n
+
2
+
y
n
2
cos
⁡
ε
.
Substituting obtained values in (21), we get
(24)
x
n
+
1
=
y
n
+
2
-
(
(
y
n
+
2
+
y
n
)
/
2
cos
⁡
ε
)
cos
⁡
ε
sin
ε
,
-
y
n
+
1
=
x
n
+
2
-
(
(
x
n
+
2
+
x
n
)
/
2
cos
⁡
ε
)
cos
⁡
ε
sin
ε
which can easily be reduced to
(25)
x
n
+
1
=
y
n
+
2
-
y
n
2
sin
ε
,
-
y
n
+
1
=
x
n
+
2
-
x
n
2
sin
ε
,
that is, to (18). The algorithm works with initial conditions: x
0
=
R, y
0
=
0, x
1
=
R
cos
⁡
ε, and y
1
=
R
sin
ε, where R
=
2
m (m is an integer). Approximating sin
ε
≈
ε one can transform (18) into another two-step algorithm:
(26)
x
n
+
2
=
x
n
-
2
·
ε
·
y
n
+
1
,
y
n
+
2
=
y
n
+
2
·
ε
·
x
n
+
1
,
with x
0
=
R, y
0
=
0, x
1
=
R
1
-
ε
2, y
1
=
R. This algorithm was first proposed in [9]; see also [10].Figure 2
Generation of 3 subsequent points along the circle.Functional diagram of interpolator which reproduces (26) is shown in Figure 3. Interpolator contains 6 registers (Register 1, Register 2, Register 3, Register 4, Output—x
n
+
2
,, and Output—y
n
+
2), 2 multipliers (Right Shift Register 1 Right Shift Register 2), subtracter and adder.Figure 3
Functional diagram for the circular interpolator (26).In this scheme, at the beginning (before the interpolator starts) values ofx
0 and y
0 are recorded in Register 2 and Register 4, respectively, while values of x
1 and y
1 are recorded in Register 1 and Register 3. Then, values of x
n
+
1 and x
n are stored in Register 1 and Register 2, and values of y
n
+
1 and y
n. are stored in Register 3 and Register 4. Right Shift Register 1 and Right Shift Register 2 perform multiplication of values y
n
+
1 and x
n
+
1 by 2
ε. Thus, at outputs of subtractor and adder parallel binary codes are formed in accordance with (26). Rounded values are recorded in Output x
n
+
2 and Output y
n
+
2.Table3 contains maximum values of absolute radial systematic error of proposed algorithm, given by (26), after full circle simulation in the case
(27)
d
=
24
,
d
=
32
,
ε
=
2
-
m
,
R
=
2
m
,
m
=
3
,
…
,
10
,
n
=
2
π
1
ε
,
where d is structure elements datapath (in bits) and n is number of iterations. Table 4 contains the same data for case R
=
2
2
m
+
3. Note that results given in Tables 1 and 2 do not include round-off errors; hence they do not depend on d.Table 3
Maximum values of the absolute  radial error (R
=
2
m) for proposed Algorithm (26).
m
Absolute radial error
d
=
24
d
=
32
Max
Min
Max
Min
3
7.712
e
-
8
-
2.567
e
-
7
4.377
e
-
10
-
9.137
e
-
10
4
1.215
e
-
7
-
6.462
e
-
7
1.284
e
-
9
-
2.666
e
-
9
5
2.458
e
-
7
-
1.267
e
-
6
9.422
e
-
10
-
5.362
e
-
9
6
4.643
e
-
7
-
2.645
e
-
6
2.223
e
-
9
-
8.292
e
-
9
7
7.045
e
-
7
-
4.942
e
-
6
3.614
e
-
9
-
1.856
e
-
8
8
1.655
e
-
6
-
9.678
e
-
6
6.936
e
-
9
-
3.794
e
-
8
9
3.258
e
-
6
-
1.898
e
-
5
1.245
e
-
8
-
7.470
e
-
8
10
6.758
e
-
6
-
3.761
e
-
5
2.510
e
-
8
-
1.456
e
-
7Table 4
Maximum values of the absolute radial error (R
=
2
2
m
+
3) for proposed Algorithm (26).
m
Absolute radial error
d
=
24
d
=
32
Max
Min
Max
Min
3
6.965
e
-
8
-
3.866
e
-
7
4.533
e
-
10
-
1.671
e
-
9
4
1.768
e
-
7
-
8.276
e
-
7
4.654
e
-
10
-
3.550
e
-
9
5
3.010
e
-
7
-
1.176
e
-
6
1.260
e
-
9
-
5.494
e
-
9
6
4.245
e
-
7
-
2.292
e
-
6
1.923
e
-
9
-
9.138
e
-
9
7
8.948
e
-
7
-
5.040
e
-
6
3.707
e
-
9
-
1.927
e
-
8
8
1.454
e
-
6
-
9.834
e
-
6
8.918
e
-
9
-
3.704
e
-
8
9
3.424
e
-
6
-
1.936
e
-
5
1.307
e
-
8
-
7.389
e
-
8
10
6.792
e
-
6
-
3.648
e
-
5
2.792
e
-
8
-
1.437
e
-
7Analogical data for algorithm given by (11) are given in Table 5 (we consider, as before, case (27)). Table 6 contains the same data for R
=
2
2
m
+
3.Table 5
Maximum values of the absolute radial error (R
=
2
m) for Algorithm (11).
m
Absolute radial error
d
=
24
d
=
32
Max
Min
Max
Min
3
1.214
e
-
5
0
1.215
e
-
5
0
4
9.641
e
-
7
-
8.277
e
-
7
7.525
e
-
7
0
5
5.496
e
-
7
-
2.386
e
-
6
4.552
e
-
8
0
6
5.508
e
-
7
-
4.660
e
-
6
5.781
e
-
9
-
1.682
e
-
8
7
1.235
e
-
6
-
9.674
e
-
6
6.069
e
-
9
-
3.579
e
-
8
8
4.009
e
-
6
-
1.858
e
-
5
1.524
e
-
8
-
7.158
e
-
8
9
6.478
e
-
6
-
3.651
e
-
5
2.034
e
-
8
-
1.484
e
-
7
10
1.268
e
-
5
-
7.285
e
-
5
5.564
e
-
8
-
2.853
e
-
7Table 6
Maximum values of the absolute radial error (R
=
2
2
m
+
3) for Algorithm (11).
m
Absolute radial error
d
=
24
d
=
32
Max
Min
Max
Min
3
7.780
e
-
4
0
7.781
e
-
4
0
4
9.650
e
-
5
0
9.632
e
-
5
0
5
1.165
e
-
5
0
1.204
e
-
5
0
6
2.469
e
-
6
-
3.933
e
-
6
1.501
e
-
6
0
7
1.731
e
-
6
-
9.497
e
-
6
1.878
e
-
7
0
8
2.359
e
-
6
-
1.942
e
-
5
3.598
e
-
8
-
6.302
e
-
8
9
5.371
e
-
6
-
3.779
e
-
5
2.655
e
-
8
-
1.421
e
-
7
10
1.057
e
-
5
-
7.478
e
-
5
5.028
e
-
8
-
2.914
e
-
7Comparing errors in Tables3 and 5, and in Tables 4 and 6, we see that the absolute radial errors of our algorithm (26) (Tables 3 and 4) are smaller than errors of well-known algorithm (11) (see Tables 5 and 6). For example, we can compare values of errors at R
=
2
m
(
m
=
10
) and d
=
24 given in Tables 3 and 5. Consequently, we will get that an absolute radial error is smaller by (
max
⁡
(
Table
5
)
)
/
(
max
⁡
(
Table
3
)
)
=
(
1.267
e
-
5
)
/
(
6.792
e
-
6
)
=
1.876 times.We also see that the smaller ism, the greater is the difference between the corresponding errors of both algorithms. For example, for (
m
=
3
) and d
=
32 we obtain that an absolute radial error is smaller by (
max
⁡
(
Table
6
)
)
/
(
max
⁡
(
Table
4
)
)
=
(
7.781
e
-
4
)
/
(
4.533
e
-
10
)
≈
1716523 times.Finally, Figures4, 5, 6, and 7 present numbers of correct bits as a function of m, for R
=
2
m and R
=
2
2
m
+
3 and for d
=
24 and d
=
32, respectively.Figure 4
Number of correct bits of result as a function ofm at R
=
2
m and d
=
24.Figure 5
Number of correct bits of result as a function ofm at R
=
2
2
m
+
3 and d
=
24.Figure 6
Number of correct bits of result as a function ofm at R
=
2
m and d
=
32.Figure 7
Number of correct bits of result as a function ofm at R
=
2
2
m
+
3 and d
=
32.Presented figures clearly show that the proposed algorithm (26) yields greater number of correct bits in comparison to algorithm (11).


## 4. Conclusion

A new circle generation DDA algorithm has been proposed. In terms of hardware implementation and precision of circle generation the proposed DDA algorithm for circle interpolator is more efficient as compared with existing algorithms. We hope that the proposed algorithm will find some applications in computer graphics (display screens), computer-controlled printing, automated control, and so forth.



---

*Source: 916539-2014-01-20.xml*