Description 
GETMIDPOINTCIRCLE return the x,y pixel coordinates of a circle
[x y] = getmidpointcircle(x0, y0, radius) returns the pixel coordinates
of the circle centered at pixel position [x0 y0] and of the given integer
radius. The midpoint circle algorithm is used for computation
(http://en.wikipedia.org/wiki/Midpoint_circle_algorithm).
This function is aimed at image processing applications, where the
integer pixel coordinates matter, and for which one pixel cannot be
missed or duplicated. In that view, using rounded trigonometric
coordinates generated using cosine calls are inadequate. The midpoint
circle algorithm is the answer.
Accent is made on performance. We compute in advance the number of point
that will be generated by the algorithm, to preallocate the coordinates
arrays. I have tried to do this using a MATLAB class implementing the
iterator pattern, to avoid computing the number of points in advance and
still be able to iterate over circle points. However, it turned out that
repeated function calls is extremely expansive, and the class version of
this function is approximately 1000 times slower. With this function, you
can get the pixel coordinates of a circle of radius 1000 in 0.16 ms, and
this time will scale linearly with increasing radius (e.g. it takes
0.16 s for a radius of 1 million).
Also, this functions ensure that sorted coordinates are returned. The
midpoint algorithm normally generates a point for the 8 circles octants
in one iteration. If they are put in an array in that order, the [x y]
points will jump from one octant to another. Here, we ensure that they
are returned in order, starting from the top point, and going clockwise.
EXAMPLE
n_circles = 20;
color_length = 100;
image_size = 128;
max_radius = 20;
I = zeros(image_size, image_size, 3, 'uint8');
colors = hsv(color_length);
for i = 1 : n_circles
x0 = round( image_size * rand);
y0 = round( image_size * rand);
radius = round( max_radius * rand );
[x y] = getmidpointcircle(x0, y0, radius);
index = 1 ;
for j = 1 : numel(x)
xp = x(j);
yp = y(j);
if ( xp < 1  yp < 1  xp > image_size  yp > image_size )
continue
end
I(xp, yp, :) = round( 255 * colors(index, :) );
index = index + 1;
if index > color_length
index = 1;
end
end
end
imshow(I, []);
