appli/examples/example_Image_Dijkstra_DanielssonCost_TwoSeedsPath.cxx

   1 #include <cstdlib>
   2 #include <iostream>
   3 #include <limits>
   4 #include <string>
   5
   6 #include <itkImage.h>
   7 #include <itkRGBAPixel.h>
   8 #include <itkImageToVTKImageFilter.h>
   9
  10 #include <itkImageFileReader.h>
  11 #include <itkMinimumMaximumImageCalculator.h>
  12 #include <itkInvertIntensityImageFilter.h>
  13 #include <itkDanielssonDistanceMapImageFilter.h>
  14 #include <itkImageFileWriter.h>
  15
  16 #include <vtkCamera.h>
  17 #include <vtkImageActor.h>
  18 #include <vtkInteractorStyleImage.h>
  19 #include <vtkPointHandleRepresentation3D.h>
  20 #include <vtkProperty.h>
  21 #include <vtkRenderer.h>
  22 #include <vtkRenderWindow.h>
  23 #include <vtkRenderWindowInteractor.h>
  24 #include <vtkSeedRepresentation.h>
  25 #include <vtkSeedWidget.h>
  26 #include <vtkSmartPointer.h>
  27
  28 #include <fpa/Image/Dijkstra.h>
  29 #include <fpa/Base/Functors/InvertCostFunction.h>
  30 #include <fpa/VTK/Image2DObserver.h>
  31
  32 // -------------------------------------------------------------------------
  33 const unsigned int Dim = 2;
  34 typedef unsigned char TInputPixel;
  35 typedef float         TOutputPixel;
  36
  37 typedef itk::Image< TInputPixel, Dim >            TInputImage;
  38 typedef itk::Image< TOutputPixel, Dim >           TOutputImage;
  39 typedef itk::ImageToVTKImageFilter< TInputImage > TVTKInputImage;
  40
  41 // -------------------------------------------------------------------------
  42 int main( int argc, char* argv[] )
  43 {
  44   if( argc < 5 )
  45   {
  46     std::cerr
  47       << "Usage: " << argv[ 0 ]
  48       << " input_image output_image neighborhood_order stop_at_one_front"
  49       << std::endl;
  50     return( 1 );
  51
  52   } // fi
  53   std::string input_image_fn = argv[ 1 ];
  54   std::string output_image_fn = argv[ 2 ];
  55   unsigned int neighborhood_order = std::atoi( argv[ 3 ] );
  56   bool stop_at_one_front = ( std::atoi( argv[ 4 ] ) != 0 );
  57
  58   // Read image
  59   itk::ImageFileReader< TInputImage >::Pointer input_image_reader =
  60     itk::ImageFileReader< TInputImage >::New( );
  61   input_image_reader->SetFileName( input_image_fn );
  62   try
  63   {
  64     input_image_reader->Update( );
  65   }
  66   catch( itk::ExceptionObject& err )
  67   {
  68     std::cerr
  69       << "Error while reading image from " << input_image_fn << ": "
  70       << err << std::endl;
  71     return( 1 );
  72
  73   } // yrt
  74   TInputImage::ConstPointer input_image = input_image_reader->GetOutput( );
  75   TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( );
  76   vtk_input_image->SetInput( input_image );
  77   vtk_input_image->Update( );
  78
  79   // VTK visualization
  80   vtkSmartPointer< vtkImageActor > actor =
  81     vtkSmartPointer< vtkImageActor >::New( );
  82   actor->SetInputData( vtk_input_image->GetOutput( ) );
  83
  84   vtkSmartPointer< vtkRenderer > renderer =
  85     vtkSmartPointer< vtkRenderer >::New( );
  86   renderer->SetBackground( 0.1, 0.2, 0.7 );
  87   renderer->AddActor( actor );
  88   vtkSmartPointer< vtkRenderWindow > window =
  89     vtkSmartPointer< vtkRenderWindow >::New( );
  90   window->SetSize( 800, 800 );
  91   window->AddRenderer( renderer );
  92
  93   // Correct camera due to the loaded image
  94   vtkCamera* camera = renderer->GetActiveCamera( );
  95   camera->SetViewUp( 0, -1, 0 );
  96   camera->SetPosition( 0, 0, -1 );
  97   camera->SetFocalPoint( 0, 0, 0 );
  98
  99   // VTK interaction
 100   vtkSmartPointer< vtkInteractorStyleImage > imageStyle =
 101     vtkSmartPointer< vtkInteractorStyleImage >::New( );
 102   vtkSmartPointer< vtkRenderWindowInteractor > interactor =
 103     vtkSmartPointer< vtkRenderWindowInteractor >::New( );
 104   interactor->SetInteractorStyle( imageStyle );
 105   window->SetInteractor( interactor );
 106
 107   // Create the widget and its representation
 108   vtkSmartPointer< vtkPointHandleRepresentation3D > handle =
 109     vtkSmartPointer< vtkPointHandleRepresentation3D >::New( );
 110   handle->GetProperty( )->SetColor( 1, 0, 0 );
 111   vtkSmartPointer< vtkSeedRepresentation > rep =
 112     vtkSmartPointer< vtkSeedRepresentation >::New( );
 113   rep->SetHandleRepresentation( handle );
 114
 115   vtkSmartPointer< vtkSeedWidget > widget =
 116     vtkSmartPointer< vtkSeedWidget >::New( );
 117   widget->SetInteractor( interactor );
 118   widget->SetRepresentation( rep );
 119
 120   // Let some interaction
 121   interactor->Initialize( );
 122   renderer->ResetCamera( );
 123   window->Render( );
 124   widget->On( );
 125   interactor->Start( );
 126
 127   // Invert input image
 128   itk::MinimumMaximumImageCalculator< TInputImage >::Pointer minmax =
 129     itk::MinimumMaximumImageCalculator< TInputImage >::New( );
 130   minmax->SetImage( input_image );
 131   minmax->Compute( );
 132
 133   itk::InvertIntensityImageFilter< TInputImage >::Pointer invert =
 134     itk::InvertIntensityImageFilter< TInputImage >::New( );
 135   invert->SetInput( input_image );
 136   invert->SetMaximum( minmax->GetMaximum( ) );
 137
 138   itk::DanielssonDistanceMapImageFilter< TInputImage, TOutputImage >::Pointer dmap =
 139     itk::DanielssonDistanceMapImageFilter< TInputImage, TOutputImage >::New( );
 140   dmap->SetInput( invert->GetOutput( ) );
 141   dmap->InputIsBinaryOn( );
 142   dmap->SquaredDistanceOn( );
 143   dmap->UseImageSpacingOn( );
 144   dmap->Update( );
 145
 146   typedef fpa::Base::Functors::InvertCostFunction< TOutputPixel > TFunction;
 147   TFunction::Pointer function = TFunction::New( );
 148
 149   // Prepare region grow filter
 150   typedef fpa::Image::Dijkstra< TOutputImage, TOutputImage > TFilter;
 151   TFilter::Pointer filter = TFilter::New( );
 152   filter->SetInput( dmap->GetOutput( ) );
 153   filter->SetConversionFunction( function );
 154   filter->SetNeighborhoodOrder( neighborhood_order );
 155   filter->SetStopAtOneFront( stop_at_one_front );
 156
 157   // Get user-given seeds
 158   for( unsigned int s = 0; s < rep->GetNumberOfSeeds( ); s++ )
 159   {
 160     double pos[ 3 ];
 161     rep->GetSeedWorldPosition( s, pos );
 162
 163     TInputImage::PointType pnt;
 164     pnt[ 0 ] = TInputImage::PointType::ValueType( pos[ 0 ] );
 165     pnt[ 1 ] = TInputImage::PointType::ValueType( pos[ 1 ] );
 166
 167     TInputImage::IndexType idx;
 168     if( input_image->TransformPhysicalPointToIndex( pnt, idx ) )
 169       filter->AddSeed( idx, 0 );
 170
 171   } // rof
 172
 173   // Prepare graphical debugger
 174   typedef fpa::VTK::Image2DObserver< TFilter, vtkRenderWindow > TDebugger;
 175   TDebugger::Pointer debugger = TDebugger::New( );
 176   debugger->SetRenderWindow( window );
 177   debugger->SetRenderPercentage( 0.001 );
 178   filter->AddObserver( itk::AnyEvent( ), debugger );
 179   filter->ThrowEventsOn( );
 180
 181   // Go!
 182   filter->Update( );
 183
 184   // Save final total cost map
 185   itk::ImageFileWriter< TOutputImage >::Pointer output_image_writer =
 186     itk::ImageFileWriter< TOutputImage >::New( );
 187   output_image_writer->SetFileName( output_image_fn );
 188   output_image_writer->SetInput( filter->GetOutput( ) );
 189   try
 190   {
 191     output_image_writer->Update( );
 192   }
 193   catch( itk::ExceptionObject& err )
 194   {
 195     std::cerr
 196       << "Error while writing image to " << output_image_fn << ": "
 197       << err << std::endl;
 198     return( 1 );
 199
 200   } // yrt
 201
 202   // Draw the path between the two seeds
 203   if( filter->GetNumberOfSeeds( ) > 1 )
 204   {
 205     TInputImage::IndexType a = filter->GetSeed( 0 );
 206     TInputImage::IndexType b =
 207       filter->GetSeed( filter->GetNumberOfSeeds( ) - 1 );
 208
 209     std::vector< TInputImage::IndexType > path;
 210     filter->GetMinimumSpanningTree( )->GetPath( path, a, b );
 211
 212     typedef itk::Image< itk::RGBAPixel< unsigned char >, Dim > TColorImage;
 213     TColorImage::Pointer path_image = TColorImage::New( );
 214     path_image->SetLargestPossibleRegion(
 215       input_image->GetLargestPossibleRegion( )
 216       );
 217     path_image->SetRequestedRegion( input_image->GetRequestedRegion( ) );
 218     path_image->SetBufferedRegion( input_image->GetBufferedRegion( ) );
 219     path_image->SetSpacing( input_image->GetSpacing( ) );
 220     path_image->SetOrigin( input_image->GetOrigin( ) );
 221     path_image->SetDirection( input_image->GetDirection( ) );
 222     path_image->Allocate( );
 223
 224     unsigned char background[ 4 ] = { 0, 0, 0, 0 };
 225     unsigned char color[ 4 ] = { 255, 0, 0, 255 };
 226     path_image->FillBuffer( itk::RGBAPixel< unsigned char >( background ) );
 227
 228     for( unsigned int pId = 0; pId < path.size( ); ++pId )
 229       path_image->SetPixel(
 230         path[ pId ], itk::RGBAPixel< unsigned char >( color )
 231         );
 232
 233     itk::ImageToVTKImageFilter< TColorImage >::Pointer vtk_path_image =
 234       itk::ImageToVTKImageFilter< TColorImage >::New( );
 235     vtk_path_image->SetInput( path_image );
 236     vtk_path_image->Update( );
 237
 238     vtkSmartPointer< vtkImageActor > path_image_actor =
 239       vtkSmartPointer< vtkImageActor >::New( );
 240     path_image_actor->SetInputData( vtk_path_image->GetOutput( ) );
 241     renderer->AddActor( path_image_actor );
 242
 243     // More interaction
 244     window->Render( );
 245     widget->Off( );
 246     interactor->Start( );
 247
 248   } // fi
 249
 250   return( 0 );
 251 }
 252
 253 // eof - $RCSfile$