examples/BronchiiInitialSegmentationWithMori.cxx

   1 #include <climits>
   2 #include <cstdlib>
   3 #include <fstream>
   4 #include <iostream>
   5 #include <string>
   6
   7 #include <itkImage.h>
   8 #include <itkImageFileReader.h>
   9 #include <itkImageFileWriter.h>
  10
  11 #include <fpa/Image/MoriFilter.h>
  12
  13 // -------------------------------------------------------------------------
  14 static const unsigned int Dim = 3;
  15 typedef short TPixel;
  16 typedef itk::Image< TPixel, Dim > TImage;
  17
  18 // -------------------------------------------------------------------------
  19 std::string GetFullPath( const std::string& filename )
  20 {
  21   /* On windows:
  22      #include <windows.h>
  23      TCHAR full_path[MAX_PATH];
  24      GetFullPathName(_T("foo.dat"), MAX_PATH, full_path, NULL);
  25   */
  26
  27   char* buffer = realpath( filename.c_str( ), NULL );
  28   std::string path = buffer;
  29   std::free( buffer );
  30   return( path );
  31 }
  32
  33 // -------------------------------------------------------------------------
  34 std::string GetFullPathToDirectory( const std::string& filename )
  35 {
  36   std::string path = GetFullPath( filename );
  37   std::size_t found = path.find_last_of( "/\\" );
  38   return( path.substr( 0, found + 1 ) );
  39 }
  40
  41 // -------------------------------------------------------------------------
  42 int main( int argc, char* argv[] )
  43 {
  44   // Command configuration
  45   if( argc < 3 )
  46   {
  47     std::cerr
  48       << "Usage: " << argv[ 0 ] << " input_image output_image"
  49       << std::endl;
  50     return( 1 );
  51
  52   } // fi
  53   std::string input_image_filename = GetFullPath( argv[ 1 ] );
  54   std::string output_image_filename = argv[ 2 ];
  55   std::string output_auxiliary_image_filename = output_image_filename + "_aux.mhd";
  56
  57   // Try to guess initial seed
  58   std::string seed_filename = GetFullPathToDirectory( argv[ 1 ] ) + "seed.txt";
  59   std::ifstream seed_str( seed_filename.c_str( ) );
  60   if( !seed_str )
  61   {
  62     std::cerr
  63       << "No \"seed.txt\" file in the same input image directory."
  64       << std::endl;
  65     return( 1 );
  66
  67   } // fi
  68   TImage::IndexType seed;
  69   seed_str >> seed[ 0 ] >> seed[ 1 ] >> seed[ 2 ];
  70   seed_str.close( );
  71
  72   // Read image
  73   typedef itk::ImageFileReader< TImage > TReader;
  74   TReader::Pointer reader = TReader::New( );
  75   reader->SetFileName( input_image_filename );
  76
  77   // Mori's algorithms
  78   typedef fpa::Image::MoriFilter< TImage, TImage > TFilter;
  79   TFilter::Pointer filter = TFilter::New( );
  80   filter->SetInput( reader->GetOutput( ) );
  81   filter->SetSeed( seed );
  82   filter->SetThresholdRange( -200, 100, 1 );
  83   filter->SetInsideValue( 1 );
  84   filter->SetOutsideValue( 0 );
  85
  86   // Write results
  87   typedef itk::ImageFileWriter< TImage > TWriter;
  88   TWriter::Pointer writer = TWriter::New( );
  89   writer->SetInput( filter->GetOutput( ) );
  90   writer->SetFileName( output_image_filename );
  91
  92   // Execute pipeline
  93   try
  94   {
  95     writer->Update( );
  96   }
  97   catch( std::exception& err )
  98   {
  99     std::cerr << "Error caught: " << err.what( ) << std::endl;
 100     return( 1 );
 101
 102   } // yrt
 103
 104   // Write auxiliary image
 105   typedef itk::ImageFileWriter< TFilter::TAuxImage > TAuxWriter;
 106   TAuxWriter::Pointer aux_writer = TAuxWriter::New( );
 107   aux_writer->SetInput( filter->GetAuxiliaryImage( ) );
 108   aux_writer->SetFileName( output_auxiliary_image_filename );
 109
 110   // Execute pipeline
 111   try
 112   {
 113     aux_writer->Update( );
 114   }
 115   catch( std::exception& err )
 116   {
 117     std::cerr << "Error caught: " << err.what( ) << std::endl;
 118     return( 1 );
 119
 120   } // yrt
 121
 122   // Write result signal
 123   std::string output_signal_filename = output_image_filename + ".csv";
 124   std::ofstream output_signal_str( output_signal_filename.c_str( ) );
 125   TFilter::TCurve curve = filter->GetCurve( );
 126   unsigned long max_count = 0;
 127   for( TFilter::TCurve::value_type d : curve )
 128   {
 129     output_signal_str << d.first << " " << d.second << std::endl;
 130     if( max_count < d.second )
 131       max_count = d.second;
 132
 133   } // rof
 134   output_signal_str.close( );
 135
 136   // Write gnuplot script
 137   std::string output_gnuplot_filename = output_image_filename + ".gnuplot";
 138   std::ofstream output_gnuplot_str( output_gnuplot_filename.c_str( ) );
 139   unsigned int thr_pos = filter->GetOptimumThreshold( );
 140   int thr = curve[ thr_pos ].first;
 141   output_gnuplot_str
 142     << "set term png" << std::endl
 143     << "set output \"" << output_image_filename << ".png\"" << std::endl
 144     << "set arrow 1 from " << thr
 145     << ",0 to " << thr << "," << max_count
 146     << std::endl
 147     << "show arrow 1" << std::endl
 148     << "plot \"" << output_signal_filename
 149     << "\" using 1:2 with linespoints title \"Evolution ("
 150     << thr << "," << thr_pos << ") " << "\""
 151     << std::endl;
 152   output_gnuplot_str.close( );
 153
 154   return( 0 );
 155 }
 156
 157 // eof - $RCSfile$