appli/CTBronchi/MoriLabelling.cxx

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