VTKで頂点距離を可視化

#include <iostream>

//VTK_MODULE_INITに必要
#include <vtkAutoInit.h>


#include <vtkSmartPointer.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>

//円筒とその表示に必要
#include <vtkPointData.h>
#include <vtkCellData.h>
#include <vtkPolyData.h>
#include <vtkActor.h>
#include <vtkPolyDataMapper.h>

#include <vtkLookupTable.h>
#include <vtkFloatArray.h>


#pragma comment(lib,"opengl32.lib")
#pragma comment(lib,"psapi.lib")
#pragma comment(lib,"dbghelp.lib")
#pragma comment(lib,"ws2_32.lib")

#include <vtkCellArray.h>
#include <vtkTriangle.h>

#include <open3d/Open3D.h>
#include <open3d/geometry/PointCloud.h>
#include <open3d/geometry/TriangleMesh.h>

// kdtree
#include <open3d/geometry/KDTreeFlann.h>

#pragma comment(lib,"Open3D.lib")

//必須
VTK_MODULE_INIT(vtkRenderingOpenGL2);
VTK_MODULE_INIT(vtkInteractionStyle);

// baseの各頂点に対して、targetの最も近い頂点までの距離を計算する
std::vector<double> calcNearestPoint(std::shared_ptr<open3d::geometry::TriangleMesh> base, std::shared_ptr<open3d::geometry::TriangleMesh> target) {

  std::vector<double> ret;
  // mymeshの頂点のkd-treeを作成
  open3d::geometry::KDTreeFlann kdtree;
  kdtree.SetGeometry(*target);
  
  // 頂点Pに最も近い頂点を探す
  std::vector<double> distances;
  std::vector<int> indices;
  for (auto& p : base->vertices_) {
    int count = kdtree.SearchKNN(p,1,indices,distances);
    double distance = sqrt(distances[0]);
    ret.push_back(distance);
  }

  return ret;
}

// テストデータを作成
void CreateMeshData_o3d_1(std::shared_ptr<open3d::geometry::TriangleMesh> mymesh, int xcount, int ycount, float width, float height,float offsetz, float depth, float density) {
  
  // 頂点の登録
  for (int i = 0; i < xcount; ++i) {
    for (int j = 0; j < ycount; ++j) {
      double x = -0.5 + 1.0 * i / (xcount - 1);
      double y = -0.5 + 1.0 * j / (ycount - 1);
      double z = sin(x * density) * depth + sin(y * density) * depth + offsetz;

      mymesh->vertices_.push_back(Eigen::Vector3d(x * width, y * height, z));

    }
  }
  // メッシュ（三角形）の作成
  for (int i = 0; i < xcount - 1; ++i) {
    for (int j = 0; j < ycount - 1; ++j) {
      // 四角形を2つの三角形に分割
      const int p1 = i * ycount + j;       // 左下
      const int p2 = (i + 1) * ycount + j; // 右下
      const int p3 = i * ycount + (j + 1); // 左上
      const int p4 = (i + 1) * ycount + (j + 1); // 右上

      // 三角形1 (p1, p2, p3)
      mymesh->triangles_.push_back(Eigen::Vector3i(p1, p2, p3));
      // 三角形2 (p2, p4, p3)
      mymesh->triangles_.push_back(Eigen::Vector3i(p2, p4, p3));

    }
  }

}

struct MeshData {
  std::shared_ptr<open3d::geometry::TriangleMesh> curve;
  std::shared_ptr<open3d::geometry::TriangleMesh> flat;
};

MeshData createData(int xcount,int ycount) {
  std::shared_ptr<open3d::geometry::TriangleMesh> mymesh1 = std::make_shared<open3d::geometry::TriangleMesh>();
  std::shared_ptr<open3d::geometry::TriangleMesh> mymesh2 = std::make_shared<open3d::geometry::TriangleMesh>();

  CreateMeshData_o3d_1(mymesh1, xcount, ycount, 2, 2, 0.0, 0.1, 10);
  CreateMeshData_o3d_1(mymesh2, xcount, ycount, 2, 2,-0.3, 0.05, 17);

  return MeshData{ mymesh1,mymesh2 };

}

      
vtkSmartPointer<vtkPolyData> CreateVTKMeshData(std::shared_ptr<open3d::geometry::TriangleMesh> tri,int xcount,int ycount) {
  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
  vtkSmartPointer<vtkCellArray> triangles = vtkSmartPointer<vtkCellArray>::New();
  vtkSmartPointer<vtkCellArray> vertices = vtkSmartPointer<vtkCellArray>::New();

  for(auto& p : tri->vertices_)
  {
    vtkIdType pid = points->InsertNextPoint(p.x(), p.y(), p.z());
    vertices->InsertNextCell(1, &pid);
  }

  // メッシュ（三角形）の作成
  for (int i = 0; i < xcount - 1; ++i) {
    for (int j = 0; j < ycount - 1; ++j) {
      // 四角形を2つの三角形に分割
      vtkIdType p1 = i * ycount + j;       // 左下
      vtkIdType p2 = (i + 1) * ycount + j; // 右下
      vtkIdType p3 = i * ycount + (j + 1); // 左上
      vtkIdType p4 = (i + 1) * ycount + (j + 1); // 右上

      // 三角形1 (p1, p2, p3)
      vtkSmartPointer<vtkTriangle> triangle1 = vtkSmartPointer<vtkTriangle>::New();
      triangle1->GetPointIds()->SetId(0, p1);
      triangle1->GetPointIds()->SetId(1, p2);
      triangle1->GetPointIds()->SetId(2, p3);
      triangles->InsertNextCell(triangle1);

      // 三角形2 (p2, p4, p3)
      vtkSmartPointer<vtkTriangle> triangle2 = vtkSmartPointer<vtkTriangle>::New();
      triangle2->GetPointIds()->SetId(0, p2);
      triangle2->GetPointIds()->SetId(1, p4);
      triangle2->GetPointIds()->SetId(2, p3);
      triangles->InsertNextCell(triangle2);
    }
  }

  // vtkPolyData の作成
  vtkSmartPointer<vtkPolyData> polyData = vtkSmartPointer<vtkPolyData>::New();
  polyData->SetPoints(points);
  polyData->SetPolys(triangles);  // 三角形セルを設定

  return polyData;
}

