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 <vtkCellArray.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 template< class I, class O >
  46 void DistanceMap(
  47   const typename I::Pointer& input, typename O::Pointer& output
  48   );
  49
  50 // -------------------------------------------------------------------------
  51 int main( int argc, char* argv[] )
  52 {
  53   if( argc < 5 )
  54   {
  55     std::cerr
  56       << "Usage: " << argv[ 0 ]
  57       << " input_image distancemap output_costmap output_labels"
  58       << std::endl;
  59     return( 1 );
  60
  61   } // fi
  62   std::string input_image_fn = argv[ 1 ];
  63   std::string distancemap_fn = argv[ 2 ];
  64   std::string output_costmap_fn = argv[ 3 ];
  65   std::string output_labels_fn = argv[ 4 ];
  66
  67   // Read image
  68   TImage::Pointer input_image;
  69   try
  70   {
  71     ReadImage< TImage >( input_image, input_image_fn );
  72   }
  73   catch( itk::ExceptionObject& err )
  74   {
  75     std::cerr
  76       << "Error caught while reading \""
  77       << input_image_fn << "\": " << err
  78       << std::endl;
  79     return( 1 );
  80
  81   } // yrt
  82
  83   TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( );
  84   vtk_input_image->SetInput( input_image );
  85   vtk_input_image->Update( );
  86
  87   // Show input image and let some interaction
  88   fpa::VTK::ImageMPR view;
  89   view.SetBackground( 0.3, 0.2, 0.8 );
  90   view.SetSize( 800, 800 );
  91   view.SetImage( vtk_input_image->GetOutput( ) );
  92
  93   vtkSmartPointer< vtkImageMarchingCubes > mc =
  94     vtkSmartPointer< vtkImageMarchingCubes >::New( );
  95   mc->SetInputData( vtk_input_image->GetOutput( ) );
  96   mc->SetValue( 0, 1e-1 );
  97   mc->Update( );
  98   view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
  99
 100   // Allow some interaction and wait for, at least, one seed
 101   view.Render( );
 102   while( view.GetNumberOfSeeds( ) == 0 )
 103     view.Start( );
 104
 105   // Get seed
 106   double p[ 3 ];
 107   view.GetSeed( 0, p );
 108   TImage::PointType pnt;
 109   TImage::IndexType seed;
 110   pnt[ 0 ] = TImage::PointType::ValueType( p[ 0 ] );
 111   pnt[ 1 ] = TImage::PointType::ValueType( p[ 1 ] );
 112   pnt[ 2 ] = TImage::PointType::ValueType( p[ 2 ] );
 113   input_image->TransformPhysicalPointToIndex( pnt, seed );
 114
 115   // Compute squared distance map
 116   TScalarImage::Pointer dmap;
 117   DistanceMap< TImage, TScalarImage >( input_image, dmap );
 118
 119   // Prepare cost conversion function
 120   typedef fpa::Base::Functors::InvertCostFunction< TScalar > TFunction;
 121   TFunction::Pointer function = TFunction::New( );
 122
 123   // Prepare Dijkstra filter
 124   typedef fpa::Image::DijkstraWithEndPointDetection< TScalarImage, TScalarImage > TFilter;
 125   TFilter::Pointer filter = TFilter::New( );
 126   filter->SetInput( dmap );
 127   filter->SetConversionFunction( function );
 128   filter->SetNeighborhoodOrder( 2 );
 129   filter->StopAtOneFrontOff( );
 130   filter->AddSeed( seed, TScalar( 0 ) );
 131
 132   // Prepare graphical debugger
 133   typedef fpa::VTK::Image3DObserver< TFilter, vtkRenderWindow > TDebugger;
 134   TDebugger::Pointer debugger = TDebugger::New( );
 135   debugger->SetRenderWindow( view.GetWindow( ) );
 136   debugger->SetRenderPercentage( 0.0001 );
 137   filter->AddObserver( itk::AnyEvent( ), debugger );
 138   filter->ThrowEventsOn( );
 139
 140   // Go!
 141   std::time_t start, end;
 142   std::time( &start );
 143   filter->Update( );
 144   std::time( &end );
 145   std::cout
 146     << "Extraction time = "
 147     << std::difftime( end, start )
 148     << " s." << std::endl;
 149
 150   // Minimum spanning tree
 151   const TFilter::TMinimumSpanningTree* mst =
 152     filter->GetMinimumSpanningTree( );
 153
 154   // Build branches from endpoints to seed
 155   vtkSmartPointer< vtkPoints > endpoints_points =
 156     vtkSmartPointer< vtkPoints >::New( );
 157   vtkSmartPointer< vtkCellArray > endpoints_cells =
 158     vtkSmartPointer< vtkCellArray >::New( );
 159   vtkSmartPointer< vtkCellArray > endpoints_vertices =
 160     vtkSmartPointer< vtkCellArray >::New( );
 161
 162   const TFilter::TUniqueVertices& endpoints = filter->GetEndPoints( );
 163   TFilter::TUniqueVertices::const_iterator eIt = endpoints.begin( );
 164   for( ; eIt != endpoints.end( ); ++eIt )
 165   {
 166     TFilter::TVertices path;
 167     mst->GetPath( path, *eIt, filter->GetSeed( 0 ) );
 168
 169     std::cout << *eIt << std::endl;
 170     TFilter::TVertices::const_iterator pIt = path.begin( );
 171     for( ; pIt != path.end( ); ++pIt )
 172     {
 173       TImage::PointType pnt;
 174       input_image->TransformIndexToPhysicalPoint( *pIt, pnt );
 175       endpoints_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 176
 177       if( pIt != path.begin( ) )
 178       {
 179         endpoints_cells->InsertNextCell( 2 );
 180         endpoints_cells->InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 2 );
 181         endpoints_cells->InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 1 );
 182       }
 183       else
 184       {
 185         std::cout << "ok" << std::endl;
 186         endpoints_vertices->InsertNextCell( 1 );
 187         endpoints_vertices->InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 1 );
 188
 189       } // fi
 190
 191
 192
 193
 194     } // rof
 195
 196   } // rof
 197   std::cout << endpoints.size( ) << std::endl;
 198
 199   vtkSmartPointer< vtkPolyData > endpoints_polydata =
 200     vtkSmartPointer< vtkPolyData >::New( );
 201   endpoints_polydata->SetPoints( endpoints_points );
 202   endpoints_polydata->SetLines( endpoints_cells );
 203   endpoints_polydata->SetVerts( endpoints_vertices );
 204
 205
 206
 207
 208
 209
 210
 211
 212
 213
 214
 215   // Bifurcations
 216   vtkSmartPointer< vtkPoints > bifurcations_points =
 217     vtkSmartPointer< vtkPoints >::New( );
 218   vtkSmartPointer< vtkCellArray > bifurcations_vertices =
 219     vtkSmartPointer< vtkCellArray >::New( );
 220
 221   const TFilter::TUniqueVertices& bifurcations = filter->GetBifurcationPoints( );
 222   TFilter::TUniqueVertices::const_iterator bIt = bifurcations.begin( );
 223   for( ; bIt != bifurcations.end( ); ++bIt )
 224   {
 225     std::cout << *bIt << std::endl;
 226
 227     TImage::PointType pnt;
 228     input_image->TransformIndexToPhysicalPoint( *bIt, pnt );
 229     bifurcations_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 230     bifurcations_vertices->InsertNextCell( 1 );
 231     bifurcations_vertices->InsertCellPoint( bifurcations_points->GetNumberOfPoints( ) - 1 );
 232
 233   } // rof
 234
 235   vtkSmartPointer< vtkPolyData > bifurcations_polydata =
 236     vtkSmartPointer< vtkPolyData >::New( );
 237   bifurcations_polydata->SetPoints( bifurcations_points );
 238   bifurcations_polydata->SetVerts( bifurcations_vertices );
 239
 240
 241
 242
 243
 244   // Build branches from endpoints to seed
 245   vtkSmartPointer< vtkPoints > branches_points =
 246     vtkSmartPointer< vtkPoints >::New( );
 247   vtkSmartPointer< vtkCellArray > branches_cells =
 248     vtkSmartPointer< vtkCellArray >::New( );
 249
 250   const TFilter::TBranches& branches = filter->GetBranches( );
 251   TFilter::TBranches::const_iterator brIt = branches.begin( );
 252   for( ; brIt != branches.end( ); ++brIt )
 253   {
 254     TFilter::TBranch::const_iterator brsIt = brIt->second.begin( );
 255     for( ; brsIt != brIt->second.end( ); ++brsIt )
 256     {
 257       TImage::PointType pnt0, pnt1;
 258       input_image->TransformIndexToPhysicalPoint( brIt->first, pnt0 );
 259       input_image->TransformIndexToPhysicalPoint( brsIt->first, pnt1 );
 260       branches_points->InsertNextPoint( pnt0[ 0 ], pnt0[ 1 ], pnt0[ 2 ] );
 261       branches_points->InsertNextPoint( pnt1[ 0 ], pnt1[ 1 ], pnt1[ 2 ] );
 262
 263       branches_cells->InsertNextCell( 2 );
 264       branches_cells->InsertCellPoint( branches_points->GetNumberOfPoints( ) - 2 );
 265       branches_cells->InsertCellPoint( branches_points->GetNumberOfPoints( ) - 1 );
 266
 267     } // rof
 268
 269   } // rof
 270
 271   vtkSmartPointer< vtkPolyData > branches_polydata =
 272     vtkSmartPointer< vtkPolyData >::New( );
 273   branches_polydata->SetPoints( branches_points );
 274   branches_polydata->SetLines( branches_cells );
 275
 276   view.AddPolyData( endpoints_polydata, 1, 0, 0, 1 );
 277   view.AddPolyData( bifurcations_polydata, 0, 0, 1, 1 );
 278   view.AddPolyData( branches_polydata, 0, 1, 0, 1 );
 279   view.Render( );
 280   view.Start( );
 281
 282   // Save final total cost map
 283   SaveImage( filter->GetOutput( ), output_costmap_fn );
 284   SaveImage( dmap.GetPointer( ), distancemap_fn );
 285   SaveImage( filter->GetLabelImage( ), output_labels_fn );
 286   return( 0 );
 287 }
 288
 289 // -------------------------------------------------------------------------
 290 template< class I >
 291 void ReadImage( typename I::Pointer& image, const std::string& filename )
 292 {
 293   typename itk::ImageFileReader< I >::Pointer reader =
 294     itk::ImageFileReader< I >::New( );
 295   reader->SetFileName( filename );
 296   reader->Update( );
 297
 298   typename itk::MinimumMaximumImageCalculator< I >::Pointer minmax =
 299     itk::MinimumMaximumImageCalculator< I >::New( );
 300   minmax->SetImage( reader->GetOutput( ) );
 301   minmax->Compute( );
 302   double vmin = double( minmax->GetMinimum( ) );
 303   double vmax = double( minmax->GetMaximum( ) );
 304
 305   typename itk::ShiftScaleImageFilter< I, I >::Pointer shift =
 306     itk::ShiftScaleImageFilter< I, I >::New( );
 307   shift->SetInput( reader->GetOutput( ) );
 308   shift->SetScale( vmax - vmin );
 309   shift->SetShift( vmin / ( vmax - vmin ) );
 310   shift->Update( );
 311
 312   image = shift->GetOutput( );
 313   image->DisconnectPipeline( );
 314 }
 315
 316 // -------------------------------------------------------------------------
 317 template< class I >
 318 void SaveImage( const I* image, const std::string& filename )
 319 {
 320   typename itk::ImageFileWriter< I >::Pointer writer =
 321     itk::ImageFileWriter< I >::New( );
 322   writer->SetInput( image );
 323   writer->SetFileName( filename );
 324   try
 325   {
 326     writer->Update( );
 327   }
 328   catch( itk::ExceptionObject& err )
 329   {
 330     std::cerr
 331       << "Error saving \"" << filename << "\": " << err
 332       << std::endl;
 333
 334   } // yrt
 335 }
 336
 337 // -------------------------------------------------------------------------
 338 template< class I, class O >
 339 void DistanceMap(
 340   const typename I::Pointer& input, typename O::Pointer& output
 341   )
 342 {
 343   typename itk::SignedMaurerDistanceMapImageFilter< I, O >::Pointer dmap =
 344     itk::SignedMaurerDistanceMapImageFilter< I, O >::New( );
 345   dmap->SetInput( input );
 346   dmap->SetBackgroundValue( ( typename I::PixelType )( 0 ) );
 347   dmap->InsideIsPositiveOn( );
 348   dmap->SquaredDistanceOn( );
 349   dmap->UseImageSpacingOn( );
 350
 351   std::time_t start, end;
 352   std::time( &start );
 353   dmap->Update( );
 354   std::time( &end );
 355   std::cout
 356     << "Distance map time = "
 357     << std::difftime( end, start )
 358     << " s." << std::endl;
 359
 360   output = dmap->GetOutput( );
 361   output->DisconnectPipeline( );
 362 }
 363
 364 // eof - $RCSfile$
 365
 366
 367
 368
 369
 370
 371
 372
 373
 374