9 #include <itkConstNeighborhoodIterator.h>
10 #include <itkNeighborhoodIterator.h>
11 #include <itkDanielssonDistanceMapImageFilter.h>
13 #include <itkImageFileReader.h>
14 #include <itkImageFileWriter.h>
15 #include <itkImageToVTKImageFilter.h>
17 #include <vtkPoints.h>
18 #include <vtkCellArray.h>
19 #include <vtkPolyData.h>
20 #include <vtkSmartPointer.h>
22 #include <fpa/Base/Functors/InvertCostFunction.h>
23 #include <fpa/Image/Dijkstra.h>
24 #include <fpa/Image/RegionGrowWithMultipleThresholds.h>
25 #include <fpa/VTK/ImageMPR.h>
26 #include <fpa/VTK/Image3DObserver.h>
28 // -------------------------------------------------------------------------
29 const unsigned int Dim = 3;
30 typedef double TPixel;
31 typedef double TScalar;
32 typedef itk::Image< TPixel, Dim > TImage;
33 typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
35 typedef itk::ImageFileReader< TImage > TImageReader;
36 typedef itk::ImageFileWriter< TImage > TImageWriter;
37 typedef fpa::Image::Dijkstra< TImage, TScalar > TBaseDijkstra;
39 typedef fpa::VTK::ImageMPR TMPR;
40 typedef fpa::VTK::Image3DObserver< TBaseDijkstra, vtkRenderWindow > TDijkstraObs;
42 // -------------------------------------------------------------------------
44 : public TBaseDijkstra
47 typedef TDijkstra Self;
48 typedef TBaseDijkstra Superclass;
49 typedef itk::SmartPointer< Self > Pointer;
50 typedef itk::SmartPointer< const Self > ConstPointer;
52 typedef Superclass::TCost TCost;
53 typedef Superclass::TVertex TVertex;
54 typedef Superclass::InputImageType TImage;
55 typedef std::deque< TVertex > TVertices;
58 typedef std::pair< TCost, TVertex > _TCandidate;
59 typedef std::multimap< TCost, TVertex > _TCandidates;
60 typedef Superclass::_TNode _TNode;
64 itkTypeMacro( TDijkstra, TBaseDijkstra );
73 virtual void _BeforeMainLoop( )
75 const TImage* img = this->GetInput( );
78 TImage::SizeType radius;
80 itk::ConstNeighborhoodIterator< TImage > it(
81 radius, img, img->GetRequestedRegion( )
83 for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
85 _TNode seed = this->m_Seeds[ s ];
86 it.SetLocation( seed.Vertex );
88 TImage::SizeType size = it.GetSize( );
90 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
92 TImage::PixelType max_value = img->GetPixel( seed.Vertex );
93 for( unsigned int i = 0; i < nN; ++i )
95 if( it.GetPixel( i ) > max_value )
97 seed.Vertex = it.GetIndex( i );
98 seed.Parent = seed.Vertex;
99 max_value = it.GetPixel( i );
104 this->m_Seeds[ s ] = seed;
108 // End initialization
109 this->Superclass::_BeforeMainLoop( );
110 this->m_Candidates.clear( );
113 virtual void _AfterMainLoop( )
115 typedef unsigned char _TMark;
116 typedef itk::Image< _TMark, TImage::ImageDimension > _TMarkImage;
118 this->Superclass::_AfterMainLoop( );
119 if( this->m_Candidates.size( ) == 0 )
122 const TImage* input = this->GetInput( );
123 TImage::SpacingType spacing = input->GetSpacing( );
125 // Prepare an object to hold marks
126 _TMarkImage::Pointer marks = _TMarkImage::New( );
127 marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
128 marks->SetRequestedRegion( input->GetRequestedRegion( ) );
129 marks->SetBufferedRegion( input->GetBufferedRegion( ) );
130 marks->SetOrigin( input->GetOrigin( ) );
131 marks->SetSpacing( spacing );
132 marks->SetDirection( input->GetDirection( ) );
134 marks->FillBuffer( _TMark( 0 ) );
136 // Iterate over the candidates, starting from the candidates that
137 // are near thin branches
138 _TCandidates::const_reverse_iterator cIt =
139 this->m_Candidates.rbegin( );
140 for( int leo = 0; leo < 1 && cIt != this->m_Candidates.rend( ); ++cIt )
142 // If pixel hasn't been visited, continue
143 TVertex v = cIt->second;
144 if( marks->GetPixel( v ) != _TMark( 0 ) )
147 // Compute nearest start candidate
148 TImage::SizeType radius;
150 itk::ConstNeighborhoodIterator< TImage > iIt(
151 radius, input, input->GetRequestedRegion( )
153 iIt.SetLocation( v );
155 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
156 nN *= iIt.GetSize( )[ d ];
157 TVertex max_vertex = iIt.GetIndex( 0 );
158 TImage::PixelType max_value = iIt.GetPixel( 0 );
159 for( unsigned int i = 1; i < nN; ++i )
161 TImage::PixelType value = iIt.GetPixel( i );
162 if( value > max_value )
165 max_vertex = iIt.GetIndex( i );
171 if( marks->GetPixel( max_vertex ) != _TMark( 0 ) )
174 std::cout << "Leaf: " << leo << " " << max_value << " " << max_vertex << std::endl;
179 double dist = double( 1.5 ) * std::sqrt( double( input->GetPixel( max_vertex ) ) );
180 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
182 ( unsigned long )( dist / spacing[ d ] ) + 1;
183 itk::NeighborhoodIterator< _TMarkImage > mIt(
184 radius, marks, marks->GetRequestedRegion( )
186 mIt.SetLocation( max_vertex );
188 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
189 nN *= mIt.GetSize( )[ d ];
190 for( unsigned int i = 0; i < nN; ++i )
191 if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
193 mIt.SetPixel( i, ( start )? 255: 100 );
198 TImage::SizeType radius;
200 mIt.SetLocation( vIt );
202 TImage::SizeType size = mIt.GetSize( );
204 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
206 for( unsigned int i = 0; i < nN; ++i )
207 if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
208 mIt.SetPixel( i, ( start )? 255: 100 );
212 // Next vertex in current path
213 max_vertex = this->_Parent( max_vertex );
215 } while( max_vertex != this->_Parent( max_vertex ) );
220 marks->SetPixel( v, _TMark( 1 ) );
224 itk::ImageFileWriter< _TMarkImage >::Pointer w =
225 itk::ImageFileWriter< _TMarkImage >::New( );
226 w->SetInput( marks );
227 w->SetFileName( "marks.mhd" );
232 this->m_Axes = vtkSmartPointer< vtkPolyData >::New( );
233 vtkSmartPointer< vtkPoints > points =
234 vtkSmartPointer< vtkPoints >::New( );
235 vtkSmartPointer< vtkCellArray > lines =
236 vtkSmartPointer< vtkCellArray >::New( );
241 std::cout << "Leaf: " << leo << " " << cIt->first << " " << vIt << std::endl;
245 double dist = double( input->GetPixel( vIt ) );
246 TImage::SizeType radius;
247 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
249 ( unsigned long )( dist / spacing[ d ] ) + 1;
250 itk::NeighborhoodIterator< _TMarkImage > mIt(
251 radius, marks, marks->GetRequestedRegion( )
254 mIt.SetLocation( vIt );
256 TImage::SizeType size = mIt.GetSize( );
258 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
260 for( unsigned int i = 0; i < nN; ++i )
261 if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
262 mIt.SetPixel( i, ( start )? 255: 100 );
267 vIt = this->_Parent( vIt );
269 } while( vIt != this->_Parent( vIt ) );
272 // std::cout << vIt << " " << this->GetSeed( 0 ) << std::endl;
276 TVertex parent = this->_Parent( vIt );
277 while( parent != vIt )
279 pS.push_front( vIt );
281 parent = this->_Parent( vIt );
283 TImage::PointType pnt;
284 this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
285 points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
286 if( points->GetNumberOfPoints( ) > 1 )
288 lines->InsertNextCell( 2 );
289 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
290 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
295 pS.push_front( vIt );
296 TImage::PointType pnt;
297 this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
298 points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
299 if( points->GetNumberOfPoints( ) > 1 )
301 lines->InsertNextCell( 2 );
302 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
303 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
307 this->m_Axes->SetPoints( points );
308 this->m_Axes->SetLines( lines );
316 virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n )
318 TCost nc = this->_Cost( nn.Vertex, n.Vertex );
319 if( TCost( 0 ) < nc )
321 nn.Cost = n.Cost + ( TCost( 1 ) / nc );
329 virtual bool _UpdateResult( _TNode& n )
331 bool ret = this->Superclass::_UpdateResult( n );
334 TCost nc = this->_Cost( n.Vertex, n.Parent );
335 if( TCost( 0 ) < nc )
337 TCost cc = n.Cost / nc;
338 this->m_Candidates.insert( _TCandidate( cc, n.Vertex ) );
340 this->GetOutput( )->SetPixel( n.Vertex, cc );
349 TDijkstra( const Self& other );
350 Self& operator=( const Self& other );
353 _TCandidates m_Candidates;
355 vtkSmartPointer< vtkPolyData > m_Axes;
358 // -------------------------------------------------------------------------
359 int main( int argc, char* argv[] )
364 << "Usage: " << argv[ 0 ]
365 << " input_image x y z output_image"
371 std::string input_image_fn = argv[ 1 ];
372 TImage::IndexType seed_idx;
373 seed_idx[ 0 ] = std::atoi( argv[ 2 ] );
374 seed_idx[ 1 ] = std::atoi( argv[ 3 ] );
375 seed_idx[ 2 ] = std::atoi( argv[ 4 ] );
376 std::string output_image_fn = argv[ 5 ];
377 bool visual_debug = false;
379 visual_debug = ( std::atoi( argv[ 6 ] ) == 1 );
382 TImageReader::Pointer input_image_reader = TImageReader::New( );
383 input_image_reader->SetFileName( input_image_fn );
386 input_image_reader->Update( );
388 catch( itk::ExceptionObject& err )
390 std::cerr << "Error caught: " << err << std::endl;
394 TImage::ConstPointer input_image = input_image_reader->GetOutput( );
397 TVTKImage::Pointer vtk_image = TVTKImage::New( );
398 vtk_image->SetInput( input_image );
399 vtk_image->Update( );
402 view.SetBackground( 0.3, 0.2, 0.8 );
403 view.SetSize( 800, 800 );
404 view.SetImage( vtk_image->GetOutput( ) );
406 // Allow some interaction
410 TDijkstra::Pointer paths = TDijkstra::New( );
411 paths->AddSeed( seed_idx, TScalar( 0 ) );
412 paths->SetInput( input_image );
413 paths->SetNeighborhoodOrder( 1 );
417 // Configure observer
418 TDijkstraObs::Pointer obs = TDijkstraObs::New( );
419 obs->SetRenderWindow( view.GetWindow( ) );
420 paths->AddObserver( itk::AnyEvent( ), obs );
421 paths->ThrowEventsOn( );
424 paths->ThrowEventsOff( );
425 std::clock_t start = std::clock( );
427 std::clock_t end = std::clock( );
428 double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
429 std::cout << "Paths extraction time = " << seconds << std::endl;
431 view.AddPolyData( paths->m_Axes, 1, 0, 0 );
434 TImageWriter::Pointer dijkstra_writer =
435 TImageWriter::New( );
436 dijkstra_writer->SetInput( paths->GetOutput( ) );
437 dijkstra_writer->SetFileName( "dijkstra.mhd" );
438 dijkstra_writer->Update( );
442 TVTKImage::Pointer output_image_vtk = TVTKImage::New( );
443 output_image_vtk->SetInput( segmentor->GetOutput( ) );
444 output_image_vtk->Update( );
446 vtkSmartPointer< vtkImageMarchingCubes > mc =
447 vtkSmartPointer< vtkImageMarchingCubes >::New( );
448 mc->SetInputData( output_image_vtk->GetOutput( ) );
451 double( segmentor->GetInsideValue( ) ) * double( 0.95 )
455 // Let some interaction and close program
456 view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 );
459 // Write resulting image
460 TImageWriter::Pointer output_image_writer = TImageWriter::New( );
461 output_image_writer->SetInput( segmentor->GetOutput( ) );
462 output_image_writer->SetFileName( output_image_fn );
465 output_image_writer->Update( );
467 catch( itk::ExceptionObject& err )
469 std::cerr << "Error caught: " << err << std::endl;