Mathematics

高次元の散布データの三角形分割と内挿

科学、 工学、統計、数学の応用分野の多くでは、凸包、VoronoiダイアグラムやDelaunay三角形分割等の構造を必要とします。Qhull [1]を使うことにより、MATLAB関数は、あらゆる次元の解析データセットを幾何学的に解析することを可能にします。

この節では、つぎの幾何学的解析のテクニックについて説明します。

• 凸包
• Delaunay三角形分割
• Voronoi ダイアグラム
• 多次元散布データの内挿

凸包

n-次元空間のデータ集合の凸包は、データ集合を含む最小凸領域として定義されます。

凸包の計算   convhulln関数は、凸包 面を構成するデータ集合内にある点のインデックスを戻します。たとえば、Xは、立方体の8つの頂点からなる8行列の行列と仮定します。Xの凸包は、その時、12の面から構成されます。

d = [-1 1];
[x,y,z] = meshgrid(d,d,d);
X = [x(:),y(:),z(:)];       % 立方体の8つの頂点
C = convhulln(X)

C =
     3     1     5
     1     2     5
     2     1     3
     7     3     5
     8     7     5
     7     8     3
     6     8     5
     2     6     5
     6     2     8
     8     4     3
     4     2     3
     2     4     8

データは3-次元なので、凸包を作る面は、三角形です。Cの12行は、12個の三角形を表わしています。Cの要素は、X内の点のインデックスです。たとえば、最初の行の 3 1 5 は、最初の三角形がX(3,:)、X(1,:)、X(5,:)の頂点をもつことを意味しています。

3次元のパッチの三角形を描いて、この凸包を表示しましょう。

figure, hold on
d = [1 2 3 1];      % C列内のインデックス
for i = 1:size(C,1) % 各三角形を描く
    j= C(i,d);      % パッチを作成するために、i番目のCを取得
    h(i)=patch(X(j,1),X(j,2),X(j,3),i,'FaceAlpha',0.9);
end                 % 'FaceAlpha'は、パッチを半透明にするために使用
hold off
view(3), axis equal, axis off
camorbit(90,-5);    % 角度を変えて表示
title('Convex hull of a cube')

Delaunay 三角形

Delaunay三角形は、データ点が、任意のシンプレックスの球の周にも含まれないようなシンプレックスの集合です。2-次元空間では、シンプレックスは三角形。3-次元空間では、シンプレックスは四面体です。

Delaunay 三角形の計算   delaunayn 関数は、データの集合体のn-次元Delaunay三角形のシンプレックスを構成するデータ集合内の点のインデックスを戻します。

この例では、 凸包 例題と同じX、つまり立方体の8つの頂点に、中心点を追加して使用します。

d = [-1 1];
[x,y,z] = meshgrid(d,d,d);
X = [x(:),y(:),z(:)];   % 立方体の8頂点
X(9,:) = [0 0 0];       % 頂点リストに中心点を追加
T = delaunayn(X)      % Delaunay 三角形を生成

T = 
     9     7     3     5
     1     9     3     5
     1     2     9     5
     4     9     7     3
     4     9     7     8
     4     1     9     3
     4     1     2     9
     6     2     9     5
     6     9     7     5
     6     9     7     8
     6     4     9     8
     6     4     2     9

Tの12行は、立方体を区切る12のシンプレックス、この場合変形四面体、を表現しています。 各行が、1つの四面体を表わし、行の要素がXの点のインデックスです。

.3-次元パッチを使って、四面体を描いて、この三角形を表示しましょう。

figure, hold on
d = [1 1 1 2; 2 2 3 3; 3 4 4 4];  % T内のインデックス
for i = 1:size(T,1)   % 各四面体の表示
   y = T(i,d);        % パッチを作るi番目のTを取得
   x1 = reshape(X(y,1),3,4);
   x2 = reshape(X(y,2),3,4);
   x3 = reshape(X(y,3),3,4);
   h(i)=patch(x1,x2,x3,(1:4)*i,'FaceAlpha',0.9);
end
hold off
view(3), axis equal
axis off
camorbit(65,120)      % 角度を変えて表示
title('Delaunay tessellation of a cube with a center point')

cameramenuを使うと任意の方向に図を回転させることができます。

Voronoi ダイアグラム

n-次元空間のmデータ点とすると、Voronoiダイアグラムはn-次元空間を、各データ点に対して1領域となるようなm個の多面体の領域に分割します。このような領域をVoronoiセルと呼びます。Voronoiセルは、全ての点が、集合内の任意のいかなる点よりもそのデータポイントに近い全ての点を含んでいるという条件を満足します。

Voronoi ダイアグラムの計算   voronoin 関数はつぎの2出力を戻します。

• Vは、n-空間内のm点のm行n列の行列です。Vの各行はVoronoi頂点を表わします。
• Cは、ベクトルのセル配列です。セル配列C内の各ベクトルは、Voronoiセルを表現します。ベクトルは、Voronoiセルの頂点であるV内の点のインデックスを含みます。

各Voronoiセルは、異なった点の数をもちます。Voronoiセルは無限にすることができるため、Vの最初の行は無限領域の点です。その後、Cの任意の無限Voronoiセルは、無限領域の点、つまりVの最初の点を含みます。

この例は、 Delaunay 例 と同じX、つまり立方体の8頂点と中心点を使用します。不規則な立方体にするために、ランダムノイズを加えます。結果としてできるVoronoiダイアグラムは、9つのVoronoiセルを持ちます。

d = [-1 1];
[x,y,z] = meshgrid(d,d,d);
X = [x(:),y(:),z(:)];      % 立方体の8頂点
X(9,:) = [0 0 0];          % 頂点リストに中心点を付加
X = X+0.01*rand(size(X));  % 立方体を不規則にする
[V,C] = voronoin(X);

V =
       Inf       Inf       Inf 
    0.0055    1.5054    0.0004
    0.0037    0.0101   -1.4990
    0.0052    0.0087   -1.4990
    0.0030    1.5054    0.0030
    0.0072    0.0072    1.4971
   -1.7912    0.0000    0.0044
   -1.4886    0.0011    0.0036
   -1.4886    0.0002    0.0045
    0.0101    0.0044    1.4971
    1.5115    0.0074    0.0033
    1.5115    0.0081    0.0040
    0.0104   -1.4846   -0.0007
    0.0026   -1.4846    0.0071

C = 
     [1x8 double] 
     [1x6 double] 
     [1x4 double] 
     [1x6 double] 
     [1x6 double] 
     [1x6 double] 
     [1x6 double] 
     [1x6 double] 
     [1x12 double]

この例では、Vは13行3列の行列で、13行は、13個のVoronoi頂点の座標です。Vの最初の行は、無限点です。Cは9行1列のセル配列です。ここで、配列内の各セルは、9個のVoronoiセルのどれか1つに対応するV内のベクトルのインデックスを含んでいます。たとえば、Voronoiダイアグラムの9番目のセルは、つぎのようになります。

C{9} = 2 3 4 5 6 7 8 9 10 11 12 13

もし、セル配列のセル内の任意のインデックスが1ならば、対応するVoronoiセルは、V内の最初の点、つまり無限の点を含んでいます。これは、Voronoiセルが無限領域であることを意味しています。有限のVoronoiセル、つまり、無限の点を含まない1点を表示するために、 convhulln関数を使ってVoronoiセルの頂点をセルを構成する面に変換することができます。たとえば、C内の9番目のセルで定義されるVoronoiセルを表示するためには、凸包で説明されている手順を使用します。

X = V(C{9},:);    % 9番目のVoronoiセルを表示
K = convhulln(X);
figure
hold on
d = [1 2 3 1];    % K内のインデックス
for i = 1:size(K,1)
   j = K(i,d);
   h(i)=patch(X(j,1),X(j,2),X(j,3),i,'FaceAlpha',0.9);
end
hold off
view(3)
axis equal
title('One cell of a Voronoi diagram')

N-次元データの内挿

多次元データ、特に散布データを内挿するために、 griddatan関数を使用します。griddatanは、データを三角形分割するために、delaunayn関数を使用し、そして三角形に基づいて内挿します。

n-個の散布点の集合を評価する関数を可視化すると仮定します。この例では、Xは点のn行3列の行列で、各行は点の1つに(x,y,z)座標をもちます。ベクトルvは、これらの点のn関数値を含んでいます。この例の関数では、原点v = x.^2 + y.^2 + z.^2 からの距離の平方です。

3-次元空間にランダムな点 n=5000点を作成し、これらの点で関数の値を計算することからはじめます。

n = 5000; 
X = 2*rand(n,3)-1; 
v = sum(X.^2,2);

次のステップで、グリッド上に関数値を計算するために内挿を使用します。グリッドを作成するために、 meshgridを使用し、内挿するために griddatanを使用します。

delta = 0.05;
d = -1:delta:1;
[x0,y0,z0] = meshgrid(d,d,d);
X0 = [x0(:), y0(:), z0(:)];
v0 = griddatan(X,v,X0);
v0 = reshape(v0, size(x0));

つぎに、関数一定値をとるような(x,y,z)値からなるサーフェスを可視化するための等値サーフェス関連関数を使用します。任意の値を選べますが、例では、0.6という値を使用しています。関数は、原点からの距離の平方根のため、一定値のサーフェスは球になります。

p = patch(isosurface(x0,y0,z0,v0,0.6)); 
isonormals(x0,y0,z0,v0,p); 
set(p,'FaceColor','red','EdgeColor','none'); 
view(3); 
camlight; 
lighting phong 
axis equal 
title('Interpolated sphere from scattered data')

注意    deltaを小さくすると、より滑らかな球が作成されますが、コンピュータの計算時間が長くなります。

散布データの三角メッシュと内挿

参考文献