segmentation/clitkExtractAirwaysTreeInfoFilter.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 CLITKEXTRACTAIRWAYSTREEINFOSFILTER_TXX
  20 #define CLITKEXTRACTAIRWAYSTREEINFOSFILTER_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 //--------------------------------------------------------------------
  39 template <class ImageType>
  40 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
  41 ExtractAirwaysTreeInfoFilter():
  42   clitk::FilterBase(),
  43   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  44   itk::ImageToImageFilter<ImageType, ImageType>()
  45 {
  46   // Default global options
  47   this->SetNumberOfRequiredInputs(1);
  48   SetBackgroundValue(0);
  49   SetForegroundValue(1);
  50 }
  51 //--------------------------------------------------------------------
  52
  53
  54 //--------------------------------------------------------------------
  55 template <class ImageType>
  56 void
  57 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
  58 SetInput(const ImageType * image)
  59 {
  60   this->SetNthInput(0, const_cast<ImageType *>(image));
  61 }
  62 //--------------------------------------------------------------------
  63
  64
  65 //--------------------------------------------------------------------
  66 template <class ImageType>
  67 template<class ArgsInfoType>
  68 void
  69 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
  70 SetArgsInfo(ArgsInfoType mArgsInfo)
  71 {
  72   SetVerboseOption_GGO(mArgsInfo);
  73   SetVerboseStep_GGO(mArgsInfo);
  74   SetWriteStep_GGO(mArgsInfo);
  75   SetVerboseWarningOff_GGO(mArgsInfo);
  76   SetAFDBFilename_GGO(mArgsInfo);
  77 }
  78 //--------------------------------------------------------------------
  79
  80
  81 //--------------------------------------------------------------------
  82 template <class ImageType>
  83 void
  84 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
  85 GenerateOutputInformation()
  86 {
  87   Superclass::GenerateOutputInformation();
  88   //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion());
  89
  90   // Get input pointer
  91   input   = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
  92
  93   // Try to load the DB
  94   try {
  95     LoadAFDB();
  96   } catch (clitk::ExceptionObject e) {
  97     // Do nothing if not found, it will be used anyway to write the result
  98   }
  99
 100
 101 }
 102 //--------------------------------------------------------------------
 103
 104
 105 //--------------------------------------------------------------------
 106 template <class ImageType>
 107 void
 108 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
 109 GenerateData()
 110 {
 111   StartNewStep("Thinning filter (skeleton)");
 112   // Extract skeleton
 113   typedef itk::BinaryThinningImageFilter3D<ImageType, ImageType> ThinningFilterType;
 114   typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
 115   thinningFilter->SetInput(input); // input = trachea
 116   thinningFilter->Update();
 117   skeleton = thinningFilter->GetOutput();
 118   StopCurrentStep<ImageType>(skeleton);
 119
 120   // Find first point for tracking
 121   StartNewStep("Find first point for tracking");
 122   typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
 123   IteratorType it(skeleton, skeleton->GetLargestPossibleRegion());
 124   it.GoToReverseBegin();
 125   while ((!it.IsAtEnd()) && (it.Get() == GetBackgroundValue())) {
 126     --it;
 127   }
 128   if (it.IsAtEnd()) {
 129     clitkExceptionMacro("first point in the trachea skeleton not found.");
 130   }
 131   typename ImageType::IndexType index = it.GetIndex();
 132   DD(index);
 133
 134   // Initialize neighborhooditerator
 135   typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
 136   typename NeighborhoodIteratorType::SizeType radius;
 137   radius.Fill(1);
 138   NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 139
 140   // Find first label number (must be different from BG and FG)
 141   typename ImageType::PixelType label = GetForegroundValue()+1;
 142   while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; }
 143   DD(label);
 144
 145   /*
 146     Tracking ?
 147     Goal : find position of C, RUL, RML, RLL, LUL, LLL bronchus
 148     Carina : ok "easy", track, slice by slice until 2 path into different label
 149     -> follow at Right
 150        -
 151     -> follow at Left
 152    */
 153
 154
 155   // Track from the first point
 156   StartNewStep("Start tracking");
 157   std::vector<BifurcationType> listOfBifurcations;
 158   m_SkeletonTree.set_head(index);
 159   TrackFromThisIndex(listOfBifurcations, skeleton, index, label, m_SkeletonTree.begin());
 160   DD("end track");
 161
 162   // Convert index into physical point coordinates
 163   for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
 164     skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index,
 165                                             listOfBifurcations[i].point);
 166   }
 167
 168   // Search for the first slice that separate the bronchus
 169   // (carina). Labelize slice by slice, stop when the two points of
 170   // the skeleton ar not in the same connected component
 171   StartNewStep("Search for carina position");
 172   typedef clitk::ExtractSliceFilter<ImageType> ExtractSliceFilterType;
 173   typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
 174   extractSliceFilter->SetInput(input);
 175   extractSliceFilter->SetDirection(2);
 176   extractSliceFilter->Update();
 177   typedef typename ExtractSliceFilterType::SliceType SliceType;
 178   std::vector<typename SliceType::Pointer> mInputSlices;
 179   extractSliceFilter->GetOutputSlices(mInputSlices);
 180   DD(mInputSlices.size());
 181
 182   bool stop = false;
 183   int slice_index = listOfBifurcations[0].index[2]; // first slice from carina in skeleton
 184   int i=0;
 185   TreeIterator firstIter = m_SkeletonTree.child(listOfBifurcations[0].treeIter, 0);
 186   TreeIterator secondIter = m_SkeletonTree.child(listOfBifurcations[0].treeIter, 1);
 187   typename SliceType::IndexType in1;
 188   typename SliceType::IndexType in2;
 189   while (!stop) {
 190     //  Labelize the current slice
 191     typename SliceType::Pointer temp = Labelize<SliceType>(mInputSlices[slice_index],
 192                                                            GetBackgroundValue(),
 193                                                            true,
 194                                                            0); // min component size=0
 195     // Check the value of the two skeleton points;
 196     in1[0] = (*firstIter)[0];
 197     in1[1] = (*firstIter)[1];
 198     typename SliceType::PixelType v1 = temp->GetPixel(in1);
 199     in2[0] = (*secondIter)[0];
 200     in2[1] = (*secondIter)[1];
 201     typename SliceType::PixelType v2 = temp->GetPixel(in2);
 202
 203     // Check the label value of the two points
 204     DD(slice_index);
 205     if (v1 != v2) {
 206       stop = true; // We found it !
 207     }
 208     else {
 209       // Check error
 210       if (slice_index == (int)(mInputSlices.size()-1)) {
 211         clitkExceptionMacro("Error while searching for carina, the two skeleton points are always in the same CC ... ???");
 212       }
 213       // Iterate
 214       i++;
 215       --slice_index;
 216       ++firstIter;
 217       ++secondIter;
 218     }
 219   }
 220   ImageIndexType carina_index; // middle position in X/Y
 221   carina_index[0] = lrint(in2[0] + in1[0])/2.0;
 222   carina_index[1] = lrint(in2[1] + in1[1])/2.0;
 223   carina_index[2] = slice_index;
 224   // Get physical coordinates
 225   ImagePointType carina_position;
 226   skeleton->TransformIndexToPhysicalPoint(carina_index,
 227                                           carina_position);
 228
 229   // Set and save Carina position
 230   if (GetVerboseStep()) {
 231     std::cout << "\t Found carina at " << carina_position << " mm" << std::endl;
 232   }
 233   GetAFDB()->SetPoint3D("Carina", carina_position);
 234
 235   // Write bifurcation (debug)
 236   for(uint i=0; i<listOfBifurcations.size(); i++) {
 237     GetAFDB()->SetPoint3D("Bif"+toString(i), listOfBifurcations[i].point);
 238   }
 239
 240   // Set the output (skeleton);
 241   this->GraftOutput(skeleton); // not SetNthOutput
 242 }
 243 //--------------------------------------------------------------------
 244
 245
 246 //--------------------------------------------------------------------
 247 template <class ImageType>
 248 void
 249 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
 250 TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
 251                    ImagePointer skeleton,
 252                    ImageIndexType index,
 253                    ImagePixelType label,
 254                    TreeIterator currentNode)
 255 {
 256   // Create NeighborhoodIterator
 257   typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
 258   typename NeighborhoodIteratorType::SizeType radius;
 259   radius.Fill(1);
 260   NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 261
 262   // Track
 263   std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
 264   bool stop = false;
 265   while (!stop) {
 266     nit.SetLocation(index);
 267     nit.SetCenterPixel(label);
 268     listOfTrackedPoint.clear();
 269     for(unsigned int i=0; i<nit.Size(); i++) {
 270       if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
 271         if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
 272           listOfTrackedPoint.push_back(nit.GetIndex(i));
 273         }
 274       }
 275     }
 276     if (listOfTrackedPoint.size() == 1) {
 277       // Add this point to the current path
 278       currentNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]);
 279       index = listOfTrackedPoint[0];
 280       skeleton->SetPixel(index, label); // change label in skeleton image
 281     }
 282     else {
 283       if (listOfTrackedPoint.size() == 2) {
 284         // m_SkeletonTree->Add(listOfTrackedPoint[0], index); // the parent is 'index'
 285         // m_SkeletonTree->Add(listOfTrackedPoint[1], index); // the parent is 'index'
 286         BifurcationType bif(index, label, label+1, label+2);
 287         bif.treeIter = currentNode;
 288         listOfBifurcations.push_back(bif);
 289         TreeIterator firstNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]);
 290         TreeIterator secondNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[1]);
 291         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1, firstNode);
 292         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2, secondNode);
 293       }
 294       else {
 295         //DD(listOfTrackedPoint.size());
 296         if (listOfTrackedPoint.size() > 2) {
 297           //clitkExceptionMacro("error while tracking trachea bifurcation. Too much bifurcation points ... ?");
 298           stop = true; // this it the end of the tracking
 299         }
 300         // Else this it the end of the tracking
 301       }
 302       stop = true;
 303     }
 304   }
 305 }
 306 //--------------------------------------------------------------------
 307
 308 /*TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
 309   ImagePointer skeleton,
 310   ImageIndexType index,
 311   ImagePixelType label)
 312   {
 313   DD("TrackFromThisIndex");
 314   DD(index);
 315   DD((int)label);
 316   // Create NeighborhoodIterator
 317   typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
 318   typename NeighborhoodIteratorType::SizeType radius;
 319   radius.Fill(1);
 320   NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 321
 322   // Track
 323   std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
 324   bool stop = false;
 325   while (!stop) {
 326   nit.SetLocation(index);
 327   // DD((int)nit.GetCenterPixel());
 328   nit.SetCenterPixel(label);
 329   listOfTrackedPoint.clear();
 330   for(unsigned int i=0; i<nit.Size(); i++) {
 331   if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
 332   //          DD(nit.GetIndex(i));
 333   if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
 334   // DD(nit.GetIndex(i));
 335   listOfTrackedPoint.push_back(nit.GetIndex(i));
 336   }
 337   }
 338   }
 339   // DD(listOfTrackedPoint.size());
 340   if (listOfTrackedPoint.size() == 1) {
 341   index = listOfTrackedPoint[0];
 342   }
 343   else {
 344   if (listOfTrackedPoint.size() == 2) {
 345   BifurcationType bif(index, label, label+1, label+2);
 346   listOfBifurcations.push_back(bif);
 347   TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1);
 348   TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2);
 349   }
 350   else {
 351   if (listOfTrackedPoint.size() > 2) {
 352   std::cerr << "too much bifurcation points ... ?" << std::endl;
 353   exit(0);
 354   }
 355   // Else this it the end of the tracking
 356   }
 357   stop = true;
 358   }
 359   }
 360   }
 361 */
 362 //--------------------------------------------------------------------
 363
 364
 365 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX