Compute Geodesic Distance on Mesh

Synopsis

Compute the geodesic distance from a provided seed vertex on a mesh.

Results

Input mesh

Output mesh

Interactive input mesh

Interactive output mesh

Code

C++

#include "itkFastMarchingQuadEdgeMeshFilterBase.h"
#include "itkQuadEdgeMeshExtendedTraits.h"
#include "itkRegularSphereMeshSource.h"
#include "itkFastMarchingThresholdStoppingCriterion.h"
#include "itkMeshFileReader.h"
#include "itkMeshFileWriter.h"

int
main(int argc, char * argv[])
{
  if (argc != 3)
  {
    std::cerr << "Usage: " << argv[0] << std::endl;
    std::cerr << " <input filename> <output filename>" << std::endl;
    return EXIT_FAILURE;
  }
  using PixelType = float;
  using CoordType = double;

  constexpr unsigned int Dimension = 3;

  using Traits = itk::QuadEdgeMeshExtendedTraits<PixelType, // type of data for vertices
                                                 Dimension, // geometrical dimension of space
                                                 2,         // Mac topological dimension of a cell
                                                 CoordType, // type for point coordinate
                                                 CoordType, // type for interpolation weight
                                                 PixelType, // type of data for cell
                                                 bool,      // type of data for primal edges
                                                 bool       // type of data for dual edges
                                                 >;

  using MeshType = itk::QuadEdgeMesh<PixelType, Dimension, Traits>;

  using ReaderType = itk::MeshFileReader<MeshType>;
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName(argv[1]);

  using FastMarchingType = itk::FastMarchingQuadEdgeMeshFilterBase<MeshType, MeshType>;

  MeshType::Pointer mesh = reader->GetOutput();

  MeshType::PointsContainerConstPointer points = mesh->GetPoints();

  MeshType::PointsContainerConstIterator pIt = points->Begin();
  MeshType::PointsContainerConstIterator pEnd = points->End();

  while (pIt != pEnd)
  {
    mesh->SetPointData(pIt->Index(), 1.);
    ++pIt;
  }

  using NodePairType = FastMarchingType::NodePairType;
  using NodePairContainerType = FastMarchingType::NodePairContainerType;

  NodePairContainerType::Pointer trial = NodePairContainerType::New();

  NodePairType nodePair(0, 0.);
  trial->push_back(nodePair);

  using CriterionType = itk::FastMarchingThresholdStoppingCriterion<MeshType, MeshType>;
  CriterionType::Pointer criterion = CriterionType::New();
  criterion->SetThreshold(100.);

  FastMarchingType::Pointer fmmFilter = FastMarchingType::New();
  fmmFilter->SetInput(mesh);
  fmmFilter->SetTrialPoints(trial);
  fmmFilter->SetStoppingCriterion(criterion);

  using WriterType = itk::MeshFileWriter<MeshType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetInput(fmmFilter->GetOutput());
  writer->SetFileName(argv[2]);

  try
  {
    writer->Update();
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

Classes demonstrated

template<typename TInput, typename TOutput>
class FastMarchingQuadEdgeMeshFilterBase : public itk::FastMarchingBase<TInput, TOutput>

Fast Marching Method on QuadEdgeMesh.

The speed function is specified by the input mesh. Data associated to each point is considered as the speed function. The speed function is set using the method SetInput().

If the speed function is constant and of value one, fast marching results is an approximate geodesic function from the initial alive points.

Implementation of this class is based on “Fast Marching Methods on Triangulated Domains”, Kimmel, R., and Sethian, J.A., Proc. Nat. Acad. Sci., 95, pp. 8341-8435, 1998.

See itk::FastMarchingQuadEdgeMeshFilterBase for additional documentation.