appli/examples/example_Image_RegionGrow_GaussianModelEstimation.cxx

   1 #include <cstdlib>
   2 #include <iostream>
   3 #include <limits>
   4 #include <string>
   5
   6 #include <itkImage.h>
   7 #include <itkImageToVTKImageFilter.h>
   8 #include <itkRGBPixel.h>
   9
  10 #include <itkImageFileReader.h>
  11
  12 #include <vtkCamera.h>
  13 #include <vtkImageActor.h>
  14 #include <vtkInteractorStyleImage.h>
  15 #include <vtkPointHandleRepresentation3D.h>
  16 #include <vtkProperty.h>
  17 #include <vtkRenderer.h>
  18 #include <vtkRenderWindow.h>
  19 #include <vtkRenderWindowInteractor.h>
  20 #include <vtkSeedRepresentation.h>
  21 #include <vtkSeedWidget.h>
  22 #include <vtkSmartPointer.h>
  23
  24 #include <cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h>
  25 #include <cpPlugins/Extensions/Algorithms/RGBToYPbPrFunction.h>
  26
  27 #include <fpa/Image/RegionGrow.h>
  28 #include <fpa/Base/Functors/TautologyFunction.h>
  29 #include <fpa/VTK/Image2DObserver.h>
  30
  31 // -------------------------------------------------------------------------
  32 const unsigned int Dim = 2;
  33 typedef unsigned char TPixel;
  34 typedef float         TScalar;
  35
  36 typedef itk::Image< itk::RGBPixel< TPixel >, Dim > TColorImage;
  37 typedef itk::Image< TPixel, Dim >                  TImage;
  38 typedef itk::ImageToVTKImageFilter< TColorImage >  TVTKImage;
  39
  40 // -------------------------------------------------------------------------
  41 template< class I, class S >
  42 class GaussianFunction
  43   : public itk::FunctionBase< typename I::PixelType, bool >
  44 {
  45 public:
  46   // Type-related and pointers
  47   typedef GaussianFunction                                 Self;
  48   typedef itk::FunctionBase< typename I::PixelType, bool > Superclass;
  49   typedef itk::SmartPointer< Self >                        Pointer;
  50   typedef itk::SmartPointer< const Self >                  ConstPointer;
  51
  52   // Superclass' types
  53   typedef cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, 3 > TEstimator;
  54   typedef cpPlugins::Extensions::Algorithms::RGBToYPbPrFunction< S > TYPbPrFunction;
  55
  56 public:
  57   itkNewMacro( Self );
  58   itkTypeMacro( GaussianFunction, itkFunctionBase );
  59
  60   itkGetConstMacro( ModelSupport, unsigned long );
  61   itkSetMacro( ModelSupport, unsigned long );
  62
  63 public:
  64   virtual bool Evaluate( const typename I::PixelType& rgb ) const
  65     {
  66       /* TODO
  67       if( !this->m_Estimating )
  68       {
  69         this->m_Estimator->AddSample( this->m_YPbPrFunction( rgb ) );
  70         if( this->m_Estimator->GetNumberOfSamples( ) == this->m_ModelSupport )
  71         {
  72           this->m_Estimating = true;
  73           this->m_Estimator->UpdateModel( );
  74           std::cout << this->m_Estimator->GetMinimumProbability( ) << std::endl;
  75           std::cout << this->m_Estimator->GetMaximumProbability( ) << std::endl;
  76
  77         } // fi
  78         return( true );
  79       }
  80       else
  81       {
  82         S p = this->m_Estimator->Probability( this->m_YPbPrFunction( rgb ) );
  83         return( p > this->m_Estimator->GetMinimumProbability( ) );
  84
  85       } // fi
  86       */
  87       return( true );
  88     }
  89
  90 protected:
  91   GaussianFunction( )
  92     : Superclass( ),
  93       m_ModelSupport( 10 )
  94     {
  95       this->m_Estimator = TEstimator::New( );
  96       this->m_Estimator->Clear( );
  97       this->m_Estimating = false;
  98     }
  99   virtual ~GaussianFunction( )
 100     { }
 101
 102 private:
 103   // Purposely not implemented
 104   GaussianFunction( const Self& );
 105   void operator=( const Self& );
 106
 107 protected:
 108   typename TEstimator::Pointer m_Estimator;
 109   TYPbPrFunction m_YPbPrFunction;
 110
 111   unsigned long m_ModelSupport;
 112   mutable bool m_Estimating;
 113 };
 114
 115 // -------------------------------------------------------------------------
 116 int main( int argc, char* argv[] )
 117 {
 118   if( argc < 4 )
 119   {
 120     std::cerr
 121       << "Usage: " << argv[ 0 ]
 122       << " input_image neighborhood_order support"
 123       << std::endl;
 124     return( 1 );
 125
 126   } // fi
 127   std::string input_image_fn = argv[ 1 ];
 128   unsigned int neighborhood_order = std::atoi( argv[ 2 ] );
 129   unsigned long support = std::atoi( argv[ 3 ] );
 130
 131   // Read image
 132   itk::ImageFileReader< TColorImage >::Pointer input_image_reader =
 133     itk::ImageFileReader< TColorImage >::New( );
 134   input_image_reader->SetFileName( input_image_fn );
 135   try
 136   {
 137     input_image_reader->Update( );
 138   }
 139   catch( itk::ExceptionObject& err )
 140   {
 141     std::cerr
 142       << "Error while reading image from " << input_image_fn << ": "
 143       << err << std::endl;
 144     return( 1 );
 145
 146   } // yrt
 147   TColorImage::ConstPointer input_image = input_image_reader->GetOutput( );
 148   TVTKImage::Pointer vtk_input_image = TVTKImage::New( );
 149   vtk_input_image->SetInput( input_image );
 150   vtk_input_image->Update( );
 151
 152   // VTK visualization
 153   vtkSmartPointer< vtkImageActor > actor =
 154     vtkSmartPointer< vtkImageActor >::New( );
 155   actor->SetInputData( vtk_input_image->GetOutput( ) );
 156
 157   vtkSmartPointer< vtkRenderer > renderer =
 158     vtkSmartPointer< vtkRenderer >::New( );
 159   renderer->SetBackground( 0.1, 0.2, 0.7 );
 160   renderer->AddActor( actor );
 161   vtkSmartPointer< vtkRenderWindow > window =
 162     vtkSmartPointer< vtkRenderWindow >::New( );
 163   window->SetSize( 800, 800 );
 164   window->AddRenderer( renderer );
 165
 166   // Correct camera due to the loaded image
 167   vtkCamera* camera = renderer->GetActiveCamera( );
 168   camera->SetViewUp( 0, -1, 0 );
 169   camera->SetPosition( 0, 0, -1 );
 170   camera->SetFocalPoint( 0, 0, 0 );
 171
 172   // VTK interaction
 173   vtkSmartPointer< vtkInteractorStyleImage > imageStyle =
 174     vtkSmartPointer< vtkInteractorStyleImage >::New( );
 175   vtkSmartPointer< vtkRenderWindowInteractor > interactor =
 176     vtkSmartPointer< vtkRenderWindowInteractor >::New( );
 177   interactor->SetInteractorStyle( imageStyle );
 178   window->SetInteractor( interactor );
 179
 180   // Create the widget and its representation
 181   vtkSmartPointer< vtkPointHandleRepresentation3D > handle =
 182     vtkSmartPointer< vtkPointHandleRepresentation3D >::New( );
 183   handle->GetProperty( )->SetColor( 1, 0, 0 );
 184   vtkSmartPointer< vtkSeedRepresentation > rep =
 185     vtkSmartPointer< vtkSeedRepresentation >::New( );
 186   rep->SetHandleRepresentation( handle );
 187
 188   vtkSmartPointer< vtkSeedWidget > widget =
 189     vtkSmartPointer< vtkSeedWidget >::New( );
 190   widget->SetInteractor( interactor );
 191   widget->SetRepresentation( rep );
 192
 193   // Let some interaction
 194   interactor->Initialize( );
 195   renderer->ResetCamera( );
 196   window->Render( );
 197   widget->On( );
 198   interactor->Start( );
 199
 200   // Prepare region grow function
 201   typedef GaussianFunction< TColorImage, TScalar > TFunction;
 202   TFunction::Pointer function = TFunction::New( );
 203   function->SetModelSupport( support );
 204
 205   // Prepare region grow filter
 206   typedef fpa::Image::RegionGrow< TColorImage, TImage > TFilter;
 207   TFilter::Pointer filter = TFilter::New( );
 208   filter->SetInput( input_image );
 209   filter->SetMembershipFunction( function );
 210   filter->SetNeighborhoodOrder( neighborhood_order );
 211   filter->SetOutsideValue( TPixel( 0 ) );
 212   filter->SetInsideValue( std::numeric_limits< TPixel >::max( ) );
 213   filter->StopAtOneFrontOff( );
 214
 215   // Get user-given seeds
 216   for( unsigned int s = 0; s < rep->GetNumberOfSeeds( ); s++ )
 217   {
 218     double pos[ 3 ];
 219     rep->GetSeedWorldPosition( s, pos );
 220
 221     TImage::PointType pnt;
 222     pnt[ 0 ] = TImage::PointType::ValueType( pos[ 0 ] );
 223     pnt[ 1 ] = TImage::PointType::ValueType( pos[ 1 ] );
 224
 225     TImage::IndexType idx;
 226     if( input_image->TransformPhysicalPointToIndex( pnt, idx ) )
 227       filter->AddSeed( idx, 0 );
 228
 229   } // rof
 230
 231   // Prepare graphical debugger
 232   typedef fpa::VTK::Image2DObserver< TFilter, vtkRenderWindow > TDebugger;
 233   TDebugger::Pointer debugger = TDebugger::New( );
 234   debugger->SetRenderWindow( window );
 235   debugger->SetRenderPercentage( 0.001 );
 236   filter->AddObserver( itk::AnyEvent( ), debugger );
 237   filter->ThrowEventsOn( );
 238
 239   // Go!
 240   filter->Update( );
 241
 242   return( 0 );
 243 }
 244
 245 // eof - $RCSfile$