基于非连续轮廓创建包含填充椭圆的矩阵

我正在尝试创建一个0值的矩阵，其中1个值填充椭圆形状 . 我的椭圆是使用minVolEllipse.m（Link 1）生成的，它返回'center form'中椭圆方程的矩阵和椭圆的中心 . 然后，我使用Ellipse_plot.m的代码片段（来自前面提到的链接）将矢量参数化为主/短轴，生成参数方程，并生成变换坐标矩阵 . 您可以看到他们的代码，看看这是如何完成的 . 结果是一个矩阵，它具有沿椭圆的点的索引位置 . 除非我将网格点数N设置为一个非常高的值，否则它不会包含沿椭圆轮廓的每个值 .

当我使用MATLAB绘图或补丁命令时，我看到确切的结果我填写了'空白' . 显然MATLAB具有此功能，但我还没有找到执行它的代码 . 我正在寻找的是类似于图像处理工具箱的bwfill工作方式（Link 2） . bwfill对我不起作用，因为我的椭圆不是连续的，所以函数返回一个完全填充1个值的矩阵 .

希望我已经很好地概述了这个问题，如果没有请注释，我可以编辑帖子以澄清 .

编辑：

我已经设计了一个策略，使用Ellipse_plot.m中的2-D X向量作为EllipseDirectFit.m（Link 3）的输入 . 该函数返回椭圆函数ax ^ 2 bxy cy ^ 2 dx dy f = 0的系数 . 使用这些系数，我计算x轴和椭圆长轴之间的角度 . 该角度与中心轴和长轴/短轴一起传递到ellipseMatrix.m（Link 4），返回填充矩阵 . 不幸的是，矩阵似乎与我想要的不一致 . 这是我的代码的一部分：

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
%  nima@seas.upenn.edu
%  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);

编辑2：

作为对@BenVoigt建议的回应，我已经为这个问题编写了一个for循环解决方案：

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

虽然这在技术上是解决问题的方法，但我对非迭代解决方案感兴趣 . 我在很多大图像上运行这个脚本，需要它非常快速和并行 .

编辑3：

感谢@ mathematical.coffee获得此解决方案：

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

但是，我相信还有更好的方法可以做到这一点 . 这是一个for循环实现，我做的比上述尝试运行得更快：

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

有没有办法实现这个目标（总图像大小仍然是[宽度高度]）没有这个for循环？这个更快的原因是因为我不必迭代整个图像来确定我的点是否在椭圆内 . 相反，我可以简单地迭代一个正方形区域，该区域是中心的长度/ - 最大的主轴 .