Hough Transform For Circle Detection - with known radius (II)

This example doing the same as the previous example, with the difference that this example does not use any "for loops" in the code. The vectorized code make the algorithm runs faster but with the trade-off of high memory consumtion. Try to compare this with the previous example:

1. Reading image and the convert to binary image. After that, the location of the value '1' is found using 'find' function. (assume the image is already preprocess, if not, the 'edge' function might help)

clear all;clc;
I = imread('pic21.bmp');
I =im2bw(I);
[y,x]=find(I);
[sy,sx]=size(I);
imshow(I);

2. Find all the require information for the transformatin. the 'totalpix' is the numbers of '1' in the image.

totalpix = length(x);

3. Preallocate memory for the Hough Matrix. Try to play around with the R, or the radius to see the different results.

HM = zeros(sy*sx,1);
R = 36;
R2 = R^2;

4. Performing Hough Transform. Notice that no "for-loop" in this portion of code.

a. Preparing all the matrices for the computation without "for-loop"

b = 1:sy;
a = zeros(sy,totalpix);

y = repmat(y',[sy,1]);
x = repmat(x',[sy,1]);

b = repmat(b',[1,totalpix]);

b. The equation for the circle

a = (round(x - sqrt(R2 - (y - b).^2)));

c. Removing all the invalid value in matrices a and b

b = b(imag(a)==0 & a>0);
a = a(imag(a)==0 & a>0);
ind = sub2ind([sy,sx],b,a);

d. Reconstruct the Hough Matrix

val = ones(length(ind),1);
data=accumarray(ind,val);
HM(1:length(data)) = data;
HM2 = reshape(HM,[sy,sx]);

5. Showing the Hough Matrix
cla
imshow(HM2,[]);

6. Finding the location of the circle with radius of R

[maxval, maxind] = max(max(HM2));
[B,A] = find(HM2==maxval);
imshow(I); hold on;
plot(mean(A),mean(B),'xr')