Compute Geodesic Distance on Mesh#

Synopsis#

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

Results#

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>;
  auto 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;

  auto trial = NodePairContainerType::New();

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

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

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

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

  try
  {
    writer->Update();
  }
  catch (const 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.