appli/CTBronchi/MoriLabelling.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 "MoriLabelling.h"
  11 #include "Functions.h"
  12
  13 // -------------------------------------------------------------------------
  14 const unsigned int Dim = 3;
  15 typedef short                                  TPixel;
  16 typedef unsigned char                          TLabel;
  17 typedef itk::NumericTraits< TPixel >::RealType TScalar;
  18 typedef itk::Image< TPixel, Dim >              TImage;
  19 typedef itk::Image< TLabel, Dim >              TLabels;
  20 typedef itk::Image< TScalar, Dim >             TScalarImage;
  21
  22 // -------------------------------------------------------------------------
  23 int main( int argc, char* argv[] )
  24 {
  25   typedef TCLAP::ValueArg< std::string > _TStringArg;
  26   typedef TCLAP::ValueArg< double > _TRealArg;
  27   typedef TCLAP::ValueArg< TPixel > _TPixelArg;
  28
  29   // Parse input line
  30   _TStringArg in( "i", "input", "Input image", true, "", "file" );
  31   _TStringArg labels( "l", "labels", "Input labels", true, "", "file" );
  32   _TStringArg vesselness( "v", "vesselness", "Input vesselness", true, "", "file" );
  33   _TStringArg out( "o", "output", "Output image", true, "", "file" );
  34   _TRealArg vThr( "a", "vesselness_thr", "Vesselness threshold", false, 0.05, "value (0.05)" );
  35   _TPixelArg uThr( "u", "upper_thr", "Upper threshold", false, -650, "value (-650)" );
  36   try
  37   {
  38     TCLAP::CmdLine cmd( "Labelling", ' ', "1.0.0" );
  39     cmd.add( uThr );
  40     cmd.add( vThr );
  41     cmd.add( out );
  42     cmd.add( vesselness );
  43     cmd.add( labels );
  44     cmd.add( in );
  45     cmd.parse( argc, argv );
  46   }
  47   catch( TCLAP::ArgException& err )
  48   {
  49     std::cerr
  50       << "===============================" << std::endl
  51       << "Error caught: " << std::endl
  52       << err.error( ) << " " << err.argId( ) << std::endl
  53       << "===============================" << std::endl
  54       << std::endl;
  55     return( 1 );
  56
  57   } // yrt
  58
  59   try
  60   {
  61     // Read input image
  62     TImage::Pointer input_image;
  63     CTBronchi::ReadImage( input_image, in.getValue( ) );
  64
  65     // Read input labels
  66     TLabels::Pointer input_labels;
  67     CTBronchi::ReadImage( input_labels, labels.getValue( ) );
  68
  69     // Read input vesselness
  70     TScalarImage::Pointer input_vesselness;
  71     CTBronchi::ReadImage( input_vesselness, vesselness.getValue( ) );
  72
  73     // Create labels
  74     typedef CTBronchi::MoriLabelling< TImage, TLabels, TScalarImage > _TLabelling;
  75     _TLabelling::Pointer labelling = _TLabelling::New( );
  76     labelling->SetInput( input_image );
  77     labelling->SetInputLabels( input_labels );
  78     labelling->SetInputVesselness( input_vesselness );
  79     labelling->SetVesselnessThreshold( vThr.getValue( ) );
  80     labelling->SetUpperThreshold( uThr.getValue( ) );
  81     labelling->SetInsideValue( 1 );
  82     labelling->SetOutsideValue( 2 );
  83     labelling->SetFillValue( 2 );
  84     double t = CTBronchi::MeasureTime( labelling );
  85     std::cout << "Labelling executed in " << t << " s" << std::endl;
  86
  87     // Write result
  88     CTBronchi::WriteImage( labelling->GetOutput( ), out.getValue( ) );
  89   }
  90   catch( std::exception& err )
  91   {
  92     std::cerr
  93       << "===============================" << std::endl
  94       << "Error caught: " << std::endl
  95       << err.what( ) << std::endl
  96       << "===============================" << std::endl
  97       << std::endl;
  98     return( 1 );
  99
 100   } // yrt
 101   return( 0 );
 102 }
 103
 104 /* TODO
 105    #include <chrono>
 106    #include <map>
 107
 108    #include <itkImage.h>
 109    #include <itkImageFileReader.h>
 110    #include <itkImageFileWriter.h>
 111    #include <itkRegionOfInterestImageFilter.h>
 112    #include <itkMinimumMaximumImageCalculator.h>
 113
 114    #include <fpa/Filters/BaseMarksInterface.h>
 115    #include <fpa/Filters/Image/SeedsFromLabelsInterface.h>
 116    #include <fpa/Filters/Image/DefaultTraits.h>
 117    #include <fpa/Filters/Image/RegionGrow.h>
 118    #include <fpa/Functors/RegionGrow/BinaryThreshold.h>
 119
 120    // -------------------------------------------------------------------------
 121    const unsigned int Dim = 3;
 122    typedef short                                  TPixel;
 123    typedef unsigned char                          TLabel;
 124    typedef itk::NumericTraits< TPixel >::RealType TScalar;
 125    typedef itk::Image< TPixel, Dim >              TImage;
 126    typedef itk::Image< TLabel, Dim >              TLabels;
 127    typedef itk::Image< TScalar, Dim >             TScalarImage;
 128
 129    class MoriLabellingTraits
 130    : public fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel >
 131    {
 132    public:
 133    typedef fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel > Superclass;
 134    typedef typename Superclass::TInternalTraits  TInternalTraits;
 135    typedef typename Superclass::TMarksImage      TMarksImage;
 136    typedef typename Superclass::TFilterInterface TFilterInterface;
 137
 138    typedef fpa::Filters::BaseMarksInterface< TInternalTraits >  TMarksInterface;
 139    typedef fpa::Filters::Image::SeedsFromLabelsInterface< TInternalTraits > TSeedsInterface;
 140    };
 141
 142    class MoriLabelling
 143    : public fpa::Filters::Image::RegionGrow< TImage, TLabels, TLabel, MoriLabellingTraits >
 144    {
 145    public:
 146    typedef MoriLabellingTraits TTraits;
 147    fpaTraitsMacro( TTraits );
 148
 149    typedef fpa::Filters::Image::RegionGrow< TImage, TLabels, TMark, TTraits > Superclass;
 150    typedef MoriLabelling                   Self;
 151    typedef itk::SmartPointer< Self >       Pointer;
 152    typedef itk::SmartPointer< const Self > ConstPointer;
 153    typedef fpa::Functors::RegionGrow::BinaryThreshold< TInputValue > TFunctor;
 154
 155    public:
 156    itkNewMacro( Self );
 157    itkTypeMacro( MoriLabelling, fpa::Filters::Image::RegionGrow );
 158
 159    itkGetConstMacro( LastThreshold, TInputValue );
 160    itkSetMacro( LastThreshold, TInputValue );
 161
 162    itkGetConstMacro( MinVertex, TVertex );
 163    itkGetConstMacro( MaxVertex, TVertex );
 164
 165    ivqITKInputMacro( InputLabels, TLabels );
 166    ivqITKInputMacro( InputVesselness, TScalarImage );
 167
 168    public:
 169    TInputValue GetUpperThreshold( ) const
 170    {
 171    return( this->m_Functor->GetUpperThreshold( ) );
 172    }
 173    void SetUpperThreshold( TInputValue t )
 174    {
 175    this->m_Functor->SetUpperThreshold( t );
 176    }
 177
 178    protected:
 179    MoriLabelling( )
 180    : Superclass( ),
 181    m_LastThreshold( TInputValue( 0 ) ),
 182    m_VesselnessThr( TScalar( 0.05 ) )
 183    {
 184    ivqITKInputConfigureMacro( InputLabels, TLabels );
 185    ivqITKInputConfigureMacro( InputVesselness, TScalarImage );
 186    this->m_Functor = TFunctor::New( );
 187    this->SetPredicate( this->m_Functor );
 188    }
 189
 190    virtual ~MoriLabelling( )
 191    {
 192    }
 193
 194    virtual const itk::DataObject* _GetReferenceInput( ) const override
 195    {
 196    return( this->GetInputLabels( ) );
 197    }
 198    virtual void _BeforeGenerateData( ) override
 199    {
 200    this->Superclass::_BeforeGenerateData( );
 201    this->m_FirstVertex = true;
 202
 203    typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
 204    _TMinMax::Pointer minMax = _TMinMax::New( );
 205    minMax->SetImage( this->GetInputVesselness( ) );
 206    minMax->Compute( );
 207    this->m_MaxVesselness = ( 1.0  - this->m_VesselnessThr ) * minMax->GetMaximum( );
 208    }
 209
 210    virtual void _PostComputeOutputValue( TNode& n ) override
 211    {
 212    this->Superclass::_PostComputeOutputValue( n );
 213    if( n.Value == this->GetInsideValue( ) )
 214    {
 215    const TImage* input = this->GetInput( );
 216    const TLabels* labels = this->GetInputLabels( );
 217    const TScalarImage* vesselness = this->GetInputVesselness( );
 218    double x = input->GetPixel( n.Vertex );
 219    if( labels->GetPixel( n.Vertex ) == 0 )
 220    {
 221    if( vesselness->GetPixel( n.Vertex ) > this->m_MaxVesselness )
 222    n.Value = this->GetInsideValue( );
 223    else
 224    n.Value = 0;
 225
 226    } // fi
 227
 228    if( !( this->m_FirstVertex ) )
 229    {
 230    for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
 231    {
 232    if( n.Vertex[ d ] < this->m_MinVertex[ d ] )
 233    this->m_MinVertex[ d ] = n.Vertex[ d ];
 234    if( this->m_MaxVertex[ d ] < n.Vertex[ d ] )
 235    this->m_MaxVertex[ d ] = n.Vertex[ d ];
 236
 237    } // rof
 238    }
 239    else
 240    {
 241    this->m_MinVertex = n.Vertex;
 242    this->m_MaxVertex = n.Vertex;
 243    this->m_FirstVertex = false;
 244
 245    } // fi
 246
 247    } // fi
 248    }
 249
 250    private:
 251    // Purposely not implemented.
 252    MoriLabelling( const Self& other );
 253    Self& operator=( const Self& other );
 254
 255    protected:
 256    TFunctor::Pointer m_Functor;
 257    TInputValue m_LastThreshold;
 258    bool m_FirstVertex;
 259    TVertex m_MinVertex;
 260    TVertex m_MaxVertex;
 261
 262    TScalar m_MaxVesselness;
 263    TScalar m_VesselnessThr;
 264    };
 265
 266    // -------------------------------------------------------------------------
 267    double MeasureTime( itk::ProcessObject* f )
 268    {
 269    std::chrono::time_point< std::chrono::high_resolution_clock > s, e;
 270    std::chrono::duration< double > t;
 271    s = std::chrono::high_resolution_clock::now( );
 272    f->Update( );
 273    e = std::chrono::high_resolution_clock::now( );
 274    t = e - s;
 275    return( t.count( ) );
 276    }
 277
 278    // -------------------------------------------------------------------------
 279    template< class _TImagePtr >
 280    void ReadImage( _TImagePtr& image, const std::string& fname )
 281    {
 282    typedef typename _TImagePtr::ObjectType _TImage;
 283    typedef itk::ImageFileReader< _TImage > _TReader;
 284    typename _TReader::Pointer reader = _TReader::New( );
 285    reader->SetFileName( fname );
 286    double t = MeasureTime( reader );
 287    std::cout << "Read " << fname << " in " << t << " s" << std::endl;
 288    image = reader->GetOutput( );
 289    image->DisconnectPipeline( );
 290    }
 291
 292    // -------------------------------------------------------------------------
 293    template< class _TImage >
 294    void WriteImage( const _TImage* image, const std::string& fname )
 295    {
 296    typedef itk::ImageFileWriter< _TImage > _TWriter;
 297    typename _TWriter::Pointer writer = _TWriter::New( );
 298    writer->SetFileName( fname );
 299    writer->SetInput( image );
 300    double t = MeasureTime( writer );
 301    std::cout << "Wrote " << fname << " in " << t << " s" << std::endl;
 302    }
 303
 304    // -------------------------------------------------------------------------
 305    template< class _TImagePtr, class _TRegion >
 306    void ROI(
 307    _TImagePtr& output,
 308    const typename _TImagePtr::ObjectType* input,
 309    const _TRegion& roi
 310    )
 311    {
 312    typedef typename _TImagePtr::ObjectType _TImage;
 313    typedef itk::RegionOfInterestImageFilter< _TImage, _TImage > _TFilter;
 314
 315    typename _TFilter::Pointer filter = _TFilter::New( );
 316    filter->SetInput( input );
 317    filter->SetRegionOfInterest( roi );
 318    double t = MeasureTime( filter );
 319    std::cout << "ROI computed in " << t << " s" << std::endl;
 320    output = filter->GetOutput( );
 321    output->DisconnectPipeline( );
 322    }
 323
 324    // -------------------------------------------------------------------------
 325    bool ParseArgs(
 326    std::map< std::string, std::string >& args, int argc, char* argv[]
 327    )
 328    {
 329    std::set< std::string > mandatory;
 330    mandatory.insert( "in" );
 331    mandatory.insert( "out" );
 332    mandatory.insert( "labels" );
 333    mandatory.insert( "vesselness" );
 334
 335    args[ "upper_threshold" ] = "-650";
 336
 337    int i = 1;
 338    while( i < argc )
 339    {
 340    std::string cmd = argv[ i ] + 1;
 341    args[ cmd ] = argv[ i + 1 ];
 342    i += 2;
 343
 344    } // elihw
 345
 346    bool complete = true;
 347    for( std::string t: mandatory )
 348    complete &= ( args.find( t ) != args.end( ) );
 349
 350    if( !complete )
 351    {
 352    std::cerr
 353    << "Usage: " << argv[ 0 ] << std::endl
 354    << "\t-in filename" << std::endl
 355    << "\t-out filename" << std::endl
 356    << "\t-labels filename" << std::endl
 357    << "\t-vesselness filename" << std::endl
 358    << "\t[-roi] filename" << std::endl
 359    << "\t[-upper_threshold value]" << std::endl;
 360    return( false );
 361
 362    } // fi
 363    return( true );
 364    }
 365
 366    // -------------------------------------------------------------------------
 367    int main( int argc, char* argv[] )
 368    {
 369    std::map< std::string, std::string > args;
 370    try
 371    {
 372    if( ParseArgs( args, argc, argv ) )
 373    {
 374    // Read input image
 375    TImage::Pointer input_image;
 376    ReadImage( input_image, args[ "in" ] );
 377
 378    // Read labels image
 379    TLabels::Pointer input_labels;
 380    ReadImage( input_labels, args[ "labels" ] );
 381
 382    // Read vesselness image
 383    TScalarImage::Pointer input_vesselness;
 384    ReadImage( input_vesselness, args[ "vesselness" ] );
 385
 386    // Mori labelling
 387    MoriLabelling::Pointer labelling = MoriLabelling::New( );
 388    labelling->SetInput( input_image );
 389    labelling->SetInputLabels( input_labels );
 390    labelling->SetInputVesselness( input_vesselness );
 391    labelling->SetOutsideValue( 2 );
 392    labelling->SetFillValue( 2 );
 393    labelling->SetInsideValue( 1 );
 394    labelling->SetUpperThreshold(
 395    TPixel( std::atof( args[ "upper_threshold" ].c_str( ) ) )
 396    );
 397    double t = MeasureTime( labelling );
 398    std::cout << "Labelling executed in " << t << " s" << std::endl;
 399
 400    std::map< std::string, std::string >::const_iterator i =
 401    args.find( "roi" );
 402    if( i != args.end( ) )
 403    {
 404    TImage::IndexType minV = labelling->GetMinVertex( );
 405    TImage::IndexType maxV = labelling->GetMaxVertex( );
 406    TImage::SizeType roiS;
 407    for( unsigned d = 0; d < TImage::ImageDimension; ++d )
 408    roiS[ d ] = maxV[ d ] - minV[ d ] + 1;
 409    TImage::RegionType roi;
 410    roi.SetIndex( minV );
 411    roi.SetSize( roiS );
 412
 413    // ROI input image
 414    TImage::Pointer input_image_roi;
 415    ROI( input_image_roi, input_image.GetPointer( ), roi );
 416    WriteImage( input_image_roi.GetPointer( ), i->second );
 417
 418    // ROI output image
 419    TLabels::Pointer output_roi;
 420    ROI( output_roi, labelling->GetOutput( ), roi );
 421    WriteImage( output_roi.GetPointer( ), args[ "out" ] );
 422    }
 423    else
 424    WriteImage( labelling->GetOutput( ), args[ "out" ] );
 425    }
 426    else
 427    return( 1 );
 428    }
 429    catch( std::exception& err )
 430    {
 431    std::cerr
 432    << "===============================" << std::endl
 433    << "Error caught: " << std::endl
 434    << err.what( )
 435    << "===============================" << std::endl
 436    << std::endl;
 437    return( 1 );
 438
 439    } // yrt
 440    return( 0 );
 441    }
 442 */
 443
 444 // eof - $RCSfile$