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>
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>
20 // -------------------------------------------------------------------------
21 const unsigned int Dim = 3;
23 typedef unsigned char TLabel;
24 typedef itk::Image< TPixel, Dim > TImage;
25 typedef itk::Image< TLabel, Dim > TLabels;
29 class MoriLabellingTraits
30 : public fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel >
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;
38 typedef fpa::Filters::BaseMarksInterface< TInternalTraits > TMarksInterface;
39 typedef fpa::Filters::Image::SeedsFromLabelsInterface< TInternalTraits > TSeedsInterface;
45 : public fpa::Filters::Image::RegionGrow< TImage, TLabels, TLabel, MoriLabellingTraits >
48 typedef MoriLabellingTraits TTraits;
49 fpaTraitsMacro( TTraits );
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;
59 itkTypeMacro( MoriLabelling, fpa::Filters::Image::RegionGrow );
61 itkGetConstMacro( LastThreshold, TInputValue );
62 itkSetMacro( LastThreshold, TInputValue );
64 itkGetConstMacro( MinVertex, TVertex );
65 itkGetConstMacro( MaxVertex, TVertex );
67 fpaFilterInputMacro( InputLabels, TLabels );
70 TInputValue GetUpperThreshold( ) const
72 return( this->m_Functor->GetUpperThreshold( ) );
74 void SetUpperThreshold( TInputValue t )
76 this->m_Functor->SetUpperThreshold( t );
82 m_LastThreshold( TInputValue( 0 ) )
84 fpaFilterInputConfigureMacro( InputLabels, TLabels );
85 this->m_Functor = TFunctor::New( );
86 this->SetPredicate( this->m_Functor );
89 virtual ~MoriLabelling( )
93 virtual const itk::DataObject* _GetReferenceInput( ) const override
95 return( this->GetInputLabels( ) );
97 virtual void _BeforeGenerateData( ) override
99 this->Superclass::_BeforeGenerateData( );
100 this->m_FirstVertex = true;
103 virtual void _PostComputeOutputValue( TNode& n ) override
105 this->Superclass::_PostComputeOutputValue( n );
106 if( n.Value == this->GetInsideValue( ) )
108 const TImage* input = this->GetInput( );
109 const TLabels* labels = this->GetInputLabels( );
110 double x = input->GetPixel( n.Vertex );
112 double a = std::fabs( x - double( this->m_Functor->GetUpperThreshold( ) ) );
113 double b = std::fabs( x - double( this->m_LastThreshold ) );
115 if( labels->GetPixel( n.Vertex ) == 0 /* && b < a*/ )
118 if( !( this->m_FirstVertex ) )
120 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
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 ];
131 this->m_MinVertex = n.Vertex;
132 this->m_MaxVertex = n.Vertex;
133 this->m_FirstVertex = false;
141 // Purposely not implemented.
142 MoriLabelling( const Self& other );
143 Self& operator=( const Self& other );
146 TFunctor::Pointer m_Functor;
147 TInputValue m_LastThreshold;
153 // -------------------------------------------------------------------------
154 double MeasureTime( itk::ProcessObject* f )
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( );
160 e = std::chrono::high_resolution_clock::now( );
162 return( t.count( ) );
165 // -------------------------------------------------------------------------
166 template< class _TImagePtr >
167 void ReadImage( _TImagePtr& image, const std::string& fname )
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( );
179 // -------------------------------------------------------------------------
180 template< class _TImage >
181 void WriteImage( const _TImage* image, const std::string& fname )
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;
191 // -------------------------------------------------------------------------
192 template< class _TImagePtr, class _TRegion >
195 const typename _TImagePtr::ObjectType* input,
199 typedef typename _TImagePtr::ObjectType _TImage;
200 typedef itk::RegionOfInterestImageFilter< _TImage, _TImage > _TFilter;
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( );
211 // -------------------------------------------------------------------------
213 std::map< std::string, std::string >& args, int argc, char* argv[]
216 std::set< std::string > mandatory;
217 mandatory.insert( "in" );
218 mandatory.insert( "out" );
219 mandatory.insert( "labels" );
221 args[ "upper_threshold" ] = "-600";
226 std::string cmd = argv[ i ] + 1;
227 args[ cmd ] = argv[ i + 1 ];
232 bool complete = true;
233 for( std::string t: mandatory )
234 complete &= ( args.find( t ) != args.end( ) );
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;
251 // -------------------------------------------------------------------------
252 int main( int argc, char* argv[] )
254 std::map< std::string, std::string > args;
257 if( ParseArgs( args, argc, argv ) )
260 TImage::Pointer input_image;
261 ReadImage( input_image, args[ "in" ] );
264 TLabels::Pointer input_labels;
265 ReadImage( input_labels, args[ "labels" ] );
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( ) ) )
277 labelling->SetLastThreshold( last_thr );
279 double t = MeasureTime( labelling );
280 std::cout << "Labelling executed in " << t << " s" << std::endl;
282 std::map< std::string, std::string >::const_iterator i =
284 if( i != args.end( ) )
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 );
296 TImage::Pointer input_image_roi;
297 ROI( input_image_roi, input_image.GetPointer( ), roi );
298 WriteImage( input_image_roi.GetPointer( ), i->second );
301 TLabels::Pointer output_roi;
302 ROI( output_roi, labelling->GetOutput( ), roi );
303 WriteImage( output_roi.GetPointer( ), args[ "out" ] );
306 WriteImage( labelling->GetOutput( ), args[ "out" ] );
311 catch( std::exception& err )
314 << "===============================" << std::endl
315 << "Error caught: " << std::endl
317 << "===============================" << std::endl