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 itk::Image< unsigned char, TImage::ImageDimension > _TMarkImage;
117 this->Superclass::_AfterMainLoop( );
118 if( this->m_Candidates.size( ) == 0 )
121 const TImage* input = this->GetInput( );
122 TImage::SpacingType spacing = input->GetSpacing( );
124 _TMarkImage::Pointer marks = _TMarkImage::New( );
125 marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
126 marks->SetRequestedRegion( input->GetRequestedRegion( ) );
127 marks->SetBufferedRegion( input->GetBufferedRegion( ) );
128 marks->SetOrigin( input->GetOrigin( ) );
129 marks->SetSpacing( spacing );
130 marks->SetDirection( input->GetDirection( ) );
132 marks->FillBuffer( 0 );
134 this->m_Axes = vtkSmartPointer< vtkPolyData >::New( );
135 vtkSmartPointer< vtkPoints > points =
136 vtkSmartPointer< vtkPoints >::New( );
137 vtkSmartPointer< vtkCellArray > lines =
138 vtkSmartPointer< vtkCellArray >::New( );
140 _TCandidates::const_reverse_iterator cIt =
141 this->m_Candidates.rbegin( );
143 for( unsigned long leo = 0; leo < 50 && cIt != this->m_Candidates.rend( ); ++cIt )
145 TVertex vIt = cIt->second;
146 if( marks->GetPixel( vIt ) != 0 )
150 std::cout << "Leaf: " << leo << " " << cIt->first << " " << vIt << std::endl;
154 double dist = double( input->GetPixel( vIt ) );
155 TImage::SizeType radius;
156 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
158 ( unsigned long )( dist / spacing[ d ] ) + 1;
159 itk::NeighborhoodIterator< _TMarkImage > mIt(
160 radius, marks, marks->GetRequestedRegion( )
163 mIt.SetLocation( vIt );
165 TImage::SizeType size = mIt.GetSize( );
167 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
169 for( unsigned int i = 0; i < nN; ++i )
170 if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
171 mIt.SetPixel( i, ( start )? 255: 100 );
176 vIt = this->_Parent( vIt );
178 } while( vIt != this->_Parent( vIt ) );
181 // std::cout << vIt << " " << this->GetSeed( 0 ) << std::endl;
185 TVertex parent = this->_Parent( vIt );
186 while( parent != vIt )
188 pS.push_front( vIt );
190 parent = this->_Parent( vIt );
192 TImage::PointType pnt;
193 this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
194 points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
195 if( points->GetNumberOfPoints( ) > 1 )
197 lines->InsertNextCell( 2 );
198 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
199 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
204 pS.push_front( vIt );
205 TImage::PointType pnt;
206 this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
207 points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
208 if( points->GetNumberOfPoints( ) > 1 )
210 lines->InsertNextCell( 2 );
211 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
212 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
216 this->m_Axes->SetPoints( points );
217 this->m_Axes->SetLines( lines );
222 itk::ImageFileWriter< _TMarkImage >::Pointer w =
223 itk::ImageFileWriter< _TMarkImage >::New( );
224 w->SetInput( marks );
225 w->SetFileName( "marks.mhd" );
230 virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n )
232 TCost nc = this->_Cost( nn.Vertex, n.Vertex );
233 if( TCost( 0 ) < nc )
235 nn.Cost = n.Cost + ( TCost( 1 ) / nc );
243 virtual bool _UpdateResult( _TNode& n )
245 bool ret = this->Superclass::_UpdateResult( n );
248 TCost nc = this->_Cost( n.Vertex, n.Parent );
249 TCost nc2 = nc * nc * nc * nc;
250 if( TCost( 0 ) < nc2 )
252 n.Result = n.Cost / nc2;
253 this->GetOutput( )->SetPixel( n.Vertex, n.Result );
254 this->m_Candidates.insert( _TCandidate( n.Result, n.Vertex ) );
263 TDijkstra( const Self& other );
264 Self& operator=( const Self& other );
267 _TCandidates m_Candidates;
269 vtkSmartPointer< vtkPolyData > m_Axes;
272 // -------------------------------------------------------------------------
273 int main( int argc, char* argv[] )
278 << "Usage: " << argv[ 0 ]
279 << " input_image x y z output_image"
285 std::string input_image_fn = argv[ 1 ];
286 TImage::IndexType seed_idx;
287 seed_idx[ 0 ] = std::atoi( argv[ 2 ] );
288 seed_idx[ 1 ] = std::atoi( argv[ 3 ] );
289 seed_idx[ 2 ] = std::atoi( argv[ 4 ] );
290 std::string output_image_fn = argv[ 5 ];
291 bool visual_debug = false;
293 visual_debug = ( std::atoi( argv[ 6 ] ) == 1 );
296 TImageReader::Pointer input_image_reader = TImageReader::New( );
297 input_image_reader->SetFileName( input_image_fn );
300 input_image_reader->Update( );
302 catch( itk::ExceptionObject& err )
304 std::cerr << "Error caught: " << err << std::endl;
308 TImage::ConstPointer input_image = input_image_reader->GetOutput( );
311 TVTKImage::Pointer vtk_image = TVTKImage::New( );
312 vtk_image->SetInput( input_image );
313 vtk_image->Update( );
316 view.SetBackground( 0.3, 0.2, 0.8 );
317 view.SetSize( 800, 800 );
318 view.SetImage( vtk_image->GetOutput( ) );
320 // Allow some interaction
324 TDijkstra::Pointer paths = TDijkstra::New( );
325 paths->AddSeed( seed_idx, TScalar( 0 ) );
326 paths->SetInput( input_image );
327 paths->SetNeighborhoodOrder( 1 );
331 // Configure observer
332 TDijkstraObs::Pointer obs = TDijkstraObs::New( );
333 obs->SetImage( input_image, view.GetWindow( ) );
334 paths->AddObserver( itk::AnyEvent( ), obs );
335 paths->ThrowEventsOn( );
338 paths->ThrowEventsOff( );
339 std::clock_t start = std::clock( );
341 std::clock_t end = std::clock( );
342 double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
343 std::cout << "Paths extraction time = " << seconds << std::endl;
345 view.AddPolyData( paths->m_Axes, 1, 0, 0 );
348 TImageWriter::Pointer dijkstra_writer =
349 TImageWriter::New( );
350 dijkstra_writer->SetInput( paths->GetOutput( ) );
351 dijkstra_writer->SetFileName( "dijkstra.mhd" );
352 dijkstra_writer->Update( );
356 TVTKImage::Pointer output_image_vtk = TVTKImage::New( );
357 output_image_vtk->SetInput( segmentor->GetOutput( ) );
358 output_image_vtk->Update( );
360 vtkSmartPointer< vtkImageMarchingCubes > mc =
361 vtkSmartPointer< vtkImageMarchingCubes >::New( );
362 mc->SetInputData( output_image_vtk->GetOutput( ) );
365 double( segmentor->GetInsideValue( ) ) * double( 0.95 )
369 // Let some interaction and close program
370 view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 );
373 // Write resulting image
374 TImageWriter::Pointer output_image_writer = TImageWriter::New( );
375 output_image_writer->SetInput( segmentor->GetOutput( ) );
376 output_image_writer->SetFileName( output_image_fn );
379 output_image_writer->Update( );
381 catch( itk::ExceptionObject& err )
383 std::cerr << "Error caught: " << err << std::endl;