8 #include <itkImageToVTKImageFilter.h>
10 #include <itkImageFileReader.h>
11 #include <itkMinimumMaximumImageCalculator.h>
12 #include <itkShiftScaleImageFilter.h>
14 #include <itkSignedMaurerDistanceMapImageFilter.h>
16 #include <itkImageFileWriter.h>
18 #include <fpa/Image/DijkstraWithEndPointDetection.h>
19 #include <fpa/Base/Functors/InvertCostFunction.h>
20 #include <fpa/VTK/ImageMPR.h>
21 #include <fpa/VTK/Image3DObserver.h>
23 #include <vtkCellArray.h>
24 #include <vtkImageMarchingCubes.h>
25 #include <vtkPoints.h>
26 #include <vtkPolyData.h>
27 #include <vtkSmartPointer.h>
29 // -------------------------------------------------------------------------
30 const unsigned int Dim = 3;
31 typedef unsigned char TPixel;
32 typedef float TScalar;
34 typedef itk::Image< TPixel, Dim > TImage;
35 typedef itk::Image< TScalar, Dim > TScalarImage;
36 typedef itk::ImageToVTKImageFilter< TImage > TVTKInputImage;
38 // -------------------------------------------------------------------------
40 void ReadImage( typename I::Pointer& image, const std::string& filename );
43 void SaveImage( const I* image, const std::string& filename );
45 template< class I, class O >
47 const typename I::Pointer& input, typename O::Pointer& output
50 // -------------------------------------------------------------------------
51 int main( int argc, char* argv[] )
56 << "Usage: " << argv[ 0 ]
57 << " input_image distancemap output_costmap output_labels"
62 std::string input_image_fn = argv[ 1 ];
63 std::string distancemap_fn = argv[ 2 ];
64 std::string output_costmap_fn = argv[ 3 ];
65 std::string output_labels_fn = argv[ 4 ];
68 TImage::Pointer input_image;
71 ReadImage< TImage >( input_image, input_image_fn );
73 catch( itk::ExceptionObject& err )
76 << "Error caught while reading \""
77 << input_image_fn << "\": " << err
83 TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( );
84 vtk_input_image->SetInput( input_image );
85 vtk_input_image->Update( );
87 // Show input image and let some interaction
88 fpa::VTK::ImageMPR view;
89 view.SetBackground( 0.3, 0.2, 0.8 );
90 view.SetSize( 800, 800 );
91 view.SetImage( vtk_input_image->GetOutput( ) );
93 vtkSmartPointer< vtkImageMarchingCubes > mc =
94 vtkSmartPointer< vtkImageMarchingCubes >::New( );
95 mc->SetInputData( vtk_input_image->GetOutput( ) );
96 mc->SetValue( 0, 1e-1 );
98 view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
100 // Allow some interaction and wait for, at least, one seed
102 while( view.GetNumberOfSeeds( ) == 0 )
107 view.GetSeed( 0, p );
108 TImage::PointType pnt;
109 TImage::IndexType seed;
110 pnt[ 0 ] = TImage::PointType::ValueType( p[ 0 ] );
111 pnt[ 1 ] = TImage::PointType::ValueType( p[ 1 ] );
112 pnt[ 2 ] = TImage::PointType::ValueType( p[ 2 ] );
113 input_image->TransformPhysicalPointToIndex( pnt, seed );
115 // Compute squared distance map
116 TScalarImage::Pointer dmap;
117 DistanceMap< TImage, TScalarImage >( input_image, dmap );
119 // Prepare cost conversion function
120 typedef fpa::Base::Functors::InvertCostFunction< TScalar > TFunction;
121 TFunction::Pointer function = TFunction::New( );
123 // Prepare Dijkstra filter
124 typedef fpa::Image::DijkstraWithEndPointDetection< TScalarImage, TScalarImage > TFilter;
125 TFilter::Pointer filter = TFilter::New( );
126 filter->SetInput( dmap );
127 filter->SetConversionFunction( function );
128 filter->SetNeighborhoodOrder( 2 );
129 filter->StopAtOneFrontOff( );
130 filter->AddSeed( seed, TScalar( 0 ) );
132 // Prepare graphical debugger
133 typedef fpa::VTK::Image3DObserver< TFilter, vtkRenderWindow > TDebugger;
134 TDebugger::Pointer debugger = TDebugger::New( );
135 debugger->SetRenderWindow( view.GetWindow( ) );
136 debugger->SetRenderPercentage( 0.0001 );
137 filter->AddObserver( itk::AnyEvent( ), debugger );
138 filter->ThrowEventsOn( );
141 std::time_t start, end;
146 << "Extraction time = "
147 << std::difftime( end, start )
148 << " s." << std::endl;
150 // Minimum spanning tree
151 const TFilter::TMinimumSpanningTree* mst =
152 filter->GetMinimumSpanningTree( );
154 // Build branches from endpoints to seed
155 vtkSmartPointer< vtkPoints > endpoints_points =
156 vtkSmartPointer< vtkPoints >::New( );
157 vtkSmartPointer< vtkCellArray > endpoints_cells =
158 vtkSmartPointer< vtkCellArray >::New( );
159 vtkSmartPointer< vtkCellArray > endpoints_vertices =
160 vtkSmartPointer< vtkCellArray >::New( );
162 const TFilter::TUniqueVertices& endpoints = filter->GetEndPoints( );
163 TFilter::TUniqueVertices::const_iterator eIt = endpoints.begin( );
164 for( ; eIt != endpoints.end( ); ++eIt )
166 TFilter::TVertices path;
167 mst->GetPath( path, *eIt, filter->GetSeed( 0 ) );
169 std::cout << *eIt << std::endl;
170 TFilter::TVertices::const_iterator pIt = path.begin( );
171 for( ; pIt != path.end( ); ++pIt )
173 TImage::PointType pnt;
174 input_image->TransformIndexToPhysicalPoint( *pIt, pnt );
175 endpoints_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
177 if( pIt != path.begin( ) )
179 endpoints_cells->InsertNextCell( 2 );
180 endpoints_cells->InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 2 );
181 endpoints_cells->InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 1 );
185 std::cout << "ok" << std::endl;
186 endpoints_vertices->InsertNextCell( 1 );
187 endpoints_vertices->InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 1 );
197 std::cout << endpoints.size( ) << std::endl;
199 vtkSmartPointer< vtkPolyData > endpoints_polydata =
200 vtkSmartPointer< vtkPolyData >::New( );
201 endpoints_polydata->SetPoints( endpoints_points );
202 endpoints_polydata->SetLines( endpoints_cells );
203 endpoints_polydata->SetVerts( endpoints_vertices );
216 vtkSmartPointer< vtkPoints > bifurcations_points =
217 vtkSmartPointer< vtkPoints >::New( );
218 vtkSmartPointer< vtkCellArray > bifurcations_vertices =
219 vtkSmartPointer< vtkCellArray >::New( );
221 const TFilter::TUniqueVertices& bifurcations = filter->GetBifurcationPoints( );
222 TFilter::TUniqueVertices::const_iterator bIt = bifurcations.begin( );
223 for( ; bIt != bifurcations.end( ); ++bIt )
225 std::cout << *bIt << std::endl;
227 TImage::PointType pnt;
228 input_image->TransformIndexToPhysicalPoint( *bIt, pnt );
229 bifurcations_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
230 bifurcations_vertices->InsertNextCell( 1 );
231 bifurcations_vertices->InsertCellPoint( bifurcations_points->GetNumberOfPoints( ) - 1 );
235 vtkSmartPointer< vtkPolyData > bifurcations_polydata =
236 vtkSmartPointer< vtkPolyData >::New( );
237 bifurcations_polydata->SetPoints( bifurcations_points );
238 bifurcations_polydata->SetVerts( bifurcations_vertices );
244 // Build branches from endpoints to seed
245 vtkSmartPointer< vtkPoints > branches_points =
246 vtkSmartPointer< vtkPoints >::New( );
247 vtkSmartPointer< vtkCellArray > branches_cells =
248 vtkSmartPointer< vtkCellArray >::New( );
250 const TFilter::TBranches& branches = filter->GetBranches( );
251 TFilter::TBranches::const_iterator brIt = branches.begin( );
252 for( ; brIt != branches.end( ); ++brIt )
254 TFilter::TBranch::const_iterator brsIt = brIt->second.begin( );
255 for( ; brsIt != brIt->second.end( ); ++brsIt )
257 TImage::PointType pnt0, pnt1;
258 input_image->TransformIndexToPhysicalPoint( brIt->first, pnt0 );
259 input_image->TransformIndexToPhysicalPoint( brsIt->first, pnt1 );
260 branches_points->InsertNextPoint( pnt0[ 0 ], pnt0[ 1 ], pnt0[ 2 ] );
261 branches_points->InsertNextPoint( pnt1[ 0 ], pnt1[ 1 ], pnt1[ 2 ] );
263 branches_cells->InsertNextCell( 2 );
264 branches_cells->InsertCellPoint( branches_points->GetNumberOfPoints( ) - 2 );
265 branches_cells->InsertCellPoint( branches_points->GetNumberOfPoints( ) - 1 );
271 vtkSmartPointer< vtkPolyData > branches_polydata =
272 vtkSmartPointer< vtkPolyData >::New( );
273 branches_polydata->SetPoints( branches_points );
274 branches_polydata->SetLines( branches_cells );
276 view.AddPolyData( endpoints_polydata, 1, 0, 0, 1 );
277 view.AddPolyData( bifurcations_polydata, 0, 0, 1, 1 );
278 view.AddPolyData( branches_polydata, 0, 1, 0, 1 );
282 // Save final total cost map
283 SaveImage( filter->GetOutput( ), output_costmap_fn );
284 SaveImage( dmap.GetPointer( ), distancemap_fn );
285 SaveImage( filter->GetLabelImage( ), output_labels_fn );
289 // -------------------------------------------------------------------------
291 void ReadImage( typename I::Pointer& image, const std::string& filename )
293 typename itk::ImageFileReader< I >::Pointer reader =
294 itk::ImageFileReader< I >::New( );
295 reader->SetFileName( filename );
298 typename itk::MinimumMaximumImageCalculator< I >::Pointer minmax =
299 itk::MinimumMaximumImageCalculator< I >::New( );
300 minmax->SetImage( reader->GetOutput( ) );
302 double vmin = double( minmax->GetMinimum( ) );
303 double vmax = double( minmax->GetMaximum( ) );
305 typename itk::ShiftScaleImageFilter< I, I >::Pointer shift =
306 itk::ShiftScaleImageFilter< I, I >::New( );
307 shift->SetInput( reader->GetOutput( ) );
308 shift->SetScale( vmax - vmin );
309 shift->SetShift( vmin / ( vmax - vmin ) );
312 image = shift->GetOutput( );
313 image->DisconnectPipeline( );
316 // -------------------------------------------------------------------------
318 void SaveImage( const I* image, const std::string& filename )
320 typename itk::ImageFileWriter< I >::Pointer writer =
321 itk::ImageFileWriter< I >::New( );
322 writer->SetInput( image );
323 writer->SetFileName( filename );
328 catch( itk::ExceptionObject& err )
331 << "Error saving \"" << filename << "\": " << err
337 // -------------------------------------------------------------------------
338 template< class I, class O >
340 const typename I::Pointer& input, typename O::Pointer& output
343 typename itk::SignedMaurerDistanceMapImageFilter< I, O >::Pointer dmap =
344 itk::SignedMaurerDistanceMapImageFilter< I, O >::New( );
345 dmap->SetInput( input );
346 dmap->SetBackgroundValue( ( typename I::PixelType )( 0 ) );
347 dmap->InsideIsPositiveOn( );
348 dmap->SquaredDistanceOn( );
349 dmap->UseImageSpacingOn( );
351 std::time_t start, end;
356 << "Distance map time = "
357 << std::difftime( end, start )
358 << " s." << std::endl;
360 output = dmap->GetOutput( );
361 output->DisconnectPipeline( );