ITK  4.13.0
Insight Segmentation and Registration Toolkit
Examples/DataRepresentation/Mesh/MeshKComplex.cxx
/*=========================================================================
*
* Copyright Insight Software Consortium
*
* 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
*
* http://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"
// 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
typedef float PixelType;
typedef itk::Mesh< PixelType, 3 > MeshType;
// 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
typedef MeshType::CellType CellType;
typedef itk::VertexCell< CellType > VertexType;
typedef itk::LineCell< CellType > LineType;
typedef itk::TriangleCell< CellType > TriangleType;
typedef itk::TetrahedronCell< CellType > TetrahedronType;
// 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
MeshType::Pointer mesh = MeshType::New();
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
typedef MeshType::PointsContainer::ConstIterator PointIterator;
PointIterator pointIterator = mesh->GetPoints()->Begin();
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
typedef MeshType::CellsContainer::ConstIterator CellIterator;
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
typedef CellType::PointIdIterator 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{GetBoundaryAssigment()}. 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;
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;
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;
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;
bool found = mesh->GetBoundaryAssignment( dimension, cellId, b1, &id );
if( found )
{
std::cout << id << std::endl;
}
}
return EXIT_SUCCESS;
}