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

0．はじめに

以前の記事では、vtkImplictitPolyDataDistanceを利用して閉曲面の内外判定を行い解析用構造格子を生成しました。
今回は別のVTKのフィルタを利用して内外判定してメッシュを作ってみたいと思います。
また、記事内のソースコードや入力データはこちらにまとめています。

なお、内外判定以外のロジックは以前作成した記事をご参照ください。

1．vtkSelectEnclosedPoitntsによる内外判定

vtkSelectEnclosedPoitntsは入力に点群と閉曲面の二つをポリゴンデータとして与えると、点群が閉曲面ポリゴン内にあるか外にあるかを判定してくれるものです。(なぜ先にこちらを説明しなかったのか)

なお、入力するポリゴンデータは閉曲面(ポリゴン中に穴がない)で且つ要素同士が向かいあっている必要があります。
この二つの条件を満たしていれば、恐らくvtkImplictitPolyDataDistanceを使った手法より高い精度で内外判定が可能です。

2．実装例

以下に実装例を示します。
ソースコード中のGenerateBaseMeshFromAABBというメソッドは以前の記事からそのままスニペットしています。

sample2.cpp

#include <iostream>

#include <vtkSelectEnclosedPoints.h>

#include <vtkGeometryFilter.h>
#include <vtkUnstructuredGrid.h>
#include <vtkCellCenters.h>
#include <vtkSphereSource.h>
#include <vtkPolyData.h>
#include <vtkPointData.h>
#include <vtkCellData.h>
#include <vtkDataArray.h>
#include <vtkThreshold.h>
#include <vtkStructuredGrid.h>
#include <vtkAppendFilter.h>

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

#include <vtkSmartPointer.h>

void GenerateBaseMeshFromAABB(vtkDataSet *out_baseMesh, vtkPolyData *in_polyData, const int in_cellDims[3]);

int main()
{
  // ① 閉曲面ポリゴンデータを作る ----
  auto reader = vtkSmartPointer<vtkDataSetReader>::New();
  reader->SetFileName("C:\\test\\ushi.vtk");
  //reader->SetFileName("C:\\test\\fillHole.vtk");
  reader->Update();
  // ----

  // ソリッドをサーフェスに変える。元々サーフェスならなんも変わらない
  auto geometry = vtkSmartPointer<vtkGeometryFilter>::New();  
  geometry->SetInputData(reader->GetOutput());
  geometry->Update();

  vtkPolyData *closedPoly = geometry->GetOutput(); // polydataのポインタで形状を管理(やらなくてもいい)
  // ----

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

  // ③ 点群データを作る ----
  // 全ボクセル中心座標をボクセル数だけ作る
  auto cellCenters = vtkSmartPointer<vtkCellCenters>::New();
  cellCenters->SetInputData(baseMesh);
  cellCenters->Update();  // 全ボクセル中心座標を取得

  vtkPolyData *polyPoints = cellCenters->GetOutput(); // polydataのポインタで形状を管理(やらなくてもいい)
  // ----

  // ④ vtkSelectEnclosedPointsで内外判定 ----
  auto selectEnclosedPts = vtkSmartPointer<vtkSelectEnclosedPoints>::New();
  selectEnclosedPts->SetInputData(polyPoints);    // 点群セット
  selectEnclosedPts->SetSurfaceData(closedPoly);  // 閉曲面セット
  selectEnclosedPts->SetTolerance(1e-5);
  selectEnclosedPts->Update();
  // ----

  // ⑤ フィルタ適用、フィルタが持つ形状にぶら下がっている
  // "SelectedPoints"という配列を参照してbaseMeshのCellDataに与える
  vtkDataArray *isInsideOrOutside = selectEnclosedPts->GetOutput()->GetPointData()->GetArray("SelectedPoints");
  baseMesh->GetCellData()->AddArray(isInsideOrOutside); // PointDataで配列を取得した後、CellDataにぶらさげるので注意

  // 形状のCellDataにぶら下がっている"SelectedPoints"という配列から情報を抽出する
  auto threshold = vtkSmartPointer<vtkThreshold>::New();
  threshold->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_CELLS, "SelectedPoints");
  threshold->SetInputData(baseMesh);
  threshold->ThresholdBetween(1, 1);  // (min, max) = (1, 1)
  threshold->Update();

  // ⑥ VTKレガシー形式で出力
  auto writer = vtkSmartPointer<vtkDataSetWriter>::New();
  writer->SetInputData(threshold->GetOutput());
  writer->SetFileName("C:\\test\\voxelMesh.vtk");
  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());
}

ソースコードを実行すると、指定ディレクトリにvoxelMesh.vtkというファイルが出力されますので、ParaViewに投げ込んでみると以下のようなボクセルボクセルメッシュが可視化されます。

3．穴が空いたポリゴンでボクセルメッシュを作りたい場合

例えば、以下のように足の部分が切られて穴のできた牛ポリゴンを使ってボクセルメッシュを作りたいとします。

穴が空いたポリゴンの場合、閉曲面を形成できないのでvtkFillHolesFilterを利用して穴を埋めてあげます。

ここで、穴埋め部分の要素の貼り方が悪いとBoundaryEdgesというものが出現し、内外判定がうまくいかない場合があります。
BoundaryEdgesは、ループを形成している線分群(要するに穴)があると出現します。穴埋めしたのに穴があるってなんだという感じですが、穴に要素を張った際に線分同士(要素間のエッジ部分)が重なるとループができてしまうわけです。

