appli/examples/example_ImageAlgorithmDijkstra_03.cxx

   1 #include <ctime>
   2 #include <iostream>
   3 #include <limits>
   4 #include <map>
   5 #include <string>
   6 #include <deque>
   7
   8 #include <itkConstNeighborhoodIterator.h>
   9 #include <itkDanielssonDistanceMapImageFilter.h>
  10 #include <itkImage.h>
  11 #include <itkImageFileReader.h>
  12 #include <itkImageFileWriter.h>
  13 #include <itkImageToVTKImageFilter.h>
  14
  15 #include <vtkPoints.h>
  16 #include <vtkCellArray.h>
  17 #include <vtkPolyData.h>
  18 #include <vtkSmartPointer.h>
  19
  20 #include <fpa/Base/Functors/InvertCostFunction.h>
  21 #include <fpa/Image/Dijkstra.h>
  22 #include <fpa/Image/RegionGrowWithMultipleThresholds.h>
  23 #include <fpa/VTK/ImageMPR.h>
  24 #include <fpa/VTK/Image3DObserver.h>
  25
  26 // -------------------------------------------------------------------------
  27 const unsigned int Dim = 3;
  28 typedef double                               TPixel;
  29 typedef double                               TScalar;
  30 typedef itk::Image< TPixel, Dim >            TImage;
  31 typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
  32
  33 typedef itk::ImageFileReader< TImage >          TImageReader;
  34 typedef itk::ImageFileWriter< TImage >          TImageWriter;
  35 typedef fpa::Image::Dijkstra< TImage, TScalar > TBaseDijkstra;
  36
  37 typedef fpa::VTK::ImageMPR TMPR;
  38 typedef fpa::VTK::Image3DObserver< TBaseDijkstra, vtkRenderWindow > TDijkstraObs;
  39
  40 // -------------------------------------------------------------------------
  41 class TDijkstra
  42   : public TBaseDijkstra
  43 {
  44 public:
  45   typedef TDijkstra                       Self;
  46   typedef TBaseDijkstra                   Superclass;
  47   typedef itk::SmartPointer< Self >       Pointer;
  48   typedef itk::SmartPointer< const Self > ConstPointer;
  49
  50   typedef Superclass::TCost          TCost;
  51   typedef Superclass::TVertex        TVertex;
  52   typedef Superclass::InputImageType TImage;
  53   typedef std::deque< TVertex >      TVertices;
  54
  55 protected:
  56   typedef std::pair< TCost, TVertex >     _TCandidate;
  57   typedef std::multimap< TCost, TVertex > _TCandidates;
  58   typedef Superclass::_TNode _TNode;
  59
  60 public:
  61   itkNewMacro( Self );
  62   itkTypeMacro( TDijkstra, TBaseDijkstra );
  63
  64 protected:
  65   TDijkstra( )
  66     : Superclass( )
  67     { }
  68   virtual ~TDijkstra( )
  69     { }
  70
  71   virtual void _BeforeMainLoop( )
  72     {
  73       const TImage* img = this->GetInput( );
  74
  75       // Correct seeds
  76       TImage::SizeType radius;
  77       radius.Fill( 3 );
  78       itk::ConstNeighborhoodIterator< TImage > it(
  79         radius, img, img->GetRequestedRegion( )
  80         );
  81       for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
  82       {
  83         _TNode seed = this->m_Seeds[ s ];
  84         it.SetLocation( seed.Vertex );
  85
  86         TImage::SizeType size = it.GetSize( );
  87         unsigned int nN = 1;
  88         for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
  89           nN *= size[ d ];
  90         TImage::PixelType max_value = img->GetPixel( seed.Vertex );
  91         for( unsigned int i = 0; i < nN; ++i )
  92         {
  93           if( it.GetPixel( i ) > max_value )
  94           {
  95             seed.Vertex = it.GetIndex( i );
  96             seed.Parent = seed.Vertex;
  97             max_value = it.GetPixel( i );
  98
  99           } // fi
 100
 101         } // rof
 102         this->m_Seeds[ s ] = seed;
 103
 104       } // rof
 105
 106       // End initialization
 107       this->Superclass::_BeforeMainLoop( );
 108       this->m_Candidates.clear( );
 109     }
 110
 111   virtual void _AfterMainLoop( )
 112     {
 113       this->Superclass::_AfterMainLoop( );
 114       if( this->m_Candidates.size( ) == 0 )
 115         return;
 116
 117       this->m_Axes = vtkSmartPointer< vtkPolyData >::New( );
 118       vtkSmartPointer< vtkPoints > points =
 119         vtkSmartPointer< vtkPoints >::New( );
 120       vtkSmartPointer< vtkCellArray > lines =
 121         vtkSmartPointer< vtkCellArray >::New( );
 122
 123       _TCandidates::const_reverse_iterator cIt =
 124         this->m_Candidates.rbegin( );
 125
 126       TVertices pS;
 127       TVertex vIt = cIt->second;
 128       TVertex parent = this->_Parent( vIt );
 129       while( parent != vIt )
 130       {
 131         pS.push_front( vIt );
 132         vIt = parent;
 133         parent = this->_Parent( vIt );
 134
 135         TImage::PointType pnt;
 136         this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
 137         points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 138         if( points->GetNumberOfPoints( ) > 1 )
 139         {
 140           lines->InsertNextCell( 2 );
 141           lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
 142           lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
 143
 144         } // fi
 145
 146       } // elihw
 147       pS.push_front( vIt );
 148       TImage::PointType pnt;
 149       this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
 150       points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 151         if( points->GetNumberOfPoints( ) > 1 )
 152         {
 153           lines->InsertNextCell( 2 );
 154           lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
 155           lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
 156
 157         } // fi
 158
 159         this->m_Axes->SetPoints( points );
 160         this->m_Axes->SetLines( lines );
 161     }
 162
 163   virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n )
 164     {
 165       TCost nc = this->_Cost( nn.Vertex, n.Vertex );
 166       if( TCost( 0 ) < nc )
 167       {
 168         nn.Cost = n.Cost + ( TCost( 1 ) / nc );
 169         nn.Result = nn.Cost;
 170         return( true );
 171       }
 172       else
 173         return( false );
 174     }
 175
 176   virtual bool _UpdateResult( _TNode& n )
 177     {
 178       bool ret = this->Superclass::_UpdateResult( n );
 179       if( ret )
 180       {
 181         TCost nc = this->_Cost( n.Vertex, n.Parent );
 182         if( TCost( 0 ) < nc )
 183         {
 184           n.Result = n.Cost / ( nc * nc );
 185           this->GetOutput( )->SetPixel( n.Vertex, n.Result );
 186           this->m_Candidates.insert( _TCandidate( n.Result, n.Vertex ) );
 187
 188         } // fi
 189
 190       } // fi
 191       return( ret );
 192     }
 193
 194 private:
 195   TDijkstra( const Self& other );
 196   Self& operator=( const Self& other );
 197
 198 protected:
 199   _TCandidates m_Candidates;
 200 public:
 201   vtkSmartPointer< vtkPolyData > m_Axes;
 202 };
 203
 204 // -------------------------------------------------------------------------
 205 int main( int argc, char* argv[] )
 206 {
 207   if( argc < 6 )
 208   {
 209     std::cerr
 210       << "Usage: " << argv[ 0 ]
 211       << " input_image x y z output_image"
 212       << " visual_debug"
 213       << std::endl;
 214     return( 1 );
 215
 216   } // fi
 217   std::string input_image_fn = argv[ 1 ];
 218   TImage::IndexType seed_idx;
 219   seed_idx[ 0 ] = std::atoi( argv[ 2 ] );
 220   seed_idx[ 1 ] = std::atoi( argv[ 3 ] );
 221   seed_idx[ 2 ] = std::atoi( argv[ 4 ] );
 222   std::string output_image_fn = argv[ 5 ];
 223   bool visual_debug = false;
 224   if( argc > 6 )
 225     visual_debug = ( std::atoi( argv[ 6 ] ) == 1 );
 226
 227   // Read image
 228   TImageReader::Pointer input_image_reader = TImageReader::New( );
 229   input_image_reader->SetFileName( input_image_fn );
 230   try
 231   {
 232     input_image_reader->Update( );
 233   }
 234   catch( itk::ExceptionObject& err )
 235   {
 236     std::cerr << "Error caught: " << err << std::endl;
 237     return( 1 );
 238
 239   } // yrt
 240   TImage::ConstPointer input_image = input_image_reader->GetOutput( );
 241
 242   // Show image
 243   TVTKImage::Pointer vtk_image = TVTKImage::New( );
 244   vtk_image->SetInput( input_image );
 245   vtk_image->Update( );
 246
 247   TMPR view;
 248   view.SetBackground( 0.3, 0.2, 0.8 );
 249   view.SetSize( 800, 800 );
 250   view.SetImage( vtk_image->GetOutput( ) );
 251
 252   // Allow some interaction
 253   view.Start( );
 254
 255   // Extract paths
 256   TDijkstra::Pointer paths = TDijkstra::New( );
 257   paths->AddSeed( seed_idx, TScalar( 0 ) );
 258   paths->SetInput( input_image );
 259   paths->SetNeighborhoodOrder( 1 );
 260
 261   if( visual_debug )
 262   {
 263     // Configure observer
 264     TDijkstraObs::Pointer obs = TDijkstraObs::New( );
 265     obs->SetImage( input_image, view.GetWindow( ) );
 266     paths->AddObserver( itk::AnyEvent( ), obs );
 267     paths->ThrowEventsOn( );
 268   }
 269   else
 270     paths->ThrowEventsOff( );
 271   std::clock_t start = std::clock( );
 272   paths->Update( );
 273   std::clock_t end = std::clock( );
 274   double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
 275   std::cout << "Paths extraction time = " << seconds << std::endl;
 276
 277   view.AddPolyData( paths->m_Axes, 1, 0, 0 );
 278   view.Start( );
 279
 280   TImageWriter::Pointer dijkstra_writer =
 281     TImageWriter::New( );
 282   dijkstra_writer->SetInput( paths->GetOutput( ) );
 283   dijkstra_writer->SetFileName( "dijkstra.mhd" );
 284   dijkstra_writer->Update( );
 285
 286   // Show result
 287   /*
 288     TVTKImage::Pointer output_image_vtk = TVTKImage::New( );
 289     output_image_vtk->SetInput( segmentor->GetOutput( ) );
 290     output_image_vtk->Update( );
 291
 292     vtkSmartPointer< vtkImageMarchingCubes > mc =
 293     vtkSmartPointer< vtkImageMarchingCubes >::New( );
 294     mc->SetInputData( output_image_vtk->GetOutput( ) );
 295     mc->SetValue(
 296     0,
 297     double( segmentor->GetInsideValue( ) ) * double( 0.95 )
 298     );
 299     mc->Update( );
 300
 301     // Let some interaction and close program
 302     view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 );
 303     view.Start( );
 304
 305     // Write resulting image
 306     TImageWriter::Pointer output_image_writer = TImageWriter::New( );
 307     output_image_writer->SetInput( segmentor->GetOutput( ) );
 308     output_image_writer->SetFileName( output_image_fn );
 309     try
 310     {
 311     output_image_writer->Update( );
 312     }
 313     catch( itk::ExceptionObject& err )
 314     {
 315     std::cerr << "Error caught: " << err << std::endl;
 316     return( 1 );
 317
 318     } // yrt
 319   */
 320   return( 0 );
 321 }
 322
 323 // eof - $RCSfile$