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)
I = imread('aaa.png'); //for this image, simply right click from the image above and save as 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.
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 = ones(sy,1).*.y;
y = ones(sy,1).*.y;
x = ones(sy,1).*.x;
b = ones(1,totalpix).*.b';
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 & real(a)>0);
a = a(imag(a)==0 & real(a)>0);
ind = sub2ind([sy,sx],b,a);
d. Reconstruct the Hough Matrix
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;
HM(1:length(data)) = data;
HM2 = matrix(HM,[sy,sx]);
5. Showing the Hough Matrix
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);
plot(mean(A),sy-mean(B),'xr');