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