Intro

The Circle Hough transform is a pretty standard algorithm for circle detection in images. It works reasonably well, and the wikipedia article explains the basic steps.

Something that bothered me from day one was that it is surprisingly brute-force. The idea is basically for every pixel on an edge, try EVERY possible circle this pixel might be a part of. As a result most of the ‘votes’ in the accumulator probably won’t be near any circles.

In this post I’ll describe a small modification to the Circle Hough algorithm, in which we will also extract the edge direction from the the image to make better guesses for where the circle centers might be.

Setup

I’ll be using python and openCV. This post will have the important (image processing) parts of the code but you can also follow along with the original iPython notebook in Jupyter.

Preparation

We load the image and convert to single-channel (grayscale). We also define a few parameters for the algorithm, such as the range of circle radii to look for, and a threshold to discard ‘unlikely’ circles.

1
2
3
4
5
6
7
8
9
10
11
12
13
# circle radii to search for
MIN_RADIUS = 18
MAX_RADIUS = 38

# accumulator must have at least this many 'votes'
ACC_THRESHOLD = 20

img = cv.imread('coins.jpg')
W = img.shape[0]
H = img.shape[1]

# convert to grayscale
gray = cv.cvtColor(img, cv.COLOR_BGR2GRAY)

This is the image I am using:
png

Finding the edges and their direction

We first apply a small blur to reduce any noise/small artifacts that we don’t care about. Then we detect the edges using a sobel filter, however, we detect the horizontal and vertical components separately. A regular hough transform at this step would just use the full edge-detect result and get only the gradient magnitude.

1
2
3
4
5
6
# blur to reduce noise
blurred = cv.blur(gray, (3, 3))

# get horizontal and vertical gradients separately
edge_x = cv.Sobel(blurred, cv.CV_32F, 1, 0)
edge_y = cv.Sobel(blurred, cv.CV_32F, 0, 1)

png

The formula for sobel’s edge magnitude should be familiar: euclidian distance between the x & y components.

Getting the edge directions isn’t anything special either, it’s actually part of a the sobel operation. Since we have the x and y components, it’s just some trigonometry:

png

1
2
edges = (np.sqrt((edge_x**2 + edge_y**2)))
direction = np.arctan2(edge_y, edge_x)

By doing the above math for every pixel, we end up with a ‘map’ of the edge directions in addition to our map of the edge strength.

png

The Accumulator

Accumulator is just a fancy word for guesses. We are basically creating a list of every circle that could conceivably be in the image, and then trying to decide which ones are definitely in the image.

The circles can be described with x, y and radius, so our accumulator is a X*Y*(MAX_RADIUS - MIN_RADIUS) array. Every time we find a pixel on an edge, we add a ‘vote’ in the accumulator to the places where the center of that circle might be.

In a regular hough transform, we would vote for every point in a circle around the edge pixel we’re looking at, and then do that for every circle radius we are looking for.

In this version, we will only vote the points that are perpendicular to the edge direction for each radius.

The difference is illustrated in this diagram below. Each red dot is a vote in the accumulator.

png

The difference in number of guesses is pretty huge! Since the regular hough has to ‘draw’ a set of circles in the accumulator, it’s a quadratic (n2) operation with respect to the number of radii we’re checking. In comparison, guessing only the 2 perpendicular points is linear (2n), and we only need to check 2 extra points per radius.

This could be advantageous for checking a large range of circle sizes. In this particular example (min_radius = 18, max_radius = 38), the regular hough would add 3284 votes to the accumulator per point, where as the modified algorithm adds 40. Checking one extra radius up (max_radius=39) would make it 3502 (+220) vs 42 (+2).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# create an accumulator for each search radius
DEPTH = MAX_RADIUS - MIN_RADIUS
accumulator = np.zeros((W, H, DEPTH))

for x,y in nonzero(edges > 128):
for z in range(DEPTH):
r = MIN_RADIUS + z
# get the direction (tangent) of the edge
dx = int(round(math.cos(direction[x, y]) * r))
dy = int(round(math.sin(direction[x, y]) * r))

# 'vote' for the two possible circle centers
# that are perpendicular to the edge and 'r' pixels away
if 0 <= x + dy < W and 0 <= y + dx < H:
accumulator[x + dy, y + dx, z] += 1
if 0 <= x - dy < W and 0 <= y - dx < H:
accumulator[x - dy, y - dx, z] += 1

As you view the different layers of the accumulator, you’ll see that the values seem to converge on a few points. These are locations where the votes from the image edges at a particular radius ‘agreed’ on a center. Those locations will have higher values, and probably be our circle centers.

Drag the slider to browse the accumulator array Radius: 18

Finding the peaks

Now we have a 3-dimensional array with all possible circles and votes for how likely each of them is. To find the final locations of the circles, we need to find the ‘peaks’ of this volume. The x & y of the each peak will be the center, and the ‘z’ will give us the radius of that circle.

First, we ‘flatten’ the accumulator into a 2d image by taking the maximum value of each x & y location along the z-axis. The maximum value is stored in our flattened accumulator, but we also store the z-value of that local maximum, which we will later look up to get the circle radius.

1
2
3
4
5
6
7
8
9
max_acc = np.zeros((W, H))
max_radius = np.zeros((W, H))
# for each location/radius, keep value with the highest number of votes
# also save the corresponding circle radius in max_radius
for x in range(W):
for y in range(H):
index = np.argmax(accumulator[x, y])
max_acc[x, y] = accumulator[x, y, index]
max_radius[x, y] = index + MIN_RADIUS

The ‘Flattened Accumulator’ shows us the how many votes the most likely radius got at each location, and the second radius ‘buffer’ shows the actual radius that had the most votes at that location

png

Think of the flattened accumulator as a probability function, where the highest values are around locations with the highest probability of being a circle center. We need to filter out the exact locations of the local maxima, and those are our cirlces.

This is a pretty simple process: For each point, we check the points around it and if it doesn’t have the highest value in its neighborhood it’s eliminated. I like to define the neighborhood as an s*s square, where s is equal to the radius of the circle proposed at that location.

It’s like asking for each point, Are there any circles within the circle that would be here that are more likely?

Before that, we blur the flattened accumulator and apply a threshold to remove any noise, as well as exclude any candidate circles that have too few votes. Im my experience, blurring the accumulator also helps to better locate the actual centers of the circles.

1
2
3
4
5
6
7
8
9
10
11
max_acc = cv.blur(max_acc, (3,3))
peaks = max_acc * (max_acc > ACC_THRESHOLD)

# non-maximum suppression
for x, y in nonzero(peaks):
s = int(max_radius[x, y]) # size of neighborhood
# left, right, top, bottom of neighborhood
l, r, t, b = max(0, x-s), min(W, x+s), max(0, y-s), min(H, y+s)

if peaks[x, y] < peaks[l:r, t:b].max():
peaks[x, y] = 0

We get local maxima of our flattened accumulator, which should correspond to the circle centers in the image:
png

Final Result

With our peaks singled out, we look up the radius from max_radius for each peak and plot the circles to see the result:

1
2
3
4
5
result = np.copy(img)
# draw the found circles
for x, y in nonzero(peaks):
radius = int(max_radius[x, y])
cv.circle(result, (y, x), radius, (0,255,0), 2)

png