8 #include <itkImageToVTKImageFilter.h>
10 #include <itkSignedMaurerDistanceMapImageFilter.h>
12 #include <vtkOutlineSource.h>
14 #include <fpa/Base/Functors/InvertCostFunction.h>
15 #include <fpa/Image/DijkstraWithEndPointDetection.h>
16 #include <fpa/VTK/Image3DObserver.h>
17 #include <fpa/IO/MinimumSpanningTreeWriter.h>
18 #include <fpa/IO/UniqueValuesContainerWriter.h>
19 #include <fpa/IO/MatrixValuesContainerWriter.h>
20 #include <fpa/VTK/UniqueVerticesToPolyDataFilter.h>
22 #include "fpa_Utility.h"
24 // -------------------------------------------------------------------------
25 const unsigned int Dim = 3;
26 typedef unsigned char TPixel;
27 typedef float TScalar;
28 typedef itk::Image< TPixel, Dim > TImage;
29 typedef itk::Image< TScalar, Dim > TScalarImage;
31 // -------------------------------------------------------------------------
32 int main( int argc, char* argv[] )
37 << "Usage: " << argv[ 0 ] << std::endl
38 << "\tinput_image" << std::endl
39 << "\toutput_labels" << std::endl
40 << "\toutput_minimum_spanning_tree" << std::endl
41 << "\toutput_endpoints" << std::endl
42 << "\toutput_bifurcations" << std::endl
43 << "\toutput_branches" << std::endl
48 std::string input_image_fn = argv[ 1 ];
49 std::string output_labels_fn = argv[ 2 ];
50 std::string output_minimum_spanning_tree_fn = argv[ 3 ];
51 std::string output_endpoints_fn = argv[ 4 ];
52 std::string output_bifurcations_fn = argv[ 5 ];
53 std::string output_branches_fn = argv[ 6 ];
56 TImage::Pointer input_image;
57 std::string err = fpa_Utility::ReadImage( input_image, input_image_fn );
60 std::cerr << err << std::endl;
65 // Show image and wait for, at least, one seed
66 itk::ImageToVTKImageFilter< TImage >::Pointer vtk_input_image =
67 itk::ImageToVTKImageFilter< TImage >::New( );
68 vtk_input_image->SetInput( input_image );
69 vtk_input_image->Update( );
71 vtkSmartPointer<vtkOutlineSource> outline =
72 vtkSmartPointer<vtkOutlineSource>::New();
73 outline->SetBounds(vtk_input_image->GetOutput( )->GetBounds());
75 vtkSmartPointer<vtkPolyDataMapper> mapper =
76 vtkSmartPointer<vtkPolyDataMapper>::New( );
77 mapper->SetInputConnection( outline->GetOutputPort( ) );
79 vtkSmartPointer<vtkActor> actor =
80 vtkSmartPointer<vtkActor>::New( );
81 actor->SetMapper( mapper );
82 vtkSmartPointer<vtkRenderer> aRenderer =
83 vtkSmartPointer<vtkRenderer>::New();
84 vtkSmartPointer<vtkRenderWindow> renWin =
85 vtkSmartPointer<vtkRenderWindow>::New();
86 renWin->AddRenderer(aRenderer);
88 vtkSmartPointer<vtkRenderWindowInteractor> iren =
89 vtkSmartPointer<vtkRenderWindowInteractor>::New();
90 iren->SetRenderWindow(renWin);
92 aRenderer->AddActor( actor );
97 fpa_Utility::Viewer2DWithSeeds viewer;
98 viewer.SetImage( vtk_input_image->GetOutput( ) );
99 while( viewer.GetNumberOfSeeds( ) == 0 )
103 // Compute squared distance map
104 itk::SignedMaurerDistanceMapImageFilter< TImage, TScalarImage >::Pointer
106 itk::SignedMaurerDistanceMapImageFilter< TImage, TScalarImage >::New( );
107 dmap->SetInput( input_image );
108 dmap->SetBackgroundValue( TPixel( 0 ) );
109 dmap->InsideIsPositiveOn( );
110 dmap->SquaredDistanceOn( );
111 dmap->UseImageSpacingOn( );
113 std::time_t start, end;
118 << "Distance map time = "
119 << std::difftime( end, start )
120 << "s." << std::endl;
122 // Prepare cost conversion function
123 typedef fpa::Base::Functors::InvertCostFunction< TScalar > TCostFunction;
124 TCostFunction::Pointer cost_f = TCostFunction::New( );
126 // Prepare Dijkstra filter
128 DijkstraWithEndPointDetection< TScalarImage, TScalarImage > TFilter;
129 TFilter::Pointer filter = TFilter::New( );
130 filter->SetInput( dmap->GetOutput( ) );
131 filter->SetConversionFunction( cost_f );
132 filter->SetNeighborhoodOrder( 2 );
133 filter->StopAtOneFrontOff( );
136 TImage::PointType pnt;
137 TImage::IndexType idx;
139 viewer.GetSeed( pnt, 0 );
142 pnt[ 1 ] = -109.36591;
143 pnt[ 2 ] = 1942.480988;
145 if( input_image->TransformPhysicalPointToIndex( pnt, idx ) )
146 filter->AddSeed( idx, 0 );
148 // Prepare graphical debugger
149 typedef fpa::VTK::Image3DObserver< TFilter, vtkRenderWindow > TDebugger;
150 TDebugger::Pointer debugger = TDebugger::New( );
151 debugger->SetRenderWindow( renWin );
152 debugger->SetRenderPercentage( 0.0001 );
153 filter->AddObserver( itk::AnyEvent( ), debugger );
154 filter->ThrowEventsOn( );
161 << "Extraction time = "
162 << std::difftime( end, start )
163 << "s." << std::endl;
166 typedef vtkSmartPointer< fpa::VTK::UniqueVerticesToPolyDataFilter< TFilter::TUniqueVertices, TFilter::TInputImage > > TVertices2PD;
168 TVertices2PD endpoints2pd = TVertices2PD::New( );
169 endpoints2pd->SetInput( filter->GetEndPoints( ) );
170 endpoints2pd->SetImage( filter->GetInput( ) );
171 endpoints2pd->Update( );
173 vtkSmartPointer< vtkPolyDataMapper > endpoints_mapper =
174 vtkSmartPointer< vtkPolyDataMapper >::New( );
175 endpoints_mapper->SetInputConnection( endpoints2pd->GetOutputPort( ) );
177 vtkSmartPointer< vtkActor > endpoints_actor =
178 vtkSmartPointer< vtkActor >::New( );
179 endpoints_actor->SetMapper( endpoints_mapper );
180 endpoints_actor->GetProperty( )->SetColor( 0, 1, 0 );
181 endpoints_actor->GetProperty( )->SetPointSize( 5 );
182 aRenderer->AddActor( endpoints_actor );
184 TVertices2PD bifurcations2pd = TVertices2PD::New( );
185 bifurcations2pd->SetInput( filter->GetBifurcations( ) );
186 bifurcations2pd->SetImage( filter->GetInput( ) );
187 bifurcations2pd->Update( );
189 vtkSmartPointer< vtkPolyDataMapper > bifurcations_mapper =
190 vtkSmartPointer< vtkPolyDataMapper >::New( );
191 bifurcations_mapper->SetInputConnection( bifurcations2pd->GetOutputPort( ) );
193 vtkSmartPointer< vtkActor > bifurcations_actor =
194 vtkSmartPointer< vtkActor >::New( );
195 bifurcations_actor->SetMapper( bifurcations_mapper );
196 bifurcations_actor->GetProperty( )->SetColor( 1, 0, 0 );
197 bifurcations_actor->GetProperty( )->SetPointSize( 5 );
198 aRenderer->AddActor( bifurcations_actor );
200 // Some more interaction and finish
205 err = fpa_Utility::SaveImage( filter->GetLabelImage( ), output_labels_fn );
207 std::cerr << err << std::endl;
209 fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::Pointer
211 fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::New( );
212 mst_writer->SetInput( filter->GetMinimumSpanningTree( ) );
213 mst_writer->SetFileName( output_minimum_spanning_tree_fn );
214 mst_writer->Update( );
216 fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer
218 fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( );
219 endpoints_writer->SetInput( filter->GetEndPoints( ) );
220 endpoints_writer->SetFileName( output_endpoints_fn );
221 endpoints_writer->Update( );
223 fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer
224 bifurcations_writer =
225 fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( );
226 bifurcations_writer->SetInput( filter->GetBifurcations( ) );
227 bifurcations_writer->SetFileName( output_bifurcations_fn );
228 bifurcations_writer->Update( );
230 fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::Pointer
232 fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::New( );
233 branches_writer->SetInput( filter->GetBranches( ) );
234 branches_writer->SetFileName( output_branches_fn );
235 branches_writer->Update( );
243 vtkSmartPointer< vtkPoints > simple_branches_points =
244 vtkSmartPointer< vtkPoints >::New( );
245 vtkSmartPointer< vtkCellArray > simple_branches_cells =
246 vtkSmartPointer< vtkCellArray >::New( );
248 vtkSmartPointer< vtkPoints > detailed_branches_points =
249 vtkSmartPointer< vtkPoints >::New( );
250 vtkSmartPointer< vtkCellArray > detailed_branches_cells =
251 vtkSmartPointer< vtkCellArray >::New( );
252 vtkSmartPointer< vtkFloatArray > detailed_branches_scalars =
253 vtkSmartPointer< vtkFloatArray >::New( );
255 TFilter::TBranches::ConstIterator brIt = branches->Begin( );
256 for( ; brIt != branches->End( ); ++brIt )
258 // Branch's first point
259 TImage::PointType first_point;
260 input_image->TransformIndexToPhysicalPoint( brIt->first, first_point );
261 unsigned long first_id = simple_branches_points->GetNumberOfPoints( );
262 simple_branches_points->InsertNextPoint(
263 first_point[ 0 ], first_point[ 1 ], first_point[ 2 ]
266 TFilter::TBranches::ConstRowIterator brRowIt = branches->Begin( brIt );
267 for( ; brRowIt != branches->End( brIt ); ++brRowIt )
269 // Branch's second point
270 TImage::PointType second_point;
272 TransformIndexToPhysicalPoint( brRowIt->first, second_point );
273 unsigned long second_id = simple_branches_points->GetNumberOfPoints( );
274 simple_branches_points->InsertNextPoint(
275 second_point[ 0 ], second_point[ 1 ], second_point[ 2 ]
277 simple_branches_cells->InsertNextCell( 2 );
278 simple_branches_cells->InsertCellPoint( first_id );
279 simple_branches_cells->InsertCellPoint( second_id );
282 double pathId = double( brRowIt->second - 1 ) / double( nBranches - 1 );
283 TFilter::TVertices path;
284 mst->GetPath( path, brIt->first, brRowIt->first );
285 TFilter::TVertices::const_iterator pIt = path.begin( );
286 for( ; pIt != path.end( ); ++pIt )
288 TImage::PointType path_point;
289 input_image->TransformIndexToPhysicalPoint( *pIt, path_point );
290 detailed_branches_points->InsertNextPoint(
291 path_point[ 0 ], path_point[ 1 ], path_point[ 2 ]
293 detailed_branches_scalars->InsertNextTuple1( pathId );
294 if( pIt != path.begin( ) )
296 unsigned long nPoints =
297 detailed_branches_points->GetNumberOfPoints( );
298 detailed_branches_cells->InsertNextCell( 2 );
299 detailed_branches_cells->InsertCellPoint( nPoints - 2 );
300 detailed_branches_cells->InsertCellPoint( nPoints - 1 );
309 vtkSmartPointer< vtkPolyData > simple_branches_polydata =
310 vtkSmartPointer< vtkPolyData >::New( );
311 simple_branches_polydata->SetPoints( simple_branches_points );
312 simple_branches_polydata->SetLines( simple_branches_cells );
313 view.AddPolyData( simple_branches_polydata, 1, 0, 1, 1 );
315 vtkSmartPointer< vtkPolyData > detailed_branches_polydata =
316 vtkSmartPointer< vtkPolyData >::New( );
317 detailed_branches_polydata->SetPoints( detailed_branches_points );
318 detailed_branches_polydata->SetLines( detailed_branches_cells );
319 detailed_branches_polydata->
320 GetPointData( )->SetScalars( detailed_branches_scalars );
321 view.AddPolyData( detailed_branches_polydata, 1 );
323 // Let some more interaction