A dijkstra example added

[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 156469ef0aeb414f70ec6803503a67448d3f22aa..f2b82794894656dd0aa8a3158e47257cf5efb49e 100644 (file)
--- a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx
+++ b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx
@@ -2,62 +2,45 @@
 #include <iostream>
 #include <string>
 
+#include <itkConstNeighborhoodIterator.h>
 #include <itkDanielssonDistanceMapImageFilter.h>
 #include <itkImage.h>
 #include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
 #include <itkImageToVTKImageFilter.h>
 
-#include <fpa/Base/Functors/InvertCostFunction.h>
-#include <fpa/Image/Dijkstra.h>
 #include <fpa/Image/RegionGrowWithMultipleThresholds.h>
 #include <fpa/VTK/ImageMPR.h>
 #include <fpa/VTK/Image3DObserver.h>
 
-/*
-  #include <limits>
-  #include <set>
-
-  #include <itkImageFileWriter.h>
-
-
-  #include <vtkImageActor.h>
-  #include <vtkImageMarchingCubes.h>
-  #include <vtkProperty.h>
-  #include <vtkRenderer.h>
-  #include <vtkRenderWindow.h>
-  #include <vtkRenderWindowInteractor.h>
-  #include <vtkSmartPointer.h>
-  #include <vtkSphereSource.h>
-
-*/
-
 // -------------------------------------------------------------------------
 const unsigned int Dim = 3;
 typedef short                                TPixel;
 typedef double                               TScalar;
 typedef itk::Image< TPixel, Dim >            TImage;
 typedef itk::Image< TScalar, Dim >           TDistanceMap;
+typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
+
+typedef itk::ImageFileReader< TImage >                       TImageReader;
+typedef itk::ImageFileWriter< TImage >                       TImageWriter;
+typedef itk::ImageFileWriter< TDistanceMap >           TDistanceMapWriter;
+typedef fpa::Image::RegionGrowWithMultipleThresholds< TImage > TSegmentor;
+typedef
+itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >
+TDistanceFilter;
 
-/*
-  typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
-  typedef itk::ImageFileReader< TImage >   TImageReader;
-  typedef itk::ImageFileWriter< TImage >   TImageWriter;
-  typedef
-  fpa::Image::RegionGrowWithMultipleThresholds< TImage >
-  TSegmentor;
-  typedef
-  
-  TObserver;
-*/
+typedef fpa::VTK::ImageMPR TMPR;
+typedef fpa::VTK::Image3DObserver< TSegmentor, vtkRenderWindow > TSegmentorObs;
 
 // -------------------------------------------------------------------------
 int main( int argc, char* argv[] )
 {
-  if( argc < 6 )
+  if( argc < 7 )
   {
     std::cerr
       << "Usage: " << argv[ 0 ]
-      << " input_image thr_0 thr_1 step output_image"
+      << " input_image thr_0 thr_1 step"
+      << " output_segmentation output_distancemap"
       << " visual_debug"
       << std::endl;
     return( 1 );
@@ -67,14 +50,14 @@ int main( int argc, char* argv[] )
   TPixel thr_0 = TPixel( std::atof( argv[ 2 ] ) );
   TPixel thr_1 = TPixel( std::atof( argv[ 3 ] ) );
   unsigned int step = std::atoi( argv[ 4 ] );
-  std::string output_image_fn = argv[ 5 ];
+  std::string output_segmentation_fn = argv[ 5 ];
+  std::string output_distancemap_fn = argv[ 6 ];
   bool visual_debug = false;
-  if( argc > 6 )
-    visual_debug = ( std::atoi( argv[ 6 ] ) == 1 );
+  if( argc > 7 )
+    visual_debug = ( std::atoi( argv[ 7 ] ) == 1 );
 
   // Read image
-  itk::ImageFileReader< TImage >::Pointer input_image_reader =
-    itk::ImageFileReader< TImage >::New( );
+  TImageReader::Pointer input_image_reader = TImageReader::New( );
   input_image_reader->SetFileName( input_image_fn );
   try
   {
@@ -89,12 +72,11 @@ int main( int argc, char* argv[] )
   TImage::ConstPointer input_image = input_image_reader->GetOutput( );
 
   // Show image
-  itk::ImageToVTKImageFilter< TImage >::Pointer vtk_image =
-    itk::ImageToVTKImageFilter< TImage >::New( );
+  TVTKImage::Pointer vtk_image = TVTKImage::New( );
   vtk_image->SetInput( input_image );
   vtk_image->Update( );
 
-  fpa::VTK::ImageMPR view;
+  TMPR view;
   view.SetBackground( 0.3, 0.2, 0.8 );
   view.SetSize( 800, 800 );
   view.SetImage( vtk_image->GetOutput( ) );
@@ -114,8 +96,7 @@ int main( int argc, char* argv[] )
   input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx );
 
   // Segment input image
-  fpa::Image::RegionGrowWithMultipleThresholds< TImage >::Pointer segmentor =
-    fpa::Image::RegionGrowWithMultipleThresholds< TImage >::New( );
+  TSegmentor::Pointer segmentor = TSegmentor::New( );
   segmentor->AddThresholds( thr_0, thr_1, step );
   segmentor->AddSeed( seed_idx, 0 );
   segmentor->SetInput( input_image );
@@ -126,8 +107,7 @@ int main( int argc, char* argv[] )
   if( visual_debug )
   {
     // Configure observer
-    fpa::VTK::Image3DObserver< fpa::Image::RegionGrowWithMultipleThresholds< TImage >, vtkRenderWindow >::Pointer obs =
-      fpa::VTK::Image3DObserver< fpa::Image::RegionGrowWithMultipleThresholds< TImage >, vtkRenderWindow >::New( );
+    TSegmentorObs::Pointer obs = TSegmentorObs::New( );
     obs->SetImage( input_image, view.GetWindow( ) );
     segmentor->AddObserver( itk::AnyEvent( ), obs );
     segmentor->ThrowEventsOn( );
@@ -141,9 +121,7 @@ int main( int argc, char* argv[] )
   std::cout << "Segmentation time = " << seconds << std::endl;
 
   // Extract distance map
-  itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >::Pointer
-    distanceMap =
-    itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >::New( );
+  TDistanceFilter::Pointer distanceMap = TDistanceFilter::New( );
   distanceMap->SetInput( segmentor->GetOutput( ) );
   distanceMap->InputIsBinaryOn( );
   start = std::clock( );
@@ -152,35 +130,19 @@ int main( int argc, char* argv[] )
   seconds = double( end - start ) / double( CLOCKS_PER_SEC );
   std::cout << "Distance map time = " << seconds << std::endl;
 
-  // Extract paths
-  fpa::Image::Dijkstra< TDistanceMap, TScalar >::Pointer paths =
-    fpa::Image::Dijkstra< TDistanceMap, TScalar >::New( );
-  paths->AddSeed( segmentor->GetSeed( 0 ), TScalar( 0 ) );
-  paths->SetInput( distanceMap->GetOutput( ) );
-  paths->SetNeighborhoodOrder( 1 );
+  std::cout << "Final seed: " << seed_idx << std::endl;
 
-  // Associate an inversion cost function
-  fpa::Base::Functors::InvertCostFunction< TScalar >::Pointer
-    cost_function =
-    fpa::Base::Functors::InvertCostFunction< TScalar >::New( );
-  paths->SetCostConversion( cost_function );
+  TDistanceMapWriter::Pointer distancemap_writer =
+    TDistanceMapWriter::New( );
+  distancemap_writer->SetInput( distanceMap->GetOutput( ) );
+  distancemap_writer->SetFileName( output_distancemap_fn );
+  distancemap_writer->Update( );
 
-  if( visual_debug )
-  {
-    // Configure observer
-    fpa::VTK::Image3DObserver< fpa::Image::Dijkstra< TDistanceMap, TScalar >, vtkRenderWindow >::Pointer obs =
-      fpa::VTK::Image3DObserver< fpa::Image::Dijkstra< TDistanceMap, TScalar >, vtkRenderWindow >::New( );
-    obs->SetImage( distanceMap->GetOutput( ), 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;
+  TImageWriter::Pointer segmentation_writer =
+    TImageWriter::New( );
+  segmentation_writer->SetInput( segmentor->GetOutput( ) );
+  segmentation_writer->SetFileName( output_segmentation_fn );
+  segmentation_writer->Update( );
 
   // Show result
   /*