function [mgn,z] = yinyang (z, opts, output_mgn)
%YINYANG  Give the signed margin of the Yin-Yang function.
%
%   Y = YINYANG(X, [ratio R orientation]) gives the output (1/-1) and
%   MGN = YINYANG(X, [ratio R], true) gives the minimal L2 distance of
%   each point in X to the boundaries. Outside of the normal Yin-Yang,
%   the boundary is the y-axis (for orientation 0).
%   X: each line is a point
%   ratio: the small disk (eye) radius vs. the big disk radius
%   R: the big disk radius
%
%   YINYANG('plot', [...]) plots the Yin-Yang symbol.

if nargin < 3, output_mgn = false; end
if nargin < 2, opts = []; end
opts = options([0.2 1 0], opts);
r = opts(1); R = opts(2); ori = opts(3);
assert(r < 0.5);
r = r * R; % radius of the eyes

if isstr(z) && strcmp(z, 'plot')
    gray = .4; % users can specify the "black" background
    if nargin >= 3, gray = output_mgn; end
    plot_yinyang(r, R, ori, gray);
    return;
end

if length(z) == 1, %% Generate random inputs
    assert(z == round(z));
    tz = 1.1 * (rand(2*z,2)*2-1);
    i = find(sum(tz.^2,2) < 1.1^2);
    assert(length(i) >= z);
    z = tz(i(1:z),:) * R;
    clear tz i
end

x = z(:,1); y = z(:,2);
if ori < 0, x = -x; ori = -ori; end
for i = 1:mod(ori,4)
    t = -x; x = y; y = t;
end

R2 = R / 2;
d1 = sqrt((x+R2).^2 + y.^2);
d2 = sqrt((x-R2).^2 + y.^2);
rd1 = r - d1;
d2r = d2 - r;
d1R2 = d1 - R2;
R2d2 = R2 - d2;

mgn = repmat(nan, length(x), 1);

i = (rd1 >= 0);                            % d1 <= r
mgn(i) = rd1(i);
i = (d2r <= 0);                            % d2 <= r
mgn(i) = d2r(i);

i = (d2r > 0) & (R2d2 >= 0) & (y <= 0);    % r < d2 & d2 <= R2 & y <= 0
mgn(i) = min(d2r(i), R2d2(i));
i = (rd1 < 0) & (d1R2 <= 0) & (y >= 0);    % r < d1 & d1 <= R2 & y >= 0
mgn(i) = max(d1R2(i), rd1(i));

c1 = -sqrt((x+R).^2 + y.^2);
c2 = sqrt((x-R).^2 + y.^2);
i = (abs(x) >= R);
c1(i) = y(i); c2(i) = y(i);

i = (d2r > 0) & (d1R2 > 0) & (y >= 0);     % r < d2 & R2 < d1 & y >= 0
mgn(i) = min(min(d2r(i), d1R2(i)), c2(i));
i = (rd1 < 0) & (R2d2 < 0) & (y <= 0);     % r < d1 & R2 < d2 & y <= 0
mgn(i) = max(max(rd1(i), R2d2(i)), c1(i));

assert(all(~isnan(mgn)));
if ~output_mgn
    mgn = sign(mgn);
end


function plot_yinyang (r, R, ori, gray)
hold on; axis equal; axis off;
gray = gray * [1 1 1];

t = [0:359 0]/180*pi;
ct = cos(t); st = sin(t);
ct2 = ct(1:(end-1)/2); st2 = st(1:(end-1)/2);

R2 = R / 2;
px = [r*ct+R2 r*ct-R2 R2-R2*ct2  R*ct2 -R2-R2*ct2 0];
py = [r*st    r*st      -R2*st2 -R*st2     R2*st2 0];
sep = length(ct);

if ori < 0, px = -px; ori = -ori; end
for i = 1:mod(ori,4)
    t = px; px = -py; py = t;
end

h = patch(px(1:sep), py(1:sep), gray);
set(h, 'EdgeColor', gray, 'LineWidth', 0.01);
h = patch(px(sep+1:end), py(sep+1:end), gray);
set(h, 'EdgeColor', gray, 'LineWidth', 0.01);
plot(R*ct, R*st, 'Color', gray);