5 #include <itkConstNeighborhoodIterator.h>
6 #include <itkDanielssonDistanceMapImageFilter.h>
8 #include <itkImageFileReader.h>
9 #include <itkImageFileWriter.h>
10 #include <itkImageToVTKImageFilter.h>
12 #include <fpa/Image/RegionGrowWithMultipleThresholds.h>
13 #include <fpa/VTK/ImageMPR.h>
14 #include <fpa/VTK/Image3DObserver.h>
16 // -------------------------------------------------------------------------
17 const unsigned int Dim = 3;
19 typedef double TScalar;
20 typedef itk::Image< TPixel, Dim > TImage;
21 typedef itk::Image< TScalar, Dim > TDistanceMap;
22 typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
24 typedef itk::ImageFileReader< TImage > TImageReader;
25 typedef itk::ImageFileWriter< TImage > TImageWriter;
26 typedef itk::ImageFileWriter< TDistanceMap > TDistanceMapWriter;
27 typedef fpa::Image::RegionGrowWithMultipleThresholds< TImage > TSegmentor;
29 itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >
32 typedef fpa::VTK::ImageMPR TMPR;
33 typedef fpa::VTK::Image3DObserver< TSegmentor, vtkRenderWindow > TSegmentorObs;
35 // -------------------------------------------------------------------------
36 int main( int argc, char* argv[] )
41 << "Usage: " << argv[ 0 ]
42 << " input_image thr_0 thr_1 step"
43 << " output_segmentation output_distancemap"
49 std::string input_image_fn = argv[ 1 ];
50 TPixel thr_0 = TPixel( std::atof( argv[ 2 ] ) );
51 TPixel thr_1 = TPixel( std::atof( argv[ 3 ] ) );
52 unsigned int step = std::atoi( argv[ 4 ] );
53 std::string output_segmentation_fn = argv[ 5 ];
54 std::string output_distancemap_fn = argv[ 6 ];
55 bool visual_debug = false;
57 visual_debug = ( std::atoi( argv[ 7 ] ) == 1 );
60 TImageReader::Pointer input_image_reader = TImageReader::New( );
61 input_image_reader->SetFileName( input_image_fn );
64 input_image_reader->Update( );
66 catch( itk::ExceptionObject& err )
68 std::cerr << "Error caught: " << err << std::endl;
72 TImage::ConstPointer input_image = input_image_reader->GetOutput( );
75 TVTKImage::Pointer vtk_image = TVTKImage::New( );
76 vtk_image->SetInput( input_image );
80 view.SetBackground( 0.3, 0.2, 0.8 );
81 view.SetSize( 800, 800 );
82 view.SetImage( vtk_image->GetOutput( ) );
84 // Wait for a seed to be given
85 while( view.GetNumberOfSeeds( ) == 0 )
90 view.GetSeed( 0, seed );
91 TImage::PointType seed_pnt;
92 seed_pnt[ 0 ] = seed[ 0 ];
93 seed_pnt[ 1 ] = seed[ 1 ];
94 seed_pnt[ 2 ] = seed[ 2 ];
95 TImage::IndexType seed_idx;
96 input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx );
98 // Segment input image
99 TSegmentor::Pointer segmentor = TSegmentor::New( );
100 segmentor->AddThresholds( thr_0, thr_1, step );
101 segmentor->AddSeed( seed_idx, 0 );
102 segmentor->SetInput( input_image );
103 segmentor->SetNeighborhoodOrder( 1 );
104 segmentor->SetDifferenceThreshold( double( 3 ) );
105 segmentor->SetInsideValue( TPixel( 0 ) ); // WARNING: This is to optimize
106 segmentor->SetOutsideValue( TPixel( 1 ) ); // distance map calculation
109 // Configure observer
110 TSegmentorObs::Pointer obs = TSegmentorObs::New( );
111 obs->SetImage( input_image, view.GetWindow( ) );
112 segmentor->AddObserver( itk::AnyEvent( ), obs );
113 segmentor->ThrowEventsOn( );
116 segmentor->ThrowEventsOff( );
117 std::clock_t start = std::clock( );
118 segmentor->Update( );
119 std::clock_t end = std::clock( );
120 double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
121 std::cout << "Segmentation time = " << seconds << std::endl;
123 // Extract distance map
124 TDistanceFilter::Pointer distanceMap = TDistanceFilter::New( );
125 distanceMap->SetInput( segmentor->GetOutput( ) );
126 distanceMap->InputIsBinaryOn( );
127 distanceMap->SquaredDistanceOn( );
128 start = std::clock( );
129 distanceMap->Update( );
131 seconds = double( end - start ) / double( CLOCKS_PER_SEC );
132 std::cout << "Distance map time = " << seconds << std::endl;
134 std::cout << "Final seed: " << seed_idx << std::endl;
136 TDistanceMapWriter::Pointer distancemap_writer =
137 TDistanceMapWriter::New( );
138 distancemap_writer->SetInput( distanceMap->GetOutput( ) );
139 distancemap_writer->SetFileName( output_distancemap_fn );
140 distancemap_writer->Update( );
142 TImageWriter::Pointer segmentation_writer =
143 TImageWriter::New( );
144 segmentation_writer->SetInput( segmentor->GetOutput( ) );
145 segmentation_writer->SetFileName( output_segmentation_fn );
146 segmentation_writer->Update( );
150 TVTKImage::Pointer output_image_vtk = TVTKImage::New( );
151 output_image_vtk->SetInput( segmentor->GetOutput( ) );
152 output_image_vtk->Update( );
154 vtkSmartPointer< vtkImageMarchingCubes > mc =
155 vtkSmartPointer< vtkImageMarchingCubes >::New( );
156 mc->SetInputData( output_image_vtk->GetOutput( ) );
159 double( segmentor->GetInsideValue( ) ) * double( 0.95 )
163 // Let some interaction and close program
164 view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 );
167 // Write resulting image
168 TImageWriter::Pointer output_image_writer = TImageWriter::New( );
169 output_image_writer->SetInput( segmentor->GetOutput( ) );
170 output_image_writer->SetFileName( output_image_fn );
173 output_image_writer->Update( );
175 catch( itk::ExceptionObject& err )
177 std::cerr << "Error caught: " << err << std::endl;