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 <itkSignedMaurerDistanceMapImageFilter.h>
  11
  12 #include <fpa/Base/Functors/InvertCostFunction.h>
  13 #include <fpa/Image/DijkstraWithEndPointDetection.h>
  14 #include <fpa/VTK/Image2DObserver.h>
  15 #include <fpa/IO/MinimumSpanningTreeWriter.h>
  16 #include <fpa/IO/UniqueValuesContainerWriter.h>
  17 #include <fpa/IO/MatrixValuesContainerWriter.h>
  18 #include <fpa/VTK/UniqueVerticesToPolyDataFilter.h>
  19
  20 #include "fpa_Utility.h"
  21
  22 // -------------------------------------------------------------------------
  23 const unsigned int Dim = 2;
  24 typedef unsigned char TPixel;
  25 typedef float         TScalar;
  26 typedef itk::Image< TPixel, Dim >  TImage;
  27 typedef itk::Image< TScalar, Dim > TScalarImage;
  28
  29 // -------------------------------------------------------------------------
  30 int main( int argc, char* argv[] )
  31 {
  32   if( argc < 7 )
  33   {
  34     std::cerr
  35       << "Usage: " << argv[ 0 ] << std::endl
  36       << "\tinput_image" << std::endl
  37       << "\toutput_labels" << std::endl
  38       << "\toutput_minimum_spanning_tree" << std::endl
  39       << "\toutput_endpoints" << std::endl
  40       << "\toutput_bifurcations" << std::endl
  41       << "\toutput_branches" << std::endl
  42       << std::endl;
  43     return( 1 );
  44
  45   } // fi
  46   std::string input_image_fn = argv[ 1 ];
  47   std::string output_labels_fn = argv[ 2 ];
  48   std::string output_minimum_spanning_tree_fn = argv[ 3 ];
  49   std::string output_endpoints_fn = argv[ 4 ];
  50   std::string output_bifurcations_fn = argv[ 5 ];
  51   std::string output_branches_fn = argv[ 6 ];
  52
  53   // Read image
  54   TImage::Pointer input_image;
  55   std::string err = fpa_Utility::ReadImage( input_image, input_image_fn );
  56   if( err != "" )
  57   {
  58     std::cerr << err << std::endl;
  59     return( 1 );
  60
  61   } // fi
  62
  63   // Show image and wait for, at least, one seed
  64   itk::ImageToVTKImageFilter< TImage >::Pointer vtk_input_image =
  65     itk::ImageToVTKImageFilter< TImage >::New( );
  66   vtk_input_image->SetInput( input_image );
  67   vtk_input_image->Update( );
  68
  69   fpa_Utility::Viewer2DWithSeeds viewer;
  70   viewer.SetImage( vtk_input_image->GetOutput( ) );
  71   while( viewer.GetNumberOfSeeds( ) == 0 )
  72     viewer.Start( );
  73
  74   // Compute squared distance map
  75   itk::SignedMaurerDistanceMapImageFilter< TImage, TScalarImage >::Pointer
  76     dmap =
  77     itk::SignedMaurerDistanceMapImageFilter< TImage, TScalarImage >::New( );
  78   dmap->SetInput( input_image );
  79   dmap->SetBackgroundValue( TPixel( 0 ) );
  80   dmap->InsideIsPositiveOn( );
  81   dmap->SquaredDistanceOn( );
  82   dmap->UseImageSpacingOn( );
  83
  84   std::time_t start, end;
  85   std::time( &start );
  86   dmap->Update( );
  87   std::time( &end );
  88   std::cout
  89     << "Distance map time = "
  90     << std::difftime( end, start )
  91     << "s." << std::endl;
  92
  93   // Prepare cost conversion function
  94   typedef fpa::Base::Functors::InvertCostFunction< TScalar > TCostFunction;
  95   TCostFunction::Pointer cost_f = TCostFunction::New( );
  96
  97   // Prepare Dijkstra filter
  98   typedef fpa::Image::
  99     DijkstraWithEndPointDetection< TScalarImage, TScalarImage > TFilter;
 100   TFilter::Pointer filter = TFilter::New( );
 101   filter->SetInput( dmap->GetOutput( ) );
 102   filter->SetConversionFunction( cost_f );
 103   filter->SetNeighborhoodOrder( 2 );
 104   filter->StopAtOneFrontOff( );
 105
 106   // Associate seed
 107   TImage::PointType pnt;
 108   TImage::IndexType idx;
 109   viewer.GetSeed( pnt, 0 );
 110   if( input_image->TransformPhysicalPointToIndex( pnt, idx ) )
 111     filter->AddSeed( idx, 0 );
 112
 113   // Prepare graphical debugger
 114   typedef fpa::VTK::Image2DObserver< TFilter, vtkRenderWindow > TDebugger;
 115   TDebugger::Pointer debugger = TDebugger::New( );
 116   debugger->SetRenderWindow( viewer.Window );
 117   debugger->SetRenderPercentage( 0.01 );
 118   filter->AddObserver( itk::AnyEvent( ), debugger );
 119   filter->ThrowEventsOn( );
 120
 121   // Go!
 122   std::time( &start );
 123   filter->Update( );
 124   std::time( &end );
 125   std::cout
 126     << "Extraction time = "
 127     << std::difftime( end, start )
 128     << "s." << std::endl;
 129
 130   // Create new actors
 131   typedef vtkSmartPointer< fpa::VTK::UniqueVerticesToPolyDataFilter< TFilter::TUniqueVertices, TFilter::TInputImage > > TVertices2PD;
 132
 133   TVertices2PD endpoints2pd = TVertices2PD::New( );
 134   endpoints2pd->SetInput( filter->GetEndPoints( ) );
 135   endpoints2pd->SetImage( filter->GetInput( ) );
 136   endpoints2pd->Update( );
 137
 138   vtkSmartPointer< vtkPolyDataMapper > endpoints_mapper =
 139     vtkSmartPointer< vtkPolyDataMapper >::New( );
 140   endpoints_mapper->SetInputConnection( endpoints2pd->GetOutputPort( ) );
 141
 142   vtkSmartPointer< vtkActor > endpoints_actor =
 143     vtkSmartPointer< vtkActor >::New( );
 144   endpoints_actor->SetMapper( endpoints_mapper );
 145   endpoints_actor->GetProperty( )->SetColor( 0, 1, 0 );
 146   endpoints_actor->GetProperty( )->SetPointSize( 5 );
 147   viewer.Renderer->AddActor( endpoints_actor );
 148
 149   TVertices2PD bifurcations2pd = TVertices2PD::New( );
 150   bifurcations2pd->SetInput( filter->GetBifurcations( ) );
 151   bifurcations2pd->SetImage( filter->GetInput( ) );
 152   bifurcations2pd->Update( );
 153
 154   vtkSmartPointer< vtkPolyDataMapper > bifurcations_mapper =
 155     vtkSmartPointer< vtkPolyDataMapper >::New( );
 156   bifurcations_mapper->SetInputConnection( bifurcations2pd->GetOutputPort( ) );
 157
 158   vtkSmartPointer< vtkActor > bifurcations_actor =
 159     vtkSmartPointer< vtkActor >::New( );
 160   bifurcations_actor->SetMapper( bifurcations_mapper );
 161   bifurcations_actor->GetProperty( )->SetColor( 1, 0, 0 );
 162   bifurcations_actor->GetProperty( )->SetPointSize( 5 );
 163   viewer.Renderer->AddActor( bifurcations_actor );
 164
 165   // Some more interaction and finish
 166   viewer.Render( );
 167   viewer.Start( );
 168
 169   // Save results
 170   err = fpa_Utility::SaveImage( filter->GetLabelImage( ), output_labels_fn );
 171   if( err != "" )
 172     std::cerr << err << std::endl;
 173
 174   fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::Pointer
 175     mst_writer =
 176     fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::New( );
 177   mst_writer->SetInput( filter->GetMinimumSpanningTree( ) );
 178   mst_writer->SetFileName( output_minimum_spanning_tree_fn );
 179   mst_writer->Update( );
 180
 181   fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer
 182     endpoints_writer =
 183     fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( );
 184   endpoints_writer->SetInput( filter->GetEndPoints( ) );
 185   endpoints_writer->SetFileName( output_endpoints_fn );
 186   endpoints_writer->Update( );
 187
 188   fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer
 189     bifurcations_writer =
 190     fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( );
 191   bifurcations_writer->SetInput( filter->GetBifurcations( ) );
 192   bifurcations_writer->SetFileName( output_bifurcations_fn );
 193   bifurcations_writer->Update( );
 194
 195   fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::Pointer
 196     branches_writer =
 197     fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::New( );
 198   branches_writer->SetInput( filter->GetBranches( ) );
 199   branches_writer->SetFileName( output_branches_fn );
 200   branches_writer->Update( );
 201
 202   return( 0 );
 203 }
 204
 205 // eof - $RCSfile$
 206
 207 /* TODO
 208    vtkSmartPointer< vtkPoints > simple_branches_points =
 209    vtkSmartPointer< vtkPoints >::New( );
 210    vtkSmartPointer< vtkCellArray > simple_branches_cells =
 211    vtkSmartPointer< vtkCellArray >::New( );
 212
 213    vtkSmartPointer< vtkPoints > detailed_branches_points =
 214    vtkSmartPointer< vtkPoints >::New( );
 215    vtkSmartPointer< vtkCellArray > detailed_branches_cells =
 216    vtkSmartPointer< vtkCellArray >::New( );
 217    vtkSmartPointer< vtkFloatArray > detailed_branches_scalars =
 218    vtkSmartPointer< vtkFloatArray >::New( );
 219
 220    TFilter::TBranches::ConstIterator brIt = branches->Begin( );
 221    for( ; brIt != branches->End( ); ++brIt )
 222    {
 223    // Branch's first point
 224    TImage::PointType first_point;
 225    input_image->TransformIndexToPhysicalPoint( brIt->first, first_point );
 226    unsigned long first_id = simple_branches_points->GetNumberOfPoints( );
 227    simple_branches_points->InsertNextPoint(
 228    first_point[ 0 ], first_point[ 1 ], first_point[ 2 ]
 229    );
 230
 231    TFilter::TBranches::ConstRowIterator brRowIt = branches->Begin( brIt );
 232    for( ; brRowIt != branches->End( brIt ); ++brRowIt )
 233    {
 234    // Branch's second point
 235    TImage::PointType second_point;
 236    input_image->
 237    TransformIndexToPhysicalPoint( brRowIt->first, second_point );
 238    unsigned long second_id = simple_branches_points->GetNumberOfPoints( );
 239    simple_branches_points->InsertNextPoint(
 240    second_point[ 0 ], second_point[ 1 ], second_point[ 2 ]
 241    );
 242    simple_branches_cells->InsertNextCell( 2 );
 243    simple_branches_cells->InsertCellPoint( first_id );
 244    simple_branches_cells->InsertCellPoint( second_id );
 245
 246    // Detailed path
 247    double pathId = double( brRowIt->second - 1 ) / double( nBranches - 1 );
 248    TFilter::TVertices path;
 249    mst->GetPath( path, brIt->first, brRowIt->first );
 250    TFilter::TVertices::const_iterator pIt = path.begin( );
 251    for( ; pIt != path.end( ); ++pIt )
 252    {
 253    TImage::PointType path_point;
 254    input_image->TransformIndexToPhysicalPoint( *pIt, path_point );
 255    detailed_branches_points->InsertNextPoint(
 256    path_point[ 0 ], path_point[ 1 ], path_point[ 2 ]
 257    );
 258    detailed_branches_scalars->InsertNextTuple1( pathId );
 259    if( pIt != path.begin( ) )
 260    {
 261    unsigned long nPoints =
 262    detailed_branches_points->GetNumberOfPoints( );
 263    detailed_branches_cells->InsertNextCell( 2 );
 264    detailed_branches_cells->InsertCellPoint( nPoints - 2 );
 265    detailed_branches_cells->InsertCellPoint( nPoints - 1 );
 266
 267    } // fi
 268
 269    } // rof
 270
 271    } // rof
 272
 273    } // rof
 274    vtkSmartPointer< vtkPolyData > simple_branches_polydata =
 275    vtkSmartPointer< vtkPolyData >::New( );
 276    simple_branches_polydata->SetPoints( simple_branches_points );
 277    simple_branches_polydata->SetLines( simple_branches_cells );
 278    view.AddPolyData( simple_branches_polydata, 1, 0, 1, 1 );
 279
 280    vtkSmartPointer< vtkPolyData > detailed_branches_polydata =
 281    vtkSmartPointer< vtkPolyData >::New( );
 282    detailed_branches_polydata->SetPoints( detailed_branches_points );
 283    detailed_branches_polydata->SetLines( detailed_branches_cells );
 284    detailed_branches_polydata->
 285    GetPointData( )->SetScalars( detailed_branches_scalars );
 286    view.AddPolyData( detailed_branches_polydata, 1 );
 287
 288    // Let some more interaction
 289    view.Render( );
 290    view.Start( );
 291
 292    return( 0 );
 293    }
 294 */
 295
 296 // eof - $RCSfile$