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