appli/examples/example_ImageAlgorithmDijkstra_03.cxx

   1 #include <cmath>
   2 #include <ctime>
   3 #include <deque>
   4 #include <iostream>
   5 #include <iomanip>
   6 #include <limits>
   7 #include <map>
   8 #include <string>
   9
  10 #include <itkConstNeighborhoodIterator.h>
  11 #include <itkNeighborhoodIterator.h>
  12 #include <itkDanielssonDistanceMapImageFilter.h>
  13 #include <itkImage.h>
  14 #include <itkImageFileReader.h>
  15 #include <itkImageFileWriter.h>
  16 #include <itkImageToVTKImageFilter.h>
  17
  18 #include <itkMinimumMaximumImageCalculator.h>
  19 #include <itkInvertIntensityImageFilter.h>
  20
  21 #include <vtkPoints.h>
  22 #include <vtkCellArray.h>
  23 #include <vtkFloatArray.h>
  24 #include <vtkPolyData.h>
  25 #include <vtkPolyDataWriter.h>
  26 #include <vtkSmartPointer.h>
  27 #include <vtkImageMarchingCubes.h>
  28 #include <vtkLookupTable.h>
  29
  30 #include <fpa/Image/DijkstraWithSphereBacktracking.h>
  31 #include <fpa/VTK/ImageMPR.h>
  32 #include <fpa/VTK/Image3DObserver.h>
  33
  34 // -------------------------------------------------------------------------
  35 const unsigned int Dim = 3;
  36 typedef float                                TPixel;
  37 typedef float                                TScalar;
  38 typedef itk::Image< TPixel, Dim >            TImage;
  39 typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
  40
  41 typedef itk::ImageFileReader< TImage >          TImageReader;
  42 typedef itk::ImageFileWriter< TImage >    TImageWriter;
  43 typedef fpa::Image::DijkstraWithSphereBacktracking< TImage, TScalar > TDijkstra;
  44
  45 typedef fpa::VTK::ImageMPR TMPR;
  46 typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs;
  47
  48 struct TBranch
  49 {
  50   double Length;
  51   TImage::PointType V1;
  52   TImage::PointType V2;
  53
  54   bool operator<( const TBranch& other ) const
  55     {
  56       return( other.Length < this->Length );
  57     }
  58 };
  59 typedef std::set< TBranch > TBranches;
  60
  61 // -------------------------------------------------------------------------
  62 int main( int argc, char* argv[] )
  63 {
  64   if( argc < 6 )
  65   {
  66     std::cerr
  67       << "Usage: " << argv[ 0 ]
  68       << " input_image x y z output_image"
  69       << " visual_debug"
  70       << std::endl;
  71     return( 1 );
  72
  73   } // fi
  74   std::string input_image_fn = argv[ 1 ];
  75   TImage::PointType seed_pnt;
  76   seed_pnt[ 0 ] = std::atof( argv[ 2 ] );
  77   seed_pnt[ 1 ] = std::atof( argv[ 3 ] );
  78   seed_pnt[ 2 ] = std::atof( argv[ 4 ] );
  79   std::string output_image_fn = argv[ 5 ];
  80   bool visual_debug = false;
  81   if( argc > 6 )
  82     visual_debug = ( std::atoi( argv[ 6 ] ) == 1 );
  83
  84   // Read image
  85   TImageReader::Pointer input_image_reader = TImageReader::New( );
  86   input_image_reader->SetFileName( input_image_fn );
  87   try
  88   {
  89     input_image_reader->Update( );
  90   }
  91   catch( itk::ExceptionObject& err )
  92   {
  93     std::cerr << "Error caught: " << err << std::endl;
  94     return( 1 );
  95
  96   } // yrt
  97   TImage::ConstPointer input_image = input_image_reader->GetOutput( );
  98
  99   // Show image
 100   TVTKImage::Pointer vtk_image = TVTKImage::New( );
 101   vtk_image->SetInput( input_image );
 102   vtk_image->Update( );
 103
 104   TMPR view;
 105   view.SetBackground( 0.3, 0.2, 0.8 );
 106   view.SetSize( 800, 800 );
 107   view.SetImage( vtk_image->GetOutput( ) );
 108
 109   /*
 110   vtkSmartPointer< vtkImageMarchingCubes > mc =
 111     vtkSmartPointer< vtkImageMarchingCubes >::New( );
 112   mc->SetInputData( vtk_image->GetOutput( ) );
 113   mc->SetValue( 0, 1e-2 );
 114   mc->Update( );
 115   view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
 116   */
 117
 118   // Allow some interaction
 119   view.Render( );
 120   view.Start( );
 121
 122   TImage::IndexType seed_idx;
 123   input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx );
 124   std::cout << seed_idx << " " << seed_pnt << std::endl;
 125   seed_idx[ 0 ] = 256;
 126   seed_idx[ 1 ] = 313;
 127   seed_idx[ 2 ] = 381;
 128
 129   // Extract paths
 130   TDijkstra::Pointer paths = TDijkstra::New( );
 131   paths->AddSeed( seed_idx, TScalar( 0 ) );
 132   paths->SetInput( input_image );
 133   paths->SetNeighborhoodOrder( 2 );
 134
 135   if( visual_debug )
 136   {
 137     // Configure observer
 138     TDijkstraObs::Pointer obs = TDijkstraObs::New( );
 139     obs->SetRenderWindow( view.GetWindow( ) );
 140     paths->AddObserver( itk::AnyEvent( ), obs );
 141     paths->ThrowEventsOn( );
 142   }
 143   else
 144     paths->ThrowEventsOff( );
 145   std::clock_t start = std::clock( );
 146   paths->Update( );
 147   std::clock_t end = std::clock( );
 148   double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
 149   std::cout << "Paths extraction time = " << seconds << std::endl;
 150
 151   // Create polydata
 152   vtkSmartPointer< vtkPoints > points =
 153     vtkSmartPointer< vtkPoints >::New( );
 154   vtkSmartPointer< vtkCellArray > cells =
 155     vtkSmartPointer< vtkCellArray >::New( );
 156   vtkSmartPointer< vtkFloatArray > scalars =
 157     vtkSmartPointer< vtkFloatArray >::New( );
 158
 159   const TDijkstra::TMarkImage* marks = paths->GetOutputMarkImage( );
 160   TDijkstra::TMark max_mark = paths->GetNumberOfBranches( );
 161   const TDijkstra::TVertices& endpoints = paths->GetEndPoints( );
 162   const TDijkstra::TVertices& bifurcations = paths->GetBifurcationPoints( );
 163   const TDijkstra::TTree& tree = paths->GetFinalTree( );
 164   TDijkstra::TVertices::const_iterator epIt = endpoints.begin( );
 165
 166   TDijkstra::TMarkImage::Pointer new_marks = TDijkstra::TMarkImage::New( );
 167   new_marks->SetLargestPossibleRegion( marks->GetLargestPossibleRegion( ) );
 168   new_marks->SetRequestedRegion( marks->GetRequestedRegion( ) );
 169   new_marks->SetBufferedRegion( marks->GetBufferedRegion( ) );
 170   new_marks->SetOrigin( marks->GetOrigin( ) );
 171   new_marks->SetSpacing( marks->GetSpacing( ) );
 172   new_marks->SetDirection( marks->GetDirection( ) );
 173   new_marks->Allocate( );
 174   new_marks->FillBuffer( 0 );
 175
 176   std::cout << max_mark << std::endl;
 177   std::cout << endpoints.size( ) << std::endl;
 178
 179   TBranches branches;
 180   TBranch branch;
 181   for( ; epIt != endpoints.end( ); ++epIt )
 182   {
 183     TDijkstra::TVertex idx = *epIt;
 184     TDijkstra::TTree::const_iterator tIt = tree.find( idx );
 185     TImage::PointType pnt, prev;
 186     input_image->TransformIndexToPhysicalPoint( idx, pnt );
 187     do
 188     {
 189       input_image->TransformIndexToPhysicalPoint( idx, pnt );
 190       branch.V2 = pnt;
 191
 192       points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 193       scalars->InsertNextTuple1( double( tIt->second.second ) );
 194       new_marks->SetPixel( idx, tIt->second.second );
 195
 196       if( idx != *epIt )
 197       {
 198         cells->InsertNextCell( 2 );
 199         cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
 200         cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
 201         branch.Length += prev.EuclideanDistanceTo( pnt );
 202
 203       } // fi
 204       prev = pnt;
 205       idx = tIt->second.first;
 206
 207       if( std::find( bifurcations.begin( ), bifurcations.end( ), idx ) != bifurcations.end( ) )
 208       {
 209         branches.insert( branch );
 210         prev = pnt;
 211         branch.V1 = pnt;
 212         branch.Length = double( 0 );
 213       }
 214
 215       tIt = tree.find( idx );
 216
 217     } while( idx != tIt->second.first );
 218
 219
 220   } // rof
 221
 222   /* TODO
 223      int i = 1;
 224      TImage::PointType ori = input_image->GetOrigin( );
 225      std::cout << ori << std::endl;
 226      for( TBranches::const_iterator bIt = branches.begin( ); bIt != branches.end( ); ++bIt, ++i )
 227      {
 228      TImage::PointType::VectorType v1 = bIt->V1 - ori;
 229      TImage::PointType::VectorType v2 = bIt->V2 - ori;
 230      std::cout
 231      << std::setiosflags( std::ios::fixed) << std::setprecision( 3 )
 232      << i << "\t1.000\t"
 233      << bIt->Length << "\t"
 234      << v1[ 0 ] << "\t"
 235      << v1[ 1 ] << "\t"
 236      << v1[ 2 ] << "\t"
 237      << v2[ 0 ] << "\t"
 238      << v2[ 1 ] << "\t"
 239      << v2[ 2 ] << "\t"
 240      << ( v2 - v1 ).GetNorm( )
 241      << std::endl;
 242
 243      } // rof
 244   */
 245
 246   vtkSmartPointer< vtkPolyData > vtk_tree =
 247     vtkSmartPointer< vtkPolyData >::New( );
 248   vtk_tree->SetPoints( points );
 249   vtk_tree->SetLines( cells );
 250   vtk_tree->GetPointData( )->SetScalars( scalars );
 251
 252   view.AddPolyData( vtk_tree );
 253   view.Render( );
 254   view.Start( );
 255
 256   vtkSmartPointer< vtkPolyDataWriter > writer =
 257     vtkSmartPointer< vtkPolyDataWriter >::New( );
 258   writer->SetInputData( vtk_tree );
 259   writer->SetFileName( output_image_fn.c_str( ) );
 260   writer->Update( );
 261
 262   itk::ImageFileWriter< TDijkstra::TMarkImage >::Pointer marks_writer =
 263     itk::ImageFileWriter< TDijkstra::TMarkImage >::New( );
 264   marks_writer->SetInput( new_marks );
 265   marks_writer->SetFileName( "marks.mhd" );
 266   marks_writer->Update( );
 267
 268   /* TODO
 269      TImageWriter::Pointer dijkstra_writer =
 270      TImageWriter::New( );
 271      dijkstra_writer->SetInput( paths->GetOutput( ) );
 272      dijkstra_writer->SetFileName( "dijkstra.mhd" );
 273      dijkstra_writer->Update( );
 274   */
 275
 276   return( 0 );
 277 }
 278
 279 // eof - $RCSfile$