8 #include <itkImageToVTKImageFilter.h>
10 #include <itkImageFileReader.h>
11 #include <fpa/VTK/ImageMPR.h>
13 #include <fpa/Image/DijkstraWithEndPointDetection.h>
14 #include <fpa/IO/MinimumSpanningTreeReader.h>
15 #include <fpa/IO/UniqueValuesContainerReader.h>
16 #include <fpa/IO/MatrixValuesContainerReader.h>
18 #include <vtkCellArray.h>
19 #include <vtkFloatArray.h>
20 #include <vtkImageMarchingCubes.h>
21 #include <vtkPoints.h>
22 #include <vtkPointData.h>
23 #include <vtkPolyData.h>
24 #include <vtkSmartPointer.h>
26 // -------------------------------------------------------------------------
27 const unsigned int Dim = 3;
28 typedef unsigned char TPixel;
30 typedef itk::Image< TPixel, Dim > TImage;
31 typedef itk::ImageToVTKImageFilter< TImage > TVTKInputImage;
33 typedef fpa::Image::DijkstraWithEndPointDetection< TImage, TImage > TFilter;
34 typedef TFilter::TMinimumSpanningTree TMinimumSpanningTree;
35 typedef TFilter::TUniqueVertices TUniqueVertices;
36 typedef TFilter::TBranches TBranches;
37 typedef TFilter::TVertices TVertices;
39 // -------------------------------------------------------------------------
40 int main( int argc, char* argv[] )
45 << "Usage: " << argv[ 0 ] << std::endl
46 << " input_image" << std::endl
47 << " input_minimum_spanning_tree" << std::endl
48 << " input_endpoints" << std::endl
49 << " input_bifurcations" << std::endl
50 << " input_branches" << std::endl
55 std::string input_image_fn = argv[ 1 ];
56 std::string mst_input_fn = argv[ 2 ];
57 std::string endpoints_input_fn = argv[ 3 ];
58 std::string bifurcations_input_fn = argv[ 4 ];
59 std::string branches_input_fn = argv[ 5 ];
62 itk::ImageFileReader< TImage >::Pointer input_image_reader =
63 itk::ImageFileReader< TImage >::New( );
64 input_image_reader->SetFileName( input_image_fn );
65 input_image_reader->Update( );
66 TImage::ConstPointer input_image = input_image_reader->GetOutput( );
68 TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( );
69 vtk_input_image->SetInput( input_image );
70 vtk_input_image->Update( );
72 // Read minimum spanning tree
73 fpa::IO::MinimumSpanningTreeReader< TMinimumSpanningTree >::Pointer
75 fpa::IO::MinimumSpanningTreeReader< TMinimumSpanningTree >::New( );
76 mst_input_reader->SetFileName( mst_input_fn );
77 mst_input_reader->Update( );
78 TMinimumSpanningTree::ConstPointer mst_input = mst_input_reader->GetOutput( );
80 fpa::IO::UniqueValuesContainerReader< TUniqueVertices >::Pointer
81 endpoints_input_reader =
82 fpa::IO::UniqueValuesContainerReader< TUniqueVertices >::New( );
83 endpoints_input_reader->SetFileName( endpoints_input_fn );
84 endpoints_input_reader->Update( );
85 TUniqueVertices::ConstPointer endpoints_input = endpoints_input_reader->GetOutput( );
87 fpa::IO::UniqueValuesContainerReader< TUniqueVertices >::Pointer
88 bifurcations_input_reader =
89 fpa::IO::UniqueValuesContainerReader< TUniqueVertices >::New( );
90 bifurcations_input_reader->SetFileName( bifurcations_input_fn );
91 bifurcations_input_reader->Update( );
92 TUniqueVertices::ConstPointer bifurcations_input =
93 bifurcations_input_reader->GetOutput( );
95 fpa::IO::MatrixValuesContainerReader< TBranches >::Pointer
96 branches_input_reader =
97 fpa::IO::MatrixValuesContainerReader< TBranches >::New( );
98 branches_input_reader->SetFileName( branches_input_fn );
99 branches_input_reader->Update( );
100 TBranches::ConstPointer branches_input = branches_input_reader->GetOutput( );
101 unsigned int nBranches = branches_input_reader->GetNumberOfLabels( );
103 // Show input image and let some interaction
104 fpa::VTK::ImageMPR view;
105 view.SetBackground( 0.3, 0.2, 0.8 );
106 view.SetSize( 800, 800 );
107 view.SetImage( vtk_input_image->GetOutput( ) );
109 vtkSmartPointer< vtkImageMarchingCubes > mc =
110 vtkSmartPointer< vtkImageMarchingCubes >::New( );
111 mc->SetInputData( vtk_input_image->GetOutput( ) );
112 mc->SetValue( 0, 1e-1 );
114 view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
117 vtkSmartPointer< vtkPoints > endpoints_points =
118 vtkSmartPointer< vtkPoints >::New( );
119 vtkSmartPointer< vtkCellArray > endpoints_cells =
120 vtkSmartPointer< vtkCellArray >::New( );
122 TUniqueVertices::ConstIterator epIt = endpoints_input->Begin( );
123 epIt != endpoints_input->End( );
127 TImage::PointType pnt;
128 input_image->TransformIndexToPhysicalPoint( *epIt, pnt );
129 endpoints_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
130 endpoints_cells->InsertNextCell( 1 );
132 InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 1 );
135 vtkSmartPointer< vtkPolyData > endpoints_polydata =
136 vtkSmartPointer< vtkPolyData >::New( );
137 endpoints_polydata->SetPoints( endpoints_points );
138 endpoints_polydata->SetVerts( endpoints_cells );
139 view.AddPolyData( endpoints_polydata, 0, 0, 1, 1 );
142 vtkSmartPointer< vtkPoints > bifurcations_points =
143 vtkSmartPointer< vtkPoints >::New( );
144 vtkSmartPointer< vtkCellArray > bifurcations_cells =
145 vtkSmartPointer< vtkCellArray >::New( );
147 TUniqueVertices::ConstIterator bfIt = bifurcations_input->Begin( );
148 bfIt != bifurcations_input->End( );
152 TImage::PointType pnt;
153 input_image->TransformIndexToPhysicalPoint( *bfIt, pnt );
154 bifurcations_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
155 bifurcations_cells->InsertNextCell( 1 );
157 InsertCellPoint( bifurcations_points->GetNumberOfPoints( ) - 1 );
160 vtkSmartPointer< vtkPolyData > bifurcations_polydata =
161 vtkSmartPointer< vtkPolyData >::New( );
162 bifurcations_polydata->SetPoints( bifurcations_points );
163 bifurcations_polydata->SetVerts( bifurcations_cells );
164 view.AddPolyData( bifurcations_polydata, 0, 1, 0, 1 );
166 // Show branches (simple and detailed)
167 vtkSmartPointer< vtkPoints > simple_branches_points =
168 vtkSmartPointer< vtkPoints >::New( );
169 vtkSmartPointer< vtkCellArray > simple_branches_cells =
170 vtkSmartPointer< vtkCellArray >::New( );
172 vtkSmartPointer< vtkPoints > detailed_branches_points =
173 vtkSmartPointer< vtkPoints >::New( );
174 vtkSmartPointer< vtkCellArray > detailed_branches_cells =
175 vtkSmartPointer< vtkCellArray >::New( );
176 vtkSmartPointer< vtkFloatArray > detailed_branches_scalars =
177 vtkSmartPointer< vtkFloatArray >::New( );
179 TBranches::ConstIterator brIt = branches_input->Begin( );
180 for( ; brIt != branches_input->End( ); ++brIt )
182 // Branch's first point
183 TImage::PointType first_point;
184 input_image->TransformIndexToPhysicalPoint( brIt->first, first_point );
185 unsigned long first_id = simple_branches_points->GetNumberOfPoints( );
186 simple_branches_points->InsertNextPoint(
187 first_point[ 0 ], first_point[ 1 ], first_point[ 2 ]
190 TBranches::ConstRowIterator brRowIt = branches_input->Begin( brIt );
191 for( ; brRowIt != branches_input->End( brIt ); ++brRowIt )
193 // Branch's second point
194 TImage::PointType second_point;
196 TransformIndexToPhysicalPoint( brRowIt->first, second_point );
197 unsigned long second_id = simple_branches_points->GetNumberOfPoints( );
198 simple_branches_points->InsertNextPoint(
199 second_point[ 0 ], second_point[ 1 ], second_point[ 2 ]
201 simple_branches_cells->InsertNextCell( 2 );
202 simple_branches_cells->InsertCellPoint( first_id );
203 simple_branches_cells->InsertCellPoint( second_id );
206 double pathId = double( brRowIt->second - 1 ) / double( nBranches - 1 );
208 mst_input->GetPath( path, brIt->first, brRowIt->first );
209 TVertices::const_iterator pIt = path.begin( );
210 for( ; pIt != path.end( ); ++pIt )
212 TImage::PointType path_point;
213 input_image->TransformIndexToPhysicalPoint( *pIt, path_point );
214 detailed_branches_points->InsertNextPoint(
215 path_point[ 0 ], path_point[ 1 ], path_point[ 2 ]
217 detailed_branches_scalars->InsertNextTuple1( pathId );
218 if( pIt != path.begin( ) )
220 unsigned long nPoints =
221 detailed_branches_points->GetNumberOfPoints( );
222 detailed_branches_cells->InsertNextCell( 2 );
223 detailed_branches_cells->InsertCellPoint( nPoints - 2 );
224 detailed_branches_cells->InsertCellPoint( nPoints - 1 );
233 vtkSmartPointer< vtkPolyData > simple_branches_polydata =
234 vtkSmartPointer< vtkPolyData >::New( );
235 simple_branches_polydata->SetPoints( simple_branches_points );
236 simple_branches_polydata->SetLines( simple_branches_cells );
237 view.AddPolyData( simple_branches_polydata, 1, 0, 1, 1 );
239 vtkSmartPointer< vtkPolyData > detailed_branches_polydata =
240 vtkSmartPointer< vtkPolyData >::New( );
241 detailed_branches_polydata->SetPoints( detailed_branches_points );
242 detailed_branches_polydata->SetLines( detailed_branches_cells );
243 detailed_branches_polydata->
244 GetPointData( )->SetScalars( detailed_branches_scalars );
245 view.AddPolyData( detailed_branches_polydata, 1 );
247 // Allow some interaction