[VTK]ポリゴンに干渉してるソリッドを抽出する

0．はじめに

今回は、二つのポリゴンとソリッド同士でお互い干渉してる部分から、ソリッド部分のみを取り出すということをやってみます。

ただし、ソリッドの要素よりもポリゴンの要素がかなり細かくないと、以下のようにうまく抽出できないこともあるので注意してください。

1．ソリッドモデルとポリゴンを準備

どんな形状でも構いませんが、VTKのUnstructuredGrid(ここではソリッドに相当)の形状とPolyData(ポリゴンに相当)の形状を用意します。
(なお、コードを若干変更する必要がありますが、ポリゴンはSTL形式のものでも構いません)

2．vtkCellLocatorでポリゴンを構成する節点に一番近いソリッドの要素を探す

vtkCellLocatorというクラスのメンバを使うと、与えた二点間に挟まれる全要素を取得したり、ある直方体に収まる要素全てを取得することなどができます。
今回は指定半径内で与点に最も近い節点と要素IDを取得するFindClosestPointWithinRadiusというメソッドを利用します。

次にVTKのvtkCellLocatorクラスを活用した実装を示します。(適当な実装なのでお目汚しすみません)

実装

sample.cpp

#include <vtkCellLocator.h>
#include <vtkIdList.h>
#include <vtkIntArray.h>
#include <vtkCellData.h>
#include <vtkExtractCells.h>

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

#include <vtkSmartPointer.h>

int main()
{
  // (1) ポリゴン読み取り ----
#if 1
  auto poly = vtkSmartPointer<vtkDataSetReader>::New();
  poly->SetFileName("C:\\test\\sphere.vtk");
#else
  auto poly = vtkSmartPointer<vtkSTLReader>::New();
  poly->SetFileName("C:\\test\\sphere.stl");
#endif
  poly->Update();
  // ----

  // (2) ソリッド読み取り ----
  auto solid = vtkSmartPointer<vtkDataSetReader>::New();
  solid->SetFileName("C:\\test\\solid.vtk");
  solid->Update();
  // ----

  // (3) vtkCellLocatorでポリゴンに干渉してるソリッドの要素IDを拾う ----
  auto cellLocator = vtkSmartPointer<vtkCellLocator>::New();
  cellLocator->SetDataSet(solid->GetOutput());
  cellLocator->BuildLocator();
  
  const double radius = 1e-6;
  double closestPoint[3];
  double closestPointDist2; 
  vtkIdType cellId; 
  int subId; 

  auto cellIdList = vtkSmartPointer<vtkIdList>::New();
  
  // ポリゴンの全節点分だけ処理
  for (vtkIdType spId = 0; spId < poly->GetOutput()->GetNumberOfPoints(); spId++)
  {
    double *surfPts = poly->GetOutput()->GetPoint(spId);
    const int isFound = cellLocator->FindClosestPointWithinRadius
      (surfPts, radius, closestPoint, cellId, subId, closestPointDist2);
    
    if (isFound)
      cellIdList->InsertUniqueId(cellId); // 重複なしで拾ったボクセルメッシュの要素IDを保存
  }
  // ----

  // (4) ポリゴンがどの要素に干渉してるかという情報を追加したソリッドデータを書き出す ----
#if 1
  // IdList -> IntArray に変換
  auto idArray = vtkSmartPointer<vtkIntArray>::New();
  idArray->SetName("IsIntersected");
  idArray->SetNumberOfComponents(1);  // 1次元配列として使う
  idArray->SetNumberOfValues(solid->GetOutput()->GetNumberOfCells());
  idArray->FillValue(0);  // とりあえず全部0で初期化
  for (vtkIdType scId = 0; scId < cellIdList->GetNumberOfIds(); scId++)
  {
    const vtkIdType bdIdOfVoxels = cellIdList->GetId(scId);
    idArray->SetValue(bdIdOfVoxels, 1); // 干渉してるセルのみ1を与える
  }
  solid->GetOutput()->GetCellData()->AddArray(idArray); // 形状データ中のCellDataにぶら下げる

  auto writer = vtkSmartPointer<vtkDataSetWriter>::New();
  writer->SetFileName("C:\\test\\addIsIntersected.vtk");
  writer->SetInputConnection(solid->GetOutputPort());
  writer->Update();
#else
  // 形状のみ取りたいならこちらで
  auto extract = vtkSmartPointer<vtkExtractCells>::New();
  extract->SetCellList(cellIdList);
  extract->SetInputData(solid->GetOutput());
  extract->Update();

  auto writer = vtkSmartPointer<vtkDataSetWriter>::New();
  writer->SetFileName("C:\\test\\extractCells.vtk");
  writer->SetInputConnection(extract->GetOutputPort());
  writer->Update();
#endif
  // ----
}

これを実行すると、指定ディレクトリに以下のようなソリッドデータが出力されます。

データ中のCellDataに"IsIntersected"という名前の配列がぶら下がっているので、ParaViewのThresholdやVTKのvtkThresholdなどを利用するとポリゴンに干渉している領域のみ抽出することができます。

3．他の組み合わせもやってみる

閉曲面ポリゴンから作った構造格子 + 閉曲面ポリゴンの部分表面

以前の記事で閉曲面ポリゴンから構造格子を生成する方法(こちらとこちら)を紹介し、またParaViewによる3Dポリゴン上の滑らかに繋がった面だけを抽出する方法(こちら)を紹介しましたが、これらを用いて作ったソリッドとポリゴンを使ってみます。

まず、こちらまたはこちらの方法から、基礎メッシュ（ここでは606060）を作り閉曲面ポリゴンからボクセルメッシュを生成します。(なお、ここで使用してる閉曲面ポリゴンは上の方で利用していたソリッドのサーフェスだけをとってきたものになります)

さらに、こちらの方法を使って部分表面を閉曲面ポリゴンから抽出します。

MeshMixerを使って部分表面をリメッシュ

ソリッドとポリゴンが用意できたので、今回紹介した方法を利用してボクセルメッシュ側に部分表面の情報を付与してあげたいところなのですが、今のままだと部分表面側の解像度が悪いので以下のように局所的にボクセルメッシュ側に反映できない部分が出てきます。

そのため、今回はMeshMixerというAutoDesk社のフリーソフトを使って部分表面に対してリメッシュをかけてあげます。

まず、[ファイル] -> [インポート] からSTLファイル(VTKの形式だと読み込めないので、ParaViewなどでSTLファイルを作ってください)を選択し、部分表面のデータを読み込ませます。
また、デフォルトではワイヤフレームの表示がないので[ビュー] -> [ワイヤフレームの表示]からワイヤフレームを表示できるようにしましょう。

データの読み込み後、Ctrl + Aで形状モデル全体を選択状態にし、[編集] -> [再メッシュ]から均等性や密度の項目などを操作した後に適用を押すとリメッシュが完了するので、[ファイル] -> [エクスポート]からSTL形式で出力してください。

リメッシュ後のポリゴン、先ほど作ったボクセルメッシュを今回作成したコードに通します。(この時、ポリゴンについてはSTLを読み込むようにし、またradiusというパラメータを少し大きめに設定するようにコードを改変してください)
指定ディレクトリにVTK形式のデータが出力されますので、ParaViewでIsIntersectedのコンターを表示すると次のようにボクセルメッシュ側にポリゴンの位置が反映された図を見ることができます。

5．おわりに

今回は、ソリッドとポリゴンを組み合わせて、ポリゴンに干渉してるソリッドを抽出するというのをやってみました。
最後に紹介した例は境界条件の設定などに使えそうなので、もう少し簡単でインタラクティブにできたら尚よさそうです。