Convert Mesh to Unstructered Grid#

Synopsis#

Convert an itk::Mesh to a vtkUnstructuredGrid.

Results#

An Output.vtu file will be generated alone with the following output.

Output:

Unstructered grid has 3 cells.

Code#

C++#

// ITK
#include "itkLineCell.h"
#include "itkMesh.h"
#include "itkTriangleCell.h"
#include "itkQuadrilateralCell.h"

// VTK
#include "vtkVersion.h"
#include <vtkCellArray.h>
#include <vtkSmartPointer.h>
#include <vtkUnstructuredGrid.h>
#include <vtkXMLUnstructuredGridWriter.h>

using MeshType = itk::Mesh<float, 3>;

// Functions
static MeshType::Pointer
CreateMeshWithEdges();
static void
ConvertMeshToUnstructuredGrid(MeshType::Pointer, vtkUnstructuredGrid *);

class VisitVTKCellsClass
{
  vtkCellArray * m_Cells;
  int *          m_LastCell;
  int *          m_TypeArray;

public:
  using CellInterfaceType = itk::CellInterface<MeshType::PixelType, MeshType::CellTraits>;

  using floatLineCell = itk::LineCell<CellInterfaceType>;
  using floatTriangleCell = itk::TriangleCell<CellInterfaceType>;
  using floatQuadrilateralCell = itk::QuadrilateralCell<CellInterfaceType>;

  // Set the vtkCellArray that will be constructed
  void
  SetCellArray(vtkCellArray * a)
  {
    m_Cells = a;
  }

  // Set the cell counter pointer
  void
  SetCellCounter(int * i)
  {
    m_LastCell = i;
  }

  // Set the type array for storing the vtk cell types
  void
  SetTypeArray(int * i)
  {
    m_TypeArray = i;
  }

  // Visit a triangle and create the VTK_TRIANGLE cell
  void
  Visit(unsigned long, floatTriangleCell * t)
  {
    m_Cells->InsertNextCell(3, (vtkIdType *)t->PointIdsBegin());
    m_TypeArray[*m_LastCell] = VTK_TRIANGLE;
    (*m_LastCell)++;
  }

  // Visit a triangle and create the VTK_QUAD cell
  void
  Visit(unsigned long, floatQuadrilateralCell * t)
  {
    m_Cells->InsertNextCell(4, (vtkIdType *)t->PointIdsBegin());
    m_TypeArray[*m_LastCell] = VTK_QUAD;
    (*m_LastCell)++;
  }

  // Visit a line and create the VTK_LINE cell
  void
  Visit(unsigned long, floatLineCell * t)
  {
    m_Cells->InsertNextCell(2, (vtkIdType *)t->PointIdsBegin());
    m_TypeArray[*m_LastCell] = VTK_LINE;
    (*m_LastCell)++;
  }
};

int
main()
{
  MeshType::Pointer mesh = CreateMeshWithEdges();

  vtkSmartPointer<vtkUnstructuredGrid> unstructuredGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
  ConvertMeshToUnstructuredGrid(mesh, unstructuredGrid);

  // Write file
  vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
  writer->SetFileName("output.vtu");
#if VTK_MAJOR_VERSION <= 5
  writer->SetInputConnection(unstructuredGrid->GetProducerPort());
#else
  writer->SetInputData(unstructuredGrid);
#endif
  writer->Write();

  return EXIT_SUCCESS;
}

MeshType::Pointer
CreateMeshWithEdges()
{

  auto mesh = MeshType::New();

  // Create 4 points and add them to the mesh
  MeshType::PointType p0, p1, p2, p3;

  p0[0] = -1.0;
  p0[1] = -1.0;
  p0[2] = 0.0;
  p1[0] = 1.0;
  p1[1] = -1.0;
  p1[2] = 0.0;
  p2[0] = 1.0;
  p2[1] = 1.0;
  p2[2] = 0.0;
  p3[0] = 1.0;
  p3[1] = 1.0;
  p3[2] = 1.0;

  mesh->SetPoint(0, p0);
  mesh->SetPoint(1, p1);
  mesh->SetPoint(2, p2);
  mesh->SetPoint(3, p3);

  // Create three lines and add them to the mesh
  using CellAutoPointer = MeshType::CellType::CellAutoPointer;
  using LineType = itk::LineCell<MeshType::CellType>;

  CellAutoPointer line0;
  line0.TakeOwnership(new LineType);
  line0->SetPointId(0, 0); // line between points 0 and 1
  line0->SetPointId(1, 1);
  mesh->SetCell(0, line0);

  CellAutoPointer line1;
  line1.TakeOwnership(new LineType);
  line1->SetPointId(0, 1); // line between points 1 and 2
  line1->SetPointId(1, 2);
  mesh->SetCell(1, line1);

  CellAutoPointer line2;
  line2.TakeOwnership(new LineType);
  line2->SetPointId(0, 2); // line between points 2 and 3
  line2->SetPointId(1, 3);
  mesh->SetCell(2, line2);

  return mesh;
}

