From: dsarrut <dsarrut>
Date: Tue, 5 Oct 2010 09:31:03 +0000 (+0000)
Subject: separate airway tracking from extract lung
X-Git-Tag: v1.2.0~360
X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=44ebb81a500e42cf19964d209e14a59812a0dd4d;p=clitk.git

separate airway tracking from extract lung
---

diff --git a/segmentation/CMakeLists.txt b/segmentation/CMakeLists.txt
index ab75b33..f97e7e5 100644
--- a/segmentation/CMakeLists.txt
+++ b/segmentation/CMakeLists.txt
@@ -28,9 +28,9 @@ IF(CLITK_BUILD_SEGMENTATION)
     ADD_EXECUTABLE(clitkExtractLung clitkExtractLung.cxx  ${clitkExtractLung_GGO_C})
     TARGET_LINK_LIBRARIES(clitkExtractLung clitkSegmentationGgoLib clitkCommon ITKIO)
 
-    WRAP_GGO(clitkExtractAirwayTreeInfo_GGO_C clitkExtractAirwayTreeInfo.ggo)
-    ADD_EXECUTABLE(clitkExtractAirwayTreeInfo clitkExtractAirwayTreeInfo.cxx ${clitkExtractAirwayTreeInfo_GGO_C})
-    TARGET_LINK_LIBRARIES(clitkExtractAirwayTreeInfo clitkCommon ITKIO)
+    WRAP_GGO(clitkExtractAirwaysTreeInfo_GGO_C clitkExtractAirwaysTreeInfo.ggo)
+    ADD_EXECUTABLE(clitkExtractAirwaysTreeInfo clitkExtractAirwaysTreeInfo.cxx ${clitkExtractAirwaysTreeInfo_GGO_C})
+    TARGET_LINK_LIBRARIES(clitkExtractAirwaysTreeInfo clitkSegmentationGgoLib clitkCommon ITKIO)
 
     WRAP_GGO(clitkExtractBones_GGO_C clitkExtractBones.ggo)
     ADD_EXECUTABLE(clitkExtractBones clitkExtractBones.cxx ${clitkExtractBones_GGO_C})
@@ -38,7 +38,7 @@ IF(CLITK_BUILD_SEGMENTATION)
 
     WRAP_GGO(clitkExtractMediastinum_GGO_C clitkExtractMediastinum.ggo)
     ADD_EXECUTABLE(clitkExtractMediastinum clitkExtractMediastinum.cxx ${clitkExtractMediastinum_GGO_C})
-    TARGET_LINK_LIBRARIES(clitkExtractMediastinum clitkCommon ITKIO )
+    TARGET_LINK_LIBRARIES(clitkExtractMediastinum clitkCommon clitkSegmentationGgoLib ITKIO)
 
     WRAP_GGO(clitkExtractLymphStations_GGO_C clitkExtractLymphStations.ggo)
     ADD_EXECUTABLE(clitkExtractLymphStations clitkExtractLymphStations.cxx clitkFilterWithAnatomicalFeatureDatabaseManagement.cxx ${clitkExtractLymphStations_GGO_C})
