function dist = Rdist(x, y)
r1 = rect_dist(x, y, 0.2, 0.1, 0.3, 0.90);
r2 = rect_dist(x, y, 0.73, 0.091, 0.76, 0.11);
c1 = bezier_dist(x, y, 0.37, 0.36, 0.30, 0.3832, 0.32, 0.46, .40, ...
    0.43, 100, 0.01); 
c2 = bezier_dist(x, y, 0.40, 0.43, 0.64, 0.34, 0.5, 0.1, 0.75, ...
    0.1, 100, 0.01); 
c3 = bezier_dist(x, y, 0.3, 0.89, 0.9, 0.90, 0.75, 0.25, 0.37, ...
    0.36, 180, 0.01); 

dist = max([r1, r2, c1, c2, c3]')';



function rd = rect_dist(x, y, x1, y1, x2, y2)
dx1 = x - x1;
dx2 = x2 - x;
dy1 = y - y1;
dy2 = y2 - y;
% preallocate
rd = dx1 - dx1 - 5;
xpos = (dx1 > 0) & (dx2 > 0);
ypos = (dy1 > 0) & (dy2 > 0);
indxs = find(xpos & ypos);
rd(indxs) = min(min(dx1(indxs), dx2(indxs)), min(dy1(indxs), ...
    dy2(indxs)));
indxs = find(xpos & ~ypos);
rd(indxs) = min(dy1(indxs), dy2(indxs));
indxs = find(ypos & ~xpos);
rd(indxs) = min(dx1(indxs), dx2(indxs));
indxs = find((~xpos) & (~ypos));
rd(indxs) = -hypot(min(abs(dx1(indxs)), abs(dx2(indxs))), ...
    min(abs(dy1(indxs)), abs(dy2(indxs))));

  
function hp = hypot(x, y)
hp = sqrt(x.*x + y.*y);


function bd = bezier_dist(x, y, x1, y1, x2, y2, x3, y3, x4, y4, ...
    steps, radius)
x1 = x1 - x;
x2 = x2 - x;
x3 = x3 - x;
x4 = x4 - x;
y1 = y1 - y;
y2 = y2 - y;
y3 = y3 - y;
y4 = y4 - y;

bd = sqdist(x1, y1);
for i = 1:steps,
  t = i / steps;
  bd = min(bd, sqdist(bez4(x1, x2, x3, x4, t), bez4(y1, y2, ...
      y3, y4, t)));
end
bd = radius - sqrt(bd);


function sd = sqdist(x, y)
sd = x.*x + y.*y;


function b4 = bez4(a, b, c, d, t)
a = linear(a,b,t);
b = linear(b,c,t);
c = linear(c,d,t);
a = linear(a,b,t);
b = linear(b,c,t);
b4 = linear(a,b,t);


function l = linear(a,b,t)
l = a + (b - a) * t;