X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLungFilter.txx;h=52276f4478bb0cd55a6040efd07a763bd291d120;hb=df51e4850f784f8f2fbaead7bb82bd6897fb9fab;hp=92b2d4f75cf2e8daa0c6633e64c23c7812886ffd;hpb=090511076dfb931319d0ec16fa5a08f2467362b6;p=clitk.git

diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx
index 92b2d4f..52276f4 100644
--- a/segmentation/clitkExtractLungFilter.txx
+++ b/segmentation/clitkExtractLungFilter.txx
@@ -24,7 +24,8 @@
 #include "clitkSetBackgroundImageFilter.h"
 #include "clitkSegmentationUtils.h"
 #include "clitkAutoCropFilter.h"
-#include "clitkExtractSliceFilter.h"
+#include "clitkCropLikeImageFilter.h"
+#include "clitkFillMaskFilter.h"
 
 // itk
 #include "itkBinaryThresholdImageFilter.h"
@@ -33,10 +34,12 @@
 #include "itkOtsuThresholdImageFilter.h"
 #include "itkBinaryThinningImageFilter3D.h"
 #include "itkImageIteratorWithIndex.h"
+#include "itkBinaryMorphologicalOpeningImageFilter.h"
+#include "itkBinaryMorphologicalClosingImageFilter.h"
 
 //--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+template <class ImageType>
+clitk::ExtractLungFilter<ImageType>::
 ExtractLungFilter():
   clitk::FilterBase(),
   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
@@ -46,7 +49,7 @@ ExtractLungFilter():
   m_MaxSeedNumber = 500;
 
   // Default global options
-  this->SetNumberOfRequiredInputs(2);
+  this->SetNumberOfRequiredInputs(1);
   SetPatientMaskBackgroundValue(0);
   SetBackgroundValue(0); // Must be zero
   SetForegroundValue(1);
@@ -84,16 +87,20 @@ ExtractLungFilter():
   p3->UseLastKeepOff();
   SetLabelizeParameters3(p3);
   
