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