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

Question

我正在尝试创建一个由 0 个值组成的矩阵，其中 1 个值填充一个椭圆形状。我的椭圆是使用 minVolEllipse.m (Link 1) 生成的，它以“中心形式”和椭圆的中心返回椭圆方程的矩阵。然后我使用 Ellipse_plot.m 中的 sn-p 代码（来自上述链接）将向量参数化为长轴/短轴，生成参数方程，并生成变换坐标矩阵。您可以查看他们的代码以了解这是如何完成的。结果是一个矩阵，该矩阵具有沿椭圆的点的索引位置。除非我将网格点的数量 N 设置为一个高得离谱的值，否则它不会包含沿椭圆轮廓的所有值。

当我使用 MATLAB plot 或 patch 命令时，我看到的正是我正在寻找的结果。但是，我希望将其表示为 0 值和 1 的矩阵，其中补丁“填充”空白。很明显 MATLAB 有这个功能，但我还没有找到执行它的代码。我正在寻找的类似于图像处理工具箱的 bwfill 工作方式（Link 2）。 bwfill 对我不起作用，因为我的椭圆不连续，因此该函数返回一个完全填充有 1 个值的矩阵。

希望我已经很好地概述了问题，如果没有，请发表评论，我可以编辑帖子以澄清。

编辑：

我设计了一个策略，使用来自 Ellipse_plot.m 的二维 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-loop 解决方案：

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

Answer 1

展开矩阵乘法，这是一个椭圆范数，给出了一个相当简单的向量化表达式：

[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;

这就是我最初评论的意图。