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(-1024) end_threshold(0) delta(1)" << std::endl
26 << " minimum_threshold(-850)" << std::endl
27 << " inside_value(255) outside_value(0)" << std::endl
28 << " signal_kernel_size(20) signal_threshold(500) signal_influence(0.5)"
30 << " [index/point] 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 std::string seed_type = argv[ 13 ];
48 TInputImage::IndexType iseed;
49 TInputImage::PointType pseed;
50 for( unsigned int i = 0; i < Dim; ++i )
52 if( seed_type == "index" )
53 iseed[ i ] = std::atoi( argv[ 14 + i ] );
55 pseed[ i ] = std::atof( argv[ 14 + i ] );
60 itk::ImageFileReader< TInputImage >::Pointer input_image_reader =
61 itk::ImageFileReader< TInputImage >::New( );
62 input_image_reader->SetFileName( input_image_filename );
65 TFilter::Pointer filter = TFilter::New( );
66 filter->SetInput( input_image_reader->GetOutput( ) );
67 if( seed_type == "index" )
68 filter->SetSeed( iseed );
70 filter->SetSeed( pseed );
71 filter->SetThresholds( init_threshold, end_threshold, delta );
72 filter->SetMinimumThreshold( minimum_threshold );
73 filter->SetInsideValue( inside_value );
74 filter->SetOutsideValue( outside_value );
75 filter->SetSignalKernelSize( signal_kernel_size );
76 filter->SetSignalThreshold( signal_threshold );
77 filter->SetSignalInfluence( signal_influence );
79 // Show some information
80 std::cout << "----------------------------------------------" << std::endl;
81 std::cout << "Image: " << input_image_filename << std::endl;
84 std::chrono::time_point< std::chrono::high_resolution_clock > ts, te;
85 std::chrono::duration< double > tr;
88 ts = std::chrono::high_resolution_clock::now( );
89 input_image_reader->Update( );
90 te = std::chrono::high_resolution_clock::now( );
92 std::cout << "Read time: " << tr.count( ) << " s" << std::endl;
94 ts = std::chrono::high_resolution_clock::now( );
96 te = std::chrono::high_resolution_clock::now( );
99 << "Mori time: " << tr.count( ) << " s" << std::endl
100 << "Optimum threshold: " << filter->GetOptimumThreshold( ) << std::endl
101 << "Number of thresholds: "
102 << filter->GetNumberOfEvaluatedThresholds( ) << std::endl;
104 catch( std::exception& err )
106 std::cerr << "Error caught: " << err.what( ) << std::endl;
112 itk::ImageFileWriter< TLabelImage >::Pointer output_image_writer =
113 itk::ImageFileWriter< TLabelImage >::New( );
114 output_image_writer->SetInput( filter->GetThresholdedOutput( ) );
115 output_image_writer->SetFileName( output_image_filename );
118 ts = std::chrono::high_resolution_clock::now( );
119 output_image_writer->Update( );
120 te = std::chrono::high_resolution_clock::now( );
122 std::cout << "Write time: " << tr.count( ) << " s" << std::endl;
124 catch( std::exception& err )
126 std::cerr << "Error caught: " << err.what( ) << std::endl;
131 // Save output signal
132 std::ofstream osignal( output_signal_filename.c_str( ) );
134 unsigned long i = 0; i < filter->GetNumberOfEvaluatedThresholds( ); ++i
139 filter->GetSignalValues( i, x, y, p );
140 osignal << x << " " << y << std::endl;
144 std::cout << "----------------------------------------------" << std::endl;