+ TVTKDistanceMap::Pointer vtk_distanceMap = TVTKDistanceMap::New( );
+ vtk_distanceMap->SetInput( distanceMap->GetOutput( ) );
+ vtk_distanceMap->Update( );
+
+ vtkSmartPointer< vtkImageMarchingCubes > mc =
+ vtkSmartPointer< vtkImageMarchingCubes >::New( );
+ mc->SetInputData( vtk_distanceMap->GetOutput( ) );
+ mc->SetValue( 0, 1e-2 );
+ mc->Update( );
+ view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
+ view.Render( );
+
+ // Extract paths
+ TDijkstra::Pointer paths = TDijkstra::New( );
+ paths->AddSeed( segmentor->GetSeed( 0 ), TScalar( 0 ) );
+ paths->SetInput( distanceMap->GetOutput( ) );
+ paths->SetNeighborhoodOrder( 1 );
+
+ if( visual_debug )
+ {
+ // Configure observer
+ TDijkstraObs::Pointer obs = TDijkstraObs::New( );
+ obs->SetRenderWindow( view.GetWindow( ) );
+ paths->AddObserver( itk::AnyEvent( ), obs );
+ paths->ThrowEventsOn( );
+ }
+ else
+ paths->ThrowEventsOff( );
+ start = std::clock( );
+ paths->Update( );
+ end = std::clock( );
+ seconds = double( end - start ) / double( CLOCKS_PER_SEC );
+ std::cout << "Paths extraction time = " << seconds << std::endl;
+
+ std::cout << "Final seed: " << paths->GetSeed( 0 ) << std::endl;
+
+ // Create polydata
+ vtkSmartPointer< vtkPoints > points =
+ vtkSmartPointer< vtkPoints >::New( );
+ vtkSmartPointer< vtkCellArray > cells =
+ vtkSmartPointer< vtkCellArray >::New( );
+ vtkSmartPointer< vtkFloatArray > scalars =
+ vtkSmartPointer< vtkFloatArray >::New( );
+
+ const TDijkstra::TVertices& endpoints = paths->GetEndPoints( );
+ const TDijkstra::TTree& tree = paths->GetFinalTree( );
+ TDijkstra::TVertices::const_iterator epIt = endpoints.begin( );
+ for( unsigned int epId = 0; epIt != endpoints.end( ); ++epIt, ++epId )
+ {
+ double pd = double( epId ) / double( endpoints.size( ) - 1 );
+
+ TDijkstra::TVertex idx = *epIt;
+ do
+ {
+ TImage::PointType pnt;
+ input_image->TransformIndexToPhysicalPoint( idx, pnt );
+
+ points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
+ scalars->InsertNextTuple1( pd );
+ if( idx != *epIt )
+ {
+ cells->InsertNextCell( 2 );
+ cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
+ cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
+
+ } // fi
+ idx = tree.find( idx )->second;
+
+ } while( idx != tree.find( idx )->second );