Skeletonization debugged

[FrontAlgorithms.git] / appli / examples / example_ImageAlgorithm_Skeletonization.cxx

diff --git a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx
index e0121bd6ba9c173d41a532a56e6d0cb0c507e1b8..7bc80a95896faac60a6972ecde2fac5f2f62b228 100644 (file)
--- a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx
+++ b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx
@@ -9,7 +9,10 @@
 #include <itkImageFileWriter.h>
 #include <itkImageToVTKImageFilter.h>
 
+#include <vtkImageMarchingCubes.h>
+
 #include <fpa/Image/RegionGrowWithMultipleThresholds.h>
+#include <fpa/Image/DijkstraWithSphereBacktracking.h>
 #include <fpa/VTK/ImageMPR.h>
 #include <fpa/VTK/Image3DObserver.h>
 
@@ -20,6 +23,7 @@ typedef double                               TScalar;
 typedef itk::Image< TPixel, Dim >            TImage;
 typedef itk::Image< TScalar, Dim >           TDistanceMap;
 typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
+typedef itk::ImageToVTKImageFilter< TDistanceMap > TVTKDistanceMap;
 
 typedef itk::ImageFileReader< TImage >                       TImageReader;
 typedef itk::ImageFileWriter< TImage >                       TImageWriter;
@@ -28,9 +32,13 @@ typedef fpa::Image::RegionGrowWithMultipleThresholds< TImage > TSegmentor;
 typedef
 itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >
 TDistanceFilter;
+typedef
+fpa::Image::DijkstraWithSphereBacktracking< TDistanceMap, TScalar >
+TDijkstra;
 
 typedef fpa::VTK::ImageMPR TMPR;
 typedef fpa::VTK::Image3DObserver< TSegmentor, vtkRenderWindow > TSegmentorObs;
+typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs;
 
 // -------------------------------------------------------------------------
 int main( int argc, char* argv[] )
@@ -131,19 +139,101 @@ int main( int argc, char* argv[] )
   seconds = double( end - start ) / double( CLOCKS_PER_SEC );
   std::cout << "Distance map time = " << seconds << std::endl;
 
-  std::cout << "Final seed: " << seed_idx << std::endl;
+  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 );
+
+  } // rof
 
-  TDistanceMapWriter::Pointer distancemap_writer =
-    TDistanceMapWriter::New( );
-  distancemap_writer->SetInput( distanceMap->GetOutput( ) );
-  distancemap_writer->SetFileName( output_distancemap_fn );
-  distancemap_writer->Update( );
+  vtkSmartPointer< vtkPolyData > vtk_tree =
+    vtkSmartPointer< vtkPolyData >::New( );
+  vtk_tree->SetPoints( points );
+  vtk_tree->SetLines( cells );
+  vtk_tree->GetPointData( )->SetScalars( scalars );
 
-  TImageWriter::Pointer segmentation_writer =
-    TImageWriter::New( );
-  segmentation_writer->SetInput( segmentor->GetOutput( ) );
-  segmentation_writer->SetFileName( output_segmentation_fn );
-  segmentation_writer->Update( );
+  view.AddPolyData( vtk_tree );
+  view.Render( );
+  view.Start( );
+
+  /* TODO
+     TDistanceMapWriter::Pointer distancemap_writer =
+     TDistanceMapWriter::New( );
+     distancemap_writer->SetInput( distanceMap->GetOutput( ) );
+     distancemap_writer->SetFileName( output_distancemap_fn );
+     distancemap_writer->Update( );
+
+     TImageWriter::Pointer segmentation_writer =
+     TImageWriter::New( );
+     segmentation_writer->SetInput( segmentor->GetOutput( ) );
+     segmentation_writer->SetFileName( output_segmentation_fn );
+     segmentation_writer->Update( );
+  */
 
   // Show result
   /*