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;
57 typedef Superclass::TEndEvent TEndEvent;
58 typedef Superclass::TBacktrackingEvent TBacktrackingEvent;
61 typedef std::pair< TCost, TVertex > _TCandidate;
62 typedef std::multimap< TCost, TVertex > _TCandidates;
63 typedef Superclass::_TNode _TNode;
67 itkTypeMacro( TDijkstra, TBaseDijkstra );
76 virtual void _BeforeMainLoop( )
78 const TImage* img = this->GetInput( );
81 TImage::SizeType radius;
83 itk::ConstNeighborhoodIterator< TImage > it(
84 radius, img, img->GetRequestedRegion( )
86 for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
88 _TNode seed = this->m_Seeds[ s ];
89 it.SetLocation( seed.Vertex );
91 TImage::SizeType size = it.GetSize( );
93 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
95 TImage::PixelType max_value = img->GetPixel( seed.Vertex );
96 for( unsigned int i = 0; i < nN; ++i )
98 if( it.GetPixel( i ) > max_value )
100 seed.Vertex = it.GetIndex( i );
101 seed.Parent = seed.Vertex;
102 max_value = it.GetPixel( i );
107 this->m_Seeds[ s ] = seed;
111 // End initialization
112 this->Superclass::_BeforeMainLoop( );
113 this->m_Candidates.clear( );
116 virtual void _AfterMainLoop( )
118 typedef unsigned char _TMark;
119 typedef itk::Image< _TMark, TImage::ImageDimension > _TMarkImage;
121 this->Superclass::_AfterMainLoop( );
122 if( this->m_Candidates.size( ) == 0 )
125 this->InvokeEvent( TEndEvent( ) );
127 const TImage* input = this->GetInput( );
128 TImage::SpacingType spacing = input->GetSpacing( );
130 // Prepare an object to hold marks
131 _TMarkImage::Pointer marks = _TMarkImage::New( );
132 marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
133 marks->SetRequestedRegion( input->GetRequestedRegion( ) );
134 marks->SetBufferedRegion( input->GetBufferedRegion( ) );
135 marks->SetOrigin( input->GetOrigin( ) );
136 marks->SetSpacing( spacing );
137 marks->SetDirection( input->GetDirection( ) );
139 marks->FillBuffer( _TMark( 0 ) );
141 // Iterate over the candidates, starting from the candidates that
142 // are near thin branches
143 _TCandidates::const_reverse_iterator cIt =
144 this->m_Candidates.rbegin( );
145 for( int backId = 0; backId < 100 && cIt != this->m_Candidates.rend( ); ++cIt )
147 // If pixel hasn't been visited, continue
148 TVertex v = cIt->second;
149 if( marks->GetPixel( v ) != _TMark( 0 ) )
152 // Compute nearest start candidate
153 TImage::SizeType radius;
155 itk::ConstNeighborhoodIterator< TImage > iIt(
156 radius, input, input->GetRequestedRegion( )
158 iIt.SetLocation( v );
160 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
161 nN *= iIt.GetSize( )[ d ];
162 TVertex max_vertex = iIt.GetIndex( 0 );
163 TImage::PixelType max_value = iIt.GetPixel( 0 );
164 for( unsigned int i = 1; i < nN; ++i )
166 TImage::PixelType value = iIt.GetPixel( i );
167 if( value > max_value )
170 max_vertex = iIt.GetIndex( i );
176 if( marks->GetPixel( max_vertex ) != _TMark( 0 ) )
179 std::cout << "Leaf: " << backId << " " << max_value << " " << max_vertex << std::endl;
184 double dist = double( 1.5 ) * std::sqrt( double( input->GetPixel( max_vertex ) ) );
185 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
187 ( unsigned long )( dist / spacing[ d ] ) + 1;
188 itk::NeighborhoodIterator< _TMarkImage > mIt(
189 radius, marks, marks->GetRequestedRegion( )
191 mIt.SetLocation( max_vertex );
193 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
194 nN *= mIt.GetSize( )[ d ];
195 for( unsigned int i = 0; i < nN; ++i )
196 if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
198 mIt.SetPixel( i, ( start )? 255: 100 );
203 TImage::SizeType radius;
205 mIt.SetLocation( vIt );
207 TImage::SizeType size = mIt.GetSize( );
209 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
211 for( unsigned int i = 0; i < nN; ++i )
212 if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
213 mIt.SetPixel( i, ( start )? 255: 100 );
217 // Next vertex in current path
218 this->InvokeEvent( TBacktrackingEvent( max_vertex, backId ) );
219 max_vertex = this->_Parent( max_vertex );
221 } while( max_vertex != this->_Parent( max_vertex ) );
226 marks->SetPixel( v, _TMark( 1 ) );
231 itk::ImageFileWriter< _TMarkImage >::Pointer w =
232 itk::ImageFileWriter< _TMarkImage >::New( );
233 w->SetInput( marks );
234 w->SetFileName( "marks.mhd" );
238 this->m_Axes = vtkSmartPointer< vtkPolyData >::New( );
239 vtkSmartPointer< vtkPoints > points =
240 vtkSmartPointer< vtkPoints >::New( );
241 vtkSmartPointer< vtkCellArray > lines =
242 vtkSmartPointer< vtkCellArray >::New( );
247 std::cout << "Leaf: " << backId << " " << cIt->first << " " << vIt << std::endl;
251 double dist = double( input->GetPixel( vIt ) );
252 TImage::SizeType radius;
253 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
255 ( unsigned long )( dist / spacing[ d ] ) + 1;
256 itk::NeighborhoodIterator< _TMarkImage > mIt(
257 radius, marks, marks->GetRequestedRegion( )
260 mIt.SetLocation( vIt );
262 TImage::SizeType size = mIt.GetSize( );
264 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
266 for( unsigned int i = 0; i < nN; ++i )
267 if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
268 mIt.SetPixel( i, ( start )? 255: 100 );
273 vIt = this->_Parent( vIt );
275 } while( vIt != this->_Parent( vIt ) );
278 // std::cout << vIt << " " << this->GetSeed( 0 ) << std::endl;
282 TVertex parent = this->_Parent( vIt );
283 while( parent != vIt )
285 pS.push_front( vIt );
287 parent = this->_Parent( vIt );
289 TImage::PointType pnt;
290 this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
291 points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
292 if( points->GetNumberOfPoints( ) > 1 )
294 lines->InsertNextCell( 2 );
295 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
296 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
301 pS.push_front( vIt );
302 TImage::PointType pnt;
303 this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
304 points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
305 if( points->GetNumberOfPoints( ) > 1 )
307 lines->InsertNextCell( 2 );
308 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
309 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
313 this->m_Axes->SetPoints( points );
314 this->m_Axes->SetLines( lines );
322 virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n )
324 TCost nc = this->_Cost( nn.Vertex, n.Vertex );
325 if( TCost( 0 ) < nc )
327 nn.Cost = n.Cost + ( TCost( 1 ) / nc );
335 virtual bool _UpdateResult( _TNode& n )
337 bool ret = this->Superclass::_UpdateResult( n );
340 TCost nc = this->_Cost( n.Vertex, n.Parent );
341 if( TCost( 0 ) < nc )
343 TCost cc = n.Cost / nc;
344 this->m_Candidates.insert( _TCandidate( cc, n.Vertex ) );
346 this->GetOutput( )->SetPixel( n.Vertex, cc );
355 TDijkstra( const Self& other );
356 Self& operator=( const Self& other );
359 _TCandidates m_Candidates;
361 vtkSmartPointer< vtkPolyData > m_Axes;
364 // -------------------------------------------------------------------------
365 int main( int argc, char* argv[] )
370 << "Usage: " << argv[ 0 ]
371 << " input_image x y z output_image"
377 std::string input_image_fn = argv[ 1 ];
378 TImage::IndexType seed_idx;
379 seed_idx[ 0 ] = std::atoi( argv[ 2 ] );
380 seed_idx[ 1 ] = std::atoi( argv[ 3 ] );
381 seed_idx[ 2 ] = std::atoi( argv[ 4 ] );
382 std::string output_image_fn = argv[ 5 ];
383 bool visual_debug = false;
385 visual_debug = ( std::atoi( argv[ 6 ] ) == 1 );
388 TImageReader::Pointer input_image_reader = TImageReader::New( );
389 input_image_reader->SetFileName( input_image_fn );
392 input_image_reader->Update( );
394 catch( itk::ExceptionObject& err )
396 std::cerr << "Error caught: " << err << std::endl;
400 TImage::ConstPointer input_image = input_image_reader->GetOutput( );
403 TVTKImage::Pointer vtk_image = TVTKImage::New( );
404 vtk_image->SetInput( input_image );
405 vtk_image->Update( );
408 view.SetBackground( 0.3, 0.2, 0.8 );
409 view.SetSize( 800, 800 );
410 view.SetImage( vtk_image->GetOutput( ) );
412 // Allow some interaction
416 TDijkstra::Pointer paths = TDijkstra::New( );
417 paths->AddSeed( seed_idx, TScalar( 0 ) );
418 paths->SetInput( input_image );
419 paths->SetNeighborhoodOrder( 1 );
423 // Configure observer
424 TDijkstraObs::Pointer obs = TDijkstraObs::New( );
425 obs->SetRenderWindow( view.GetWindow( ) );
426 paths->AddObserver( itk::AnyEvent( ), obs );
427 paths->ThrowEventsOn( );
430 paths->ThrowEventsOff( );
431 std::clock_t start = std::clock( );
433 std::clock_t end = std::clock( );
434 double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
435 std::cout << "Paths extraction time = " << seconds << std::endl;
437 view.AddPolyData( paths->m_Axes, 1, 0, 0 );
440 TImageWriter::Pointer dijkstra_writer =
441 TImageWriter::New( );
442 dijkstra_writer->SetInput( paths->GetOutput( ) );
443 dijkstra_writer->SetFileName( "dijkstra.mhd" );
444 dijkstra_writer->Update( );
448 TVTKImage::Pointer output_image_vtk = TVTKImage::New( );
449 output_image_vtk->SetInput( segmentor->GetOutput( ) );
450 output_image_vtk->Update( );
452 vtkSmartPointer< vtkImageMarchingCubes > mc =
453 vtkSmartPointer< vtkImageMarchingCubes >::New( );
454 mc->SetInputData( output_image_vtk->GetOutput( ) );
457 double( segmentor->GetInsideValue( ) ) * double( 0.95 )
461 // Let some interaction and close program
462 view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 );
465 // Write resulting image
466 TImageWriter::Pointer output_image_writer = TImageWriter::New( );
467 output_image_writer->SetInput( segmentor->GetOutput( ) );
468 output_image_writer->SetFileName( output_image_fn );
471 output_image_writer->Update( );
473 catch( itk::ExceptionObject& err )
475 std::cerr << "Error caught: " << err << std::endl;