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 #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 template <class TInputImageType>
  85 template<class ArgsInfoType>
  86 void
  87 clitk::ExtractBonesFilter<TInputImageType>::
  88 SetArgsInfo(ArgsInfoType mArgsInfo)
  89 {
  90   SetVerboseOption_GGO(mArgsInfo);
  91   SetVerboseStep_GGO(mArgsInfo);
  92   SetWriteStep_GGO(mArgsInfo);
  93   SetVerboseWarningOff_GGO(mArgsInfo);
  94
  95   SetAFDBFilename_GGO(mArgsInfo);
  96   SetOutputBonesFilename_GGO(mArgsInfo);
  97
  98   SetInitialSmoothing_GGO(mArgsInfo);
  99   SetSmoothingConductanceParameter_GGO(mArgsInfo);
 100   SetSmoothingNumberOfIterations_GGO(mArgsInfo);
 101   SetSmoothingTimeStep_GGO(mArgsInfo);
 102   SetSmoothingUseImageSpacing_GGO(mArgsInfo);
 103
 104   SetMinimalComponentSize_GGO(mArgsInfo);
 105   SetUpperThreshold1_GGO(mArgsInfo);
 106   SetLowerThreshold1_GGO(mArgsInfo);
 107   SetFullConnectivity_GGO(mArgsInfo);
 108
 109   SetUpperThreshold2_GGO(mArgsInfo);
 110   SetLowerThreshold2_GGO(mArgsInfo);
 111   SetRadius2_GGO(mArgsInfo);
 112   SetSampleRate2_GGO(mArgsInfo);
 113   SetAutoCrop_GGO(mArgsInfo);
 114   SetFillHoles_GGO(mArgsInfo);
 115 }
 116 //--------------------------------------------------------------------
 117
 118
 119 //--------------------------------------------------------------------
 120 template <class TInputImageType>
 121 void
 122 clitk::ExtractBonesFilter<TInputImageType>::
 123 GenerateOutputInformation() {
 124
 125   // Get input pointers
 126   InputImagePointer input   = dynamic_cast<TInputImageType*>(itk::ProcessObject::GetInput(0));
 127   Superclass::GenerateOutputInformation();
 128   MaskImagePointer outputImage = this->GetOutput(0);
 129   outputImage->SetRegions(input->GetLargestPossibleRegion());
 130
 131   // Read DB
 132   LoadAFDB();
 133
 134   // typedefs
 135   typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
 136   typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> BinarizeFilterType;
 137   typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
 138   typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
 139   typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
 140   typedef itk::CastImageFilter<InternalImageType,MaskImageType> CastImageFilterType;
 141   typedef itk::ImageFileWriter<MaskImageType> WriterType;
 142
 143   //---------------------------------
 144   // Smoothing [Optional]
 145   //---------------------------------
 146   if (GetInitialSmoothing()) {
 147     StartNewStep("Initial Smoothing");
 148     typedef itk::CurvatureAnisotropicDiffusionImageFilter<InputImageType, InputImageType> FilterType;
 149     typename FilterType::Pointer df = FilterType::New();
 150     df->SetConductanceParameter(GetSmoothingConductanceParameter());
 151     df->SetNumberOfIterations(GetSmoothingNumberOfIterations());
 152     df->SetTimeStep(GetSmoothingTimeStep());
 153     df->SetUseImageSpacing(GetSmoothingUseImageSpacing());
 154     df->SetInput(input);
 155     df->Update();
 156     filtered_input = df->GetOutput();
 157     StopCurrentStep<InputImageType>(filtered_input);
 158   }
 159   else {
 160     filtered_input = input;
 161   }
 162
 163   //--------------------------------------------------------------------
 164   //--------------------------------------------------------------------
 165   StartNewStep("Initial Labeling");
 166   typename InternalImageType::Pointer firstLabelImage;
 167
 168   //---------------------------------
 169   // Binarize the image
 170   //---------------------------------
 171   typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 172   binarizeFilter->SetInput(filtered_input);
 173   binarizeFilter->SetLowerThreshold(GetLowerThreshold1());
 174   binarizeFilter->SetUpperThreshold(GetUpperThreshold1());
 175   binarizeFilter->SetInsideValue(this->GetForegroundValue());
 176   binarizeFilter->SetOutsideValue(this->GetBackgroundValue());
 177
 178   //---------------------------------
 179   // Label the connected components
 180   //---------------------------------
 181   typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
 182   connectFilter->SetInput(binarizeFilter->GetOutput());
 183   connectFilter->SetBackgroundValue(this->GetBackgroundValue());
 184   connectFilter->SetFullyConnected(GetFullConnectivity());
 185
 186   //---------------------------------
 187   // Sort the labels according to size
 188   //---------------------------------
 189   typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
 190   relabelFilter->SetInput(connectFilter->GetOutput());
 191   relabelFilter->SetMinimumObjectSize(GetMinimalComponentSize());
 192
 193   //---------------------------------
 194   // Keep the label
 195   //---------------------------------
 196   typename BinarizeFilterType::Pointer binarizeFilter2=BinarizeFilterType::New();
 197   binarizeFilter2->SetInput(relabelFilter->GetOutput());
 198   binarizeFilter2->SetLowerThreshold(1);
 199   binarizeFilter2->SetUpperThreshold(1);
 200   binarizeFilter2->SetInsideValue(this->GetForegroundValue());
 201   binarizeFilter2->SetOutsideValue(this->GetBackgroundValue());
 202   binarizeFilter2->Update();
 203
 204   firstLabelImage = binarizeFilter2->GetOutput();
 205   StopCurrentStep<InternalImageType>(firstLabelImage);
 206
 207   //--------------------------------------------------------------------
 208   //--------------------------------------------------------------------
 209   StartNewStep("Neighborhood connected filter");
 210   typename InternalImageType::Pointer secondLabelImage;
 211
 212   //---------------------------------
 213   //Neighborhood connected RG
 214   //---------------------------------
 215   typedef itk::NeighborhoodConnectedImageFilter<InputImageType, InternalImageType>
 216     NeighborhoodConnectedImageFilterType;
 217   typename NeighborhoodConnectedImageFilterType::Pointer neighborhoodConnectedImageFilter=
 218     NeighborhoodConnectedImageFilterType::New();
 219
 220   // thresholds
 221   neighborhoodConnectedImageFilter->SetLower(GetLowerThreshold2());
 222   neighborhoodConnectedImageFilter->SetUpper(GetUpperThreshold2());
 223   neighborhoodConnectedImageFilter->SetReplaceValue(this->GetForegroundValue());
 224   neighborhoodConnectedImageFilter->SetRadius(GetRadius2());
 225   neighborhoodConnectedImageFilter->SetInput(filtered_input);
 226
 227   // Seeds from label image
 228   typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
 229   IteratorType it(firstLabelImage, firstLabelImage->GetLargestPossibleRegion());
 230   typename InputImageType::IndexType index;
 231   unsigned int counter=0;
 232   while (!it.IsAtEnd())
 233     {
 234       if (it.Get()==this->GetForegroundValue())
 235         {
 236           counter++;
 237           index=it.GetIndex();
 238           neighborhoodConnectedImageFilter->AddSeed(index);
 239           ++it;
 240           unsigned int i=0;
 241           while (!it.IsAtEnd()  &&  i< (unsigned int) GetSampleRate2())
 242             {
 243               ++it;
 244               i++;
 245             }
 246         }
 247       else ++it;
 248     }
 249
 250   neighborhoodConnectedImageFilter->Update();
 251   secondLabelImage = neighborhoodConnectedImageFilter->GetOutput();
 252
 253   //--------------------------------------------------------------------
 254   //--------------------------------------------------------------------
 255   StartNewStep("Combine the images");
 256   typedef clitk::SetBackgroundImageFilter<InternalImageType, InternalImageType, InternalImageType>
 257     SetBackgroundImageFilterType;
 258   typename SetBackgroundImageFilterType::Pointer setBackgroundFilter=SetBackgroundImageFilterType::New();
 259   setBackgroundFilter->SetInput(firstLabelImage);
 260   setBackgroundFilter->SetInput2(secondLabelImage);
 261   setBackgroundFilter->SetMaskValue(this->GetForegroundValue());
 262   setBackgroundFilter->SetOutsideValue(this->GetForegroundValue());
 263   setBackgroundFilter->Update();
 264
 265   output = setBackgroundFilter->GetOutput();
 266
 267   //--------------------------------------------------------------------
 268   //--------------------------------------------------------------------
 269   // Fill Bones
 270   if (GetFillHoles()) {
 271     StartNewStep("Fill Holes");
 272     typedef clitk::FillMaskFilter<InternalImageType> FillMaskFilterType;
 273     typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New();
 274     fillMaskFilter->SetInput(output);
 275     fillMaskFilter->AddDirection(2);
 276     fillMaskFilter->Update();
 277     output = fillMaskFilter->GetOutput();
 278     StopCurrentStep<InternalImageType>(output);
 279   }
 280
 281   //--------------------------------------------------------------------
 282   //--------------------------------------------------------------------
 283   // [Optional]
 284   if (GetAutoCrop()) {
 285     StartNewStep("AutoCrop");
 286     typedef clitk::AutoCropFilter<InternalImageType> CropFilterType;
 287     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 288     cropFilter->SetInput(output);
 289     cropFilter->SetBackgroundValue(GetBackgroundValue());
 290     cropFilter->Update();
 291     output = cropFilter->GetOutput();
 292     StopCurrentStep<InternalImageType>(output);
 293     outputImage->SetRegions(output->GetLargestPossibleRegion());
 294   }
 295
 296 }
 297 //--------------------------------------------------------------------
 298
 299
 300 //--------------------------------------------------------------------
 301 template <class TInputImageType>
 302 void
 303 clitk::ExtractBonesFilter<TInputImageType>::
 304 GenerateData() {
 305
 306   //--------------------------------------------------------------------
 307   //--------------------------------------------------------------------
 308   // Final Cast
 309   typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 310   typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 311   caster->SetInput(output);
 312   caster->Update();
 313   this->GraftOutput(caster->GetOutput());
 314
 315   // Store image filenames into AFDB
 316   GetAFDB()->SetImageFilename("bones", this->GetOutputBonesFilename());
 317   WriteAFDB();
 318   return;
 319 }
 320 //--------------------------------------------------------------------
 321
 322
 323 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX