Creating a matrix containing a filled ellipse based on a non-contiguous outline

Question

I'm trying to create a matrix of 0 values, with 1 values filling a ellipse shape. My ellipse was generated using minVolEllipse.m (Link 1) which returns a matrix of the ellipse equation in the 'center form' and the center of the ellipse. I then use a snippet of code from Ellipse_plot.m (from the aforementioned link) to parameterize the vector into major/minor axes, generate a parametric equation, and generate a matrix of transformed coordinates. You can see their code to see how this is done. The result is a matrix that has index locations for points along the ellipse. It does not encompass every value along the outline of the ellipse unless I set the number of grid points, N, to a ridiculously high value.

When I use the MATLAB plot or patch commands I see exactly the result I'm looking for. However, I want this represented as a matrix of 0 values with 1s where patch 'fills in' the blanks. It is apparent that MATLAB has this functionality, but I have yet to find the code to execute it. What I am looking for is similar to how bwfill of the image processing toolbox works (Link 2). bwfill does not work for me because my ellipse is not contiguous, so the function returns a matrix filled completely with 1 values.

Hopefully I have outlined the problem well enough, if not please comment and I can edit the post to clarify.

EDIT:

I have devised a strategy using the 2-D X vector from Ellipse_plot.m as an input to EllipseDirectFit.m (Link 3). This function returns the coefficients for the ellipse function ax^2+bxy+cy^2+dx+dy+f=0. Using these coefficients I calculate the angle between the x-axis and the major axis of the ellipse. This angle, along with the center and major/minor axes are passed into ellipseMatrix.m (Link 4), which returns a filled matrix. Unfortunately, the matrix appears to be out of rotation from what I want. Here is the portion of my code:

N = 20; %Number of grid points in ellipse
ellipsepoints = clusterpoints(clusterpoints(:,1)==i,2:3)';
[A,C]         = minVolEllipse(ellipsepoints,0.001,N);
%%%%%%%%%%%%%%
%
%Adapted from:
%  Ellipse_plot.m
%  Nima Moshtagh
%  [email protected]
%  University of Pennsylvania
%  Feb 1, 2007
%  Updated: Feb 3, 2007
%%%%%%%%%%%%%%
%
%
% "singular value decomposition" to extract the orientation and the
% axes of the ellipsoid
%------------------------------------
[U D V] = svd(A);

%
% get the major and minor axes
%------------------------------------
a = 1/sqrt(D(1,1))
b = 1/sqrt(D(2,2))

%theta values
theta = [0:1/N:2*pi+1/N];

%
% Parametric equation of the ellipse
%----------------------------------------
state(1,:) = a*cos(theta); 
state(2,:) = b*sin(theta);

%
% Coordinate transform 
%----------------------------------------
X = V * state;
X(1,:) = X(1,:) + C(1);
X(2,:) = X(2,:) + C(2);

% Output: Elip_Eq = [a b c d e f]' is the vector of algebraic 
%       parameters of the fitting ellipse:
Elip_Eq = EllipseDirectFit(X')
% http://mathworld.wolfram.com/Ellipse.html gives the equation for finding the angle theta (teta).
% The coefficients from EllipseDirectFit are rescaled to match what is expected in the wolfram link.
Elip_Eq(2) = Elip_Eq(2)/2;
Elip_Eq(4) = Elip_Eq(4)/2;
Elip_Eq(5) = Elip_Eq(5)/2;

if Elip_Eq(2)==0
    if Elip_Eq(1) < Elip_Eq(3)
        teta = 0;
    else
        teta = (1/2)*pi;
    endif
else
    tetap = (1/2)*acot((Elip_Eq(1)-Elip_Eq(3))/(Elip_Eq(2)));
    if Elip_Eq(1) < Elip_Eq(3)
        teta = tetap;
    else
        teta = (pi/2)+tetap;
    endif
endif

blank_mask = zeros([height width]);

if teta < 0
    teta = pi+teta;
endif

%I may need to switch a and b, depending on which is larger (so that the fist is the major axis)
filled_mask1 = ellipseMatrix(C(2),C(1),b,a,teta,blank_mask,1);

EDIT 2:

As a response to the suggestion from @BenVoigt, I have written a for-loop solution to the problem, here:

N = 20; %Number of grid points in ellipse
ellipsepoints = clusterpoints(clusterpoints(:,1)==i,2:3)';
[A,C]         = minVolEllipse(ellipsepoints,0.001,N);

filled_mask = zeros([height width]);
for y=0:1:height
   for x=0:1:width
      point = ([x;y]-C)'*A*([x;y]-C);
      if point < 1
         filled_mask(y,x) = 1;
      endif
   endfor
endfor

Although this is technically a solution to the problem, I am interested in a non-iterative solution. I'm running this script over many large images, and need it to be very fast and parallel.

EDIT 3:

Thanks @mathematical.coffee for this solution:

[X,Y] = meshgrid(0:width,0:height);
fill_mask=arrayfun(@(x,y) ([x;y]-C)'*A*([x;y]-C),X,Y) < 1;

However, I believe there is yet a better way to do this. Here is a for-loop implementation that I did that runs faster than both above attempts:

ellip_mask = zeros([height width]);
[U D V] = svd(A);
a = 1/sqrt(D(1,1));
b = 1/sqrt(D(2,2));
maxab = ceil(max(a,b));

xstart = round(max(C(1)-maxab,1));
xend = round(min(C(1)+maxab,width));
ystart = round(max(C(2)-maxab,1));
yend = round(min(C(2)+maxab,height));

for y = ystart:1:yend
   for x = xstart:1:xend
      point = ([x;y]-C)'*A*([x;y]-C);
      if point < 1
         ellip_mask(y,x) = 1;
      endif
   endfor
endfor

Is there a way to accomplish this goal (the total image size is still [width height]) without this for-loop? The reason this is faster is because I don't have to iterate over the entire image to determine if my point is within the ellipse. Instead, I can simply iterate over a square region that is the length of the center +/- the largest principle axis.

Answer 1

Expanding the matrix multiply, which is an elliptic norm, gives a fairly simply vectorized expression:

[X,Y] = meshgrid(0:width,0:height);
X = X - C(1);
Y = Y - C(2);
fill_mask = X.^2 * A(1,1) + X.*Y * (A(1,2) + A(2,1)) + Y.^2 * A(2,2) < 1;

This is what I intended by my original comment.