segmentation/clitkLevelSetSegmentationGenericFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   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
   8
   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.
  12
  13   It is distributed under dual licence
  14
  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
  20
  21 /* =================================================
  22  * @file   clitkLevelSetSegmentationGenericFilter.txx
  23  * @author Jef Vandemeulebroucke <jef@creatis.insa-lyon.fr>
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30
  31 namespace clitk
  32 {
  33
  34   //-------------------------------------------------------------------
  35   // Update with the number of dimensions
  36   //-------------------------------------------------------------------
  37   template<unsigned int Dimension>
  38   void
  39   LevelSetSegmentationGenericFilter::UpdateWithDim(std::string PixelType)
  40   {
  41     if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
  42
  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>();
  46     //     }
  47
  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>();
  51     //     }
  52
  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>();
  56     //     }
  57
  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>();
  61     //     }
  62     //else {
  63     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
  64     UpdateWithDimAndPixelType<Dimension, float>();
  65     // }
  66   }
  67
  68
  69   //-------------------------------------------------------------------
  70   // Update with the number of dimensions and the pixeltype
  71   //-------------------------------------------------------------------
  72   template <unsigned int Dimension, class  PixelType>
  73   void
  74   LevelSetSegmentationGenericFilter::UpdateWithDimAndPixelType()
  75   {
  76
  77     // ImageTypes
  78     typedef itk::Image<float, Dimension> InputImageType;
  79     typedef itk::Image<float, Dimension> FeatureImageType;
  80     typedef itk::Image<float, Dimension> OutputImageType;
  81
  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();
  87
  88     // Feature image
  89     typename FeatureImageType::Pointer featureImage;
  90     if( m_ArgsInfo.feature_given)
  91       {
  92         // Read it
  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();
  98
  99         // Edge preserving smoothing
 100         if (m_ArgsInfo.smooth_flag)
 101           {
 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();
 110
 111           }
 112
 113         // Recursive gaussian gradient magnitude
 114         if (m_ArgsInfo.gradMag_flag)
 115           {
 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();
 121
 122           }
 123
 124         // Recursive gaussian gradient magnitude
 125         if (m_ArgsInfo.gradMagGauss_flag)
 126           {
 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();
 133           }
 134
 135         // Sigmoid
 136         if (m_ArgsInfo.sigmoid_flag)
 137           {
 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  );
 143             sigmoid->Update();
 144             featureImage=sigmoid->GetOutput();
 145           }
 146
 147       }
 148
 149     // Filter
 150     typename FeatureImageType::Pointer output;
 151     if (m_ArgsInfo.GAC_flag)
 152       {
 153
 154         // Create the filter
 155         typedef  itk::GeodesicActiveContourLevelSetImageFilter< InputImageType, FeatureImageType >    GeodesicActiveContourFilterType;
 156         typename GeodesicActiveContourFilterType::Pointer geodesicActiveContour = GeodesicActiveContourFilterType::New();
 157
 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);
 165
 166         // Monitor
 167         unsigned int totalNumberOfIterations=0;
 168         if(m_ArgsInfo.monitorIm_given)
 169           {
 170             geodesicActiveContour->SetNumberOfIterations( m_ArgsInfo.monitorIt_arg );
 171             while (true)
 172               {
 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;
 180               }
 181           }
 182         else
 183           {
 184             geodesicActiveContour->SetNumberOfIterations( m_ArgsInfo.iter_arg );
 185             geodesicActiveContour->Update();
 186             totalNumberOfIterations=geodesicActiveContour->GetElapsedIterations();
 187           }
 188
 189
 190         // Print
 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;
 197
 198         output = geodesicActiveContour->GetOutput();
 199       }
 200
 201     // Write levelset
 202     if (m_ArgsInfo.levelSet_given)
 203       {
 204         typedef itk::ImageFileWriter<FeatureImageType> WriterType;
 205         typename WriterType::Pointer writer = WriterType::New();
 206         writer->SetFileName(m_ArgsInfo.levelSet_arg);
 207         writer->SetInput(output);
 208         writer->Update();
 209       }
 210
 211     // Threshold
 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();
 221
 222     // Output
 223     typedef itk::ImageFileWriter<FeatureImageType> WriterType;
 224     typename WriterType::Pointer writer = WriterType::New();
 225     writer->SetFileName(m_ArgsInfo.output_arg);
 226     writer->SetInput(output);
 227     writer->Update();
 228
 229   }
 230
 231
 232 }//end clitk
 233
 234 #endif //#define clitkLevelSetSegmentationGenericFilter_txx