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 itk::Image< unsigned char, TImage::ImageDimension > _TMarkImage;
 116
 117       this->Superclass::_AfterMainLoop( );
 118       if( this->m_Candidates.size( ) == 0 )
 119         return;
 120
 121       const TImage* input = this->GetInput( );
 122       TImage::SpacingType spacing = input->GetSpacing( );
 123
 124       _TMarkImage::Pointer marks = _TMarkImage::New( );
 125       marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
 126       marks->SetRequestedRegion( input->GetRequestedRegion( ) );
 127       marks->SetBufferedRegion( input->GetBufferedRegion( ) );
 128       marks->SetOrigin( input->GetOrigin( ) );
 129       marks->SetSpacing( spacing );
 130       marks->SetDirection( input->GetDirection( ) );
 131       marks->Allocate( );
 132       marks->FillBuffer( 0 );
 133
 134       this->m_Axes = vtkSmartPointer< vtkPolyData >::New( );
 135       vtkSmartPointer< vtkPoints > points =
 136         vtkSmartPointer< vtkPoints >::New( );
 137       vtkSmartPointer< vtkCellArray > lines =
 138         vtkSmartPointer< vtkCellArray >::New( );
 139
 140       _TCandidates::const_reverse_iterator cIt =
 141         this->m_Candidates.rbegin( );
 142
 143       for( unsigned long leo = 0; leo < 50 && cIt != this->m_Candidates.rend( ); ++cIt )
 144       {
 145         TVertex vIt = cIt->second;
 146         if( marks->GetPixel( vIt ) != 0 )
 147           continue;
 148
 149         leo++;
 150         std::cout << "Leaf: " << leo << " " << cIt->first << " " << vIt << std::endl;
 151         bool start = true;
 152         do
 153         {
 154           double dist = double( input->GetPixel( vIt ) );
 155         TImage::SizeType radius;
 156           for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 157             radius[ d ] =
 158               ( unsigned long )( dist / spacing[ d ] ) + 1;
 159           itk::NeighborhoodIterator< _TMarkImage > mIt(
 160             radius, marks, marks->GetRequestedRegion( )
 161             );
 162           mIt.GoToBegin( );
 163           mIt.SetLocation( vIt );
 164
 165           TImage::SizeType size = mIt.GetSize( );
 166           unsigned int nN = 1;
 167           for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 168             nN *= size[ d ];
 169           for( unsigned int i = 0; i < nN; ++i )
 170             if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
 171               mIt.SetPixel( i, ( start )? 255: 100 );
 172
 173           start = false;
 174
 175           // Next vertex
 176           vIt = this->_Parent( vIt );
 177
 178         } while( vIt != this->_Parent( vIt ) );
 179
 180         // Last vertex
 181         // std::cout << vIt << " " << this->GetSeed( 0 ) << std::endl;
 182
 183         /* TODO
 184            TVertices pS;
 185            TVertex parent = this->_Parent( vIt );
 186            while( parent != vIt )
 187            {
 188            pS.push_front( vIt );
 189            vIt = parent;
 190            parent = this->_Parent( vIt );
 191
 192            TImage::PointType pnt;
 193            this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
 194            points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 195            if( points->GetNumberOfPoints( ) > 1 )
 196            {
 197            lines->InsertNextCell( 2 );
 198            lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
 199            lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
 200
 201            } // fi
 202
 203            } // elihw
 204            pS.push_front( vIt );
 205            TImage::PointType pnt;
 206            this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt );
 207            points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 208            if( points->GetNumberOfPoints( ) > 1 )
 209            {
 210            lines->InsertNextCell( 2 );
 211            lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
 212            lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
 213
 214            } // fi
 215
 216            this->m_Axes->SetPoints( points );
 217            this->m_Axes->SetLines( lines );
 218         */
 219
 220       } // rof
 221
 222       itk::ImageFileWriter< _TMarkImage >::Pointer w =
 223         itk::ImageFileWriter< _TMarkImage >::New( );
 224       w->SetInput( marks );
 225       w->SetFileName( "marks.mhd" );
 226       w->Update( );
 227
 228     }
 229
 230   virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n )
 231     {
 232       TCost nc = this->_Cost( nn.Vertex, n.Vertex );
 233       if( TCost( 0 ) < nc )
 234       {
 235         nn.Cost = n.Cost + ( TCost( 1 ) / nc );
 236         nn.Result = nn.Cost;
 237         return( true );
 238       }
 239       else
 240         return( false );
 241     }
 242
 243   virtual bool _UpdateResult( _TNode& n )
 244     {
 245       bool ret = this->Superclass::_UpdateResult( n );
 246       if( ret )
 247       {
 248         TCost nc = this->_Cost( n.Vertex, n.Parent );
 249         TCost nc2 = nc * nc * nc * nc;
 250         if( TCost( 0 ) < nc2 )
 251         {
 252           n.Result = n.Cost / nc2;
 253           this->GetOutput( )->SetPixel( n.Vertex, n.Result );
 254           this->m_Candidates.insert( _TCandidate( n.Result, n.Vertex ) );
 255
 256         } // fi
 257
 258       } // fi
 259       return( ret );
 260     }
 261
 262 private:
 263   TDijkstra( const Self& other );
 264   Self& operator=( const Self& other );
 265
 266 protected:
 267   _TCandidates m_Candidates;
 268 public:
 269   vtkSmartPointer< vtkPolyData > m_Axes;
 270 };
 271
 272 // -------------------------------------------------------------------------
 273 int main( int argc, char* argv[] )
 274 {
 275   if( argc < 6 )
 276   {
 277     std::cerr
 278       << "Usage: " << argv[ 0 ]
 279       << " input_image x y z output_image"
 280       << " visual_debug"
 281       << std::endl;
 282     return( 1 );
 283
 284   } // fi
 285   std::string input_image_fn = argv[ 1 ];
 286   TImage::IndexType seed_idx;
 287   seed_idx[ 0 ] = std::atoi( argv[ 2 ] );
 288   seed_idx[ 1 ] = std::atoi( argv[ 3 ] );
 289   seed_idx[ 2 ] = std::atoi( argv[ 4 ] );
 290   std::string output_image_fn = argv[ 5 ];
 291   bool visual_debug = false;
 292   if( argc > 6 )
 293     visual_debug = ( std::atoi( argv[ 6 ] ) == 1 );
 294
 295   // Read image
 296   TImageReader::Pointer input_image_reader = TImageReader::New( );
 297   input_image_reader->SetFileName( input_image_fn );
 298   try
 299   {
 300     input_image_reader->Update( );
 301   }
 302   catch( itk::ExceptionObject& err )
 303   {
 304     std::cerr << "Error caught: " << err << std::endl;
 305     return( 1 );
 306
 307   } // yrt
 308   TImage::ConstPointer input_image = input_image_reader->GetOutput( );
 309
 310   // Show image
 311   TVTKImage::Pointer vtk_image = TVTKImage::New( );
 312   vtk_image->SetInput( input_image );
 313   vtk_image->Update( );
 314
 315   TMPR view;
 316   view.SetBackground( 0.3, 0.2, 0.8 );
 317   view.SetSize( 800, 800 );
 318   view.SetImage( vtk_image->GetOutput( ) );
 319
 320   // Allow some interaction
 321   view.Start( );
 322
 323   // Extract paths
 324   TDijkstra::Pointer paths = TDijkstra::New( );
 325   paths->AddSeed( seed_idx, TScalar( 0 ) );
 326   paths->SetInput( input_image );
 327   paths->SetNeighborhoodOrder( 1 );
 328
 329   if( visual_debug )
 330   {
 331     // Configure observer
 332     TDijkstraObs::Pointer obs = TDijkstraObs::New( );
 333     obs->SetImage( input_image, view.GetWindow( ) );
 334     paths->AddObserver( itk::AnyEvent( ), obs );
 335     paths->ThrowEventsOn( );
 336   }
 337   else
 338     paths->ThrowEventsOff( );
 339   std::clock_t start = std::clock( );
 340   paths->Update( );
 341   std::clock_t end = std::clock( );
 342   double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
 343   std::cout << "Paths extraction time = " << seconds << std::endl;
 344
 345   view.AddPolyData( paths->m_Axes, 1, 0, 0 );
 346   view.Start( );
 347
 348   TImageWriter::Pointer dijkstra_writer =
 349     TImageWriter::New( );
 350   dijkstra_writer->SetInput( paths->GetOutput( ) );
 351   dijkstra_writer->SetFileName( "dijkstra.mhd" );
 352   dijkstra_writer->Update( );
 353
 354   // Show result
 355   /*
 356     TVTKImage::Pointer output_image_vtk = TVTKImage::New( );
 357     output_image_vtk->SetInput( segmentor->GetOutput( ) );
 358     output_image_vtk->Update( );
 359
 360     vtkSmartPointer< vtkImageMarchingCubes > mc =
 361     vtkSmartPointer< vtkImageMarchingCubes >::New( );
 362     mc->SetInputData( output_image_vtk->GetOutput( ) );
 363     mc->SetValue(
 364     0,
 365     double( segmentor->GetInsideValue( ) ) * double( 0.95 )
 366     );
 367     mc->Update( );
 368
 369     // Let some interaction and close program
 370     view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 );
 371     view.Start( );
 372
 373     // Write resulting image
 374     TImageWriter::Pointer output_image_writer = TImageWriter::New( );
 375     output_image_writer->SetInput( segmentor->GetOutput( ) );
 376     output_image_writer->SetFileName( output_image_fn );
 377     try
 378     {
 379     output_image_writer->Update( );
 380     }
 381     catch( itk::ExceptionObject& err )
 382     {
 383     std::cerr << "Error caught: " << err << std::endl;
 384     return( 1 );
 385
 386     } // yrt
 387   */
 388   return( 0 );
 389 }
 390
 391 // eof - $RCSfile$