8 #include <itkConstNeighborhoodIterator.h>
9 #include <itkDanielssonDistanceMapImageFilter.h>
11 #include <itkImageFileReader.h>
12 #include <itkImageFileWriter.h>
13 #include <itkImageToVTKImageFilter.h>
15 #include <vtkPoints.h>
16 #include <vtkCellArray.h>
17 #include <vtkPolyData.h>
18 #include <vtkSmartPointer.h>
20 #include <fpa/Base/Functors/InvertCostFunction.h>
21 #include <fpa/Image/Dijkstra.h>
22 #include <fpa/Image/RegionGrowWithMultipleThresholds.h>
23 #include <fpa/VTK/ImageMPR.h>
24 #include <fpa/VTK/Image3DObserver.h>
26 // -------------------------------------------------------------------------
27 const unsigned int Dim = 3;
28 typedef double TPixel;
29 typedef double TScalar;
30 typedef itk::Image< TPixel, Dim > TImage;
31 typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
33 typedef itk::ImageFileReader< TImage > TImageReader;
34 typedef itk::ImageFileWriter< TImage > TImageWriter;
35 typedef fpa::Image::Dijkstra< TImage, TScalar > TBaseDijkstra;
37 typedef fpa::VTK::ImageMPR TMPR;
38 typedef fpa::VTK::Image3DObserver< TBaseDijkstra, vtkRenderWindow > TDijkstraObs;
40 // -------------------------------------------------------------------------
42 : public TBaseDijkstra
45 typedef TDijkstra Self;
46 typedef TBaseDijkstra Superclass;
47 typedef itk::SmartPointer< Self > Pointer;
48 typedef itk::SmartPointer< const Self > ConstPointer;
50 typedef Superclass::TCost TCost;
51 typedef Superclass::TVertex TVertex;
52 typedef Superclass::InputImageType TImage;
53 typedef std::deque< TVertex > TVertices;
56 typedef std::pair< TCost, TVertex > _TCandidate;
57 typedef std::multimap< TCost, TVertex > _TCandidates;
58 typedef Superclass::_TNode _TNode;
62 itkTypeMacro( TDijkstra, TBaseDijkstra );
71 virtual void _BeforeMainLoop( )
73 const TImage* img = this->GetInput( );
76 TImage::SizeType radius;
78 itk::ConstNeighborhoodIterator< TImage > it(
79 radius, img, img->GetRequestedRegion( )
81 for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
83 _TNode seed = this->m_Seeds[ s ];
84 it.SetLocation( seed.Vertex );
86 TImage::SizeType size = it.GetSize( );
88 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
90 TImage::PixelType max_value = img->GetPixel( seed.Vertex );
91 for( unsigned int i = 0; i < nN; ++i )
93 if( it.GetPixel( i ) > max_value )
95 seed.Vertex = it.GetIndex( i );
96 seed.Parent = seed.Vertex;
97 max_value = it.GetPixel( i );
102 this->m_Seeds[ s ] = seed;
106 // End initialization
107 this->Superclass::_BeforeMainLoop( );
108 this->m_Candidates.clear( );
111 virtual void _AfterMainLoop( )
113 this->Superclass::_AfterMainLoop( );
114 if( this->m_Candidates.size( ) == 0 )
117 this->m_Axes = vtkSmartPointer< vtkPolyData >::New( );
118 vtkSmartPointer< vtkPoints > points =
119 vtkSmartPointer< vtkPoints >::New( );
120 vtkSmartPointer< vtkCellArray > lines =
121 vtkSmartPointer< vtkCellArray >::New( );
123 _TCandidates::const_reverse_iterator cIt =
124 this->m_Candidates.rbegin( );
127 TVertex vIt = cIt->second;
128 TVertex parent = this->_Parent( vIt );
129 while( parent != vIt )
131 pS.push_front( vIt );
133 parent = this->_Parent( vIt );
135 TImage::PointType pnt;
136 this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
137 points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
138 if( points->GetNumberOfPoints( ) > 1 )
140 lines->InsertNextCell( 2 );
141 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
142 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
147 pS.push_front( vIt );
148 TImage::PointType pnt;
149 this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
150 points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
151 if( points->GetNumberOfPoints( ) > 1 )
153 lines->InsertNextCell( 2 );
154 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
155 lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
159 this->m_Axes->SetPoints( points );
160 this->m_Axes->SetLines( lines );
163 virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n )
165 TCost nc = this->_Cost( nn.Vertex, n.Vertex );
166 if( TCost( 0 ) < nc )
168 nn.Cost = n.Cost + ( TCost( 1 ) / nc );
176 virtual bool _UpdateResult( _TNode& n )
178 bool ret = this->Superclass::_UpdateResult( n );
181 TCost nc = this->_Cost( n.Vertex, n.Parent );
182 if( TCost( 0 ) < nc )
184 n.Result = n.Cost / ( nc * nc );
185 this->GetOutput( )->SetPixel( n.Vertex, n.Result );
186 this->m_Candidates.insert( _TCandidate( n.Result, n.Vertex ) );
195 TDijkstra( const Self& other );
196 Self& operator=( const Self& other );
199 _TCandidates m_Candidates;
201 vtkSmartPointer< vtkPolyData > m_Axes;
204 // -------------------------------------------------------------------------
205 int main( int argc, char* argv[] )
210 << "Usage: " << argv[ 0 ]
211 << " input_image x y z output_image"
217 std::string input_image_fn = argv[ 1 ];
218 TImage::IndexType seed_idx;
219 seed_idx[ 0 ] = std::atoi( argv[ 2 ] );
220 seed_idx[ 1 ] = std::atoi( argv[ 3 ] );
221 seed_idx[ 2 ] = std::atoi( argv[ 4 ] );
222 std::string output_image_fn = argv[ 5 ];
223 bool visual_debug = false;
225 visual_debug = ( std::atoi( argv[ 6 ] ) == 1 );
228 TImageReader::Pointer input_image_reader = TImageReader::New( );
229 input_image_reader->SetFileName( input_image_fn );
232 input_image_reader->Update( );
234 catch( itk::ExceptionObject& err )
236 std::cerr << "Error caught: " << err << std::endl;
240 TImage::ConstPointer input_image = input_image_reader->GetOutput( );
243 TVTKImage::Pointer vtk_image = TVTKImage::New( );
244 vtk_image->SetInput( input_image );
245 vtk_image->Update( );
248 view.SetBackground( 0.3, 0.2, 0.8 );
249 view.SetSize( 800, 800 );
250 view.SetImage( vtk_image->GetOutput( ) );
252 // Allow some interaction
256 TDijkstra::Pointer paths = TDijkstra::New( );
257 paths->AddSeed( seed_idx, TScalar( 0 ) );
258 paths->SetInput( input_image );
259 paths->SetNeighborhoodOrder( 1 );
263 // Configure observer
264 TDijkstraObs::Pointer obs = TDijkstraObs::New( );
265 obs->SetImage( input_image, view.GetWindow( ) );
266 paths->AddObserver( itk::AnyEvent( ), obs );
267 paths->ThrowEventsOn( );
270 paths->ThrowEventsOff( );
271 std::clock_t start = std::clock( );
273 std::clock_t end = std::clock( );
274 double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
275 std::cout << "Paths extraction time = " << seconds << std::endl;
277 view.AddPolyData( paths->m_Axes, 1, 0, 0 );
280 TImageWriter::Pointer dijkstra_writer =
281 TImageWriter::New( );
282 dijkstra_writer->SetInput( paths->GetOutput( ) );
283 dijkstra_writer->SetFileName( "dijkstra.mhd" );
284 dijkstra_writer->Update( );
288 TVTKImage::Pointer output_image_vtk = TVTKImage::New( );
289 output_image_vtk->SetInput( segmentor->GetOutput( ) );
290 output_image_vtk->Update( );
292 vtkSmartPointer< vtkImageMarchingCubes > mc =
293 vtkSmartPointer< vtkImageMarchingCubes >::New( );
294 mc->SetInputData( output_image_vtk->GetOutput( ) );
297 double( segmentor->GetInsideValue( ) ) * double( 0.95 )
301 // Let some interaction and close program
302 view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 );
305 // Write resulting image
306 TImageWriter::Pointer output_image_writer = TImageWriter::New( );
307 output_image_writer->SetInput( segmentor->GetOutput( ) );
308 output_image_writer->SetFileName( output_image_fn );
311 output_image_writer->Update( );
313 catch( itk::ExceptionObject& err )
315 std::cerr << "Error caught: " << err << std::endl;