appli/examples/example_Image_Dijkstra_EndPointDetection_WithoutVTK.cxx

   1 #include <cstdlib>
   2 #include <ctime>
   3 #include <iostream>
   4 #include <limits>
   5 #include <string>
   6
   7 #include <itkImage.h>
   8
   9 #include <itkImageFileReader.h>
  10 #include <itkMinimumMaximumImageCalculator.h>
  11 #include <itkShiftScaleImageFilter.h>
  12
  13 #include <itkSignedMaurerDistanceMapImageFilter.h>
  14
  15 #include <itkImageFileWriter.h>
  16
  17 #include <fpa/Image/DijkstraWithEndPointDetection.h>
  18 #include <fpa/Base/Functors/InvertCostFunction.h>
  19
  20 #include <fpa/IO/MinimumSpanningTreeWriter.h>
  21 #include <fpa/IO/UniqueValuesContainerWriter.h>
  22 #include <fpa/IO/MatrixValuesContainerWriter.h>
  23
  24 // -------------------------------------------------------------------------
  25 const unsigned int Dim = 3;
  26 typedef unsigned char TPixel;
  27 typedef float         TScalar;
  28
  29 typedef itk::Image< TPixel, Dim >  TImage;
  30 typedef itk::Image< TScalar, Dim > TScalarImage;
  31
  32 // -------------------------------------------------------------------------
  33 template< class I >
  34 void ReadImage( typename I::Pointer& image, const std::string& filename );
  35
  36 template< class I >
  37 void SaveImage( const I* image, const std::string& filename );
  38
  39 template< class I, class O >
  40 void DistanceMap(
  41   const typename I::Pointer& input, typename O::Pointer& output
  42   );
  43
  44 // -------------------------------------------------------------------------
  45 int main( int argc, char* argv[] )
  46 {
  47   if( argc < 12 )
  48   {
  49     std::cerr
  50       << "Usage: " << argv[ 0 ] << std::endl
  51       << " input_image" << std::endl
  52       << " seed_x seed_y seed_z" << std::endl
  53       << " output_distancemap" << std::endl
  54       << " output_costmap" << std::endl
  55       << " output_labels" << std::endl
  56       << " output_minimum_spanning_tree" << std::endl
  57       << " output_endpoints" << std::endl
  58       << " output_bifurcations" << std::endl
  59       << " output_branches" << std::endl
  60       << std::endl;
  61     return( 1 );
  62
  63   } // fi
  64   std::string input_image_fn = argv[ 1 ];
  65   double seed_x = std::atof( argv[ 2 ] );
  66   double seed_y = std::atof( argv[ 3 ] );
  67   double seed_z = std::atof( argv[ 4 ] );
  68   std::string distancemap_fn = argv[ 5 ];
  69   std::string output_costmap_fn = argv[ 6 ];
  70   std::string output_labels_fn = argv[ 7 ];
  71   std::string mst_output_fn = argv[ 8 ];
  72   std::string endpoints_output_fn = argv[ 9 ];
  73   std::string bifurcations_output_fn = argv[ 10 ];
  74   std::string branches_output_fn = argv[ 11 ];
  75
  76   // Read image
  77   TImage::Pointer input_image;
  78   try
  79   {
  80     ReadImage< TImage >( input_image, input_image_fn );
  81   }
  82   catch( itk::ExceptionObject& err )
  83   {
  84     std::cerr
  85       << "Error caught while reading \""
  86       << input_image_fn << "\": " << err
  87       << std::endl;
  88     return( 1 );
  89
  90   } // yrt
  91
  92   // Get seed
  93   double p[ 3 ];
  94   TImage::PointType pnt;
  95   TImage::IndexType seed;
  96   pnt[ 0 ] = seed_x;
  97   pnt[ 1 ] = seed_y;
  98   pnt[ 2 ] = seed_z;
  99   input_image->TransformPhysicalPointToIndex( pnt, seed );
 100
 101   // Compute squared distance map
 102   TScalarImage::Pointer dmap;
 103   DistanceMap< TImage, TScalarImage >( input_image, dmap );
 104
 105   // Prepare cost conversion function
 106   typedef fpa::Base::Functors::InvertCostFunction< TScalar > TFunction;
 107   TFunction::Pointer function = TFunction::New( );
 108
 109   // Prepare Dijkstra filter
 110   typedef fpa::Image::DijkstraWithEndPointDetection< TScalarImage, TScalarImage > TFilter;
 111   TFilter::Pointer filter = TFilter::New( );
 112   filter->SetInput( dmap );
 113   filter->SetConversionFunction( function );
 114   filter->SetNeighborhoodOrder( 2 );
 115   filter->StopAtOneFrontOff( );
 116   filter->AddSeed( seed, TScalar( 0 ) );
 117
 118   // Go!
 119   std::time_t start, end;
 120   std::time( &start );
 121   filter->Update( );
 122   std::time( &end );
 123   std::cout
 124     << "Extraction time = "
 125     << std::difftime( end, start )
 126     << " s." << std::endl;
 127
 128   // Outputs
 129   const TFilter::TOutputImage* accumulated_costs = filter->GetOutput( );
 130   const TFilter::TLabelImage* labeled_image = filter->GetLabelImage( );
 131   const TFilter::TMinimumSpanningTree* mst = filter->GetMinimumSpanningTree( );
 132   const TFilter::TUniqueVertices* endpoints = filter->GetEndPoints( );
 133   const TFilter::TUniqueVertices* bifurcations = filter->GetBifurcations( );
 134   const TFilter::TBranches* branches = filter->GetBranches( );
 135   unsigned long nBranches = filter->GetNumberOfBranches( );
 136
 137   // Save outputs
 138   SaveImage( dmap.GetPointer( ), distancemap_fn );
 139   SaveImage( accumulated_costs, output_costmap_fn );
 140   SaveImage( labeled_image, output_labels_fn );
 141
 142   fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::Pointer
 143     mst_writer =
 144     fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::New( );
 145   mst_writer->SetInput( mst );
 146   mst_writer->SetFileName( mst_output_fn );
 147   mst_writer->Update( );
 148
 149   fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer
 150     endpoints_writer =
 151     fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( );
 152   endpoints_writer->SetInput( endpoints );
 153   endpoints_writer->SetFileName( endpoints_output_fn );
 154   endpoints_writer->Update( );
 155
 156   fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer
 157     bifurcations_writer =
 158     fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( );
 159   bifurcations_writer->SetInput( bifurcations );
 160   bifurcations_writer->SetFileName( bifurcations_output_fn );
 161   bifurcations_writer->Update( );
 162
 163   fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::Pointer
 164     branches_writer =
 165     fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::New( );
 166   branches_writer->SetInput( branches );
 167   branches_writer->SetFileName( branches_output_fn );
 168   branches_writer->Update( );
 169
 170   return( 0 );
 171 }
 172
 173 // -------------------------------------------------------------------------
 174 template< class I >
 175 void ReadImage( typename I::Pointer& image, const std::string& filename )
 176 {
 177   typename itk::ImageFileReader< I >::Pointer reader =
 178     itk::ImageFileReader< I >::New( );
 179   reader->SetFileName( filename );
 180   reader->Update( );
 181
 182   typename itk::MinimumMaximumImageCalculator< I >::Pointer minmax =
 183     itk::MinimumMaximumImageCalculator< I >::New( );
 184   minmax->SetImage( reader->GetOutput( ) );
 185   minmax->Compute( );
 186   double vmin = double( minmax->GetMinimum( ) );
 187   double vmax = double( minmax->GetMaximum( ) );
 188
 189   typename itk::ShiftScaleImageFilter< I, I >::Pointer shift =
 190     itk::ShiftScaleImageFilter< I, I >::New( );
 191   shift->SetInput( reader->GetOutput( ) );
 192   shift->SetScale( vmax - vmin );
 193   shift->SetShift( vmin / ( vmax - vmin ) );
 194   shift->Update( );
 195
 196   image = shift->GetOutput( );
 197   image->DisconnectPipeline( );
 198 }
 199
 200 // -------------------------------------------------------------------------
 201 template< class I >
 202 void SaveImage( const I* image, const std::string& filename )
 203 {
 204   typename itk::ImageFileWriter< I >::Pointer writer =
 205     itk::ImageFileWriter< I >::New( );
 206   writer->SetInput( image );
 207   writer->SetFileName( filename );
 208   try
 209   {
 210     writer->Update( );
 211   }
 212   catch( itk::ExceptionObject& err )
 213   {
 214     std::cerr
 215       << "Error saving \"" << filename << "\": " << err
 216       << std::endl;
 217
 218   } // yrt
 219 }
 220
 221 // -------------------------------------------------------------------------
 222 template< class I, class O >
 223 void DistanceMap(
 224   const typename I::Pointer& input, typename O::Pointer& output
 225   )
 226 {
 227   typename itk::SignedMaurerDistanceMapImageFilter< I, O >::Pointer dmap =
 228     itk::SignedMaurerDistanceMapImageFilter< I, O >::New( );
 229   dmap->SetInput( input );
 230   dmap->SetBackgroundValue( ( typename I::PixelType )( 0 ) );
 231   dmap->InsideIsPositiveOn( );
 232   dmap->SquaredDistanceOn( );
 233   dmap->UseImageSpacingOn( );
 234
 235   std::time_t start, end;
 236   std::time( &start );
 237   dmap->Update( );
 238   std::time( &end );
 239   std::cout
 240     << "Distance map time = "
 241     << std::difftime( end, start )
 242     << " s." << std::endl;
 243
 244   output = dmap->GetOutput( );
 245   output->DisconnectPipeline( );
 246 }
 247
 248 // eof - $RCSfile$
 249