Principles:

1) Edge detection (gradient based) necessary for reduced hough transform computation

2) Initializating hough paramter space, i.e., 3D accumulator matrix A (x, y , r)

3) Looping through all circles centered at the foreground pixels to obtain circle point indices

4) Finalizing the accumulator matrix consisting of voting number (futher normalization is helpful)

5) Finding the hough peaks (local maxima of the accumulator matrix)

The positions of the local maxima are corresponding to the circle centers in the original space.

Implementation:

[circles, A] = chough(bw, 75:100, 'normalize', 1, 'npeaks', 1, 'display', 1); %self-defined function

Comments： hough transform is computationally expensive, especially in case of a wide range of radii to be looped over, and finding the local maximum in 3D space is one of rating-limiting steps.

Reference: https://en.wikipedia.org/wiki/Circle_Hough_Transform