appli/examples/example_ImageAlgorithmDijkstra_03.cxx

   1 #include <cmath>
   2 #include <ctime>
   3 #include <deque>
   4 #include <iostream>
   5 #include <limits>
   6 #include <map>
   7 #include <string>
   8
   9 #include <itkConstNeighborhoodIterator.h>
  10 #include <itkNeighborhoodIterator.h>
  11 #include <itkDanielssonDistanceMapImageFilter.h>
  12 #include <itkImage.h>
  13 #include <itkImageFileReader.h>
  14 #include <itkImageFileWriter.h>
  15 #include <itkImageToVTKImageFilter.h>
  16
  17 #include <vtkPoints.h>
  18 #include <vtkCellArray.h>
  19 #include <vtkFloatArray.h>
  20 #include <vtkPolyData.h>
  21 #include <vtkSmartPointer.h>
  22
  23 #include <fpa/Image/DijkstraWithSphereBacktracking.h>
  24 #include <fpa/VTK/ImageMPR.h>
  25 #include <fpa/VTK/Image3DObserver.h>
  26
  27 // -------------------------------------------------------------------------
  28 const unsigned int Dim = 3;
  29 typedef double                               TPixel;
  30 typedef double                               TScalar;
  31 typedef itk::Image< TPixel, Dim >            TImage;
  32 typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
  33
  34 typedef itk::ImageFileReader< TImage >          TImageReader;
  35 typedef itk::ImageFileWriter< TImage >          TImageWriter;
  36 typedef fpa::Image::DijkstraWithSphereBacktracking< TImage, TScalar > TDijkstra;
  37
  38 typedef fpa::VTK::ImageMPR TMPR;
  39 typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs;
  40
  41 // -------------------------------------------------------------------------
  42 int main( int argc, char* argv[] )
  43 {
  44   if( argc < 4 )
  45   {
  46     std::cerr
  47       << "Usage: " << argv[ 0 ]
  48       << " input_image output_image"
  49       << " neighborhood_order"
  50       << " visual_debug"
  51       << std::endl;
  52     return( 1 );
  53
  54   } // fi
  55   std::string input_image_fn = argv[ 1 ];
  56   std::string output_image_fn = argv[ 2 ];
  57   unsigned int neighborhood_order = std::atoi( argv[ 3 ] );
  58   bool visual_debug = false;
  59   if( argc > 4 )
  60     visual_debug = ( std::atoi( argv[ 4 ] ) == 1 );
  61
  62   // Read image
  63   TImageReader::Pointer input_image_reader = TImageReader::New( );
  64   input_image_reader->SetFileName( input_image_fn );
  65   try
  66   {
  67     input_image_reader->Update( );
  68   }
  69   catch( itk::ExceptionObject& err )
  70   {
  71     std::cerr << "Error caught: " << err << std::endl;
  72     return( 1 );
  73
  74   } // yrt
  75   TImage::ConstPointer input_image = input_image_reader->GetOutput( );
  76
  77   // Show image
  78   TVTKImage::Pointer vtk_image = TVTKImage::New( );
  79   vtk_image->SetInput( input_image );
  80   vtk_image->Update( );
  81
  82   TMPR view;
  83   view.SetBackground( 0.3, 0.2, 0.8 );
  84   view.SetSize( 800, 800 );
  85   view.SetImage( vtk_image->GetOutput( ) );
  86
  87   // Wait for a seed to be given
  88   while( view.GetNumberOfSeeds( ) == 0 )
  89     view.Start( );
  90
  91   // Get seed from user
  92   double seed[ 3 ];
  93   view.GetSeed( 0, seed );
  94   TImage::PointType seed_pnt;
  95   seed_pnt[ 0 ] = seed[ 0 ];
  96   seed_pnt[ 1 ] = seed[ 1 ];
  97   seed_pnt[ 2 ] = seed[ 2 ];
  98   TImage::IndexType seed_idx;
  99   input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx );
 100
 101   // Extract paths
 102   TDijkstra::Pointer paths = TDijkstra::New( );
 103   paths->AddSeed( seed_idx, TScalar( 0 ) );
 104   paths->SetInput( input_image );
 105   paths->SetNeighborhoodOrder( neighborhood_order );
 106
 107   if( visual_debug )
 108   {
 109     // Configure observer
 110     TDijkstraObs::Pointer obs = TDijkstraObs::New( );
 111     obs->SetRenderWindow( view.GetWindow( ) );
 112     paths->AddObserver( itk::AnyEvent( ), obs );
 113     paths->ThrowEventsOn( );
 114   }
 115   else
 116     paths->ThrowEventsOff( );
 117   std::clock_t start = std::clock( );
 118   paths->Update( );
 119   std::clock_t end = std::clock( );
 120   double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
 121   std::cout << "Paths extraction time = " << seconds << std::endl;
 122
 123   // Create polydata
 124   vtkSmartPointer< vtkPoints > points =
 125     vtkSmartPointer< vtkPoints >::New( );
 126   vtkSmartPointer< vtkCellArray > cells =
 127     vtkSmartPointer< vtkCellArray >::New( );
 128   vtkSmartPointer< vtkFloatArray > scalars =
 129     vtkSmartPointer< vtkFloatArray >::New( );
 130
 131   const TDijkstra::TVertices& endpoints = paths->GetEndPoints( );
 132   const TDijkstra::TTree& tree = paths->GetFinalTree( );
 133   TDijkstra::TVertices::const_iterator epIt = endpoints.begin( );
 134   for( unsigned int epId = 0; epIt != endpoints.end( ); ++epIt, ++epId )
 135   {
 136     double pd = double( epId ) / double( endpoints.size( ) - 1 );
 137
 138     TDijkstra::TVertex idx = *epIt;
 139     do
 140     {
 141       TImage::PointType pnt;
 142       input_image->TransformIndexToPhysicalPoint( idx, pnt );
 143
 144       points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 145       scalars->InsertNextTuple1( pd );
 146       if( idx != *epIt )
 147       {
 148         cells->InsertNextCell( 2 );
 149         cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
 150         cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
 151
 152       } // fi
 153       idx = tree.find( idx )->second;
 154
 155     } while( idx != tree.find( idx )->second );
 156
 157   } // rof
 158
 159   vtkSmartPointer< vtkPolyData > vtk_tree =
 160     vtkSmartPointer< vtkPolyData >::New( );
 161   vtk_tree->SetPoints( points );
 162   vtk_tree->SetLines( cells );
 163   vtk_tree->GetPointData( )->SetScalars( scalars );
 164
 165   view.AddPolyData( vtk_tree );
 166   view.Render( );
 167   view.Start( );
 168
 169   /* TODO
 170      TImageWriter::Pointer dijkstra_writer =
 171      TImageWriter::New( );
 172      dijkstra_writer->SetInput( paths->GetOutput( ) );
 173      dijkstra_writer->SetFileName( "dijkstra.mhd" );
 174      dijkstra_writer->Update( );
 175   */
 176
 177   return( 0 );
 178 }
 179
 180 // eof - $RCSfile$