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