void
ConvertMeshToUnstructuredGrid(MeshType::Pointer mesh, vtkUnstructuredGrid * unstructuredGrid)
{
  // Get the number of points in the mesh
  int numPoints = mesh->GetNumberOfPoints();
  if (numPoints == 0)
  {
    mesh->Print(std::cerr);
    std::cerr << "no points in Grid " << std::endl;
    exit(-1);
  }

  // Create the vtkPoints object and set the number of points
  vtkPoints * vpoints = vtkPoints::New();
  vpoints->SetNumberOfPoints(numPoints);
  // Iterate over all the points in the itk mesh filling in
  // the vtkPoints object as we go
  MeshType::PointsContainer::Pointer points = mesh->GetPoints();

  // In ITK the point container is not necessarily a vector, but in VTK it is
  vtkIdType                VTKId = 0;
  std::map<vtkIdType, int> IndexMap;

  for (MeshType::PointsContainer::Iterator i = points->Begin(); i != points->End(); ++i, VTKId++)
  {
    // Get the point index from the point container iterator
    IndexMap[VTKId] = i->Index();

    // Set the vtk point at the index with the the coord array from itk
    // itk returns a const pointer, but vtk is not const correct, so
    // we have to use a const cast to get rid of the const
    vpoints->SetPoint(VTKId, const_cast<float *>(i->Value().GetDataPointer()));
  }

  // Set the points on the vtk grid
  unstructuredGrid->SetPoints(vpoints);

  // Setup some VTK things
  int    vtkCellCount = 0; // running counter for current cell being inserted into vtk
  int    numCells = mesh->GetNumberOfCells();
  auto * types = new int[numCells]; // type array for vtk
  // create vtk cells and estimate the size
  vtkCellArray * cells = vtkCellArray::New();
  cells->EstimateSize(numCells, 4);

  // Setup the line visitor
  using LineVisitor = itk::CellInterfaceVisitorImplementation<
    float,
    MeshType::CellTraits,
    itk::LineCell<itk::CellInterface<MeshType::PixelType, MeshType::CellTraits>>,
    VisitVTKCellsClass>;
  auto lv = LineVisitor::New();
  lv->SetTypeArray(types);
  lv->SetCellCounter(&vtkCellCount);
  lv->SetCellArray(cells);

  // Setup the triangle visitor
  using TriangleVisitor = itk::CellInterfaceVisitorImplementation<
    float,
    MeshType::CellTraits,
    itk::TriangleCell<itk::CellInterface<MeshType::PixelType, MeshType::CellTraits>>,
    VisitVTKCellsClass>;
  auto tv = TriangleVisitor::New();
  tv->SetTypeArray(types);
  tv->SetCellCounter(&vtkCellCount);
  tv->SetCellArray(cells);

  // Setup the quadrilateral visitor
  using QuadrilateralVisitor = itk::CellInterfaceVisitorImplementation<
    float,
    MeshType::CellTraits,
    itk::QuadrilateralCell<itk::CellInterface<MeshType::PixelType, MeshType::CellTraits>>,
    VisitVTKCellsClass>;
  auto qv = QuadrilateralVisitor::New();
  qv->SetTypeArray(types);
  qv->SetCellCounter(&vtkCellCount);
  qv->SetCellArray(cells);

  // Add the visitors to a multivisitor

  MeshType::CellType::MultiVisitor::Pointer mv = MeshType::CellType::MultiVisitor::New();

  mv->AddVisitor(tv);
  mv->AddVisitor(qv);
  mv->AddVisitor(lv);

  // Now ask the mesh to accept the multivisitor which
  // will Call Visit for each cell in the mesh that matches the
  // cell types of the visitors added to the MultiVisitor
  mesh->Accept(mv);

  // Now set the cells on the vtk grid with the type array and cell array
  unstructuredGrid->SetCells(types, cells);
  std::cout << "Unstructured grid has " << unstructuredGrid->GetNumberOfCells() << " cells." << std::endl;

  // Clean up vtk objects
  cells->Delete();
  vpoints->Delete();
}

Classes demonstrated#

template<typename TPixelType, unsigned int VDimension = 3, typename TMeshTraits = DefaultStaticMeshTraits<TPixelType, VDimension, VDimension>> class Mesh : public itk::PointSet<TPixelType, VDimension, TMeshTraits>

Implements the N-dimensional mesh structure.

Mesh is an adaptive, evolving structure. Typically points and cells are created, with the cells referring to their defining points. If additional topological information is required, then BuildCellLinks() is called and links from the points back to the cells that use them are created. This allows implicit topological information about the faces and edges of the cells to be determined. (For example, a “face” neighbor to a cell can be determined by intersection the sets of cells that use the points defining the face. This is an inherent assumption on the manifold relationship of the cells in the mesh.) In some cases, either because the mesh is non-manifold, because we wish to explicitly store information with the faces and edges of the mesh, or because performance requirements demand that boundaries are explicitly represented (the set intersection does not need to be performed); then Mesh can be further extended by adding explicit boundary assignments.

Overview: Mesh implements the N-dimensional mesh structure for ITK. It provides an API to perform operations on points, cells, boundaries, etc., but does not tie down the underlying implementation and storage. A “MeshTraits” structure is used to define the container and identifier types that will be used to access the mesh. See DefaultStaticMeshTraits for the set of type definitions needed. All types that are defined in the “MeshTraits” structure will have duplicate type alias in the resulting mesh itself.

One of the most important parts of using this mesh is how to create cells to insert into it. The cells for the mesh take two template parameters. The first is the pixel type, and should correspond exactly to that type given to the mesh. The second is a “CellTraits” which holds a sub-set of the “MeshTraits” structure definitions, and is also a member of them. Any cell which is to be inserted to a mesh should have MeshTraits::CellTraits as its second template parameter.

Usage: Mesh has three template parameters. The first is the pixel type, or the type of data stored (optionally) with points, cells, and/or boundaries. The second is the geometric dimension of the points defining the mesh. This also limits the maximum topological dimension of the cells that can be inserted. The third template parameter is the “MeshTraits” structure controlling type information for the mesh. Most users will be happy with the defaults, and will not have to worry about this third argument.

Template parameters for Mesh:

TPixelType = The type stored as data for an entity (cell, point, or boundary).

TMeshTraits = Type information structure for the mesh.

References

No reference information is available.

ITK Sphinx Examples: