[VTK]閉曲面ポリゴンから解析用構造格子を生成する

1. はじめに

---

今回は、VTKを使用して閉曲面から解析用構造格子を生成する方法について書きます。
やや複雑な形状を解析したい場合などに利用できるかもしれません。
また本記事の中で、VTKのクラスをいくつか使用しておりますが、具体的な使用方法に関して知りたい場合はGoogleなどで「(クラス名) example」などと検索してサンプルを探してみてください。

さらに、VTKの導入を考えている方は前回の記事を見ていただければ幸甚です。
なお、使用するVTKのバージョンは8.1.2であり、実装したプログラムはC++で書いてます。

1. 目標

---

この記事では、次のような閉曲面を作っているポリゴンデータを3次元構造格子(ボクセルメッシュ)に変換することを目指します。

具体的にはこういった閉曲面ポリゴンを

こうします。

1. ベースとなる構造格子を作る

---

まず、vtkPolyDataで表現される閉曲面に対して、それを包括するAABB(Axis-Aligned Bounding Box)を作り、さらにその中に構造格子(今回は等間隔格子)を切ります。

この構造格子に対してvtkCellCentersを使ってボクセル中心に節点を打ち込むと下準備は完了です。

1. 内外判定

---

ベースの構造格子における各ボクセルが閉曲面内にあるかどうか確認したいため、ボクセル中心点の内外判定を行います。
VTKによる内外判定の方法は色々あり、vtkOBBTreeやvtkCellLocatorを使う方法などもありますが、今回はvtkImplicitPolyDataDistanceを利用します。

vtkImplicitPolyDataDistanceは入力したポリゴンを陰関数とした時に、与点と一番近いポリゴン表面までの符号付き距離を測ってくれるものです。

閉曲面を構成しているポリゴンデータの質にもよりますが、大体の場合において、閉曲面内側に節点を与えたときに負の距離を返します。

ここで、負の距離が返されたときのボクセルIDを保存しておくと、ベースの構造格子から欲しい部分を抜き出すことができます。

1. ボクセルの抽出

---

vtkExtractCellsを利用し、3. 内外判定で作成したボクセルIDのリストからベースの構造格子から抽出したいボクセル群を決定します。

1. 実装

---

