5 #include <itkConstNeighborhoodIterator.h>
6 #include <itkDanielssonDistanceMapImageFilter.h>
8 #include <itkImageFileReader.h>
9 #include <itkImageFileWriter.h>
10 #include <itkImageToVTKImageFilter.h>
12 #include <vtkImageMarchingCubes.h>
14 #include <fpa/Image/RegionGrowWithMultipleThresholds.h>
15 #include <fpa/Image/DijkstraWithSphereBacktracking.h>
16 #include <fpa/VTK/ImageMPR.h>
17 #include <fpa/VTK/Image3DObserver.h>
19 // -------------------------------------------------------------------------
20 const unsigned int Dim = 3;
22 typedef float TScalar;
23 typedef itk::Image< TPixel, Dim > TImage;
24 typedef itk::Image< TScalar, Dim > TDistanceMap;
25 typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
26 typedef itk::ImageToVTKImageFilter< TDistanceMap > TVTKDistanceMap;
28 typedef itk::ImageFileReader< TImage > TImageReader;
29 typedef itk::ImageFileWriter< TImage > TImageWriter;
30 typedef itk::ImageFileWriter< TDistanceMap > TDistanceMapWriter;
31 typedef fpa::Image::RegionGrowWithMultipleThresholds< TImage > TSegmentor;
33 itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >
36 fpa::Image::DijkstraWithSphereBacktracking< TDistanceMap, TScalar >
39 typedef fpa::VTK::ImageMPR TMPR;
40 typedef fpa::VTK::Image3DObserver< TSegmentor, vtkRenderWindow > TSegmentorObs;
41 typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs;
43 // -------------------------------------------------------------------------
44 int main( int argc, char* argv[] )
49 << "Usage: " << argv[ 0 ]
50 << " input_image thr_0 thr_1 step"
51 << " output_segmentation output_distancemap output_dijkstra"
57 std::string input_image_fn = argv[ 1 ];
58 TPixel thr_0 = TPixel( std::atof( argv[ 2 ] ) );
59 TPixel thr_1 = TPixel( std::atof( argv[ 3 ] ) );
60 unsigned int step = std::atoi( argv[ 4 ] );
61 std::string output_segmentation_fn = argv[ 5 ];
62 std::string output_distancemap_fn = argv[ 6 ];
63 std::string output_dijkstra_fn = argv[ 7 ];
64 bool visual_debug = false;
66 visual_debug = ( std::atoi( argv[ 8 ] ) == 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 // Wait for a seed to be given
94 while( view.GetNumberOfSeeds( ) == 0 )
99 view.GetSeed( 0, seed );
100 TImage::PointType seed_pnt;
101 seed_pnt[ 0 ] = seed[ 0 ];
102 seed_pnt[ 1 ] = seed[ 1 ];
103 seed_pnt[ 2 ] = seed[ 2 ];
104 TImage::IndexType seed_idx;
105 input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx );
107 // Segment input image
108 TSegmentor::Pointer segmentor = TSegmentor::New( );
109 segmentor->AddThresholds( thr_0, thr_1, step );
110 segmentor->AddSeed( seed_idx, 0 );
111 segmentor->SetInput( input_image );
112 segmentor->SetNeighborhoodOrder( 1 );
113 segmentor->SetDifferenceThreshold( double( 3 ) );
114 segmentor->SetInsideValue( TPixel( 0 ) ); // WARNING: This is to optimize
115 segmentor->SetOutsideValue( TPixel( 1 ) ); // distance map calculation
118 // Configure observer
119 TSegmentorObs::Pointer obs = TSegmentorObs::New( );
120 obs->SetRenderWindow( view.GetWindow( ) );
121 segmentor->AddObserver( itk::AnyEvent( ), obs );
122 segmentor->ThrowEventsOn( );
125 segmentor->ThrowEventsOff( );
126 std::clock_t start = std::clock( );
127 segmentor->Update( );
128 std::clock_t end = std::clock( );
129 double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
130 std::cout << "Segmentation time = " << seconds << std::endl;
132 // Extract distance map
133 TDistanceFilter::Pointer distanceMap = TDistanceFilter::New( );
134 distanceMap->SetInput( segmentor->GetOutput( ) );
135 distanceMap->InputIsBinaryOn( );
136 distanceMap->SquaredDistanceOn( );
137 distanceMap->UseImageSpacingOn( );
138 start = std::clock( );
139 distanceMap->Update( );
141 seconds = double( end - start ) / double( CLOCKS_PER_SEC );
142 std::cout << "Distance map time = " << seconds << std::endl;
144 TVTKDistanceMap::Pointer vtk_distanceMap = TVTKDistanceMap::New( );
145 vtk_distanceMap->SetInput( distanceMap->GetOutput( ) );
146 vtk_distanceMap->Update( );
148 vtkSmartPointer< vtkImageMarchingCubes > mc =
149 vtkSmartPointer< vtkImageMarchingCubes >::New( );
150 mc->SetInputData( vtk_distanceMap->GetOutput( ) );
151 mc->SetValue( 0, 1e-2 );
153 view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
157 TDijkstra::Pointer paths = TDijkstra::New( );
158 paths->AddSeed( segmentor->GetSeed( 0 ), TScalar( 0 ) );
159 paths->SetInput( distanceMap->GetOutput( ) );
160 paths->SetNeighborhoodOrder( 2 );
164 // Configure observer
165 TDijkstraObs::Pointer obs = TDijkstraObs::New( );
166 obs->SetRenderWindow( view.GetWindow( ) );
167 paths->AddObserver( itk::AnyEvent( ), obs );
168 paths->ThrowEventsOn( );
171 paths->ThrowEventsOff( );
172 start = std::clock( );
175 seconds = double( end - start ) / double( CLOCKS_PER_SEC );
176 std::cout << "Paths extraction time = " << seconds << std::endl;
178 std::cout << "Final seed: " << paths->GetSeed( 0 ) << std::endl;
181 vtkSmartPointer< vtkPoints > points =
182 vtkSmartPointer< vtkPoints >::New( );
183 vtkSmartPointer< vtkCellArray > cells =
184 vtkSmartPointer< vtkCellArray >::New( );
185 vtkSmartPointer< vtkFloatArray > scalars =
186 vtkSmartPointer< vtkFloatArray >::New( );
188 const TDijkstra::TVertices& endpoints = paths->GetEndPoints( );
189 const TDijkstra::TTree& tree = paths->GetFinalTree( );
190 TDijkstra::TVertices::const_iterator epIt = endpoints.begin( );
191 for( unsigned int epId = 0; epIt != endpoints.end( ); ++epIt, ++epId )
193 double pd = double( epId ) / double( endpoints.size( ) - 1 );
195 TDijkstra::TVertex idx = *epIt;
198 TImage::PointType pnt;
199 input_image->TransformIndexToPhysicalPoint( idx, pnt );
201 points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
202 scalars->InsertNextTuple1( pd );
205 cells->InsertNextCell( 2 );
206 cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
207 cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
210 idx = tree.find( idx )->second;
212 } while( idx != tree.find( idx )->second );
216 vtkSmartPointer< vtkPolyData > vtk_tree =
217 vtkSmartPointer< vtkPolyData >::New( );
218 vtk_tree->SetPoints( points );
219 vtk_tree->SetLines( cells );
220 vtk_tree->GetPointData( )->SetScalars( scalars );
222 view.AddPolyData( vtk_tree );
226 itk::ImageFileWriter< TImage >::Pointer segmentation_writer =
227 itk::ImageFileWriter< TImage >::New( );
228 segmentation_writer->SetInput( segmentor->GetOutput( ) );
229 segmentation_writer->SetFileName( output_segmentation_fn );
231 itk::ImageFileWriter< TDistanceMap >::Pointer dmap_writer =
232 itk::ImageFileWriter< TDistanceMap >::New( );
233 dmap_writer->SetInput( distanceMap->GetOutput( ) );
234 dmap_writer->SetFileName( output_distancemap_fn );
236 itk::ImageFileWriter< TDistanceMap >::Pointer dijk_writer =
237 itk::ImageFileWriter< TDistanceMap >::New( );
238 dijk_writer->SetInput( paths->GetOutput( ) );
239 dijk_writer->SetFileName( output_dijkstra_fn );
243 segmentation_writer->Update( );
244 dmap_writer->Update( );
245 dijk_writer->Update( );
247 catch( itk::ExceptionObject& err )
249 std::cerr << "Error caught: " << err << std::endl;