vtkでトーラスの陰関数を計算する


方法1 vtkTorus

(自分で陰関数を実装するのでなければ)vtkの陰関数クラスの命名方法から推定するとvtkTorusというクラスがありそうなものなので、googleで検索すると Sebastien Jourdain の作成したvtkTorus.cxxのコミットに関するmlでの議論がヒットする。これは結局バグトラッカーへ送られるが2016年に新システムへmoveされたあとcloseされている。しかしながら旧バグトラッカーの添付ファイル
https://www.vtk.org/Bug/file/7287/vtkTorus-contribution.zip
は生存しているのでこれをビルドすれば他の陰関数クラス(ex. vtkSphere)等と同様に使用可能。

方法2 vtkSuperquadric

Jourdain氏はその後kitwareへ就職しているにも関わらず件のコミットが放置されているというのは不自然、というかこれはおそらく既に他にTorus陰関数を得る方法が現存しているため、と考えた方が自然なので、vtkの中からそれらしいクラスを探してみるとvtkSuperquadricというものが見つかる。

kitware webでのexampleではアステロイド状の物体の表示例となっており、トーラスを描画するためにはパラメータを変更する必要があるが、vtkSuperquadricクラスのパラメータが何を意味するかわかりづらい。そこでvtkSuperquadric.hのコメントにある参考文献 "Rigid physically based superquadrics", A. H. Barr,in "Graphics Gems III", David Kirk, ed., Academic Press, 1992.(以下"文献"と略)にあたってみる。つまり、

  • 文献での定式化パラメータとvtkSuperquadricのパラメータの関係
  • いわゆるトーラスのパラメータ(大半径、小半径等)とそれらの関係

を整理してvtkSuperquadricでの(楕円)トーラス描画を実行する。

文献での扱い

文献では超楕円体(superquadric)は6個のパラメータで定義されている。即ち

  1. $n$: roundness/squareness paremeter in north-south direction
  2. $e$: roundness/squareness paremeter in east-west direction
  3. $a_1$: length
  4. $a_2$: width
  5. $a_3$: depth
  6. $\alpha$: hole diameter (環状体の場合。) > 1

文献 p.148で与えられている環状体(toroid)の陰関数では、$n$はポロイダル方向の、$e$はトロイダル方向の丸み(roundness)であるので、軸対称トーラス(いわゆるドーナッツ型)の場合、n=1, e=1を代入し、さらに軸対称なので$a_1 = a_2$とすると

f(x, y, z) = 
\left| \sqrt{\left(\frac{x}{a_1}\right)^2 
+ \left(\frac{y}{a_1}\right)^2} -\alpha \right|^2
+ \left(\frac{z}{a_3}\right)^2 -1

となる。

トーラスパラメータ

大半径$R$, トロイダル軸方向小半径$r_1$, トロイダル径方向小半径$r_2$とすると、軸対称楕円トーラスの陰関数は

f(x, y, z) =
\left(
\sqrt{\left(\frac{x}{r_2}\right)^2 + \left(\frac{y}{r_2}\right)^2} 
      -\frac{R}{r_2}
 \right)^2
+\left(\frac{z}{r_1}\right)^2 -1

であるので文献の陰関数のパラメータとの関係は

  • $\alpha =R/r_2$
  • $a_1 = a_2 = r_2$
  • $a_3 = r_1$

となる。

vtkでの扱い

vtkSuperquadricでは制御パラメータとして
1. Toroidal(bool = false)
2. Thickness(double = 0.3333)
3. PhiRoundness(double = 1)
4. ThetaRoundness(double = 1)
5. Center(double[3]={0, 0, 0})
6. Scale(double[3]={1, 1, 1})
7. Size(double =0.5)
の7個(括弧内はデフォルト値)が用いられている。ソースを見ると、これらと$n, e, a_{123},\alpha$のトーラスの場合の関係は

  • Thickness = $1/\alpha = r_2/R$
  • PhiRoundness = $n$
  • ThetaRoundness = $e$
  • Scale[i]*Size = $(\alpha+1)$ {${a_1, a_3, a_2}$} = $(\alpha+1)${$r_2, r_1, r_2$}※

※vtkSuperquadricはトロイダル軸方向をy方向にしている
となっているので、例えば$R=8.5, r_1=5, r_2=1.5$なる楕円トーラスは

auto iFunc = vtkSmartPointer<vtkImplicitFunction>::New();
// PhiRoundness, ThetaRoundnessはデフォルトで1なので設定不要
iFunc->SetToroidal(true); // 環状体(troid)の場合on
iFunc->SetSize(1); // 実サイズはScaleで設定するのでここでは1で固定する

const double R=8.5, r1=5, r2=1.5;
const double alpha = R/r2;
iFunc->SetThickness(1.0/alpha); // -> hole半径は10
iFunc->SetScale((alpha+1)*r2, (alpha+1)*r1, (alpha+1)*r2);

とすればよい。vtk exampleのsample functionの陰関数部分を上記の例で書き換えるとこのように表示される。