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 <CTBronchi/MoriLabelling.h>
19 // -------------------------------------------------------------------------
20 const unsigned int Dim = 3;
21 typedef short TInputPixel;
22 typedef unsigned char TLabelPixel;
23 typedef itk::Image< TInputPixel, Dim > TInputImage;
24 typedef itk::Image< TLabelPixel, Dim > TLabelImage;
25 typedef std::map< std::string, std::string > TMap;
27 // -------------------------------------------------------------------------
29 TInputImage::PointType global_seed;
30 TLabelPixel inside_value = std::numeric_limits< TLabelPixel >::max( );
31 TLabelPixel outside_value = TLabelPixel( 0 );
33 // -------------------------------------------------------------------------
34 double MeasureTime( itk::ProcessObject* f )
36 std::chrono::time_point< std::chrono::high_resolution_clock > s, e;
37 std::chrono::duration< double > t;
38 s = std::chrono::high_resolution_clock::now( );
40 e = std::chrono::high_resolution_clock::now( );
45 // -------------------------------------------------------------------------
46 template< class _TImagePtr >
47 void ReadImage( _TImagePtr& image, const std::string& fname )
49 typedef typename _TImagePtr::ObjectType _TImage;
50 typedef itk::ImageFileReader< _TImage > _TReader;
51 typename _TReader::Pointer reader = _TReader::New( );
52 reader->SetFileName( fname );
53 double t = MeasureTime( reader );
54 std::cout << "Read " << fname << " in " << t << " s" << std::endl;
55 image = reader->GetOutput( );
56 image->DisconnectPipeline( );
59 // -------------------------------------------------------------------------
60 template< class _TImagePtr >
61 void WriteImage( const _TImagePtr& image, const std::string& fname )
63 typedef typename _TImagePtr::ObjectType _TImage;
64 typedef itk::ImageFileWriter< _TImage > _TWriter;
65 typename _TWriter::Pointer writer = _TWriter::New( );
66 writer->SetFileName( fname );
67 writer->SetInput( image );
68 double t = MeasureTime( writer );
69 std::cout << "Wrote " << fname << " in " << t << " s" << std::endl;
72 // -------------------------------------------------------------------------
73 template< class _TInputPtr, class _TOutputPtr >
74 void Mori( _TOutputPtr& output, const _TInputPtr& input )
76 typedef typename _TInputPtr::ObjectType _TInput;
77 typedef typename _TOutputPtr::ObjectType _TOutput;
78 typedef fpa::Filters::Image::Mori< _TInput, _TOutput > _TMori;
80 typename _TMori::Pointer mori = _TMori::New( );
81 mori->SetInput( input );
82 mori->SetSeed( global_seed );
83 mori->SetInsideValue( inside_value );
84 mori->SetOutsideValue( outside_value );
85 mori->SetMinimumThreshold(
86 TInputPixel( std::atof( Args[ "mori_minimum_threshold" ].c_str( ) ) )
88 mori->SetSignalKernelSize(
89 std::atoi( Args[ "mori_signal_kernel_size" ].c_str( ) )
91 mori->SetSignalThreshold(
92 std::atof( Args[ "mori_signal_threshold" ].c_str( ) )
94 mori->SetSignalInfluence(
95 std::atof( Args[ "mori_signal_influence" ].c_str( ) )
98 TInputPixel( std::atof( Args[ "mori_lower_threshold" ].c_str( ) ) ),
99 TInputPixel( std::atof( Args[ "mori_upper_threshold" ].c_str( ) ) ),
100 TInputPixel( std::atof( Args[ "mori_delta_threshold" ].c_str( ) ) )
102 double t = MeasureTime( mori );
103 std::cout << "Mori executed in " << t << " s" << std::endl;
104 output = mori->GetOutput( );
106 TMap::const_iterator i = Args.find( "out_mori" );
107 if( i != Args.end( ) )
108 WriteImage( output, i->second );
110 i = Args.find( "out_signal" );
111 if( i != Args.end( ) )
113 std::stringstream signal;
114 unsigned long nthr = mori->GetNumberOfEvaluatedThresholds( );
115 signal << "# nThr = " << nthr << std::endl;
116 signal << "# Opt = " << mori->GetOptimumThreshold( ) << std::endl;
117 for( unsigned long j = 0; j < nthr; ++j )
119 typename _TMori::TPeak p;
121 mori->GetSignalValues( j, x, y, p );
122 signal << x << " " << y << std::endl;
126 std::ofstream signals_str( i->second.c_str( ) );
127 signals_str << signal.str( );
128 signals_str.close( );
131 output->DisconnectPipeline( );
134 // -------------------------------------------------------------------------
135 template< class _TInputPtr, class _TOutputPtr >
136 void Label( _TOutputPtr& output, const _TInputPtr& input )
140 // -------------------------------------------------------------------------
141 bool ParseArgs( int argc, char* argv[] )
143 std::set< std::string > mandatory;
144 mandatory.insert( "in" );
145 mandatory.insert( "out_segmentation" );
146 mandatory.insert( "seed_x" );
147 mandatory.insert( "seed_y" );
148 mandatory.insert( "seed_z" );
150 Args[ "mori_minimum_threshold" ] = "-850";
151 Args[ "mori_signal_kernel_size" ] = "20";
152 Args[ "mori_signal_threshold" ] = "100";
153 Args[ "mori_signal_influence" ] = "0.5";
154 Args[ "mori_lower_threshold" ] = "-1024";
155 Args[ "mori_upper_threshold" ] = "0";
156 Args[ "mori_delta_threshold" ] = "1";
157 for( int i = 1; i < argc; i += 2 )
158 Args[ argv[ i ] + 1 ] = argv[ i + 1 ];
160 bool complete = true;
161 for( std::string t: mandatory )
162 complete &= ( Args.find( t ) != Args.end( ) );
167 << "Usage: " << argv[ 0 ] << std::endl
168 << "\t-in filename" << std::endl
169 << "\t-out_segmentation filename" << std::endl
170 << "\t-seed_x value -seed_y value -seed_z value" << std::endl
171 << "\t[-out_mori filename]" << std::endl
172 << "\t[-out_signal filename]" << std::endl
173 << "\t[-out_labels filename]" << std::endl
174 << "\t[-mori_minimum_threshold value]" << std::endl
175 << "\t[-mori_signal_kernel_size value]" << std::endl
176 << "\t[-mori_signal_threshold value]" << std::endl
177 << "\t[-mori_signal_influence value]" << std::endl
178 << "\t[-mori_lower_threshold value]" << std::endl
179 << "\t[-mori_upper_threshold value]" << std::endl
180 << "\t[-mori_delta_threshold value]" << std::endl;
187 // -------------------------------------------------------------------------
188 int main( int argc, char* argv[] )
192 if( ParseArgs( argc, argv ) )
195 global_seed[ 0 ] = std::atof( Args[ "seed_x" ].c_str( ) );
196 global_seed[ 1 ] = std::atof( Args[ "seed_y" ].c_str( ) );
197 global_seed[ 2 ] = std::atof( Args[ "seed_z" ].c_str( ) );
200 TInputImage::Pointer input_image;
201 ReadImage( input_image, Args[ "in" ] );
204 TLabelImage::Pointer mori;
205 Mori( mori, input_image );
208 TLabelImage::Pointer labels;
209 Label( labels, mori );
216 catch( std::exception& err )
219 << "===============================" << std::endl
220 << "Error caught: " << std::endl
222 << "===============================" << std::endl
228 if [ -z "$label_upper_threshold" ]; then label_upper_threshold="-600"; fi
229 if [ -z "$label_inside" ]; then label_inside="1"; fi
230 if [ -z "$label_outside" ]; then label_outside="2"; fi
231 if [ -z "$random_walker_beta" ]; then random_walker_beta="20"; fi