9 #include <itkConstNeighborhoodIterator.h>
10 #include <itkNeighborhoodIterator.h>
11 #include <itkDanielssonDistanceMapImageFilter.h>
13 #include <itkImageFileReader.h>
14 #include <itkImageFileWriter.h>
15 #include <itkImageToVTKImageFilter.h>
17 #include <itkMinimumMaximumImageCalculator.h>
18 #include <itkInvertIntensityImageFilter.h>
20 #include <vtkPoints.h>
21 #include <vtkCellArray.h>
22 #include <vtkFloatArray.h>
23 #include <vtkPolyData.h>
24 #include <vtkSmartPointer.h>
25 #include <vtkImageMarchingCubes.h>
27 #include <fpa/Image/DijkstraWithSphereBacktracking.h>
28 #include <fpa/VTK/ImageMPR.h>
29 #include <fpa/VTK/Image3DObserver.h>
31 // -------------------------------------------------------------------------
32 const unsigned int Dim = 3;
34 typedef float TScalar;
35 typedef itk::Image< TPixel, Dim > TImage;
36 typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
38 typedef itk::ImageFileReader< TImage > TImageReader;
39 typedef itk::ImageFileWriter< TImage > TImageWriter;
40 typedef fpa::Image::DijkstraWithSphereBacktracking< TImage, TScalar > TDijkstra;
42 typedef fpa::VTK::ImageMPR TMPR;
43 typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs;
45 // -------------------------------------------------------------------------
46 int main( int argc, char* argv[] )
51 << "Usage: " << argv[ 0 ]
52 << " input_image x y z output_image"
58 std::string input_image_fn = argv[ 1 ];
59 TImage::PointType seed_pnt;
60 seed_pnt[ 0 ] = std::atof( argv[ 2 ] );
61 seed_pnt[ 1 ] = std::atof( argv[ 3 ] );
62 seed_pnt[ 2 ] = std::atof( argv[ 4 ] );
63 std::string output_image_fn = argv[ 5 ];
64 bool visual_debug = false;
66 visual_debug = ( std::atoi( argv[ 6 ] ) == 1 );
69 TImageReader::Pointer input_image_reader = TImageReader::New( );
70 input_image_reader->SetFileName( input_image_fn );
73 input_image_reader->Update( );
75 catch( itk::ExceptionObject& err )
77 std::cerr << "Error caught: " << err << std::endl;
81 TImage::ConstPointer input_image = input_image_reader->GetOutput( );
84 TVTKImage::Pointer vtk_image = TVTKImage::New( );
85 vtk_image->SetInput( input_image );
89 view.SetBackground( 0.3, 0.2, 0.8 );
90 view.SetSize( 800, 800 );
91 view.SetImage( vtk_image->GetOutput( ) );
93 vtkSmartPointer< vtkImageMarchingCubes > mc =
94 vtkSmartPointer< vtkImageMarchingCubes >::New( );
95 mc->SetInputData( vtk_image->GetOutput( ) );
96 mc->SetValue( 0, 1e-2 );
98 view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
100 // Allow some interaction
104 TImage::IndexType seed_idx;
105 input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx );
106 std::cout << seed_idx << " " << seed_pnt << std::endl;
112 TDijkstra::Pointer paths = TDijkstra::New( );
113 paths->AddSeed( seed_idx, TScalar( 0 ) );
114 paths->SetInput( input_image );
115 paths->SetNeighborhoodOrder( 2 );
119 // Configure observer
120 TDijkstraObs::Pointer obs = TDijkstraObs::New( );
121 obs->SetRenderWindow( view.GetWindow( ) );
122 paths->AddObserver( itk::AnyEvent( ), obs );
123 paths->ThrowEventsOn( );
126 paths->ThrowEventsOff( );
127 std::clock_t start = std::clock( );
129 std::clock_t end = std::clock( );
130 double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
131 std::cout << "Paths extraction time = " << seconds << std::endl;
134 vtkSmartPointer< vtkPoints > points =
135 vtkSmartPointer< vtkPoints >::New( );
136 vtkSmartPointer< vtkCellArray > cells =
137 vtkSmartPointer< vtkCellArray >::New( );
138 vtkSmartPointer< vtkFloatArray > scalars =
139 vtkSmartPointer< vtkFloatArray >::New( );
141 const TDijkstra::TVertices& endpoints = paths->GetEndPoints( );
142 const TDijkstra::TTree& tree = paths->GetFinalTree( );
143 TDijkstra::TVertices::const_iterator epIt = endpoints.begin( );
144 for( unsigned int epId = 0; epIt != endpoints.end( ); ++epIt, ++epId )
146 double pd = double( epId ) / double( endpoints.size( ) - 1 );
148 TDijkstra::TVertex idx = *epIt;
151 TImage::PointType pnt;
152 input_image->TransformIndexToPhysicalPoint( idx, pnt );
154 points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
155 scalars->InsertNextTuple1( pd );
158 cells->InsertNextCell( 2 );
159 cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
160 cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
163 idx = tree.find( idx )->second;
165 } while( idx != tree.find( idx )->second );
169 vtkSmartPointer< vtkPolyData > vtk_tree =
170 vtkSmartPointer< vtkPolyData >::New( );
171 vtk_tree->SetPoints( points );
172 vtk_tree->SetLines( cells );
173 vtk_tree->GetPointData( )->SetScalars( scalars );
175 view.AddPolyData( vtk_tree );
180 TImageWriter::Pointer dijkstra_writer =
181 TImageWriter::New( );
182 dijkstra_writer->SetInput( paths->GetOutput( ) );
183 dijkstra_writer->SetFileName( "dijkstra.mhd" );
184 dijkstra_writer->Update( );