3 #include <itkImageFileReader.h>
4 #include <itkImageFileWriter.h>
5 #include <fpa/Image/Mori.h>
7 // -------------------------------------------------------------------------
8 const unsigned int Dim = 3;
10 typedef unsigned short TLabel;
12 typedef itk::Image< TPixel, Dim > TInputImage;
13 typedef itk::Image< TLabel, Dim > TLabelImage;
14 typedef fpa::Image::Mori< TInputImage, TLabelImage > TFilter;
16 // -------------------------------------------------------------------------
17 int main( int argc, char* argv[] )
23 << "Usage: " << argv[ 0 ] << std::endl
24 << " input_image output_image output_signal" << std::endl
25 << " init_threshold end_threshold delta" << std::endl
26 << " minimum_threshold" << std::endl
27 << " inside_value outside_value" << std::endl
28 << " signal_kernel_size signal_threshold signal_influence"
30 << " seed_x seed_y seed_z"
35 std::string input_image_filename = argv[ 1 ];
36 std::string output_image_filename = argv[ 2 ];
37 std::string output_signal_filename = argv[ 3 ];
38 TPixel init_threshold = std::atoi( argv[ 4 ] );
39 TPixel end_threshold = std::atoi( argv[ 5 ] );
40 TPixel delta = std::atoi( argv[ 6 ] );
41 TPixel minimum_threshold = std::atoi( argv[ 7 ] );
42 TLabel inside_value = std::atoi( argv[ 8 ] );
43 TLabel outside_value = std::atoi( argv[ 9 ] );
44 unsigned long signal_kernel_size = std::atoi( argv[ 10 ] );
45 double signal_threshold = std::atof( argv[ 11 ] );
46 double signal_influence = std::atof( argv[ 12 ] );
47 TInputImage::PointType pseed;
48 for( int i = 0; i < Dim; ++i )
49 pseed[ i ] = std::atof( argv[ 13 + i ] );
52 itk::ImageFileReader< TInputImage >::Pointer input_image_reader =
53 itk::ImageFileReader< TInputImage >::New( );
54 input_image_reader->SetFileName( input_image_filename );
57 TFilter::Pointer filter = TFilter::New( );
58 filter->SetInput( input_image_reader->GetOutput( ) );
59 filter->SetSeed( pseed );
60 filter->SetThresholds( init_threshold, end_threshold, delta );
61 filter->SetMinimumThreshold( minimum_threshold );
62 filter->SetInsideValue( inside_value );
63 filter->SetOutsideValue( outside_value );
64 filter->SetSignalKernelSize( signal_kernel_size );
65 filter->SetSignalThreshold( signal_threshold );
66 filter->SetSignalInfluence( signal_influence );
68 // Show some information
69 std::cout << "----------------------------------------------" << std::endl;
70 std::cout << "Image: " << input_image_filename << std::endl;
73 std::chrono::time_point< std::chrono::high_resolution_clock > ts, te;
74 std::chrono::duration< double > tr;
77 ts = std::chrono::high_resolution_clock::now( );
78 input_image_reader->Update( );
79 te = std::chrono::high_resolution_clock::now( );
81 std::cout << "Read time: " << tr.count( ) << " s" << std::endl;
83 ts = std::chrono::high_resolution_clock::now( );
85 te = std::chrono::high_resolution_clock::now( );
88 << "Mori time: " << tr.count( ) << " s" << std::endl
89 << "Optimum threshold: " << filter->GetOptimumThreshold( ) << std::endl
90 << "Number of thresholds: "
91 << filter->GetNumberOfEvaluatedThresholds( ) << std::endl;
93 catch( std::exception& err )
95 std::cerr << "Error caught: " << err.what( ) << std::endl;
101 itk::ImageFileWriter< TLabelImage >::Pointer output_image_writer =
102 itk::ImageFileWriter< TLabelImage >::New( );
103 output_image_writer->SetInput( filter->GetThresholdedOutput( ) );
104 output_image_writer->SetFileName( output_image_filename );
107 ts = std::chrono::high_resolution_clock::now( );
108 output_image_writer->Update( );
109 te = std::chrono::high_resolution_clock::now( );
111 std::cout << "Write time: " << tr.count( ) << " s" << std::endl;
113 catch( std::exception& err )
115 std::cerr << "Error caught: " << err.what( ) << std::endl;
120 // Save output signal
121 std::ofstream osignal( output_signal_filename.c_str( ) );
123 unsigned long i = 0; i < filter->GetNumberOfEvaluatedThresholds( ); ++i
128 filter->GetSignalValues( i, x, y, p );
129 osignal << x << " " << y << std::endl;
133 std::cout << "----------------------------------------------" << std::endl;