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