vtkSmartPointer<vtkActor> createActor(vtkSmartPointer<vtkPolyData> polyData, vtkSmartPointer<vtkLookupTable> table,double min_,double max_ ) {
  vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
  mapper->SetInputData(polyData);

  if( table != nullptr)
  {
    mapper->SetLookupTable(table);
    mapper->SetScalarRange(min_, max_);
    mapper->SetInterpolateScalarsBeforeMapping(true);
  }

  vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
  actor->SetMapper(mapper);

  return actor;
}

//! @brief サーモグラフィー風のカラールックアップテーブルを生成する
//! @param [in] scalar_min スカラー値の最小値
//! @param [in] scalar_max スカラー値の最大値
//! @param [in] division スカラー値の範囲を何分割するか
//! @return カラールックアップテーブル
vtkSmartPointer<vtkLookupTable> ThermographyColorTable(const double scalar_min, const double scalar_max, const int division) {

  vtkSmartPointer<vtkLookupTable> lut = vtkSmartPointer<vtkLookupTable>::New();
  lut->SetNumberOfTableValues(division); // テーブルの要素数を設定
  lut->SetTableRange(scalar_min, scalar_max); // スカラー値の範囲を設定
  lut->Build();

  // テーブルの i 番目の要素に対応する色を設定
  for (int i = 0; i < division; i++)
  {
    double t = static_cast<double>(i) / float(division - 1);
    double r, g, b;
    if (t < 0.33)
    {
      // 青から緑への補間
      double factor = t / 0.33;
      r = 0.0;
      g = factor;
      b = 1.0 - factor;
    }
    else if (t < 0.66)
    {
      // 緑から黄への補間
      double factor = (t - 0.33) / 0.33;
      r = factor;
      g = 1.0;
      b = 0.0;
    }
    else
    {
      // 黄から赤への補間
      double factor = (t - 0.66) / 0.34;
      r = 1.0;
      g = 1.0 - factor;
      b = 0.0;
    }

    // テーブルの i 番目の要素に色を設定
    lut->SetTableValue(
      i  // テーブルの要素番号
      , r, g, b, 1.0
    );
  }
  return lut;
}

int main(int /*argc*/, char** /*argv*/)
{
  auto meshes = createData(100, 100);
  std::vector<double>  scalar = calcNearestPoint(meshes.curve, meshes.flat);

  vtkSmartPointer<vtkPolyData> curve = CreateVTKMeshData(meshes.curve, 100, 100);
  vtkSmartPointer<vtkPolyData> flat = CreateVTKMeshData(meshes.flat, 100, 100);

  // scalarをvtkFloatArrayに変換
  vtkSmartPointer<vtkFloatArray> scalars = vtkSmartPointer<vtkFloatArray>::New();
  scalars->SetName("MyScalarValues");
  double _min = *std::min_element(scalar.begin(), scalar.end());
  double _max = *std::max_element(scalar.begin(), scalar.end());
  for (auto& s : scalar) {
    scalars->InsertNextValue(s);
  }
  curve->GetPointData()->AddArray(scalars);
  curve->GetPointData()->SetActiveScalars("MyScalarValues");

  vtkSmartPointer<vtkLookupTable> colortable = ThermographyColorTable(_min, _max, 100);
  vtkSmartPointer<vtkActor> actor1 = createActor(curve, colortable,_min,_max);
  vtkSmartPointer<vtkActor> actor2 = createActor(flat,nullptr,0,0);

  //////////////////////////////////////
  auto renderer = vtkSmartPointer<vtkRenderer>::New();
  renderer->AddActor(actor1);
  renderer->AddActor(actor2);
  renderer->ResetCamera();

  //////////////////////////////////////
  auto interactor = vtkSmartPointer<vtkRenderWindowInteractor>::New();

  //////////////////////////////////////
  auto renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
  renderWindow->AddRenderer(renderer);
  renderWindow->SetInteractor(interactor);
  renderWindow->Render();

  interactor->Start(); //イベントループへ入る

  return 0;
}