segmentation/clitkExtractAirwayTreeInfoFilter.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 CLITKEXTRACTAIRWAYTREEINFOSFILTER_TXX
  20 #define CLITKEXTRACTAIRWAYTREEINFOSFILTER_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 #include "itkBinaryThinningImageFilter3D.h"
  34 #include "itkImageIteratorWithIndex.h"
  35
  36
  37 //--------------------------------------------------------------------
  38 template <class ImageType>
  39 clitk::ExtractAirwayTreeInfoFilter<ImageType>::
  40 ExtractAirwayTreeInfoFilter():
  41   clitk::FilterBase(),
  42   itk::ImageToImageFilter<ImageType, ImageType>()
  43 {
  44   // Default global options
  45   this->SetNumberOfRequiredInputs(1);
  46   SetBackgroundValue(0);
  47   SetForegroundValue(1);
  48 }
  49 //--------------------------------------------------------------------
  50
  51
  52 //--------------------------------------------------------------------
  53 template <class ImageType>
  54 void
  55 clitk::ExtractAirwayTreeInfoFilter<ImageType>::
  56 SetInput(const ImageType * image)
  57 {
  58   this->SetNthInput(0, const_cast<ImageType *>(image));
  59 }
  60 //--------------------------------------------------------------------
  61
  62
  63 //--------------------------------------------------------------------
  64 template <class ImageType>
  65 template<class ArgsInfoType>
  66 void
  67 clitk::ExtractAirwayTreeInfoFilter<ImageType>::
  68 SetArgsInfo(ArgsInfoType mArgsInfo)
  69 {
  70   SetVerboseOption_GGO(mArgsInfo);
  71   SetVerboseStep_GGO(mArgsInfo);
  72   SetWriteStep_GGO(mArgsInfo);
  73   SetVerboseWarningOff_GGO(mArgsInfo);
  74 }
  75 //--------------------------------------------------------------------
  76
  77
  78 //--------------------------------------------------------------------
  79 template <class ImageType>
  80 void
  81 clitk::ExtractAirwayTreeInfoFilter<ImageType>::
  82 GenerateOutputInformation()
  83 {
  84   Superclass::GenerateOutputInformation();
  85   //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion());
  86
  87   // Get input pointer
  88   input   = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
  89
  90   StartNewStepOrStop("Find bronchial bifurcations");
  91   // Step 1 : extract skeleton
  92   // Define the thinning filter
  93   typedef itk::BinaryThinningImageFilter3D<ImageType, ImageType> ThinningFilterType;
  94   typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
  95   thinningFilter->SetInput(input);
  96   thinningFilter->Update();
  97   skeleton = thinningFilter->GetOutput();
  98   writeImage<ImageType>(skeleton, "skeleton.mhd");
  99
 100   // Step 2 : tracking
 101   DD("tracking");
 102
 103   // Step 2.1 : find first point for tracking
 104   typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
 105   IteratorType it(skeleton, skeleton->GetLargestPossibleRegion());
 106   it.GoToReverseBegin();
 107   while ((!it.IsAtEnd()) && (it.Get() == GetBackgroundValue())) {
 108     --it;
 109   }
 110   if (it.IsAtEnd()) {
 111     this->SetLastError("ERROR: first point in the skeleton not found ! Abort");
 112     return;
 113   }
 114   DD(skeleton->GetLargestPossibleRegion().GetIndex());
 115   typename ImageType::IndexType index = it.GetIndex();
 116   skeleton->TransformIndexToPhysicalPoint(index, m_FirstTracheaPoint);
 117   DD(index);
 118   DD(m_FirstTracheaPoint);
 119
 120   // Step 2.2 : initialize neighborhooditerator
 121   typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
 122   typename NeighborhoodIteratorType::SizeType radius;
 123   radius.Fill(1);
 124   NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 125   DD(nit.GetSize());
 126   DD(nit.Size());
 127
 128   // Find first label number (must be different from BG and FG)
 129   typename ImageType::PixelType label = GetForegroundValue()+1;
 130   while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; }
 131   DD(label);
 132
 133   // Track from the first point
 134   std::vector<BifurcationType> listOfBifurcations;
 135   TrackFromThisIndex(listOfBifurcations, skeleton, index, label);
 136   DD("end track");
 137   DD(listOfBifurcations.size());
 138   writeImage<ImageType>(skeleton, "skeleton2.mhd");
 139
 140   for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
 141     typename ImageType::PointType p;
 142     skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index, p);
 143     DD(p);
 144   }
 145
 146 }
 147 //--------------------------------------------------------------------
 148
 149
 150 //--------------------------------------------------------------------
 151 template <class ImageType>
 152 void
 153 clitk::ExtractAirwayTreeInfoFilter<ImageType>::
 154 GenerateData()
 155 {
 156   // Do not put some "startnewstep" here, because the object if
 157   // modified and the filter's pipeline it do two times. But it is
 158   // required to quit if MustStop was set before.
 159   if (GetMustStop()) return;
 160
 161   // If everything goes well, set the output
 162   this->GraftOutput(skeleton); // not SetNthOutput
 163 }
 164 //--------------------------------------------------------------------
 165
 166
 167 //--------------------------------------------------------------------
 168 template <class ImageType>
 169 void
 170 clitk::ExtractAirwayTreeInfoFilter<ImageType>::
 171 TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
 172                    ImagePointer skeleton,
 173                    ImageIndexType index,
 174                    ImagePixelType label)
 175 {
 176   DD("TrackFromThisIndex");
 177   DD(index);
 178   DD((int)label);
 179   // Create NeighborhoodIterator
 180   typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
 181   typename NeighborhoodIteratorType::SizeType radius;
 182   radius.Fill(1);
 183   NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 184
 185   // Track
 186   std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
 187   bool stop = false;
 188   while (!stop) {
 189     nit.SetLocation(index);
 190     // DD((int)nit.GetCenterPixel());
 191     nit.SetCenterPixel(label);
 192     listOfTrackedPoint.clear();
 193     for(unsigned int i=0; i<nit.Size(); i++) {
 194       if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
 195         //          DD(nit.GetIndex(i));
 196         if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
 197           // DD(nit.GetIndex(i));
 198           listOfTrackedPoint.push_back(nit.GetIndex(i));
 199         }
 200       }
 201     }
 202     // DD(listOfTrackedPoint.size());
 203     if (listOfTrackedPoint.size() == 1) {
 204       index = listOfTrackedPoint[0];
 205     }
 206     else {
 207       if (listOfTrackedPoint.size() == 2) {
 208         BifurcationType bif(index, label, label+1, label+2);
 209         listOfBifurcations.push_back(bif);
 210         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1);
 211         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2);
 212       }
 213       else {
 214         if (listOfTrackedPoint.size() > 2) {
 215           std::cerr << "too much bifurcation points ... ?" << std::endl;
 216           exit(0);
 217         }
 218         // Else this it the end of the tracking
 219       }
 220       stop = true;
 221     }
 222   }
 223 }
 224 //--------------------------------------------------------------------
 225
 226
 227 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX