1 // =========================================================================
2 // @author Leonardo Florez Valencia
3 // @email florez-l@javeriana.edu.co
4 // =========================================================================
8 #include <tclap/CmdLine.h>
10 #include "MoriLabelling.h"
11 #include "Functions.h"
13 // -------------------------------------------------------------------------
14 const unsigned int Dim = 3;
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;
22 // -------------------------------------------------------------------------
23 int main( int argc, char* argv[] )
25 typedef TCLAP::ValueArg< std::string > _TStringArg;
26 typedef TCLAP::ValueArg< double > _TRealArg;
27 typedef TCLAP::ValueArg< TPixel > _TPixelArg;
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)" );
38 TCLAP::CmdLine cmd( "Labelling", ' ', "1.0.0" );
42 cmd.add( vesselness );
45 cmd.parse( argc, argv );
47 catch( TCLAP::ArgException& err )
50 << "===============================" << std::endl
51 << "Error caught: " << std::endl
52 << err.error( ) << " " << err.argId( ) << std::endl
53 << "===============================" << std::endl
62 TImage::Pointer input_image;
63 CTBronchi::ReadImage( input_image, in.getValue( ) );
66 TLabels::Pointer input_labels;
67 CTBronchi::ReadImage( input_labels, labels.getValue( ) );
69 // Read input vesselness
70 TScalarImage::Pointer input_vesselness;
71 CTBronchi::ReadImage( input_vesselness, vesselness.getValue( ) );
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;
88 CTBronchi::WriteImage( labelling->GetOutput( ), out.getValue( ) );
90 catch( std::exception& err )
93 << "===============================" << std::endl
94 << "Error caught: " << std::endl
95 << err.what( ) << std::endl
96 << "===============================" << std::endl
108 #include <itkImage.h>
109 #include <itkImageFileReader.h>
110 #include <itkImageFileWriter.h>
111 #include <itkRegionOfInterestImageFilter.h>
112 #include <itkMinimumMaximumImageCalculator.h>
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>
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;
129 class MoriLabellingTraits
130 : public fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel >
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;
138 typedef fpa::Filters::BaseMarksInterface< TInternalTraits > TMarksInterface;
139 typedef fpa::Filters::Image::SeedsFromLabelsInterface< TInternalTraits > TSeedsInterface;
143 : public fpa::Filters::Image::RegionGrow< TImage, TLabels, TLabel, MoriLabellingTraits >
146 typedef MoriLabellingTraits TTraits;
147 fpaTraitsMacro( TTraits );
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;
157 itkTypeMacro( MoriLabelling, fpa::Filters::Image::RegionGrow );
159 itkGetConstMacro( LastThreshold, TInputValue );
160 itkSetMacro( LastThreshold, TInputValue );
162 itkGetConstMacro( MinVertex, TVertex );
163 itkGetConstMacro( MaxVertex, TVertex );
165 ivqITKInputMacro( InputLabels, TLabels );
166 ivqITKInputMacro( InputVesselness, TScalarImage );
169 TInputValue GetUpperThreshold( ) const
171 return( this->m_Functor->GetUpperThreshold( ) );
173 void SetUpperThreshold( TInputValue t )
175 this->m_Functor->SetUpperThreshold( t );
181 m_LastThreshold( TInputValue( 0 ) ),
182 m_VesselnessThr( TScalar( 0.05 ) )
184 ivqITKInputConfigureMacro( InputLabels, TLabels );
185 ivqITKInputConfigureMacro( InputVesselness, TScalarImage );
186 this->m_Functor = TFunctor::New( );
187 this->SetPredicate( this->m_Functor );
190 virtual ~MoriLabelling( )
194 virtual const itk::DataObject* _GetReferenceInput( ) const override
196 return( this->GetInputLabels( ) );
198 virtual void _BeforeGenerateData( ) override
200 this->Superclass::_BeforeGenerateData( );
201 this->m_FirstVertex = true;
203 typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
204 _TMinMax::Pointer minMax = _TMinMax::New( );
205 minMax->SetImage( this->GetInputVesselness( ) );
207 this->m_MaxVesselness = ( 1.0 - this->m_VesselnessThr ) * minMax->GetMaximum( );
210 virtual void _PostComputeOutputValue( TNode& n ) override
212 this->Superclass::_PostComputeOutputValue( n );
213 if( n.Value == this->GetInsideValue( ) )
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 )
221 if( vesselness->GetPixel( n.Vertex ) > this->m_MaxVesselness )
222 n.Value = this->GetInsideValue( );
228 if( !( this->m_FirstVertex ) )
230 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
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 ];
241 this->m_MinVertex = n.Vertex;
242 this->m_MaxVertex = n.Vertex;
243 this->m_FirstVertex = false;
251 // Purposely not implemented.
252 MoriLabelling( const Self& other );
253 Self& operator=( const Self& other );
256 TFunctor::Pointer m_Functor;
257 TInputValue m_LastThreshold;
262 TScalar m_MaxVesselness;
263 TScalar m_VesselnessThr;
266 // -------------------------------------------------------------------------
267 double MeasureTime( itk::ProcessObject* f )
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( );
273 e = std::chrono::high_resolution_clock::now( );
275 return( t.count( ) );
278 // -------------------------------------------------------------------------
279 template< class _TImagePtr >
280 void ReadImage( _TImagePtr& image, const std::string& fname )
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( );
292 // -------------------------------------------------------------------------
293 template< class _TImage >
294 void WriteImage( const _TImage* image, const std::string& fname )
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;
304 // -------------------------------------------------------------------------
305 template< class _TImagePtr, class _TRegion >
308 const typename _TImagePtr::ObjectType* input,
312 typedef typename _TImagePtr::ObjectType _TImage;
313 typedef itk::RegionOfInterestImageFilter< _TImage, _TImage > _TFilter;
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( );
324 // -------------------------------------------------------------------------
326 std::map< std::string, std::string >& args, int argc, char* argv[]
329 std::set< std::string > mandatory;
330 mandatory.insert( "in" );
331 mandatory.insert( "out" );
332 mandatory.insert( "labels" );
333 mandatory.insert( "vesselness" );
335 args[ "upper_threshold" ] = "-650";
340 std::string cmd = argv[ i ] + 1;
341 args[ cmd ] = argv[ i + 1 ];
346 bool complete = true;
347 for( std::string t: mandatory )
348 complete &= ( args.find( t ) != args.end( ) );
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;
366 // -------------------------------------------------------------------------
367 int main( int argc, char* argv[] )
369 std::map< std::string, std::string > args;
372 if( ParseArgs( args, argc, argv ) )
375 TImage::Pointer input_image;
376 ReadImage( input_image, args[ "in" ] );
379 TLabels::Pointer input_labels;
380 ReadImage( input_labels, args[ "labels" ] );
382 // Read vesselness image
383 TScalarImage::Pointer input_vesselness;
384 ReadImage( input_vesselness, args[ "vesselness" ] );
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( ) ) )
397 double t = MeasureTime( labelling );
398 std::cout << "Labelling executed in " << t << " s" << std::endl;
400 std::map< std::string, std::string >::const_iterator i =
402 if( i != args.end( ) )
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 );
414 TImage::Pointer input_image_roi;
415 ROI( input_image_roi, input_image.GetPointer( ), roi );
416 WriteImage( input_image_roi.GetPointer( ), i->second );
419 TLabels::Pointer output_roi;
420 ROI( output_roi, labelling->GetOutput( ), roi );
421 WriteImage( output_roi.GetPointer( ), args[ "out" ] );
424 WriteImage( labelling->GetOutput( ), args[ "out" ] );
429 catch( std::exception& err )
432 << "===============================" << std::endl
433 << "Error caught: " << std::endl
435 << "===============================" << std::endl