appli/examples/example_ShowSkeleton.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 <fpa/VTK/ImageMPR.h>
  12
  13 #include <fpa/Image/DijkstraWithEndPointDetection.h>
  14 #include <fpa/IO/MinimumSpanningTreeReader.h>
  15 #include <fpa/IO/UniqueValuesContainerReader.h>
  16 #include <fpa/IO/MatrixValuesContainerReader.h>
  17
  18 #include <vtkCellArray.h>
  19 #include <vtkFloatArray.h>
  20 #include <vtkImageMarchingCubes.h>
  21 #include <vtkPoints.h>
  22 #include <vtkPointData.h>
  23 #include <vtkPolyData.h>
  24 #include <vtkSmartPointer.h>
  25
  26 // -------------------------------------------------------------------------
  27 const unsigned int Dim = 3;
  28 typedef unsigned char TPixel;
  29
  30 typedef itk::Image< TPixel, Dim >            TImage;
  31 typedef itk::ImageToVTKImageFilter< TImage > TVTKInputImage;
  32
  33 typedef fpa::Image::DijkstraWithEndPointDetection< TImage, TImage > TFilter;
  34 typedef TFilter::TMinimumSpanningTree TMinimumSpanningTree;
  35 typedef TFilter::TUniqueVertices TUniqueVertices;
  36 typedef TFilter::TBranches TBranches;
  37 typedef TFilter::TVertices TVertices;
  38
  39 // -------------------------------------------------------------------------
  40 int main( int argc, char* argv[] )
  41 {
  42   if( argc < 6 )
  43   {
  44     std::cerr
  45       << "Usage: " << argv[ 0 ] << std::endl
  46       << " input_image" << std::endl
  47       << " input_minimum_spanning_tree" << std::endl
  48       << " input_endpoints" << std::endl
  49       << " input_bifurcations" << std::endl
  50       << " input_branches" << std::endl
  51       << std::endl;
  52     return( 1 );
  53
  54   } // fi
  55   std::string input_image_fn = argv[ 1 ];
  56   std::string mst_input_fn = argv[ 2 ];
  57   std::string endpoints_input_fn = argv[ 3 ];
  58   std::string bifurcations_input_fn = argv[ 4 ];
  59   std::string branches_input_fn = argv[ 5 ];
  60
  61   // Read image
  62   typename itk::ImageFileReader< TImage >::Pointer input_image_reader =
  63     itk::ImageFileReader< TImage >::New( );
  64   input_image_reader->SetFileName( input_image_fn );
  65   input_image_reader->Update( );
  66   TImage::ConstPointer input_image = input_image_reader->GetOutput( );
  67
  68   TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( );
  69   vtk_input_image->SetInput( input_image );
  70   vtk_input_image->Update( );
  71
  72   // Read minimum spanning tree
  73   fpa::IO::MinimumSpanningTreeReader< TMinimumSpanningTree >::Pointer
  74     mst_input_reader =
  75     fpa::IO::MinimumSpanningTreeReader< TMinimumSpanningTree >::New( );
  76   mst_input_reader->SetFileName( mst_input_fn );
  77   mst_input_reader->Update( );
  78   TMinimumSpanningTree::ConstPointer mst_input = mst_input_reader->GetOutput( );
  79
  80   fpa::IO::UniqueValuesContainerReader< TUniqueVertices >::Pointer
  81     endpoints_input_reader =
  82     fpa::IO::UniqueValuesContainerReader< TUniqueVertices >::New( );
  83   endpoints_input_reader->SetFileName( endpoints_input_fn );
  84   endpoints_input_reader->Update( );
  85   TUniqueVertices::ConstPointer endpoints_input = endpoints_input_reader->GetOutput( );
  86
  87   fpa::IO::UniqueValuesContainerReader< TUniqueVertices >::Pointer
  88     bifurcations_input_reader =
  89     fpa::IO::UniqueValuesContainerReader< TUniqueVertices >::New( );
  90   bifurcations_input_reader->SetFileName( bifurcations_input_fn );
  91   bifurcations_input_reader->Update( );
  92   TUniqueVertices::ConstPointer bifurcations_input =
  93     bifurcations_input_reader->GetOutput( );
  94
  95   fpa::IO::MatrixValuesContainerReader< TBranches >::Pointer
  96     branches_input_reader =
  97     fpa::IO::MatrixValuesContainerReader< TBranches >::New( );
  98   branches_input_reader->SetFileName( branches_input_fn );
  99   branches_input_reader->Update( );
 100   TBranches::ConstPointer branches_input = branches_input_reader->GetOutput( );
 101   unsigned int nBranches = branches_input_reader->GetNumberOfLabels( );
 102
 103   // Show input image and let some interaction
 104   fpa::VTK::ImageMPR view;
 105   view.SetBackground( 0.3, 0.2, 0.8 );
 106   view.SetSize( 800, 800 );
 107   view.SetImage( vtk_input_image->GetOutput( ) );
 108
 109   vtkSmartPointer< vtkImageMarchingCubes > mc =
 110     vtkSmartPointer< vtkImageMarchingCubes >::New( );
 111   mc->SetInputData( vtk_input_image->GetOutput( ) );
 112   mc->SetValue( 0, 1e-1 );
 113   mc->Update( );
 114   view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
 115
 116   // Show endpoints
 117   vtkSmartPointer< vtkPoints > endpoints_points =
 118     vtkSmartPointer< vtkPoints >::New( );
 119   vtkSmartPointer< vtkCellArray > endpoints_cells =
 120     vtkSmartPointer< vtkCellArray >::New( );
 121   for(
 122     TUniqueVertices::ConstIterator epIt = endpoints_input->Begin( );
 123     epIt != endpoints_input->End( );
 124     ++epIt
 125     )
 126   {
 127     TImage::PointType pnt;
 128     input_image->TransformIndexToPhysicalPoint( *epIt, pnt );
 129     endpoints_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 130     endpoints_cells->InsertNextCell( 1 );
 131     endpoints_cells->
 132       InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 1 );
 133
 134   } // rof
 135   vtkSmartPointer< vtkPolyData > endpoints_polydata =
 136     vtkSmartPointer< vtkPolyData >::New( );
 137   endpoints_polydata->SetPoints( endpoints_points );
 138   endpoints_polydata->SetVerts( endpoints_cells );
 139   view.AddPolyData( endpoints_polydata, 0, 0, 1, 1 );
 140
 141   // Show bifurcations
 142   vtkSmartPointer< vtkPoints > bifurcations_points =
 143     vtkSmartPointer< vtkPoints >::New( );
 144   vtkSmartPointer< vtkCellArray > bifurcations_cells =
 145     vtkSmartPointer< vtkCellArray >::New( );
 146   for(
 147     TUniqueVertices::ConstIterator bfIt = bifurcations_input->Begin( );
 148     bfIt != bifurcations_input->End( );
 149     ++bfIt
 150     )
 151   {
 152     TImage::PointType pnt;
 153     input_image->TransformIndexToPhysicalPoint( *bfIt, pnt );
 154     bifurcations_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 155     bifurcations_cells->InsertNextCell( 1 );
 156     bifurcations_cells->
 157       InsertCellPoint( bifurcations_points->GetNumberOfPoints( ) - 1 );
 158
 159   } // rof
 160   vtkSmartPointer< vtkPolyData > bifurcations_polydata =
 161     vtkSmartPointer< vtkPolyData >::New( );
 162   bifurcations_polydata->SetPoints( bifurcations_points );
 163   bifurcations_polydata->SetVerts( bifurcations_cells );
 164   view.AddPolyData( bifurcations_polydata, 0, 1, 0, 1 );
 165
 166   // Show branches (simple and detailed)
 167   vtkSmartPointer< vtkPoints > simple_branches_points =
 168     vtkSmartPointer< vtkPoints >::New( );
 169   vtkSmartPointer< vtkCellArray > simple_branches_cells =
 170     vtkSmartPointer< vtkCellArray >::New( );
 171
 172   vtkSmartPointer< vtkPoints > detailed_branches_points =
 173     vtkSmartPointer< vtkPoints >::New( );
 174   vtkSmartPointer< vtkCellArray > detailed_branches_cells =
 175     vtkSmartPointer< vtkCellArray >::New( );
 176   vtkSmartPointer< vtkFloatArray > detailed_branches_scalars =
 177     vtkSmartPointer< vtkFloatArray >::New( );
 178
 179   TBranches::ConstIterator brIt = branches_input->Begin( );
 180   for( ; brIt != branches_input->End( ); ++brIt )
 181   {
 182     // Branch's first point
 183     TImage::PointType first_point;
 184     input_image->TransformIndexToPhysicalPoint( brIt->first, first_point );
 185     unsigned long first_id = simple_branches_points->GetNumberOfPoints( );
 186     simple_branches_points->InsertNextPoint(
 187       first_point[ 0 ], first_point[ 1 ], first_point[ 2 ]
 188       );
 189
 190     TBranches::ConstRowIterator brRowIt = branches_input->Begin( brIt );
 191     for( ; brRowIt != branches_input->End( brIt ); ++brRowIt )
 192     {
 193       // Branch's second point
 194       TImage::PointType second_point;
 195       input_image->
 196         TransformIndexToPhysicalPoint( brRowIt->first, second_point );
 197       unsigned long second_id = simple_branches_points->GetNumberOfPoints( );
 198       simple_branches_points->InsertNextPoint(
 199         second_point[ 0 ], second_point[ 1 ], second_point[ 2 ]
 200         );
 201       simple_branches_cells->InsertNextCell( 2 );
 202       simple_branches_cells->InsertCellPoint( first_id );
 203       simple_branches_cells->InsertCellPoint( second_id );
 204
 205       // Detailed path
 206       double pathId = double( brRowIt->second - 1 ) / double( nBranches - 1 );
 207       TVertices path;
 208       mst_input->GetPath( path, brIt->first, brRowIt->first );
 209       TVertices::const_iterator pIt = path.begin( );
 210       for( ; pIt != path.end( ); ++pIt )
 211       {
 212         TImage::PointType path_point;
 213         input_image->TransformIndexToPhysicalPoint( *pIt, path_point );
 214         detailed_branches_points->InsertNextPoint(
 215           path_point[ 0 ], path_point[ 1 ], path_point[ 2 ]
 216           );
 217         detailed_branches_scalars->InsertNextTuple1( pathId );
 218         if( pIt != path.begin( ) )
 219         {
 220           unsigned long nPoints =
 221             detailed_branches_points->GetNumberOfPoints( );
 222           detailed_branches_cells->InsertNextCell( 2 );
 223           detailed_branches_cells->InsertCellPoint( nPoints - 2 );
 224           detailed_branches_cells->InsertCellPoint( nPoints - 1 );
 225
 226         } // fi
 227
 228       } // rof
 229
 230     } // rof
 231
 232   } // rof
 233   vtkSmartPointer< vtkPolyData > simple_branches_polydata =
 234     vtkSmartPointer< vtkPolyData >::New( );
 235   simple_branches_polydata->SetPoints( simple_branches_points );
 236   simple_branches_polydata->SetLines( simple_branches_cells );
 237   view.AddPolyData( simple_branches_polydata, 1, 0, 1, 1 );
 238
 239   vtkSmartPointer< vtkPolyData > detailed_branches_polydata =
 240     vtkSmartPointer< vtkPolyData >::New( );
 241   detailed_branches_polydata->SetPoints( detailed_branches_points );
 242   detailed_branches_polydata->SetLines( detailed_branches_cells );
 243   detailed_branches_polydata->
 244     GetPointData( )->SetScalars( detailed_branches_scalars );
 245   view.AddPolyData( detailed_branches_polydata, 1 );
 246
 247   // Allow some interaction
 248   view.Render( );
 249   view.Start( );
 250
 251   return( 0 );
 252 }
 253
 254 // eof - $RCSfile$