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

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

##方法2 vtkSuperquadric
Jourdain氏はその後kitwareへ[就職している](https://www.kitware.com/sebastien-jourdain/)にも関わらず件のコミットが放置されているというのは不自然、というかこれはおそらく既に他にTorus陰関数を得る方法が現存しているため、と考えた方が自然なので、vtkの中からそれらしいクラスを探してみると[vtkSuperquadric](https://vtk.org/doc/nightly/html/classvtkSuperquadric.html)というものが見つかる。

kitware webでの[example](https://lorensen.github.io/VTKExamples/site/Java/ImplicitFunctions/SampleFunction/)ではアステロイド状の物体の表示例となっており、トーラスを描画するためにはパラメータを変更する必要があるが、vtkSuperquadricクラスのパラメータが何を意味するかわかりづらい。そこでvtkSuperquadric.hのコメントにある参考文献 ["Rigid physically based superquadrics", A. H. Barr,in "Graphics Gems III", David Kirk, ed., Academic Press, 1992.](https://www.oreilly.com/library/view/graphics-gems-iii/9780124096738/xhtml/B9780080507552500385.htm)(以下"文献"と略)にあたってみる。つまり、

- 文献での定式化パラメータと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$とすると

```math
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$とすると、軸対称楕円トーラスの陰関数は

```math
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$なる楕円トーラスは

```code
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](https://lorensen.github.io/VTKExamples/site/Cxx/ImplicitFunctions/SampleFunction/)の陰関数部分を上記の例で書き換えるとこのように表示される。
![elliptictorus.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/488299/f3d0a14c-8684-0ee6-11d2-b223ffa67ef0.png)