1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
18 #ifndef clitkLevelSetSegmentationGenericFilter_txx
19 #define clitkLevelSetSegmentationGenericFilter_txx
21 /* =================================================
22 * @file clitkLevelSetSegmentationGenericFilter.txx
23 * @author Jef Vandemeulebroucke <jef@creatis.insa-lyon.fr>
28 ===================================================*/
34 //-------------------------------------------------------------------
35 // Update with the number of dimensions
36 //-------------------------------------------------------------------
37 template<unsigned int Dimension>
39 LevelSetSegmentationGenericFilter::UpdateWithDim(std::string PixelType)
41 if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
43 // if(PixelType == "short"){
44 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
45 // UpdateWithDimAndPixelType<Dimension, signed short>();
48 // else if(PixelType == "unsigned_short"){
49 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
50 // UpdateWithDimAndPixelType<Dimension, unsigned short>();
53 // else if (PixelType == "unsigned_char"){
54 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
55 // UpdateWithDimAndPixelType<Dimension, unsigned char>();
58 // else if (PixelType == "char"){
59 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
60 // UpdateWithDimAndPixelType<Dimension, signed char>();
63 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
64 UpdateWithDimAndPixelType<Dimension, float>();
69 //-------------------------------------------------------------------
70 // Update with the number of dimensions and the pixeltype
71 //-------------------------------------------------------------------
72 template <unsigned int Dimension, class PixelType>
74 LevelSetSegmentationGenericFilter::UpdateWithDimAndPixelType()
78 typedef itk::Image<float, Dimension> InputImageType;
79 typedef itk::Image<float, Dimension> FeatureImageType;
80 typedef itk::Image<float, Dimension> OutputImageType;
82 // Read the input (initial level set)
83 typedef itk::ImageFileReader<InputImageType> InputReaderType;
84 typename InputReaderType::Pointer reader = InputReaderType::New();
85 reader->SetFileName( m_ArgsInfo.input_arg);
86 typename InputImageType::Pointer input= reader->GetOutput();
89 typename FeatureImageType::Pointer featureImage;
90 if( m_ArgsInfo.feature_given)
93 typedef itk::ImageFileReader<FeatureImageType> FeatureReaderType;
94 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
95 featureReader->SetFileName(m_ArgsInfo.feature_arg);
96 featureReader->Update();
97 featureImage=featureReader->GetOutput();
99 // Edge preserving smoothing
100 if (m_ArgsInfo.smooth_flag)
102 typedef itk::CurvatureAnisotropicDiffusionImageFilter<FeatureImageType, FeatureImageType > SmoothingFilterType;
103 typename SmoothingFilterType::Pointer smoothingFilter = SmoothingFilterType::New();
104 smoothingFilter->SetInput(featureImage);
105 smoothingFilter->SetTimeStep( m_ArgsInfo.timeStep_arg );
106 smoothingFilter->SetNumberOfIterations( m_ArgsInfo.iterSmooth_arg );
107 smoothingFilter->SetConductanceParameter( m_ArgsInfo.cond_arg );
108 smoothingFilter->Update();
109 featureImage=smoothingFilter->GetOutput();
113 // Recursive gaussian gradient magnitude
114 if (m_ArgsInfo.gradMag_flag)
116 typedef itk::GradientMagnitudeImageFilter< FeatureImageType,FeatureImageType> GradientFilterType;
117 typename GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();
118 gradientMagnitude->SetInput(featureImage);
119 gradientMagnitude->Update();
120 featureImage=gradientMagnitude->GetOutput();
124 // Recursive gaussian gradient magnitude
125 if (m_ArgsInfo.gradMagGauss_flag)
127 typedef itk::GradientMagnitudeRecursiveGaussianImageFilter< FeatureImageType,FeatureImageType> GradientFilterType;
128 typename GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();
129 gradientMagnitude->SetInput(featureImage);
130 gradientMagnitude->SetSigma( m_ArgsInfo.sigma_arg );
131 gradientMagnitude->Update();
132 featureImage=gradientMagnitude->GetOutput();
136 if (m_ArgsInfo.sigmoid_flag)
138 typedef itk::SigmoidImageFilter<FeatureImageType, FeatureImageType > SigmoidFilterType;
139 typename SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
140 sigmoid->SetInput(featureImage);
141 sigmoid->SetAlpha( m_ArgsInfo.alpha_arg );
142 sigmoid->SetBeta( m_ArgsInfo.beta_arg );
144 featureImage=sigmoid->GetOutput();
150 typename FeatureImageType::Pointer output;
151 if (m_ArgsInfo.GAC_flag)
155 typedef itk::GeodesicActiveContourLevelSetImageFilter< InputImageType, FeatureImageType > GeodesicActiveContourFilterType;
156 typename GeodesicActiveContourFilterType::Pointer geodesicActiveContour = GeodesicActiveContourFilterType::New();
158 geodesicActiveContour->SetPropagationScaling( m_ArgsInfo.propScale_arg );
159 geodesicActiveContour->SetCurvatureScaling( m_ArgsInfo.curveScale_arg );
160 geodesicActiveContour->SetAdvectionScaling( m_ArgsInfo.advectionScale_arg );
161 geodesicActiveContour->SetMaximumRMSError( m_ArgsInfo.maxRMS_arg );
162 geodesicActiveContour->SetInput( input );
163 geodesicActiveContour->SetFeatureImage( featureImage );
164 geodesicActiveContour->SetUseImageSpacing(true);
167 unsigned int totalNumberOfIterations=0;
168 if(m_ArgsInfo.monitorIm_given)
170 geodesicActiveContour->SetNumberOfIterations( m_ArgsInfo.monitorIt_arg );
173 geodesicActiveContour->Update();
174 totalNumberOfIterations+=geodesicActiveContour->GetElapsedIterations();
175 if(m_Verbose) std::cout <<"Writing image after "<< totalNumberOfIterations<<"..."<<std::endl;
176 writeImage<InputImageType>(geodesicActiveContour->GetOutput(), m_ArgsInfo.monitorIm_arg);
177 geodesicActiveContour->SetInput(geodesicActiveContour->GetOutput());
178 geodesicActiveContour->SetNumberOfIterations( std::min( (m_ArgsInfo.iter_arg-totalNumberOfIterations) ,(unsigned int) m_ArgsInfo.monitorIt_arg ) );
179 if (totalNumberOfIterations> (unsigned int) m_ArgsInfo.iter_arg) break;
184 geodesicActiveContour->SetNumberOfIterations( m_ArgsInfo.iter_arg );
185 geodesicActiveContour->Update();
186 totalNumberOfIterations=geodesicActiveContour->GetElapsedIterations();
191 std::cout << std::endl;
192 std::cout << "Max. no. iterations: " << m_ArgsInfo.iter_arg << std::endl;
193 std::cout << "Max. RMS error: " << geodesicActiveContour->GetMaximumRMSError() << std::endl;
194 std::cout << std::endl;
195 std::cout << "No. elpased iterations: " << totalNumberOfIterations << std::endl;
196 std::cout << "RMS change: " << geodesicActiveContour->GetRMSChange() << std::endl;
198 output = geodesicActiveContour->GetOutput();
202 if (m_ArgsInfo.levelSet_given)
204 typedef itk::ImageFileWriter<FeatureImageType> WriterType;
205 typename WriterType::Pointer writer = WriterType::New();
206 writer->SetFileName(m_ArgsInfo.levelSet_arg);
207 writer->SetInput(output);
212 typedef itk::BinaryThresholdImageFilter< FeatureImageType,FeatureImageType > ThresholdingFilterType;
213 typename ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
214 thresholder->SetLowerThreshold( -1000.0 );
215 thresholder->SetUpperThreshold( 0.0 );
216 thresholder->SetOutsideValue( 0 );
217 thresholder->SetInsideValue( 1 );
218 thresholder->SetInput( output );
219 thresholder->Update();
220 output=thresholder->GetOutput();
223 typedef itk::ImageFileWriter<FeatureImageType> WriterType;
224 typename WriterType::Pointer writer = WriterType::New();
225 writer->SetFileName(m_ArgsInfo.output_arg);
226 writer->SetInput(output);
234 #endif //#define clitkLevelSetSegmentationGenericFilter_txx