8 #include <itkImageFileReader.h>
9 #include <itkImageFileWriter.h>
11 #include <fpa/Image/MoriFilter.h>
13 // -------------------------------------------------------------------------
14 static const unsigned int Dim = 3;
16 typedef itk::Image< TPixel, Dim > TImage;
18 // -------------------------------------------------------------------------
19 std::string GetFullPath( const std::string& filename )
23 TCHAR full_path[MAX_PATH];
24 GetFullPathName(_T("foo.dat"), MAX_PATH, full_path, NULL);
27 char* buffer = realpath( filename.c_str( ), NULL );
28 std::string path = buffer;
33 // -------------------------------------------------------------------------
34 std::string GetFullPathToDirectory( const std::string& filename )
36 std::string path = GetFullPath( filename );
37 std::size_t found = path.find_last_of( "/\\" );
38 return( path.substr( 0, found + 1 ) );
41 // -------------------------------------------------------------------------
42 int main( int argc, char* argv[] )
44 // Command configuration
48 << "Usage: " << argv[ 0 ] << " input_image output_image"
53 std::string input_image_filename = GetFullPath( argv[ 1 ] );
54 std::string output_image_filename = argv[ 2 ];
55 std::string output_auxiliary_image_filename = output_image_filename + "_aux.mhd";
57 // Try to guess initial seed
58 std::string seed_filename = GetFullPathToDirectory( argv[ 1 ] ) + "seed.txt";
59 std::ifstream seed_str( seed_filename.c_str( ) );
63 << "No \"seed.txt\" file in the same input image directory."
68 TImage::IndexType seed;
69 seed_str >> seed[ 0 ] >> seed[ 1 ] >> seed[ 2 ];
73 typedef itk::ImageFileReader< TImage > TReader;
74 TReader::Pointer reader = TReader::New( );
75 reader->SetFileName( input_image_filename );
78 typedef fpa::Image::MoriFilter< TImage, TImage > TFilter;
79 TFilter::Pointer filter = TFilter::New( );
80 filter->SetInput( reader->GetOutput( ) );
81 filter->SetSeed( seed );
82 filter->SetThresholdRange( -200, 100, 1 );
83 filter->SetInsideValue( 1 );
84 filter->SetOutsideValue( 0 );
87 typedef itk::ImageFileWriter< TImage > TWriter;
88 TWriter::Pointer writer = TWriter::New( );
89 writer->SetInput( filter->GetOutput( ) );
90 writer->SetFileName( output_image_filename );
97 catch( std::exception& err )
99 std::cerr << "Error caught: " << err.what( ) << std::endl;
104 // Write auxiliary image
105 typedef itk::ImageFileWriter< TFilter::TAuxImage > TAuxWriter;
106 TAuxWriter::Pointer aux_writer = TAuxWriter::New( );
107 aux_writer->SetInput( filter->GetAuxiliaryImage( ) );
108 aux_writer->SetFileName( output_auxiliary_image_filename );
113 aux_writer->Update( );
115 catch( std::exception& err )
117 std::cerr << "Error caught: " << err.what( ) << std::endl;
122 // Write result signal
123 std::string output_signal_filename = output_image_filename + ".csv";
124 std::ofstream output_signal_str( output_signal_filename.c_str( ) );
125 TFilter::TCurve curve = filter->GetCurve( );
126 unsigned long max_count = 0;
127 for( TFilter::TCurve::value_type d : curve )
129 output_signal_str << d.first << " " << d.second << std::endl;
130 if( max_count < d.second )
131 max_count = d.second;
134 output_signal_str.close( );
136 // Write gnuplot script
137 std::string output_gnuplot_filename = output_image_filename + ".gnuplot";
138 std::ofstream output_gnuplot_str( output_gnuplot_filename.c_str( ) );
139 unsigned int thr_pos = filter->GetOptimumThreshold( );
140 int thr = curve[ thr_pos ].first;
142 << "set term png" << std::endl
143 << "set output \"" << output_image_filename << ".png\"" << std::endl
144 << "set arrow 1 from " << thr
145 << ",0 to " << thr << "," << max_count
147 << "show arrow 1" << std::endl
148 << "plot \"" << output_signal_filename
149 << "\" using 1:2 with linespoints title \"Evolution ("
150 << thr << "," << thr_pos << ") " << "\""
152 output_gnuplot_str.close( );