Circle to be found
Figure 6.6. Original data in (x,y) domain.
Again it is possible to reduce the amount of work by using the gradient direction to indicate
the likely arc within which the circle centre is expected to lie. Figure 6.7 illustrates this
technique.
It is possible to look for the following types of circles:
different radii plot in (a,b,R) space
different radii, same vertical centres plot in (b,R) space
Four cicles coincide here
only
different radii, same horizontal centres plot in (a,R) space
Figure 6.7 Illustration of Hough circle transform (looking for circles radius 1/(cid:214) 2).
Corresponding accumulator circles in (a,b) domain.
If the circle radius is known to be one of three values, say, then (a,b,R) space can be three
planes of (a,b) arrays.
The following points are important:
1. As the number of unknown parameters increases, the amount of processing
increases exponentially.
2. The Hough technique above can be used to discover any edge that can be expressed
as a simple identity.
3. The generalized Hough transform can also be used to discover shapes that can not
be represented by simple mathematical identities. This is described below.
Technique 6.3. The generalized Hough transform.
in its most general form-of any size or orientation in an USE. Find a known shape (cid:190)
image. In practice it is best to go for a known size and orientation.
OPERATION. Some preparation is needed prior to the analysis of the image. Given
the object boundary, and assuming that the object in the image is of the same size and
orientation (otherwise a number of accumulator arrays have to beset up for different
sizes and orientations), a ‘centre’ (x,y) is chosen somewhere within the boundary of
the object.
The boundary is then traversed and after every step d alone the boundary the angle of
the boundary tangent with respect to horizontal is noted, and the x difference and y
difference of the boundary position from the centre point are also noted.
For every pixel I(x, y) in the edge-detected image, the gradient direction is found. The
accumulator array (same size as the image) is then incremented by 1 for each such
element.
Finally, the highest-valued elements in the accumulator array point to the possible
centres of the object in the image.
6.3 Bresenham’s Algorithm
Bresenham’s line algorithm is an efficient method for scan-converting straight lines in
that it uses only integer addition, subtraction, and multiplication by 2. As a very well known
fact, the computer can perform the operations of integer addition and subtraction very rapidly.
The computer is also time-efficient when performing integer multiplication and division by
powers of 2.
U
yk+1
d2
y
d1
D
yk
xk
xk + 1
The algorithm described in the following is a modified version of the Bresenham algorithm. It is commonly referred to as the midpoint line algorithm.
Figure 6.8 Midpoint algorithm
The equation of a straight line in 2-dimensional space can be written in an implicit form as
F(x, y) = ax + by + c = 0
= y +
x B From the slope-intercept form
dy
dx
we can bring it to the implicit form as
(cid:215) - (cid:215) dy x +
dx y Bdx = 0
So a = dy, b = - dx, c = Bdx
2 ). We have
2 ) + c
Suppose that point (xi, yi) has been plotted. We move xi to xi + 1. The problem is to
select between two pixels, U(xi + 1, yi + 1) and D(xi + 1, yi). For this purpose, we
consider the middle pixel M(xi + 1, yi + 1
2 )
d = F(M) = a(xi + 1) + b( yi + 1
choose U
If
choose D
choose either U or D, so choose U. d > 0 ,
d < 0 ,
d = 0 ,
- When D is chosen, M is incremented one step in the x direction. So
2 ) + c
dnew = F(xi +2, yi + 1
= a(xi + 2) + b(yi + 1
while
2 ) = a (xi + 1) + b (yi + 1
2 ) + c
dold = F(xi + 1, yi + 1
So the increment in d (denoted dD) is
dD = dnew - dold = a = dy
2 )
2 ) + c
- When U (xi + 1, yi + 1) is chosen, M is incremented one step in both directions:
dnew = F (xi +2, yi + 3
= a (xi + 2) + b( yi + 3
= dold + a + b
So the increment in d (denoted dU ) is
dU = a + b = dy - dx
In summary, at each step, the algorithm chooses between two pixels based on the sign of d. It
updates d by adding dD or dU to the old value.
2 ) and
First, we have the point (x1, y1). So M (x1 +1, y1 + 1
2 ) + c
F(M) = a(x1 + 1) + b (y1 + 1
= F(x1, y1 ) + a + b/2
Since F (x1 , y1) = 0, we have
d = d1 = dy - dx/2
In order to avoid a division by 2, we use 2d1 instead. Afterward, 2d is used. So, with d used in
place of 2d, we have
First set d1 = 2dy - dx
0 then If di ‡ xi+1 = xi + 1, yi+1 = yi + 1 and
di+1 = di + 2 (dy - dx)
If di < 0 then xi+1 = xi + 1, yi+1 = yi
di+1 = di + 2dy
The algorithm can be summarized as follows:
Midpoint Line Algorithm [Scan-convert the line between (x1, y1) and (x2, y2)]
dx = x2 - x1;
dy = y2 - y1;
d = 2*dy - dx; /* initial value of d */
dD = 2*dy; /* increment used to move D */
dU = 2*(dy - dx); /* increment used to move U */
x = x1;
y = y1 ;
Plot Point (x, y); /* the first pixel */
While (x < x1)
if d <0 then
d = d + dD; / * choose D */
x = x + 1;
else
d = d + dU; /* choose U */
x = x + 1;
y = y + 1;
endif
Plot Point (x, y); /* the selected pixel closest to the line */
EndWhile
Remark The described algorithm works only for those lines with slope between 0 and 1. It is
generalized to lines with arbitrary slope by considering the symmetry between the
various octants and quadrants of the xy-plane.
Example. Scan-convert the line between (5, 8) and (9, 11).
Since for the points, x < y, consequently the algorithm can apply. Here dy = 11 - 8 = 3, dx =
9 - 5 = 4
First d1 = 2dy - dx = 6 - 4 = 2 > 0 So the new point is (6, 9) and
d2 = d1 + 2 (dy - dx) = 2 + 2(- 1) = 0
(cid:222) the chosen pixel is (7, 10) and
d3 = d2 + 2 (dy - dx) = 0 +2(- 1) = - 2 < 0
the chosen pixel is (8, 10), then
d4 = d3 + 2dy = - 1 +6 = 5 > 0
The chosen pixel is (9, 11).
6.3.2 Circle incrementation
A circle is a symmetrical figure. Any circle-generating algorithm can take advantage
of the circle’s symmetry to plot eight points for each value that the algorithm calculates.
Eight-way symmetry is used by reflecting each calculated point around each 45(cid:176) axis. For
example, if point 1 in Figure 6.9 were calculated with a circle algorithm, seven more points
could be found by reflection. The reflection is accomplished by reversing the x, y coordinates
as in point 2, reversing the x, y coordinates and reflecting about the y axis as in point 3,
reflecting about the y
y
(-2, 8)
(-y, x)
9
(2, 8)
(y, x)
› fi
(8, 2)
(x, y)
(-8, 2)
(-x, y)
x
9
‹ fl
(x, -y)
(8, -2)
(-x, -y)
(-8, -2)
(cid:176) ‡
(y, -x)
(2, -8)
(-y, -x)
(-2, -8)
† –
Figure 6.9 Eight-way symmetry of a circle.
axis as in point 4, switching the signs of x and y as in point 5, reversing the x, y coordinates,
reflecting about the y axis and reflecting about the x axis as in point 6, reversing the x, y
coordinates and reflecting about the y axis as in point 7, and reflecting about the x axis as in
point 8.
To summarize:
P1 = (x, y)
P2 = (y, x)
P3 = (- y, x)
P4 = (- x, y) P5 = (- y, - x)
P1 = (- y, - x)
P7 = (y, - x)
P8 = (x, - y)
(i) Defining a Circle
There are two standard methods of mathematically defining a circle centered at the
origin. The first method defines a circle with the second-order polynomial equation (see
Figure 6.10).
y2 = r2 - x2
2
2
where x =
y =
r = the x coordinate
the y coordinate
the circle radius
- x r
With this method, each x coordinate in the sector, from 90 to 45(cid:176) , is found by stepping
for each step of x.
x from 0 to r / 2 , and each y coordinate is found by evaluating
This is a very inefficient method, however, because for each point both x and r must be
squared and subtracted from each other; then the square root of the result must be found.
y
y
2
2
The second method of defining a circle makes use of trigonometric functions (see Figure 6.11):
P=(r cos q , r sin q )
y
r sin q
r
q
= - P ( ,
x r x )
x
r cos q
x
x
Fig. 6.11 Circle defined with trigonometric Fig. 6.10 Circle defined with a second-
degree polynomial equation. functions.
x = r cosq y = r sinq
where q =
r =
x =
y = current angle
circle radius
x coordinate
y coordinate
By this method, q is stepped from q to p
/ 4, and each value of x and y is calculated.
However, computation of the values of sinq and cosq is even more time-consuming than the
calculations required by the first method.
(ii) Bresenham’s Circle Algorithm
If a circle is to be plotted efficiently, the use of trigonometric and power functions
must be avoided. And as with the generation of a straight line, it is also desirable to perform
the calculations necessary to find the scan-converted points with only integer addition,
subtraction, and multiplication by powers of 2. Bresenham’s circle algorithm allows these
goals to be met.
Scan-converting a circle using Bresenham’s algorithm works are follows. If the eight-
way symmetry of a circle is used to generate a circle, points will only have to be generated
through a 45(cid:176) angle. And, if points are generated from 90 to 45(cid:176) , moves will be made only in
the +x and -y directions (see Figure 6.12).
-y
45(cid:176)
+x
Figure 6.12 Circle scan-converted with Bresenham’s algorithm.
The best approximation of the true circle will be described by those pixels in the raster
that fall the least distance from the true circle. Examine Figures 6.13(a) and 6.13(b). Notice
that if points are generated from 90 and 45(cid:176) , each new point closest to the true circle can be
found by taking either of two actions: (1) move in the x direction one unit or (2) move in the x
direction one unit and move in the negative y direction one unit. Therefore, a method of
selecting between these two choices is all that is necessary to find the points closest to the true
circle.
Due to the 8-way symmetry, we need to concentrate only on the are from (0, r) to
,
(r / 2 r / 2 ) . Here we assume r to be an integer.
Suppose that P(xi, yi) has been selected as closest to the circle. The choice of the next pixel is between U and D (Fig.2.8).
Let F(x, y) =
F(x, y) =
x2 + y2 - r2. We know that
0
> 0
< 0 then (x, y) lies on the circle
then (x, y) is outside the circle
then (x, y) is inside the circle
Let M be the midpoint of DU. If M is outside then pixel D is closer to the circle, and
if M is inside, pixel U is closer to the circle.
2 )2 - r2
Let dold
= F(xi+1, yi - 1
2 )
= (xi + 1)2 + (yi - 1
* If dold < 0, then U (xi+1, yi) is chosen and the next midpoint will be one increment over x.
Thus
dnew = F(xi+2, yi - 1
2 )
= dold + 2xi + 3
The increment in d is
dU = dnew - dold = 2xi + 3
* 0, M is outside the circle and D is chosen. The new midpoint will be one If dold ‡
increment over x and one increment down in y:
dnew = F (xi + 2, yi - 3
2 )
= dold + 2xi - 2yi + 5
The increment in d is therefore
dD = dnew - dold = 2(xi - yi ) + 5
Since the increments dU and dD are functions of (xi , yi), we call point P(xi, yi) the point of
evaluation.
2 ) and so
Initial point : (0, r). The next midpoint lies at (1, r- 1
4 - r
2 ) = 1 + (r - 1
2 )2 - r2 = 5
1
4 . So the initials value of h is 1 - r
To avoid the fractional initialization of d, we take h = d - 1
and the comparison d < 0 becomes h < -
4 . However, since h starts out with an integer value
and is incremented with integer values (dU and dD), we can change the comparison to h < 0.
Thus we have an integer algorithm in terms of h. It is summarized as follows:
(0, r)
U(xi + 1, yi )
P(xi, yi)
(r /
,
2 r /
2
)
·
M
O
D(xi +1, yi - 1)
(b)
(a)
F(1, r - 1
Figure 6.13 Bresenham’s Circle Algorithm (Midpoint algorithm)
Bresenham Midpoint Circle Algorithm
h = 1 - r ; /*initialization */
x = 0;
y = r;
Plot Point (x, y);
While y > x
if h < 0 then /* Select U */
dU = 2*x + 3;
h = h + dU;
x = x + 1;
else /* Select D */
dD = 2*(x - y) + 5;
h = h - dD;
x = x + 1;
y = y - 1;
endif
End While
(iii) Second-order differences
If U is chosen in the current iteration, the point of evaluation moves from (xi, yi ) to
(xi+1, yi ). The first-order difference has been calculated as
++
(2
3)1
=¢
Ud
ix
dU = 2xi + 3
. Thus the second-order difference is At point (xi + 1, yi ), this will be
U
d
2=
¢=
d
U
U
- D
¢d D = 2(xi +1- yi ) + 5. Thus the Similarly, dD at (xi, yi ) is 2(xi - yi )+5 and at (xi +1, yi ) is
second-order difference is
D
D
- = ¢
D D d d = 2
If D is chosen in the current iteration, the point of evaluation moves from (xi, yi ) to (xi
+1, yi -1). The first-order differences are
2(
d
+
5)
-
d
y
i
-+
(1
)]1
=+
5
(
2
++
54)
y
i
y
i
x
i
=
D
=¢
[ 2
D
=
x
i
x
i
+
d
x
U
2
=
3
++
d
i
2(
x
3)1
- -
i
U
¢
Thus the second-order differences are
= D D U 2, D = 4
So the revised algorithm using the second-order differences is as follows:
D U = 3, D D = 5 - 2r, plot point (x, y) (1) h = 1 - r, x = 0 , y = r ,
(initial point) (2) Test if the condition y = x is reached.
It not then
If (3) h < 0 : select U
= x + 1
x
= h + D U
h
D U = D U + 2
D D = D D + 2
else : select D
x
= x + 1
= y - 1
y
= h + D D
h
D U = D U + 2
D D = D D + 4
end if
plot point (x, y)
6.4 Using interest point
The previous chapter described how interest points might be discovered from an image. From
these, it is possible to determine whether the object being viewed is a “known” object. Here
the two-dimensional problem, without occlusion (objects being covered up by other objects),
is considered. Assume that the interest points from the known two dimensional shape are held
on file in some way and that the two-dimensional shape to be identified has been processed by
the same interest points that now have to be compared with a known shape. We further
assume that the shape may be have been related, scaled, and/or translated from the original
known shape. Hence it is necessary to determine a matrix that satisfies:
discovered interest point = known shape interest point · M
or D = KM
where M is two-dimensional transformation matrix of the form
(cid:246) (cid:230) a b 0 (cid:247) (cid:231) (cid:247) (cid:231) c d 0 (cid:247) (cid:231) e f 1 ł Ł
and the interest point sets are of the form
2
...
2
...
n
n
(cid:246) (cid:230) 1 x
1 y
1 (cid:247) (cid:231) (cid:247) (cid:231) x y 1 (cid:247) (cid:231) ... (cid:247) (cid:231) (cid:247) (cid:231) x y 1 ł Ł
The matrix M described above does not allow for sheering transformations because this is
essentially a three-dimensional transformation of an original shape.
There is usually some error in the calculations of interest point positions so that
D = K M + e
2
and the purpose is to find M with the largest error and then determine whether that error is
small enough to indicate that the match is correct or not. A good approach is to use a least-
squares approximation to determine M and the errors, i.e. minimize
2 + y1
F(D-KM) where F(Z) = x1
2
This gives the following normal equations:
(cid:246) (cid:230) (cid:246) (cid:230) (cid:246) (cid:230) x x a xX (cid:229) (cid:229) (cid:229) (cid:229) (cid:247) (cid:231) (cid:247) (cid:231) (cid:247) (cid:231) xy
2 = = · (cid:247) (cid:231) (cid:247) (cid:231) (cid:247) (cid:231) xy y y c yX or La (cid:229) (cid:229) (cid:229) (cid:229) s
1 (cid:247) (cid:231) (cid:247) (cid:231) (cid:247) (cid:231) x y n e X (cid:229) (cid:229) (cid:229) ł Ł ł Ł ł Ł
2
and
2
(cid:246) (cid:230) (cid:246) (cid:230) (cid:246) (cid:230) x x b xY (cid:229) (cid:229) (cid:229) (cid:229) (cid:247) (cid:231) (cid:247) (cid:231) (cid:247) (cid:231) xy
2 = = · (cid:247) (cid:231) (cid:247) (cid:231) (cid:247) (cid:231) y or Lb s (cid:229) (cid:229) (cid:229) (cid:229) (cid:247) (cid:231) (cid:247) (cid:231) (cid:247) (cid:231) xy
x y
y n d
f yY
Y (cid:229) (cid:229) (cid:229) ł Ł ł Ł ł Ł
If the inverse of the square L matrix is calculated, then the values for a to f can be evaluated
and the error determinated. This is calculate as
L-1L a = L-1 s1 and L-1L b = L-1 s2
Resulting in
a = L-1s1 and b = L-1s2.
6.5 Problems
There are some problems with interest point. First, coordinates must be paired beforehand.
That is, there are known library coordinates, each of which must correspond to correct
unknown coordinate for a match to occur. This can be done by extensive searching, i.e. by
matching each known coordinate with each captured coordinate, all possible permutations
have to be considered. For example, consider an interest point algorithm that delivers five
interest points for a known objects. Also let there be N images, each containing an unknown
object, the purpose of the exercise being to identify if any or all of the images contain the
known object.
A reduction on the search can be done by eliminating all those images that do not have five
interest points. If this leaves n images there will be b x 5! = 120n possible permutations to
search. One search reduction method is to order the interest points. The interest operator itself
may give a value which can place that interest point at a particular position in the list.
Alternatively, a simple sum of the brightness of the surrounding pixels can be used to give a
position. Either way, if the order is known, the searches are reduced from 0(n x i!) to 0(n),
where i is the number of interest points in the image. The second problem is that the system
cannot deal with occlusion or part views of objects, nor can it deal with three-dimensional
objects in different orientations.
6.6 Exercises
6.6.1 Using standard graph paper, perform a straight line Hough transform on the binary
pixels array shown in the following figure transforming into (m,c) space.
Figure 6.8 Binary array
6.6.2 A library object has the following ordered interest point classification
{(0,0), (3,0), (1,0), (2,4)}
Identify, using the above technique, which of the following two sets of interest points
represent a transition, rotation, and/or scaling of the above object:
{(1,1), (6,12), (2,5), (12,23)}
{(1,3), (1,12), (-1,8), (3,6)}
Check your answer by showing that a final point maps near to its corresponding known point.