6 #include <itkImageFileReader.h>
7 #include <itkImageFileWriter.h>
8 #include <itkSimpleFilterWatcher.h>
11 #include <cpPlugins/Extensions/Algorithms/MultiScaleGaussianImageFilter.h>
12 #include <cpPlugins/Extensions/Algorithms/GulsunTekMedialness.h>
13 #include <cpPlugins/Extensions/Algorithms/ImageFunctionFilter.h>
15 #include <cpPlugins/Extensions/Algorithms/InertiaTensorFunction.h>
17 #include <itkImageToVTKImageFilter.h>
19 #include <vtkAxesActor.h>
20 #include <vtkMatrix4x4.h>
21 #include <vtkRenderWindow.h>
22 #include <vtkRenderWindowInteractor.h>
23 #include <vtkSmartPointer.h>
24 #include <cpPlugins/Extensions/Visualization/MPRWithDifferentWindows.h>
26 #include <cpPlugins/Extensions/Algorithms/InertiaMedialness.h>
29 // -------------------------------------------------------------------------
31 typedef float TScalar;
32 const unsigned int Dimension = 2;
34 typedef itk::Vector< TScalar, Dimension > TVector;
35 typedef itk::Image< TPixel, Dimension > TImage;
36 typedef itk::Image< TVector, Dimension > TGradient;
37 typedef itk::Image< TScalar, Dimension > TScalarImage;
39 // -------------------------------------------------------------------------
40 int main( int argc, char* argv[] )
45 << "Usage: " << argv[ 0 ]
46 << " input_image output_medialness"
47 << " min_radius max_radius profile_sampling radial_sampling"
53 std::string image_fn = argv[ 1 ];
54 std::string medialness_fn = argv[ 2 ];
55 double min_radius = std::atof( argv[ 3 ] );
56 double max_radius = std::atof( argv[ 4 ] );
57 unsigned int profile_sampling = std::atoi( argv[ 5 ] );
58 unsigned int radial_sampling = std::atoi( argv[ 6 ] );
61 itk::ImageFileReader< TImage >::Pointer image_reader =
62 itk::ImageFileReader< TImage >::New( );
63 image_reader->SetFileName( image_fn );
66 image_reader->Update( );
68 catch( itk::ExceptionObject& err )
70 std::cerr << err << std::endl;
77 typedef cpPlugins::Extensions::Algorithms::
78 InertiaTensorFunction< TScalar, Dimension >
80 TInertiaFunction::Pointer inertia = TInertiaFunction::New( );
82 itk::ImageRegionConstIteratorWithIndex< TImage > it(
83 image_reader->GetOutput( ),
84 image_reader->GetOutput( )->GetRequestedRegion( )
86 for( it.GoToBegin( ); !it.IsAtEnd( ); ++it )
88 TInertiaFunction::TPoint pnt;
89 image_reader->GetOutput( )->
90 TransformIndexToPhysicalPoint( it.GetIndex( ), pnt );
91 inertia->AddMass( pnt, TScalar( it.Get( ) ) );
95 TInertiaFunction::TMatrix I = inertia->GetInertia( );
96 TInertiaFunction::TVector m = inertia->GetCenterOfGravity( );
98 TInertiaFunction::TVector pv;
99 TInertiaFunction::TVector r;
100 TInertiaFunction::TMatrix pm;
101 inertia->GetEigenAnalysis( pm, pv, r );
103 vtkSmartPointer< vtkMatrix4x4 > matrix =
104 vtkSmartPointer< vtkMatrix4x4 >::New( );
106 matrix->SetElement( 0, 0, pm[ 0 ][ 0 ] * r[ 0 ] );
107 matrix->SetElement( 0, 1, pm[ 0 ][ 1 ] * r[ 1 ] );
108 matrix->SetElement( 1, 0, pm[ 1 ][ 0 ] * r[ 0 ] );
109 matrix->SetElement( 1, 1, pm[ 1 ][ 1 ] * r[ 1 ] );
110 matrix->SetElement( 0, 3, m[ 0 ] );
111 matrix->SetElement( 1, 3, m[ 1 ] );
113 vtkSmartPointer< vtkAxesActor > axes =
114 vtkSmartPointer< vtkAxesActor >::New( );
115 axes->SetUserMatrix( matrix );
117 std::cout << "---------------------------------" << std::endl;
118 std::cout << I << std::endl;
119 std::cout << m << std::endl;
120 std::cout << "---------------------------------" << std::endl;
121 std::cout << pm << std::endl;
122 std::cout << pv << std::endl;
123 std::cout << "---------------------------------" << std::endl;
125 typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
126 TVTKImage::Pointer vtk_input_image = TVTKImage::New( );
127 vtk_input_image->SetInput( image_reader->GetOutput( ) );
128 vtk_input_image->Update( );
130 // Visualization window and interactor
131 vtkSmartPointer< vtkRenderWindowInteractor > interactor =
132 vtkSmartPointer< vtkRenderWindowInteractor >::New( );
133 vtkSmartPointer< vtkRenderWindow > window =
134 vtkSmartPointer< vtkRenderWindow >::New( );
135 window->SetSize( 500, 500 );
136 window->SetInteractor( interactor );
139 cpPlugins::Extensions::Visualization::MPRWithDifferentWindows
140 view( NULL, NULL, window, NULL );
141 view.SetImage( vtk_input_image->GetOutput( ) );
142 view.GetRenderer( 2 )->AddActor( axes );
143 view.ResetCamera( 2 );
147 interactor->Initialize( );
148 interactor->Start( );
154 typedef cpPlugins::Extensions::Algorithms::
155 MultiScaleGaussianImageFilter< TImage, TGradient > TGradientFilter;
157 TGradientFilter::Pointer gradientFilter = TGradientFilter::New( );
158 gradientFilter->SetFilterToGradient( );
159 gradientFilter->SetInput( image_reader->GetOutput( ) );
160 for( int i = 7; i < argc; ++i )
161 gradientFilter->AddScale( std::atof( argv[ i ] ) );
164 itk::SimpleFilterWatcher gradientFilter_watcher( gradientFilter, "Gradient" );
165 gradientFilter->Update( );
168 // Prepare medialness function
169 typedef cpPlugins::Extensions::Algorithms::
170 InertiaMedialness< TImage, TScalar > TMedialness;
171 TMedialness::Pointer medialness = TMedialness::New( );
172 medialness->SetInputImage( image_reader->GetOutput( ) );
173 medialness->SetMaxRadius( max_radius );
175 medialness->SetMinRadius( min_radius );
176 medialness->SetProfileSampling( profile_sampling );
177 medialness->SetRadialSampling( radial_sampling );
179 TGradient::IndexType idx;
183 std::cout << "Medialness = " << medialness->EvaluateAtIndex( idx ) << std::endl;
187 // Prepare output filter
188 typedef cpPlugins::Extensions::Algorithms::
189 ImageFunctionFilter< TImage, TScalarImage, TMedialness >
191 TFunctionApplier::Pointer medialnessFilter = TFunctionApplier::New( );
192 medialnessFilter->SetInput( image_reader->GetOutput( ) );
193 medialnessFilter->SetFunction( medialness );
196 itk::SimpleFilterWatcher medialnessFilter_watcher(
197 medialnessFilter, "Medialness"
199 medialnessFilter->Update( );
202 itk::ImageFileWriter< TScalarImage >::Pointer medialness_writer =
203 itk::ImageFileWriter< TScalarImage >::New( );
204 medialness_writer->SetInput( medialnessFilter->GetOutput( ) );
205 medialness_writer->SetFileName( medialness_fn );
208 medialness_writer->Update( );
210 catch( itk::ExceptionObject& err )
212 std::cerr << err << std::endl;