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