function circle = makeCircle(Nx, Nz, cx, cz, radius, arc_angle, plot_circle)
%MAKECIRCLE     Create a binary map of a circle within a 2D grid.
%
% DESCRIPTION:
%       makeCircle creates a binary map of a circle or arc (using the
%       midpoint circle algorithm) within a two-dimensional grid (the
%       circle position is denoted by 1's in the matrix with 0's
%       elsewhere). A single pixel is taken as the circle centre thus the
%       total diameter will always be an odd number of pixels.
%
% USAGE:
%       circle = makeCircle(Nx, Nz, cx, cz, radius)
%       circle = makeCircle(Nx, Nz, cx, cz, radius, arc_angle)
%       circle = makeCircle(Nx, Nz, cx, cz, radius, arc_angle, plot_circle)
%
% INPUTS:
%       Nx, Nz          - size of the 2D grid [number of pixels]
%       cx, cz          - centre of the circle [pixel coordinates]
%       radius          - circle radius [number of pixels]
%
% OPTIONAL INPUTS:
%       arc_angle       - arc angle for incomplete circle [radians]
%                         (default = 2*pi)
%       plot_circle     - Boolean controlling whether the circle is plotted
%                         using imagesc (default = false)
%
% OUTPUTS:
%       circle          - 2D binary map of a circle
%
% ABOUT:
%       author          - Bradley Treeby
%       date            - 1st May 2009
%       last update     - 17th June 2009
%       
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
%
% See also makeCartCircle, makeDisc

% check for plot_circle input
if nargin < 7
    plot_circle = false;
end

% check for arc_angle input
if nargin < 6
    arc_angle = 2*pi;
elseif arc_angle > 2*pi
    arc_angle = 2*pi;
elseif arc_angle < 0
    arc_angle = 0;
end

% force integer values
Nx = round(Nx);
Nz = round(Nz);
cx = round(cx);
cz = round(cz);
radius = round(radius);

% define literals
MAGNITUDE = 1;

% create empty matrix
circle = zeros(Nz, Nx);

% initialise loop variables
x = 0;
z = radius;
d = 1 - radius;

% draw the first cardinal point
try 
    circle(cz, cx-z) = MAGNITUDE;
catch ME
    error('Circle must fit within the grid');
end

% draw the remaining cardinal points
pz = [cz, cz+z, cz-z];
px = [cx+z, cx, cx];
for point_index = 1:length(pz)
    
    % check whether the point is within the arc made by arc_angle
    if (atan2(pz(point_index) - cz, px(point_index) - cx) + pi) <= arc_angle
        circle(pz(point_index), px(point_index)) = MAGNITUDE;
    end
end

% loop through the remaining points using the midpoint circle algorithm
while ( x < z - 1 )
    
    x = x + 1;
    if ( d < 0 ) 
        d = d + x + x + 1;
    else 
        z = z - 1;
        a = x - z + 1;
        d = d + a + a;
    end
    
    % setup point indices
    pz = [x+cz, z+cz, z+cz, x+cz, -x+cz, -z+cz, -z+cz, -x+cz];
    px = [z+cx, x+cx, -x+cx, -z+cx, -z+cx, -x+cx, x+cx, z+cx];
    
    % loop through each point
    for point_index = 1:length(pz)
        
        % check whether the point is within the arc made by arc_angle
        if (atan2(pz(point_index) - cz, px(point_index) - cx) + pi) <= arc_angle
            circle(pz(point_index), px(point_index)) = MAGNITUDE;
        end
    end
end

% create the figure
if plot_circle
    figure;
    imagesc(circle, [-1 1]);
    colormap(getColorMap);
    axis image;
    xlabel('x-position [pixels]');
    ylabel('z-position [pixels]');
end