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