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