circle_hough.ipynb

Download .ipynb

import math
import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt
from ipywidgets import interact

def nonzero(arr):
    return np.transpose(np.nonzero(arr))

# 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)
plt.imshow(gray, cmap='gray')
plt.title('Original Image'), plt.axis('off')
plt.show()

# 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)

plt.subplot(1, 2, 1)
plt.imshow(edge_x, cmap='gray')
plt.title('Horizontal Edges'), plt.axis('off')
plt.subplot(1, 2, 2)
plt.imshow(edge_y, cmap='gray')
plt.title('Vertical Edges'), plt.axis('off')
plt.show()

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

plt.subplot(1, 2, 1)
plt.imshow(edges, cmap='gray')
plt.title('magnitude'), plt.axis('off')
plt.subplot(1, 2, 2)
# mask directions that are not on an edge for clarity
plt.imshow(direction * (edges > 128)) 
plt.title('direction'), plt.axis('off')
plt.show()

# 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
                    
def view_acc(radius):
    plt.imshow(accumulator[:, :, radius - MIN_RADIUS], cmap='gray')
    plt.axis('off')
    plt.savefig('acc'+str(radius)+'.jpg', bbox_inches='tight', pad_inches=0)
    plt.show()
interact(view_acc, radius=(MIN_RADIUS, MAX_RADIUS -1))

<function __main__.view_acc>

max_acc = np.zeros((W, H))
max_radius = np.zeros((W, H))
# for each location/radius, keep the accumulator radius with the 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

max_acc = cv.blur(max_acc, (3,3))
        
plt.subplot(1, 2, 1); plt.imshow(max_acc, cmap='gray')
plt.title('Flattened Accumulator'), plt.axis('off')
plt.subplot(1, 2, 2); plt.imshow(max_radius)
plt.title('Radius'); plt.axis('off')
plt.show()

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

plt.imshow(peaks > 0, cmap='gray')
plt.title('Accumulator Peaks'), plt.axis('off')
plt.show()

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)
  
plt.imshow(result)
plt.title('Result'); plt.axis('off')
plt.show()