More than 3 years have passed since last update.

C++Advent Calendar 2018Day 14

VTKライブラリ

Last updated at 2019-12-14Posted at 2018-12-13

はじめに

この記事は、C++ Advent Calendar 2018の14日目です。
前日の記事は@tyanmahouさんでした。
明日の記事は@yumetodoさんです。

アメリカにいると15時間の時差のせいで13日に記事を投稿しないければならず、なんとなく損した気分になります笑（実際はそんなことないけど）。

VTKライブラリとは

Visualization Toolkit（VTK）ライブラリは（科学）データを操作・表示するためのC++で書かれたオープンソースライブラリです。アカデミック分野でよく使われているデータ解析・グラフィックス表示ソフトウェアParaViewは、このVTKライブラリを元に作られています。ParaViewはノートパソコンからスパコンまでの多様なコンピュータ上で動き、ペタサイズの巨大なデータセットすら取り扱うことができます。

VTKを使ったグラフィックスの例

シミュレーション結果をParaViewで表示するためには、ParaViewがサポートするファイル形式でデータを出力する必要があります。ParaViewは多様な形式をサポートしていますが、通常はVTK形式（VTK file format）を使います。

VTKライブラリにVTK形式で入出力するためのクラスが用意されているのですが、巨大なライブラリであることと、オブジェクト指向に基づいて継承を多用して実装しているため、C++初心者にとって使いやすいライブラリとは言い難いです。英語を理解できるならば、オンラインレファレンスや使用例を見るよりも、まずは無料でダウンロードできるユーザーガイドや教科書から必要な個所を拾い読みする方がいいでしょう。

本記事では、VTKライブラリの基本について解説します。VTK形式で入出力する処理をC++で書きたいと考えているアカデミック分野の学生の手助けになれば幸いです。

（本当はVTK形式の入出力を簡単にするラッパークラスを作って公開しようと思っていたのですが、時間が足りず断念しました…）

VTKライブラリ

概要

VTKライブラリは、3Dグラフィックスを表示するための機能をオブジェクト指向に基づいて実装したライブラリとして、1993年に作られました。

VTK | About | History

ほとんど全てのVTKライブラリのクラスはvtkObjectを継承して作られ、必ずヒープ領域にメモリが確保されます。vtkObjectは自身を参照するオブジェクトを数えており（参照カウント方式）、後述するVTKライブラリの提供するスマートポインタ（vtkNew/vtkSmartPointer/vtkWeakPointer）と共に使うことで、std::unique_ptr/std::shared_ptr/std::weak_ptrに似た機能を実現しています。

VTKライブラリのクラスの生成は必ずこれらのスマートポインタを通しておこなうため、メンバ関数はアロー演算子を通して呼び出します。25年前に作られたライブラリで、現代ではあまり多用されない設計に基づいているため、使い方に慣れるまで時間がかかると思います。

VTKライブラリのスマートポインタ

vtkライブラリの提供するスマートポインタ

vtkNew
vtkSmartPointer
vtkWeakPointer

については、VTKライブラリの開発元であるKitwareのブログ「A Tour of VTK Pointer Classes」でわかりやすく解説されています。簡単にまとめると、

vtkNewは、オブジェクト生成時にメモリを確保し、スコープを抜けると自動的に破棄します。vtkNewオブジェクト生成後に他のオブジェクトへの参照を代入することはできません。
vtkSmartPointerはvtkObjectへの参照を保持します。メモリ確保は各クラスのstatic関数New()を通して行います。生成した後で参照するオブジェクトを変更することが可能です。
vtkWeakPointerはダングリングを避けるためのスマートポインタです。

基本的に、後から参照を変更したり割り当てたい場合はvtkSmartPointerを使い、そうでない場合はvtkNewを使いましょう。

VTK形式（VTK File Format）

VTK形式は、幾何形状を規定するデータ（Grid Data）および頂点とセルに付随するデータ（Field Data）に分けられます。後述するLegacy形式では、Grid Dataとして

Structured Points
Structured Grid
Rectilinear Grid
Polygonal Grid
Unstructured Grid

の5つが定義されています。

VTK形式はその中でLegacy/XMLに分けられ、それぞれASCII/Binaryを選択できます。つまり下の表の(1)~(4)の組み合わせがあります。

	ASCII	Binary
Legacy	(1)	(2)
XML	(3)	(4)

XML形式は並列入出力をサポートしているため、大規模なデータセットの場合はXML形式を使用しましょう。

Legacy形式の場合、VTKファイルの拡張子はvtkのみですが、XML形式の場合はSerial/Parallelおよびデータ型ごとに拡張子が分けられています。詳しくはVTK file formatのXML File Formatsの節を参照してください。

Grid Dataは、頂点（Point）のリストとセル（Cell）のリストからなります（セルが必要ないグリッド形式もあります）。頂点は三次元座標(x y z)で表され、セルはセルを構成する頂点の番号のリスト(v1 v2 ...)として表されます。

Field Dataは、頂点またはセルに紐付けられたものと、幾何形状とは全く関連のないものに分けられます。スカラー・ベクトル・対称テンソル・テンソルがサポートされています。各成分の順番は、ベクトルは(x y z)、対称テンソルは(xx yy zz xy yz xz)、テンソルは(xx xy xz yx yy yz zx zy zz)です。

VTK形式の入出力

グリッドデータの入手力

例えば以下のようなオリジナルのメッシュクラスのデータをUnstructured GridとしてVTK形式に入出力したいとします。

struct point {
  double x;
  double y;
  double z;
};

using tetra_cell = std::array<std::ptrdiff_t, 4>; 

struct tetra_mesh {
  std::vector<point> points;
  std::vector<tetra_cell> cells;

  auto n_points() const noexcept { return points.size(); }
  auto n_cells() const noexcept { return cells.size(); }
};

VTKファイルからこのメッシュクラスにデータをコピーするには以下のようにします。

// ファイルを読み込み
vtkNew<vtkUnstructuredGridReader> reader;
const std::string filename = "mesh.vtk";
reader->SetFileName(filename.c_str());
reader->Update();                          // ファイルを読み込み
auto grid = reader->GetOutput();           // gridの型はvtkUnstructuredGrid*

tetra_mesh mesh;

const auto npoints = grid->GetNumberOfPoints(); // npointsの方はvtkIdType（整数型）
auto &points = mesh.points;
points.resize(npoints);
// 頂点データをコピー
for (vtkIdType i = 0; i < npoints; ++i) {
  auto p = grid->GetPoint(i); // pの型はdouble*
  auto &pt = points[i];
  pt.x = p[0];
  pt.y = p[1];
  pt.z = p[2];
}

const auto ncells = grid->GetNumberOfCells();
auto &cells = mesh.cells;
cells.resize(ncells);
// セルデータをコピー
for (vtkIdType i = 0; i < ncells; ++i) {
  auto cellpts = grid->GetCell(i)->GetPointIds(); // cellptsの型はvtkIdList*
  auto &cell = cells[i];
  for (std::size_t j = 0; j < cell.size(); ++j) {
    cell[j] = cellpts->GetId(j);
  }
}

他のグリッド形式の場合は、Legacy形式はvtk<GridName>Readerを、XML形式はvtkXML<GridName>Readerを使ってVTK File Formatを読み込みます。

逆にグリッドデータをVTKファイルに書き込む場合は以下のようにします。

vtkNew<vtkUnstructuredGrid> grid;
vtkNew<vtkPoints> vtkpoints;
vtkpoints->SetDataTypeToDouble();
vtkpoints->SetNumberOfPoints(mesh.n_points());

const auto& points = mesh.points;
// 頂点データをコピー
for (std::size_t i = 0; i < points.size(); ++i) {
  const auto& p = points[i];
  vtkpoints->SetPoint(i, p.x, p.y, p.z);
}
grid->SetPoints(vtkpoints);

const auto& cells = mesh.cells;

vtkNew<vtkCellArray> vtkcells;
constexpr int ncellpts = 4;
vtkcells->Allocate(vtkcells->EstimateSize(mesh.n_cells(), ncellpts));

vtkNew<vtkIdList> vtkcellpts;
vtkcellpts->SetNumberOfIds(4);
// セルデータをコピー
for (std::size_t i = 0; i < cells.size(); ++i) {
  const auto& cell = cells[i];
  for (std::size_t j = 0; j < cell.size(); ++i) {
    vtkcellpts->SetId(j, cell[j]);
  }
  vtkcells->InsertNextCell(vtkcellpts.GetPointer());
}
// セルデータをグリッドに渡す
// VTK_TETRAはセル形式を表すマクロ定数（vtkCellType.hに定義されている）
grid->SetCells(VTK_TETRA, vtkcells);

// ファイルへ書き込み
vtkNew<vtkUnstructuredGridWriter> writer;
writer->SetInputData(grid);
const std::string filename = "mesh.vtk";
writer->SetFileName(filename.c_str());
writer->SetFileTypeToBinary(); //  ASCIIの場合はSetFileTypeToASCII()を使う
writer->Write();

フィールドデータの入出力

既にファイルを読み込んでいるものとします。フィールドデータを取得するには以下のようにします。

// セルに紐付けられたフィールドデータを取得
auto cell_data = grid->GetCellData(); // cell_dataの型はvtkCellData*

const std::string fieldname = "pressure";
// セルに紐付けられたフィールドデータを名前から探す
if (cell_data->HasArray(fieldname.c_str())) {
  // 全てのフィールドデータはvtkAbstractArrayへのポインタとして保持されているため、
  // 要素にアクセスする際は必ず各具体的なarrayへダウンキャストが必要。
  auto array  = cell_data->GetAbstractArray(fieldname.c_str()); // arrayの型はvtkAbstractArray*
  auto darray = vtkDoubleArray::SafeDownCast(array);            // darrayの方はvtkDoubleArray*
  // あとはvtkDoubleArrayに用意されているメンバ関数を使って要素を取得するだけ

  // フィールドデータのサイズは成分数(number of components)x要素数(number of tuples)
  // セルに紐づけられたフィールドデータの要素数は当然ながらセル数と等しい
  const auto ntuples = darray->GetNumberOfTuples();      // = number of cells
  // 成分数は、スカラ:1、ベクトル:3、対称テンソル:6、テンソル:9
  const auto ncomp   = darray->GetNumberOfComponents();

  // 各要素の最初の成分へのポインタ
  for (int i = 0; i < ntuples; ++i) {
    auto p = darray->GetTuple(i);  // pはdouble*型
    // ...
  }
}

// 頂点に紐付けられたデータ
auto point_data = grid->GetPointData();
// 頂点に紐付けられたフィールドデータも、セルの場合と同様に探すことができる。
// 要素数は当然頂点の数と等しい。
// ...

C++以外の言語からの利用

Julia/FortranにはVTK形式を入出力するオープンソースライブラリがあります。

Pythonならpip install vtkしてimport vtkすればVTKライブラリ自体を使うことができます。

最後に

もっと色々書こうと思ったのですが、力尽きました…。気が向いたら加筆します。
VTK形式の入出力はVTKライブラリが提供する機能のほんの一部です。非常に多くの機能を提供していますので、興味を持ったら自分でコードを書いて遊んでみてください。

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up