diff --git a/segmentation/clitkAnatomicalFeatureDatabase.cxx b/segmentation/clitkAnatomicalFeatureDatabase.cxx
index a132f58..afd082e 100644
--- a/segmentation/clitkAnatomicalFeatureDatabase.cxx
+++ b/segmentation/clitkAnatomicalFeatureDatabase.cxx
@@ -56,7 +56,7 @@ void clitk::AnatomicalFeatureDatabase::Load()
   // load line by line string/string
   while (!is.fail()) {
     std::string tag;
-    is >> tag; DD(tag);
+    is >> tag; 
     std::string value;
     std::getline(is,value,'\n');
     m_MapOfTag[tag] = value;
diff --git a/segmentation/clitkExtractLung.ggo b/segmentation/clitkExtractLung.ggo
index e91a2ea..37eada2 100644
--- a/segmentation/clitkExtractLung.ggo
+++ b/segmentation/clitkExtractLung.ggo
@@ -15,8 +15,7 @@ section "I/O"
 
 option "input"         i	"Input CT image filename"	  string  	yes
 option "patient"       p	"Input patient mask filename"	  string  	yes
-option "patientBG"   -	"Patient Background"              int           default="0" no
-option "afdb"          a	"Output Anatomical Feature DB (Carina position)"    string  	no
+option "patientBG"     -	"Patient Background"              int           default="0" no
 option "output"        o 	"Output lungs mask filename"      string 	yes
 option "outputTrachea" t 	"Output trachea mask filename"    string 	no
 
@@ -51,3 +50,6 @@ option "remove3"     -	"Labels not to keep in lungs mask (trachea)"    int		no
 option "firstKeep3"  -	"First label to keep"			  int		no  	default="1"
 option "lastKeep3"   -	"Last label to keep"		          int		no  	default="2"
 
+section "Step 5 : openclose"
+option "openclose"		-  "Final OpenClose"		flag off
+option "opencloseRadius"   	-  "Final OpenClose radius"	int no default="1"
diff --git a/segmentation/clitkExtractLungFilter.h b/segmentation/clitkExtractLungFilter.h
index cf94b07..6ec698d 100644
--- a/segmentation/clitkExtractLungFilter.h
+++ b/segmentation/clitkExtractLungFilter.h
@@ -24,12 +24,9 @@
 #include "clitkDecomposeAndReconstructImageFilter.h"
 #include "clitkExplosionControlledThresholdConnectedImageFilter.h"
 #include "clitkSegmentationUtils.h"
-#include "clitkFilterWithAnatomicalFeatureDatabaseManagement.h"
-#include "tree.hh"
 
 // itk
 #include "itkStatisticsImageFilter.h"
-#include "itkTreeContainer.h"
 
 namespace clitk {
   
@@ -57,39 +54,10 @@ namespace clitk {
   */
   //--------------------------------------------------------------------
   
-
-  //--------------------------------------------------------------------
-
-class Bifurcation
-{
-public:
-  typedef itk::Index<3> IndexType;
-  typedef itk::Point<double, 3> PointType;
-  typedef double PixelType;
-  Bifurcation(IndexType _index, PixelType _l, PixelType _l1, PixelType _l2) {
-    index = _index;
-    _l = l;
-    _l1 = l1;
-    _l2 = l2;
-  }
-  IndexType index;
-  PointType point;
-  PixelType l;
-  PixelType l1;
-  PixelType l2;
-  typedef itk::Index<3> NodeType;
-  typedef tree<NodeType> TreeType;
-  typedef TreeType::iterator TreeIterator;
-  TreeIterator treeIter;
-};
-  //--------------------------------------------------------------------
-
-
   //--------------------------------------------------------------------
   template <class TImageType, class TMaskImageType>
   class ITK_EXPORT ExtractLungFilter: 
     public virtual clitk::FilterBase, 
-    public clitk::FilterWithAnatomicalFeatureDatabaseManagement,
     public itk::ImageToImageFilter<TImageType, TMaskImageType> 
   {
     
@@ -133,11 +101,6 @@ public:
     typedef typename InternalImageType::IndexType                    InternalIndexType;
     typedef LabelizeParameters<InternalPixelType>                    LabelParamType;
     
-    typedef Bifurcation BifurcationType;
-    typedef MaskImageIndexType NodeType;
-    typedef tree<NodeType> TreeType;
-    typedef typename TreeType::iterator TreeIterator;
-
     /** Connect inputs */
     void SetInput(const ImageType * image);
     void SetInputPatientMask(MaskImageType * mask, MaskImagePixelType BG);
@@ -216,15 +179,15 @@ public:
     itkGetConstMacro(LabelizeParameters3, LabelParamType*);
     GGO_DefineOption_LabelParam(3, SetLabelizeParameters3, LabelParamType);
 
-    // Step 5 options LungSeparation
-    //     itkSetMacro(FinalOpenClose, bool);
-    //     itkGetConstMacro(FinalOpenClose, bool);
-    //     itkBooleanMacro(FinalOpenClose);
+    // Step 5 final openclose
+    itkSetMacro(FinalOpenClose, bool);
+    itkGetConstMacro(FinalOpenClose, bool);
+    itkBooleanMacro(FinalOpenClose);
+    GGO_DefineOption_Flag(openclose, SetFinalOpenClose);
 
-    // Bronchial bifurcations
-    itkSetMacro(FindBronchialBifurcations, bool);
-    itkGetConstMacro(FindBronchialBifurcations, bool);
-    itkBooleanMacro(FindBronchialBifurcations);
+    itkSetMacro(FinalOpenCloseRadius, int);
+    itkGetConstMacro(FinalOpenCloseRadius, int);
+    GGO_DefineOption(opencloseRadius, SetFinalOpenCloseRadius, int);
 
   protected:
     ExtractLungFilter();
@@ -270,21 +233,14 @@ public:
     LabelParamType* m_LabelizeParameters3;
 
     // Step 5
-    //     bool m_FinalOpenClose;    
-    bool m_FindBronchialBifurcations;
-    
+    bool m_FinalOpenClose;    
+    int m_FinalOpenCloseRadius;
+
+    // Main functions
     virtual void GenerateOutputInformation();
     virtual void GenerateData();
-
-    TreeType m_SkeletonTree;
-
-    void TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations, 
-                            MaskImagePointer skeleton, 
-                            MaskImageIndexType index,
-                            MaskImagePixelType label, 
-			    TreeIterator currentNode);
-        
-
+    
+    // Functions for trachea extraction
     bool SearchForTracheaSeed(int skip);
     void SearchForTrachea();
     void TracheaRegionGrowing();
diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx
index 92b2d4f..3de027c 100644
--- a/segmentation/clitkExtractLungFilter.txx
+++ b/segmentation/clitkExtractLungFilter.txx
@@ -24,7 +24,6 @@
 #include "clitkSetBackgroundImageFilter.h"
 #include "clitkSegmentationUtils.h"
 #include "clitkAutoCropFilter.h"
-#include "clitkExtractSliceFilter.h"
 
 // itk
 #include "itkBinaryThresholdImageFilter.h"
@@ -33,13 +32,14 @@
 #include "itkOtsuThresholdImageFilter.h"
 #include "itkBinaryThinningImageFilter3D.h"
 #include "itkImageIteratorWithIndex.h"
+#include "itkBinaryMorphologicalOpeningImageFilter.h"
+#include "itkBinaryMorphologicalClosingImageFilter.h"
 
 //--------------------------------------------------------------------
 template <class ImageType, class MaskImageType>
 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 ExtractLungFilter():
   clitk::FilterBase(),
-  clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
   itk::ImageToImageFilter<ImageType, MaskImageType>()
 {
   SetNumberOfSteps(10);
@@ -84,8 +84,9 @@ ExtractLungFilter():
   p3->UseLastKeepOff();
   SetLabelizeParameters3(p3);
   
-  // Step 5 : find bronchial bifurcations
-  FindBronchialBifurcationsOn();
+  // Step 5
+  FinalOpenCloseOff();
+  SetFinalOpenCloseRadius(1);
 }
 //--------------------------------------------------------------------
 
@@ -158,7 +159,8 @@ SetArgsInfo(ArgsInfoType mArgsInfo)
   SetRadiusForTrachea_GGO(mArgsInfo);
   SetLabelizeParameters3_GGO(mArgsInfo);
   
-  SetAFDBFilename_GGO(mArgsInfo);
+  SetFinalOpenCloseRadius_GGO(mArgsInfo);
+  SetFinalOpenClose_GGO(mArgsInfo);
 }
 //--------------------------------------------------------------------
 
@@ -170,7 +172,6 @@ clitk::ExtractLungFilter<ImageType, MaskImageType>::
 GenerateOutputInformation() 
 { 
   Superclass::GenerateOutputInformation();
-  //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion());
 
   // Get input pointers
   patient = dynamic_cast<const MaskImageType*>(itk::ProcessObject::GetInput(1));
@@ -335,6 +336,37 @@ GenerateOutputInformation()
   working_image = cropFilter2->GetOutput();
   StopCurrentStep<InternalImageType>(working_image);
 
+  //--------------------------------------------------------------------
+  //--------------------------------------------------------------------
+  // Final OpenClose
+  if (GetFinalOpenClose()) {
+    StartNewStep("Open/Close"); 
+
+    // Structuring element
+    typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
+    KernelType structuringElement;
+    structuringElement.SetRadius(GetFinalOpenCloseRadius());
+    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();
+  }
+
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
   StartNewStep("Separate Left/Right lungs");
@@ -395,125 +427,7 @@ 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);
-    }
-  }
 }
 //--------------------------------------------------------------------
 
