function circle_image = circle_mask( N, r2 )
% Create a NxN image where
% all points less than radius sqrt(r2) away from the center are set to 1,
% everywhere outside is set to 0.
% For odd N, the circle is centered on the center pixel;
% for even N, the circle is centered on the center 4 pixels.
% N is assumed to be integral; r2 can be any double.
% Pixels exactly radius r away from the center
% are considered to be "outside" the circle
% and are zeroed.
% Even N causes the diameter to always be even.
% Odd N causes the diameter to always be odd.
% For example, with r2=1, N=odd, the single center pixel is turned on.
% For example, with r2=1, N=even, the central 2x2 square is turned on.
% For example, with r2=2, N=even, the central 2x2 square is turned on.
% For example, with r2=2, N=3, the cross
%   [0 1 0;
%    1 1 1;
%    0 1 0]
%   is returned
%   even though the corner pixels are at (radius^2)=2 exactly.
% For example, for circle_mask(4, 4.5) returns
%   [0 1 1 0;
%    1 1 1 1;
%    1 1 1 1;
%    0 1 1 1].
% If there are no pixels less than sqrt(r2) away from the center --
% -- i.e., N = even and r2 <= 1/2, or N=odd and r2 <= 0 --
% -- then this returns the same matrix as zeros(N).
% Works for both even and odd N.
%
% See also CIRCLE.


% Change Log
% 1998-09-16:DAV: documented; tweaked implementation to make
%   documentation simpler.
% 1998-03-01:DAV: wrote David Cary <d.cary@ieee.org>

% FIXME: Probably using "meshgrid()" would make this a little simpler.

even = rem(N,2) == 0;
if(even)
    % N is even
    % get quarter mask centered just above and to left of upper-left pixel
    
    % make vector of [0.5 1.5 2.5 3.5 ....]
    index = 0.5 + (0:(N-1));
else
    % N is odd
    % get quarter mask centered on upper-left pixel
    index = 0:(N-1);
end
onerow = ones(1,N);
x = onerow' * index;
y = x';
% find distance from center
distance2 = (y.^2) + (x.^2);
quarter_mask = distance2 < r2;
% shift so it is centered in grid
% (centered on center pixel for odd N,
% centered halfway between the 2 center columns and rows for even N)
quarter_mask = fftshift(quarter_mask);
half_mask = quarter_mask + fliplr(quarter_mask);
full_mask = half_mask + flipud(half_mask);

% if size is odd, the center pixel now has value "4".
% Convert back to binary 0 and 1.
circle_image = 0 < full_mask;

% end circle_mask