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 <itkMinimumMaximumImageCalculator.h>
  18 #include <itkInvertIntensityImageFilter.h>
  19
  20 #include <vtkPoints.h>
  21 #include <vtkCellArray.h>
  22 #include <vtkFloatArray.h>
  23 #include <vtkPolyData.h>
  24 #include <vtkPolyDataWriter.h>
  25 #include <vtkSmartPointer.h>
  26 #include <vtkImageMarchingCubes.h>
  27 #include <vtkLookupTable.h>
  28
  29 #include <fpa/Image/DijkstraWithSphereBacktracking.h>
  30 #include <fpa/VTK/ImageMPR.h>
  31 #include <fpa/VTK/Image3DObserver.h>
  32
  33 // -------------------------------------------------------------------------
  34 const unsigned int Dim = 3;
  35 typedef float                                TPixel;
  36 typedef float                                TScalar;
  37 typedef itk::Image< TPixel, Dim >            TImage;
  38 typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
  39
  40 typedef itk::ImageFileReader< TImage >          TImageReader;
  41 typedef itk::ImageFileWriter< TImage >    TImageWriter;
  42 typedef fpa::Image::DijkstraWithSphereBacktracking< TImage, TScalar > TDijkstra;
  43
  44 typedef fpa::VTK::ImageMPR TMPR;
  45 typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs;
  46
  47 // -------------------------------------------------------------------------
  48 int main( int argc, char* argv[] )
  49 {
  50   if( argc < 6 )
  51   {
  52     std::cerr
  53       << "Usage: " << argv[ 0 ]
  54       << " input_image x y z output_image"
  55       << " visual_debug"
  56       << std::endl;
  57     return( 1 );
  58
  59   } // fi
  60   std::string input_image_fn = argv[ 1 ];
  61   TImage::PointType seed_pnt;
  62   seed_pnt[ 0 ] = std::atof( argv[ 2 ] );
  63   seed_pnt[ 1 ] = std::atof( argv[ 3 ] );
  64   seed_pnt[ 2 ] = std::atof( argv[ 4 ] );
  65   std::string output_image_fn = argv[ 5 ];
  66   bool visual_debug = false;
  67   if( argc > 6 )
  68     visual_debug = ( std::atoi( argv[ 6 ] ) == 1 );
  69
  70   // Read image
  71   TImageReader::Pointer input_image_reader = TImageReader::New( );
  72   input_image_reader->SetFileName( input_image_fn );
  73   try
  74   {
  75     input_image_reader->Update( );
  76   }
  77   catch( itk::ExceptionObject& err )
  78   {
  79     std::cerr << "Error caught: " << err << std::endl;
  80     return( 1 );
  81
  82   } // yrt
  83   TImage::ConstPointer input_image = input_image_reader->GetOutput( );
  84
  85   // Show image
  86   TVTKImage::Pointer vtk_image = TVTKImage::New( );
  87   vtk_image->SetInput( input_image );
  88   vtk_image->Update( );
  89
  90   TMPR view;
  91   view.SetBackground( 0.3, 0.2, 0.8 );
  92   view.SetSize( 800, 800 );
  93   view.SetImage( vtk_image->GetOutput( ) );
  94
  95   vtkSmartPointer< vtkImageMarchingCubes > mc =
  96     vtkSmartPointer< vtkImageMarchingCubes >::New( );
  97   mc->SetInputData( vtk_image->GetOutput( ) );
  98   mc->SetValue( 0, 1e-2 );
  99   mc->Update( );
 100   view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
 101
 102   // Allow some interaction
 103   view.Render( );
 104   view.Start( );
 105
 106   TImage::IndexType seed_idx;
 107   input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx );
 108   std::cout << seed_idx << " " << seed_pnt << std::endl;
 109   seed_idx[ 0 ] = 256;
 110   seed_idx[ 1 ] = 313;
 111   seed_idx[ 2 ] = 381;
 112
 113   // Extract paths
 114   TDijkstra::Pointer paths = TDijkstra::New( );
 115   paths->AddSeed( seed_idx, TScalar( 0 ) );
 116   paths->SetInput( input_image );
 117   paths->SetNeighborhoodOrder( 2 );
 118
 119   if( visual_debug )
 120   {
 121     // Configure observer
 122     TDijkstraObs::Pointer obs = TDijkstraObs::New( );
 123     obs->SetRenderWindow( view.GetWindow( ) );
 124     paths->AddObserver( itk::AnyEvent( ), obs );
 125     paths->ThrowEventsOn( );
 126   }
 127   else
 128     paths->ThrowEventsOff( );
 129   std::clock_t start = std::clock( );
 130   paths->Update( );
 131   std::clock_t end = std::clock( );
 132   double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
 133   std::cout << "Paths extraction time = " << seconds << std::endl;
 134
 135   // Create polydata
 136   vtkSmartPointer< vtkPoints > points =
 137     vtkSmartPointer< vtkPoints >::New( );
 138   vtkSmartPointer< vtkCellArray > cells =
 139     vtkSmartPointer< vtkCellArray >::New( );
 140   vtkSmartPointer< vtkFloatArray > scalars =
 141     vtkSmartPointer< vtkFloatArray >::New( );
 142
 143   const TDijkstra::TMarkImage* marks = paths->GetOutputMarkImage( );
 144   TDijkstra::TMark max_mark = paths->GetNumberOfBranches( );
 145   const TDijkstra::TVertices& endpoints = paths->GetEndPoints( );
 146   const TDijkstra::TTree& tree = paths->GetFinalTree( );
 147   TDijkstra::TVertices::const_iterator epIt = endpoints.begin( );
 148   for( unsigned int epId = 0; epIt != endpoints.end( ); ++epIt, ++epId )
 149   {
 150     TDijkstra::TVertex idx = *epIt;
 151     do
 152     {
 153       TImage::PointType pnt;
 154       input_image->TransformIndexToPhysicalPoint( idx, pnt );
 155
 156       TDijkstra::TMark mark = marks->GetPixel( idx );
 157       double pd = double( mark );
 158
 159       points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 160       scalars->InsertNextTuple1( pd );
 161       if( idx != *epIt )
 162       {
 163         cells->InsertNextCell( 2 );
 164         cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
 165         cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
 166
 167       } // fi
 168       idx = tree.find( idx )->second;
 169
 170     } while( idx != tree.find( idx )->second );
 171
 172   } // rof
 173
 174   vtkSmartPointer< vtkPolyData > vtk_tree =
 175     vtkSmartPointer< vtkPolyData >::New( );
 176   vtk_tree->SetPoints( points );
 177   vtk_tree->SetLines( cells );
 178   vtk_tree->GetPointData( )->SetScalars( scalars );
 179
 180   vtkSmartPointer<vtkLookupTable> lut =
 181     vtkSmartPointer<vtkLookupTable>::New( );
 182   lut->SetNumberOfTableValues( max_mark + 1 );
 183   lut->SetTableRange( 0, max_mark );
 184   lut->Build( );
 185
 186   view.AddPolyData( vtk_tree, lut );
 187   view.Render( );
 188   view.Start( );
 189
 190   vtkSmartPointer< vtkPolyDataWriter > writer =
 191     vtkSmartPointer< vtkPolyDataWriter >::New( );
 192   writer->SetInputData( vtk_tree );
 193   writer->SetFileName( output_image_fn.c_str( ) );
 194   writer->Update( );
 195
 196   itk::ImageFileWriter< TDijkstra::TMarkImage >::Pointer marks_writer =
 197     itk::ImageFileWriter< TDijkstra::TMarkImage >::New( );
 198   marks_writer->SetInput( marks );
 199   marks_writer->SetFileName( "marks.mhd" );
 200   marks_writer->Update( );
 201
 202   /* TODO
 203      TImageWriter::Pointer dijkstra_writer =
 204      TImageWriter::New( );
 205      dijkstra_writer->SetInput( paths->GetOutput( ) );
 206      dijkstra_writer->SetFileName( "dijkstra.mhd" );
 207      dijkstra_writer->Update( );
 208   */
 209
 210   return( 0 );
 211 }
 212
 213 // eof - $RCSfile$