demoPR.m

%
% demoPR - using PR with random overlapping ellipses 
%
% NOTE - for working with ground truth, etc., see EllipseSimulation folder
% in this project.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (c) 2016, Drexel University
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
% 
% * Redistributions of source code must retain the above copyright notice, this
%   list of conditions and the following disclaimer.
% 
% * Redistributions in binary form must reproduce the above copyright notice,
%   this list of conditions and the following disclaimer in the documentation
%   and/or other materials provided with the distribution.
% 
% * Neither the name of PixelRep nor the names of its
%   contributors may be used to endorse or promote products derived from
%   this software without specific prior written permission.
% 
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
% FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

K=3; % the number of  ellipses (cells) we're fitting
cmap=hsv(K);
% generate a random image. you can use your own image here instead...
bw=GetRandomEllipseImage(K); 

figure(1);
hold off
imshow(bw)
RMAX=K*30;
rLim = [500-RMAX,500+RMAX];
xlim(rLim); ylim(rLim)
hold on

% NOTE -- important -- bw only has 1 connected component! If bw has multiple
% connected components, process each one separately:
% [L num]=bwlabel(bw);
% for n=1:num
%   ptsReplicated = PixelReplicate(L,n);
%   objPR = fitgmdist(ptsReplicated,K,'replicates',5);
%   ....
% end

% For more than a million replicate points use sampling at 10% instead.
numPoints = sum(round(bwdist(~bw)));
if ( numPoints < 1e6  )
    ptsReplicated = PixelReplicate(bw);
else
    ptsReplicated = PixelReplicateSample(bw,0.1);
end

% fit gmm to PR points
% NOTE - more replicates is more accurate fit, but takes longer. you can
% spmd this, or adjust as needed...
objPR = fitgmdist(ptsReplicated, K, 'replicates',5);

% draw the clustering
[y,x]=find(bw);
idx = objPR.cluster([x,y]);
for kk=1:K
    idxK=find(idx==kk);
    plot(x(idxK),y(idxK),'.','color',cmap(kk,:))
end

%draw the gaussian isocontours for the triangular elliptical approximation
%(see online methods). these are exact for circles, approximate for
%ellipses
for kk=1:K
    ptsEdge=genEllipse(objPR.mu(kk,:),objPR.Sigma(:,:,kk),sqrt(20/3),0);
    plot(ptsEdge(:,1),ptsEdge(:,2),'-','color',cmap(kk,:),'linewidth',2);
end