7 #include <itkRGBAPixel.h>
8 #include <itkImageToVTKImageFilter.h>
10 #include <itkImageFileReader.h>
11 #include <itkMinimumMaximumImageCalculator.h>
12 #include <itkInvertIntensityImageFilter.h>
13 #include <itkDanielssonDistanceMapImageFilter.h>
14 #include <itkImageFileWriter.h>
16 #include <vtkCamera.h>
17 #include <vtkImageActor.h>
18 #include <vtkInteractorStyleImage.h>
19 #include <vtkPointHandleRepresentation3D.h>
20 #include <vtkProperty.h>
21 #include <vtkRenderer.h>
22 #include <vtkRenderWindow.h>
23 #include <vtkRenderWindowInteractor.h>
24 #include <vtkSeedRepresentation.h>
25 #include <vtkSeedWidget.h>
26 #include <vtkSmartPointer.h>
28 #include <fpa/Image/Dijkstra.h>
29 #include <fpa/Base/Functors/InvertCostFunction.h>
30 #include <fpa/VTK/Image2DObserver.h>
32 // -------------------------------------------------------------------------
33 const unsigned int Dim = 2;
34 typedef unsigned char TInputPixel;
35 typedef float TOutputPixel;
37 typedef itk::Image< TInputPixel, Dim > TInputImage;
38 typedef itk::Image< TOutputPixel, Dim > TOutputImage;
39 typedef itk::ImageToVTKImageFilter< TInputImage > TVTKInputImage;
41 // -------------------------------------------------------------------------
42 int main( int argc, char* argv[] )
47 << "Usage: " << argv[ 0 ]
48 << " input_image output_image neighborhood_order stop_at_one_front"
53 std::string input_image_fn = argv[ 1 ];
54 std::string output_image_fn = argv[ 2 ];
55 unsigned int neighborhood_order = std::atoi( argv[ 3 ] );
56 bool stop_at_one_front = ( std::atoi( argv[ 4 ] ) != 0 );
59 itk::ImageFileReader< TInputImage >::Pointer input_image_reader =
60 itk::ImageFileReader< TInputImage >::New( );
61 input_image_reader->SetFileName( input_image_fn );
64 input_image_reader->Update( );
66 catch( itk::ExceptionObject& err )
69 << "Error while reading image from " << input_image_fn << ": "
74 TInputImage::ConstPointer input_image = input_image_reader->GetOutput( );
75 TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( );
76 vtk_input_image->SetInput( input_image );
77 vtk_input_image->Update( );
80 vtkSmartPointer< vtkImageActor > actor =
81 vtkSmartPointer< vtkImageActor >::New( );
82 actor->SetInputData( vtk_input_image->GetOutput( ) );
84 vtkSmartPointer< vtkRenderer > renderer =
85 vtkSmartPointer< vtkRenderer >::New( );
86 renderer->SetBackground( 0.1, 0.2, 0.7 );
87 renderer->AddActor( actor );
88 vtkSmartPointer< vtkRenderWindow > window =
89 vtkSmartPointer< vtkRenderWindow >::New( );
90 window->SetSize( 800, 800 );
91 window->AddRenderer( renderer );
93 // Correct camera due to the loaded image
94 vtkCamera* camera = renderer->GetActiveCamera( );
95 camera->SetViewUp( 0, -1, 0 );
96 camera->SetPosition( 0, 0, -1 );
97 camera->SetFocalPoint( 0, 0, 0 );
100 vtkSmartPointer< vtkInteractorStyleImage > imageStyle =
101 vtkSmartPointer< vtkInteractorStyleImage >::New( );
102 vtkSmartPointer< vtkRenderWindowInteractor > interactor =
103 vtkSmartPointer< vtkRenderWindowInteractor >::New( );
104 interactor->SetInteractorStyle( imageStyle );
105 window->SetInteractor( interactor );
107 // Create the widget and its representation
108 vtkSmartPointer< vtkPointHandleRepresentation3D > handle =
109 vtkSmartPointer< vtkPointHandleRepresentation3D >::New( );
110 handle->GetProperty( )->SetColor( 1, 0, 0 );
111 vtkSmartPointer< vtkSeedRepresentation > rep =
112 vtkSmartPointer< vtkSeedRepresentation >::New( );
113 rep->SetHandleRepresentation( handle );
115 vtkSmartPointer< vtkSeedWidget > widget =
116 vtkSmartPointer< vtkSeedWidget >::New( );
117 widget->SetInteractor( interactor );
118 widget->SetRepresentation( rep );
120 // Let some interaction
121 interactor->Initialize( );
122 renderer->ResetCamera( );
125 interactor->Start( );
127 // Invert input image
128 itk::MinimumMaximumImageCalculator< TInputImage >::Pointer minmax =
129 itk::MinimumMaximumImageCalculator< TInputImage >::New( );
130 minmax->SetImage( input_image );
133 itk::InvertIntensityImageFilter< TInputImage >::Pointer invert =
134 itk::InvertIntensityImageFilter< TInputImage >::New( );
135 invert->SetInput( input_image );
136 invert->SetMaximum( minmax->GetMaximum( ) );
138 itk::DanielssonDistanceMapImageFilter< TInputImage, TOutputImage >::Pointer dmap =
139 itk::DanielssonDistanceMapImageFilter< TInputImage, TOutputImage >::New( );
140 dmap->SetInput( invert->GetOutput( ) );
141 dmap->InputIsBinaryOn( );
142 dmap->SquaredDistanceOn( );
143 dmap->UseImageSpacingOn( );
146 typedef fpa::Base::Functors::InvertCostFunction< TOutputPixel > TFunction;
147 TFunction::Pointer function = TFunction::New( );
149 // Prepare region grow filter
150 typedef fpa::Image::Dijkstra< TOutputImage, TOutputImage > TFilter;
151 TFilter::Pointer filter = TFilter::New( );
152 filter->SetInput( dmap->GetOutput( ) );
153 filter->SetConversionFunction( function );
154 filter->SetNeighborhoodOrder( neighborhood_order );
155 filter->SetStopAtOneFront( stop_at_one_front );
157 // Get user-given seeds
158 for( unsigned int s = 0; s < rep->GetNumberOfSeeds( ); s++ )
161 rep->GetSeedWorldPosition( s, pos );
163 TInputImage::PointType pnt;
164 pnt[ 0 ] = TInputImage::PointType::ValueType( pos[ 0 ] );
165 pnt[ 1 ] = TInputImage::PointType::ValueType( pos[ 1 ] );
167 TInputImage::IndexType idx;
168 if( input_image->TransformPhysicalPointToIndex( pnt, idx ) )
169 filter->AddSeed( idx, 0 );
173 // Prepare graphical debugger
174 typedef fpa::VTK::Image2DObserver< TFilter, vtkRenderWindow > TDebugger;
175 TDebugger::Pointer debugger = TDebugger::New( );
176 debugger->SetRenderWindow( window );
177 debugger->SetRenderPercentage( 0.001 );
178 filter->AddObserver( itk::AnyEvent( ), debugger );
179 filter->ThrowEventsOn( );
184 // Save final total cost map
185 itk::ImageFileWriter< TOutputImage >::Pointer output_image_writer =
186 itk::ImageFileWriter< TOutputImage >::New( );
187 output_image_writer->SetFileName( output_image_fn );
188 output_image_writer->SetInput( filter->GetOutput( ) );
191 output_image_writer->Update( );
193 catch( itk::ExceptionObject& err )
196 << "Error while writing image to " << output_image_fn << ": "
202 // Draw the path between the two seeds
203 if( filter->GetNumberOfSeeds( ) > 1 )
205 TInputImage::IndexType a = filter->GetSeed( 0 );
206 TInputImage::IndexType b =
207 filter->GetSeed( filter->GetNumberOfSeeds( ) - 1 );
209 std::vector< TInputImage::IndexType > path;
210 filter->GetMinimumSpanningTree( )->GetPath( path, a, b );
212 typedef itk::Image< itk::RGBAPixel< unsigned char >, Dim > TColorImage;
213 TColorImage::Pointer path_image = TColorImage::New( );
214 path_image->SetLargestPossibleRegion(
215 input_image->GetLargestPossibleRegion( )
217 path_image->SetRequestedRegion( input_image->GetRequestedRegion( ) );
218 path_image->SetBufferedRegion( input_image->GetBufferedRegion( ) );
219 path_image->SetSpacing( input_image->GetSpacing( ) );
220 path_image->SetOrigin( input_image->GetOrigin( ) );
221 path_image->SetDirection( input_image->GetDirection( ) );
222 path_image->Allocate( );
224 unsigned char background[ 4 ] = { 0, 0, 0, 0 };
225 unsigned char color[ 4 ] = { 255, 0, 0, 255 };
226 path_image->FillBuffer( itk::RGBAPixel< unsigned char >( background ) );
228 for( unsigned int pId = 0; pId < path.size( ); ++pId )
229 path_image->SetPixel(
230 path[ pId ], itk::RGBAPixel< unsigned char >( color )
233 itk::ImageToVTKImageFilter< TColorImage >::Pointer vtk_path_image =
234 itk::ImageToVTKImageFilter< TColorImage >::New( );
235 vtk_path_image->SetInput( path_image );
236 vtk_path_image->Update( );
238 vtkSmartPointer< vtkImageActor > path_image_actor =
239 vtkSmartPointer< vtkImageActor >::New( );
240 path_image_actor->SetInputData( vtk_path_image->GetOutput( ) );
241 renderer->AddActor( path_image_actor );
246 interactor->Start( );