-  // Step 5 : find bronchial bifurcations
-  FindBronchialBifurcationsOn();
+  // Step 5
+  OpenCloseOff();
+  SetOpenCloseRadius(1);
+  
+  // Step 6
+  FillHolesOn();
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
 void 
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
 SetInput(const ImageType * image) 
 {
   this->SetNthInput(0, const_cast<ImageType *>(image));
@@ -102,21 +109,9 @@ SetInput(const ImageType * image)
 
 
 //--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
 void 
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
-SetInputPatientMask(MaskImageType * image, MaskImagePixelType bg ) 
-{
-  this->SetNthInput(1, const_cast<MaskImageType *>(image));
-  SetPatientMaskBackgroundValue(bg);
-}
-//--------------------------------------------------------------------
-
-
-//--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
-void 
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
 AddSeed(InternalIndexType s) 
 { 
   m_Seeds.push_back(s);
@@ -125,10 +120,10 @@ AddSeed(InternalIndexType s)
 
 
 //--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
 template<class ArgsInfoType>
 void 
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
 SetArgsInfo(ArgsInfoType mArgsInfo)
 {
   SetVerboseOption_GGO(mArgsInfo);
@@ -136,6 +131,10 @@ SetArgsInfo(ArgsInfoType mArgsInfo)
   SetWriteStep_GGO(mArgsInfo);
   SetVerboseWarningOff_GGO(mArgsInfo);
 
+  SetAFDBFilename_GGO(mArgsInfo);
+  SetOutputLungFilename_GGO(mArgsInfo);
+  SetOutputTracheaFilename_GGO(mArgsInfo);
+
   SetUpperThreshold_GGO(mArgsInfo);
   SetLowerThreshold_GGO(mArgsInfo);
   SetNumberOfSlicesToSkipBeforeSearchingSeed_GGO(mArgsInfo);
@@ -158,34 +157,47 @@ SetArgsInfo(ArgsInfoType mArgsInfo)
   SetRadiusForTrachea_GGO(mArgsInfo);
   SetLabelizeParameters3_GGO(mArgsInfo);
   
-  SetAFDBFilename_GGO(mArgsInfo);
+  SetOpenCloseRadius_GGO(mArgsInfo);
+  SetOpenClose_GGO(mArgsInfo);
+  
+  SetFillHoles_GGO(mArgsInfo);
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
 void 
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
 GenerateOutputInformation() 
 { 
   Superclass::GenerateOutputInformation();
-  //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion());
+  
+  // Read DB
+  LoadAFDB();
 
   // Get input pointers
-  patient = dynamic_cast<const MaskImageType*>(itk::ProcessObject::GetInput(1));
   input   = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
+  patient = GetAFDB()->template GetImage <MaskImageType>("patient");  
 
-  // Check image
-  if (!HaveSameSizeAndSpacing<ImageType, MaskImageType>(input, patient)) {
-    clitkExceptionMacro("the 'input' and 'patient' masks must have the same size & spacing.");
-  }
-  
+  //--------------------------------------------------------------------
+  //--------------------------------------------------------------------
+  // Crop input like patient image (must have the same spacing)
+  StartNewStep("Crop input image to 'patient' extends");
+  typedef clitk::CropLikeImageFilter<ImageType> CropImageFilter;
+  typename CropImageFilter::Pointer cropFilter = CropImageFilter::New();
+  cropFilter->SetInput(input);
+  cropFilter->SetCropLikeImage(patient);
+  cropFilter->Update();
+  working_input = cropFilter->GetOutput();
+  DD(working_input->GetLargestPossibleRegion());
+  StopCurrentStep<ImageType>(working_input);
+ 
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
   StartNewStep("Set background to initial image");
   working_input = SetBackground<ImageType, MaskImageType>
-    (input, patient, GetPatientMaskBackgroundValue(), -1000);
+    (working_input, patient, GetPatientMaskBackgroundValue(), -1000);
   StopCurrentStep<ImageType>(working_input);
 
   //--------------------------------------------------------------------
@@ -312,15 +324,15 @@ GenerateOutputInformation()
 
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
-  typedef clitk::AutoCropFilter<InternalImageType> CropFilterType;
-  typename CropFilterType::Pointer cropFilter = CropFilterType::New();
+  typedef clitk::AutoCropFilter<InternalImageType> AutoCropFilterType;
+  typename AutoCropFilterType::Pointer autocropFilter = AutoCropFilterType::New();
   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
-    StartNewStep("Croping trachea");
-    cropFilter->SetInput(trachea_tmp);
-    cropFilter->Update(); // Needed
+    StartNewStep("Cropping trachea");
+    autocropFilter->SetInput(trachea_tmp);
+    autocropFilter->Update(); // Needed
     typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
     typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
-    caster->SetInput(cropFilter->GetOutput());
+    caster->SetInput(autocropFilter->GetOutput());
     caster->Update();   
     trachea = caster->GetOutput();
     StopCurrentStep<MaskImageType>(trachea);  
@@ -328,13 +340,60 @@ GenerateOutputInformation()
 
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
-  StartNewStep("Croping lung");
-  typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline
-  cropFilter2->SetInput(working_image);
-  cropFilter2->Update();   
-  working_image = cropFilter2->GetOutput();
+  StartNewStep("Cropping lung");
+  typename AutoCropFilterType::Pointer autocropFilter2 = AutoCropFilterType::New(); // Needed to reset pipeline
+  autocropFilter2->SetInput(working_image);
+  autocropFilter2->Update();   
+  working_image = autocropFilter2->GetOutput();
+  DD(working_image->GetLargestPossibleRegion());
   StopCurrentStep<InternalImageType>(working_image);
 
+  //--------------------------------------------------------------------
+  //--------------------------------------------------------------------
+  // Final OpenClose
+  if (GetOpenClose()) {
+    StartNewStep("Open/Close"); 
+
+    // Structuring element
+    typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
+    KernelType structuringElement;
+    structuringElement.SetRadius(GetOpenCloseRadius());
+    structuringElement.CreateStructuringElement();
+	
+    // Open
+    typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType, KernelType> OpenFilterType;
+    typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
+    openFilter->SetInput(working_image);
+    openFilter->SetBackgroundValue(GetBackgroundValue());
+    openFilter->SetForegroundValue(GetForegroundValue());
+    openFilter->SetKernel(structuringElement);
+	
+    // Close
+    typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType, KernelType> CloseFilterType;
+    typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
+    closeFilter->SetInput(openFilter->GetOutput());
+    closeFilter->SetSafeBorder(true);
+    closeFilter->SetForegroundValue(GetForegroundValue());
+    closeFilter->SetKernel(structuringElement);
+    closeFilter->Update();
+    working_image = closeFilter->GetOutput();
+  }
+
+  //--------------------------------------------------------------------
+  //--------------------------------------------------------------------
+  // Fill Lungs
+  if (GetFillHoles()) {
+    StartNewStep("Fill Holes");
+    typedef clitk::FillMaskFilter<InternalImageType> FillMaskFilterType;
+    typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New();
+    fillMaskFilter->SetInput(working_image);
+    fillMaskFilter->AddDirection(2);
+    //fillMaskFilter->AddDirection(1);
+    fillMaskFilter->Update();   
+    working_image = fillMaskFilter->GetOutput();
+    StopCurrentStep<InternalImageType>(working_image);
+  }
+
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
   StartNewStep("Separate Left/Right lungs");
@@ -355,13 +414,14 @@ GenerateOutputInformation()
   // Decompose the first label
   static const unsigned int Dim = ImageType::ImageDimension;
   if (initialNumberOfLabels<2) {
+    DD(initialNumberOfLabels);
     // Structuring element radius
     typename ImageType::SizeType radius;
     for (unsigned int i=0;i<Dim;i++) radius[i]=1;
     typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> DecomposeAndReconstructFilterType;
     typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
     decomposeAndReconstructFilter->SetInput(working_image);
-    decomposeAndReconstructFilter->SetVerbose(false);
+    decomposeAndReconstructFilter->SetVerbose(true);
     decomposeAndReconstructFilter->SetRadius(radius);
     decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
     decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
@@ -394,206 +454,30 @@ GenerateOutputInformation()
 
   // Update output info
   this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
-
-  // Try to extract bifurcation in the trachea (bronchi)
-  if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
-
-    if (GetFindBronchialBifurcations()) {
-      StartNewStep("Find bronchial bifurcations");
-      // Step 1 : extract skeleton
-      typedef itk::BinaryThinningImageFilter3D<MaskImageType, MaskImageType> ThinningFilterType;
-      typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
-      thinningFilter->SetInput(trachea);
-      thinningFilter->Update();
-      typename MaskImageType::Pointer skeleton = thinningFilter->GetOutput();
-
-      // Step 2.1 : find first point for tracking
-      typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> IteratorType;
-      IteratorType it(skeleton, skeleton->GetLargestPossibleRegion());
-      it.GoToReverseBegin();
-      while ((!it.IsAtEnd()) && (it.Get() == GetBackgroundValue())) { 
-        --it;
-      }
-      if (it.IsAtEnd()) {
-        clitkExceptionMacro("first point in the trachea skeleton not found.");
-      }
-      typename MaskImageType::IndexType index = it.GetIndex();
-    
-      // Step 2.2 : initialize neighborhooditerator
-      typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
-      typename NeighborhoodIteratorType::SizeType radius;
-      radius.Fill(1);
-      NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
-    
-      // Find first label number (must be different from BG and FG)
-      typename MaskImageType::PixelType label = GetForegroundValue()+1;
-      while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; }
-
-      // Track from the first point
-      std::vector<BifurcationType> listOfBifurcations;
-      m_SkeletonTree.set_head(index);
-      TrackFromThisIndex(listOfBifurcations, skeleton, index, label, m_SkeletonTree.begin());
-      DD("end track");
-      DD(listOfBifurcations.size());
-      DD(m_SkeletonTree.size());
-      
-      for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
-        skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index, 
-                                                listOfBifurcations[i].point);
-      }
-
-      // Search for the first slice that separate the bronchus (carena)
-      typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
-      typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
-      extractSliceFilter->SetInput(trachea);
-      extractSliceFilter->SetDirection(2);
-      extractSliceFilter->Update();
-      typedef typename ExtractSliceFilterType::SliceType SliceType;
-      std::vector<typename SliceType::Pointer> mInputSlices;
-      extractSliceFilter->GetOutputSlices(mInputSlices);
-      
-      bool stop = false;
-      DD(listOfBifurcations[0].index);
-      DD(listOfBifurcations[1].index);
-      int slice_index = listOfBifurcations[0].index[2]; // first slice from carena in skeleton
-      int i=0;
-      TreeIterator firstIter = m_SkeletonTree.child(listOfBifurcations[1].treeIter, 0);
-      TreeIterator secondIter = m_SkeletonTree.child(listOfBifurcations[1].treeIter, 1);
-      typename SliceType::IndexType in1;
-      typename SliceType::IndexType in2;
-      while (!stop) {
-        DD(slice_index);
-
-	//  Labelize the current slice
-        typename SliceType::Pointer temp = Labelize<SliceType>(mInputSlices[slice_index],
-                                                               GetBackgroundValue(), 
-                                                               true, 
-                                                               GetMinimalComponentSize());
-	// Check the value of the two skeleton points;
-        in1[0] = (*firstIter)[0];
-        in1[1] = (*firstIter)[1];
-	typename SliceType::PixelType v1 = temp->GetPixel(in1);
-	DD(in1);
-	DD((int)v1);
-        in2[0] = (*secondIter)[0];
-        in2[1] = (*secondIter)[1];
-	typename SliceType::PixelType v2 = temp->GetPixel(in2);
-	DD(in2);
-	DD((int)v2);
-
-	// TODO IF NOT FOUND ????
-
-        if (v1 != v2) {
-	  stop = true;
-	}
-	else {
-	  i++;
-	  --slice_index;
-	  ++firstIter;
-	  ++secondIter;
-	}
-      }
-      MaskImageIndexType carena_index;
-      carena_index[0] = lrint(in2[0] + in1[0])/2.0;
-      carena_index[1] = lrint(in2[1] + in1[1])/2.0;
-      carena_index[2] = slice_index;
-      MaskImagePointType carena_position;
-      DD(carena_index);
-      skeleton->TransformIndexToPhysicalPoint(carena_index,
-					      carena_position);
-      DD(carena_position);
-
-      // Set and save Carina position      
-      StartNewStep("Save carina position");
-      // Try to load the DB
-      try {
-        LoadAFDB();
-      } catch (clitk::ExceptionObject e) {
-        // Do nothing if not found, it will be used anyway to write the result
-      }
-      GetAFDB()->SetPoint3D("Carena", carena_position);
-    }
-  }
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
 void 
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
 GenerateData() 
 {
-  // If everything goes well, set the output
+  // Set the output
   this->GraftOutput(output); // not SetNthOutput
+  // Store image filenames into AFDB 
+  GetAFDB()->SetImageFilename("lungs", this->GetOutputLungFilename());  
+  GetAFDB()->SetImageFilename("trachea", this->GetOutputTracheaFilename());  
+  WriteAFDB();
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
-void 
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
-TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations, 
-                   MaskImagePointer skeleton, 
-                   MaskImageIndexType index,
-                   MaskImagePixelType label, 
-		   TreeIterator currentNode) 
-{
-  // Create NeighborhoodIterator 
-  typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
-  typename NeighborhoodIteratorType::SizeType radius;
-  radius.Fill(1);
-  NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
-      
-  // Track
-  std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
-  bool stop = false;
-  while (!stop) {
-    nit.SetLocation(index);
-    nit.SetCenterPixel(label);
-    listOfTrackedPoint.clear();
-    for(unsigned int i=0; i<nit.Size(); i++) {
-      if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
-        if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
-          listOfTrackedPoint.push_back(nit.GetIndex(i));
-        }
-      }
-    }
-    if (listOfTrackedPoint.size() == 1) {
-      // Add this point to the current path
-      currentNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]);
-      index = listOfTrackedPoint[0];
-    }
-    else {
-      if (listOfTrackedPoint.size() == 2) {
-        // m_SkeletonTree->Add(listOfTrackedPoint[0], index); // the parent is 'index'
-        // m_SkeletonTree->Add(listOfTrackedPoint[1], index); // the parent is 'index'
-        BifurcationType bif(index, label, label+1, label+2);
-	bif.treeIter = currentNode;
-        listOfBifurcations.push_back(bif);
-	TreeIterator firstNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]);
-        TreeIterator secondNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[1]);
-        TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1, firstNode);
-        TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2, secondNode);
-      }
-      else {
-        if (listOfTrackedPoint.size() > 2) {
-          clitkExceptionMacro("error while tracking trachea bifurcation. Too much bifurcation points ... ?");
-        }
-        // Else this it the end of the tracking
-      }
-      stop = true;
-    }
-  }
-}
-//--------------------------------------------------------------------
-
-
-//--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
 bool 
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
 SearchForTracheaSeed(int skip)
 {
   if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)    
