appli/examples/example_ImageAlgorithm_Skeletonization.cxx

   1 #include <ctime>
   2 #include <iostream>
   3 #include <string>
   4
   5 #include <itkConstNeighborhoodIterator.h>
   6 #include <itkDanielssonDistanceMapImageFilter.h>
   7 #include <itkImage.h>
   8 #include <itkImageFileReader.h>
   9 #include <itkImageFileWriter.h>
  10 #include <itkImageToVTKImageFilter.h>
  11
  12 #include <vtkImageMarchingCubes.h>
  13
  14 #include <fpa/Image/RegionGrowWithMultipleThresholds.h>
  15 #include <fpa/Image/DijkstraWithSphereBacktracking.h>
  16 #include <fpa/VTK/ImageMPR.h>
  17 #include <fpa/VTK/Image3DObserver.h>
  18
  19 // -------------------------------------------------------------------------
  20 const unsigned int Dim = 3;
  21 typedef short                                TPixel;
  22 typedef double                               TScalar;
  23 typedef itk::Image< TPixel, Dim >            TImage;
  24 typedef itk::Image< TScalar, Dim >           TDistanceMap;
  25 typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
  26 typedef itk::ImageToVTKImageFilter< TDistanceMap > TVTKDistanceMap;
  27
  28 typedef itk::ImageFileReader< TImage >                       TImageReader;
  29 typedef itk::ImageFileWriter< TImage >                       TImageWriter;
  30 typedef itk::ImageFileWriter< TDistanceMap >           TDistanceMapWriter;
  31 typedef fpa::Image::RegionGrowWithMultipleThresholds< TImage > TSegmentor;
  32 typedef
  33 itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >
  34 TDistanceFilter;
  35 typedef
  36 fpa::Image::DijkstraWithSphereBacktracking< TDistanceMap, TScalar >
  37 TDijkstra;
  38
  39 typedef fpa::VTK::ImageMPR TMPR;
  40 typedef fpa::VTK::Image3DObserver< TSegmentor, vtkRenderWindow > TSegmentorObs;
  41 typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs;
  42
  43 // -------------------------------------------------------------------------
  44 int main( int argc, char* argv[] )
  45 {
  46   if( argc < 7 )
  47   {
  48     std::cerr
  49       << "Usage: " << argv[ 0 ]
  50       << " input_image thr_0 thr_1 step"
  51       << " output_segmentation output_distancemap"
  52       << " visual_debug"
  53       << std::endl;
  54     return( 1 );
  55
  56   } // fi
  57   std::string input_image_fn = argv[ 1 ];
  58   TPixel thr_0 = TPixel( std::atof( argv[ 2 ] ) );
  59   TPixel thr_1 = TPixel( std::atof( argv[ 3 ] ) );
  60   unsigned int step = std::atoi( argv[ 4 ] );
  61   std::string output_segmentation_fn = argv[ 5 ];
  62   std::string output_distancemap_fn = argv[ 6 ];
  63   bool visual_debug = false;
  64   if( argc > 7 )
  65     visual_debug = ( std::atoi( argv[ 7 ] ) == 1 );
  66
  67   // Read image
  68   TImageReader::Pointer input_image_reader = TImageReader::New( );
  69   input_image_reader->SetFileName( input_image_fn );
  70   try
  71   {
  72     input_image_reader->Update( );
  73   }
  74   catch( itk::ExceptionObject& err )
  75   {
  76     std::cerr << "Error caught: " << err << std::endl;
  77     return( 1 );
  78
  79   } // yrt
  80   TImage::ConstPointer input_image = input_image_reader->GetOutput( );
  81
  82   // Show image
  83   TVTKImage::Pointer vtk_image = TVTKImage::New( );
  84   vtk_image->SetInput( input_image );
  85   vtk_image->Update( );
  86
  87   TMPR view;
  88   view.SetBackground( 0.3, 0.2, 0.8 );
  89   view.SetSize( 800, 800 );
  90   view.SetImage( vtk_image->GetOutput( ) );
  91
  92   // Wait for a seed to be given
  93   while( view.GetNumberOfSeeds( ) == 0 )
  94     view.Start( );
  95
  96   // Get seed from user
  97   double seed[ 3 ];
  98   view.GetSeed( 0, seed );
  99   TImage::PointType seed_pnt;
 100   seed_pnt[ 0 ] = seed[ 0 ];
 101   seed_pnt[ 1 ] = seed[ 1 ];
 102   seed_pnt[ 2 ] = seed[ 2 ];
 103   TImage::IndexType seed_idx;
 104   input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx );
 105
 106   // Segment input image
 107   TSegmentor::Pointer segmentor = TSegmentor::New( );
 108   segmentor->AddThresholds( thr_0, thr_1, step );
 109   segmentor->AddSeed( seed_idx, 0 );
 110   segmentor->SetInput( input_image );
 111   segmentor->SetNeighborhoodOrder( 1 );
 112   segmentor->SetDifferenceThreshold( double( 3 ) );
 113   segmentor->SetInsideValue( TPixel( 0 ) );  // WARNING: This is to optimize
 114   segmentor->SetOutsideValue( TPixel( 1 ) ); // distance map calculation
 115   if( visual_debug )
 116   {
 117     // Configure observer
 118     TSegmentorObs::Pointer obs = TSegmentorObs::New( );
 119     obs->SetRenderWindow( view.GetWindow( ) );
 120     segmentor->AddObserver( itk::AnyEvent( ), obs );
 121     segmentor->ThrowEventsOn( );
 122   }
 123   else
 124     segmentor->ThrowEventsOff( );
 125   std::clock_t start = std::clock( );
 126   segmentor->Update( );
 127   std::clock_t end = std::clock( );
 128   double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
 129   std::cout << "Segmentation time = " << seconds << std::endl;
 130
 131   // Extract distance map
 132   TDistanceFilter::Pointer distanceMap = TDistanceFilter::New( );
 133   distanceMap->SetInput( segmentor->GetOutput( ) );
 134   distanceMap->InputIsBinaryOn( );
 135   distanceMap->SquaredDistanceOn( );
 136   start = std::clock( );
 137   distanceMap->Update( );
 138   end = std::clock( );
 139   seconds = double( end - start ) / double( CLOCKS_PER_SEC );
 140   std::cout << "Distance map time = " << seconds << std::endl;
 141
 142   TVTKDistanceMap::Pointer vtk_distanceMap = TVTKDistanceMap::New( );
 143   vtk_distanceMap->SetInput( distanceMap->GetOutput( ) );
 144   vtk_distanceMap->Update( );
 145
 146   vtkSmartPointer< vtkImageMarchingCubes > mc =
 147     vtkSmartPointer< vtkImageMarchingCubes >::New( );
 148   mc->SetInputData( vtk_distanceMap->GetOutput( ) );
 149   mc->SetValue( 0, 1e-2 );
 150   mc->Update( );
 151   view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
 152   view.Render( );
 153
 154   // Extract paths
 155   TDijkstra::Pointer paths = TDijkstra::New( );
 156   paths->AddSeed( segmentor->GetSeed( 0 ), TScalar( 0 ) );
 157   paths->SetInput( distanceMap->GetOutput( ) );
 158   paths->SetNeighborhoodOrder( 1 );
 159
 160   if( visual_debug )
 161   {
 162     // Configure observer
 163     TDijkstraObs::Pointer obs = TDijkstraObs::New( );
 164     obs->SetRenderWindow( view.GetWindow( ) );
 165     paths->AddObserver( itk::AnyEvent( ), obs );
 166     paths->ThrowEventsOn( );
 167   }
 168   else
 169     paths->ThrowEventsOff( );
 170   start = std::clock( );
 171   paths->Update( );
 172   end = std::clock( );
 173   seconds = double( end - start ) / double( CLOCKS_PER_SEC );
 174   std::cout << "Paths extraction time = " << seconds << std::endl;
 175
 176   std::cout << "Final seed: " << paths->GetSeed( 0 ) << std::endl;
 177
 178   // Create polydata
 179   vtkSmartPointer< vtkPoints > points =
 180     vtkSmartPointer< vtkPoints >::New( );
 181   vtkSmartPointer< vtkCellArray > cells =
 182     vtkSmartPointer< vtkCellArray >::New( );
 183   vtkSmartPointer< vtkFloatArray > scalars =
 184     vtkSmartPointer< vtkFloatArray >::New( );
 185
 186   const TDijkstra::TVertices& endpoints = paths->GetEndPoints( );
 187   const TDijkstra::TTree& tree = paths->GetFinalTree( );
 188   TDijkstra::TVertices::const_iterator epIt = endpoints.begin( );
 189   for( unsigned int epId = 0; epIt != endpoints.end( ); ++epIt, ++epId )
 190   {
 191     double pd = double( epId ) / double( endpoints.size( ) - 1 );
 192
 193     TDijkstra::TVertex idx = *epIt;
 194     do
 195     {
 196       TImage::PointType pnt;
 197       input_image->TransformIndexToPhysicalPoint( idx, pnt );
 198
 199       points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 200       scalars->InsertNextTuple1( pd );
 201       if( idx != *epIt )
 202       {
 203         cells->InsertNextCell( 2 );
 204         cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
 205         cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
 206
 207       } // fi
 208       idx = tree.find( idx )->second;
 209
 210     } while( idx != tree.find( idx )->second );
 211
 212   } // rof
 213
 214   vtkSmartPointer< vtkPolyData > vtk_tree =
 215     vtkSmartPointer< vtkPolyData >::New( );
 216   vtk_tree->SetPoints( points );
 217   vtk_tree->SetLines( cells );
 218   vtk_tree->GetPointData( )->SetScalars( scalars );
 219
 220   view.AddPolyData( vtk_tree );
 221   view.Render( );
 222   view.Start( );
 223
 224   /* TODO
 225      TDistanceMapWriter::Pointer distancemap_writer =
 226      TDistanceMapWriter::New( );
 227      distancemap_writer->SetInput( distanceMap->GetOutput( ) );
 228      distancemap_writer->SetFileName( output_distancemap_fn );
 229      distancemap_writer->Update( );
 230
 231      TImageWriter::Pointer segmentation_writer =
 232      TImageWriter::New( );
 233      segmentation_writer->SetInput( segmentor->GetOutput( ) );
 234      segmentation_writer->SetFileName( output_segmentation_fn );
 235      segmentation_writer->Update( );
 236   */
 237
 238   // Show result
 239   /*
 240     TVTKImage::Pointer output_image_vtk = TVTKImage::New( );
 241     output_image_vtk->SetInput( segmentor->GetOutput( ) );
 242     output_image_vtk->Update( );
 243
 244     vtkSmartPointer< vtkImageMarchingCubes > mc =
 245     vtkSmartPointer< vtkImageMarchingCubes >::New( );
 246     mc->SetInputData( output_image_vtk->GetOutput( ) );
 247     mc->SetValue(
 248     0,
 249     double( segmentor->GetInsideValue( ) ) * double( 0.95 )
 250     );
 251     mc->Update( );
 252
 253     // Let some interaction and close program
 254     view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 );
 255     view.Start( );
 256
 257     // Write resulting image
 258     TImageWriter::Pointer output_image_writer = TImageWriter::New( );
 259     output_image_writer->SetInput( segmentor->GetOutput( ) );
 260     output_image_writer->SetFileName( output_image_fn );
 261     try
 262     {
 263     output_image_writer->Update( );
 264     }
 265     catch( itk::ExceptionObject& err )
 266     {
 267     std::cerr << "Error caught: " << err << std::endl;
 268     return( 1 );
 269
 270     } // yrt
 271   */
 272   return( 0 );
 273 }
 274
 275 // eof - $RCSfile$