function [ellipse_cells, boundaries] = ellipse_fit_on_cells(cells)
%ELLIPSE_FIT_ON_CELLS Ellipse fitting for cells or other binary objects
%   Fits ellipses on segmented cells modelled as connected components. 
%   Each cell (or any other object) is replaced by an equivalent ellipse 
%   computed considering the cells's morphological features. 
%
% Author: Salim Arslan
% Last updated: 28/06/2014
%
% Input
%  cells: Binary or segmented image of cells
%
% Output
%  ellipse_cells: Labeled ellipses fitted on each cell in the input image
%  boundaries: Label map containing only the perimeter pixels of ellipses
%

    cells = bwlabeln(cells);
    ellipse_cells = zeros(size(cells));
    boundaries = zeros(size(cells));

    for i = 1 : max(max(cells))
        cell = cells == i;
        [ellipse_filled, ellipse_boundary] = ellipse_fit(cell);
        ellipse_cells(ellipse_filled == 1) = i;
        boundaries(ellipse_boundary == 1) = i;
    end

end

function [ellipse_filled, ellipse_boundary] = ellipse_fit(cell)

    phi = linspace(0,2*pi,100);
    cosphi = cos(phi);
    sinphi = sin(phi);

    s = regionprops(cell, 'Centroid', ...
                          'MajorAxisLength', ...
                          'MinorAxisLength', ...
                          'Orientation');
    xbar = s.Centroid(1);
    ybar = s.Centroid(2);

    a = s.MajorAxisLength/2;
    b = s.MinorAxisLength/2;

    theta = pi*s.Orientation/180;
    R = [ cos(theta)   sin(theta)
         -sin(theta)   cos(theta)];

    xy = [a*cosphi; b*sinphi];
    xy = R*xy;

    x = xy(1,:) + xbar;
    y = xy(2,:) + ybar;

    ellipse_filled = roipoly(zeros(size(cell)), x, y);
    ellipse_boundary = bwperim(ellipse_filled);

end