11 #include <itkConstNeighborhoodIterator.h>
12 #include <itkNeighborhoodIterator.h>
13 #include <itkDanielssonDistanceMapImageFilter.h>
15 #include <itkImageFileReader.h>
16 #include <itkImageFileWriter.h>
17 #include <itkImageToVTKImageFilter.h>
19 #include <itkMinimumMaximumImageCalculator.h>
20 #include <itkInvertIntensityImageFilter.h>
22 #include <vtkPoints.h>
23 #include <vtkCellArray.h>
24 #include <vtkFloatArray.h>
25 #include <vtkPolyData.h>
26 #include <vtkPolyDataWriter.h>
27 #include <vtkSmartPointer.h>
28 #include <vtkImageMarchingCubes.h>
29 #include <vtkLookupTable.h>
31 #include <fpa/Image/DijkstraWithSphereBacktracking.h>
32 #include <fpa/VTK/ImageMPR.h>
33 #include <fpa/VTK/Image3DObserver.h>
35 // -------------------------------------------------------------------------
36 const unsigned int Dim = 3;
38 typedef float TScalar;
39 typedef itk::Image< TPixel, Dim > TImage;
40 typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
42 typedef itk::ImageFileReader< TImage > TImageReader;
43 typedef itk::ImageFileWriter< TImage > TImageWriter;
44 typedef fpa::Image::DijkstraWithSphereBacktracking< TImage, TScalar > TDijkstra;
46 typedef fpa::VTK::ImageMPR TMPR;
47 typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs;
52 TImage::PointType::VectorType V1;
53 TImage::PointType::VectorType V2;
55 bool operator<( const TBranch& other ) const
57 return( other.Length < this->Length );
60 typedef std::multiset< TBranch > TBranches;
62 // -------------------------------------------------------------------------
63 int main( int argc, char* argv[] )
68 << "Usage: " << argv[ 0 ]
69 << " input_image x y z output_image output_imagej output_vtk"
75 std::string input_image_fn = argv[ 1 ];
76 TImage::PointType seed_pnt;
77 seed_pnt[ 0 ] = std::atof( argv[ 2 ] );
78 seed_pnt[ 1 ] = std::atof( argv[ 3 ] );
79 seed_pnt[ 2 ] = std::atof( argv[ 4 ] );
80 std::string output_image_fn = argv[ 5 ];
81 std::string output_imagej_fn = argv[ 6 ];
82 std::string output_vtk_fn = argv[ 7 ];
83 bool visual_debug = false;
85 visual_debug = ( std::atoi( argv[ 8 ] ) == 1 );
88 TImageReader::Pointer input_image_reader = TImageReader::New( );
89 input_image_reader->SetFileName( input_image_fn );
92 input_image_reader->Update( );
94 catch( itk::ExceptionObject& err )
96 std::cerr << "Error caught: " << err << std::endl;
100 TImage::ConstPointer input_image = input_image_reader->GetOutput( );
103 TVTKImage::Pointer vtk_image = TVTKImage::New( );
104 vtk_image->SetInput( input_image );
105 vtk_image->Update( );
108 view.SetBackground( 0.3, 0.2, 0.8 );
109 view.SetSize( 800, 800 );
110 view.SetImage( vtk_image->GetOutput( ) );
113 vtkSmartPointer< vtkImageMarchingCubes > mc =
114 vtkSmartPointer< vtkImageMarchingCubes >::New( );
115 mc->SetInputData( vtk_image->GetOutput( ) );
116 mc->SetValue( 0, 1e-2 );
118 view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
121 // Allow some interaction
126 TImage::IndexType seed_idx;
127 input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx );
129 std::cout << seed_idx << " " << seed_pnt << std::endl;
136 TDijkstra::Pointer paths = TDijkstra::New( );
137 paths->AddSeed( seed_idx, TScalar( 0 ) );
138 paths->SetInput( input_image );
139 paths->SetNeighborhoodOrder( 2 );
143 // Configure observer
144 TDijkstraObs::Pointer obs = TDijkstraObs::New( );
145 obs->SetRenderWindow( view.GetWindow( ) );
146 paths->AddObserver( itk::AnyEvent( ), obs );
147 paths->ThrowEventsOn( );
150 paths->ThrowEventsOff( );
151 std::clock_t start = std::clock( );
153 std::clock_t end = std::clock( );
154 double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
155 std::cout << "Paths extraction time = " << seconds << std::endl;
158 vtkSmartPointer< vtkPoints > points =
159 vtkSmartPointer< vtkPoints >::New( );
160 vtkSmartPointer< vtkCellArray > cells =
161 vtkSmartPointer< vtkCellArray >::New( );
162 vtkSmartPointer< vtkFloatArray > scalars =
163 vtkSmartPointer< vtkFloatArray >::New( );
165 TDijkstra::TMarkImage::Pointer new_marks = TDijkstra::TMarkImage::New( );
166 new_marks->SetLargestPossibleRegion( input_image->GetLargestPossibleRegion( ) );
167 new_marks->SetRequestedRegion( input_image->GetRequestedRegion( ) );
168 new_marks->SetBufferedRegion( input_image->GetBufferedRegion( ) );
169 new_marks->SetOrigin( input_image->GetOrigin( ) );
170 new_marks->SetSpacing( input_image->GetSpacing( ) );
171 new_marks->SetDirection( input_image->GetDirection( ) );
172 new_marks->Allocate( );
173 new_marks->FillBuffer( 0 );
175 const TDijkstra::TMarkImage* marks = paths->GetOutputMarkImage( );
176 TDijkstra::TMark max_mark = paths->GetNumberOfBranches( );
177 const TDijkstra::TVertices& endpoints = paths->GetEndPoints( );
178 const TDijkstra::TVertices& bifurcations = paths->GetBifurcationPoints( );
179 const TDijkstra::TTree& tree = paths->GetFullTree( );
180 const TDijkstra::TTree& reduced_tree = paths->GetReducedTree( );
181 TDijkstra::TVertices::const_iterator epIt = endpoints.begin( );
184 TImage::PointType ori = input_image->GetOrigin( );
186 TDijkstra::TTree::const_iterator rIt = reduced_tree.begin( );
187 for( ; rIt != reduced_tree.end( ); ++rIt )
189 TDijkstra::TTree::const_iterator tIt = tree.find( rIt->first );
191 // Prepare branch "a la ImageJ"
193 TImage::PointType p1, p2;
194 input_image->TransformIndexToPhysicalPoint( rIt->second.first, p1 );
195 input_image->TransformIndexToPhysicalPoint( rIt->first, p2 );
196 branch.V1 = p1 - ori;
197 branch.V2 = p2 - ori;
198 branch.Length = double( 0 );
201 ( double( tIt->second.second ) - double( 1 ) ) /
202 ( double( max_mark ) - double( 1 ) );
204 TImage::PointType prev;
207 TDijkstra::TVertex idx = tIt->first;
208 new_marks->SetPixel( idx, 255 );
209 TImage::PointType pnt;
210 input_image->TransformIndexToPhysicalPoint( idx, pnt );
211 points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
212 scalars->InsertNextTuple1( pd );
215 cells->InsertNextCell( 2 );
216 cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
217 cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
218 branch.Length += prev.EuclideanDistanceTo( pnt );
223 tIt = tree.find( tIt->second.first );
225 } while( tIt->first != rIt->second.first );
228 TDijkstra::TVertex idx = tIt->first;
229 new_marks->SetPixel( idx, 255 );
230 TImage::PointType pnt;
231 input_image->TransformIndexToPhysicalPoint( idx, pnt );
232 points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
233 scalars->InsertNextTuple1( pd );
236 cells->InsertNextCell( 2 );
237 cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
238 cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
239 branch.Length += prev.EuclideanDistanceTo( pnt );
242 branches.insert( branch );
246 std::ofstream bout( output_imagej_fn.c_str( ) );
249 std::cerr << "Error: could not open " << output_imagej_fn << std::endl;
254 for( TBranches::const_iterator bIt = branches.begin( ); bIt != branches.end( ); ++bIt, ++i )
257 << std::setiosflags( std::ios::fixed) << std::setprecision( 3 )
259 << bIt->Length << "\t"
260 << bIt->V1[ 0 ] << "\t"
261 << bIt->V1[ 1 ] << "\t"
262 << bIt->V1[ 2 ] << "\t"
263 << bIt->V2[ 0 ] << "\t"
264 << bIt->V2[ 1 ] << "\t"
265 << bIt->V2[ 2 ] << "\t"
266 << ( bIt->V2 - bIt->V1 ).GetNorm( )
272 vtkSmartPointer< vtkPolyData > vtk_tree =
273 vtkSmartPointer< vtkPolyData >::New( );
274 vtk_tree->SetPoints( points );
275 vtk_tree->SetLines( cells );
276 vtk_tree->GetPointData( )->SetScalars( scalars );
278 view.AddPolyData( vtk_tree );
283 vtkSmartPointer< vtkPolyDataWriter > writer =
284 vtkSmartPointer< vtkPolyDataWriter >::New( );
285 writer->SetInputData( vtk_tree );
286 writer->SetFileName( output_vtk_fn.c_str( ) );
289 itk::ImageFileWriter< TDijkstra::TMarkImage >::Pointer marks_writer =
290 itk::ImageFileWriter< TDijkstra::TMarkImage >::New( );
291 marks_writer->SetInput( new_marks );
292 marks_writer->SetFileName( output_image_fn );
293 marks_writer->Update( );
296 TImageWriter::Pointer dijkstra_writer =
297 TImageWriter::New( );
298 dijkstra_writer->SetInput( paths->GetOutput( ) );
299 dijkstra_writer->SetFileName( "dijkstra.mhd" );
300 dijkstra_writer->Update( );