**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.**

HM = zeros(sy,sx,50);

R = 1:50;

R2 = R.^2;

sz = sy*sx;

**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.**

for cnt = 1:totalpix

for cntR = 1:50

b = 1:sy;

a = (round(x(cnt) - sqrt(R2(cntR) - (y(cnt) - [1:sy]).^2)));

b = b(imag(a)==0 & real(a)>0);

a = a(imag(a)==0 & real(a)>0);

ind = sub2ind([sy,sx],b,a);

HM(sz*(cntR-1)+ind) = HM(sz*(cntR-1)+ind) + 1;

end

end

**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(max(HM(:,:,cnt)));

end

plot(H,'*-');