function [simData, X, simSeg] = demo_simulateDonuts(P2D, optsin)
imSize = [1280 1280];

%% PARSE INPUTS
opts.nFrames = 1; tooltips.nFrames = 'number of frames to simulate';
opts.rOuter = 30; tooltips.rOuter = 'outer donut radius';
opts.rInner = 18; tooltips.rInner = 'inner donut radius';
opts.brightness = 1000; tooltips.brightness = 'photons/donut';
opts.nDonuts = 50; tooltips.nDonuts = 'number of donuts';
opts.darkRate = 1; tooltips.darkRate= 'dark noise level';
%Update with user inputs
if nargin>1 %if the user specified options
    for field = fieldnames(optsin)'
         opts.(field{1}) = optsin.(field{1});
    end
else %no user options; use a GUI
      opts = optionsGUI(opts, tooltips);
end
opts.isSimulation = true;

%% INITIALIZE
SEo= strel('disk', opts.rOuter, 0);
SEi= strel('disk', opts.rInner, 0);
censor = zeros(4*opts.rOuter+1); censor(2*opts.rOuter+1,2*opts.rOuter+1) = true;
censor = imdilate(imdilate(censor,SEo),SEo);
[xq,yq] = meshgrid(-639:640);
SLM = xq.^2+yq.^2 < 635^2;
valid = imerode(imerode(SLM, SEo), SEo);


%% SIMULATION
%Generate Donuts
segIM = zeros(imSize);
inds = find(valid(:));
inds = inds(randperm(length(inds)));
n = 0; ix = 1;
while n<opts.nDonuts && ix<length(inds)
    if valid(inds(ix))
        [rr,cc] = ind2sub(imSize, inds(ix));
        segIM(rr,cc) = n+1;
        valid(rr+[-2*opts.rOuter:2*opts.rOuter],cc+[-2*opts.rOuter:2*opts.rOuter]) = valid(rr+[-2*opts.rOuter:2*opts.rOuter],cc+[-2*opts.rOuter:2*opts.rOuter]) & ~censor;
        n=n+1;
    end
    ix = ix+1;
end
segIM = imdilate(segIM, SEo) .* double(~imdilate(segIM,SEi));

%make the sparse segmentation matrix
ixs = find(segIM(:)>0);
v = segIM(ixs);
simSeg = sparse(ixs,v, ones(size(v)), prod(imSize), opts.nDonuts);

%simulate activity
X = 1+4*rand(opts.nDonuts,opts.nFrames);

%Simulate Imaging
bFac = mean(sum(P2D*simSeg,1));
rates = opts.brightness/bFac*full(P2D*simSeg*X)+opts.darkRate;
simRef = nan([imSize opts.nFrames]);
for frame = 1:opts.nFrames
    simRef(:,:,frame) = reshape(simSeg*X(:,frame), imSize);
    simData.frames(frame).pmtData = poissrnd(rates(:,frame));
end
end