@@ -641,9 +525,9 @@ SearchForTracheaSeed(int skip)
 
   
 //--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
 void 
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
 TracheaRegionGrowing()
 {
   // Explosion controlled region growing
@@ -654,7 +538,8 @@ TracheaRegionGrowing()
   f->SetLower(-2000);
   f->SetUpper(GetUpperThresholdForTrachea());
   f->SetMinimumLowerThreshold(-2000);
-  f->SetMaximumUpperThreshold(0);
+  //  f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
+  f->SetMaximumUpperThreshold(-800); // MAYBE TO CHANGE ???
   f->SetAdaptLowerBorder(false);
   f->SetAdaptUpperBorder(true);
   f->SetMinimumSize(5000); 
@@ -662,12 +547,15 @@ TracheaRegionGrowing()
   f->SetMultiplier(GetMultiplierForTrachea());
   f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
   f->SetMinimumThresholdStepSize(1);
+  f->VerboseOn();
   for(unsigned int i=0; i<m_Seeds.size();i++) {
     f->AddSeed(m_Seeds[i]);
     // DD(m_Seeds[i]);
   }  
   f->Update();
 
+  writeImage<InternalImageType>(f->GetOutput(), "trg.mhd");
+
   // take first (main) connected component
   trachea_tmp = Labelize<InternalImageType>(f->GetOutput(), 
 					    GetBackgroundValue(), 
@@ -682,9 +570,9 @@ TracheaRegionGrowing()
 
 
 //--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
 double 
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
 ComputeTracheaVolume()
 {
   typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
@@ -703,9 +591,9 @@ ComputeTracheaVolume()
 
 
 //--------------------------------------------------------------------
-template <class ImageType, class MaskImageType>
+template <class ImageType>
 void 
-clitk::ExtractLungFilter<ImageType, MaskImageType>::
+clitk::ExtractLungFilter<ImageType>::
 SearchForTrachea()
 {
   // Search for seed among n slices, skip some slices before starting