appli/examples/example_Image_Dijkstra_LabelSkeleton.cxx

   1 #include <cstdlib>
   2 #include <ctime>
   3 #include <iostream>
   4 #include <limits>
   5 #include <string>
   6
   7 #include <itkImage.h>
   8 #include <itkImageToVTKImageFilter.h>
   9
  10 #include <itkImageFileReader.h>
  11 #include <itkImageFileWriter.h>
  12
  13 #include <fpa/Image/Functors/ImageCostFunction.h>
  14 #include <fpa/Image/DijkstraWithEndPointDetection.h>
  15 #include <fpa/VTK/ImageMPR.h>
  16 #include <fpa/VTK/Image3DObserver.h>
  17
  18 #include <fpa/IO/MinimumSpanningTreeWriter.h>
  19 #include <fpa/IO/UniqueValuesContainerWriter.h>
  20 #include <fpa/IO/MatrixValuesContainerWriter.h>
  21
  22 #include <vtkCellArray.h>
  23 #include <vtkFloatArray.h>
  24 #include <vtkImageMarchingCubes.h>
  25 #include <vtkPoints.h>
  26 #include <vtkPolyData.h>
  27 #include <vtkSmartPointer.h>
  28
  29 // -------------------------------------------------------------------------
  30 const unsigned int Dim = 3;
  31 typedef unsigned char TPixel;
  32 typedef float         TScalar;
  33
  34 typedef itk::Image< TPixel, Dim >            TImage;
  35 typedef itk::Image< TScalar, Dim >           TScalarImage;
  36 typedef itk::ImageToVTKImageFilter< TImage > TVTKInputImage;
  37
  38 // -------------------------------------------------------------------------
  39 template< class I >
  40 void ReadImage( typename I::Pointer& image, const std::string& filename );
  41
  42 template< class I >
  43 void SaveImage( const I* image, const std::string& filename );
  44
  45 // -------------------------------------------------------------------------
  46 template< class F >
  47 class SkeletonCostFunction
  48   : public fpa::Image::Functors::ImageCostFunction< typename F::TInputImage, typename F::TResult >
  49 {
  50 public:
  51   typedef fpa::Image::Functors::ImageCostFunction< typename F::TInputImage, typename F::TResult > Superclass;
  52   typedef SkeletonCostFunction            Self;
  53   typedef itk::SmartPointer< Self >       Pointer;
  54   typedef itk::SmartPointer< const Self > ConstPointer;
  55
  56   typedef typename Superclass::TInputImage TInputImage;
  57   typedef typename Superclass::TResult     TResult;
  58   typedef typename Superclass::TIndex      TIndex;
  59
  60 public:
  61   itkNewMacro( Self );
  62   itkTypeMacro( SkeletonCostFunction, fpa_Image_Functors_ImageCostFunction );
  63
  64 public:
  65   virtual void SetInputImage( const TInputImage* img )
  66     {
  67       this->Superclass::SetInputImage( img );
  68
  69       typename TInputImage::SpacingType spac = img->GetSpacing( );
  70       typename TInputImage::SpacingType::ValueType min_spac = spac[ 0 ];
  71       for( unsigned int d = 1; d < TInputImage::ImageDimension; ++d )
  72         min_spac = ( spac[ d ] < min_spac )? spac[ d ]: min_spac;
  73       this->m_Distance = TResult( min_spac );
  74     }
  75
  76   virtual TResult Evaluate( const TIndex& v, const TIndex& p ) const
  77     {
  78       TResult res = this->Superclass::Evaluate( v, p );
  79       if( res > TResult( 0 ) )
  80         return( this->m_Distance );
  81       else
  82         return( TResult( -1 ) );
  83     }
  84
  85 protected:
  86   SkeletonCostFunction( )
  87     : Superclass( )
  88     { }
  89   virtual ~SkeletonCostFunction( )
  90     { }
  91
  92 private:
  93   // Purposely not implemented
  94   SkeletonCostFunction( const Self& );
  95   void operator=( const Self& );
  96
  97 protected:
  98   TResult m_Distance;
  99 };
 100
 101 // -------------------------------------------------------------------------
 102 int main( int argc, char* argv[] )
 103 {
 104   if( argc < 11 )
 105   {
 106     std::cerr
 107       << "Usage: " << argv[ 0 ] << std::endl
 108       << " input_image" << std::endl
 109       << " seed_x seed_y seed_z" << std::endl
 110       << " output_costmap" << std::endl
 111       << " output_labels" << std::endl
 112       << " output_minimum_spanning_tree" << std::endl
 113       << " output_endpoints" << std::endl
 114       << " output_bifurcations" << std::endl
 115       << " output_branches" << std::endl
 116       << std::endl;
 117     return( 1 );
 118
 119   } // fi
 120   std::string input_image_fn = argv[ 1 ];
 121   std::string output_costmap_fn = argv[ 5 ];
 122   std::string output_labels_fn = argv[ 6 ];
 123   std::string mst_output_fn = argv[ 7 ];
 124   std::string endpoints_output_fn = argv[ 8 ];
 125   std::string bifurcations_output_fn = argv[ 9 ];
 126   std::string branches_output_fn = argv[ 10 ];
 127
 128   // Get seed
 129   TImage::PointType pnt;
 130   pnt[ 0 ] = TImage::PointType::ValueType( std::atof( argv[ 2 ] ) );
 131   pnt[ 1 ] = TImage::PointType::ValueType( std::atof( argv[ 3 ] ) );
 132   pnt[ 2 ] = TImage::PointType::ValueType( std::atof( argv[ 4 ] ) );
 133
 134   // Read image
 135   TImage::Pointer input_image;
 136   try
 137   {
 138     ReadImage< TImage >( input_image, input_image_fn );
 139   }
 140   catch( itk::ExceptionObject& err )
 141   {
 142     std::cerr
 143       << "Error caught while reading \""
 144       << input_image_fn << "\": " << err
 145       << std::endl;
 146     return( 1 );
 147
 148   } // yrt
 149
 150   TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( );
 151   vtk_input_image->SetInput( input_image );
 152   vtk_input_image->Update( );
 153
 154   // Show input image and let some interaction
 155   fpa::VTK::ImageMPR view;
 156   view.SetBackground( 0.3, 0.2, 0.8 );
 157   view.SetSize( 500, 500 );
 158   view.SetImage( vtk_input_image->GetOutput( ) );
 159
 160   vtkSmartPointer< vtkImageMarchingCubes > mc =
 161     vtkSmartPointer< vtkImageMarchingCubes >::New( );
 162   mc->SetInputData( vtk_input_image->GetOutput( ) );
 163   mc->SetValue( 0, 1e-1 );
 164   mc->Update( );
 165   view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
 166
 167   // Allow some interaction and wait for, at least, one seed
 168   view.Render( );
 169   view.Start( );
 170
 171   // Transform seed
 172   TImage::IndexType seed;
 173   input_image->TransformPhysicalPointToIndex( pnt, seed );
 174
 175   // Prepare Dijkstra filter
 176   typedef fpa::Image::DijkstraWithEndPointDetection< TImage, TScalarImage > TFilter;
 177   typedef SkeletonCostFunction< TFilter > TFunction;
 178
 179   TFunction::Pointer function = TFunction::New( );
 180
 181   TFilter::Pointer filter = TFilter::New( );
 182   filter->SetInput( input_image );
 183   filter->SetCostFunction( function );
 184   filter->SetNeighborhoodOrder( 2 );
 185   filter->SetSafetyNeighborhoodSize( 0 );
 186   filter->StopAtOneFrontOff( );
 187   filter->CorrectSeedsOff( );
 188   filter->CorrectEndPointsOff( );
 189   filter->AddSeed( seed, TScalar( 0 ) );
 190
 191   // Prepare graphical debugger
 192   typedef fpa::VTK::Image3DObserver< TFilter, vtkRenderWindow > TDebugger;
 193   TDebugger::Pointer debugger = TDebugger::New( );
 194   debugger->SetRenderWindow( view.GetWindow( ) );
 195   debugger->SetRenderPercentage( 0.0000001 );
 196   filter->AddObserver( itk::AnyEvent( ), debugger );
 197   filter->ThrowEventsOn( );
 198
 199   // Go!
 200   std::time_t start, end;
 201   std::time( &start );
 202   filter->Update( );
 203   std::time( &end );
 204   std::cout
 205     << "Extraction time = "
 206     << std::difftime( end, start )
 207     << " s." << std::endl;
 208
 209   // Outputs
 210   const TFilter::TOutputImage* accumulated_costs = filter->GetOutput( );
 211   const TFilter::TLabelImage* labeled_image = filter->GetLabelImage( );
 212   const TFilter::TMinimumSpanningTree* mst = filter->GetMinimumSpanningTree( );
 213   const TFilter::TUniqueVertices* endpoints = filter->GetEndPoints( );
 214   const TFilter::TUniqueVertices* bifurcations = filter->GetBifurcations( );
 215   const TFilter::TBranches* branches = filter->GetBranches( );
 216   unsigned long nBranches = filter->GetNumberOfBranches( );
 217
 218   // Save outputs
 219   SaveImage( accumulated_costs, output_costmap_fn );
 220   SaveImage( labeled_image, output_labels_fn );
 221
 222   fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::Pointer
 223     mst_writer =
 224     fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::New( );
 225   mst_writer->SetInput( mst );
 226   mst_writer->SetFileName( mst_output_fn );
 227   mst_writer->Update( );
 228
 229   fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer
 230     endpoints_writer =
 231     fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( );
 232   endpoints_writer->SetInput( endpoints );
 233   endpoints_writer->SetFileName( endpoints_output_fn );
 234   endpoints_writer->Update( );
 235
 236   fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer
 237     bifurcations_writer =
 238     fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( );
 239   bifurcations_writer->SetInput( bifurcations );
 240   bifurcations_writer->SetFileName( bifurcations_output_fn );
 241   bifurcations_writer->Update( );
 242
 243   fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::Pointer
 244     branches_writer =
 245     fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::New( );
 246   branches_writer->SetInput( branches );
 247   branches_writer->SetFileName( branches_output_fn );
 248   branches_writer->Update( );
 249
 250   // Show endpoints
 251   vtkSmartPointer< vtkPoints > endpoints_points =
 252     vtkSmartPointer< vtkPoints >::New( );
 253   vtkSmartPointer< vtkCellArray > endpoints_cells =
 254     vtkSmartPointer< vtkCellArray >::New( );
 255   for(
 256     TFilter::TUniqueVertices::ConstIterator epIt = endpoints->Begin( );
 257     epIt != endpoints->End( );
 258     ++epIt
 259     )
 260   {
 261     TImage::PointType pnt;
 262     input_image->TransformIndexToPhysicalPoint( *epIt, pnt );
 263     endpoints_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 264     endpoints_cells->InsertNextCell( 1 );
 265     endpoints_cells->
 266       InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 1 );
 267
 268   } // rof
 269   vtkSmartPointer< vtkPolyData > endpoints_polydata =
 270     vtkSmartPointer< vtkPolyData >::New( );
 271   endpoints_polydata->SetPoints( endpoints_points );
 272   endpoints_polydata->SetVerts( endpoints_cells );
 273   view.AddPolyData( endpoints_polydata, 0, 0, 1, 1 );
 274
 275   // Show bifurcations
 276   vtkSmartPointer< vtkPoints > bifurcations_points =
 277     vtkSmartPointer< vtkPoints >::New( );
 278   vtkSmartPointer< vtkCellArray > bifurcations_cells =
 279     vtkSmartPointer< vtkCellArray >::New( );
 280   for(
 281     TFilter::TUniqueVertices::ConstIterator bfIt = bifurcations->Begin( );
 282     bfIt != bifurcations->End( );
 283     ++bfIt
 284     )
 285   {
 286     TImage::PointType pnt;
 287     input_image->TransformIndexToPhysicalPoint( *bfIt, pnt );
 288     bifurcations_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 289     bifurcations_cells->InsertNextCell( 1 );
 290     bifurcations_cells->
 291       InsertCellPoint( bifurcations_points->GetNumberOfPoints( ) - 1 );
 292
 293   } // rof
 294   vtkSmartPointer< vtkPolyData > bifurcations_polydata =
 295     vtkSmartPointer< vtkPolyData >::New( );
 296   bifurcations_polydata->SetPoints( bifurcations_points );
 297   bifurcations_polydata->SetVerts( bifurcations_cells );
 298   view.AddPolyData( bifurcations_polydata, 0, 1, 0, 1 );
 299
 300   // Show branches (simple and detailed)
 301   vtkSmartPointer< vtkPoints > simple_branches_points =
 302     vtkSmartPointer< vtkPoints >::New( );
 303   vtkSmartPointer< vtkCellArray > simple_branches_cells =
 304     vtkSmartPointer< vtkCellArray >::New( );
 305
 306   vtkSmartPointer< vtkPoints > detailed_branches_points =
 307     vtkSmartPointer< vtkPoints >::New( );
 308   vtkSmartPointer< vtkCellArray > detailed_branches_cells =
 309     vtkSmartPointer< vtkCellArray >::New( );
 310   vtkSmartPointer< vtkFloatArray > detailed_branches_scalars =
 311     vtkSmartPointer< vtkFloatArray >::New( );
 312
 313   TFilter::TBranches::ConstIterator brIt = branches->Begin( );
 314   for( ; brIt != branches->End( ); ++brIt )
 315   {
 316     // Branch's first point
 317     TImage::PointType first_point;
 318     input_image->TransformIndexToPhysicalPoint( brIt->first, first_point );
 319     unsigned long first_id = simple_branches_points->GetNumberOfPoints( );
 320     simple_branches_points->InsertNextPoint(
 321       first_point[ 0 ], first_point[ 1 ], first_point[ 2 ]
 322       );
 323
 324     TFilter::TBranches::ConstRowIterator brRowIt = branches->Begin( brIt );
 325     for( ; brRowIt != branches->End( brIt ); ++brRowIt )
 326     {
 327       // Branch's second point
 328       TImage::PointType second_point;
 329       input_image->
 330         TransformIndexToPhysicalPoint( brRowIt->first, second_point );
 331       unsigned long second_id = simple_branches_points->GetNumberOfPoints( );
 332       simple_branches_points->InsertNextPoint(
 333         second_point[ 0 ], second_point[ 1 ], second_point[ 2 ]
 334         );
 335       simple_branches_cells->InsertNextCell( 2 );
 336       simple_branches_cells->InsertCellPoint( first_id );
 337       simple_branches_cells->InsertCellPoint( second_id );
 338
 339       // Detailed path
 340       double pathId = double( brRowIt->second - 1 ) / double( nBranches - 1 );
 341       TFilter::TVertices path;
 342       mst->GetPath( path, brIt->first, brRowIt->first );
 343       TFilter::TVertices::const_iterator pIt = path.begin( );
 344       for( ; pIt != path.end( ); ++pIt )
 345       {
 346         TImage::PointType path_point;
 347         input_image->TransformIndexToPhysicalPoint( *pIt, path_point );
 348         detailed_branches_points->InsertNextPoint(
 349           path_point[ 0 ], path_point[ 1 ], path_point[ 2 ]
 350           );
 351         detailed_branches_scalars->InsertNextTuple1( pathId );
 352         if( pIt != path.begin( ) )
 353         {
 354           unsigned long nPoints =
 355             detailed_branches_points->GetNumberOfPoints( );
 356           detailed_branches_cells->InsertNextCell( 2 );
 357           detailed_branches_cells->InsertCellPoint( nPoints - 2 );
 358           detailed_branches_cells->InsertCellPoint( nPoints - 1 );
 359
 360         } // fi
 361
 362       } // rof
 363
 364     } // rof
 365
 366   } // rof
 367   vtkSmartPointer< vtkPolyData > simple_branches_polydata =
 368     vtkSmartPointer< vtkPolyData >::New( );
 369   simple_branches_polydata->SetPoints( simple_branches_points );
 370   simple_branches_polydata->SetLines( simple_branches_cells );
 371   view.AddPolyData( simple_branches_polydata, 1, 0, 1, 1 );
 372
 373   vtkSmartPointer< vtkPolyData > detailed_branches_polydata =
 374     vtkSmartPointer< vtkPolyData >::New( );
 375   detailed_branches_polydata->SetPoints( detailed_branches_points );
 376   detailed_branches_polydata->SetLines( detailed_branches_cells );
 377   detailed_branches_polydata->
 378     GetPointData( )->SetScalars( detailed_branches_scalars );
 379   view.AddPolyData( detailed_branches_polydata, 1 );
 380
 381   // Let some more interaction
 382   view.Render( );
 383   view.Start( );
 384
 385   return( 0 );
 386 }
 387
 388 // -------------------------------------------------------------------------
 389 template< class I >
 390 void ReadImage( typename I::Pointer& image, const std::string& filename )
 391 {
 392   typename itk::ImageFileReader< I >::Pointer reader =
 393     itk::ImageFileReader< I >::New( );
 394   reader->SetFileName( filename );
 395   reader->Update( );
 396   image = reader->GetOutput( );
 397   image->DisconnectPipeline( );
 398 }
 399
 400 // -------------------------------------------------------------------------
 401 template< class I >
 402 void SaveImage( const I* image, const std::string& filename )
 403 {
 404   typename itk::ImageFileWriter< I >::Pointer writer =
 405     itk::ImageFileWriter< I >::New( );
 406   writer->SetInput( image );
 407   writer->SetFileName( filename );
 408   try
 409   {
 410     writer->Update( );
 411   }
 412   catch( itk::ExceptionObject& err )
 413   {
 414     std::cerr
 415       << "Error saving \"" << filename << "\": " << err
 416       << std::endl;
 417
 418   } // yrt
 419 }
 420
 421 // eof - $RCSfile$