appli/examples/example_2DGulsunTekMedialness.cxx

   1 #include <cstdlib>
   2 #include <iostream>
   3 #include <string>
   4
   5 #include <itkImage.h>
   6 #include <itkImageFileReader.h>
   7 #include <itkImageFileWriter.h>
   8 #include <itkSimpleFilterWatcher.h>
   9 #include <itkVector.h>
  10
  11 #include <cpPlugins/Extensions/Algorithms/MultiScaleGaussianImageFilter.h>
  12 #include <cpPlugins/Extensions/Algorithms/GulsunTekMedialness.h>
  13 #include <cpPlugins/Extensions/Algorithms/ImageFunctionFilter.h>
  14
  15 #include <cpPlugins/Extensions/Algorithms/InertiaTensorFunction.h>
  16
  17 #include <itkImageToVTKImageFilter.h>
  18
  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>
  25
  26 #include <cpPlugins/Extensions/Algorithms/InertiaMedialness.h>
  27
  28
  29 // -------------------------------------------------------------------------
  30 typedef short        TPixel;
  31 typedef float        TScalar;
  32 const   unsigned int Dimension = 2;
  33
  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;
  38
  39 // -------------------------------------------------------------------------
  40 int main( int argc, char* argv[] )
  41 {
  42   if( argc < 8 )
  43   {
  44     std::cout
  45       << "Usage: " << argv[ 0 ]
  46       << " input_image output_medialness"
  47       << " min_radius max_radius profile_sampling radial_sampling"
  48       << " s0 s1 ..."
  49       << std::endl;
  50     return( 1 );
  51
  52   } // fi
  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 ] );
  59
  60   // Read image
  61   itk::ImageFileReader< TImage >::Pointer image_reader =
  62     itk::ImageFileReader< TImage >::New( );
  63   image_reader->SetFileName( image_fn );
  64   try
  65   {
  66     image_reader->Update( );
  67   }
  68   catch( itk::ExceptionObject& err )
  69   {
  70     std::cerr << err << std::endl;
  71     return( 1 );
  72
  73   } // yrt
  74
  75   // Inertia tensor
  76   /* TODO
  77   typedef cpPlugins::Extensions::Algorithms::
  78     InertiaTensorFunction< TScalar, Dimension >
  79     TInertiaFunction;
  80   TInertiaFunction::Pointer inertia = TInertiaFunction::New( );
  81
  82   itk::ImageRegionConstIteratorWithIndex< TImage > it(
  83     image_reader->GetOutput( ),
  84     image_reader->GetOutput( )->GetRequestedRegion( )
  85     );
  86   for( it.GoToBegin( ); !it.IsAtEnd( ); ++it )
  87   {
  88     TInertiaFunction::TPoint pnt;
  89     image_reader->GetOutput( )->
  90       TransformIndexToPhysicalPoint( it.GetIndex( ), pnt );
  91     inertia->AddMass( pnt, TScalar( it.Get( ) ) );
  92
  93   } // rof
  94
  95   TInertiaFunction::TMatrix I = inertia->GetInertia( );
  96   TInertiaFunction::TVector m = inertia->GetCenterOfGravity( );
  97
  98   TInertiaFunction::TVector pv;
  99   TInertiaFunction::TVector r;
 100   TInertiaFunction::TMatrix pm;
 101   inertia->GetEigenAnalysis( pm, pv, r );
 102
 103   vtkSmartPointer< vtkMatrix4x4 > matrix =
 104     vtkSmartPointer< vtkMatrix4x4 >::New( );
 105   matrix->Identity( );
 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 ] );
 112
 113   vtkSmartPointer< vtkAxesActor > axes =
 114     vtkSmartPointer< vtkAxesActor >::New( );
 115   axes->SetUserMatrix( matrix );
 116
 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;
 124
 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( );
 129
 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 );
 137
 138   // Viewer
 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 );
 144   view.Render( 2 );
 145
 146   // Start
 147   interactor->Initialize( );
 148   interactor->Start( );
 149
 150   */
 151
 152   /*
 153   // Prepare filter
 154   typedef cpPlugins::Extensions::Algorithms::
 155     MultiScaleGaussianImageFilter< TImage, TGradient > TGradientFilter;
 156
 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 ] ) );
 162
 163   // Execute filter
 164   itk::SimpleFilterWatcher gradientFilter_watcher( gradientFilter, "Gradient" );
 165   gradientFilter->Update( );
 166   */
 167
 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 );
 174   /*
 175     medialness->SetMinRadius( min_radius );
 176     medialness->SetProfileSampling( profile_sampling );
 177     medialness->SetRadialSampling( radial_sampling );
 178
 179     TGradient::IndexType idx;
 180     idx[ 0 ] = 288;
 181     idx[ 1 ] = 690;
 182
 183     std::cout << "Medialness = " << medialness->EvaluateAtIndex( idx ) << std::endl;
 184     // return( 0 );
 185     */
 186
 187   // Prepare output filter
 188   typedef cpPlugins::Extensions::Algorithms::
 189     ImageFunctionFilter< TImage, TScalarImage, TMedialness >
 190     TFunctionApplier;
 191   TFunctionApplier::Pointer medialnessFilter = TFunctionApplier::New( );
 192   medialnessFilter->SetInput( image_reader->GetOutput( ) );
 193   medialnessFilter->SetFunction( medialness );
 194
 195   // Execute filter
 196   itk::SimpleFilterWatcher medialnessFilter_watcher(
 197     medialnessFilter, "Medialness"
 198     );
 199   medialnessFilter->Update( );
 200
 201   // Write image
 202   itk::ImageFileWriter< TScalarImage >::Pointer medialness_writer =
 203     itk::ImageFileWriter< TScalarImage >::New( );
 204   medialness_writer->SetInput( medialnessFilter->GetOutput( ) );
 205   medialness_writer->SetFileName( medialness_fn );
 206   try
 207   {
 208     medialness_writer->Update( );
 209   }
 210   catch( itk::ExceptionObject& err )
 211   {
 212     std::cerr << err << std::endl;
 213     return( 1 );
 214
 215   } // yrt
 216
 217   return( 0 );
 218 }
 219
 220 // eof - $RCSfile$