Examples/DataRepresentation/Mesh/PointSetWithCovariantVectors.cxx

/*=========================================================================
 *
 *  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
//
//  It is common to represent geometric objects by using points on their
//  surfaces and normals associated with those points.  This structure can be
//  easily instantiated with the \doxygen{PointSet} class.
//
//  The natural class for representing normals to surfaces and
//  gradients of functions is the \doxygen{CovariantVector}. A
//  covariant vector differs from a vector in the way it behaves
//  under affine transforms, in particular under anisotropic
//  scaling. If a covariant vector represents the gradient of a
//  function, the transformed covariant vector will still be the valid
//  gradient of the transformed function, a property which would not
//  hold with a regular vector.
//
//  \index{itk::PointSet!itk::CovariantVector}
//  \index{itk::CovariantVector!itk::PointSet}
//
//  The following example demonstrates how a \code{CovariantVector} can
//  be used as the \code{PixelType} for the \code{PointSet} class.  The
//  example illustrates how a deformable model could move under
//  the influence of the gradient of a potential function.
//
//  In order to use the CovariantVector class it is necessary to
//  include its header file along with the header of the point set.
//
//  \index{itk::CovariantVector!Header}
//
//  Software Guide : EndLatex
 
// Software Guide : BeginCodeSnippet
#include "itkCovariantVector.h"
#include "itkPointSet.h"
// Software Guide : EndCodeSnippet
 
int
main(int, char *[])
{
  //  Software Guide : BeginLatex
  //
  //  The CovariantVector class is templated over the type used to
  //  represent the spatial coordinates and over the space dimension.  Since
  //  the PixelType is independent of the PointType, we are free to select any
  //  dimension for the covariant vectors to be used as pixel type. However,
  //  we want to illustrate here the spirit of a deformable model. It is then
  //  required for the vectors representing gradients to be of the same
  //  dimension as the points in space.
  //
  //  \index{itk::CovariantVector!Instantiation}
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  constexpr unsigned int Dimension = 3;
  using PixelType = itk::CovariantVector<float, Dimension>;
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  Then we use the PixelType (which are actually CovariantVectors) to
  //  instantiate the PointSet type and subsequently create a PointSet object.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using PointSetType = itk::PointSet<PixelType, Dimension>;
  auto pointSet = PointSetType::New();
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The following code generates a circle and assigns gradient values to
  //  the points. The components of the CovariantVectors in this example are
  //  computed to represent the normals to the circle.
  //
  //  \index{itk::PointSet!SetPoint()}
  //  \index{itk::PointSet!SetPointData()}
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  PointSetType::PixelType gradient;
  PointSetType::PointType point;
 
  unsigned int     pointId = 0;
  constexpr double radius = 300.0;
 
  for (unsigned int i = 0; i < 360; ++i)
  {
    const double angle = i * std::atan(1.0) / 45.0;
    point[0] = radius * std::sin(angle);
    point[1] = radius * std::cos(angle);
    point[2] = 1.0; // flat on the Z plane
    gradient[0] = std::sin(angle);
    gradient[1] = std::cos(angle);
    gradient[2] = 0.0; // flat on the Z plane
    pointSet->SetPoint(pointId, point);
    pointSet->SetPointData(pointId, gradient);
    pointId++;
  }
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  We can now visit all the points and use the vector on the pixel values
  //  to apply a deformation on the points by following the gradient of the
  //  function. This is along the spirit of what a deformable model could do
  //  at each one of its iterations. To be more formal we should use the
  //  function gradients as forces and multiply them by local stress tensors
  //  in order to obtain local deformations.  The resulting deformations
  //  would finally be used to apply displacements on the points.  However,
  //  to shorten the example, we will ignore this complexity for the moment.
  //
  //  \index{itk::PointSet!PointDataIterator}
  //
  //  Software Guide : EndLatex
 
 
  // Software Guide : BeginCodeSnippet
  using PointDataIterator = PointSetType::PointDataContainer::ConstIterator;
  PointDataIterator       pixelIterator = pointSet->GetPointData()->Begin();
  const PointDataIterator pixelEnd = pointSet->GetPointData()->End();
 
  using PointIterator = PointSetType::PointsContainer::Iterator;
  PointIterator pointIterator = pointSet->GetPoints()->Begin();
  PointIterator pointEnd = pointSet->GetPoints()->End();
 
  while (pixelIterator != pixelEnd && pointIterator != pointEnd)
  {
    point = pointIterator.Value();
    gradient = pixelIterator.Value();
    for (unsigned int i = 0; i < Dimension; ++i)
    {
      point[i] += gradient[i];
    }
    pointIterator.Value() = point;
    ++pixelIterator;
    ++pointIterator;
  }
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The CovariantVector class does not overload the \code{+}
  //  operator with the \doxygen{Point}. In other words, CovariantVectors can
  //  not be added to points in order to get new points. Further, since we
  //  are ignoring physics in the example, we are also forced to do the
  //  illegal addition manually between the components of the gradient and
  //  the coordinates of the points.
  //
  //  Note that the absence of some basic operators on the ITK geometry
  //  classes is completely intentional with the aim of preventing the
  //  incorrect use of the mathematical concepts they represent.
  //
  //  \index{itk::CovariantVector}
  //
  //  Software Guide : EndLatex
  //
 
 
  //  We can finally visit all the points and print out the new values.
  //
  pointIterator = pointSet->GetPoints()->Begin();
  pointEnd = pointSet->GetPoints()->End();
  while (pointIterator != pointEnd)
  {
    std::cout << pointIterator.Value() << std::endl;
    ++pointIterator;
  }
 
 
  return EXIT_SUCCESS;
}