マーチングキューブ法 (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（）を呼び出し、マーチングキューブ法の処理を行っている。

marching_cubes.pde

 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種類の中から一つを定めている(図２)。

marching_cubes.pde

   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)$

marching_cubes.pde

 fill(H,255,255)

実装

今回はProcessingでの実装となるため，コードの変更が必要である．
以下のように論理演算の処理部分を変える必要がある．

marching_cubes.pde

   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]);

以下のように論理演算を四則演算へと変えて対応した．

marching_cubes.pde

   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)を参考にして学んだ。