7 #include <itkImageToVTKImageFilter.h>
9 #include <itkImageFileReader.h>
10 #include <itkMinimumMaximumImageCalculator.h>
11 #include <itkInvertIntensityImageFilter.h>
12 #include <itkDanielssonDistanceMapImageFilter.h>
13 #include <itkImageFileWriter.h>
15 #include <vtkCamera.h>
16 #include <vtkImageActor.h>
17 #include <vtkInteractorStyleImage.h>
18 #include <vtkPointHandleRepresentation3D.h>
19 #include <vtkProperty.h>
20 #include <vtkRenderer.h>
21 #include <vtkRenderWindow.h>
22 #include <vtkRenderWindowInteractor.h>
23 #include <vtkSeedRepresentation.h>
24 #include <vtkSeedWidget.h>
25 #include <vtkSmartPointer.h>
27 #include <fpa/Image/Dijkstra.h>
28 #include <fpa/Base/Functors/InvertCostFunction.h>
29 #include <fpa/VTK/Image2DObserver.h>
31 // -------------------------------------------------------------------------
32 const unsigned int Dim = 2;
33 typedef unsigned char TInputPixel;
34 typedef float TOutputPixel;
36 typedef itk::Image< TInputPixel, Dim > TInputImage;
37 typedef itk::Image< TOutputPixel, Dim > TOutputImage;
38 typedef itk::ImageToVTKImageFilter< TInputImage > TVTKInputImage;
40 // -------------------------------------------------------------------------
41 int main( int argc, char* argv[] )
46 << "Usage: " << argv[ 0 ]
47 << " input_image output_image neighborhood_order stop_at_one_front"
52 std::string input_image_fn = argv[ 1 ];
53 std::string output_image_fn = argv[ 2 ];
54 unsigned int neighborhood_order = std::atoi( argv[ 3 ] );
55 bool stop_at_one_front = ( std::atoi( argv[ 4 ] ) != 0 );
58 itk::ImageFileReader< TInputImage >::Pointer input_image_reader =
59 itk::ImageFileReader< TInputImage >::New( );
60 input_image_reader->SetFileName( input_image_fn );
63 input_image_reader->Update( );
65 catch( itk::ExceptionObject& err )
68 << "Error while reading image from " << input_image_fn << ": "
73 TInputImage::ConstPointer input_image = input_image_reader->GetOutput( );
74 TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( );
75 vtk_input_image->SetInput( input_image );
76 vtk_input_image->Update( );
79 vtkSmartPointer< vtkImageActor > actor =
80 vtkSmartPointer< vtkImageActor >::New( );
81 actor->SetInputData( vtk_input_image->GetOutput( ) );
83 vtkSmartPointer< vtkRenderer > renderer =
84 vtkSmartPointer< vtkRenderer >::New( );
85 renderer->SetBackground( 0.1, 0.2, 0.7 );
86 renderer->AddActor( actor );
87 vtkSmartPointer< vtkRenderWindow > window =
88 vtkSmartPointer< vtkRenderWindow >::New( );
89 window->SetSize( 800, 800 );
90 window->AddRenderer( renderer );
92 // Correct camera due to the loaded image
93 vtkCamera* camera = renderer->GetActiveCamera( );
94 camera->SetViewUp( 0, -1, 0 );
95 camera->SetPosition( 0, 0, -1 );
96 camera->SetFocalPoint( 0, 0, 0 );
99 vtkSmartPointer< vtkInteractorStyleImage > imageStyle =
100 vtkSmartPointer< vtkInteractorStyleImage >::New( );
101 vtkSmartPointer< vtkRenderWindowInteractor > interactor =
102 vtkSmartPointer< vtkRenderWindowInteractor >::New( );
103 interactor->SetInteractorStyle( imageStyle );
104 window->SetInteractor( interactor );
106 // Create the widget and its representation
107 vtkSmartPointer< vtkPointHandleRepresentation3D > handle =
108 vtkSmartPointer< vtkPointHandleRepresentation3D >::New( );
109 handle->GetProperty( )->SetColor( 1, 0, 0 );
110 vtkSmartPointer< vtkSeedRepresentation > rep =
111 vtkSmartPointer< vtkSeedRepresentation >::New( );
112 rep->SetHandleRepresentation( handle );
114 vtkSmartPointer< vtkSeedWidget > widget =
115 vtkSmartPointer< vtkSeedWidget >::New( );
116 widget->SetInteractor( interactor );
117 widget->SetRepresentation( rep );
119 // Let some interaction
120 interactor->Initialize( );
121 renderer->ResetCamera( );
124 interactor->Start( );
126 // Invert input image
127 itk::MinimumMaximumImageCalculator< TInputImage >::Pointer minmax =
128 itk::MinimumMaximumImageCalculator< TInputImage >::New( );
129 minmax->SetImage( input_image );
132 itk::InvertIntensityImageFilter< TInputImage >::Pointer invert =
133 itk::InvertIntensityImageFilter< TInputImage >::New( );
134 invert->SetInput( input_image );
135 invert->SetMaximum( minmax->GetMaximum( ) );
137 itk::DanielssonDistanceMapImageFilter< TInputImage, TOutputImage >::Pointer dmap =
138 itk::DanielssonDistanceMapImageFilter< TInputImage, TOutputImage >::New( );
139 dmap->SetInput( invert->GetOutput( ) );
140 dmap->InputIsBinaryOn( );
141 dmap->SquaredDistanceOn( );
142 dmap->UseImageSpacingOn( );
145 typedef fpa::Base::Functors::InvertCostFunction< TOutputPixel > TFunction;
146 TFunction::Pointer function = TFunction::New( );
148 // Prepare region grow filter
149 typedef fpa::Image::Dijkstra< TOutputImage, TOutputImage > TFilter;
150 TFilter::Pointer filter = TFilter::New( );
151 filter->SetInput( dmap->GetOutput( ) );
152 filter->SetConversionFunction( function );
153 filter->SetNeighborhoodOrder( neighborhood_order );
154 filter->SetStopAtOneFront( stop_at_one_front );
156 // Get user-given seeds
157 for( unsigned int s = 0; s < rep->GetNumberOfSeeds( ); s++ )
160 rep->GetSeedWorldPosition( s, pos );
162 TInputImage::PointType pnt;
163 pnt[ 0 ] = TInputImage::PointType::ValueType( pos[ 0 ] );
164 pnt[ 1 ] = TInputImage::PointType::ValueType( pos[ 1 ] );
166 TInputImage::IndexType idx;
167 if( input_image->TransformPhysicalPointToIndex( pnt, idx ) )
168 filter->AddSeed( idx, 0 );
172 // Prepare graphical debugger
173 typedef fpa::VTK::Image2DObserver< TFilter, vtkRenderWindow > TDebugger;
174 TDebugger::Pointer debugger = TDebugger::New( );
175 debugger->SetRenderWindow( window );
176 debugger->SetRenderPercentage( 0.001 );
177 filter->AddObserver( itk::AnyEvent( ), debugger );
178 filter->ThrowEventsOn( );
183 // Save final total cost map
184 itk::ImageFileWriter< TOutputImage >::Pointer output_image_writer =
185 itk::ImageFileWriter< TOutputImage >::New( );
186 output_image_writer->SetFileName( output_image_fn );
187 output_image_writer->SetInput( filter->GetOutput( ) );
190 output_image_writer->Update( );
192 catch( itk::ExceptionObject& err )
195 << "Error while writing image to " << output_image_fn << ": "