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 fromt 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 < 100 && 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         TImage::PixelType min_value = iIt.GetPixel( 0 );
 158         TVertex min_vertex = iIt.GetIndex( 0 );
 159         if( !( min_value > TImage::PixelType( 0 ) ) )
 160           min_value = std::numeric_limits< TImage::PixelType >::max( );
 161         for( unsigned int i = 1; i < nN; ++i )
 162         {
 163           TImage::PixelType value = iIt.GetPixel( i );
 164           if( !( value > TImage::PixelType( 0 ) ) )
 165             value = std::numeric_limits< TImage::PixelType >::max( );
 166           if( value < min_value )
 167           {
 168             min_value = value;
 169             min_vertex = iIt.GetIndex( i );
 170
 171           } // fi
 172
 173         } // rof
 174
 175         if( min_value < std::numeric_limits< TImage::PixelType >::max( ) )
 176         {
 177           if( marks->GetPixel( min_vertex ) != _TMark( 0 ) )
 178             continue;
 179           leo++;
 180           std::cout << "Leaf: " << leo << " " << min_value << " " << min_vertex << std::endl;
 181
 182           bool start = true;
 183           do
 184           {
 185             double dist = double( 1.5 ) * std::sqrt( double( input->GetPixel( min_vertex ) ) );
 186             for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 187               radius[ d ] =
 188                 ( unsigned long )( dist / spacing[ d ] ) + 1;
 189             itk::NeighborhoodIterator< _TMarkImage > mIt(
 190               radius, marks, marks->GetRequestedRegion( )
 191               );
 192             mIt.SetLocation( min_vertex );
 193             nN = 1;
 194             for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 195               nN *= mIt.GetSize( )[ d ];
 196             for( unsigned int i = 0; i < nN; ++i )
 197               if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
 198               {
 199                 mIt.SetPixel( i, ( start )? 255: 100 );
 200                 start = false;
 201               }
 202
 203             /*
 204               TImage::SizeType radius;
 205               mIt.GoToBegin( );
 206               mIt.SetLocation( vIt );
 207
 208               TImage::SizeType size = mIt.GetSize( );
 209               unsigned int nN = 1;
 210               for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 211               nN *= size[ d ];
 212               for( unsigned int i = 0; i < nN; ++i )
 213               if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
 214               mIt.SetPixel( i, ( start )? 255: 100 );
 215
 216               start = false;
 217             */
 218             // Next vertex in current path
 219             min_vertex = this->_Parent( min_vertex );
 220
 221           } while( min_vertex != this->_Parent( min_vertex ) );
 222         }
 223         else
 224           marks->SetPixel( v, _TMark( 1 ) );
 225
 226       } // rof
 227
 228       itk::ImageFileWriter< _TMarkImage >::Pointer w =
 229         itk::ImageFileWriter< _TMarkImage >::New( );
 230       w->SetInput( marks );
 231       w->SetFileName( "marks.mhd" );
 232       w->Update( );
 233       std::exit( 1 );
 234
 235       /*
 236
 237         this->m_Axes = vtkSmartPointer< vtkPolyData >::New( );
 238         vtkSmartPointer< vtkPoints > points =
 239         vtkSmartPointer< vtkPoints >::New( );
 240         vtkSmartPointer< vtkCellArray > lines =
 241         vtkSmartPointer< vtkCellArray >::New( );
 242
 243         {
 244
 245         leo++;
 246         std::cout << "Leaf: " << leo << " " << cIt->first << " " << vIt << std::endl;
 247         bool start = true;
 248         do
 249         {
 250         double dist = double( input->GetPixel( vIt ) );
 251         TImage::SizeType radius;
 252         for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 253         radius[ d ] =
 254         ( unsigned long )( dist / spacing[ d ] ) + 1;
 255         itk::NeighborhoodIterator< _TMarkImage > mIt(
 256         radius, marks, marks->GetRequestedRegion( )
 257         );
 258         mIt.GoToBegin( );
 259         mIt.SetLocation( vIt );
 260
 261         TImage::SizeType size = mIt.GetSize( );
 262         unsigned int nN = 1;
 263         for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 264         nN *= size[ d ];
 265         for( unsigned int i = 0; i < nN; ++i )
 266         if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
 267         mIt.SetPixel( i, ( start )? 255: 100 );
 268
 269         start = false;
 270
 271         // Next vertex
 272         vIt = this->_Parent( vIt );
 273
 274         } while( vIt != this->_Parent( vIt ) );
 275
 276         // Last vertex
 277         // std::cout << vIt << " " << this->GetSeed( 0 ) << std::endl;
 278         */
 279       /* TODO
 280          TVertices pS;
 281          TVertex parent = this->_Parent( vIt );
 282          while( parent != vIt )
 283          {
 284          pS.push_front( vIt );
 285          vIt = parent;
 286          parent = this->_Parent( vIt );
 287
 288          TImage::PointType pnt;
 289          this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
 290          points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 291          if( points->GetNumberOfPoints( ) > 1 )
 292          {
 293          lines->InsertNextCell( 2 );
 294          lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
 295          lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
 296
 297          } // fi
 298
 299          } // elihw
 300          pS.push_front( vIt );
 301          TImage::PointType pnt;
 302          this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
 303          points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 304          if( points->GetNumberOfPoints( ) > 1 )
 305          {
 306          lines->InsertNextCell( 2 );
 307          lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
 308          lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
 309
 310          } // fi
 311
 312          this->m_Axes->SetPoints( points );
 313          this->m_Axes->SetLines( lines );
 314       */
 315
 316       /*
 317         } // rof
 318       */
 319     }
 320
 321   virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n )
 322     {
 323       TCost nc = this->_Cost( nn.Vertex, n.Vertex );
 324       if( TCost( 0 ) < nc )
 325       {
 326         nn.Cost = n.Cost + ( TCost( 1 ) / nc );
 327         nn.Result = nn.Cost;
 328         return( true );
 329       }
 330       else
 331         return( false );
 332     }
 333
 334   virtual bool _UpdateResult( _TNode& n )
 335     {
 336       bool ret = this->Superclass::_UpdateResult( n );
 337       if( ret )
 338       {
 339         TCost nc = this->_Cost( n.Vertex, n.Parent );
 340         if( TCost( 0 ) < nc )
 341         {
 342           TCost cc = n.Cost / nc;
 343           this->m_Candidates.insert( _TCandidate( cc, n.Vertex ) );
 344           /* TODO: debug code
 345              this->GetOutput( )->SetPixel( n.Vertex, cc );
 346           */
 347         } // fi
 348
 349       } // fi
 350       return( ret );
 351     }
 352
 353 private:
 354   TDijkstra( const Self& other );
 355   Self& operator=( const Self& other );
 356
 357 protected:
 358   _TCandidates m_Candidates;
 359 public:
 360   vtkSmartPointer< vtkPolyData > m_Axes;
 361 };
 362
 363 // -------------------------------------------------------------------------
 364 int main( int argc, char* argv[] )
 365 {
 366   if( argc < 6 )
 367   {
 368     std::cerr
 369       << "Usage: " << argv[ 0 ]
 370       << " input_image x y z output_image"
 371       << " visual_debug"
 372       << std::endl;
 373     return( 1 );
 374
 375   } // fi
 376   std::string input_image_fn = argv[ 1 ];
 377   TImage::IndexType seed_idx;
 378   seed_idx[ 0 ] = std::atoi( argv[ 2 ] );
 379   seed_idx[ 1 ] = std::atoi( argv[ 3 ] );
 380   seed_idx[ 2 ] = std::atoi( argv[ 4 ] );
 381   std::string output_image_fn = argv[ 5 ];
 382   bool visual_debug = false;
 383   if( argc > 6 )
 384     visual_debug = ( std::atoi( argv[ 6 ] ) == 1 );
 385
 386   // Read image
 387   TImageReader::Pointer input_image_reader = TImageReader::New( );
 388   input_image_reader->SetFileName( input_image_fn );
 389   try
 390   {
 391     input_image_reader->Update( );
 392   }
 393   catch( itk::ExceptionObject& err )
 394   {
 395     std::cerr << "Error caught: " << err << std::endl;
 396     return( 1 );
 397
 398   } // yrt
 399   TImage::ConstPointer input_image = input_image_reader->GetOutput( );
 400
 401   // Show image
 402   TVTKImage::Pointer vtk_image = TVTKImage::New( );
 403   vtk_image->SetInput( input_image );
 404   vtk_image->Update( );
 405
 406   TMPR view;
 407   view.SetBackground( 0.3, 0.2, 0.8 );
 408   view.SetSize( 800, 800 );
 409   view.SetImage( vtk_image->GetOutput( ) );
 410
 411   // Allow some interaction
 412   view.Start( );
 413
 414   // Extract paths
 415   TDijkstra::Pointer paths = TDijkstra::New( );
 416   paths->AddSeed( seed_idx, TScalar( 0 ) );
 417   paths->SetInput( input_image );
 418   paths->SetNeighborhoodOrder( 1 );
 419
 420   if( visual_debug )
 421   {
 422     // Configure observer
 423     TDijkstraObs::Pointer obs = TDijkstraObs::New( );
 424     obs->SetImage( input_image, view.GetWindow( ) );
 425     paths->AddObserver( itk::AnyEvent( ), obs );
 426     paths->ThrowEventsOn( );
 427   }
 428   else
 429     paths->ThrowEventsOff( );
 430   std::clock_t start = std::clock( );
 431   paths->Update( );
 432   std::clock_t end = std::clock( );
 433   double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
 434   std::cout << "Paths extraction time = " << seconds << std::endl;
 435
 436   view.AddPolyData( paths->m_Axes, 1, 0, 0 );
 437   view.Start( );
 438
 439   TImageWriter::Pointer dijkstra_writer =
 440     TImageWriter::New( );
 441   dijkstra_writer->SetInput( paths->GetOutput( ) );
 442   dijkstra_writer->SetFileName( "dijkstra.mhd" );
 443   dijkstra_writer->Update( );
 444
 445   // Show result
 446   /*
 447     TVTKImage::Pointer output_image_vtk = TVTKImage::New( );
 448     output_image_vtk->SetInput( segmentor->GetOutput( ) );
 449     output_image_vtk->Update( );
 450
 451     vtkSmartPointer< vtkImageMarchingCubes > mc =
 452     vtkSmartPointer< vtkImageMarchingCubes >::New( );
 453     mc->SetInputData( output_image_vtk->GetOutput( ) );
 454     mc->SetValue(
 455     0,
 456     double( segmentor->GetInsideValue( ) ) * double( 0.95 )
 457     );
 458     mc->Update( );
 459
 460     // Let some interaction and close program
 461     view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 );
 462     view.Start( );
 463
 464     // Write resulting image
 465     TImageWriter::Pointer output_image_writer = TImageWriter::New( );
 466     output_image_writer->SetInput( segmentor->GetOutput( ) );
 467     output_image_writer->SetFileName( output_image_fn );
 468     try
 469     {
 470     output_image_writer->Update( );
 471     }
 472     catch( itk::ExceptionObject& err )
 473     {
 474     std::cerr << "Error caught: " << err << std::endl;
 475     return( 1 );
 476
 477     } // yrt
 478   */
 479   return( 0 );
 480 }
 481
 482 // eof - $RCSfile$