demoPR_3D.m

%
% demoPR_3D - 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,pts_rc]=GetRandomEllipseImage3D(K); 

hold off
subplot(1,3,1)
imshow(max(bw,[],3))
subplot(1,3,2)
imshow(squeeze(max(bw,[],2)))
subplot(1,3,3)
imshow(squeeze(max(bw,[],1)))
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

ptsReplicated_xy = PixelReplicate(bw);

% 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_xy, K, 'replicates',5);

% draw the clustering
ind=find(bw);
coord_xy = SwapXY_RC(IndToCoord(size(bw),ind));
idx = objPR.cluster(coord_xy);
for kk=1:K
    idxK=find(idx==kk);
    
    subplot(1,3,1)
    hold on
    plot(coord_xy(idxK,1),coord_xy(idxK,2),'.','color',cmap(kk,:))
    subplot(1,3,2)
    hold on
    plot(coord_xy(idxK,3),coord_xy(idxK,2),'.','color',cmap(kk,:))
    subplot(1,3,3)
    hold on
    plot(coord_xy(idxK,3),coord_xy(idxK,1),'.','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
    idxK=find(idx==kk);
    
    subplot(1,3,1)
    hold on
    curCoord = coord_xy(idxK,[1,2]);
    objPR = fitgmdist(curCoord,1,'replicates',5);
    ptsEdge=genEllipse(objPR.mu,objPR.Sigma,sqrt(20/3),0);
    plot(ptsEdge(:,1),ptsEdge(:,2),'-','color',cmap(kk,:),'linewidth',2);
    title('Z projection');
    
    subplot(1,3,2)
    hold on
    curCoord = coord_xy(idxK,[3,2]);
    objPR = fitgmdist(curCoord,1,'replicates',5);
    ptsEdge=genEllipse(objPR.mu,objPR.Sigma,sqrt(20/3),0);
    plot(ptsEdge(:,1),ptsEdge(:,2),'-','color',cmap(kk,:),'linewidth',2);
    title('X projection');
    
    subplot(1,3,3)
    hold on
    curCoord = coord_xy(idxK,[3,1]);
    objPR = fitgmdist(curCoord,1,'replicates',5);
    ptsEdge=genEllipse(objPR.mu,objPR.Sigma,sqrt(20/3),0);
    plot(ptsEdge(:,1),ptsEdge(:,2),'-','color',cmap(kk,:),'linewidth',2);
    title('Y projection');
end