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

The detection of circle with unknown radius has been shown in previous post - Hough Transform For Circle Detection - with unknown radius (I) , by using a lot of "for" loops.This example shows the same function with "vectorized" code, which run 25% faster than the previous code, in a Pentium M 1.5 GHz computer with 512 M RAM.Let's see how to perform this:


1. Reading image and the convert to binary image.

I = imread('aaa.png');
I =im2bw(double(I),0.5);
[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.

stacksize('max');
HM = zeros(sy*sx*50,1);
R = 1:50;
R2 = zeros(1,1,length(R)); // Extra line for scilab
R2(1,1,:) = R.^2;
//R2 = repmat(R2,[sy,totalpix,1]);
R2 = ones(sy,totalpix,1).*.R2;
//y = repmat(y',[sy,1,50]);
y = ones(sy,1,50).*.y;
//x = repmat(x',[sy,1,50]);
x = ones(sy,1,50).*.x;

b = 1:sy;
//b = repmat(b',[1,totalpix,50]);
b = ones(1,totalpix,50).*.b';

a = zeros(sy,totalpix);

4. Performing Hough Transform. Notice the accumulator is located in the inner for loop. This portion of codes will map the original image to the a-b domain.

a. Find the a parameter

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

b. Find the location of ‘a’ and ‘b’, as well as the valid value both.

[xx,yy,zz]=find((imag(a)==0 & real(a)>0));
b = b(imag(a)==0 & real(a)>0);
a = a(imag(a)==0 & real(a)>0);

c. Find the indices and remap it to different layers

ind = sub2ind([sy,sx],b,a);
ind = ind + ((ceil(zz)-1)*sx*sy)'; // interesting

d. Construct the Hough Matrix of 50 layers

val = ones(length(ind),1);

data = zeros(max(real(ind)),1);
ind = real(ind);
for cnt = 1:length(ind)
data(ind(cnt)) = data(ind(cnt)) + 1;
end

HM(1:length(data)) = data;
HM2 = matrix(HM,[sy,sx,50]);

5. Find for the maximum value for each layer, or in other words, the layer with maximum value will indicate the correspond R for the circle.

for cnt = 1:50
H(cnt) = (max(HM2(:,:,cnt)));
end
plot(H,'*-');

6. Extract the information from the layer with maximum value, and overlap with the original image.

[maxval, maxind] = max(H);
[B,A] = find(HM2(:,:,maxind)==maxval);
imshow(double(I));
plot(mean(A),sy-mean(B),'xr')
xnumb(mean(A),sy-mean(B),maxind,100)
h = gce();
h.font_foreground = 255;