@@ -524,72 +438,12 @@ void
 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 GenerateData() 
 {
-  // If everything goes well, set the output
+  // Set the output
   this->GraftOutput(output); // not SetNthOutput
 }
 //--------------------------------------------------------------------
 
 
-//--------------------------------------------------------------------
-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>
 bool 
diff --git a/segmentation/clitkExtractLungGenericFilter.txx b/segmentation/clitkExtractLungGenericFilter.txx
index e942dcf..74317ed 100644
--- a/segmentation/clitkExtractLungGenericFilter.txx
+++ b/segmentation/clitkExtractLungGenericFilter.txx
@@ -95,7 +95,6 @@ void clitk::ExtractLungGenericFilter<ArgsInfoType>::UpdateWithInputImageType()
   typename OutputImageType::Pointer output = filter->GetOutput();
   this->template SetNextOutput<OutputImageType>(output); 
   this->template SetNextOutput<typename FilterType::MaskImageType>(filter->GetTracheaImage()); 
-  filter->WriteAFDB();
 }
 //--------------------------------------------------------------------
 
diff --git a/segmentation/clitkExtractLymphStationsFilter.txx b/segmentation/clitkExtractLymphStationsFilter.txx
index 919d480..b81ae07 100644
--- a/segmentation/clitkExtractLymphStationsFilter.txx
+++ b/segmentation/clitkExtractLymphStationsFilter.txx
@@ -117,7 +117,7 @@ GenerateOutputInformation() {
   if (!m_CarinaZPositionInMMIsSet) {
     ImagePointType carina;
     LoadAFDB();
-    GetAFDB()->GetPoint3D("Carena", carina);
+    GetAFDB()->GetPoint3D("Carina", carina);
     DD(carina);
     m_CarinaZPositionInMM = carina[2];
   }
@@ -125,16 +125,16 @@ GenerateOutputInformation() {
 
   // ----------------------------------------------------------------
   // ----------------------------------------------------------------
-  // Superior limit = carena
-  // Inferior limit = origine middle lobe bronchus
-  StartNewStep("Inf/Sup mediastinum limits with carena/bronchus");
+  // Superior limit = carina
+  // Inferior limit = origin right middle lobe bronchus
+  StartNewStep("Inf/Sup mediastinum limits with carina/bronchus");
   ImageRegionType region = m_mediastinum->GetLargestPossibleRegion(); DD(region);
   ImageSizeType size = region.GetSize();
   ImageIndexType index = region.GetIndex();
   DD(m_CarinaZPositionInMM);
   DD(m_MiddleLobeBronchusZPositionInMM);
   index[2] = floor((m_MiddleLobeBronchusZPositionInMM - m_mediastinum->GetOrigin()[2]) / m_mediastinum->GetSpacing()[2]);
-  size[2] = ceil((m_CarinaZPositionInMM-m_MiddleLobeBronchusZPositionInMM) / m_mediastinum->GetSpacing()[2]);
+  size[2] = 1+ceil((m_CarinaZPositionInMM-m_MiddleLobeBronchusZPositionInMM) / m_mediastinum->GetSpacing()[2]); // +1 to 
   region.SetSize(size);
   region.SetIndex(index);
   DD(region);
@@ -152,7 +152,7 @@ GenerateOutputInformation() {
   // ----------------------------------------------------------------
   // ----------------------------------------------------------------
   // Separate trachea in two CCL
-  StartNewStep("Separate trachea under carena");
+  StartNewStep("Separate trachea under carina");
   // DD(region);
   ImageRegionType trachea_region = m_trachea->GetLargestPossibleRegion();
   for(int i=0; i<3; i++) {
diff --git a/segmentation/clitkExtractMediastinum.ggo b/segmentation/clitkExtractMediastinum.ggo
index 92bff77..4a42029 100644
--- a/segmentation/clitkExtractMediastinum.ggo
+++ b/segmentation/clitkExtractMediastinum.ggo
@@ -23,6 +23,7 @@ option "lungLeft"      -	"Lung Left value"                 int           default
 option "lungRight"     -	"Lung Right value"                int           default="2" no
 option "trachea"       t	"Input trachea mask filename"	  string  	yes
 option "tracheaBG"     -	"Trachea Background"              int           default="0" no
+option "afdb"          a	"Input Anatomical Feature DB"     string  	no
 
 option "output"        o 	"Output lungs mask filename"      string 	yes
 
diff --git a/segmentation/clitkExtractMediastinumFilter.h b/segmentation/clitkExtractMediastinumFilter.h
index 0c4e444..4b64d13 100644
--- a/segmentation/clitkExtractMediastinumFilter.h
+++ b/segmentation/clitkExtractMediastinumFilter.h
@@ -20,6 +20,7 @@
 #define CLITKEXTRACTMEDIASTINUMFILTER_H
 
 #include "clitkFilterBase.h"
+#include "clitkFilterWithAnatomicalFeatureDatabaseManagement.h"
 
 namespace clitk {
   
@@ -36,7 +37,8 @@ namespace clitk {
   
   template <class TImageType>
   class ITK_EXPORT ExtractMediastinumFilter: 
-    public clitk::FilterBase, 
+    public virtual clitk::FilterBase, 
+    public clitk::FilterWithAnatomicalFeatureDatabaseManagement,
     public itk::ImageToImageFilter<TImageType, TImageType> 
   {
 
diff --git a/segmentation/clitkExtractMediastinumFilter.txx b/segmentation/clitkExtractMediastinumFilter.txx
index 1f0d5f2..7f5ca18 100644
--- a/segmentation/clitkExtractMediastinumFilter.txx
+++ b/segmentation/clitkExtractMediastinumFilter.txx
@@ -24,10 +24,12 @@
 #include "clitkExtractMediastinumFilter.h"
 #include "clitkAddRelativePositionConstraintToLabelImageFilter.h"
 #include "clitkSegmentationUtils.h"
-#include "clitkExtractAirwayTreeInfoFilter.h"
+#include "clitkExtractAirwaysTreeInfoFilter.h"
 
-// itk
+// std
 #include <deque>
+
+// itk
 #include "itkStatisticsLabelMapFilter.h"
 #include "itkLabelImageToStatisticsLabelMapFilter.h"
 #include "itkRegionOfInterestImageFilter.h"
@@ -41,6 +43,7 @@ template <class ImageType>
 clitk::ExtractMediastinumFilter<ImageType>::
 ExtractMediastinumFilter():
   clitk::FilterBase(),
+  clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
   itk::ImageToImageFilter<ImageType, ImageType>()
 {
   this->SetNumberOfRequiredInputs(4);
@@ -256,17 +259,12 @@ GenerateData()
 
   // Find ant-post coordinate of trachea (assume the carena position is a
   // good estimation of the ant-post position of the trachea)
-  typedef clitk::ExtractAirwayTreeInfoFilter<ImageType> AirwayFilter;
-  typename AirwayFilter::Pointer airwayfilter = AirwayFilter::New();
-  airwayfilter->SetVerboseStep(false);
-  airwayfilter->SetWriteStep(false);
-  airwayfilter->SetInput(trachea);
-  airwayfilter->Update();
-  DD(airwayfilter->GetFirstTracheaPoint());
-  ImagePointType point_trachea = airwayfilter->GetCarinaPoint();
+  ImagePointType carina;
+  LoadAFDB();
+  GetAFDB()->GetPoint3D("Carina", carina);
+  DD(carina);
   ImageIndexType index_trachea;
-  bones->TransformPhysicalPointToIndex(point_trachea, index_trachea);
-  DD(point_trachea);
+  bones->TransformPhysicalPointToIndex(carina, index_trachea);
   DD(index_trachea);
   
   // Split bone image first into two parts (ant and post)
diff --git a/segmentation/clitkFilterWithAnatomicalFeatureDatabaseManagement.h b/segmentation/clitkFilterWithAnatomicalFeatureDatabaseManagement.h
index fd9ecc1..9bcd9d9 100644
--- a/segmentation/clitkFilterWithAnatomicalFeatureDatabaseManagement.h
+++ b/segmentation/clitkFilterWithAnatomicalFeatureDatabaseManagement.h
@@ -42,19 +42,26 @@ namespace clitk {
     itkTypeMacro(FilterWithAnatomicalFeatureDatabaseManagement, Object);
 
     // Set/Get filename 
+    itkBooleanMacro(AFDBFilenameGivenFlag);
+    itkSetMacro(AFDBFilenameGivenFlag, bool);
+    itkGetConstMacro(AFDBFilenameGivenFlag, bool);
+    GGO_DefineOption_Flag(afdb, SetAFDBFilenameGivenFlag);
+
     itkSetMacro(AFDBFilename, std::string);
     itkGetConstMacro(AFDBFilename, std::string);
-    GGO_DefineOption(afdb, SetAFDBFilename, std::string);
+    GGO_DefineOption_WithTest(afdb, SetAFDBFilename, std::string, AFDBFilenameGivenFlag);
 
     void WriteAFDB();
     void LoadAFDB();
     AnatomicalFeatureDatabase * GetAFDB();
+    void SetAFDB(AnatomicalFeatureDatabase * a) { m_AFDB = a; }
 
   protected:
     FilterWithAnatomicalFeatureDatabaseManagement();
     virtual ~FilterWithAnatomicalFeatureDatabaseManagement() {}    
     
     std::string m_AFDBFilename;
+    bool m_AFDBFilenameGivenFlag;
     clitk::AnatomicalFeatureDatabase * m_AFDB;
 
   private: