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   // But the lungs can touch the outside air. In such case, use primaryOpeningRadius with 1 kernel
 121   typedef itk::ConstantPadImageFilter<InputImageType, InputImageType> PadFilterType;
 122   typename PadFilterType::Pointer padFilter = PadFilterType::New();
 123   padFilter->SetInput(input);
 124   padFilter->SetConstant(GetUpperThreshold() - 1);
 125   typename InputImageType::SizeType bounds;
 126   for (unsigned i = 0; i < Dim - 1; ++i)
 127     bounds[i] = 1;
 128   bounds[Dim - 1] = 0;
 129   padFilter->SetPadLowerBound(bounds);
 130   padFilter->SetPadUpperBound(bounds);
 131   padFilter->Update();
 132
 133   typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> BinarizeFilterType;
 134   typename BinarizeFilterType::Pointer binarizeFilter=BinarizeFilterType::New();
 135   binarizeFilter->SetInput(padFilter->GetOutput());
 136   if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(GetLowerThreshold());
 137   binarizeFilter->SetUpperThreshold(GetUpperThreshold());
 138   binarizeFilter ->SetInsideValue(this->GetForegroundValue());
 139   binarizeFilter ->SetOutsideValue(this->GetBackgroundValue());
 140   padFilter->GetOutput()->ReleaseData();
 141   working_image = binarizeFilter->GetOutput();
 142
 143   typedef itk::BinaryBallStructuringElement<InternalPixelType,Dim> KernelType;
 144   unsigned int radius = this->GetPrimaryOpeningRadius();
 145   if (radius > 0)
 146   {
 147     if (this->GetVerboseOptionFlag()) std::cout << ("Opening after threshold; R = ") << radius << std::endl;
 148     KernelType kernel;
 149     kernel.SetRadius(radius);
 150
 151     typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType2;
 152     typename OpenFilterType2::Pointer openFilter2 = OpenFilterType2::New();
 153     openFilter2->SetInput(working_image);
 154     openFilter2->SetBackgroundValue(0);
 155     openFilter2->SetForegroundValue(1);
 156     openFilter2->SetKernel(kernel);
 157     openFilter2->Update();
 158     working_image->ReleaseData();
 159     working_image = openFilter2->GetOutput();
 160   }
 161
 162   if (this->GetVerboseOptionFlag()) std::cout << ("Labelling") << std::endl;
 163   // Connected component labeling
 164   typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
 165   typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
 166   connectFilter->SetInput(working_image);
 167   connectFilter->SetBackgroundValue(this->GetBackgroundValue());
 168   connectFilter->SetFullyConnected(false);
 169   connectFilter->Update();
 170   working_image->ReleaseData();
 171   working_image = connectFilter->GetOutput();
 172
 173   if (this->GetVerboseOptionFlag()) std::cout << ("RelabelComponentImageFilter") << std::endl;
 174   // Sort labels according to size
 175   typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
 176   typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
 177   relabelFilter->InPlaceOn();
 178   relabelFilter->SetInput(connectFilter->GetOutput());
 179   relabelFilter->Update();
 180   working_image->ReleaseData();
 181   working_image = relabelFilter->GetOutput();
 182
 183   // End
 184   StopCurrentStep<InternalImageType>(working_image);
 185
 186   //--------------------------------------------------------------------
 187   //--------------------------------------------------------------------
 188   // [Optional]
 189   if (GetDecomposeAndReconstructDuringFirstStep()) {
 190     StartNewStep("First Decompose & Reconstruct step");
 191     typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> FilterType;
 192     typename FilterType::Pointer f = FilterType::New();
 193     f->SetInput(working_image);
 194     // f->SetVerbose(m_Verbose);
 195     f->SetRadius(GetRadius1());
 196     f->SetMaximumNumberOfLabels(GetMaximumNumberOfLabels1());
 197     f->SetBackgroundValue(this->GetBackgroundValue());
 198     f->SetForegroundValue(this->GetForegroundValue());
 199     f->SetFullyConnected(true);
 200     f->SetNumberOfNewLabels(GetNumberOfNewLabels1());
 201     f->Update();
 202     working_image->ReleaseData();
 203     working_image = f->GetOutput();
 204     StopCurrentStep<InternalImageType>(working_image);
 205   }
 206
 207   //--------------------------------------------------------------------
 208   //--------------------------------------------------------------------
 209   if (this->GetVerboseOptionFlag()) std::cout << ("Remove the air (largest area)") << std::endl;
 210   StartNewStep("Remove the air (largest area)");
 211   typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> iBinarizeFilterType;
 212   typename iBinarizeFilterType::Pointer binarizeFilter2 = iBinarizeFilterType::New();
 213   binarizeFilter2->SetInput(working_image);
 214   binarizeFilter2->SetLowerThreshold(GetFirstKeep());
 215   binarizeFilter2->SetUpperThreshold(GetLastKeep());
 216   binarizeFilter2 ->SetInsideValue(0);
 217   binarizeFilter2 ->SetOutsideValue(1);
 218   binarizeFilter2 ->Update();
 219   working_image->ReleaseData();
 220   working_image = binarizeFilter2->GetOutput();
 221
 222   typename ConnectFilterType::Pointer connectFilter2 = ConnectFilterType::New();
 223   connectFilter2->SetInput(working_image);
 224   connectFilter2->SetBackgroundValue(this->GetBackgroundValue());
 225   connectFilter2->SetFullyConnected(false);
 226   connectFilter2->Update();
 227   working_image->ReleaseData();
 228   working_image = connectFilter2->GetOutput();
 229
 230   typename RelabelFilterType::Pointer relabelFilter2 = RelabelFilterType::New();
 231   relabelFilter2->SetInput(working_image);
 232   relabelFilter2->Update();
 233   working_image->ReleaseData();
 234   working_image = relabelFilter2->GetOutput();
 235
 236   // Keep main label
 237   working_image = KeepLabels<InternalImageType>
 238     (working_image, GetBackgroundValue(), GetForegroundValue(), 1, 1, true);
 239   StopCurrentStep<InternalImageType>(working_image);
 240
 241   //--------------------------------------------------------------------
 242   //--------------------------------------------------------------------
 243   // [Optional]
 244   if (GetDecomposeAndReconstructDuringSecondStep()) {
 245     StartNewStep("Second Decompose & Reconstruct step");
 246     typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> FilterType;
 247     typename FilterType::Pointer f = FilterType::New();
 248     f->SetInput(working_image);
 249     // f->SetVerbose(m_Verbose);
 250     f->SetRadius(GetRadius2());
 251     f->SetMaximumNumberOfLabels(GetMaximumNumberOfLabels2());
 252     f->SetBackgroundValue(this->GetBackgroundValue());
 253     f->SetForegroundValue(this->GetForegroundValue());
 254     f->SetFullyConnected(true);
 255     f->SetNumberOfNewLabels(GetNumberOfNewLabels2());
 256     f->Update();
 257     working_image->ReleaseData();
 258     working_image = f->GetOutput();
 259     StopCurrentStep<InternalImageType>(working_image);
 260   }
 261
 262   //--------------------------------------------------------------------
 263   //--------------------------------------------------------------------
 264   // [Optional]
 265   if (GetFinalOpenClose()) {
 266     StartNewStep("Final OpenClose");
 267     // Open
 268     KernelType structuringElement;
 269     structuringElement.SetRadius(1);
 270     structuringElement.CreateStructuringElement();
 271     typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType;
 272     typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
 273     openFilter->SetInput(working_image);
 274     openFilter->SetBackgroundValue(this->GetBackgroundValue());
 275     openFilter->SetForegroundValue(this->GetForegroundValue());
 276     openFilter->SetKernel(structuringElement);
 277     // Close
 278     typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType , KernelType> CloseFilterType;
 279     typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
 280     closeFilter->SetInput(openFilter->GetOutput());
 281     closeFilter->SetSafeBorder(true);
 282     closeFilter->SetForegroundValue(this->GetForegroundValue());
 283     //  closeFilter->SetBackgroundValue(SetBackgroundValue());
 284     closeFilter->SetKernel(structuringElement);
 285     closeFilter->Update();
 286     working_image->ReleaseData();
 287     working_image = closeFilter->GetOutput();
 288     StopCurrentStep<InternalImageType>(working_image);
 289   }
 290
 291   //--------------------------------------------------------------------
 292   //--------------------------------------------------------------------
 293   // Final Cast
 294   typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 295   typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 296   caster->SetInput(working_image);
 297   caster->Update();
 298   working_image->ReleaseData();
 299   output = caster->GetOutput();
 300
 301   //--------------------------------------------------------------------
 302   //--------------------------------------------------------------------
 303   // [Optional]
 304   if (GetAutoCrop()) {
 305   StartNewStep("AutoCrop");
 306     typedef clitk::AutoCropFilter<MaskImageType> CropFilterType;
 307     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 308     cropFilter->SetInput(output);
 309     cropFilter->SetBackgroundValue(GetBackgroundValue());
 310     cropFilter->Update();
 311     output->ReleaseData();
 312     output = cropFilter->GetOutput();
 313     StopCurrentStep<MaskImageType>(output);
 314   }
 315   else
 316   {
 317     // Remove Padding region
 318     typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
 319     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 320     cropFilter->SetInput(output);
 321     cropFilter->SetLowerBoundaryCropSize(bounds);
 322     cropFilter->SetUpperBoundaryCropSize(bounds);
 323     cropFilter->Update();
 324     output->ReleaseData();
 325     output = cropFilter->GetOutput();
 326   }
 327 }
 328 //--------------------------------------------------------------------
 329
 330
 331 //--------------------------------------------------------------------
 332 template <class TInputImageType>
 333 void
 334 clitk::ExtractPatientFilter<TInputImageType>::
 335 GenerateData() {
 336   // Final Graft
 337   this->GraftOutput(output);
 338   // Store image filename into AFDB
 339   GetAFDB()->SetImageFilename("Patient", this->GetOutputPatientFilename());
 340   WriteAFDB();
 341 }
 342 //--------------------------------------------------------------------
 343
 344
 345 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX