Doxygen/html/Examples_2DataRepresentation_2Mesh_2MeshKComplex_8cxx-example.html

/*=========================================================================

 *

 *  Copyright NumFOCUS

 *

 *  Licensed under the Apache License, Version 2.0 (the "License");

 *  you may not use this file except in compliance with the License.

 *  You may obtain a copy of the License at

 *

 *         https://www.apache.org/licenses/LICENSE-2.0.txt

 *

 *  Unless required by applicable law or agreed to in writing, software

 *  distributed under the License is distributed on an "AS IS" BASIS,

 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.

 *  See the License for the specific language governing permissions and

 *  limitations under the License.

 *

 *=========================================================================*/


//  Software Guide : BeginLatex

//

//  The \doxygen{Mesh} class supports the representation of formal

//  topologies. In particular the concept of \emph{K-Complex} can be

//  correctly represented in the Mesh. An informal definition of K-Complex

//  may be as follows: a K-Complex is a topological structure in which for

//  every cell of dimension $N$, its boundary faces (which are cells of

//  dimension $N-1$) also belong to the structure.

//

//  This section illustrates how to instantiate a K-Complex structure using

//  the mesh. The example structure is composed of one tetrahedron, its four

//  triangle faces, its six line edges and its four vertices.

//

//  \index{itk::Mesh!K-Complex}

//

//  Software Guide : EndLatex


//  Software Guide : BeginLatex

//

//  The header files of all the cell types involved should be loaded along

//  with the header file of the mesh class.

//

//  \index{itk::LineCell!header}

//  \index{itk::VertexCell!header}

//  \index{itk::TriangleCell!header}

//  \index{itk::TetrahedronCell!header}

//

//  Software Guide : EndLatex


// Software Guide : BeginCodeSnippet

#include "itkMesh.h"

#include "itkLineCell.h"

#include "itkTetrahedronCell.h"

// Software Guide : EndCodeSnippet


int

main(int, char *[])

{

  //  Software Guide : BeginLatex

  //

  //  Then the PixelType is defined and the mesh type is instantiated with it.

  //  Note that the dimension of the space is three in this case.

  //

  //  \index{itk::Mesh!Instantiation}

  //  \index{itk::Mesh!PixelType}

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  using PixelType = float;

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

  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex

  //

  //  The cell type can now be instantiated using the traits

  //  taken from the Mesh.

  //

  //  \index{itk::LineCell!Instantiation}

  //  \index{itk::VertexCell!Instantiation}

  //  \index{itk::TriangleCell!Instantiation}

  //  \index{itk::TetrahedronCell!Instantiation}

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  using CellType = MeshType::CellType;

  using VertexType = itk::VertexCell<CellType>;

  using LineType = itk::LineCell<CellType>;

  using TriangleType = itk::TriangleCell<CellType>;

  using TetrahedronType = itk::TetrahedronCell<CellType>;

  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex

  //

  //  The mesh is created and the points associated with the vertices are

  //  inserted.  Note that there is an important distinction between the

  //  points in the mesh and the \doxygen{VertexCell} concept. A VertexCell

  //  is a cell of dimension zero. Its main difference as compared to a point

  //  is that the cell can be aware of neighborhood relationships with other

  //  cells. Points are not aware of the existence of cells. In fact, from

  //  the pure topological point of view, the coordinates of points in the

  //  mesh are completely irrelevant.  They may as well be absent from the

  //  mesh structure altogether.  VertexCells on the other hand are necessary

  //  to represent the full set of neighborhood relationships on the

  //  K-Complex.

  //

  //  The geometrical coordinates of the nodes of a regular tetrahedron can be

  //  obtained by taking every other node from a regular cube.

  //

  //  \index{itk::Mesh!New()}

  //  \index{itk::Mesh!SetPoint()}

  //  \index{itk::Mesh!PointType}

  //  \index{itk::Mesh!Pointer}

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  auto mesh = MeshType::New();


  MeshType::PointType point0;

  MeshType::PointType point1;

  MeshType::PointType point2;

  MeshType::PointType point3;


  point0[0] = -1;

  point0[1] = -1;

  point0[2] = -1;

  point1[0] = 1;

  point1[1] = 1;

  point1[2] = -1;

  point2[0] = 1;

  point2[1] = -1;

  point2[2] = 1;

  point3[0] = -1;

  point3[1] = 1;

  point3[2] = 1;


  mesh->SetPoint(0, point0);

  mesh->SetPoint(1, point1);

  mesh->SetPoint(2, point2);

  mesh->SetPoint(3, point3);

  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex

  //

  //  We proceed now to create the cells, associate them with the points and

  //  insert them on the mesh. Starting with the tetrahedron we write the

  //  following code.

  //

  //  \index{itk::AutoPointer!TakeOwnership()}

  //  \index{CellAutoPointer!TakeOwnership()}

  //  \index{CellType!creation}

  //  \index{itk::Mesh!SetCell()}

  //  \index{itk::TetrahedronCell!Instantiation}

  //  \index{itk::TetrahedronCell!SetPointId()}

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  CellType::CellAutoPointer cellpointer;


  cellpointer.TakeOwnership(new TetrahedronType);

  cellpointer->SetPointId(0, 0);

  cellpointer->SetPointId(1, 1);

  cellpointer->SetPointId(2, 2);

  cellpointer->SetPointId(3, 3);

  mesh->SetCell(0, cellpointer);

  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex

  //

  //  Four triangular faces are created and associated with the mesh now.

  //  The first triangle connects points {0,1,2}.

  //

  //  \index{itk::TriangleCell!Instantiation}

  //  \index{itk::TriangleCell!SetPointId()}

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  cellpointer.TakeOwnership(new TriangleType);

  cellpointer->SetPointId(0, 0);

  cellpointer->SetPointId(1, 1);

  cellpointer->SetPointId(2, 2);

  mesh->SetCell(1, cellpointer);

  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex

  //

  //  The second triangle connects points { 0, 2, 3 }.

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  cellpointer.TakeOwnership(new TriangleType);

  cellpointer->SetPointId(0, 0);

  cellpointer->SetPointId(1, 2);

  cellpointer->SetPointId(2, 3);

  mesh->SetCell(2, cellpointer);

  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex

  //

  //  The third triangle connects points { 0, 3, 1 }.

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  cellpointer.TakeOwnership(new TriangleType);

  cellpointer->SetPointId(0, 0);

  cellpointer->SetPointId(1, 3);

  cellpointer->SetPointId(2, 1);

  mesh->SetCell(3, cellpointer);

  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex

  //

  //  The fourth triangle connects points { 3, 2, 1 }.

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  cellpointer.TakeOwnership(new TriangleType);

  cellpointer->SetPointId(0, 3);

  cellpointer->SetPointId(1, 2);

  cellpointer->SetPointId(2, 1);

  mesh->SetCell(4, cellpointer);

  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex

  //

  //  Note how the \code{CellAutoPointer} is reused every time. Reminder: the

  //  \doxygen{AutoPointer} loses ownership of the cell when it is passed as

  //  an argument of the \code{SetCell()} method. The AutoPointer is attached

  //  to a new cell by using the \code{TakeOwnership()} method.

  //

  //  The construction of the K-Complex continues now with the creation of the

  //  six lines on the tetrahedron edges.

  //

  //  \index{itk::LineCell!Instantiation}

  //  \index{itk::LineCell!SetPointId()}

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  cellpointer.TakeOwnership(new LineType);

  cellpointer->SetPointId(0, 0);

  cellpointer->SetPointId(1, 1);

  mesh->SetCell(5, cellpointer);


  cellpointer.TakeOwnership(new LineType);

  cellpointer->SetPointId(0, 1);

  cellpointer->SetPointId(1, 2);

  mesh->SetCell(6, cellpointer);


  cellpointer.TakeOwnership(new LineType);

  cellpointer->SetPointId(0, 2);

  cellpointer->SetPointId(1, 0);

  mesh->SetCell(7, cellpointer);


  cellpointer.TakeOwnership(new LineType);

  cellpointer->SetPointId(0, 1);

  cellpointer->SetPointId(1, 3);

  mesh->SetCell(8, cellpointer);


  cellpointer.TakeOwnership(new LineType);

  cellpointer->SetPointId(0, 3);

  cellpointer->SetPointId(1, 2);

  mesh->SetCell(9, cellpointer);


  cellpointer.TakeOwnership(new LineType);

  cellpointer->SetPointId(0, 3);

  cellpointer->SetPointId(1, 0);

  mesh->SetCell(10, cellpointer);

  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex

  //

  //  Finally the zero dimensional cells represented by the

  //  \doxygen{VertexCell} are created and inserted in the mesh.

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  cellpointer.TakeOwnership(new VertexType);

  cellpointer->SetPointId(0, 0);

  mesh->SetCell(11, cellpointer);


  cellpointer.TakeOwnership(new VertexType);

  cellpointer->SetPointId(0, 1);

  mesh->SetCell(12, cellpointer);


  cellpointer.TakeOwnership(new VertexType);

  cellpointer->SetPointId(0, 2);

  mesh->SetCell(13, cellpointer);


  cellpointer.TakeOwnership(new VertexType);

  cellpointer->SetPointId(0, 3);

  mesh->SetCell(14, cellpointer);

  // Software Guide : EndCodeSnippet


  // Print out the number of points and the number of cells.

  std::cout << "# Points= " << mesh->GetNumberOfPoints() << std::endl;

  std::cout << "# Cell  = " << mesh->GetNumberOfCells() << std::endl;


  //  Software Guide : BeginLatex

  //

  //  At this point the Mesh contains four points and fifteen cells enumerated

  //  from 0 to 14.  The points can be visited using PointContainer iterators.

  //

  // \index{itk::Mesh!PointsContainer}

  // \index{itk::Mesh!PointsIterators}

  // \index{itk::Mesh!GetPoints()}

  // \index{PointsContainer!Begin()}

  // \index{PointsContainer!End()}

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  using PointIterator = MeshType::PointsContainer::ConstIterator;

  PointIterator       pointIterator = mesh->GetPoints()->Begin();

  const PointIterator pointEnd = mesh->GetPoints()->End();


  while (pointIterator != pointEnd)

  {

    std::cout << pointIterator.Value() << std::endl;

    ++pointIterator;

  }

  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex

  //

  //  The cells can be visited using CellsContainer iterators.

  //

  // \index{itk::Mesh!CellsContainer}

  // \index{itk::Mesh!CellsIterators}

  // \index{itk::Mesh!GetCells()}

  // \index{CellsContainer!Begin()}

  // \index{CellsContainer!End()}

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  using CellIterator = MeshType::CellsContainer::ConstIterator;


  CellIterator cellIterator = mesh->GetCells()->Begin();

  CellIterator cellEnd = mesh->GetCells()->End();


  while (cellIterator != cellEnd)

  {

    CellType * cell = cellIterator.Value();

    std::cout << cell->GetNumberOfPoints() << std::endl;

    ++cellIterator;

  }

  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex

  //

  //  Note that cells are stored as pointer to a generic cell type that is the

  //  base class of all the specific cell classes. This means that at this

  //  level we can only have access to the virtual methods defined in the

  //  \code{CellType}.

  //

  //  The point identifiers to which the cells have been associated can be

  //  visited using iterators defined in the \code{CellType} trait. The

  //  following code illustrates the use of the PointIdIterators. The

  //  \code{PointIdsBegin()} method returns the iterator to the first

  //  point-identifier in the cell.  The \code{PointIdsEnd()} method returns

  //  the iterator to the past-end point-identifier in the cell.

  //

  //  \index{CellType!PointIdsBegin()}

  //  \index{CellType!PointIdsEnd()}

  //  \index{CellType!PointIdIterator}

  //  \index{PointIdIterator}

  //  \index{PointIdsBegin()}

  //  \index{PointIdsEnd()}

  //

  //  Software Guide : EndLatex


  cellIterator = mesh->GetCells()->Begin();

  cellEnd = mesh->GetCells()->End();


  while (cellIterator != cellEnd)

  {

    CellType * cell = cellIterator.Value();


    std::cout << "cell with " << cell->GetNumberOfPoints();

    std::cout << " points   " << std::endl;


    // Software Guide : BeginCodeSnippet

    using PointIdIterator = CellType::PointIdIterator;


    PointIdIterator pointIditer = cell->PointIdsBegin();

    PointIdIterator pointIdend = cell->PointIdsEnd();


    while (pointIditer != pointIdend)

    {

      std::cout << *pointIditer << std::endl;

      ++pointIditer;

    }

    // Software Guide : EndCodeSnippet


    ++cellIterator;

  }


  //  Software Guide : BeginLatex

  //

  //  Note that the point-identifier is obtained from the iterator using the

  //  more traditional \code{*iterator} notation instead the \code{Value()}

  //  notation used by cell-iterators.

  //

  //  Software Guide : EndLatex


  //  Software Guide : BeginLatex

  //

  //  Up to here, the topology of the K-Complex is not completely defined

  //  since we have only introduced the cells. ITK allows the user to define

  //  explicitly the neighborhood relationships between cells. It is clear

  //  that a clever exploration of the point identifiers could have allowed a

  //  user to figure out the neighborhood relationships. For example, two

  //  triangle cells sharing the same two point identifiers will probably be

  //  neighbor cells. Some of the drawbacks on this implicit discovery of

  //  neighborhood relationships is that it takes computing time and that some

  //  applications may not accept the same assumptions. A specific case is

  //  surgery simulation. This application typically simulates bistoury cuts

  //  in a mesh representing an organ. A small cut in the surface may be made

  //  by specifying that two triangles are not considered to be neighbors any

  //  more.

  //

  //  Neighborhood relationships are represented in the mesh by the

  //  notion of \emph{BoundaryFeature}. Every cell has an internal list of

  //  cell-identifiers pointing to other cells that are considered to be its

  //  neighbors. Boundary features are classified by dimension. For example, a

  //  line will have two boundary features of dimension zero corresponding to

  //  its two vertices. A tetrahedron will have boundary features of dimension

  //  zero, one and two, corresponding to its four vertices, six edges and

  //  four triangular faces. It is up to the user to specify the connections

  //  between the cells.

  //

  //  \index{BoundaryFeature}

  //  \index{CellBoundaryFeature}

  //  \index{itk::Mesh!BoundaryFeature}

  //

  //  Let's take in our current example the tetrahedron cell that was

  //  associated with the cell-identifier \code{0} and assign to it the four

  //  vertices as boundaries of dimension zero. This is done by invoking the

  //  \code{SetBoundaryAssignment()} method on the Mesh class.

  //

  //  \index{itk::Mesh!SetBoundaryAssignment()}

  //  \index{SetBoundaryAssignment()!itk::Mesh}

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  MeshType::CellIdentifier cellId = 0; // the tetrahedron


  int dimension = 0; // vertices


  MeshType::CellFeatureIdentifier featureId = 0;


  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 11);

  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 12);

  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 13);

  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 14);

  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex

  //

  //  The \code{featureId} is simply a number associated with the sequence of

  //  the boundary cells of the same dimension in a specific cell. For

  //  example, the zero-dimensional features of a tetrahedron are its four

  //  vertices. Then the zero-dimensional feature-Ids for this cell will range

  //  from zero to three. The one-dimensional features of the tetrahedron are

  //  its six edges, hence its one-dimensional feature-Ids will range from

  //  zero to five. The two-dimensional features of the tetrahedron are its

  //  four triangular faces. The two-dimensional feature ids will then range

  //  from zero to three. The following table summarizes the use on indices

  //  for boundary assignments.

  //

  //  \begin{center}

  //  \begin{tabular}{ | c || c | c | c | }

  //  \hline

  //  Dimension & CellType & FeatureId range & Cell Ids \\ \hline\hline

  //  0 & VertexCell & [0:3] & \{11,12,13,14\} \\   \hline

  //  1 & LineCell & [0:5] & \{5,6,7,8,9,10\} \\   \hline

  //  2 & TriangleCell & [0:3] & \{1,2,3,4\} \\ \hline

  //  \end{tabular}

  //  \end{center}

  //

  //  In the code example above, the values of featureId range from zero to

  //  three. The cell identifiers of the triangle cells in this example are

  //  the numbers \{1,2,3,4\}, while the cell identifiers of the vertex cells

  //  are the numbers \{11,12,13,14\}.


  //

  //  Let's now assign one-dimensional boundary features of the tetrahedron.

  //  Those are the line cells with identifiers  \{5,6,7,8,9,10\}. Note that

  //  the feature identifier is reinitialized to zero since the count is

  //  independent for each dimension.

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  cellId = 0;    // still the tetrahedron

  dimension = 1; // one-dimensional features = edges

  featureId = 0; // reinitialize the count


  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 5);

  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 6);

  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 7);

  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 8);

  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 9);

  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 10);

  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex

  //

  //  Finally we assign the two-dimensional boundary features of the

  //  tetrahedron. These are the four triangular cells with identifiers

  //  \{1,2,3,4\}. The featureId is reset to zero since feature-Ids are

  //  independent on each dimension.

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  cellId = 0;    // still the tetrahedron

  dimension = 2; // two-dimensional features = triangles

  featureId = 0; // reinitialize the count


  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 1);

  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 2);

  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 3);

  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 4);

  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex

  //

  //  At this point we can query the tetrahedron cell for information about

  //  its boundary features. For example, the number of boundary features of

  //  each dimension can be obtained with the method

  //  \code{GetNumberOfBoundaryFeatures()}.

  //

  //  \index{itk::Mesh!CellFeatureCount}

  //  \index{itk::Mesh!GetNumberOfBoundaryFeatures()}

  //  \index{GetNumberOfBoundaryFeatures()!itk::Mesh}

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  cellId = 0; // still the tetrahedron


  MeshType::CellFeatureCount n0; // number of zero-dimensional features

  MeshType::CellFeatureCount n1; // number of  one-dimensional features

  MeshType::CellFeatureCount n2; // number of  two-dimensional features


  n0 = mesh->GetNumberOfCellBoundaryFeatures(0, cellId);

  n1 = mesh->GetNumberOfCellBoundaryFeatures(1, cellId);

  n2 = mesh->GetNumberOfCellBoundaryFeatures(2, cellId);

  // Software Guide : EndCodeSnippet


  std::cout << "Number of boundary features of the cellId= " << cellId

            << std::endl;

  std::cout << "Dimension 0 = " << n0 << std::endl;

  std::cout << "Dimension 1 = " << n1 << std::endl;

  std::cout << "Dimension 2 = " << n2 << std::endl;


  //  Software Guide : BeginLatex

  //

  //  The boundary assignments can be recovered with the method

  //  \code{GetBoundaryAssignment()}. For example, the zero-dimensional

  //  features of the tetrahedron can be obtained with the following code.

  //

  // \index{itk::Mesh!GetBoundaryAssignment()}

  // \index{GetBoundaryAssignment()!itk::Mesh}

  //

  //  Software Guide : EndLatex


  std::cout << "Boundary features of dimension 0 " << std::endl;


  // Software Guide : BeginCodeSnippet

  dimension = 0;

  for (unsigned int b0 = 0; b0 < n0; ++b0)

  {

    MeshType::CellIdentifier id;

    const bool               found =

      mesh->GetBoundaryAssignment(dimension, cellId, b0, &id);

    if (found)

      std::cout << id << std::endl;

  }

  // Software Guide : EndCodeSnippet


  dimension = 1;

  std::cout << "Boundary features of dimension " << dimension << std::endl;

  for (unsigned int b1 = 0; b1 < n1; ++b1)

  {

    MeshType::CellIdentifier id;

    const bool               found =

      mesh->GetBoundaryAssignment(dimension, cellId, b1, &id);

    if (found)

    {

      std::cout << id << std::endl;

    }

  }


  dimension = 2;

  std::cout << "Boundary features of dimension " << dimension << std::endl;

  for (unsigned int b2 = 0; b2 < n2; ++b2)

  {

    MeshType::CellIdentifier id;

    const bool               found =

      mesh->GetBoundaryAssignment(dimension, cellId, b2, &id);

    if (found)

    {

      std::cout << id << std::endl;

    }

  }


  //  Software Guide : BeginLatex

  //

  //  The following code illustrates how to set the edge boundaries for one of

  //  the triangular faces.

  //

  //  Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet

  cellId = 2;    // one of the triangles

  dimension = 1; // boundary edges

  featureId = 0; // start the count of features


  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 7);

  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 9);

  mesh->SetBoundaryAssignment(dimension, cellId, featureId++, 10);

  // Software Guide : EndCodeSnippet


  std::cout << "In cell Id = " << cellId << std::endl;

  std::cout << "Boundary features of dimension " << dimension;

  n1 = mesh->GetNumberOfCellBoundaryFeatures(dimension, cellId);

  std::cout << " = " << n1 << std::endl;

  for (unsigned int b1 = 0; b1 < n1; ++b1)

  {

    MeshType::CellIdentifier id;

    const bool               found =

      mesh->GetBoundaryAssignment(dimension, cellId, b1, &id);

    if (found)

    {

      std::cout << id << std::endl;

    }

  }


  return EXIT_SUCCESS;

}