以上の方法に基づいてC++コードを書いてみました。なお、readerで読み取っている閉曲面データは、こちらのサンプル(http://www.vtk.org/files/release/5.0/vtkdata-5.0.2.zip) を利用しています。(コードがキレイでないのはすみません)

また、今回はVTKの閉曲面ポリゴンを入力していますが、vtkSTLReaderを利用するとSTLデータを入力することもできます。

sample.cpp

#include <iostream>
#include <vtkImplicitPolyDataDistance.h>

#include <vtkGeometryFilter.h>
#include <vtkPolyData.h>
#include <vtkPoints.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>
#include <vtkAppendFilter.h>
#include <vtkCellCenters.h>
#include <vtkIdList.h>
#include <vtkExtractCells.h>

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

#include <vtkSmartPointer.h>

#define WRITE_XML

using namespace std;

void GenerateBaseMeshFromAABB(vtkDataSet *out_baseMesh, vtkPolyData *in_polyData, const int in_cellDims[3]);
void ExtractVoxelMeshWithinClosedSurface(vtkDataSet *out_voxelMesh, vtkPoints *in_voxelCenterPoints, vtkDataSet *in_baseGrid, vtkPolyData *in_closedSurface);
const bool IsInsideOrOutside(const double *in_points, vtkPolyData *in_closedSurface, vtkImplicitPolyDataDistance *in_implicitPolyDataDistance);

int main()
{
  // 閉曲面読み取り
  auto reader = vtkSmartPointer<vtkDataSetReader>::New();
  reader->SetFileName("C:\\testGenerateVoxelMesh\\IDGH.vtk");
  //reader->SetFileName("C:\\testGenerateVoxelMesh\\rbc_001.vtk");
  reader->Update();

  // PolyDataに変換
  auto geometry = vtkSmartPointer<vtkGeometryFilter>::New();
  geometry->SetInputData(reader->GetOutput());
  geometry->Update();

  auto baseMesh = vtkSmartPointer<vtkUnstructuredGrid>::New();  // 基礎メッシュ
  
  const int cellDims[] = { 50, 50, 50 };
  GenerateBaseMeshFromAABB(baseMesh, geometry->GetOutput(), cellDims);

  {
#ifndef WRITE_XML
    auto writer = vtkSmartPointer<vtkDataSetWriter>::New();
    writer->SetFileName("C:\\testGenerateVoxelMesh\\BaseMesh.vtk");
#else
    auto writer = vtkSmartPointer<vtkXMLDataSetWriter>::New();
    writer->SetFileName("C:\\testGenerateVoxelMesh\\BaseMesh.vtu");
#endif
    writer->SetInputData(baseMesh);
    writer->Update();
  }

  // 全ボクセル中心座標を取得
  auto cellCenters = vtkSmartPointer<vtkCellCenters>::New();
  cellCenters->SetInputData(baseMesh);
  cellCenters->Update();

  vtkPoints *voxelCenterPoints = cellCenters->GetOutput()->GetPoints();

  // 閉曲面に囲まれたボクセル群の形状
  auto voxelsWithinClosedSurface = vtkSmartPointer<vtkUnstructuredGrid>::New(); 

  ExtractVoxelMeshWithinClosedSurface(voxelsWithinClosedSurface, voxelCenterPoints, baseMesh, geometry->GetOutput());

  {
#ifndef WRITE_XML
    auto writer = vtkSmartPointer<vtkDataSetWriter>::New();
    writer->SetFileName("C:\\testGenerateVoxelMesh\\VoxelsWithinClosedSurface.vtk");
#else
    auto writer = vtkSmartPointer<vtkXMLDataSetWriter>::New();
    writer->SetFileName("C:\\testGenerateVoxelMesh\\VoxelsWithinClosedSurface.vtu");
#endif
    writer->SetInputData(voxelsWithinClosedSurface);
    writer->Update();
  }
}
// 閉曲面のAABB(Axis-Aligned Bounding Box)から基礎メッシュを作成
// 第1引数に空の形状データを与えること
void GenerateBaseMeshFromAABB(vtkDataSet *out_baseMesh, vtkPolyData *in_polyData, const int in_cellDims[3])
{
  // 諸元
  auto bounds = in_polyData->GetBounds();
  const double meshPitches[3] =
  {
    (bounds[1] - bounds[0]) / in_cellDims[0],
    (bounds[3] - bounds[2]) / in_cellDims[1],
    (bounds[5] - bounds[4]) / in_cellDims[2]
  };
  const double mins[3] = { bounds[0], bounds[2], bounds[4] };

  // 等間隔格子点生成
  auto points = vtkSmartPointer<vtkPoints>::New();
  for (int zId = 0; zId < in_cellDims[2] + 1; zId++)
  {
    for (int yId = 0; yId < in_cellDims[1] + 1; yId++)
    {
      for (int xId = 0; xId < in_cellDims[0] + 1; xId++)
      {
        const double x = xId * meshPitches[0] + mins[0];
        const double y = yId * meshPitches[1] + mins[1];
        const double z = zId * meshPitches[2] + mins[2];
        points->InsertNextPoint(x, y, z);
      }
    }
  }

  // 閉曲面を包括するAABB内に構造格子を生成
  auto baseMesh = vtkSmartPointer<vtkStructuredGrid>::New();
  baseMesh->SetExtent(0, in_cellDims[0], 0, in_cellDims[1], 0, in_cellDims[2]);
  baseMesh->SetPoints(points);

  // 何かと使いにくいことがあるので、vtkStructuredGridをvtkUnstructuredGridに変換
  auto append = vtkSmartPointer<vtkAppendFilter>::New();
  append->AddInputData(baseMesh);
  append->Update();

  // 第1引数にコピー
  out_baseMesh->DeepCopy(append->GetOutput());
}

// 基礎メッシュ内のボクセルのうち、閉曲面以内にあるものだけを抽出
// 第1引数に空の形状データを与えること
void ExtractVoxelMeshWithinClosedSurface(vtkDataSet *out_voxelMesh, vtkPoints *in_voxelCenterPoints, vtkDataSet *in_baseGrid, vtkPolyData *in_closedSurface)
{
  // 符号付き距離を測る準備
  auto implicitPolyDataDistance =
    vtkSmartPointer<vtkImplicitPolyDataDistance>::New();
  implicitPolyDataDistance->SetInput(in_closedSurface);

  // 基礎メッシュ中の全ボクセルを参照
  auto cellIdList = vtkSmartPointer<vtkIdList>::New();
  //for (vtkIdType bgId = 0; bgId < in_baseGrid->GetNumberOfCells(); bgId++)
  for (vtkIdType bgId = 0; bgId < in_voxelCenterPoints->GetNumberOfPoints(); bgId++)
  {
    const double *currentCenter = in_voxelCenterPoints->GetPoint(bgId);
    if (IsInsideOrOutside(currentCenter, in_closedSurface, implicitPolyDataDistance))
      cellIdList->InsertNextId(bgId);
  }

  // cellIdListに登録したボクセルのみ抽出する
  auto extractCells = vtkSmartPointer<vtkExtractCells>::New();
  extractCells->SetInputData(in_baseGrid);
  extractCells->SetCellList(cellIdList);
  extractCells->Update();

  out_voxelMesh->DeepCopy(extractCells->GetOutput());
}

// 内外判定
const bool IsInsideOrOutside(const double *in_points, vtkPolyData *in_closedSurface, vtkImplicitPolyDataDistance *in_implicitPolyDataDistance)
{
  double distance = in_implicitPolyDataDistance->FunctionValue(in_points);

  if (distance <= 0)
    return true;
  else
    return false;
}

指定ディレクトリにReaderで入力している平曲面VTKを置き、プログラムを実行しうまくいくと"VoxelsWithinClosedSurface.vtk"という名前で目的の構造格子データが生成されます。

赤血球

また、こちらから"rbc_001.vtk"というデータを拾って同様に試してみます。
元形状は以下のような赤血球です。

これに対して構造格子を生成すると、以下のようになります。

ちゃんと赤血球がボクセル化されてます。

ウサギさん

元の形状

粗いウサギさんになりました。

1. おわりに

---

数値計算で使えそうなVTKの使い方の一例を紹介してみました。このようにVTKは割りと便利なので、よろしかったら皆さんも勉強してみてください。

P.S.
以上の方法でうまくボクセルメッシュが生成されない場合は、閉曲面の各要素に定義されている法線ベクトルの一部に、閉曲面外向きを正と設定できていないものがあるかもしれません。
法線ベクトルに問題がある場合、ベースメッシュの分割数によっては生成した構造格子の周辺にゴミが漂う場合があるので、ParaViewやVTKのConnectivityFilter等を利用して除去してください。

また、内外判定の精度を高めたい場合はボクセル内にさらに点を打つといいです。（どの点がどのボクセルに対応するかは構造体などを使って管理するといいかもしれません）