1 // =========================================================================
2 // @author Leonardo Florez Valencia
3 // @email florez-l@javeriana.edu.co
4 // =========================================================================
9 #include <itkImageFileReader.h>
10 #include <itkImageFileWriter.h>
11 #include <fpa/Filters/Image/Mori.h>
13 // -------------------------------------------------------------------------
14 const unsigned int Dim = 3;
16 typedef unsigned char TLabel;
17 typedef itk::Image< TPixel, Dim > TImage;
18 typedef itk::Image< TLabel, Dim > TLabels;
20 // -------------------------------------------------------------------------
21 double MeasureTime( itk::ProcessObject* f )
23 std::chrono::time_point< std::chrono::high_resolution_clock > s, e;
24 std::chrono::duration< double > t;
25 s = std::chrono::high_resolution_clock::now( );
27 e = std::chrono::high_resolution_clock::now( );
32 // -------------------------------------------------------------------------
33 template< class _TImagePtr >
34 void ReadImage( _TImagePtr& image, const std::string& fname )
36 typedef typename _TImagePtr::ObjectType _TImage;
37 typedef itk::ImageFileReader< _TImage > _TReader;
38 typename _TReader::Pointer reader = _TReader::New( );
39 reader->SetFileName( fname );
40 double t = MeasureTime( reader );
41 std::cout << "Read " << fname << " in " << t << " s" << std::endl;
42 image = reader->GetOutput( );
43 image->DisconnectPipeline( );
46 // -------------------------------------------------------------------------
47 template< class _TImage >
48 void WriteImage( const _TImage* image, const std::string& fname )
50 typedef itk::ImageFileWriter< _TImage > _TWriter;
51 typename _TWriter::Pointer writer = _TWriter::New( );
52 writer->SetFileName( fname );
53 writer->SetInput( image );
54 double t = MeasureTime( writer );
55 std::cout << "Wrote " << fname << " in " << t << " s" << std::endl;
58 // -------------------------------------------------------------------------
60 std::map< std::string, std::string >& args, int argc, char* argv[]
63 std::set< std::string > mandatory;
64 mandatory.insert( "in" );
65 mandatory.insert( "out" );
66 mandatory.insert( "seed" );
68 args[ "minimum_threshold" ] = "-850";
69 args[ "signal_kernel_size" ] = "20";
70 args[ "signal_threshold" ] = "100";
71 args[ "signal_influence" ] = "0.5";
72 args[ "lower_threshold" ] = "-1024";
73 args[ "upper_threshold" ] = "0";
74 args[ "delta_threshold" ] = "1";
79 std::string cmd = argv[ i ] + 1;
82 std::stringstream seed;
84 << argv[ i + 1 ] << ";"
85 << argv[ i + 2 ] << ";"
87 args[ cmd ] = seed.str( );
92 args[ cmd ] = argv[ i + 1 ];
100 for( std::string t: mandatory )
101 complete &= ( args.find( t ) != args.end( ) );
106 << "Usage: " << argv[ 0 ] << std::endl
107 << "\t-in filename" << std::endl
108 << "\t-out filename" << std::endl
109 << "\t-seed x y z" << std::endl
110 << "\t[-out_signal filename]" << std::endl
111 << "\t[-minimum_threshold value]" << std::endl
112 << "\t[-signal_kernel_size value]" << std::endl
113 << "\t[-signal_threshold value]" << std::endl
114 << "\t[-signal_influence value]" << std::endl
115 << "\t[-lower_threshold value]" << std::endl
116 << "\t[-upper_threshold value]" << std::endl
117 << "\t[-delta_threshold value]" << std::endl;
124 // -------------------------------------------------------------------------
125 int main( int argc, char* argv[] )
127 std::map< std::string, std::string > args;
130 if( ParseArgs( args, argc, argv ) )
133 TImage::PointType seed;
134 char* str = new char[ args[ "seed" ].size( ) + 1 ];
135 strcpy( str, args[ "seed" ].c_str( ) );
136 seed[ 0 ] = std::atof( strtok( str, ";" ) );
137 seed[ 1 ] = std::atof( strtok( NULL, ";" ) );
138 seed[ 2 ] = std::atof( strtok( NULL, ";" ) );
141 TImage::Pointer input_image;
142 ReadImage( input_image, args[ "in" ] );
145 typedef fpa::Filters::Image::Mori< TImage, TLabels > TMori;
146 TMori::Pointer mori = TMori::New( );
147 mori->SetInput( input_image );
148 mori->SetSeed( seed );
149 mori->SetInsideValue( 1 );
150 mori->SetOutsideValue( 0 );
151 mori->SetMinimumThreshold(
152 TPixel( std::atof( args[ "minimum_threshold" ].c_str( ) ) )
154 mori->SetSignalKernelSize(
155 std::atoi( args[ "signal_kernel_size" ].c_str( ) )
157 mori->SetSignalThreshold(
158 std::atof( args[ "signal_threshold" ].c_str( ) )
160 mori->SetSignalInfluence(
161 std::atof( args[ "signal_influence" ].c_str( ) )
164 TPixel( std::atof( args[ "lower_threshold" ].c_str( ) ) ),
165 TPixel( std::atof( args[ "upper_threshold" ].c_str( ) ) ),
166 TPixel( std::atof( args[ "delta_threshold" ].c_str( ) ) )
168 double t = MeasureTime( mori );
169 std::cout << "Mori executed in " << t << " s" << std::endl;
170 WriteImage( mori->GetOutput( ), args[ "out" ] );
172 std::map< std::string, std::string >::const_iterator i =
173 args.find( "out_signal" );
174 if( i != args.end( ) )
176 std::stringstream signal;
177 unsigned long nthr = mori->GetNumberOfEvaluatedThresholds( );
178 signal << "# nThr = " << nthr << std::endl;
179 signal << "# Opt = " << mori->GetOptimumThreshold( ) << std::endl;
180 for( unsigned long j = 0; j < nthr; ++j )
182 typename TMori::TPeak p;
184 mori->GetSignalValues( j, x, y, p );
185 signal << x << " " << y << std::endl;
188 std::ofstream signals_str( i->second.c_str( ) );
189 signals_str << signal.str( );
190 signals_str.close( );
197 catch( std::exception& err )
200 << "===============================" << std::endl
201 << "Error caught: " << std::endl
203 << "===============================" << std::endl