1 // =========================================================================
2 // @author Leonardo Florez Valencia
3 // @email florez-l@javeriana.edu.co
4 // =========================================================================
10 #include <itkImageFileReader.h>
11 #include <itkImageFileWriter.h>
12 #include <itkRegionOfInterestImageFilter.h>
13 #include <itkMinimumMaximumImageCalculator.h>
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>
21 // -------------------------------------------------------------------------
22 const unsigned int Dim = 3;
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;
32 class MoriLabellingTraits
33 : public fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel >
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;
41 typedef fpa::Filters::BaseMarksInterface< TInternalTraits > TMarksInterface;
42 typedef fpa::Filters::Image::SeedsFromLabelsInterface< TInternalTraits > TSeedsInterface;
48 : public fpa::Filters::Image::RegionGrow< TImage, TLabels, TLabel, MoriLabellingTraits >
51 typedef MoriLabellingTraits TTraits;
52 fpaTraitsMacro( TTraits );
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;
62 itkTypeMacro( MoriLabelling, fpa::Filters::Image::RegionGrow );
64 itkGetConstMacro( LastThreshold, TInputValue );
65 itkSetMacro( LastThreshold, TInputValue );
67 itkGetConstMacro( MinVertex, TVertex );
68 itkGetConstMacro( MaxVertex, TVertex );
70 ivqITKInputMacro( InputLabels, TLabels );
71 ivqITKInputMacro( InputVesselness, TScalarImage );
74 TInputValue GetUpperThreshold( ) const
76 return( this->m_Functor->GetUpperThreshold( ) );
78 void SetUpperThreshold( TInputValue t )
80 this->m_Functor->SetUpperThreshold( t );
86 m_LastThreshold( TInputValue( 0 ) ),
87 m_VesselnessThr( TScalar( 0.05 ) )
89 ivqITKInputConfigureMacro( InputLabels, TLabels );
90 ivqITKInputConfigureMacro( InputVesselness, TScalarImage );
91 this->m_Functor = TFunctor::New( );
92 this->SetPredicate( this->m_Functor );
95 virtual ~MoriLabelling( )
99 virtual const itk::DataObject* _GetReferenceInput( ) const override
101 return( this->GetInputLabels( ) );
103 virtual void _BeforeGenerateData( ) override
105 this->Superclass::_BeforeGenerateData( );
106 this->m_FirstVertex = true;
108 typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
109 _TMinMax::Pointer minMax = _TMinMax::New( );
110 minMax->SetImage( this->GetInputVesselness( ) );
112 this->m_MaxVesselness = ( 1.0 - this->m_VesselnessThr ) * minMax->GetMaximum( );
115 virtual void _PostComputeOutputValue( TNode& n ) override
117 this->Superclass::_PostComputeOutputValue( n );
118 if( n.Value == this->GetInsideValue( ) )
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 );
125 double a = std::fabs( x - double( this->m_Functor->GetUpperThreshold( ) ) );
126 double b = std::fabs( x - double( this->m_LastThreshold ) );
128 if( labels->GetPixel( n.Vertex ) == 0 )
130 if( vesselness->GetPixel( n.Vertex ) > this->m_MaxVesselness )
131 n.Value = this->GetInsideValue( );
137 if( !( this->m_FirstVertex ) )
139 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
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 ];
150 this->m_MinVertex = n.Vertex;
151 this->m_MaxVertex = n.Vertex;
152 this->m_FirstVertex = false;
160 // Purposely not implemented.
161 MoriLabelling( const Self& other );
162 Self& operator=( const Self& other );
165 TFunctor::Pointer m_Functor;
166 TInputValue m_LastThreshold;
171 TScalar m_MaxVesselness;
172 TScalar m_VesselnessThr;
175 // -------------------------------------------------------------------------
176 double MeasureTime( itk::ProcessObject* f )
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( );
182 e = std::chrono::high_resolution_clock::now( );
184 return( t.count( ) );
187 // -------------------------------------------------------------------------
188 template< class _TImagePtr >
189 void ReadImage( _TImagePtr& image, const std::string& fname )
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( );
201 // -------------------------------------------------------------------------
202 template< class _TImage >
203 void WriteImage( const _TImage* image, const std::string& fname )
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;
213 // -------------------------------------------------------------------------
214 template< class _TImagePtr, class _TRegion >
217 const typename _TImagePtr::ObjectType* input,
221 typedef typename _TImagePtr::ObjectType _TImage;
222 typedef itk::RegionOfInterestImageFilter< _TImage, _TImage > _TFilter;
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( );
233 // -------------------------------------------------------------------------
235 std::map< std::string, std::string >& args, int argc, char* argv[]
238 std::set< std::string > mandatory;
239 mandatory.insert( "in" );
240 mandatory.insert( "out" );
241 mandatory.insert( "labels" );
242 mandatory.insert( "vesselness" );
244 args[ "upper_threshold" ] = "-650";
249 std::string cmd = argv[ i ] + 1;
250 args[ cmd ] = argv[ i + 1 ];
255 bool complete = true;
256 for( std::string t: mandatory )
257 complete &= ( args.find( t ) != args.end( ) );
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;
275 // -------------------------------------------------------------------------
276 int main( int argc, char* argv[] )
278 std::map< std::string, std::string > args;
281 if( ParseArgs( args, argc, argv ) )
284 TImage::Pointer input_image;
285 ReadImage( input_image, args[ "in" ] );
288 TLabels::Pointer input_labels;
289 ReadImage( input_labels, args[ "labels" ] );
291 // Read vesselness image
292 TScalarImage::Pointer input_vesselness;
293 ReadImage( input_vesselness, args[ "vesselness" ] );
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( ) ) )
307 labelling->SetLastThreshold( last_thr );
309 double t = MeasureTime( labelling );
310 std::cout << "Labelling executed in " << t << " s" << std::endl;
312 std::map< std::string, std::string >::const_iterator i =
314 if( i != args.end( ) )
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 );
326 TImage::Pointer input_image_roi;
327 ROI( input_image_roi, input_image.GetPointer( ), roi );
328 WriteImage( input_image_roi.GetPointer( ), i->second );
331 TLabels::Pointer output_roi;
332 ROI( output_roi, labelling->GetOutput( ), roi );
333 WriteImage( output_roi.GetPointer( ), args[ "out" ] );
336 WriteImage( labelling->GetOutput( ), args[ "out" ] );
341 catch( std::exception& err )
344 << "===============================" << std::endl
345 << "Error caught: " << std::endl
347 << "===============================" << std::endl