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 /*
  24   #include <itkMinimumMaximumImageCalculator.h>
  25   #include <itkInvertIntensityImageFilter.h>
  26   #include <itkDanielssonDistanceMapImageFilter.h>
  27   #include <itkImageFileWriter.h>
  28
  29   #include <vtkCamera.h>
  30   #include <vtkImageActor.h>
  31   #include <vtkInteractorStyleImage.h>
  32   #include <vtkPointHandleRepresentation3D.h>
  33   #include <vtkProperty.h>
  34   #include <vtkRenderer.h>
  35   #include <vtkRenderWindow.h>
  36   #include <vtkRenderWindowInteractor.h>
  37   #include <vtkSeedRepresentation.h>
  38   #include <vtkSeedWidget.h>
  39   #include <vtkSmartPointer.h>
  40
  41   #include <fpa/Image/Dijkstra.h>
  42 */
  43
  44 // -------------------------------------------------------------------------
  45 const unsigned int Dim = 3;
  46 typedef unsigned char TPixel;
  47 typedef float         TScalar;
  48
  49 typedef itk::Image< TPixel, Dim >            TImage;
  50 typedef itk::Image< TScalar, Dim >           TScalarImage;
  51 typedef itk::ImageToVTKImageFilter< TImage > TVTKInputImage;
  52
  53 // -------------------------------------------------------------------------
  54 template< class I >
  55 void ReadImage( typename I::Pointer& image, const std::string& filename );
  56
  57 template< class I, class O >
  58 void DistanceMap(
  59   const typename I::Pointer& input, typename O::Pointer& output
  60   );
  61
  62 // -------------------------------------------------------------------------
  63 int main( int argc, char* argv[] )
  64 {
  65   if( argc < 2 )
  66   {
  67     std::cerr
  68       << "Usage: " << argv[ 0 ]
  69       << " input_image"
  70       << std::endl;
  71     return( 1 );
  72
  73   } // fi
  74   std::string input_image_fn = argv[ 1 ];
  75
  76   // Read image
  77   TImage::Pointer input_image;
  78   try
  79   {
  80     ReadImage< TImage >( input_image, input_image_fn );
  81   }
  82   catch( itk::ExceptionObject& err )
  83   {
  84     std::cerr
  85       << "Error caught while reading \""
  86       << input_image_fn << "\": " << err
  87       << std::endl;
  88     return( 1 );
  89
  90   } // yrt
  91   TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( );
  92   vtk_input_image->SetInput( input_image );
  93   vtk_input_image->Update( );
  94
  95   // Show input image and let some interaction
  96   fpa::VTK::ImageMPR view;
  97   view.SetBackground( 0.3, 0.2, 0.8 );
  98   view.SetSize( 800, 800 );
  99   view.SetImage( vtk_input_image->GetOutput( ) );
 100
 101   // Allow some interaction and wait for, at least, one seed
 102   view.Render( );
 103   while( view.GetNumberOfSeeds( ) == 0 )
 104     view.Start( );
 105
 106   // Get seed
 107   double p[ 3 ];
 108   view.GetSeed( 0, p );
 109   TImage::PointType pnt;
 110   TImage::IndexType seed;
 111   pnt[ 0 ] = TImage::PointType::ValueType( p[ 0 ] );
 112   pnt[ 1 ] = TImage::PointType::ValueType( p[ 1 ] );
 113   pnt[ 2 ] = TImage::PointType::ValueType( p[ 2 ] );
 114   input_image->TransformPhysicalPointToIndex( pnt, seed );
 115
 116   // Compute squared distance map
 117   TScalarImage::Pointer dmap;
 118   DistanceMap< TImage, TScalarImage >( input_image, dmap );
 119
 120   // Prepare cost conversion function
 121   typedef fpa::Base::Functors::InvertCostFunction< TScalar > TFunction;
 122   TFunction::Pointer function = TFunction::New( );
 123
 124   // Prepare Dijkstra filter
 125   typedef fpa::Image::DijkstraWithEndPointDetection< TScalarImage, TScalarImage > TFilter;
 126   TFilter::Pointer filter = TFilter::New( );
 127   filter->SetInput( dmap );
 128   filter->SetConversionFunction( function );
 129   filter->SetNeighborhoodOrder( 2 );
 130   filter->StopAtOneFrontOff( );
 131   filter->AddSeed( seed, TScalar( 0 ) );
 132
 133   // Prepare graphical debugger
 134   typedef fpa::VTK::Image3DObserver< TFilter, vtkRenderWindow > TDebugger;
 135   TDebugger::Pointer debugger = TDebugger::New( );
 136   debugger->SetRenderWindow( view.GetWindow( ) );
 137   debugger->SetRenderPercentage( 0.0001 );
 138   filter->AddObserver( itk::AnyEvent( ), debugger );
 139   filter->ThrowEventsOn( );
 140
 141   // Go!
 142   std::clock_t start = std::clock( );
 143   filter->Update( );
 144   std::clock_t end = std::clock( );
 145   std::cout
 146     << "Extraction time = "
 147     << ( double( end - start ) / double( CLOCKS_PER_SEC ) )
 148     << " s." << std::endl;
 149
 150   /* TODO
 151   // Save final total cost map
 152   itk::ImageFileWriter< TOutputImage >::Pointer output_image_writer =
 153     itk::ImageFileWriter< TOutputImage >::New( );
 154   output_image_writer->SetFileName( output_image_fn );
 155   output_image_writer->SetInput( filter->GetOutput( ) );
 156   try
 157   {
 158     output_image_writer->Update( );
 159   }
 160   catch( itk::ExceptionObject& err )
 161   {
 162     std::cerr
 163       << "Error while writing image to " << output_image_fn << ": "
 164       << err << std::endl;
 165     return( 1 );
 166
 167   } // yrt
 168   */
 169   return( 0 );
 170 }
 171
 172 // -------------------------------------------------------------------------
 173 template< class I >
 174 void ReadImage( typename I::Pointer& image, const std::string& filename )
 175 {
 176   typename itk::ImageFileReader< I >::Pointer reader =
 177     itk::ImageFileReader< I >::New( );
 178   reader->SetFileName( filename );
 179   reader->Update( );
 180
 181   typename itk::MinimumMaximumImageCalculator< I >::Pointer minmax =
 182     itk::MinimumMaximumImageCalculator< I >::New( );
 183   minmax->SetImage( reader->GetOutput( ) );
 184   minmax->Compute( );
 185   double vmin = double( minmax->GetMinimum( ) );
 186   double vmax = double( minmax->GetMaximum( ) );
 187
 188   typename itk::ShiftScaleImageFilter< I, I >::Pointer shift =
 189     itk::ShiftScaleImageFilter< I, I >::New( );
 190   shift->SetInput( reader->GetOutput( ) );
 191   shift->SetScale( vmax - vmin );
 192   shift->SetShift( vmin / ( vmax - vmin ) );
 193   shift->Update( );
 194
 195   image = shift->GetOutput( );
 196   image->DisconnectPipeline( );
 197 }
 198
 199 // -------------------------------------------------------------------------
 200 template< class I, class O >
 201 void DistanceMap(
 202   const typename I::Pointer& input, typename O::Pointer& output
 203   )
 204 {
 205   typename itk::SignedMaurerDistanceMapImageFilter< I, O >::Pointer dmap =
 206     itk::SignedMaurerDistanceMapImageFilter< I, O >::New( );
 207   dmap->SetInput( input );
 208   dmap->SetBackgroundValue( ( typename I::PixelType )( 0 ) );
 209   dmap->InsideIsPositiveOn( );
 210   dmap->SquaredDistanceOn( );
 211   dmap->UseImageSpacingOn( );
 212
 213   std::clock_t start = std::clock( );
 214   dmap->Update( );
 215   std::clock_t end = std::clock( );
 216   std::cout
 217     << "Distance map time = "
 218     << ( double( end - start ) / double( CLOCKS_PER_SEC ) )
 219     << " s." << std::endl;
 220
 221   output = dmap->GetOutput( );
 222   output->DisconnectPipeline( );
 223 }
 224
 225 // eof - $RCSfile$