segmentation/clitkExtractLungFilter.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 CLITKEXTRACTLUNGSFILTER_TXX
  20 #define CLITKEXTRACTLUNGSFILTER_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 "itkOtsuThresholdImageFilter.h"
  33
  34 //--------------------------------------------------------------------
  35 template <class TInputImageType, class TMaskImageType>
  36 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
  37 ExtractLungFilter():
  38   clitk::FilterBase(),
  39   itk::ImageToImageFilter<TInputImageType, TMaskImageType>()
  40 {
  41   // Default global options
  42   this->SetNumberOfRequiredInputs(2);
  43   SetPatientMaskBackgroundValue(0);
  44   SetBackgroundValue(0); // Must be zero
  45   SetForegroundValue(1);
  46
  47   // Step 1 default values
  48   SetUpperThreshold(-300);
  49   SetLowerThreshold(-1000);
  50   UseLowerThresholdOff();
  51   LabelParamType * p1 = new LabelParamType;
  52   p1->SetFirstKeep(1);
  53   p1->UseLastKeepOff();
  54   p1->AddLabelToRemove(2);
  55   SetLabelizeParameters1(p1);
  56
  57   // Step 2 default values
  58   SetUpperThresholdForTrachea(-900);
  59   SetMultiplierForTrachea(5);
  60   SetThresholdStepSizeForTrachea(64);
  61
  62   // Step 3 default values
  63   SetNumberOfHistogramBins(500);
  64   LabelParamType * p2 = new LabelParamType;
  65   p2->SetFirstKeep(1);
  66   p2->UseLastKeepOff();
  67   // p->AddLabelToRemove(2); // No label to remove by default
  68   SetLabelizeParameters2(p2);
  69
  70   // Step 4 default values
  71   SetRadiusForTrachea(1);
  72   LabelParamType * p3 = new LabelParamType;
  73   p3->SetFirstKeep(1);
  74   p3->SetLastKeep(2);
  75   p3->UseLastKeepOff();
  76   SetLabelizeParameters3(p3);
  77 }
  78 //--------------------------------------------------------------------
  79
  80
  81 //--------------------------------------------------------------------
  82 template <class TInputImageType, class TMaskImageType>
  83 void
  84 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
  85 SetInput(const TInputImageType * image)
  86 {
  87   this->SetNthInput(0, const_cast<TInputImageType *>(image));
  88 }
  89 //--------------------------------------------------------------------
  90
  91
  92 //--------------------------------------------------------------------
  93 template <class TInputImageType, class TMaskImageType>
  94 void
  95 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
  96 SetInputPatientMask(TMaskImageType * image, MaskImagePixelType bg )
  97 {
  98   this->SetNthInput(1, const_cast<TMaskImageType *>(image));
  99   SetPatientMaskBackgroundValue(bg);
 100 }
 101 //--------------------------------------------------------------------
 102
 103
 104 //--------------------------------------------------------------------
 105 template <class TInputImageType, class TMaskImageType>
 106 void
 107 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
 108 AddSeed(InternalIndexType s)
 109 {
 110   m_Seeds.push_back(s);
 111 }
 112 //--------------------------------------------------------------------
 113
 114
 115 //--------------------------------------------------------------------
 116 template <class TInputImageType, class TMaskImageType>
 117 template<class ArgsInfoType>
 118 void
 119 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
 120 SetArgsInfo(ArgsInfoType mArgsInfo)
 121 {
 122   SetVerboseOption_GGO(mArgsInfo);
 123   SetVerboseStep_GGO(mArgsInfo);
 124   SetWriteStep_GGO(mArgsInfo);
 125   SetVerboseWarningOff_GGO(mArgsInfo);
 126
 127   SetUpperThreshold_GGO(mArgsInfo);
 128   SetLowerThreshold_GGO(mArgsInfo);
 129   SetLabelizeParameters1_GGO(mArgsInfo);
 130   if (!mArgsInfo.remove1_given) {
 131     GetLabelizeParameters1()->AddLabelToRemove(2);
 132     if (GetVerboseOption()) VerboseOption("remove1", 2);
 133   }
 134
 135   SetUpperThresholdForTrachea_GGO(mArgsInfo);
 136   SetMultiplierForTrachea_GGO(mArgsInfo);
 137   SetThresholdStepSizeForTrachea_GGO(mArgsInfo);
 138   AddSeed_GGO(mArgsInfo);
 139
 140   SetNumberOfHistogramBins_GGO(mArgsInfo);
 141   SetLabelizeParameters2_GGO(mArgsInfo);
 142
 143   SetRadiusForTrachea_GGO(mArgsInfo);
 144   SetLabelizeParameters3_GGO(mArgsInfo);
 145 }
 146 //--------------------------------------------------------------------
 147
 148
 149 //--------------------------------------------------------------------
 150 template <class TInputImageType, class TMaskImageType>
 151 void
 152 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
 153 GenerateOutputInformation()
 154 {
 155   input   = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
 156   Superclass::GenerateOutputInformation();
 157 //   MaskImagePointer output = this->GetOutput(0);
 158
 159   // Get input pointers
 160   input   = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
 161   patient = dynamic_cast<const TMaskImageType*>(itk::ProcessObject::GetInput(1));
 162
 163   // Check image
 164   if (!HaveSameSizeAndSpacing<TInputImageType, TMaskImageType>(input, patient)) {
 165     this->SetLastError("* ERROR * the images (input and patient mask) must have the same size & spacing");
 166     return;
 167   }
 168
 169   //--------------------------------------------------------------------
 170   //--------------------------------------------------------------------
 171   StartNewStep("Set background to initial image");
 172   working_input = SetBackground<TInputImageType, TMaskImageType>
 173     (input, patient, GetPatientMaskBackgroundValue(), -1000);
 174   StopCurrentStep<InputImageType>(working_input);
 175
 176   //--------------------------------------------------------------------
 177   //--------------------------------------------------------------------
 178   StartNewStep("Remove Air");
 179   // Threshold to get air
 180   typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
 181   typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 182   binarizeFilter->SetInput(working_input);
 183   if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
 184   binarizeFilter->SetUpperThreshold(m_UpperThreshold);
 185   binarizeFilter ->SetInsideValue(this->GetForegroundValue());
 186   binarizeFilter ->SetOutsideValue(this->GetBackgroundValue());
 187   binarizeFilter->Update();
 188   working_image = binarizeFilter->GetOutput();
 189
 190   // Labelize and keep right labels
 191   working_image = Labelize<InternalImageType>(working_image, GetBackgroundValue(), true, GetMinimalComponentSize());
 192   working_image = RemoveLabels<InternalImageType>
 193     (working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
 194   typename InternalImageType::Pointer air = KeepLabels<InternalImageType>
 195     (working_image,
 196      GetBackgroundValue(),
 197      GetForegroundValue(),
 198      GetLabelizeParameters1()->GetFirstKeep(),
 199      GetLabelizeParameters1()->GetLastKeep(),
 200      GetLabelizeParameters1()->GetUseLastKeep());
 201
 202   // Set Air to BG
 203   working_input = SetBackground<TInputImageType, InternalImageType>
 204     (working_input, air, this->GetForegroundValue(), this->GetBackgroundValue());
 205
 206   // End
 207   StopCurrentStep<InputImageType>(working_input);
 208
 209   //--------------------------------------------------------------------
 210   //--------------------------------------------------------------------
 211   StartNewStep("Find the trachea");
 212   //DD(m_Seeds.size());
 213   if (m_Seeds.size() == 0) { // try to find seed
 214     // Search seed (parameters = UpperThresholdForTrachea)
 215     static const unsigned int Dim = InputImageType::ImageDimension;
 216     typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
 217     typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
 218     typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
 219     sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-5;
 220     sliceRegionSize[Dim-1]=5;
 221     sliceRegion.SetSize(sliceRegionSize);
 222     sliceRegion.SetIndex(sliceRegionIndex);
 223     //DD(GetUpperThresholdForTrachea());
 224     //DD(sliceRegion);
 225     typedef  itk::ImageRegionConstIterator<InputImageType> IteratorType;
 226     IteratorType it(working_input, sliceRegion);
 227     it.GoToBegin();
 228     while (!it.IsAtEnd()) {
 229       if(it.Get() < GetUpperThresholdForTrachea() ) {
 230         AddSeed(it.GetIndex());
 231         //      DD(it.GetIndex());
 232       }
 233       ++it;
 234     }
 235   }
 236
 237   //DD(m_Seeds.size());
 238   if (m_Seeds.size() != 0) {
 239     // Explosion controlled region growing
 240     typedef clitk::ExplosionControlledThresholdConnectedImageFilter<InputImageType, InternalImageType> ImageFilterType;
 241     typename ImageFilterType::Pointer f= ImageFilterType::New();
 242     f->SetInput(working_input);
 243     f->SetVerbose(false);
 244     f->SetLower(-2000);
 245     f->SetUpper(GetUpperThresholdForTrachea());
 246     f->SetMinimumLowerThreshold(-2000);
 247     f->SetMaximumUpperThreshold(0);
 248     f->SetAdaptLowerBorder(false);
 249     f->SetAdaptUpperBorder(true);
 250     f->SetMinimumSize(5000);
 251     f->SetReplaceValue(1);
 252     f->SetMultiplier(GetMultiplierForTrachea());
 253     f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
 254     f->SetMinimumThresholdStepSize(1);
 255     for(unsigned int i=0; i<m_Seeds.size();i++) {
 256       //    std::cout<<"Adding seed " <<m_Seeds[i]<<"..."<<std::endl;
 257       f->AddSeed(m_Seeds[i]);
 258     }
 259     f->Update();
 260     trachea_tmp = f->GetOutput();
 261     // Set output
 262     StopCurrentStep<InternalImageType>(trachea_tmp);
 263
 264   }
 265   else { // Trachea not found
 266     this->SetWarning("* WARNING * No seed found for trachea.");
 267     // Set output
 268     StopCurrentStep();
 269   }
 270
 271   //--------------------------------------------------------------------
 272   //--------------------------------------------------------------------
 273   StartNewStep("Extract the lung with Otsu filter");
 274   // Automated Otsu thresholding and relabeling
 275   typedef itk::OtsuThresholdImageFilter<InputImageType,InternalImageType> OtsuThresholdImageFilterType;
 276   typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
 277   otsuFilter->SetInput(working_input);
 278   otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
 279   otsuFilter->SetInsideValue(this->GetForegroundValue());
 280   otsuFilter->SetOutsideValue(this->GetBackgroundValue());
 281   otsuFilter->Update();
 282   working_image =  otsuFilter->GetOutput();
 283
 284   // Set output
 285   StopCurrentStep<InternalImageType>(working_image);
 286
 287   //--------------------------------------------------------------------
 288   //--------------------------------------------------------------------
 289   StartNewStep("Select labels");
 290   // Keep right labels
 291   working_image = LabelizeAndSelectLabels<InternalImageType>
 292     (working_image,
 293      GetBackgroundValue(),
 294      GetForegroundValue(),
 295      false,
 296      GetMinimalComponentSize(),
 297      GetLabelizeParameters2());
 298
 299   // Set output
 300   StopCurrentStep<InternalImageType>(working_image);
 301
 302   //--------------------------------------------------------------------
 303   //--------------------------------------------------------------------
 304   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 305     StartNewStep("Remove the trachea");
 306     // Set the trachea
 307     working_image = SetBackground<InternalImageType, InternalImageType>
 308       (working_image, trachea_tmp, 1, -1);
 309
 310    // Dilate the trachea
 311     static const unsigned int Dim = InputImageType::ImageDimension;
 312     typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
 313     KernelType structuringElement;
 314     structuringElement.SetRadius(GetRadiusForTrachea());
 315     structuringElement.CreateStructuringElement();
 316     typedef clitk::ConditionalBinaryDilateImageFilter<InternalImageType, InternalImageType, KernelType> ConditionalBinaryDilateImageFilterType;
 317     typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
 318     dilateFilter->SetBoundaryToForeground(false);
 319     dilateFilter->SetKernel(structuringElement);
 320     dilateFilter->SetBackgroundValue (1);
 321     dilateFilter->SetForegroundValue (-1);
 322     dilateFilter->SetInput (working_image);
 323     dilateFilter->Update();
 324     working_image = dilateFilter->GetOutput();
 325
 326     // Set trachea with dilatation
 327     trachea_tmp = SetBackground<InternalImageType, InternalImageType>
 328       (trachea_tmp, working_image, -1, this->GetForegroundValue());
 329
 330     // Remove the trachea
 331     working_image = SetBackground<InternalImageType, InternalImageType>
 332       (working_image, working_image, -1, this->GetBackgroundValue());
 333
 334     // Label
 335     working_image = LabelizeAndSelectLabels<InternalImageType>
 336       (working_image,
 337        GetBackgroundValue(),
 338        GetForegroundValue(),
 339        false,
 340        GetMinimalComponentSize(),
 341        GetLabelizeParameters3());
 342
 343     // Set output
 344     StopCurrentStep<InternalImageType>(working_image);
 345   }
 346
 347
 348   //--------------------------------------------------------------------
 349   //--------------------------------------------------------------------
 350   typedef clitk::AutoCropFilter<InternalImageType> CropFilterType;
 351   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 352   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 353     StartNewStep("Croping trachea");
 354     cropFilter->SetInput(trachea_tmp);
 355     cropFilter->Update(); // Needed
 356     typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 357     typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 358     caster->SetInput(cropFilter->GetOutput());
 359     caster->Update();
 360     trachea = caster->GetOutput();
 361     StopCurrentStep<MaskImageType>(trachea);
 362   }
 363
 364   //--------------------------------------------------------------------
 365   //--------------------------------------------------------------------
 366   StartNewStep("Croping lung");
 367   typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline
 368   cropFilter2->SetInput(working_image);
 369   cropFilter2->Update();
 370   working_image = cropFilter2->GetOutput();
 371   StopCurrentStep<InternalImageType>(working_image);
 372
 373   //--------------------------------------------------------------------
 374   //--------------------------------------------------------------------
 375   StartNewStep("Separate Left/Right lungs");
 376   // Initial label
 377   working_image = Labelize<InternalImageType>(working_image,
 378                                               GetBackgroundValue(),
 379                                               false,
 380                                               GetMinimalComponentSize());
 381
 382   // Count the labels
 383   typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
 384   typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
 385   statisticsImageFilter->SetInput(working_image);
 386   statisticsImageFilter->Update();
 387   unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
 388   working_image = statisticsImageFilter->GetOutput();
 389
 390   // Decompose the first label
 391   static const unsigned int Dim = InputImageType::ImageDimension;
 392   if (initialNumberOfLabels<2) {
 393     // Structuring element radius
 394     typename InputImageType::SizeType radius;
 395     for (unsigned int i=0;i<Dim;i++) radius[i]=1;
 396     typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> DecomposeAndReconstructFilterType;
 397     typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
 398     decomposeAndReconstructFilter->SetInput(working_image);
 399     decomposeAndReconstructFilter->SetVerbose(false);
 400     decomposeAndReconstructFilter->SetRadius(radius);
 401     decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
 402     decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
 403     decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
 404     decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
 405     decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
 406     decomposeAndReconstructFilter->SetFullyConnected(true);
 407     decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
 408     decomposeAndReconstructFilter->Update();
 409     working_image = decomposeAndReconstructFilter->GetOutput();
 410   }
 411
 412   // Retain labels (lungs)
 413   typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
 414   typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
 415   thresholdFilter->SetInput(working_image);
 416   thresholdFilter->ThresholdAbove(2);
 417   thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
 418   thresholdFilter->Update();
 419   working_image = thresholdFilter->GetOutput();
 420   StopCurrentStep<InternalImageType> (working_image);
 421 }
 422 //--------------------------------------------------------------------
 423
 424
 425 //--------------------------------------------------------------------
 426 template <class TInputImageType, class TMaskImageType>
 427 void
 428 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
 429 GenerateData() {
 430   // Final Cast
 431   typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 432   typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 433   caster->SetInput(working_image);
 434   caster->Update();
 435   // Set output
 436   //this->SetNthOutput(0, caster->GetOutput()); // -> no because redo filter otherwise
 437   this->GraftOutput(caster->GetOutput());
 438   return;
 439 }
 440 //--------------------------------------------------------------------
 441
 442
 443 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX