]> Creatis software - FrontAlgorithms.git/blobdiff - appli/examples/example_Image_Dijkstra_EndPointDetection.cxx
CMake updated. Some other filters added.
[FrontAlgorithms.git] / appli / examples / example_Image_Dijkstra_EndPointDetection.cxx
index a334de3ae881a372578e7e999cc01d9315695d6f..aea0fe504f203100b5f04ab8123c5af82332133d 100644 (file)
 #include <itkImage.h>
 #include <itkImageToVTKImageFilter.h>
 
-#include <itkImageFileReader.h>
-#include <itkMinimumMaximumImageCalculator.h>
-#include <itkShiftScaleImageFilter.h>
-
 #include <itkSignedMaurerDistanceMapImageFilter.h>
 
-#include <itkImageFileWriter.h>
-
-#include <fpa/Image/DijkstraWithEndPointDetection.h>
 #include <fpa/Base/Functors/InvertCostFunction.h>
-#include <fpa/VTK/ImageMPR.h>
-#include <fpa/VTK/Image3DObserver.h>
-
+#include <fpa/Image/DijkstraWithEndPointDetection.h>
+#include <fpa/VTK/Image2DObserver.h>
 #include <fpa/IO/MinimumSpanningTreeWriter.h>
 #include <fpa/IO/UniqueValuesContainerWriter.h>
 #include <fpa/IO/MatrixValuesContainerWriter.h>
+#include <fpa/VTK/UniqueVerticesToPolyDataFilter.h>
 
-#include <vtkCellArray.h>
-#include <vtkFloatArray.h>
-#include <vtkImageMarchingCubes.h>
-#include <vtkPoints.h>
-#include <vtkPolyData.h>
-#include <vtkSmartPointer.h>
+#include "fpa_Utility.h"
 
 // -------------------------------------------------------------------------
-const unsigned int Dim = 3;
+const unsigned int Dim = 2;
 typedef unsigned char TPixel;
 typedef float         TScalar;
-
-typedef itk::Image< TPixel, Dim >            TImage;
-typedef itk::Image< TScalar, Dim >           TScalarImage;
-typedef itk::ImageToVTKImageFilter< TImage > TVTKInputImage;
-
-// -------------------------------------------------------------------------
-template< class I >
-void ReadImage( typename I::Pointer& image, const std::string& filename );
-
-template< class I >
-void SaveImage( const I* image, const std::string& filename );
-
-template< class I, class O >
-void DistanceMap(
-  const typename I::Pointer& input, typename O::Pointer& output
-  );
+typedef itk::Image< TPixel, Dim >  TImage;
+typedef itk::Image< TScalar, Dim > TScalarImage;
 
 // -------------------------------------------------------------------------
 int main( int argc, char* argv[] )
 {
-  if( argc < 9 )
+  if( argc < 7 )
   {
     std::cerr
       << "Usage: " << argv[ 0 ] << std::endl
-      << " input_image" << std::endl
-      << " output_distancemap" << std::endl
-      << " output_costmap" << std::endl
-      << " output_labels" << std::endl
-      << " output_minimum_spanning_tree" << std::endl
-      << " output_endpoints" << std::endl
-      << " output_bifurcations" << std::endl
-      << " output_branches" << std::endl
+      << "\tinput_image" << std::endl
+      << "\toutput_labels" << std::endl
+      << "\toutput_minimum_spanning_tree" << std::endl
+      << "\toutput_endpoints" << std::endl
+      << "\toutput_bifurcations" << std::endl
+      << "\toutput_branches" << std::endl
       << std::endl;
     return( 1 );
 
   } // fi
   std::string input_image_fn = argv[ 1 ];
-  std::string distancemap_fn = argv[ 2 ];
-  std::string output_costmap_fn = argv[ 3 ];
-  std::string output_labels_fn = argv[ 4 ];
-  std::string mst_output_fn = argv[ 5 ];
-  std::string endpoints_output_fn = argv[ 6 ];
-  std::string bifurcations_output_fn = argv[ 7 ];
-  std::string branches_output_fn = argv[ 8 ];
+  std::string output_labels_fn = argv[ 2 ];
+  std::string output_minimum_spanning_tree_fn = argv[ 3 ];
+  std::string output_endpoints_fn = argv[ 4 ];
+  std::string output_bifurcations_fn = argv[ 5 ];
+  std::string output_branches_fn = argv[ 6 ];
 
   // Read image
   TImage::Pointer input_image;
-  try
+  std::string err = fpa_Utility::ReadImage( input_image, input_image_fn );
+  if( err != "" )
   {
-    ReadImage< TImage >( input_image, input_image_fn );
-  }
-  catch( itk::ExceptionObject& err )
-  {
-    std::cerr
-      << "Error caught while reading \""
-      << input_image_fn << "\": " << err
-      << std::endl;
+    std::cerr << err << std::endl;
     return( 1 );
 
-  } // yrt
+  } // fi
 
-  TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( );
+  // Show image and wait for, at least, one seed
+  itk::ImageToVTKImageFilter< TImage >::Pointer vtk_input_image =
+    itk::ImageToVTKImageFilter< TImage >::New( );
   vtk_input_image->SetInput( input_image );
   vtk_input_image->Update( );
 
-  // Show input image and let some interaction
-  fpa::VTK::ImageMPR view;
-  view.SetBackground( 0.3, 0.2, 0.8 );
-  view.SetSize( 800, 800 );
-  view.SetImage( vtk_input_image->GetOutput( ) );
-
-  vtkSmartPointer< vtkImageMarchingCubes > mc =
-    vtkSmartPointer< vtkImageMarchingCubes >::New( );
-  mc->SetInputData( vtk_input_image->GetOutput( ) );
-  mc->SetValue( 0, 1e-1 );
-  mc->Update( );
-  view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
-
-  // Allow some interaction and wait for, at least, one seed
-  view.Render( );
-  while( view.GetNumberOfSeeds( ) == 0 )
-    view.Start( );
-
-  // Get seed
-  double p[ 3 ];
-  view.GetSeed( 0, p );
-  TImage::PointType pnt;
-  TImage::IndexType seed;
-  pnt[ 0 ] = TImage::PointType::ValueType( p[ 0 ] );
-  pnt[ 1 ] = TImage::PointType::ValueType( p[ 1 ] );
-  pnt[ 2 ] = TImage::PointType::ValueType( p[ 2 ] );
-  input_image->TransformPhysicalPointToIndex( pnt, seed );
+  fpa_Utility::Viewer2DWithSeeds viewer;
+  viewer.SetImage( vtk_input_image->GetOutput( ) );
+  while( viewer.GetNumberOfSeeds( ) == 0 )
+    viewer.Start( );
 
   // Compute squared distance map
-  TScalarImage::Pointer dmap;
-  DistanceMap< TImage, TScalarImage >( input_image, dmap );
+  typename
+    itk::SignedMaurerDistanceMapImageFilter< TImage, TScalarImage >::Pointer
+    dmap =
+    itk::SignedMaurerDistanceMapImageFilter< TImage, TScalarImage >::New( );
+  dmap->SetInput( input_image );
+  dmap->SetBackgroundValue( TPixel( 0 ) );
+  dmap->InsideIsPositiveOn( );
+  dmap->SquaredDistanceOn( );
+  dmap->UseImageSpacingOn( );
+
+  std::time_t start, end;
+  std::time( &start );
+  dmap->Update( );
+  std::time( &end );
+  std::cout
+    << "Distance map time = "
+    << std::difftime( end, start )
+    << "s." << std::endl;
 
   // Prepare cost conversion function
-  typedef fpa::Base::Functors::InvertCostFunction< TScalar > TFunction;
-  TFunction::Pointer function = TFunction::New( );
+  typedef fpa::Base::Functors::InvertCostFunction< TScalar > TCostFunction;
+  TCostFunction::Pointer cost_f = TCostFunction::New( );
 
   // Prepare Dijkstra filter
-  typedef fpa::Image::DijkstraWithEndPointDetection< TScalarImage, TScalarImage > TFilter;
+  typedef fpa::Image::
+    DijkstraWithEndPointDetection< TScalarImage, TScalarImage > TFilter;
   TFilter::Pointer filter = TFilter::New( );
-  filter->SetInput( dmap );
-  filter->SetConversionFunction( function );
+  filter->SetInput( dmap->GetOutput( ) );
+  filter->SetConversionFunction( cost_f );
   filter->SetNeighborhoodOrder( 2 );
   filter->StopAtOneFrontOff( );
-  filter->AddSeed( seed, TScalar( 0 ) );
+
+  // Associate seed
+  TImage::PointType pnt;
+  TImage::IndexType idx;
+  viewer.GetSeed( pnt, 0 );
+  if( input_image->TransformPhysicalPointToIndex( pnt, idx ) )
+    filter->AddSeed( idx, 0 );
 
   // Prepare graphical debugger
-  typedef fpa::VTK::Image3DObserver< TFilter, vtkRenderWindow > TDebugger;
+  typedef fpa::VTK::Image2DObserver< TFilter, vtkRenderWindow > TDebugger;
   TDebugger::Pointer debugger = TDebugger::New( );
-  debugger->SetRenderWindow( view.GetWindow( ) );
-  debugger->SetRenderPercentage( 0.0001 );
+  debugger->SetRenderWindow( viewer.Window );
+  debugger->SetRenderPercentage( 0.01 );
   filter->AddObserver( itk::AnyEvent( ), debugger );
   filter->ThrowEventsOn( );
 
   // Go!
-  std::time_t start, end;
   std::time( &start );
   filter->Update( );
   std::time( &end );
   std::cout
     << "Extraction time = "
     << std::difftime( end, start )
-    << " s." << std::endl;
-
-  // Outputs
-  const TFilter::TOutputImage* accumulated_costs = filter->GetOutput( );
-  const TFilter::TLabelImage* labeled_image = filter->GetLabelImage( );
-  const TFilter::TMinimumSpanningTree* mst = filter->GetMinimumSpanningTree( );
-  const TFilter::TUniqueVertices* endpoints = filter->GetEndPoints( );
-  const TFilter::TUniqueVertices* bifurcations = filter->GetBifurcations( );
-  const TFilter::TBranches* branches = filter->GetBranches( );
-  unsigned long nBranches = filter->GetNumberOfBranches( );
-
-  // Save outputs
-  SaveImage( dmap.GetPointer( ), distancemap_fn );
-  SaveImage( accumulated_costs, output_costmap_fn );
-  SaveImage( labeled_image, output_labels_fn );
+    << "s." << std::endl;
+
+  // Create new actors
+  typedef vtkSmartPointer< fpa::VTK::UniqueVerticesToPolyDataFilter< TFilter::TUniqueVertices, TFilter::TInputImage > > TVertices2PD;
+
+  TVertices2PD endpoints2pd = TVertices2PD::New( );
+  endpoints2pd->SetInput( filter->GetEndPoints( ) );
+  endpoints2pd->SetImage( filter->GetInput( ) );
+  endpoints2pd->Update( );
+
+  vtkSmartPointer< vtkPolyDataMapper > endpoints_mapper =
+    vtkSmartPointer< vtkPolyDataMapper >::New( );
+  endpoints_mapper->SetInputConnection( endpoints2pd->GetOutputPort( ) );
+
+  vtkSmartPointer< vtkActor > endpoints_actor =
+    vtkSmartPointer< vtkActor >::New( );
+  endpoints_actor->SetMapper( endpoints_mapper );
+  endpoints_actor->GetProperty( )->SetColor( 0, 1, 0 );
+  endpoints_actor->GetProperty( )->SetPointSize( 5 );
+  viewer.Renderer->AddActor( endpoints_actor );
+
+  TVertices2PD bifurcations2pd = TVertices2PD::New( );
+  bifurcations2pd->SetInput( filter->GetBifurcations( ) );
+  bifurcations2pd->SetImage( filter->GetInput( ) );
+  bifurcations2pd->Update( );
+
+  vtkSmartPointer< vtkPolyDataMapper > bifurcations_mapper =
+    vtkSmartPointer< vtkPolyDataMapper >::New( );
+  bifurcations_mapper->SetInputConnection( bifurcations2pd->GetOutputPort( ) );
+
+  vtkSmartPointer< vtkActor > bifurcations_actor =
+    vtkSmartPointer< vtkActor >::New( );
+  bifurcations_actor->SetMapper( bifurcations_mapper );
+  bifurcations_actor->GetProperty( )->SetColor( 1, 0, 0 );
+  bifurcations_actor->GetProperty( )->SetPointSize( 5 );
+  viewer.Renderer->AddActor( bifurcations_actor );
+
+  // Some more interaction and finish
+  viewer.Render( );
+  viewer.Start( );
+
+  // Save results
+  err = fpa_Utility::SaveImage( filter->GetLabelImage( ), output_labels_fn );
+  if( err != "" )
+    std::cerr << err << std::endl;
 
   fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::Pointer
     mst_writer =
     fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::New( );
-  mst_writer->SetInput( mst );
-  mst_writer->SetFileName( mst_output_fn );
+  mst_writer->SetInput( filter->GetMinimumSpanningTree( ) );
+  mst_writer->SetFileName( output_minimum_spanning_tree_fn );
   mst_writer->Update( );
 
   fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer
     endpoints_writer =
     fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( );
-  endpoints_writer->SetInput( endpoints );
-  endpoints_writer->SetFileName( endpoints_output_fn );
+  endpoints_writer->SetInput( filter->GetEndPoints( ) );
+  endpoints_writer->SetFileName( output_endpoints_fn );
   endpoints_writer->Update( );
 
   fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer
     bifurcations_writer =
     fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( );
-  bifurcations_writer->SetInput( bifurcations );
-  bifurcations_writer->SetFileName( bifurcations_output_fn );
+  bifurcations_writer->SetInput( filter->GetBifurcations( ) );
+  bifurcations_writer->SetFileName( output_bifurcations_fn );
   bifurcations_writer->Update( );
 
   fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::Pointer
     branches_writer =
     fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::New( );
-  branches_writer->SetInput( branches );
-  branches_writer->SetFileName( branches_output_fn );
+  branches_writer->SetInput( filter->GetBranches( ) );
+  branches_writer->SetFileName( output_branches_fn );
   branches_writer->Update( );
 
-  // Show endpoints
-  vtkSmartPointer< vtkPoints > endpoints_points =
-    vtkSmartPointer< vtkPoints >::New( );
-  vtkSmartPointer< vtkCellArray > endpoints_cells =
-    vtkSmartPointer< vtkCellArray >::New( );
-  for(
-    TFilter::TUniqueVertices::ConstIterator epIt = endpoints->Begin( );
-    epIt != endpoints->End( );
-    ++epIt
-    )
-  {
-    TImage::PointType pnt;
-    input_image->TransformIndexToPhysicalPoint( *epIt, pnt );
-    endpoints_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
-    endpoints_cells->InsertNextCell( 1 );
-    endpoints_cells->
-      InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 1 );
-
-  } // rof
-  vtkSmartPointer< vtkPolyData > endpoints_polydata =
-    vtkSmartPointer< vtkPolyData >::New( );
-  endpoints_polydata->SetPoints( endpoints_points );
-  endpoints_polydata->SetVerts( endpoints_cells );
-  view.AddPolyData( endpoints_polydata, 0, 0, 1, 1 );
-
-  // Show bifurcations
-  vtkSmartPointer< vtkPoints > bifurcations_points =
-    vtkSmartPointer< vtkPoints >::New( );
-  vtkSmartPointer< vtkCellArray > bifurcations_cells =
-    vtkSmartPointer< vtkCellArray >::New( );
-  for(
-    TFilter::TUniqueVertices::ConstIterator bfIt = bifurcations->Begin( );
-    bfIt != bifurcations->End( );
-    ++bfIt
-    )
-  {
-    TImage::PointType pnt;
-    input_image->TransformIndexToPhysicalPoint( *bfIt, pnt );
-    bifurcations_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
-    bifurcations_cells->InsertNextCell( 1 );
-    bifurcations_cells->
-      InsertCellPoint( bifurcations_points->GetNumberOfPoints( ) - 1 );
-
-  } // rof
-  vtkSmartPointer< vtkPolyData > bifurcations_polydata =
-    vtkSmartPointer< vtkPolyData >::New( );
-  bifurcations_polydata->SetPoints( bifurcations_points );
-  bifurcations_polydata->SetVerts( bifurcations_cells );
-  view.AddPolyData( bifurcations_polydata, 0, 1, 0, 1 );
-
-  // Show branches (simple and detailed)
-  vtkSmartPointer< vtkPoints > simple_branches_points =
-    vtkSmartPointer< vtkPoints >::New( );
-  vtkSmartPointer< vtkCellArray > simple_branches_cells =
-    vtkSmartPointer< vtkCellArray >::New( );
-
-  vtkSmartPointer< vtkPoints > detailed_branches_points =
-    vtkSmartPointer< vtkPoints >::New( );
-  vtkSmartPointer< vtkCellArray > detailed_branches_cells =
-    vtkSmartPointer< vtkCellArray >::New( );
-  vtkSmartPointer< vtkFloatArray > detailed_branches_scalars =
-    vtkSmartPointer< vtkFloatArray >::New( );
-
-  TFilter::TBranches::ConstIterator brIt = branches->Begin( );
-  for( ; brIt != branches->End( ); ++brIt )
-  {
-    // Branch's first point
-    TImage::PointType first_point;
-    input_image->TransformIndexToPhysicalPoint( brIt->first, first_point );
-    unsigned long first_id = simple_branches_points->GetNumberOfPoints( );
-    simple_branches_points->InsertNextPoint(
-      first_point[ 0 ], first_point[ 1 ], first_point[ 2 ]
-      );
-
-    TFilter::TBranches::ConstRowIterator brRowIt = branches->Begin( brIt );
-    for( ; brRowIt != branches->End( brIt ); ++brRowIt )
-    {
-      // Branch's second point
-      TImage::PointType second_point;
-      input_image->
-        TransformIndexToPhysicalPoint( brRowIt->first, second_point );
-      unsigned long second_id = simple_branches_points->GetNumberOfPoints( );
-      simple_branches_points->InsertNextPoint(
-        second_point[ 0 ], second_point[ 1 ], second_point[ 2 ]
-        );
-      simple_branches_cells->InsertNextCell( 2 );
-      simple_branches_cells->InsertCellPoint( first_id );
-      simple_branches_cells->InsertCellPoint( second_id );
-
-      // Detailed path
-      double pathId = double( brRowIt->second - 1 ) / double( nBranches - 1 );
-      TFilter::TVertices path;
-      mst->GetPath( path, brIt->first, brRowIt->first );
-      TFilter::TVertices::const_iterator pIt = path.begin( );
-      for( ; pIt != path.end( ); ++pIt )
-      {
-        TImage::PointType path_point;
-        input_image->TransformIndexToPhysicalPoint( *pIt, path_point );
-        detailed_branches_points->InsertNextPoint(
-          path_point[ 0 ], path_point[ 1 ], path_point[ 2 ]
-          );
-        detailed_branches_scalars->InsertNextTuple1( pathId );
-        if( pIt != path.begin( ) )
-        {
-          unsigned long nPoints =
-            detailed_branches_points->GetNumberOfPoints( );
-          detailed_branches_cells->InsertNextCell( 2 );
-          detailed_branches_cells->InsertCellPoint( nPoints - 2 );
-          detailed_branches_cells->InsertCellPoint( nPoints - 1 );
-
-        } // fi
-
-      } // rof
-
-    } // rof
-
-  } // rof
-  vtkSmartPointer< vtkPolyData > simple_branches_polydata =
-    vtkSmartPointer< vtkPolyData >::New( );
-  simple_branches_polydata->SetPoints( simple_branches_points );
-  simple_branches_polydata->SetLines( simple_branches_cells );
-  view.AddPolyData( simple_branches_polydata, 1, 0, 1, 1 );
-
-  vtkSmartPointer< vtkPolyData > detailed_branches_polydata =
-    vtkSmartPointer< vtkPolyData >::New( );
-  detailed_branches_polydata->SetPoints( detailed_branches_points );
-  detailed_branches_polydata->SetLines( detailed_branches_cells );
-  detailed_branches_polydata->
-    GetPointData( )->SetScalars( detailed_branches_scalars );
-  view.AddPolyData( detailed_branches_polydata, 1 );
-
-  // Let some more interaction
-  view.Render( );
-  view.Start( );
-
   return( 0 );
 }
 
-// -------------------------------------------------------------------------
-template< class I >
-void ReadImage( typename I::Pointer& image, const std::string& filename )
-{
-  typename itk::ImageFileReader< I >::Pointer reader =
-    itk::ImageFileReader< I >::New( );
-  reader->SetFileName( filename );
-  reader->Update( );
-
-  typename itk::MinimumMaximumImageCalculator< I >::Pointer minmax =
-    itk::MinimumMaximumImageCalculator< I >::New( );
-  minmax->SetImage( reader->GetOutput( ) );
-  minmax->Compute( );
-  double vmin = double( minmax->GetMinimum( ) );
-  double vmax = double( minmax->GetMaximum( ) );
-
-  typename itk::ShiftScaleImageFilter< I, I >::Pointer shift =
-    itk::ShiftScaleImageFilter< I, I >::New( );
-  shift->SetInput( reader->GetOutput( ) );
-  shift->SetScale( vmax - vmin );
-  shift->SetShift( vmin / ( vmax - vmin ) );
-  shift->Update( );
-
-  image = shift->GetOutput( );
-  image->DisconnectPipeline( );
-}
-
-// -------------------------------------------------------------------------
-template< class I >
-void SaveImage( const I* image, const std::string& filename )
-{
-  typename itk::ImageFileWriter< I >::Pointer writer =
-    itk::ImageFileWriter< I >::New( );
-  writer->SetInput( image );
-  writer->SetFileName( filename );
-  try
-  {
-    writer->Update( );
-  }
-  catch( itk::ExceptionObject& err )
-  {
-    std::cerr
-      << "Error saving \"" << filename << "\": " << err
-      << std::endl;
-
-  } // yrt
-}
-
-// -------------------------------------------------------------------------
-template< class I, class O >
-void DistanceMap(
-  const typename I::Pointer& input, typename O::Pointer& output
-  )
-{
-  typename itk::SignedMaurerDistanceMapImageFilter< I, O >::Pointer dmap =
-    itk::SignedMaurerDistanceMapImageFilter< I, O >::New( );
-  dmap->SetInput( input );
-  dmap->SetBackgroundValue( ( typename I::PixelType )( 0 ) );
-  dmap->InsideIsPositiveOn( );
-  dmap->SquaredDistanceOn( );
-  dmap->UseImageSpacingOn( );
-
-  std::time_t start, end;
-  std::time( &start );
-  dmap->Update( );
-  std::time( &end );
-  std::cout
-    << "Distance map time = "
-    << std::difftime( end, start )
-    << " s." << std::endl;
+// eof - $RCSfile$
 
-  output = dmap->GetOutput( );
-  output->DisconnectPipeline( );
-}
+/* TODO
+   vtkSmartPointer< vtkPoints > simple_branches_points =
+   vtkSmartPointer< vtkPoints >::New( );
+   vtkSmartPointer< vtkCellArray > simple_branches_cells =
+   vtkSmartPointer< vtkCellArray >::New( );
+
+   vtkSmartPointer< vtkPoints > detailed_branches_points =
+   vtkSmartPointer< vtkPoints >::New( );
+   vtkSmartPointer< vtkCellArray > detailed_branches_cells =
+   vtkSmartPointer< vtkCellArray >::New( );
+   vtkSmartPointer< vtkFloatArray > detailed_branches_scalars =
+   vtkSmartPointer< vtkFloatArray >::New( );
+
+   TFilter::TBranches::ConstIterator brIt = branches->Begin( );
+   for( ; brIt != branches->End( ); ++brIt )
+   {
+   // Branch's first point
+   TImage::PointType first_point;
+   input_image->TransformIndexToPhysicalPoint( brIt->first, first_point );
+   unsigned long first_id = simple_branches_points->GetNumberOfPoints( );
+   simple_branches_points->InsertNextPoint(
+   first_point[ 0 ], first_point[ 1 ], first_point[ 2 ]
+   );
+
+   TFilter::TBranches::ConstRowIterator brRowIt = branches->Begin( brIt );
+   for( ; brRowIt != branches->End( brIt ); ++brRowIt )
+   {
+   // Branch's second point
+   TImage::PointType second_point;
+   input_image->
+   TransformIndexToPhysicalPoint( brRowIt->first, second_point );
+   unsigned long second_id = simple_branches_points->GetNumberOfPoints( );
+   simple_branches_points->InsertNextPoint(
+   second_point[ 0 ], second_point[ 1 ], second_point[ 2 ]
+   );
+   simple_branches_cells->InsertNextCell( 2 );
+   simple_branches_cells->InsertCellPoint( first_id );
+   simple_branches_cells->InsertCellPoint( second_id );
+
+   // Detailed path
+   double pathId = double( brRowIt->second - 1 ) / double( nBranches - 1 );
+   TFilter::TVertices path;
+   mst->GetPath( path, brIt->first, brRowIt->first );
+   TFilter::TVertices::const_iterator pIt = path.begin( );
+   for( ; pIt != path.end( ); ++pIt )
+   {
+   TImage::PointType path_point;
+   input_image->TransformIndexToPhysicalPoint( *pIt, path_point );
+   detailed_branches_points->InsertNextPoint(
+   path_point[ 0 ], path_point[ 1 ], path_point[ 2 ]
+   );
+   detailed_branches_scalars->InsertNextTuple1( pathId );
+   if( pIt != path.begin( ) )
+   {
+   unsigned long nPoints =
+   detailed_branches_points->GetNumberOfPoints( );
+   detailed_branches_cells->InsertNextCell( 2 );
+   detailed_branches_cells->InsertCellPoint( nPoints - 2 );
+   detailed_branches_cells->InsertCellPoint( nPoints - 1 );
+
+   } // fi
+
+   } // rof
+
+   } // rof
+
+   } // rof
+   vtkSmartPointer< vtkPolyData > simple_branches_polydata =
+   vtkSmartPointer< vtkPolyData >::New( );
+   simple_branches_polydata->SetPoints( simple_branches_points );
+   simple_branches_polydata->SetLines( simple_branches_cells );
+   view.AddPolyData( simple_branches_polydata, 1, 0, 1, 1 );
+
+   vtkSmartPointer< vtkPolyData > detailed_branches_polydata =
+   vtkSmartPointer< vtkPolyData >::New( );
+   detailed_branches_polydata->SetPoints( detailed_branches_points );
+   detailed_branches_polydata->SetLines( detailed_branches_cells );
+   detailed_branches_polydata->
+   GetPointData( )->SetScalars( detailed_branches_scalars );
+   view.AddPolyData( detailed_branches_polydata, 1 );
+
+   // Let some more interaction
+   view.Render( );
+   view.Start( );
+
+   return( 0 );
+   }
+*/
 
 // eof - $RCSfile$