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
 183   DD("REDO !!!!!!!!!!!!");
 184   /**
 185      => chercher la bif qui a les plus important sous-arbres
 186    **/
 187
 188   bool stop = false;
 189   int slice_index = listOfBifurcations[0].index[2]; // first slice from carina in skeleton
 190   int i=0;
 191   TreeIterator firstIter = m_SkeletonTree.child(listOfBifurcations[0].treeIter, 0);
 192   TreeIterator secondIter = m_SkeletonTree.child(listOfBifurcations[0].treeIter, 1);
 193   DD(firstIter.number_of_children());
 194   DD(secondIter.number_of_children());
 195   typename SliceType::IndexType in1;
 196   typename SliceType::IndexType in2;
 197   while (!stop) {
 198     //  Labelize the current slice
 199     typename SliceType::Pointer temp = Labelize<SliceType>(mInputSlices[slice_index],
 200                                                            GetBackgroundValue(),
 201                                                            true,
 202                                                            0); // min component size=0
 203     DD(*firstIter);
 204     DD(*secondIter);
 205     // Check the value of the two skeleton points;
 206     in1[0] = (*firstIter)[0];
 207     in1[1] = (*firstIter)[1];
 208     typename SliceType::PixelType v1 = temp->GetPixel(in1);
 209     in2[0] = (*secondIter)[0];
 210     in2[1] = (*secondIter)[1];
 211     typename SliceType::PixelType v2 = temp->GetPixel(in2);
 212
 213     // Check the label value of the two points
 214     DD(slice_index);
 215     if (v1 != v2) {
 216       stop = true; // We found it !
 217     }
 218     else {
 219       // Check error
 220       if (slice_index == (int)(mInputSlices.size()-1)) {
 221         clitkExceptionMacro("Error while searching for carina, the two skeleton points are always in the same CC ... ???");
 222       }
 223       // Iterate
 224       i++;
 225       --slice_index;
 226       ++firstIter;
 227       ++secondIter;
 228     }
 229   }
 230   ImageIndexType carina_index; // middle position in X/Y
 231   carina_index[0] = lrint(in2[0] + in1[0])/2.0;
 232   carina_index[1] = lrint(in2[1] + in1[1])/2.0;
 233   carina_index[2] = slice_index;
 234   // Get physical coordinates
 235   ImagePointType carina_position;
 236   skeleton->TransformIndexToPhysicalPoint(carina_index,
 237                                           carina_position);
 238
 239   // Set and save Carina position
 240   if (GetVerboseStep()) {
 241     std::cout << "\t Found carina at " << carina_position << " mm" << std::endl;
 242   }
 243   GetAFDB()->SetPoint3D("Carina", carina_position);
 244
 245   // Write bifurcation (debug)
 246   for(uint i=0; i<listOfBifurcations.size(); i++) {
 247     GetAFDB()->SetPoint3D("Bif"+toString(i), listOfBifurcations[i].point);
 248   }
 249
 250   // Set the output (skeleton);
 251   this->GraftOutput(skeleton); // not SetNthOutput
 252 }
 253 //--------------------------------------------------------------------
 254
 255
 256 //--------------------------------------------------------------------
 257 template <class ImageType>
 258 void
 259 clitk::ExtractAirwaysTreeInfoFilter<ImageType>::
 260 TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
 261                    ImagePointer skeleton,
 262                    ImageIndexType index,
 263                    ImagePixelType label,
 264                    TreeIterator currentNode)
 265 {
 266   // Create NeighborhoodIterator
 267   typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
 268   typename NeighborhoodIteratorType::SizeType radius;
 269   radius.Fill(1);
 270   NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 271
 272   // Track
 273   std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
 274   bool stop = false;
 275   while (!stop) {
 276     nit.SetLocation(index);
 277     nit.SetCenterPixel(label);
 278     listOfTrackedPoint.clear();
 279     for(unsigned int i=0; i<nit.Size(); i++) {
 280       if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
 281         if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
 282           listOfTrackedPoint.push_back(nit.GetIndex(i));
 283         }
 284       }
 285     }
 286     if (listOfTrackedPoint.size() == 1) {
 287       // Add this point to the current path
 288       currentNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]);
 289       index = listOfTrackedPoint[0];
 290       skeleton->SetPixel(index, label); // change label in skeleton image
 291     }
 292     else {
 293       if (listOfTrackedPoint.size() == 2) {
 294         // m_SkeletonTree->Add(listOfTrackedPoint[0], index); // the parent is 'index'
 295         // m_SkeletonTree->Add(listOfTrackedPoint[1], index); // the parent is 'index'
 296         DD("BifurcationType");
 297         DD(listOfTrackedPoint[0]);
 298         DD(listOfTrackedPoint[1]);
 299         BifurcationType bif(index, label, label+1, label+2);
 300         bif.treeIter = currentNode;
 301         listOfBifurcations.push_back(bif);
 302         TreeIterator firstNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]);
 303         TreeIterator secondNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[1]);
 304         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1, firstNode);
 305         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2, secondNode);
 306       }
 307       else {
 308         //DD(listOfTrackedPoint.size());
 309         if (listOfTrackedPoint.size() > 2) {
 310           //clitkExceptionMacro("error while tracking trachea bifurcation. Too much bifurcation points ... ?");
 311           stop = true; // this it the end of the tracking
 312         }
 313         // Else this it the end of the tracking
 314       }
 315       stop = true;
 316     }
 317   }
 318 }
 319 //--------------------------------------------------------------------
 320
 321 /*TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
 322   ImagePointer skeleton,
 323   ImageIndexType index,
 324   ImagePixelType label)
 325   {
 326   DD("TrackFromThisIndex");
 327   DD(index);
 328   DD((int)label);
 329   // Create NeighborhoodIterator
 330   typedef itk::NeighborhoodIterator<ImageType> NeighborhoodIteratorType;
 331   typename NeighborhoodIteratorType::SizeType radius;
 332   radius.Fill(1);
 333   NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 334
 335   // Track
 336   std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
 337   bool stop = false;
 338   while (!stop) {
 339   nit.SetLocation(index);
 340   // DD((int)nit.GetCenterPixel());
 341   nit.SetCenterPixel(label);
 342   listOfTrackedPoint.clear();
 343   for(unsigned int i=0; i<nit.Size(); i++) {
 344   if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
 345   //          DD(nit.GetIndex(i));
 346   if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
 347   // DD(nit.GetIndex(i));
 348   listOfTrackedPoint.push_back(nit.GetIndex(i));
 349   }
 350   }
 351   }
 352   // DD(listOfTrackedPoint.size());
 353   if (listOfTrackedPoint.size() == 1) {
 354   index = listOfTrackedPoint[0];
 355   }
 356   else {
 357   if (listOfTrackedPoint.size() == 2) {
 358   BifurcationType bif(index, label, label+1, label+2);
 359   listOfBifurcations.push_back(bif);
 360   TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1);
 361   TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2);
 362   }
 363   else {
 364   if (listOfTrackedPoint.size() > 2) {
 365   std::cerr << "too much bifurcation points ... ?" << std::endl;
 366   exit(0);
 367   }
 368   // Else this it the end of the tracking
 369   }
 370   stop = true;
 371   }
 372   }
 373   }
 374 */
 375 //--------------------------------------------------------------------
 376
 377
 378 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX