マーチングキューブ法 (Marching cubes algorithm)を用いた簡易的なメタボール
はじめに
この記事はコンピュータグラフィックスで扱われるアルゴリズム手法の一つマーチングキューブ法を実装した記事である．マーチングキューブ法とは3次元ボクセルデータをポリゴンに変換するアルゴリズムである．三次元表示や実体模型化と相性が良くCTやMRIなどに利用されている．
日本語で掲載されている記事はこれが理解しやすい．
改めて，今回はマーチングキューブ法を用いて、3次元空間におけるメタボールを表現した。マーチングキューブ法のコードはweb上に公開されているもの (Polygonising a scalar field, Paul Bourke, May 1994)を引用し、それをprocessing仕様にコードに変更を加えて実装した。
プログラムの仕様
まず、3次元空間上の(x,y,z)=(0,0,0)の点を中心に一辺の長さが50の立方体内を描画の範囲に設定した。そして、その立方体内における、一辺の長さが1の立方体を、一つのボクセルとして考え、一辺の長さが50の立方体内に計50×50×50=1250000のボクセルが入っていると考える。のちにマーチングキューブ法を用いて描画を行うために、void setupにて、立方体内にあるすべてのボクセルの8つの頂点すべてに値を持たせている。各頂点の値(=gridvalue)の持たせ方としては、点(10,10,10)からの距離aと点(-10,-10,-10)からの距離bを用いて以下のような関数により設定した。
gridvalue = r(\frac{1}{a}+\frac{1}{b})
rの値は任意の値で、rの値を変えることでメタボールの加減を変更することができる(今回は180に設定)。これにより、閾値(=isolevel)=20より低いボクセルの点をもとに、マーチングキューブ法で描画を行い、実行時のメタボールを表現した。
まず、void setupにてボクセルの８つの頂点それぞれに値(=gridvalue)を持たせた後、void draw内にて、図1のように各ボクセルに対し関数marchingcubes（）を呼び出し、マーチングキューブ法の処理を行っている。
for(int x=m; x<n; x++){
for(int y=m; y<n; y++){
for(int z=m; z<n; z++){
float d=dist(x,y,z,a,a,a)+dist(x,y,z,b,b,b);
fill(r/d, 255, 255);
marchingcubes(isolevel, index, edgeTable, triTable, gridvalue);
index+=1;
}
}
}
マーチングキューブ法の処理を行う関数を呼び出している部分
変数indexがボクセルの番号である。関数marchingcubes()の中では、まず初めに、ボクセルの8つの頂点それぞれが持つ値(=gridvalue)がisolevelより低いか、高いかを判定することにより、ボクセルを描画する方法全256種類の中から一つを定めている(図２)。
cubeindex = 0;
if (gridvalue[index][0] < isolevel) cubeindex += 1;
if (gridvalue[index][1] < isolevel) cubeindex += 2;
if (gridvalue[index][2] < isolevel) cubeindex += 4;
if (gridvalue[index][3] < isolevel) cubeindex += 8;
if (gridvalue[index][4] < isolevel) cubeindex += 16;
if (gridvalue[index][5] < isolevel) cubeindex += 32;
if (gridvalue[index][6] < isolevel) cubeindex += 64;
if (gridvalue[index][7] < isolevel) cubeindex += 128;
全２５６種類の描画方法から一つの描画方法を定める処理を行う部分
この変数cubeindexが持つ値が、描画方法の「番号」である。これ以降ではその描画方法に合わせて、描画に用いる点(詳しくは、三角形を描画する際の三つの点を、算出するための点)を配列vertlistに格納し、その後配列verttlistに格納されている点をもとに三角形を0~複数個描画することでポリゴンを作成する処理を行っている。
色相を付ける
メタボールの色相に変化を加えて実装した。点(10,10,10)からの距離aと点(-10,-10,-10)からの距離bを加算したものに重みrをつけ、以下のようにそれぞれの点からの距離によって色相Hが変化する処理を施している。
$H=r(a+b)$
fill(H,255,255)
実装
今回はProcessingでの実装となるため，コードの変更が必要である．
以下のように論理演算の処理部分を変える必要がある．
cubeindex = 0;
if (grid.val[0] < isolevel) cubeindex |= 1;
if (grid.val[1] < isolevel) cubeindex |= 2;
if (grid.val[2] < isolevel) cubeindex |= 4;
if (grid.val[3] < isolevel) cubeindex |= 8;
if (grid.val[4] < isolevel) cubeindex |= 16;
if (grid.val[5] < isolevel) cubeindex |= 32;
if (grid.val[6] < isolevel) cubeindex |= 64;
if (grid.val[7] < isolevel) cubeindex |= 128;
/* Cube is entirely in/out of the surface */
if (edgeTable[cubeindex] == 0)
return(0);
/* Find the vertices where the surface intersects the cube */
if (edgeTable[cubeindex] & 1)
vertlist[0] =
VertexInterp(isolevel,grid.p[0],grid.p[1],grid.val[0],grid.val[1]);
if (edgeTable[cubeindex] & 2)
vertlist[1] =
VertexInterp(isolevel,grid.p[1],grid.p[2],grid.val[1],grid.val[2]);
if (edgeTable[cubeindex] & 4)
vertlist[2] =
VertexInterp(isolevel,grid.p[2],grid.p[3],grid.val[2],grid.val[3]);
if (edgeTable[cubeindex] & 8)
vertlist[3] =
VertexInterp(isolevel,grid.p[3],grid.p[0],grid.val[3],grid.val[0]);
if (edgeTable[cubeindex] & 16)
vertlist[4] =
VertexInterp(isolevel,grid.p[4],grid.p[5],grid.val[4],grid.val[5]);
if (edgeTable[cubeindex] & 32)
vertlist[5] =
VertexInterp(isolevel,grid.p[5],grid.p[6],grid.val[5],grid.val[6]);
if (edgeTable[cubeindex] & 64)
vertlist[6] =
VertexInterp(isolevel,grid.p[6],grid.p[7],grid.val[6],grid.val[7]);
if (edgeTable[cubeindex] & 128)
vertlist[7] =
VertexInterp(isolevel,grid.p[7],grid.p[4],grid.val[7],grid.val[4]);
if (edgeTable[cubeindex] & 256)
vertlist[8] =
VertexInterp(isolevel,grid.p[0],grid.p[4],grid.val[0],grid.val[4]);
if (edgeTable[cubeindex] & 512)
vertlist[9] =
VertexInterp(isolevel,grid.p[1],grid.p[5],grid.val[1],grid.val[5]);
if (edgeTable[cubeindex] & 1024)
vertlist[10] =
VertexInterp(isolevel,grid.p[2],grid.p[6],grid.val[2],grid.val[6]);
if (edgeTable[cubeindex] & 2048)
vertlist[11] =
VertexInterp(isolevel,grid.p[3],grid.p[7],grid.val[3],grid.val[7]);
以下のように論理演算を四則演算へと変えて対応した．
i=0;
while(edgeTable[cubeindex][i] != -1){
if(edgeTable[cubeindex][i]<7){
n=edgeTable[cubeindex][i];
vertlist[n] =
VertexInterp(isolevel, xyz[index][n], xyz[index][n+1], gridvalue[index][n], gridvalue[index][n+1]);
}
else if(edgeTable[cubeindex][i]==7){
vertlist[7] =
VertexInterp(isolevel,xyz[index][7],xyz[index][4],gridvalue[index][7],gridvalue[index][4]);
}
else if(edgeTable[cubeindex][i]==8){
vertlist[8] =
VertexInterp(isolevel,xyz[index][0],xyz[index][4],gridvalue[index][0],gridvalue[index][4]);
}
else if(edgeTable[cubeindex][i]==9){
vertlist[9] =
VertexInterp(isolevel,xyz[index][1],xyz[index][5],gridvalue[index][1],gridvalue[index][5]);
}
else if(edgeTable[cubeindex][i]==10){
vertlist[10] =
VertexInterp(isolevel,xyz[index][2],xyz[index][6],gridvalue[index][2],gridvalue[index][6]);
}
else if(edgeTable[cubeindex][i]==11){
vertlist[11] =
VertexInterp(isolevel,xyz[index][3],xyz[index][7],gridvalue[index][3],gridvalue[index][7]);
}
i+=1;
}
自分が実装したコードはここにあるのでぜひ試してほしい．
まとめ
今回はコンピュータグラフィックスの世界で活用されているアルゴリズム，マーチングキューブ法を実装した。マーチングキューブ法の理解にはweb上に公開されている(ソースコード)[http://paulbourke.net/geometry/polygonise/]と日本語でマーチングキューブ法を説明しているweb上の資料から学んだ。また、メタボールについても勉強した。メタボールに関してはThe Coding Trainがyoutubeに公開している動画(Coding Challenge #28: Metaballs)を参考にして学んだ。