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 #include "clitkExtractSliceFilter.h"
  28
  29 // itk
  30 #include "itkBinaryThresholdImageFilter.h"
  31 #include "itkConnectedComponentImageFilter.h"
  32 #include "itkRelabelComponentImageFilter.h"
  33 #include "itkOtsuThresholdImageFilter.h"
  34 #include "itkBinaryThinningImageFilter3D.h"
  35 #include "itkImageIteratorWithIndex.h"
  36
  37 //--------------------------------------------------------------------
  38 template <class ImageType, class MaskImageType>
  39 clitk::ExtractLungFilter<ImageType, MaskImageType>::
  40 ExtractLungFilter():
  41   clitk::FilterBase(),
  42   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  43   itk::ImageToImageFilter<ImageType, MaskImageType>()
  44 {
  45   SetNumberOfSteps(10);
  46   m_MaxSeedNumber = 500;
  47
  48   // Default global options
  49   this->SetNumberOfRequiredInputs(2);
  50   SetPatientMaskBackgroundValue(0);
  51   SetBackgroundValue(0); // Must be zero
  52   SetForegroundValue(1);
  53   SetMinimalComponentSize(100);
  54
  55   // Step 1 default values
  56   SetUpperThreshold(-300);
  57   SetLowerThreshold(-1000);
  58   UseLowerThresholdOff();
  59   LabelParamType * p1 = new LabelParamType;
  60   p1->SetFirstKeep(1);
  61   p1->UseLastKeepOff();
  62   p1->AddLabelToRemove(2);
  63   SetLabelizeParameters1(p1);
  64
  65   // Step 2 default values
  66   SetUpperThresholdForTrachea(-900);
  67   SetMultiplierForTrachea(5);
  68   SetThresholdStepSizeForTrachea(64);
  69   SetNumberOfSlicesToSkipBeforeSearchingSeed(0);
  70
  71   // Step 3 default values
  72   SetNumberOfHistogramBins(500);
  73   LabelParamType * p2 = new LabelParamType;
  74   p2->SetFirstKeep(1);
  75   p2->UseLastKeepOff();
  76   // p->AddLabelToRemove(2); // No label to remove by default
  77   SetLabelizeParameters2(p2);
  78
  79   // Step 4 default values
  80   SetRadiusForTrachea(1);
  81   LabelParamType * p3 = new LabelParamType;
  82   p3->SetFirstKeep(1);
  83   p3->SetLastKeep(2);
  84   p3->UseLastKeepOff();
  85   SetLabelizeParameters3(p3);
  86
  87   // Step 5 : find bronchial bifurcations
  88   FindBronchialBifurcationsOn();
  89 }
  90 //--------------------------------------------------------------------
  91
  92
  93 //--------------------------------------------------------------------
  94 template <class ImageType, class MaskImageType>
  95 void
  96 clitk::ExtractLungFilter<ImageType, MaskImageType>::
  97 SetInput(const ImageType * image)
  98 {
  99   this->SetNthInput(0, const_cast<ImageType *>(image));
 100 }
 101 //--------------------------------------------------------------------
 102
 103
 104 //--------------------------------------------------------------------
 105 template <class ImageType, class MaskImageType>
 106 void
 107 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 108 SetInputPatientMask(MaskImageType * image, MaskImagePixelType bg )
 109 {
 110   this->SetNthInput(1, const_cast<MaskImageType *>(image));
 111   SetPatientMaskBackgroundValue(bg);
 112 }
 113 //--------------------------------------------------------------------
 114
 115
 116 //--------------------------------------------------------------------
 117 template <class ImageType, class MaskImageType>
 118 void
 119 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 120 AddSeed(InternalIndexType s)
 121 {
 122   m_Seeds.push_back(s);
 123 }
 124 //--------------------------------------------------------------------
 125
 126
 127 //--------------------------------------------------------------------
 128 template <class ImageType, class MaskImageType>
 129 template<class ArgsInfoType>
 130 void
 131 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 132 SetArgsInfo(ArgsInfoType mArgsInfo)
 133 {
 134   SetVerboseOption_GGO(mArgsInfo);
 135   SetVerboseStep_GGO(mArgsInfo);
 136   SetWriteStep_GGO(mArgsInfo);
 137   SetVerboseWarningOff_GGO(mArgsInfo);
 138
 139   SetUpperThreshold_GGO(mArgsInfo);
 140   SetLowerThreshold_GGO(mArgsInfo);
 141   SetNumberOfSlicesToSkipBeforeSearchingSeed_GGO(mArgsInfo);
 142   SetLabelizeParameters1_GGO(mArgsInfo);
 143   if (!mArgsInfo.remove1_given) {
 144     GetLabelizeParameters1()->AddLabelToRemove(2);
 145     if (GetVerboseOption()) VerboseOption("remove1", 2);
 146   }
 147
 148   SetUpperThresholdForTrachea_GGO(mArgsInfo);
 149   SetMultiplierForTrachea_GGO(mArgsInfo);
 150   SetThresholdStepSizeForTrachea_GGO(mArgsInfo);
 151   AddSeed_GGO(mArgsInfo);
 152
 153   SetMinimalComponentSize_GGO(mArgsInfo);
 154
 155   SetNumberOfHistogramBins_GGO(mArgsInfo);
 156   SetLabelizeParameters2_GGO(mArgsInfo);
 157
 158   SetRadiusForTrachea_GGO(mArgsInfo);
 159   SetLabelizeParameters3_GGO(mArgsInfo);
 160
 161   SetAFDBFilename_GGO(mArgsInfo);
 162 }
 163 //--------------------------------------------------------------------
 164
 165
 166 //--------------------------------------------------------------------
 167 template <class ImageType, class MaskImageType>
 168 void
 169 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 170 GenerateOutputInformation()
 171 {
 172   Superclass::GenerateOutputInformation();
 173   //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion());
 174
 175   // Get input pointers
 176   patient = dynamic_cast<const MaskImageType*>(itk::ProcessObject::GetInput(1));
 177   input   = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
 178
 179   // Check image
 180   if (!HaveSameSizeAndSpacing<ImageType, MaskImageType>(input, patient)) {
 181     clitkExceptionMacro("the 'input' and 'patient' masks must have the same size & spacing.");
 182   }
 183
 184   //--------------------------------------------------------------------
 185   //--------------------------------------------------------------------
 186   StartNewStep("Set background to initial image");
 187   working_input = SetBackground<ImageType, MaskImageType>
 188     (input, patient, GetPatientMaskBackgroundValue(), -1000);
 189   StopCurrentStep<ImageType>(working_input);
 190
 191   //--------------------------------------------------------------------
 192   //--------------------------------------------------------------------
 193   StartNewStep("Remove Air");
 194   // Check threshold
 195   if (m_UseLowerThreshold) {
 196     if (m_LowerThreshold > m_UpperThreshold) {
 197       clitkExceptionMacro("lower threshold cannot be greater than upper threshold.");
 198     }
 199   }
 200   // Threshold to get air
 201   typedef itk::BinaryThresholdImageFilter<ImageType, InternalImageType> InputBinarizeFilterType;
 202   typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 203   binarizeFilter->SetInput(working_input);
 204   if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
 205   binarizeFilter->SetUpperThreshold(m_UpperThreshold);
 206   binarizeFilter ->SetInsideValue(this->GetForegroundValue());
 207   binarizeFilter ->SetOutsideValue(this->GetBackgroundValue());
 208   binarizeFilter->Update();
 209   working_image = binarizeFilter->GetOutput();
 210
 211   // Labelize and keep right labels
 212   working_image = Labelize<InternalImageType>(working_image, GetBackgroundValue(), true, GetMinimalComponentSize());
 213
 214   working_image = RemoveLabels<InternalImageType>
 215     (working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
 216
 217   typename InternalImageType::Pointer air = KeepLabels<InternalImageType>
 218     (working_image,
 219      GetBackgroundValue(),
 220      GetForegroundValue(),
 221      GetLabelizeParameters1()->GetFirstKeep(),
 222      GetLabelizeParameters1()->GetLastKeep(),
 223      GetLabelizeParameters1()->GetUseLastKeep());
 224
 225   // Set Air to BG
 226   working_input = SetBackground<ImageType, InternalImageType>
 227     (working_input, air, this->GetForegroundValue(), this->GetBackgroundValue());
 228
 229   // End
 230   StopCurrentStep<ImageType>(working_input);
 231
 232   //--------------------------------------------------------------------
 233   //--------------------------------------------------------------------
 234   StartNewStep("Search for the trachea");
 235   SearchForTrachea();
 236
 237   //--------------------------------------------------------------------
 238   //--------------------------------------------------------------------
 239   StartNewStep("Extract the lung with Otsu filter");
 240   // Automated Otsu thresholding and relabeling
 241   typedef itk::OtsuThresholdImageFilter<ImageType,InternalImageType> OtsuThresholdImageFilterType;
 242   typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
 243   otsuFilter->SetInput(working_input);
 244   otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
 245   otsuFilter->SetInsideValue(this->GetForegroundValue());
 246   otsuFilter->SetOutsideValue(this->GetBackgroundValue());
 247   otsuFilter->Update();
 248   working_image =  otsuFilter->GetOutput();
 249
 250   // Set output
 251   StopCurrentStep<InternalImageType>(working_image);
 252
 253   //--------------------------------------------------------------------
 254   //--------------------------------------------------------------------
 255   StartNewStep("Select labels");
 256   // Keep right labels
 257   working_image = LabelizeAndSelectLabels<InternalImageType>
 258     (working_image,
 259      GetBackgroundValue(),
 260      GetForegroundValue(),
 261      false,
 262      GetMinimalComponentSize(),
 263      GetLabelizeParameters2());
 264
 265   // Set output
 266   StopCurrentStep<InternalImageType>(working_image);
 267
 268   //--------------------------------------------------------------------
 269   //--------------------------------------------------------------------
 270   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 271     StartNewStep("Remove the trachea");
 272     // Set the trachea
 273     working_image = SetBackground<InternalImageType, InternalImageType>
 274       (working_image, trachea_tmp, 1, -1);
 275
 276     // Dilate the trachea
 277     static const unsigned int Dim = ImageType::ImageDimension;
 278     typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
 279     KernelType structuringElement;
 280     structuringElement.SetRadius(GetRadiusForTrachea());
 281     structuringElement.CreateStructuringElement();
 282     typedef clitk::ConditionalBinaryDilateImageFilter<InternalImageType, InternalImageType, KernelType> ConditionalBinaryDilateImageFilterType;
 283     typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
 284     dilateFilter->SetBoundaryToForeground(false);
 285     dilateFilter->SetKernel(structuringElement);
 286     dilateFilter->SetBackgroundValue (1);
 287     dilateFilter->SetForegroundValue (-1);
 288     dilateFilter->SetInput (working_image);
 289     dilateFilter->Update();
 290     working_image = dilateFilter->GetOutput();
 291
 292     // Set trachea with dilatation
 293     trachea_tmp = SetBackground<InternalImageType, InternalImageType>
 294       (trachea_tmp, working_image, -1, this->GetForegroundValue());
 295
 296     // Remove the trachea
 297     working_image = SetBackground<InternalImageType, InternalImageType>
 298       (working_image, working_image, -1, this->GetBackgroundValue());
 299
 300     // Label
 301     working_image = LabelizeAndSelectLabels<InternalImageType>
 302       (working_image,
 303        GetBackgroundValue(),
 304        GetForegroundValue(),
 305        false,
 306        GetMinimalComponentSize(),
 307        GetLabelizeParameters3());
 308
 309     // Set output
 310     StopCurrentStep<InternalImageType>(working_image);
 311   }
 312
 313   //--------------------------------------------------------------------
 314   //--------------------------------------------------------------------
 315   typedef clitk::AutoCropFilter<InternalImageType> CropFilterType;
 316   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 317   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 318     StartNewStep("Croping trachea");
 319     cropFilter->SetInput(trachea_tmp);
 320     cropFilter->Update(); // Needed
 321     typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 322     typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 323     caster->SetInput(cropFilter->GetOutput());
 324     caster->Update();
 325     trachea = caster->GetOutput();
 326     StopCurrentStep<MaskImageType>(trachea);
 327   }
 328
 329   //--------------------------------------------------------------------
 330   //--------------------------------------------------------------------
 331   StartNewStep("Croping lung");
 332   typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline
 333   cropFilter2->SetInput(working_image);
 334   cropFilter2->Update();
 335   working_image = cropFilter2->GetOutput();
 336   StopCurrentStep<InternalImageType>(working_image);
 337
 338   //--------------------------------------------------------------------
 339   //--------------------------------------------------------------------
 340   StartNewStep("Separate Left/Right lungs");
 341   // Initial label
 342   working_image = Labelize<InternalImageType>(working_image,
 343                                               GetBackgroundValue(),
 344                                               false,
 345                                               GetMinimalComponentSize());
 346
 347   // Count the labels
 348   typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
 349   typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
 350   statisticsImageFilter->SetInput(working_image);
 351   statisticsImageFilter->Update();
 352   unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
 353   working_image = statisticsImageFilter->GetOutput();
 354
 355   // Decompose the first label
 356   static const unsigned int Dim = ImageType::ImageDimension;
 357   if (initialNumberOfLabels<2) {
 358     // Structuring element radius
 359     typename ImageType::SizeType radius;
 360     for (unsigned int i=0;i<Dim;i++) radius[i]=1;
 361     typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> DecomposeAndReconstructFilterType;
 362     typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
 363     decomposeAndReconstructFilter->SetInput(working_image);
 364     decomposeAndReconstructFilter->SetVerbose(false);
 365     decomposeAndReconstructFilter->SetRadius(radius);
 366     decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
 367     decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
 368     decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
 369     decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
 370     decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
 371     decomposeAndReconstructFilter->SetFullyConnected(true);
 372     decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
 373     decomposeAndReconstructFilter->Update();
 374     working_image = decomposeAndReconstructFilter->GetOutput();
 375   }
 376
 377   // Retain labels ('1' is largset lung, so right. '2' is left)
 378   typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
 379   typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
 380   thresholdFilter->SetInput(working_image);
 381   thresholdFilter->ThresholdAbove(2);
 382   thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
 383   thresholdFilter->Update();
 384   working_image = thresholdFilter->GetOutput();
 385   StopCurrentStep<InternalImageType> (working_image);
 386
 387   // Final Cast
 388   StartNewStep("Cast the lung mask");
 389   typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 390   typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 391   caster->SetInput(working_image);
 392   caster->Update();
 393   output = caster->GetOutput();
 394
 395   // Update output info
 396   this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
 397
 398   // Try to extract bifurcation in the trachea (bronchi)
 399   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 400
 401     if (GetFindBronchialBifurcations()) {
 402       StartNewStep("Find bronchial bifurcations");
 403       // Step 1 : extract skeleton
 404       typedef itk::BinaryThinningImageFilter3D<MaskImageType, MaskImageType> ThinningFilterType;
 405       typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
 406       thinningFilter->SetInput(trachea);
 407       thinningFilter->Update();
 408       typename MaskImageType::Pointer skeleton = thinningFilter->GetOutput();
 409
 410       // Step 2.1 : find first point for tracking
 411       typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> IteratorType;
 412       IteratorType it(skeleton, skeleton->GetLargestPossibleRegion());
 413       it.GoToReverseBegin();
 414       while ((!it.IsAtEnd()) && (it.Get() == GetBackgroundValue())) {
 415         --it;
 416       }
 417       if (it.IsAtEnd()) {
 418         clitkExceptionMacro("first point in the trachea skeleton not found.");
 419       }
 420       typename MaskImageType::IndexType index = it.GetIndex();
 421
 422       // Step 2.2 : initialize neighborhooditerator
 423       typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
 424       typename NeighborhoodIteratorType::SizeType radius;
 425       radius.Fill(1);
 426       NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 427
 428       // Find first label number (must be different from BG and FG)
 429       typename MaskImageType::PixelType label = GetForegroundValue()+1;
 430       while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; }
 431
 432       // Track from the first point
 433       std::vector<BifurcationType> listOfBifurcations;
 434       m_SkeletonTree.set_head(index);
 435       TrackFromThisIndex(listOfBifurcations, skeleton, index, label, m_SkeletonTree.begin());
 436       DD("end track");
 437       DD(listOfBifurcations.size());
 438       DD(m_SkeletonTree.size());
 439
 440       for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
 441         skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index,
 442                                                 listOfBifurcations[i].point);
 443       }
 444
 445       // Search for the first slice that separate the bronchus (carena)
 446       typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
 447       typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
 448       extractSliceFilter->SetInput(trachea);
 449       extractSliceFilter->SetDirection(2);
 450       extractSliceFilter->Update();
 451       typedef typename ExtractSliceFilterType::SliceType SliceType;
 452       std::vector<typename SliceType::Pointer> mInputSlices;
 453       extractSliceFilter->GetOutputSlices(mInputSlices);
 454
 455       bool stop = false;
 456       DD(listOfBifurcations[0].index);
 457       DD(listOfBifurcations[1].index);
 458       int slice_index = listOfBifurcations[0].index[2]; // first slice from carena in skeleton
 459       int i=0;
 460       TreeIterator firstIter = m_SkeletonTree.child(listOfBifurcations[1].treeIter, 0);
 461       TreeIterator secondIter = m_SkeletonTree.child(listOfBifurcations[1].treeIter, 1);
 462       typename SliceType::IndexType in1;
 463       typename SliceType::IndexType in2;
 464       while (!stop) {
 465         DD(slice_index);
 466
 467         //  Labelize the current slice
 468         typename SliceType::Pointer temp = Labelize<SliceType>(mInputSlices[slice_index],
 469                                                                GetBackgroundValue(),
 470                                                                true,
 471                                                                GetMinimalComponentSize());
 472         // Check the value of the two skeleton points;
 473         in1[0] = (*firstIter)[0];
 474         in1[1] = (*firstIter)[1];
 475         typename SliceType::PixelType v1 = temp->GetPixel(in1);
 476         DD(in1);
 477         DD((int)v1);
 478         in2[0] = (*secondIter)[0];
 479         in2[1] = (*secondIter)[1];
 480         typename SliceType::PixelType v2 = temp->GetPixel(in2);
 481         DD(in2);
 482         DD((int)v2);
 483
 484         // TODO IF NOT FOUND ????
 485
 486         if (v1 != v2) {
 487           stop = true;
 488         }
 489         else {
 490           i++;
 491           --slice_index;
 492           ++firstIter;
 493           ++secondIter;
 494         }
 495       }
 496       MaskImageIndexType carena_index;
 497       carena_index[0] = lrint(in2[0] + in1[0])/2.0;
 498       carena_index[1] = lrint(in2[1] + in1[1])/2.0;
 499       carena_index[2] = slice_index;
 500       MaskImagePointType carena_position;
 501       DD(carena_index);
 502       skeleton->TransformIndexToPhysicalPoint(carena_index,
 503                                               carena_position);
 504       DD(carena_position);
 505
 506       // Set and save Carina position
 507       StartNewStep("Save carina position");
 508       // Try to load the DB
 509       try {
 510         LoadAFDB();
 511       } catch (clitk::ExceptionObject e) {
 512         // Do nothing if not found, it will be used anyway to write the result
 513       }
 514       GetAFDB()->SetPoint3D("Carena", carena_position);
 515     }
 516   }
 517 }
 518 //--------------------------------------------------------------------
 519
 520
 521 //--------------------------------------------------------------------
 522 template <class ImageType, class MaskImageType>
 523 void
 524 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 525 GenerateData()
 526 {
 527   // If everything goes well, set the output
 528   this->GraftOutput(output); // not SetNthOutput
 529 }
 530 //--------------------------------------------------------------------
 531
 532
 533 //--------------------------------------------------------------------
 534 template <class ImageType, class MaskImageType>
 535 void
 536 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 537 TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
 538                    MaskImagePointer skeleton,
 539                    MaskImageIndexType index,
 540                    MaskImagePixelType label,
 541                    TreeIterator currentNode)
 542 {
 543   // Create NeighborhoodIterator
 544   typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
 545   typename NeighborhoodIteratorType::SizeType radius;
 546   radius.Fill(1);
 547   NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 548
 549   // Track
 550   std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
 551   bool stop = false;
 552   while (!stop) {
 553     nit.SetLocation(index);
 554     nit.SetCenterPixel(label);
 555     listOfTrackedPoint.clear();
 556     for(unsigned int i=0; i<nit.Size(); i++) {
 557       if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
 558         if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
 559           listOfTrackedPoint.push_back(nit.GetIndex(i));
 560         }
 561       }
 562     }
 563     if (listOfTrackedPoint.size() == 1) {
 564       // Add this point to the current path
 565       currentNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]);
 566       index = listOfTrackedPoint[0];
 567     }
 568     else {
 569       if (listOfTrackedPoint.size() == 2) {
 570         // m_SkeletonTree->Add(listOfTrackedPoint[0], index); // the parent is 'index'
 571         // m_SkeletonTree->Add(listOfTrackedPoint[1], index); // the parent is 'index'
 572         BifurcationType bif(index, label, label+1, label+2);
 573         bif.treeIter = currentNode;
 574         listOfBifurcations.push_back(bif);
 575         TreeIterator firstNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]);
 576         TreeIterator secondNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[1]);
 577         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1, firstNode);
 578         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2, secondNode);
 579       }
 580       else {
 581         if (listOfTrackedPoint.size() > 2) {
 582           clitkExceptionMacro("error while tracking trachea bifurcation. Too much bifurcation points ... ?");
 583         }
 584         // Else this it the end of the tracking
 585       }
 586       stop = true;
 587     }
 588   }
 589 }
 590 //--------------------------------------------------------------------
 591
 592
 593 //--------------------------------------------------------------------
 594 template <class ImageType, class MaskImageType>
 595 bool
 596 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 597 SearchForTracheaSeed(int skip)
 598 {
 599   if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
 600     // Restart the search until a seed is found, skipping more and more slices
 601     bool stop = false;
 602     while (!stop) {
 603       // Search seed (parameters = UpperThresholdForTrachea)
 604       static const unsigned int Dim = ImageType::ImageDimension;
 605       typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
 606       typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
 607       typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
 608       sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
 609       sliceRegionSize[Dim-1]=5;
 610       sliceRegion.SetSize(sliceRegionSize);
 611       sliceRegion.SetIndex(sliceRegionIndex);
 612       typedef  itk::ImageRegionConstIterator<ImageType> IteratorType;
 613       IteratorType it(working_input, sliceRegion);
 614       it.GoToBegin();
 615       while (!it.IsAtEnd()) {
 616         if(it.Get() < GetUpperThresholdForTrachea() ) {
 617           AddSeed(it.GetIndex());
 618           // DD(it.GetIndex());
 619         }
 620         ++it;
 621       }
 622
 623       // if we do not found : restart
 624       stop = (m_Seeds.size() != 0);
 625       if (!stop) {
 626         if (GetVerboseStep()) {
 627           std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
 628         }
 629         if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
 630           // we want to skip more than a half of the image, it is probably a bug
 631           std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
 632           stop = true;
 633         }
 634         skip += 5;
 635       }
 636     }
 637   }
 638   return (m_Seeds.size() != 0);
 639 }
 640 //--------------------------------------------------------------------
 641
 642
 643 //--------------------------------------------------------------------
 644 template <class ImageType, class MaskImageType>
 645 void
 646 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 647 TracheaRegionGrowing()
 648 {
 649   // Explosion controlled region growing
 650   typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, InternalImageType> ImageFilterType;
 651   typename ImageFilterType::Pointer f= ImageFilterType::New();
 652   f->SetInput(working_input);
 653   f->SetVerbose(false);
 654   f->SetLower(-2000);
 655   f->SetUpper(GetUpperThresholdForTrachea());
 656   f->SetMinimumLowerThreshold(-2000);
 657   f->SetMaximumUpperThreshold(0);
 658   f->SetAdaptLowerBorder(false);
 659   f->SetAdaptUpperBorder(true);
 660   f->SetMinimumSize(5000);
 661   f->SetReplaceValue(1);
 662   f->SetMultiplier(GetMultiplierForTrachea());
 663   f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
 664   f->SetMinimumThresholdStepSize(1);
 665   for(unsigned int i=0; i<m_Seeds.size();i++) {
 666     f->AddSeed(m_Seeds[i]);
 667     // DD(m_Seeds[i]);
 668   }
 669   f->Update();
 670
 671   // take first (main) connected component
 672   trachea_tmp = Labelize<InternalImageType>(f->GetOutput(),
 673                                             GetBackgroundValue(),
 674                                             true,
 675                                             GetMinimalComponentSize());
 676   trachea_tmp = KeepLabels<InternalImageType>(trachea_tmp,
 677                                               GetBackgroundValue(),
 678                                               GetForegroundValue(),
 679                                               1, 1, false);
 680 }
 681 //--------------------------------------------------------------------
 682
 683
 684 //--------------------------------------------------------------------
 685 template <class ImageType, class MaskImageType>
 686 double
 687 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 688 ComputeTracheaVolume()
 689 {
 690   typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
 691   IteratorType iter(trachea_tmp, trachea_tmp->GetLargestPossibleRegion());
 692   iter.GoToBegin();
 693   double volume = 0.0;
 694   while (!iter.IsAtEnd()) {
 695     if (iter.Get() == this->GetForegroundValue()) volume++;
 696     ++iter;
 697   }
 698
 699   double voxelsize = trachea_tmp->GetSpacing()[0]*trachea_tmp->GetSpacing()[1]*trachea_tmp->GetSpacing()[2];
 700   return volume*voxelsize;
 701 }
 702 //--------------------------------------------------------------------
 703
 704
 705 //--------------------------------------------------------------------
 706 template <class ImageType, class MaskImageType>
 707 void
 708 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 709 SearchForTrachea()
 710 {
 711   // Search for seed among n slices, skip some slices before starting
 712   // if not found -> skip more and restart
 713
 714   // when seed found : perform region growing
 715   // compute trachea volume
 716   // if volume not plausible  -> skip more slices and restart
 717
 718   bool stop = false;
 719   double volume = 0.0;
 720   int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
 721   while (!stop) {
 722     stop = SearchForTracheaSeed(skip);
 723     if (stop) {
 724       TracheaRegionGrowing();
 725       volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
 726       if ((volume > 10) && (volume < 55 )) { // it is ok
 727         // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
 728         if (GetVerboseStep()) {
 729           std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
 730         }
 731         stop = true;
 732       }
 733       else {
 734         if (GetVerboseStep()) {
 735           std::cout << "\t The volume of the trachea (" << volume
 736                     << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
 737                     << std::endl;
 738         }
 739         skip += 5;
 740         stop = false;
 741         // empty the list of seed
 742         m_Seeds.clear();
 743       }
 744     }
 745   }
 746
 747   if (volume != 0.0) {
 748     // Set output
 749     StopCurrentStep<InternalImageType>(trachea_tmp);
 750   }
 751   else { // Trachea not found
 752     this->SetWarning("* WARNING * No seed found for trachea.");
 753     StopCurrentStep();
 754   }
 755 }
 756 //--------------------------------------------------------------------
 757
 758 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX