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

Question

我正在尝试创建一个包含 0 个值的矩阵，其中 1 个值填充一个椭圆形。我的椭圆是使用 minVolEllipse.m ( Link 1 ) 生成的，它返回“中心形式”和椭圆中心的椭圆方程矩阵。然后，我使用 Ellipse_plot.m(来自上述链接)中的一段代码将向量参数化为长轴/短轴，生成参数方程，并生成变换坐标矩阵。您可以查看他们的代码以了解这是如何完成的。结果是一个矩阵，该矩阵具有沿椭圆的点的索引位置。它不会包含沿着椭圆轮廓的每个值，除非我将网格点数 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循环的情况下实现这个目标(总图像大小仍然是[width height])？这更快的原因是因为我不必遍历整个图像来确定我的点是否在椭圆内。相反，我可以简单地遍历一个方形区域，该区域的长度为中心的长度 +/- 最大主轴。

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;

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