8 #include <itkImageToVTKImageFilter.h>
10 #include <itkImageFileReader.h>
11 #include <itkMinimumMaximumImageCalculator.h>
12 #include <itkShiftScaleImageFilter.h>
14 #include <itkSignedMaurerDistanceMapImageFilter.h>
16 #include <itkImageFileWriter.h>
18 #include <fpa/Image/DijkstraWithEndPointDetection.h>
19 #include <fpa/Base/Functors/InvertCostFunction.h>
20 #include <fpa/VTK/ImageMPR.h>
21 #include <fpa/VTK/Image3DObserver.h>
24 #include <itkMinimumMaximumImageCalculator.h>
25 #include <itkInvertIntensityImageFilter.h>
26 #include <itkDanielssonDistanceMapImageFilter.h>
27 #include <itkImageFileWriter.h>
29 #include <vtkCamera.h>
30 #include <vtkImageActor.h>
31 #include <vtkInteractorStyleImage.h>
32 #include <vtkPointHandleRepresentation3D.h>
33 #include <vtkProperty.h>
34 #include <vtkRenderer.h>
35 #include <vtkRenderWindow.h>
36 #include <vtkRenderWindowInteractor.h>
37 #include <vtkSeedRepresentation.h>
38 #include <vtkSeedWidget.h>
39 #include <vtkSmartPointer.h>
41 #include <fpa/Image/Dijkstra.h>
44 // -------------------------------------------------------------------------
45 const unsigned int Dim = 3;
46 typedef unsigned char TPixel;
47 typedef float TScalar;
49 typedef itk::Image< TPixel, Dim > TImage;
50 typedef itk::Image< TScalar, Dim > TScalarImage;
51 typedef itk::ImageToVTKImageFilter< TImage > TVTKInputImage;
53 // -------------------------------------------------------------------------
55 void ReadImage( typename I::Pointer& image, const std::string& filename );
57 template< class I, class O >
59 const typename I::Pointer& input, typename O::Pointer& output
62 // -------------------------------------------------------------------------
63 int main( int argc, char* argv[] )
68 << "Usage: " << argv[ 0 ]
74 std::string input_image_fn = argv[ 1 ];
77 TImage::Pointer input_image;
80 ReadImage< TImage >( input_image, input_image_fn );
82 catch( itk::ExceptionObject& err )
85 << "Error caught while reading \""
86 << input_image_fn << "\": " << err
91 TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( );
92 vtk_input_image->SetInput( input_image );
93 vtk_input_image->Update( );
95 // Show input image and let some interaction
96 fpa::VTK::ImageMPR view;
97 view.SetBackground( 0.3, 0.2, 0.8 );
98 view.SetSize( 800, 800 );
99 view.SetImage( vtk_input_image->GetOutput( ) );
101 // Allow some interaction and wait for, at least, one seed
103 while( view.GetNumberOfSeeds( ) == 0 )
108 view.GetSeed( 0, p );
109 TImage::PointType pnt;
110 TImage::IndexType seed;
111 pnt[ 0 ] = TImage::PointType::ValueType( p[ 0 ] );
112 pnt[ 1 ] = TImage::PointType::ValueType( p[ 1 ] );
113 pnt[ 2 ] = TImage::PointType::ValueType( p[ 2 ] );
114 input_image->TransformPhysicalPointToIndex( pnt, seed );
116 // Compute squared distance map
117 TScalarImage::Pointer dmap;
118 DistanceMap< TImage, TScalarImage >( input_image, dmap );
120 // Prepare cost conversion function
121 typedef fpa::Base::Functors::InvertCostFunction< TScalar > TFunction;
122 TFunction::Pointer function = TFunction::New( );
124 // Prepare Dijkstra filter
125 typedef fpa::Image::DijkstraWithEndPointDetection< TScalarImage, TScalarImage > TFilter;
126 TFilter::Pointer filter = TFilter::New( );
127 filter->SetInput( dmap );
128 filter->SetConversionFunction( function );
129 filter->SetNeighborhoodOrder( 2 );
130 filter->StopAtOneFrontOff( );
131 filter->AddSeed( seed, TScalar( 0 ) );
133 // Prepare graphical debugger
134 typedef fpa::VTK::Image3DObserver< TFilter, vtkRenderWindow > TDebugger;
135 TDebugger::Pointer debugger = TDebugger::New( );
136 debugger->SetRenderWindow( view.GetWindow( ) );
137 debugger->SetRenderPercentage( 0.0001 );
138 filter->AddObserver( itk::AnyEvent( ), debugger );
139 filter->ThrowEventsOn( );
142 std::time_t start, end;
147 << "Extraction time = "
148 << std::difftime( end, start )
149 << " s." << std::endl;
152 // Save final total cost map
153 itk::ImageFileWriter< TOutputImage >::Pointer output_image_writer =
154 itk::ImageFileWriter< TOutputImage >::New( );
155 output_image_writer->SetFileName( output_image_fn );
156 output_image_writer->SetInput( filter->GetOutput( ) );
159 output_image_writer->Update( );
161 catch( itk::ExceptionObject& err )
164 << "Error while writing image to " << output_image_fn << ": "
173 // -------------------------------------------------------------------------
175 void ReadImage( typename I::Pointer& image, const std::string& filename )
177 typename itk::ImageFileReader< I >::Pointer reader =
178 itk::ImageFileReader< I >::New( );
179 reader->SetFileName( filename );
182 typename itk::MinimumMaximumImageCalculator< I >::Pointer minmax =
183 itk::MinimumMaximumImageCalculator< I >::New( );
184 minmax->SetImage( reader->GetOutput( ) );
186 double vmin = double( minmax->GetMinimum( ) );
187 double vmax = double( minmax->GetMaximum( ) );
189 typename itk::ShiftScaleImageFilter< I, I >::Pointer shift =
190 itk::ShiftScaleImageFilter< I, I >::New( );
191 shift->SetInput( reader->GetOutput( ) );
192 shift->SetScale( vmax - vmin );
193 shift->SetShift( vmin / ( vmax - vmin ) );
196 image = shift->GetOutput( );
197 image->DisconnectPipeline( );
200 // -------------------------------------------------------------------------
201 template< class I, class O >
203 const typename I::Pointer& input, typename O::Pointer& output
206 typename itk::SignedMaurerDistanceMapImageFilter< I, O >::Pointer dmap =
207 itk::SignedMaurerDistanceMapImageFilter< I, O >::New( );
208 dmap->SetInput( input );
209 dmap->SetBackgroundValue( ( typename I::PixelType )( 0 ) );
210 dmap->InsideIsPositiveOn( );
211 dmap->SquaredDistanceOn( );
212 dmap->UseImageSpacingOn( );
214 std::time_t start, end;
219 << "Distance map time = "
220 << std::difftime( end, start )
221 << " s." << std::endl;
223 output = dmap->GetOutput( );
224 output->DisconnectPipeline( );