8 #include <itkImageToVTKImageFilter.h>
10 #include <itkSignedMaurerDistanceMapImageFilter.h>
12 #include <fpa/Base/Functors/InvertCostFunction.h>
13 #include <fpa/Image/DijkstraWithEndPointDetection.h>
14 #include <fpa/VTK/Image2DObserver.h>
15 #include <fpa/IO/MinimumSpanningTreeWriter.h>
16 #include <fpa/IO/UniqueValuesContainerWriter.h>
17 #include <fpa/IO/MatrixValuesContainerWriter.h>
18 #include <fpa/VTK/UniqueVerticesToPolyDataFilter.h>
20 #include "fpa_Utility.h"
22 // -------------------------------------------------------------------------
23 const unsigned int Dim = 2;
24 typedef unsigned char TPixel;
25 typedef float TScalar;
26 typedef itk::Image< TPixel, Dim > TImage;
27 typedef itk::Image< TScalar, Dim > TScalarImage;
29 // -------------------------------------------------------------------------
30 int main( int argc, char* argv[] )
35 << "Usage: " << argv[ 0 ] << std::endl
36 << "\tinput_image" << std::endl
37 << "\toutput_labels" << std::endl
38 << "\toutput_minimum_spanning_tree" << std::endl
39 << "\toutput_endpoints" << std::endl
40 << "\toutput_bifurcations" << std::endl
41 << "\toutput_branches" << std::endl
46 std::string input_image_fn = argv[ 1 ];
47 std::string output_labels_fn = argv[ 2 ];
48 std::string output_minimum_spanning_tree_fn = argv[ 3 ];
49 std::string output_endpoints_fn = argv[ 4 ];
50 std::string output_bifurcations_fn = argv[ 5 ];
51 std::string output_branches_fn = argv[ 6 ];
54 TImage::Pointer input_image;
55 std::string err = fpa_Utility::ReadImage( input_image, input_image_fn );
58 std::cerr << err << std::endl;
63 // Show image and wait for, at least, one seed
64 itk::ImageToVTKImageFilter< TImage >::Pointer vtk_input_image =
65 itk::ImageToVTKImageFilter< TImage >::New( );
66 vtk_input_image->SetInput( input_image );
67 vtk_input_image->Update( );
69 fpa_Utility::Viewer2DWithSeeds viewer;
70 viewer.SetImage( vtk_input_image->GetOutput( ) );
71 while( viewer.GetNumberOfSeeds( ) == 0 )
74 // Compute squared distance map
76 itk::SignedMaurerDistanceMapImageFilter< TImage, TScalarImage >::Pointer
78 itk::SignedMaurerDistanceMapImageFilter< TImage, TScalarImage >::New( );
79 dmap->SetInput( input_image );
80 dmap->SetBackgroundValue( TPixel( 0 ) );
81 dmap->InsideIsPositiveOn( );
82 dmap->SquaredDistanceOn( );
83 dmap->UseImageSpacingOn( );
85 std::time_t start, end;
90 << "Distance map time = "
91 << std::difftime( end, start )
94 // Prepare cost conversion function
95 typedef fpa::Base::Functors::InvertCostFunction< TScalar > TCostFunction;
96 TCostFunction::Pointer cost_f = TCostFunction::New( );
98 // Prepare Dijkstra filter
100 DijkstraWithEndPointDetection< TScalarImage, TScalarImage > TFilter;
101 TFilter::Pointer filter = TFilter::New( );
102 filter->SetInput( dmap->GetOutput( ) );
103 filter->SetConversionFunction( cost_f );
104 filter->SetNeighborhoodOrder( 2 );
105 filter->StopAtOneFrontOff( );
108 TImage::PointType pnt;
109 TImage::IndexType idx;
110 viewer.GetSeed( pnt, 0 );
111 if( input_image->TransformPhysicalPointToIndex( pnt, idx ) )
112 filter->AddSeed( idx, 0 );
114 // Prepare graphical debugger
115 typedef fpa::VTK::Image2DObserver< TFilter, vtkRenderWindow > TDebugger;
116 TDebugger::Pointer debugger = TDebugger::New( );
117 debugger->SetRenderWindow( viewer.Window );
118 debugger->SetRenderPercentage( 0.01 );
119 filter->AddObserver( itk::AnyEvent( ), debugger );
120 filter->ThrowEventsOn( );
127 << "Extraction time = "
128 << std::difftime( end, start )
129 << "s." << std::endl;
132 typedef vtkSmartPointer< fpa::VTK::UniqueVerticesToPolyDataFilter< TFilter::TUniqueVertices, TFilter::TInputImage > > TVertices2PD;
134 TVertices2PD endpoints2pd = TVertices2PD::New( );
135 endpoints2pd->SetInput( filter->GetEndPoints( ) );
136 endpoints2pd->SetImage( filter->GetInput( ) );
137 endpoints2pd->Update( );
139 vtkSmartPointer< vtkPolyDataMapper > endpoints_mapper =
140 vtkSmartPointer< vtkPolyDataMapper >::New( );
141 endpoints_mapper->SetInputConnection( endpoints2pd->GetOutputPort( ) );
143 vtkSmartPointer< vtkActor > endpoints_actor =
144 vtkSmartPointer< vtkActor >::New( );
145 endpoints_actor->SetMapper( endpoints_mapper );
146 endpoints_actor->GetProperty( )->SetColor( 0, 1, 0 );
147 endpoints_actor->GetProperty( )->SetPointSize( 5 );
148 viewer.Renderer->AddActor( endpoints_actor );
150 TVertices2PD bifurcations2pd = TVertices2PD::New( );
151 bifurcations2pd->SetInput( filter->GetBifurcations( ) );
152 bifurcations2pd->SetImage( filter->GetInput( ) );
153 bifurcations2pd->Update( );
155 vtkSmartPointer< vtkPolyDataMapper > bifurcations_mapper =
156 vtkSmartPointer< vtkPolyDataMapper >::New( );
157 bifurcations_mapper->SetInputConnection( bifurcations2pd->GetOutputPort( ) );
159 vtkSmartPointer< vtkActor > bifurcations_actor =
160 vtkSmartPointer< vtkActor >::New( );
161 bifurcations_actor->SetMapper( bifurcations_mapper );
162 bifurcations_actor->GetProperty( )->SetColor( 1, 0, 0 );
163 bifurcations_actor->GetProperty( )->SetPointSize( 5 );
164 viewer.Renderer->AddActor( bifurcations_actor );
166 // Some more interaction and finish
171 err = fpa_Utility::SaveImage( filter->GetLabelImage( ), output_labels_fn );
173 std::cerr << err << std::endl;
175 fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::Pointer
177 fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::New( );
178 mst_writer->SetInput( filter->GetMinimumSpanningTree( ) );
179 mst_writer->SetFileName( output_minimum_spanning_tree_fn );
180 mst_writer->Update( );
182 fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer
184 fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( );
185 endpoints_writer->SetInput( filter->GetEndPoints( ) );
186 endpoints_writer->SetFileName( output_endpoints_fn );
187 endpoints_writer->Update( );
189 fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer
190 bifurcations_writer =
191 fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( );
192 bifurcations_writer->SetInput( filter->GetBifurcations( ) );
193 bifurcations_writer->SetFileName( output_bifurcations_fn );
194 bifurcations_writer->Update( );
196 fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::Pointer
198 fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::New( );
199 branches_writer->SetInput( filter->GetBranches( ) );
200 branches_writer->SetFileName( output_branches_fn );
201 branches_writer->Update( );
209 vtkSmartPointer< vtkPoints > simple_branches_points =
210 vtkSmartPointer< vtkPoints >::New( );
211 vtkSmartPointer< vtkCellArray > simple_branches_cells =
212 vtkSmartPointer< vtkCellArray >::New( );
214 vtkSmartPointer< vtkPoints > detailed_branches_points =
215 vtkSmartPointer< vtkPoints >::New( );
216 vtkSmartPointer< vtkCellArray > detailed_branches_cells =
217 vtkSmartPointer< vtkCellArray >::New( );
218 vtkSmartPointer< vtkFloatArray > detailed_branches_scalars =
219 vtkSmartPointer< vtkFloatArray >::New( );
221 TFilter::TBranches::ConstIterator brIt = branches->Begin( );
222 for( ; brIt != branches->End( ); ++brIt )
224 // Branch's first point
225 TImage::PointType first_point;
226 input_image->TransformIndexToPhysicalPoint( brIt->first, first_point );
227 unsigned long first_id = simple_branches_points->GetNumberOfPoints( );
228 simple_branches_points->InsertNextPoint(
229 first_point[ 0 ], first_point[ 1 ], first_point[ 2 ]
232 TFilter::TBranches::ConstRowIterator brRowIt = branches->Begin( brIt );
233 for( ; brRowIt != branches->End( brIt ); ++brRowIt )
235 // Branch's second point
236 TImage::PointType second_point;
238 TransformIndexToPhysicalPoint( brRowIt->first, second_point );
239 unsigned long second_id = simple_branches_points->GetNumberOfPoints( );
240 simple_branches_points->InsertNextPoint(
241 second_point[ 0 ], second_point[ 1 ], second_point[ 2 ]
243 simple_branches_cells->InsertNextCell( 2 );
244 simple_branches_cells->InsertCellPoint( first_id );
245 simple_branches_cells->InsertCellPoint( second_id );
248 double pathId = double( brRowIt->second - 1 ) / double( nBranches - 1 );
249 TFilter::TVertices path;
250 mst->GetPath( path, brIt->first, brRowIt->first );
251 TFilter::TVertices::const_iterator pIt = path.begin( );
252 for( ; pIt != path.end( ); ++pIt )
254 TImage::PointType path_point;
255 input_image->TransformIndexToPhysicalPoint( *pIt, path_point );
256 detailed_branches_points->InsertNextPoint(
257 path_point[ 0 ], path_point[ 1 ], path_point[ 2 ]
259 detailed_branches_scalars->InsertNextTuple1( pathId );
260 if( pIt != path.begin( ) )
262 unsigned long nPoints =
263 detailed_branches_points->GetNumberOfPoints( );
264 detailed_branches_cells->InsertNextCell( 2 );
265 detailed_branches_cells->InsertCellPoint( nPoints - 2 );
266 detailed_branches_cells->InsertCellPoint( nPoints - 1 );
275 vtkSmartPointer< vtkPolyData > simple_branches_polydata =
276 vtkSmartPointer< vtkPolyData >::New( );
277 simple_branches_polydata->SetPoints( simple_branches_points );
278 simple_branches_polydata->SetLines( simple_branches_cells );
279 view.AddPolyData( simple_branches_polydata, 1, 0, 1, 1 );
281 vtkSmartPointer< vtkPolyData > detailed_branches_polydata =
282 vtkSmartPointer< vtkPolyData >::New( );
283 detailed_branches_polydata->SetPoints( detailed_branches_points );
284 detailed_branches_polydata->SetLines( detailed_branches_cells );
285 detailed_branches_polydata->
286 GetPointData( )->SetScalars( detailed_branches_scalars );
287 view.AddPolyData( detailed_branches_polydata, 1 );
289 // Let some more interaction