appli/examples/example_ImageAlgorithmDijkstra_03.cxx

   1 #include <cmath>
   2 #include <ctime>
   3 #include <deque>
   4 #include <iostream>
   5 #include <limits>
   6 #include <map>
   7 #include <string>
   8
   9 #include <itkConstNeighborhoodIterator.h>
  10 #include <itkNeighborhoodIterator.h>
  11 #include <itkDanielssonDistanceMapImageFilter.h>
  12 #include <itkImage.h>
  13 #include <itkImageFileReader.h>
  14 #include <itkImageFileWriter.h>
  15 #include <itkImageToVTKImageFilter.h>
  16
  17 #include <vtkPoints.h>
  18 #include <vtkCellArray.h>
  19 #include <vtkPolyData.h>
  20 #include <vtkSmartPointer.h>
  21
  22 #include <fpa/Base/Functors/InvertCostFunction.h>
  23 #include <fpa/Image/Dijkstra.h>
  24 #include <fpa/Image/RegionGrowWithMultipleThresholds.h>
  25 #include <fpa/VTK/ImageMPR.h>
  26 #include <fpa/VTK/Image3DObserver.h>
  27
  28 // -------------------------------------------------------------------------
  29 const unsigned int Dim = 3;
  30 typedef double                               TPixel;
  31 typedef double                               TScalar;
  32 typedef itk::Image< TPixel, Dim >            TImage;
  33 typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
  34
  35 typedef itk::ImageFileReader< TImage >          TImageReader;
  36 typedef itk::ImageFileWriter< TImage >          TImageWriter;
  37 typedef fpa::Image::Dijkstra< TImage, TScalar > TBaseDijkstra;
  38
  39 typedef fpa::VTK::ImageMPR TMPR;
  40 typedef fpa::VTK::Image3DObserver< TBaseDijkstra, vtkRenderWindow > TDijkstraObs;
  41
  42 // -------------------------------------------------------------------------
  43 class TDijkstra
  44   : public TBaseDijkstra
  45 {
  46 public:
  47   typedef TDijkstra                       Self;
  48   typedef TBaseDijkstra                   Superclass;
  49   typedef itk::SmartPointer< Self >       Pointer;
  50   typedef itk::SmartPointer< const Self > ConstPointer;
  51
  52   typedef Superclass::TCost          TCost;
  53   typedef Superclass::TVertex        TVertex;
  54   typedef Superclass::InputImageType TImage;
  55   typedef std::deque< TVertex >      TVertices;
  56
  57   typedef Superclass::TEndEvent          TEndEvent;
  58   typedef Superclass::TBacktrackingEvent TBacktrackingEvent;
  59
  60 protected:
  61   typedef std::pair< TCost, TVertex >     _TCandidate;
  62   typedef std::multimap< TCost, TVertex > _TCandidates;
  63   typedef Superclass::_TNode _TNode;
  64
  65 public:
  66   itkNewMacro( Self );
  67   itkTypeMacro( TDijkstra, TBaseDijkstra );
  68
  69 protected:
  70   TDijkstra( )
  71     : Superclass( )
  72     { }
  73   virtual ~TDijkstra( )
  74     { }
  75
  76   virtual void _BeforeMainLoop( )
  77     {
  78       const TImage* img = this->GetInput( );
  79
  80       // Correct seeds
  81       TImage::SizeType radius;
  82       radius.Fill( 3 );
  83       itk::ConstNeighborhoodIterator< TImage > it(
  84         radius, img, img->GetRequestedRegion( )
  85         );
  86       for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
  87       {
  88         _TNode seed = this->m_Seeds[ s ];
  89         it.SetLocation( seed.Vertex );
  90
  91         TImage::SizeType size = it.GetSize( );
  92         unsigned int nN = 1;
  93         for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
  94           nN *= size[ d ];
  95         TImage::PixelType max_value = img->GetPixel( seed.Vertex );
  96         for( unsigned int i = 0; i < nN; ++i )
  97         {
  98           if( it.GetPixel( i ) > max_value )
  99           {
 100             seed.Vertex = it.GetIndex( i );
 101             seed.Parent = seed.Vertex;
 102             max_value = it.GetPixel( i );
 103
 104           } // fi
 105
 106         } // rof
 107         this->m_Seeds[ s ] = seed;
 108
 109       } // rof
 110
 111       // End initialization
 112       this->Superclass::_BeforeMainLoop( );
 113       this->m_Candidates.clear( );
 114     }
 115
 116   virtual void _AfterMainLoop( )
 117     {
 118       typedef unsigned char                                _TMark;
 119       typedef itk::Image< _TMark, TImage::ImageDimension > _TMarkImage;
 120
 121       this->Superclass::_AfterMainLoop( );
 122       if( this->m_Candidates.size( ) == 0 )
 123         return;
 124
 125       this->InvokeEvent( TEndEvent( ) );
 126
 127       const TImage* input = this->GetInput( );
 128       TImage::SpacingType spacing = input->GetSpacing( );
 129
 130       // Prepare an object to hold marks
 131       _TMarkImage::Pointer marks = _TMarkImage::New( );
 132       marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
 133       marks->SetRequestedRegion( input->GetRequestedRegion( ) );
 134       marks->SetBufferedRegion( input->GetBufferedRegion( ) );
 135       marks->SetOrigin( input->GetOrigin( ) );
 136       marks->SetSpacing( spacing );
 137       marks->SetDirection( input->GetDirection( ) );
 138       marks->Allocate( );
 139       marks->FillBuffer( _TMark( 0 ) );
 140
 141       // Iterate over the candidates, starting from the candidates that
 142       // are near thin branches
 143       _TCandidates::const_reverse_iterator cIt =
 144         this->m_Candidates.rbegin( );
 145       for( int backId = 0; backId < 100 && cIt != this->m_Candidates.rend( ); ++cIt )
 146       {
 147         // If pixel hasn't been visited, continue
 148         TVertex v = cIt->second;
 149         if( marks->GetPixel( v ) != _TMark( 0 ) )
 150           continue;
 151
 152         // Compute nearest start candidate
 153         TImage::SizeType radius;
 154         radius.Fill( 3 );
 155         itk::ConstNeighborhoodIterator< TImage > iIt(
 156           radius, input, input->GetRequestedRegion( )
 157           );
 158         iIt.SetLocation( v );
 159         unsigned int nN = 1;
 160         for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 161           nN *= iIt.GetSize( )[ d ];
 162         TVertex max_vertex = iIt.GetIndex( 0 );
 163         TImage::PixelType max_value = iIt.GetPixel( 0 );
 164         for( unsigned int i = 1; i < nN; ++i )
 165         {
 166           TImage::PixelType value = iIt.GetPixel( i );
 167           if( value > max_value )
 168           {
 169             max_value = value;
 170             max_vertex = iIt.GetIndex( i );
 171
 172           } // fi
 173
 174         } // rof
 175
 176         if( marks->GetPixel( max_vertex ) != _TMark( 0 ) )
 177           continue;
 178         backId++;
 179         std::cout << "Leaf: " << backId << " " << max_value << " " << max_vertex << std::endl;
 180
 181         bool start = true;
 182         do
 183         {
 184           double dist = double( 1.5 ) * std::sqrt( double( input->GetPixel( max_vertex ) ) );
 185           for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 186             radius[ d ] =
 187               ( unsigned long )( dist / spacing[ d ] ) + 1;
 188           itk::NeighborhoodIterator< _TMarkImage > mIt(
 189             radius, marks, marks->GetRequestedRegion( )
 190             );
 191           mIt.SetLocation( max_vertex );
 192           nN = 1;
 193           for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 194             nN *= mIt.GetSize( )[ d ];
 195           for( unsigned int i = 0; i < nN; ++i )
 196             if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
 197             {
 198               mIt.SetPixel( i, ( start )? 255: 100 );
 199               start = false;
 200             }
 201
 202           /*
 203             TImage::SizeType radius;
 204             mIt.GoToBegin( );
 205             mIt.SetLocation( vIt );
 206
 207             TImage::SizeType size = mIt.GetSize( );
 208             unsigned int nN = 1;
 209             for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 210             nN *= size[ d ];
 211             for( unsigned int i = 0; i < nN; ++i )
 212             if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
 213             mIt.SetPixel( i, ( start )? 255: 100 );
 214
 215             start = false;
 216           */
 217           // Next vertex in current path
 218           this->InvokeEvent( TBacktrackingEvent( max_vertex, backId ) );
 219           max_vertex = this->_Parent( max_vertex );
 220
 221         } while( max_vertex != this->_Parent( max_vertex ) );
 222
 223       } // rof
 224       /*
 225         else
 226         marks->SetPixel( v, _TMark( 1 ) );
 227         } // rof
 228       */
 229
 230       /*
 231         itk::ImageFileWriter< _TMarkImage >::Pointer w =
 232         itk::ImageFileWriter< _TMarkImage >::New( );
 233         w->SetInput( marks );
 234         w->SetFileName( "marks.mhd" );
 235         w->Update( );
 236
 237
 238         this->m_Axes = vtkSmartPointer< vtkPolyData >::New( );
 239         vtkSmartPointer< vtkPoints > points =
 240         vtkSmartPointer< vtkPoints >::New( );
 241         vtkSmartPointer< vtkCellArray > lines =
 242         vtkSmartPointer< vtkCellArray >::New( );
 243
 244         {
 245
 246         backId++;
 247         std::cout << "Leaf: " << backId << " " << cIt->first << " " << vIt << std::endl;
 248         bool start = true;
 249         do
 250         {
 251         double dist = double( input->GetPixel( vIt ) );
 252         TImage::SizeType radius;
 253         for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 254         radius[ d ] =
 255         ( unsigned long )( dist / spacing[ d ] ) + 1;
 256         itk::NeighborhoodIterator< _TMarkImage > mIt(
 257         radius, marks, marks->GetRequestedRegion( )
 258         );
 259         mIt.GoToBegin( );
 260         mIt.SetLocation( vIt );
 261
 262         TImage::SizeType size = mIt.GetSize( );
 263         unsigned int nN = 1;
 264         for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 265         nN *= size[ d ];
 266         for( unsigned int i = 0; i < nN; ++i )
 267         if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
 268         mIt.SetPixel( i, ( start )? 255: 100 );
 269
 270         start = false;
 271
 272         // Next vertex
 273         vIt = this->_Parent( vIt );
 274
 275         } while( vIt != this->_Parent( vIt ) );
 276
 277         // Last vertex
 278         // std::cout << vIt << " " << this->GetSeed( 0 ) << std::endl;
 279         */
 280       /* TODO
 281          TVertices pS;
 282          TVertex parent = this->_Parent( vIt );
 283          while( parent != vIt )
 284          {
 285          pS.push_front( vIt );
 286          vIt = parent;
 287          parent = this->_Parent( vIt );
 288
 289          TImage::PointType pnt;
 290          this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
 291          points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 292          if( points->GetNumberOfPoints( ) > 1 )
 293          {
 294          lines->InsertNextCell( 2 );
 295          lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
 296          lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
 297
 298          } // fi
 299
 300          } // elihw
 301          pS.push_front( vIt );
 302          TImage::PointType pnt;
 303          this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
 304          points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 305          if( points->GetNumberOfPoints( ) > 1 )
 306          {
 307          lines->InsertNextCell( 2 );
 308          lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
 309          lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
 310
 311          } // fi
 312
 313          this->m_Axes->SetPoints( points );
 314          this->m_Axes->SetLines( lines );
 315       */
 316
 317       /*
 318         } // rof
 319       */
 320     }
 321
 322   virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n )
 323     {
 324       TCost nc = this->_Cost( nn.Vertex, n.Vertex );
 325       if( TCost( 0 ) < nc )
 326       {
 327         nn.Cost = n.Cost + ( TCost( 1 ) / nc );
 328         nn.Result = nn.Cost;
 329         return( true );
 330       }
 331       else
 332         return( false );
 333     }
 334
 335   virtual bool _UpdateResult( _TNode& n )
 336     {
 337       bool ret = this->Superclass::_UpdateResult( n );
 338       if( ret )
 339       {
 340         TCost nc = this->_Cost( n.Vertex, n.Parent );
 341         if( TCost( 0 ) < nc )
 342         {
 343           TCost cc = n.Cost / nc;
 344           this->m_Candidates.insert( _TCandidate( cc, n.Vertex ) );
 345           /* TODO: debug code
 346              this->GetOutput( )->SetPixel( n.Vertex, cc );
 347           */
 348         } // fi
 349
 350       } // fi
 351       return( ret );
 352     }
 353
 354 private:
 355   TDijkstra( const Self& other );
 356   Self& operator=( const Self& other );
 357
 358 protected:
 359   _TCandidates m_Candidates;
 360 public:
 361   vtkSmartPointer< vtkPolyData > m_Axes;
 362 };
 363
 364 // -------------------------------------------------------------------------
 365 int main( int argc, char* argv[] )
 366 {
 367   if( argc < 6 )
 368   {
 369     std::cerr
 370       << "Usage: " << argv[ 0 ]
 371       << " input_image x y z output_image"
 372       << " visual_debug"
 373       << std::endl;
 374     return( 1 );
 375
 376   } // fi
 377   std::string input_image_fn = argv[ 1 ];
 378   TImage::IndexType seed_idx;
 379   seed_idx[ 0 ] = std::atoi( argv[ 2 ] );
 380   seed_idx[ 1 ] = std::atoi( argv[ 3 ] );
 381   seed_idx[ 2 ] = std::atoi( argv[ 4 ] );
 382   std::string output_image_fn = argv[ 5 ];
 383   bool visual_debug = false;
 384   if( argc > 6 )
 385     visual_debug = ( std::atoi( argv[ 6 ] ) == 1 );
 386
 387   // Read image
 388   TImageReader::Pointer input_image_reader = TImageReader::New( );
 389   input_image_reader->SetFileName( input_image_fn );
 390   try
 391   {
 392     input_image_reader->Update( );
 393   }
 394   catch( itk::ExceptionObject& err )
 395   {
 396     std::cerr << "Error caught: " << err << std::endl;
 397     return( 1 );
 398
 399   } // yrt
 400   TImage::ConstPointer input_image = input_image_reader->GetOutput( );
 401
 402   // Show image
 403   TVTKImage::Pointer vtk_image = TVTKImage::New( );
 404   vtk_image->SetInput( input_image );
 405   vtk_image->Update( );
 406
 407   TMPR view;
 408   view.SetBackground( 0.3, 0.2, 0.8 );
 409   view.SetSize( 800, 800 );
 410   view.SetImage( vtk_image->GetOutput( ) );
 411
 412   // Allow some interaction
 413   view.Start( );
 414
 415   // Extract paths
 416   TDijkstra::Pointer paths = TDijkstra::New( );
 417   paths->AddSeed( seed_idx, TScalar( 0 ) );
 418   paths->SetInput( input_image );
 419   paths->SetNeighborhoodOrder( 1 );
 420
 421   if( visual_debug )
 422   {
 423     // Configure observer
 424     TDijkstraObs::Pointer obs = TDijkstraObs::New( );
 425     obs->SetRenderWindow( view.GetWindow( ) );
 426     paths->AddObserver( itk::AnyEvent( ), obs );
 427     paths->ThrowEventsOn( );
 428   }
 429   else
 430     paths->ThrowEventsOff( );
 431   std::clock_t start = std::clock( );
 432   paths->Update( );
 433   std::clock_t end = std::clock( );
 434   double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
 435   std::cout << "Paths extraction time = " << seconds << std::endl;
 436
 437   view.AddPolyData( paths->m_Axes, 1, 0, 0 );
 438   view.Start( );
 439
 440   TImageWriter::Pointer dijkstra_writer =
 441     TImageWriter::New( );
 442   dijkstra_writer->SetInput( paths->GetOutput( ) );
 443   dijkstra_writer->SetFileName( "dijkstra.mhd" );
 444   dijkstra_writer->Update( );
 445
 446   // Show result
 447   /*
 448     TVTKImage::Pointer output_image_vtk = TVTKImage::New( );
 449     output_image_vtk->SetInput( segmentor->GetOutput( ) );
 450     output_image_vtk->Update( );
 451
 452     vtkSmartPointer< vtkImageMarchingCubes > mc =
 453     vtkSmartPointer< vtkImageMarchingCubes >::New( );
 454     mc->SetInputData( output_image_vtk->GetOutput( ) );
 455     mc->SetValue(
 456     0,
 457     double( segmentor->GetInsideValue( ) ) * double( 0.95 )
 458     );
 459     mc->Update( );
 460
 461     // Let some interaction and close program
 462     view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 );
 463     view.Start( );
 464
 465     // Write resulting image
 466     TImageWriter::Pointer output_image_writer = TImageWriter::New( );
 467     output_image_writer->SetInput( segmentor->GetOutput( ) );
 468     output_image_writer->SetFileName( output_image_fn );
 469     try
 470     {
 471     output_image_writer->Update( );
 472     }
 473     catch( itk::ExceptionObject& err )
 474     {
 475     std::cerr << "Error caught: " << err << std::endl;
 476     return( 1 );
 477
 478     } // yrt
 479   */
 480   return( 0 );
 481 }
 482
 483 // eof - $RCSfile$