この場合、vtkCleanPolyDataを使ってや重複している線分を消去してあげる必要があります。

以下に、実装例を示します。

testFillHoles.cpp

#include <iostream>

#include <vtkFillHolesFilter.h>

#include <vtkGeometryFilter.h>
#include <vtkPolyData.h>
#include <vtkCleanPolyData.h>
#include <vtkPolyDataNormals.h>

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

#include <vtkSmartPointer.h>

int main()
{
  // ① 穴があるポリゴンデータを読み込む ----
  auto reader = vtkSmartPointer<vtkDataSetReader>::New();
  reader->SetFileName("C:\\test\\ushiThatCutLeg.vtk");
  reader->Update();

  // vtkUnstructuredGrid -> vtkPolyData
  auto geometry = vtkSmartPointer<vtkGeometryFilter>::New();
  geometry->SetInputData(reader->GetOutput());
  geometry->Update();

  vtkPolyData *closedPoly = geometry->GetOutput();
  // ----

  // ② 穴を埋める ----
  auto fillHole = vtkSmartPointer<vtkFillHolesFilter>::New();
  fillHole->SetInputData(closedPoly);
  fillHole->SetHoleSize(100.0); // 適当なサイズを入力してください
  fillHole->Update();
  // ----

  // ③ 重複節点を消す ----
  // 形状の大きさからトレランスを求める
  double tol = 0.0;
  double *bds = fillHole->GetOutput()->GetBounds();
  for (int id = 0; id < 3; id++)
    tol += abs(bds[id * 2 + 1] - bds[id * 2]);
  tol = tol / 3.0 * 0.001; // 形状の大きさの平均の0.1%

  auto clean = vtkSmartPointer<vtkCleanPolyData>::New();
  clean->SetInputData(fillHole->GetOutput());
  clean->SetTolerance(tol);
  clean->Update();
  // ----

  // ④ VTKレガシー形式で出力
  auto writer = vtkSmartPointer<vtkDataSetWriter>::New();
  writer->SetInputData(clean->GetOutput());
  writer->SetFileName("C:\\test\\fillHole.vtk");
  writer->Update();
}

上記のコードを実行すると、指定した大きさの穴を埋めてくれます。
なお、コード中の "トレランス" とは、点群がお互い非常に近い位置にある場合に節点同士の重なりをどれだけ許すかというパラメータです。
この値が大きいと、重複してると判断したい点群が一つにならない可能性があるので注意しましょう。

穴を埋めると以下のような具合になります。

4．内外判定がうまくいかない場合

4.1．ボクセルメッシュの周辺にゴミが漂う

sample2.cpp中のselectEnclosedPts->SetTolerance(1e-5);という命令をコメントアウトして実行すると、以下のようにボクセルメッシュ周辺にゴミが漂う現象が起こります。

上記のように、閉曲面ポリゴンを入力したにも関わらず内外判定がうまくいかない場合、vtkSelectEnclosedPointsで設定されているToleranceに問題がある可能性が考えられます。

vtkSelectEnclosedPoints内ではメンバにToleranceという値(デフォルトでは $ 10^{-3} $ )が定義されており、この値によってはうまく内外判定ができないことがあります。

ここでのToleranceは恐らく前述のvtkCleanPolyDataのものと同様と思われますが、この値によって判定に不具合が出るのは、以下のようなアルゴリズムを使用しているためと推測されます。

vtkSelectEnclosedPoints内の内外判定アルゴリズム

2019.01.26追記：vtkSelectEnclosedPoints内部の実装を読み直したら考察に誤りがあったので修正しました。

vtkSelectEnclosedPointsの内部では、「与点及び、与点から十分離れたランダムな位置にある点がなす線分がポリゴンに何度交差したかをカウントし内外判定を行う」という処理を行っています。

ここでToleranceの値が大きいと、線分が閉曲面スレスレに生成されるときなどに恐らく内外判定がうまくできなくなるものと考えられます。

実際に、SetToleranceメソッドを使ってデフォルト( $ 10^{-3} $ )から $ 10^{-5} $ に設定してみると上手く内外判定できます。

---

4.2．そもそもボクセルメッシュが生成されない

このケースは、閉曲面中にBoundaryEdges又はNonManifoldEdgesが存在している可能性が挙げられます。

まず、BoundaryEdgesは前述の3．穴が空いたポリゴンでボクセルメッシュを作りたい場合で説明した穴空いた部分になります。

次にNonManifoldEdgesですが、これは以下のように、ある一つのエッジを共有するセルが3つ以上の時に出る線分です。

この`NonManifoldEdges`が閉曲面上にある場合、VTKのフィルタで除去するのはなかなか難しかったりするので、CAD等で該当部分を消去するといった泥臭い作業が必要になってきます。

なお、NonManifoldEdgesがあるかどうか確認したい場合、ParaView上で [Filter] -> [Alphabetical] -> ([Extract Surface] -> ) [Feature Edges] を選択し、[Properties] 内の [Non-ManifoldEdges] にチェックを入れた状態で [Apply] を選択すると [Information] またはビューアー上で確認できます。BoundaryEdgesに関しても同様です。

5．おわりに

前回の記事とは別の方法で解析用構造格子を作ってみました。

そのうち、前回と今回の手法を両方取り扱えるフィルタをモジュール化して、VTKライクに使えるものを作りたいと思います。