1 // =========================================================================
2 // @author Leonardo Florez Valencia
3 // @email florez-l@javeriana.edu.co
4 // =========================================================================
14 #include <itkImageFileReader.h>
15 #include <itkImageFileWriter.h>
16 #include <fpa/Filters/Image/Mori.h>
17 #include <fpa/Filters/Image/RandomWalker.h>
18 #include <fpa/Functors/Dijkstra/Image/Gaussian.h>
19 #include <CTBronchi/MoriLabelling.h>
21 // -------------------------------------------------------------------------
22 const unsigned int Dim = 3;
23 typedef short TInputPixel;
24 typedef unsigned char TLabelPixel;
25 typedef float TScalar;
26 typedef itk::Image< TInputPixel, Dim > TInputImage;
27 typedef itk::Image< TLabelPixel, Dim > TLabelImage;
28 typedef itk::Image< TScalar, Dim > TScalarImage;
29 typedef std::map< std::string, std::string > TMap;
31 // -------------------------------------------------------------------------
33 TInputImage::PointType global_seed;
34 TLabelPixel inside_value = TLabelPixel( 1 );
35 TLabelPixel outside_value = TLabelPixel( 0 );
36 TLabelPixel inside_label = TLabelPixel( 1 );
37 TLabelPixel outside_label = TLabelPixel( 2 );
39 // -------------------------------------------------------------------------
40 double MeasureTime( itk::ProcessObject* f )
42 std::chrono::time_point< std::chrono::high_resolution_clock > s, e;
43 std::chrono::duration< double > t;
44 s = std::chrono::high_resolution_clock::now( );
46 e = std::chrono::high_resolution_clock::now( );
51 // -------------------------------------------------------------------------
52 template< class _TImagePtr >
53 void ReadImage( _TImagePtr& image, const std::string& fname )
55 typedef typename _TImagePtr::ObjectType _TImage;
56 typedef itk::ImageFileReader< _TImage > _TReader;
57 typename _TReader::Pointer reader = _TReader::New( );
58 reader->SetFileName( fname );
59 double t = MeasureTime( reader );
60 std::cout << "Read " << fname << " in " << t << " s" << std::endl;
61 image = reader->GetOutput( );
62 image->DisconnectPipeline( );
65 // -------------------------------------------------------------------------
66 template< class _TImagePtr >
67 void WriteImage( const _TImagePtr& image, const std::string& fname )
69 typedef typename _TImagePtr::ObjectType _TImage;
70 typedef itk::ImageFileWriter< _TImage > _TWriter;
71 typename _TWriter::Pointer writer = _TWriter::New( );
72 writer->SetFileName( fname );
73 writer->SetInput( image );
74 double t = MeasureTime( writer );
75 std::cout << "Wrote " << fname << " in " << t << " s" << std::endl;
78 // -------------------------------------------------------------------------
79 template< class _TInputPtr, class _TOutputPtr >
80 void Mori( _TOutputPtr& output, const _TInputPtr& input )
82 typedef typename _TInputPtr::ObjectType _TInput;
83 typedef typename _TOutputPtr::ObjectType _TOutput;
84 typedef fpa::Filters::Image::Mori< _TInput, _TOutput > _TMori;
86 typename _TMori::Pointer mori = _TMori::New( );
87 mori->SetInput( input );
88 mori->SetSeed( global_seed );
89 mori->SetInsideValue( inside_value );
90 mori->SetOutsideValue( outside_value );
91 mori->SetMinimumThreshold(
92 TInputPixel( std::atof( Args[ "mori_minimum_threshold" ].c_str( ) ) )
94 mori->SetSignalKernelSize(
95 std::atoi( Args[ "mori_signal_kernel_size" ].c_str( ) )
97 mori->SetSignalThreshold(
98 std::atof( Args[ "mori_signal_threshold" ].c_str( ) )
100 mori->SetSignalInfluence(
101 std::atof( Args[ "mori_signal_influence" ].c_str( ) )
104 TInputPixel( std::atof( Args[ "mori_lower_threshold" ].c_str( ) ) ),
105 TInputPixel( std::atof( Args[ "mori_upper_threshold" ].c_str( ) ) ),
106 TInputPixel( std::atof( Args[ "mori_delta_threshold" ].c_str( ) ) )
108 double t = MeasureTime( mori );
109 std::cout << "Mori executed in " << t << " s" << std::endl;
110 output = mori->GetOutput( );
112 TMap::const_iterator i = Args.find( "out_mori" );
113 if( i != Args.end( ) )
114 WriteImage( output, i->second );
116 i = Args.find( "out_signal" );
117 if( i != Args.end( ) )
119 std::stringstream signal;
120 unsigned long nthr = mori->GetNumberOfEvaluatedThresholds( );
121 signal << "# nThr = " << nthr << std::endl;
122 signal << "# Opt = " << mori->GetOptimumThreshold( ) << std::endl;
123 for( unsigned long j = 0; j < nthr; ++j )
125 typename _TMori::TPeak p;
127 mori->GetSignalValues( j, x, y, p );
128 signal << x << " " << y << std::endl;
132 std::ofstream signals_str( i->second.c_str( ) );
133 signals_str << signal.str( );
134 signals_str.close( );
137 output->DisconnectPipeline( );
140 // -------------------------------------------------------------------------
141 template< class _TInputPtr, class _TOutputPtr >
142 void Label( _TOutputPtr& output, const _TInputPtr& input, const _TOutputPtr& labels )
144 typedef typename _TInputPtr::ObjectType _TInput;
145 typedef typename _TOutputPtr::ObjectType _TOutput;
146 typedef CTBronchi::MoriLabelling< _TInput, _TOutput > _TLabelling;
148 typename _TLabelling::Pointer labelling = _TLabelling::New( );
149 labelling->SetInput( input );
150 labelling->SetInputLabels( labels );
151 labelling->SetOutsideValue( outside_label );
152 labelling->SetInsideValue( inside_label );
153 labelling->SetInputInsideValue( inside_value );
154 labelling->SetUpperThreshold(
155 TInputPixel( std::atof( Args[ "labelling_upper_threshold" ].c_str( ) ) )
157 double t = MeasureTime( labelling );
158 std::cout << "Labelling executed in " << t << " s" << std::endl;
159 output = labelling->GetOutput( );
160 TMap::const_iterator i = Args.find( "out_labels" );
161 if( i != Args.end( ) )
162 WriteImage( output, i->second );
163 output->DisconnectPipeline( );
166 // -------------------------------------------------------------------------
167 template< class _TRawPtr, class _TLabelPtr >
168 void RandomWalker( _TLabelPtr& output, const _TRawPtr& raw, const _TLabelPtr& labels )
170 typedef typename _TRawPtr::ObjectType _TRaw;
171 typedef typename _TLabelPtr::ObjectType _TLabel;
174 // -------------------------------------------------------------------------
175 bool ParseArgs( int argc, char* argv[] )
177 std::set< std::string > mandatory;
178 mandatory.insert( "in" );
179 mandatory.insert( "out_segmentation" );
180 mandatory.insert( "seed_x" );
181 mandatory.insert( "seed_y" );
182 mandatory.insert( "seed_z" );
184 Args[ "mori_minimum_threshold" ] = "-850";
185 Args[ "mori_signal_kernel_size" ] = "20";
186 Args[ "mori_signal_threshold" ] = "100";
187 Args[ "mori_signal_influence" ] = "0.5";
188 Args[ "mori_lower_threshold" ] = "-1024";
189 Args[ "mori_upper_threshold" ] = "0";
190 Args[ "mori_delta_threshold" ] = "1";
191 Args[ "labelling_upper_threshold" ] = "-400";
193 for( int i = 1; i < argc; i += 2 )
194 Args[ argv[ i ] + 1 ] = argv[ i + 1 ];
196 bool complete = true;
197 for( std::string t: mandatory )
198 complete &= ( Args.find( t ) != Args.end( ) );
203 << "Usage: " << argv[ 0 ] << std::endl
204 << "\t-in filename" << std::endl
205 << "\t-out_segmentation filename" << std::endl
206 << "\t-seed_x value -seed_y value -seed_z value" << std::endl
207 << "\t[-out_mori filename]" << std::endl
208 << "\t[-out_signal filename]" << std::endl
209 << "\t[-out_labels filename]" << std::endl
210 << "\t[-mori_minimum_threshold value]" << std::endl
211 << "\t[-mori_signal_kernel_size value]" << std::endl
212 << "\t[-mori_signal_threshold value]" << std::endl
213 << "\t[-mori_signal_influence value]" << std::endl
214 << "\t[-mori_lower_threshold value]" << std::endl
215 << "\t[-mori_upper_threshold value]" << std::endl
216 << "\t[-mori_delta_threshold value]" << std::endl
217 << "\t[-labelling_upper_threshold value]" << std::endl;
224 // -------------------------------------------------------------------------
225 int main( int argc, char* argv[] )
229 if( ParseArgs( argc, argv ) )
232 global_seed[ 0 ] = std::atof( Args[ "seed_x" ].c_str( ) );
233 global_seed[ 1 ] = std::atof( Args[ "seed_y" ].c_str( ) );
234 global_seed[ 2 ] = std::atof( Args[ "seed_z" ].c_str( ) );
237 TInputImage::Pointer input_image;
238 ReadImage( input_image, Args[ "in" ] );
241 TLabelImage::Pointer mori;
242 Mori( mori, input_image );
245 TLabelImage::Pointer labels;
246 Label( labels, input_image, mori );
253 catch( std::exception& err )
256 << "===============================" << std::endl
257 << "Error caught: " << std::endl
259 << "===============================" << std::endl
265 if [ -z "$label_upper_threshold" ]; then label_upper_threshold="-600"; fi
266 if [ -z "$label_inside" ]; then label_inside="1"; fi
267 if [ -z "$label_outside" ]; then label_outside="2"; fi
268 if [ -z "$random_walker_beta" ]; then random_walker_beta="20"; fi