appli/CTBronchi/Skeleton.cxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia
   3 // @email florez-l@javeriana.edu.co
   4 // =========================================================================
   5
   6 #include <sstream>
   7 #include <string>
   8 #include <tclap/CmdLine.h>
   9 #include <itkImage.h>
  10 #include <ivq/ITK/ImageSkeletonWriter.h>
  11 #include <fpa/Filters/Image/ExtractSkeleton.h>
  12 #include "Functions.h"
  13
  14 // -------------------------------------------------------------------------
  15 const unsigned int Dim = 3;
  16 typedef unsigned char             TPixel;
  17 typedef itk::Image< TPixel, Dim > TImage;
  18
  19 // -------------------------------------------------------------------------
  20 int main( int argc, char* argv[] )
  21 {
  22   typedef TCLAP::ValueArg< std::string >   _TStringArg;
  23   typedef TCLAP::ValueArg< double >        _TRealArg;
  24   typedef TCLAP::ValueArg< unsigned long > _TUIntArg;
  25   typedef TCLAP::ValueArg< TPixel >        _TPixelArg;
  26
  27   // Parse input line
  28   _TStringArg in( "i", "input", "Input image", true, "", "file" );
  29   _TStringArg out( "o", "output", "Output skeleton", true, "", "file" );
  30   _TStringArg oute( "e", "end_points", "Output end points", false, "", "file" );
  31   _TStringArg aSeed( "p", "seed", "Input seed", true, "", "\"x y z\"" );
  32   try
  33   {
  34     TCLAP::CmdLine cmd( "Skeleton", ' ', "1.0.0" );
  35     cmd.add( aSeed );
  36     cmd.add( oute );
  37     cmd.add( out );
  38     cmd.add( in );
  39     cmd.parse( argc, argv );
  40   }
  41   catch( TCLAP::ArgException& err )
  42   {
  43     std::cerr
  44       << "===============================" << std::endl
  45       << "Error caught: " << std::endl
  46       << err.error( ) << " " << err.argId( )
  47       << "===============================" << std::endl
  48       << std::endl;
  49     return( 1 );
  50
  51   } // yrt
  52
  53   // Get seed
  54   std::istringstream sSeed( aSeed.getValue( ) );
  55   TImage::PointType pnt;
  56   sSeed >> pnt[ 0 ] >> pnt[ 1 ] >> pnt[ 2 ];
  57
  58   try
  59   {
  60     // Read image
  61     TImage::Pointer input_image;
  62     CTBronchi::ReadImage( input_image, in.getValue( ) );
  63
  64     TImage::IndexType seed;
  65     input_image->TransformPhysicalPointToIndex( pnt, seed );
  66
  67     // Skeleton extraction
  68     typedef fpa::Filters::Image::ExtractSkeleton< TImage > TFilter;
  69     TFilter::Pointer filter = TFilter::New( );
  70     filter->SetInput( input_image );
  71     filter->SeedFromMaximumDistanceOff( );
  72     filter->SetSeed( seed );
  73     filter->GetDistanceMap( )->InsideIsPositiveOn( );
  74     filter->GetDistanceMap( )->SquaredDistanceOff( );
  75     filter->GetDistanceMap( )->UseImageSpacingOn( );
  76     double t = CTBronchi::MeasureTime( filter );
  77     std::cout << "Skeleton executed in " << t << " s" << std::endl;
  78
  79     // Save results
  80     typedef ivq::ITK::ImageSkeletonWriter< TFilter::TSkeleton > TWriter;
  81     TWriter::Pointer skeleton_writer = TWriter::New( );
  82     skeleton_writer->SetInput( filter->GetOutput( ) );
  83     skeleton_writer->SetFileName( out.getValue( ) );
  84     skeleton_writer->Update( );
  85
  86     // Write end points
  87     if( oute.getValue( ) != "" )
  88     {
  89       std::stringstream eStr;
  90       auto pnts = filter->GetEndPoints( );
  91       for( unsigned int i = 0; i < 10; ++i )
  92       {
  93         eStr
  94           << pnts[ i ][ 0 ] << " "
  95           << pnts[ i ][ 1 ] << " "
  96           << pnts[ i ][ 2 ] << " " << i % 3 << std::endl;
  97
  98       } // rof
  99       std::ofstream oStrE( oute.getValue( ).c_str( ) );
 100       oStrE << eStr.str( );
 101       oStrE.close( );
 102
 103     } // fi
 104   }
 105   catch( std::exception& err )
 106   {
 107     std::cerr
 108       << "===============================" << std::endl
 109       << "Error caught: " << std::endl
 110       << err.what( )
 111       << "===============================" << std::endl
 112       << std::endl;
 113     return( 1 );
 114
 115   } // yrt
 116   return( 0 );
 117 }
 118
 119 // eof - $RCSfile$