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