function img = BuildCircularityTestImage(ver,options)
%function img = BuildCircularityTestImage(ver,options)
%
if(nargin < 1), ver = 1;, end
if(nargin < 2), options = [];, end
d.RandSeed = ceil(rand*100000);
options = ResolveMissingOptions(options,d);
seedrand(options.RandSeed);
d.bUseCircularBoundary = 1;
d.NoiseType1 = 2; % 0 for none, 1 for white, 2 for 1/f noise
d.NoiseType2 = 2;% 0 for none, 1 for white, 2 for 1/f noise
d.fNoiseParam1 = -3; % -2 is hard -3 is medium -4 is easy
d.fNoiseParam2 = -3;% -2 is hard -3 is medium -4 is easy
d.xc = [];
d.yc = [];
d.r = [];
seedrand(options.RandSeed + 1);
d.Bias1 = 0;
seedrand(options.RandSeed + 2);
d.Bias2 = 0;
d.Dim = 256;
options = ResolveMissingOptions(options,d);

switch ver
  case 1
  switch options.NoiseType1
    case 0
    im1 = zeros(options.Dim);
    case 1
    im1 = rand(options.Dim);
    case 2
    im1 = Noise1OverF(options.Dim,options.fNoiseParam1);
  end
  switch options.NoiseType2
    case 0
    im2 = zeros(options.Dim);
    case 1
    im2 = rand(options.Dim);
    case 2
    im2 = Noise1OverF(options.Dim,options.fNoiseParam1);
  end
  im1 = (im1 - min(im1(:))) / (max(im1(:)) - min(im1(:)));
  im2 = (im2 - min(im2(:))) / (max(im2(:)) - min(im2(:)));
  im1 = im1 + options.Bias1;
  im2 = im2 + options.Bias2;
  Mask = BuildMask(options);
  img = zeros(options.Dim);
  img(find(Mask)) = im1(find(Mask));
  img(find(~Mask)) = im2(find(~Mask));
end

function mask = BuildMask(options)
  mask = zeros(options.Dim);
  if(options.bUseCircularBoundary)
    % pick a center
    if(isempty(options.xc))
      seedrand(options.RandSeed + 3);
      xc = rand* options.Dim;
    else
      xc = options.xc;
    end
    if(isempty(options.yc))
      seedrand(options.RandSeed + 4);
      yc = rand* options.Dim;
    else
      yc = options.yc;
    end
    % pick a radius
    if(isempty(options.r))
      seedrand(options.RandSeed + 5);
      r = (0.8*rand+.2) * options.Dim / 2;
    else
      r = options.r;
    end
    % build the mask
    for i = 1:size(mask,1);
      for j = 1:size(mask,2);
        d = sqrt((i - xc)^2 + (j - yc)^2);
	if(d < r)
	  mask(i,j) = 1;
	end
      end
    end
  else  %linear boundary
    % pick a point
    seedrand(options.RandSeed + 6);
    xc = (0.8*rand+.1)* options.Dim;
    seedrand(options.RandSeed + 7);
    yc = (0.8*rand+.1)* options.Dim;
    % pick a slope
    seedrand(options.RandSeed + 8);
    th = atan2(rand-.5,rand-.5);
    % build the mask
    for i = 1:size(mask,1);
      for j = 1:size(mask,2);
        mask(i,j) = atan2((i-xc),(j-yc));
      end
    end
    mask = abs(mod(mask - th,2*pi)) < (pi);
  end
  
