function [pts, dists] = randomPtsOnR(num)
% produce NUM random points on the R
pts = zeros(num, 2);
dists = zeros(num, 1);
idx = 1;

r1 = rect_perim(0.2, 0.1, 0.3, 0.90);
r2 = rect_perim(0.73, 0.091, 0.76, 0.11);
c1 = 2 * bezier_length(0.37, 0.36, 0.30, 0.3832, 0.32, 0.46, .40, ...
    0.43, 100); 
c2 = 2 * bezier_length(0.40, 0.43, 0.64, 0.34, 0.5, 0.1, 0.75, ...
    0.1, 100); 
c3 = 2 * bezier_length(0.3, 0.89, 0.9, 0.90, 0.75, 0.25, 0.37, ...
    0.36, 180); 

tot = sum([r1 r2 c1 c2 c3]);

for idx = 1:num,
  p = [];
  loc = tot * rand;
  if (loc < r1),
    p = rand_rect(0.2, 0.1, 0.3, 0.90);
  end
  loc = loc - r1;
  if (isempty(p) & (loc < r2)),
    p = rand_rect(0.73, 0.091, 0.76, 0.11);
  end
  loc = loc - r2;
  if (isempty(p) & (loc < c1)),
    p = rand_bezier(0.37, 0.36, 0.30, 0.3832, 0.32, 0.46, .40, ...
        0.43, 100, 0.01); 
  end
  loc = loc - c1;
  if (isempty(p) & (loc < c2)),
    p = rand_bezier(0.40, 0.43, 0.64, 0.34, 0.5, 0.1, 0.75, ...
        0.1, 100, 0.01); 
  end
  loc = loc - c2;
  if isempty(p),
    if loc > c3,
      'foo!'
    end
    p = rand_bezier(0.3, 0.89, 0.9, 0.90, 0.75, 0.25, 0.37, ...
        0.36, 180, 0.01); 
  end

  
  pts(idx,:) = p;
  dists(idx) = Rdist(p);
end


function p = rect_perim(x1, y1, x2, y2)
p = 2 * (abs(x1 - x2) + abs(y1 - y2));


function l = bezier_length(x1, y1, x2, y2, x3, y3, x4, y4, steps)
p_prev = [x1 y1];
l = 0;
for i = 1:steps,
  t = i / steps;
  p_next = [bez4(x1, x2, x3, x4, t) bez4(y1, y2, y3, y4, t)];
  l = l + norm(p_prev - p_next);
  p_prev = p_next;
end


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;


function p = rand_rect(x1, y1, x2, y2)
loc = rand * (abs(x1 - x2) + abs(y1 - y2));
if (loc < abs(y1 - y2)),
  if (rand > 0.5),
    p = [x1 linear(y1, y2, rand)];
  else
    p = [x2 linear(y1, y2, rand)];
  end
else
  if (rand > 0.5),
    p = [linear(x1, x2, rand) y1];
  else
    p = [linear(x1, x2, rand) y2];
  end
end


function p = rand_bezier(x1, y1, x2, y2, x3, y3, x4, y4, ...
    steps, radius)
t1 = rand;
t2 = t1 + 1 / steps;
p1 = [bez4(x1, x2, x3, x4, t1) bez4(y1, y2, y3, y4, t1)];
p2 = [bez4(x1, x2, x3, x4, t2) bez4(y1, y2, y3, y4, t2)];
if (rand < 0.5),
  temp = p1;
  p1 = p2;
  p2 = temp;
end
dxdt = p2 - p1;
% rotate 90 degrees to the right
dxdt = [dxdt(2) -dxdt(1)];
p = p1 + radius * dxdt / norm(dxdt);