segmentation/clitkExtractPatientFilter.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 CLITKEXTRACTPATIENTFILTER_TXX
  20 #define CLITKEXTRACTPATIENTFILTER_TXX
  21
  22 // clitk
  23 #include "clitkImageCommon.h"
  24 #include "clitkSetBackgroundImageFilter.h"
  25 #include "clitkDecomposeAndReconstructImageFilter.h"
  26 #include "clitkAutoCropFilter.h"
  27 #include "clitkMemoryUsage.h"
  28 #include "clitkSegmentationUtils.h"
  29
  30 // itk
  31 #include "itkBinaryThresholdImageFilter.h"
  32 #include "itkConnectedComponentImageFilter.h"
  33 #include "itkRelabelComponentImageFilter.h"
  34 #include "itkBinaryMorphologicalClosingImageFilter.h"
  35 #include "itkBinaryMorphologicalOpeningImageFilter.h"
  36 #include "itkBinaryBallStructuringElement.h"
  37 #include "itkCastImageFilter.h"
  38 #include "itkConstantPadImageFilter.h"
  39
  40 //--------------------------------------------------------------------
  41 template <class TInputImageType>
  42 clitk::ExtractPatientFilter<TInputImageType>::
  43 ExtractPatientFilter():
  44   clitk::FilterBase(),
  45   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  46   itk::ImageToImageFilter<TInputImageType, MaskImageType>()
  47 {
  48   this->SetNumberOfRequiredInputs(1);
  49   SetBackgroundValue(0); // Must be zero
  50   SetForegroundValue(1);
  51   SetPrimaryOpeningRadius(0);
  52
  53   // Step 1: Threshold + CC + sort (Find low density areas)
  54   SetUpperThreshold(-300);
  55   SetLowerThreshold(-1000);
  56   UseLowerThresholdOff();
  57
  58   // Step 2: DecomposeAndReconstructImageFilter (optional)
  59   DecomposeAndReconstructDuringFirstStepOff();
  60   InternalImageSizeType r;
  61   r.Fill(1);
  62   SetRadius1(r);
  63   SetMaximumNumberOfLabels1(2);
  64   SetNumberOfNewLabels1(1);
  65
  66   // Step 3: Remove the air (largest area).
  67
  68   // Step 4: 2nd DecomposeAndReconstructImageFilter
  69   DecomposeAndReconstructDuringSecondStepOff();
  70   SetRadius2(r);
  71   SetMaximumNumberOfLabels2(2);
  72   SetNumberOfNewLabels2(1);
  73
  74   // Step 5: Only keep label corresponding (Keep patient's labels)
  75   SetFirstKeep(1);
  76   SetLastKeep(1);
  77
  78   // Step 4: OpenClose (option)
  79   FinalOpenCloseOff();
  80   AutoCropOn();
  81 }
  82 //--------------------------------------------------------------------
  83
  84
  85 //--------------------------------------------------------------------
  86 template <class TInputImageType>
  87 void
  88 clitk::ExtractPatientFilter<TInputImageType>::
  89 SetInput(const TInputImageType * image)
  90 {
  91   this->SetNthInput(0, const_cast<TInputImageType *>(image));
  92 }
  93 //--------------------------------------------------------------------
  94
  95
  96 //--------------------------------------------------------------------
  97 template <class TInputImageType>
  98 void
  99 clitk::ExtractPatientFilter<TInputImageType>::
 100 GenerateOutputInformation() {
 101
 102   clitk::PrintMemory(GetVerboseMemoryFlag(), "Initial memory"); // OK
 103
 104   Superclass::GenerateOutputInformation();
 105   input = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
 106
 107   // MaskImagePointer outputImage = this->GetOutput(0);
 108 //   outputImage->SetRegions(input->GetLargestPossibleRegion());
 109
 110   // Get input pointers
 111   static const unsigned int Dim = InputImageType::ImageDimension;
 112   //input = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
 113
 114   //--------------------------------------------------------------------
 115   //--------------------------------------------------------------------
 116   // Step 1:
 117   StartNewStep("Find low densities areas");
 118
 119   // Pad images with air to prevent patient touching the image border
 120   typedef itk::ConstantPadImageFilter<InputImageType, InputImageType> PadFilterType;
 121   typename PadFilterType::Pointer padFilter = PadFilterType::New();
 122   padFilter->SetInput(input);
 123   padFilter->SetConstant(GetUpperThreshold() - 1);
 124   typename InputImageType::SizeType bounds;
 125   for (unsigned i = 0; i < Dim - 1; ++i)
 126     bounds[i] = 1;
 127   bounds[Dim - 1] = 0;
 128   padFilter->SetPadLowerBound(bounds);
 129   padFilter->SetPadUpperBound(bounds);
 130
 131   typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> BinarizeFilterType;
 132   typename BinarizeFilterType::Pointer binarizeFilter=BinarizeFilterType::New();
 133   binarizeFilter->SetInput(padFilter->GetOutput());
 134   if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(GetLowerThreshold());
 135   binarizeFilter->SetUpperThreshold(GetUpperThreshold());
 136   binarizeFilter ->SetInsideValue(this->GetForegroundValue());
 137   binarizeFilter ->SetOutsideValue(this->GetBackgroundValue());
 138   working_image = binarizeFilter->GetOutput();
 139
 140   typedef itk::BinaryBallStructuringElement<InternalPixelType,Dim> KernelType;
 141   unsigned int radius = this->GetPrimaryOpeningRadius();
 142   if (radius > 0)
 143   {
 144     if (this->GetVerboseOptionFlag()) std::cout << ("Opening after threshold; R = ") << radius << std::endl;
 145     KernelType kernel;
 146     kernel.SetRadius(radius);
 147
 148     typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType2;
 149     typename OpenFilterType2::Pointer openFilter2 = OpenFilterType2::New();
 150     openFilter2->SetInput(working_image);
 151     openFilter2->SetBackgroundValue(0);
 152     openFilter2->SetForegroundValue(1);
 153     openFilter2->SetKernel(kernel);
 154     openFilter2->Update();
 155     working_image = openFilter2->GetOutput();
 156   }
 157
 158   if (this->GetVerboseOptionFlag()) std::cout << ("Labelling") << std::endl;
 159   // Connected component labeling
 160   typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
 161   typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
 162   connectFilter->SetInput(working_image);
 163   connectFilter->SetBackgroundValue(this->GetBackgroundValue());
 164   connectFilter->SetFullyConnected(false);
 165   connectFilter->Update();
 166
 167   if (this->GetVerboseOptionFlag()) std::cout << ("RelabelComponentImageFilter") << std::endl;
 168   // Sort labels according to size
 169   typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
 170   typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
 171   relabelFilter->InPlaceOn();
 172   relabelFilter->SetInput(connectFilter->GetOutput());
 173   relabelFilter->Update();
 174   working_image = relabelFilter->GetOutput();
 175
 176   // End
 177   StopCurrentStep<InternalImageType>(working_image);
 178
 179   //--------------------------------------------------------------------
 180   //--------------------------------------------------------------------
 181   // [Optional]
 182   if (GetDecomposeAndReconstructDuringFirstStep()) {
 183     StartNewStep("First Decompose & Reconstruct step");
 184     typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> FilterType;
 185     typename FilterType::Pointer f = FilterType::New();
 186     f->SetInput(working_image);
 187     // f->SetVerbose(m_Verbose);
 188     f->SetRadius(GetRadius1());
 189     f->SetMaximumNumberOfLabels(GetMaximumNumberOfLabels1());
 190     f->SetBackgroundValue(this->GetBackgroundValue());
 191     f->SetForegroundValue(this->GetForegroundValue());
 192     f->SetFullyConnected(true);
 193     f->SetNumberOfNewLabels(GetNumberOfNewLabels1());
 194     f->Update();
 195     working_image = f->GetOutput();
 196     StopCurrentStep<InternalImageType>(working_image);
 197   }
 198
 199   //--------------------------------------------------------------------
 200   //--------------------------------------------------------------------
 201   if (this->GetVerboseOptionFlag()) std::cout << ("Remove the air (largest area)") << std::endl;
 202   StartNewStep("Remove the air (largest area)");
 203   typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> iBinarizeFilterType;
 204   typename iBinarizeFilterType::Pointer binarizeFilter2 = iBinarizeFilterType::New();
 205   binarizeFilter2->SetInput(working_image);
 206   binarizeFilter2->SetLowerThreshold(GetFirstKeep());
 207   binarizeFilter2->SetUpperThreshold(GetLastKeep());
 208   binarizeFilter2 ->SetInsideValue(0);
 209   binarizeFilter2 ->SetOutsideValue(1);
 210   //  binarizeFilter2 ->Update(); // NEEDED ?
 211
 212   typename ConnectFilterType::Pointer connectFilter2 = ConnectFilterType::New();
 213   connectFilter2->SetInput(binarizeFilter2->GetOutput());
 214   connectFilter2->SetBackgroundValue(this->GetBackgroundValue());
 215   connectFilter2->SetFullyConnected(false);
 216
 217   typename RelabelFilterType::Pointer relabelFilter2 = RelabelFilterType::New();
 218   relabelFilter2->SetInput(connectFilter2->GetOutput());
 219   relabelFilter2->Update();
 220   working_image = relabelFilter2->GetOutput();
 221
 222   // Keep main label
 223   working_image = KeepLabels<InternalImageType>
 224     (working_image, GetBackgroundValue(), GetForegroundValue(), 1, 1, true);
 225   StopCurrentStep<InternalImageType>(working_image);
 226
 227   //--------------------------------------------------------------------
 228   //--------------------------------------------------------------------
 229   // [Optional]
 230   if (GetDecomposeAndReconstructDuringSecondStep()) {
 231     StartNewStep("Second Decompose & Reconstruct step");
 232     typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> FilterType;
 233     typename FilterType::Pointer f = FilterType::New();
 234     f->SetInput(working_image);
 235     // f->SetVerbose(m_Verbose);
 236     f->SetRadius(GetRadius2());
 237     f->SetMaximumNumberOfLabels(GetMaximumNumberOfLabels2());
 238     f->SetBackgroundValue(this->GetBackgroundValue());
 239     f->SetForegroundValue(this->GetForegroundValue());
 240     f->SetFullyConnected(true);
 241     f->SetNumberOfNewLabels(GetNumberOfNewLabels2());
 242     f->Update();
 243     working_image = f->GetOutput();
 244     StopCurrentStep<InternalImageType>(working_image);
 245   }
 246
 247   //--------------------------------------------------------------------
 248   //--------------------------------------------------------------------
 249   // [Optional]
 250   if (GetFinalOpenClose()) {
 251     StartNewStep("Final OpenClose");
 252     // Open
 253     KernelType structuringElement;
 254     structuringElement.SetRadius(1);
 255     structuringElement.CreateStructuringElement();
 256     typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType;
 257     typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
 258     openFilter->SetInput(working_image);
 259     openFilter->SetBackgroundValue(this->GetBackgroundValue());
 260     openFilter->SetForegroundValue(this->GetForegroundValue());
 261     openFilter->SetKernel(structuringElement);
 262     // Close
 263     typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType , KernelType> CloseFilterType;
 264     typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
 265     closeFilter->SetInput(openFilter->GetOutput());
 266     closeFilter->SetSafeBorder(true);
 267     closeFilter->SetForegroundValue(this->GetForegroundValue());
 268     //  closeFilter->SetBackgroundValue(SetBackgroundValue());
 269     closeFilter->SetKernel(structuringElement);
 270     closeFilter->Update();
 271     working_image = closeFilter->GetOutput();
 272     StopCurrentStep<InternalImageType>(working_image);
 273   }
 274
 275   //--------------------------------------------------------------------
 276   //--------------------------------------------------------------------
 277   // Final Cast
 278   typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 279   typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 280   caster->SetInput(working_image);
 281   caster->Update();
 282   output = caster->GetOutput();
 283
 284   //--------------------------------------------------------------------
 285   //--------------------------------------------------------------------
 286   // [Optional]
 287   if (GetAutoCrop()) {
 288   StartNewStep("AutoCrop");
 289     typedef clitk::AutoCropFilter<MaskImageType> CropFilterType;
 290     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 291     cropFilter->SetInput(output);
 292     cropFilter->SetBackgroundValue(GetBackgroundValue());
 293     cropFilter->Update();
 294     output = cropFilter->GetOutput();
 295     StopCurrentStep<MaskImageType>(output);
 296   }
 297   else
 298   {
 299     // Remove Padding region
 300     typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
 301     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 302     cropFilter->SetInput(output);
 303     cropFilter->SetLowerBoundaryCropSize(bounds);
 304     cropFilter->SetUpperBoundaryCropSize(bounds);
 305     cropFilter->Update();
 306     output = cropFilter->GetOutput();
 307   }
 308 }
 309 //--------------------------------------------------------------------
 310
 311
 312 //--------------------------------------------------------------------
 313 template <class TInputImageType>
 314 void
 315 clitk::ExtractPatientFilter<TInputImageType>::
 316 GenerateData() {
 317   // Final Graft
 318   this->GraftOutput(output);
 319   // Store image filename into AFDB
 320   GetAFDB()->SetImageFilename("Patient", this->GetOutputPatientFilename());
 321   WriteAFDB();
 322 }
 323 //--------------------------------------------------------------------
 324
 325
 326 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX