segmentation/clitkExtractBonesFilter.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
  19 #ifndef CLITKEXTRACTBONESSFILTER_TXX
  20 #define CLITKEXTRACTBONESSFILTER_TXX
  21
  22 // clitk
  23 #include "clitkImageCommon.h"
  24 #include "clitkSetBackgroundImageFilter.h"
  25 #include "clitkSegmentationUtils.h"
  26 #include "clitkAutoCropFilter.h"
  27
  28 // itk
  29 #include "itkBinaryThresholdImageFilter.h"
  30 #include "itkConnectedComponentImageFilter.h"
  31 #include "itkRelabelComponentImageFilter.h"
  32 #include "itkNeighborhoodConnectedImageFilter.h"
  33 #include "itkCurvatureAnisotropicDiffusionImageFilter.h"
  34
  35 //--------------------------------------------------------------------
  36 template <class TInputImageType>
  37 clitk::ExtractBonesFilter<TInputImageType>::
  38 ExtractBonesFilter():
  39   clitk::FilterBase(),
  40   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  41   itk::ImageToImageFilter<TInputImageType, MaskImageType>()
  42 {
  43   // Default global options
  44   this->SetNumberOfRequiredInputs(1);
  45   SetBackgroundValue(0); // Must be zero
  46   SetForegroundValue(1);
  47   SetInitialSmoothing(false);
  48
  49   SetSmoothingConductanceParameter(3.0);
  50   SetSmoothingNumberOfIterations(5);
  51   SetSmoothingTimeStep(0.0625);
  52   SetSmoothingUseImageSpacing(false);
  53
  54   SetMinimalComponentSize(100);
  55   SetUpperThreshold1(1500);
  56   SetLowerThreshold1(100);
  57   SetFullConnectivity(false);
  58
  59   SetUpperThreshold2(1500);
  60   SetLowerThreshold2(10);
  61   InputImageSizeType s;
  62   s.Fill(1);
  63   SetRadius2(s);
  64   SetSampleRate2(0);
  65   AutoCropOn();
  66 }
  67 //--------------------------------------------------------------------
  68
  69
  70 //--------------------------------------------------------------------
  71 template <class TInputImageType>
  72 void
  73 clitk::ExtractBonesFilter<TInputImageType>::
  74 SetInput(const TInputImageType * image)
  75 {
  76   this->SetNthInput(0, const_cast<TInputImageType *>(image));
  77 }
  78 //--------------------------------------------------------------------
  79
  80
  81 //--------------------------------------------------------------------
  82 template <class TInputImageType>
  83 template<class ArgsInfoType>
  84 void
  85 clitk::ExtractBonesFilter<TInputImageType>::
  86 SetArgsInfo(ArgsInfoType mArgsInfo)
  87 {
  88   SetVerboseOption_GGO(mArgsInfo);
  89   SetVerboseStep_GGO(mArgsInfo);
  90   SetWriteStep_GGO(mArgsInfo);
  91   SetVerboseWarningOff_GGO(mArgsInfo);
  92
  93   SetOutputBonesFilename_GGO(mArgsInfo);
  94
  95   SetInitialSmoothing_GGO(mArgsInfo);
  96   SetSmoothingConductanceParameter_GGO(mArgsInfo);
  97   SetSmoothingNumberOfIterations_GGO(mArgsInfo);
  98   SetSmoothingTimeStep_GGO(mArgsInfo);
  99   SetSmoothingUseImageSpacing_GGO(mArgsInfo);
 100
 101   SetMinimalComponentSize_GGO(mArgsInfo);
 102   SetUpperThreshold1_GGO(mArgsInfo);
 103   SetLowerThreshold1_GGO(mArgsInfo);
 104   SetFullConnectivity_GGO(mArgsInfo);
 105
 106   SetUpperThreshold2_GGO(mArgsInfo);
 107   SetLowerThreshold2_GGO(mArgsInfo);
 108   SetRadius2_GGO(mArgsInfo);
 109   SetSampleRate2_GGO(mArgsInfo);
 110   SetAutoCrop_GGO(mArgsInfo);
 111
 112   SetAFDBFilename_GGO(mArgsInfo);
 113 }
 114 //--------------------------------------------------------------------
 115
 116
 117 //--------------------------------------------------------------------
 118 template <class TInputImageType>
 119 void
 120 clitk::ExtractBonesFilter<TInputImageType>::
 121 GenerateOutputInformation() {
 122
 123   // Get input pointers
 124   InputImagePointer input   = dynamic_cast<TInputImageType*>(itk::ProcessObject::GetInput(0));
 125   Superclass::GenerateOutputInformation();
 126   MaskImagePointer outputImage = this->GetOutput(0);
 127   outputImage->SetRegions(input->GetLargestPossibleRegion());
 128
 129   // Read DB
 130   LoadAFDB();
 131
 132   // typedefs
 133   typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
 134   typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> BinarizeFilterType;
 135   typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
 136   typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
 137   typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
 138   typedef itk::CastImageFilter<InternalImageType,MaskImageType> CastImageFilterType;
 139   typedef itk::ImageFileWriter<MaskImageType> WriterType;
 140
 141   //---------------------------------
 142   // Smoothing [Optional]
 143   //---------------------------------
 144   if (GetInitialSmoothing()) {
 145     StartNewStep("Initial Smoothing");
 146     typedef itk::CurvatureAnisotropicDiffusionImageFilter<InputImageType, InputImageType> FilterType;
 147     typename FilterType::Pointer df = FilterType::New();
 148     df->SetConductanceParameter(GetSmoothingConductanceParameter());
 149     df->SetNumberOfIterations(GetSmoothingNumberOfIterations());
 150     df->SetTimeStep(GetSmoothingTimeStep());
 151     df->SetUseImageSpacing(GetSmoothingUseImageSpacing());
 152     df->SetInput(input);
 153     df->Update();
 154     filtered_input = df->GetOutput();
 155     StopCurrentStep<InputImageType>(filtered_input);
 156   }
 157   else {
 158     filtered_input = input;
 159   }
 160
 161   //--------------------------------------------------------------------
 162   //--------------------------------------------------------------------
 163   StartNewStep("Initial Labeling");
 164   typename InternalImageType::Pointer firstLabelImage;
 165
 166   //---------------------------------
 167   // Binarize the image
 168   //---------------------------------
 169   typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 170   binarizeFilter->SetInput(filtered_input);
 171   binarizeFilter->SetLowerThreshold(GetLowerThreshold1());
 172   binarizeFilter->SetUpperThreshold(GetUpperThreshold1());
 173   binarizeFilter->SetInsideValue(this->GetForegroundValue());
 174   binarizeFilter->SetOutsideValue(this->GetBackgroundValue());
 175
 176   //---------------------------------
 177   // Label the connected components
 178   //---------------------------------
 179   typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
 180   connectFilter->SetInput(binarizeFilter->GetOutput());
 181   connectFilter->SetBackgroundValue(this->GetBackgroundValue());
 182   connectFilter->SetFullyConnected(GetFullConnectivity());
 183
 184   //---------------------------------
 185   // Sort the labels according to size
 186   //---------------------------------
 187   typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
 188   relabelFilter->SetInput(connectFilter->GetOutput());
 189   relabelFilter->SetMinimumObjectSize(GetMinimalComponentSize());
 190
 191   //---------------------------------
 192   // Keep the label
 193   //---------------------------------
 194   typename BinarizeFilterType::Pointer binarizeFilter2=BinarizeFilterType::New();
 195   binarizeFilter2->SetInput(relabelFilter->GetOutput());
 196   binarizeFilter2->SetLowerThreshold(1);
 197   binarizeFilter2->SetUpperThreshold(1);
 198   binarizeFilter2->SetInsideValue(this->GetForegroundValue());
 199   binarizeFilter2->SetOutsideValue(this->GetBackgroundValue());
 200   binarizeFilter2->Update();
 201
 202   firstLabelImage = binarizeFilter2->GetOutput();
 203   StopCurrentStep<InternalImageType>(firstLabelImage);
 204
 205   //--------------------------------------------------------------------
 206   //--------------------------------------------------------------------
 207   StartNewStep("Neighborhood connected filter");
 208   typename InternalImageType::Pointer secondLabelImage;
 209
 210   //---------------------------------
 211   //Neighborhood connected RG
 212   //---------------------------------
 213   typedef itk::NeighborhoodConnectedImageFilter<InputImageType, InternalImageType>
 214     NeighborhoodConnectedImageFilterType;
 215   typename NeighborhoodConnectedImageFilterType::Pointer neighborhoodConnectedImageFilter=
 216     NeighborhoodConnectedImageFilterType::New();
 217
 218   // thresholds
 219   neighborhoodConnectedImageFilter->SetLower(GetLowerThreshold2());
 220   neighborhoodConnectedImageFilter->SetUpper(GetUpperThreshold2());
 221   neighborhoodConnectedImageFilter->SetReplaceValue(this->GetForegroundValue());
 222   neighborhoodConnectedImageFilter->SetRadius(GetRadius2());
 223   neighborhoodConnectedImageFilter->SetInput(filtered_input);
 224
 225   // Seeds from label image
 226   typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
 227   IteratorType it(firstLabelImage, firstLabelImage->GetLargestPossibleRegion());
 228   typename InputImageType::IndexType index;
 229   unsigned int counter=0;
 230   while (!it.IsAtEnd())
 231     {
 232       if (it.Get()==this->GetForegroundValue())
 233         {
 234           counter++;
 235           index=it.GetIndex();
 236           neighborhoodConnectedImageFilter->AddSeed(index);
 237           ++it;
 238           unsigned int i=0;
 239           while (!it.IsAtEnd()  &&  i< (unsigned int) GetSampleRate2())
 240             {
 241               ++it;
 242               i++;
 243             }
 244         }
 245       else ++it;
 246     }
 247
 248   neighborhoodConnectedImageFilter->Update();
 249   secondLabelImage = neighborhoodConnectedImageFilter->GetOutput();
 250
 251   //--------------------------------------------------------------------
 252   //--------------------------------------------------------------------
 253   StartNewStep("Combine the images");
 254   typedef clitk::SetBackgroundImageFilter<InternalImageType, InternalImageType, InternalImageType>
 255     SetBackgroundImageFilterType;
 256   typename SetBackgroundImageFilterType::Pointer setBackgroundFilter=SetBackgroundImageFilterType::New();
 257   setBackgroundFilter->SetInput(firstLabelImage);
 258   setBackgroundFilter->SetInput2(secondLabelImage);
 259   setBackgroundFilter->SetMaskValue(this->GetForegroundValue());
 260   setBackgroundFilter->SetOutsideValue(this->GetForegroundValue());
 261   setBackgroundFilter->Update();
 262
 263   output = setBackgroundFilter->GetOutput();
 264
 265   //--------------------------------------------------------------------
 266   //--------------------------------------------------------------------
 267   // [Optional]
 268   if (GetAutoCrop()) {
 269     StartNewStep("AutoCrop");
 270     typedef clitk::AutoCropFilter<InternalImageType> CropFilterType;
 271     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 272     cropFilter->SetInput(output);
 273     cropFilter->SetBackgroundValue(GetBackgroundValue());
 274     cropFilter->Update();
 275     output = cropFilter->GetOutput();
 276     StopCurrentStep<InternalImageType>(output);
 277     outputImage->SetRegions(output->GetLargestPossibleRegion());
 278   }
 279
 280 }
 281 //--------------------------------------------------------------------
 282
 283
 284 //--------------------------------------------------------------------
 285 template <class TInputImageType>
 286 void
 287 clitk::ExtractBonesFilter<TInputImageType>::
 288 GenerateData() {
 289
 290   //--------------------------------------------------------------------
 291   //--------------------------------------------------------------------
 292   // Final Cast
 293   typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 294   typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 295   caster->SetInput(output);
 296   caster->Update();
 297   this->GraftOutput(caster->GetOutput());
 298
 299   // Store image filenames into AFDB
 300   GetAFDB()->SetImageFilename("bones", this->GetOutputBonesFilename());
 301   WriteAFDB();
 302   return;
 303 }
 304 //--------------------------------------------------------------------
 305
 306
 307 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX