3

0個の値のマトリックスを作成しようとしています。1個の値が楕円形を塗りつぶしています。私の楕円は、minVolEllipse.m（リンク1）を使用して生成されました。これは、楕円方程式の行列を「中心形式」と楕円の中心で返します。次に、Ellipse_plot.m（前述のリンクから）のコードスニペットを使用して、ベクトルを長軸/短軸にパラメーター化し、パラメトリック方程式を生成し、変換された座標の行列を生成します。あなたはこれがどのように行われるかを見るために彼らのコードを見ることができます。結果は、楕円に沿ったポイントのインデックス位置を持つ行列です。グリッドポイントの数Nを途方もなく高い値に設定しない限り、楕円の輪郭に沿ったすべての値が含まれるわけではありません。

MATLABのプロットまたはパッチコマンドを使用すると、探している結果が正確に表示されます。ただし、これを0の値と1の行列として表現し、パッチが空白を埋める必要があります。MATLABにこの機能があることは明らかですが、それを実行するためのコードをまだ見つけていません。私が探しているのは、画像処理ツールボックスのbwfillがどのように機能するかに似ています（リンク2）。楕円が連続していないため、bwfillは機能しません。そのため、関数は1つの値で完全に満たされた行列を返します。

うまくいけば、私は問題の概要を十分に説明しました。そうでない場合はコメントしてください。投稿を編集して明確にすることができます。

編集：

Ellipse_plot.mからの2-DXベクトルをEllipseDirectFit.mへの入力として使用する戦略を考案しました（リンク3）。この関数は、楕円関数ax ^ 2 + bxy + cy ^ 2 + dx + dy + f=0の係数を返します。これらの係数を使用して、楕円のx軸と主軸の間の角度を計算します。この角度は、中心軸と長軸/短軸とともにellipseMatrix.m（リンク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ループなしでこの目標を達成する方法はありますか（合計画像サイズは[幅の高さ]のままです）？これが高速である理由は、ポイントが楕円内にあるかどうかを判断するために画像全体を反復処理する必要がないためです。代わりに、中心の長さ+/-最大の主軸である正方形の領域を単純に反復することができます。

4

1 に答える 1

2

楕円ノルムである行列の乗算を展開すると、かなり単純なベクトル化された式が得られます。

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

これは私が最初のコメントで意図したことです。

于 2012-01-24T16:16:01.007 に答える