[VTK]数値計算の前処理などで利用できそうなVTKの使用方法を紹介する

# はじめに
この記事では、自分が数値計算関連の前処理などにおいて、VTK(Visualization Toolkit)を利用する際に便利と感じた内容について紹介したいと思います。

なお、VTKはPythonやJavaなどでも利用が可能ですが、本記事ではC++で利用することを前提にして書きたいと考えております。(コードが汚いのはお許しください…)

(本記事を読んで、「これ違うよ！」っていう部分があればご指摘いただけると助かります！)

# VTKとは
VTK(Visualization Toolkit)とは、可視化や3DCGの技術をパイプラインの組み合わせで実現できるOSSです。
OpenGLなどよりも比較的単純に利用でき、またグラフィックスの知識が少なくとも形状処理のプログラムを作ることができるので、「数値計算用のメッシュをちょっといじりたいなー」という時などに使えます。

また、「本気で前処理をしたい！」という方は、VTK以外にも色々な形状処理ツールがありますので、そちらも調べてみるといいでしょう。


# VTKの導入
VTKの導入については[こちらの記事](https://qiita.com/MusicScience37/items/0d197b4fed062b880db1)を参照すると良いと思われます。なお、この記事の説明はバージョン8.1.1に対応するものであり、最新のもの(8.1.2)とは異なりますが特に問題はありません。パス名などに注意を払いながら、記事通りに作業を進めていけば問題なく導入ができます。

### 導入時の注意点
- 開発環境によっては途中で「～.dllが足りません！」などとしつこく言われるかもしれませんが、自分の場合は無視しても大丈夫でした。
- プロジェクトの「プロパティ」にて必要事項を入力する前に、ビルドのモードが適切なものになっているか(Releaseや64bitモードになっているかそうでないか)などに気を付けてください。
- コントロールパネルの「システム詳細設定」から環境変数に手を加えたくない人は、Visual Studio 2017のプロジェクトの [プロパティ] -> [デバッグ] -> [環境] を選択し、以下のように「PATH=」の後にDLLまでの絶対/相対パスを入力してください。

![DLLのPath設定.png](https://qiita-image-store.s3.amazonaws.com/0/304390/dc3cbc98-cc5f-eed1-974b-65f4e656b2e4.png)

# VTKでの形状処理例
## フィルタの利用
VTKでの形状処理はフィルタと呼ばれるものを利用します。以下に、数値計算の前処理などで使えそうな例を示します。

### 具体例
- ソリッドメッシュをサーフェスメッシュに変える([vtkGeometryFilter](https://www.vtk.org/Wiki/VTK/Examples/Cxx/PolyData/GeometryFilter))
- 繋がっていない形状を調べる([vtkConnectivityFilter](https://www.vtk.org/Wiki/VTK/Examples/Cxx/Filtering/ConnectivityFilter))
- 複数の独立した形状を一つの形状にまとめる([vtkAppendFilter](https://www.vtk.org/Wiki/VTK/Examples/Cxx/Filtering/AppendFilter))
- 重複節点を除去する([vtkCleanPolyData](https://www.vtk.org/Wiki/VTK/Examples/Cxx/PolyData/CleanPolyData))
- (節点/セルごとに)スカラー値を持つ形状から、指定した値を持つ領域のみ抽出する([vtkThreshold](https://www.vtk.org/Wiki/VTK/Examples/Cxx/PolyData/ThresholdCells))
- 凹凸の多い形状を滑らかにする([vtkSmoothPolyDataFilter](https://www.vtk.org/Wiki/VTK/Examples/Cxx/PolyData/SmoothPolyDataFilter))

特に、vtkConnectivityFilterとvtkThresholdは組み合わせて使うことが多いです。ちなみにこれらはParaViewの [Filter] からも利用ができますので、簡単なテストをしたい場合は自分でプログラムを作るよりParaViewを使うことをお勧めします。

### 実装例(vtkGeometryFilterを使ってみる)
※ソースコード中の入力ファイル「testUnstructuredGrid.vtk」は、ParaViewの [Sources] -> [Unstructured Cell Types]を選択し、[Properties] 内の [Apply] を選択することで生成した、ソリッド要素のボクセルメッシュ(以下、ソリッドメッシュとします)をvtkの形式で保存した物を利用しています。


``` SolidToSurface.cpp
#include <iostream>

#include <vtkDataSetReader.h>
#include <vtkGeometryFilter.h>
#include <vtkDataSetWriter.h>

int main(int, char *argv[])
{
  // vtk legacy format の形状(ソリッド)を読み込む
  vtkSmartPointer<vtkDataSetReader> reader =
    vtkSmartPointer<vtkDataSetReader>::New();
  reader->SetFileName("C:\\testVTK\\testUnstructuredGrid.vtk");
  reader->Update();

  // ソリッド → サーフェス
  vtkSmartPointer<vtkGeometryFilter> geometry =
    vtkSmartPointer<vtkGeometryFilter>::New();
  geometry->SetInputConnection(reader->GetOutputPort());
  geometry->Update();

  // vtk legacy format の形状(サーフェス)を書き込む
  vtkSmartPointer<vtkDataSetWriter> writer =
    vtkSmartPointer<vtkDataSetWriter>::New();
  writer->SetInputConnection(geometry->GetOutputPort());
  writer->SetFileName("C:\\testVTK\\testSurface.vtk");
  writer->Update();

  return 0;
}
```
上記のコードを実行すると、指定パスにソリッドメッシュの表面部分だけを抽出したサーフェスメッシュのvtkファイルが生成されます。
このファイルをParaViewにドラッグ&ドロップで放り込むと、フィルタ後の形状を詳しく見ることができます。
![paraviewで作った形状を見る.png](https://qiita-image-store.s3.amazonaws.com/0/304390/0b723505-06c7-66c8-6124-de841552eb66.png)

なお、この例ではvtkGeometryFilter適用前後でモデル中の要素のタイプ(CellType)がHexa(三次元)からQuad(二次元)に変化しています。

## フィルタ前後での節点/セルが持つ情報を見てみる

- フィルタ前のvtkの形状データに対して予め節点やセルにスカラー値などを設定しておくと、その値が形状処理後にも残ってくれます。
- 例えば、ソリッドメッシュからサーフェスメッシュに変換した形状の節点/セルに、フィルタ前の形状が持っていた情報を保持するといったことができます。


### 実装例(フィルタ前後でのセルが持つ情報の変化を見る)
``` SolidIdsToSurfaceIds.cpp
#include <iostream>

#include <vtkDataSetReader.h>
#include <vtkUnstructuredGrid.h>
#include <vtkIdTypeArray.h>
#include <vtkCellData.h>
#include <vtkGeometryFilter.h>
#include <vtkDataSetWriter.h>

void SetDSCellIds(vtkDataSet *inDS, const char *inArrayName);

int main(int, char *argv[])
{
  // vtk legacy format の形状(ソリッド)を読み込む
  vtkSmartPointer<vtkDataSetReader> reader =
    vtkSmartPointer<vtkDataSetReader>::New();
  reader->SetFileName("C:\\testVTK\\testUnstructuredGrid.vtk");
  reader->Update();

  vtkSmartPointer<vtkUnstructuredGrid> unsGrid =
    vtkSmartPointer<vtkUnstructuredGrid>::New();
  unsGrid->DeepCopy(reader->GetOutput());   // readerの持つ形状をコピー
  //dataSet->NewInstance();
  // セルの通し番号を登録
  SetDSCellIds(unsGrid, "OriginalCellIds");

  // ソリッド → サーフェス
  vtkSmartPointer<vtkGeometryFilter> geometry =
    vtkSmartPointer<vtkGeometryFilter>::New();
  geometry->SetInputData(unsGrid);
  geometry->Update();

  // vtk legacy format の形状(サーフェス)を書き込む
  vtkSmartPointer<vtkDataSetWriter> writer =
    vtkSmartPointer<vtkDataSetWriter>::New();
  writer->SetInputConnection(geometry->GetOutputPort());
  writer->SetFileName("C:\\testVTK\\testSurface.vtk");
  writer->Update();

  return 0;
}

void SetDSCellIds(vtkDataSet *inDS, const char *inArrayName)  // 汎用性を考慮して、引数は抽象クラスで受け取るようにする
{
  // 空のvtkDataArrayに形状データのセル通し番号を登録
  vtkSmartPointer<vtkIdTypeArray> idArray =
    vtkSmartPointer<vtkIdTypeArray>::New();
  idArray->SetName(inArrayName);
  idArray->SetNumberOfComponents(1); // 1次元の配列を定義

  for (int cId = 0; cId < inDS->GetNumberOfCells(); cId++)
  {
    idArray->InsertNextTuple1(cId);
  }

  // 形状データのCellDataにDataArrayを登録する
  inDS->GetCellData()->AddArray(idArray);
}
```
上記のコードから生成したサーフェスメッシュのVTKファイル(vtk legacy format)の中身を見てみると、「CELL_DATA」内の「OriginalCellIds」以下にソリッドメッシュに登録させたセルIDが残っています。

![フィルタ前後での通し番号の変化.png](https://qiita-image-store.s3.amazonaws.com/0/304390/29fc91d0-a9eb-80ab-a0db-a66ecf11d7a2.png)

また、このセルIDは次のように、サーフェスメッシュのセル数分だけ0から順に参照することで取得できます。

``` sampleMethod.cpp
// サーフェスメッシュに登録されたソリッドメッシュのセルIDを取得
void ConvertSolidIdsToSurfaceIds(vtkDataSet *inSurfaceDS, vtkIdTypeArray *inSolidIdsOnSurface) // 第2引数は空の配列を渡す
{
  //サーフェスメッシュの持つセルの分だけ処理
  vtkIdTypeArray *solidIdsOnSurface 
    = vtkIdTypeArray::SafeDownCast(inSurfaceDS->GetCellData()->GetAbstractArray("OriginalCellIds"));

  for (int cId = 0; cId < inSurfaceDS->GetNumberOfCells(); cId++)
  {
    inSolidIdsOnSurface->InsertNextTuple1(solidIdsOnSurface->GetTuple1(cId));
  }
}
```

以上のように形状データに予め配列をぶら下げることで、フィルタ前後で元の形状のセルIDがどのように変化したかを調べることができます。
これの何が嬉しいかと言いますと、例えばオリジナルのモデルの表面に何らかの境界条件を付したい時に、vtkGeometryFilter適用後の形状のセルIDを拾ってくるだけで設定したいソリッドモデルのセルを参照出来るので非常に便利です。

もちろん、節点ID対しても同様のことができますし、IDではなく特定のスカラー値をぶら下げることも可能です。

# VTKファイルのフォーマットと形状データ構造の種類
- VTKを使って生成した形状データは、Legacy FormatとXML Formatの二つで管理されます。
- 本記事ではフォーマットについての説明は省きますが、特にLegacy Formatの内部構造は多少理解しておくと色々便利なので、[こちら](https://qiita.com/shohirose/items/27de6e3f1005cbb67abe)や[こちら](https://qiita.com/kaityo256/items/661833e9e2bfbac31d4b)、[こちら](https://qiita.com/shohirose/items/dd2bdc4201644455e3a8#vtk%E5%BD%A2%E5%BC%8Fvtk-file-format)などを読んでおくことをお勧めします。
- 特に、UnstructuredGridやPolyDataなどの形状に関するデータ構造は頻繁に利用するので慣れておくと無難です。

# VTKに興味が沸いてきたけど、メソッドの使用方法などが分からんという人
### フィルタの使用手順
基本的にVTKのフィルタは以下のような手順で使用することが殆どです。

    1. filter->SetInputData(preFilter->GetOutput());　// または filter->SetInputConnection(preFilter->GetOutputPort);
    2. filter->SetHoge(piyo);    // 何らかの諸元を設定
    3. filter->Update();   // 設定を反映

### 公式資料
(英語ですが)VTKには公式の[Reference](https://www.vtk.org/doc/nightly/html/annotated.html)や豊富な[Examples](https://lorensen.github.io/VTKExamples/site/Cxx/)が用意されており、WEB検索をすると色々な実装例が出てきます。

実装例を少しいじれば自分のやりたいことを実現できることもあるので、はじめの内はExamplesを参考にしながら実装すると良いと思われます。