]> Creatis software - clitk.git/blobdiff - segmentation/clitkExtractLymphStation_8.txx
Merge branch 'master' of /home/dsarrut/clitk3.server
[clitk.git] / segmentation / clitkExtractLymphStation_8.txx
index 469ff2f836f38a053739151f642c79dc461c8163..d4bbfc389217188c9ba4470d63a56dec39b33a65 100644 (file)
@@ -2,6 +2,48 @@
 #include <itkBinaryDilateImageFilter.h>
 #include <itkMirrorPadImageFilter.h>
 
+//--------------------------------------------------------------------
+template <class ImageType>
+void 
+clitk::ExtractLymphStationsFilter<ImageType>::
+ExtractStation_8_SetDefaultValues()
+{
+  MaskImagePointType p;
+  p[0] = 15; p[1] = 2; p[2] = 1;
+  SetEsophagusDiltationForAnt(p);
+  p[0] = 5; p[1] = 10; p[2] = 1;
+  SetEsophagusDiltationForRight(p);
+  SetFuzzyThreshold("8", "Esophagus", 0.5);
+  SetInjectedThresholdForS8(150);
+}
+//--------------------------------------------------------------------
+
+//--------------------------------------------------------------------
+template <class TImageType>
+void 
+clitk::ExtractLymphStationsFilter<TImageType>::
+ExtractStation_8() 
+{
+  if (!CheckForStation("8")) return;
+
+  StartNewStep("Station 8");
+  StartSubStep();   
+  ExtractStation_8_SI_Limits();         // OK, validated
+  ExtractStation_8_Ant_Limits();        // OK, validated
+  ExtractStation_8_Left_Sup_Limits();   // OK, validated
+  ExtractStation_8_Left_Inf_Limits();   // OK, validated
+  ExtractStation_8_Single_CCL_Limits(); // OK, validated
+  ExtractStation_8_Remove_Structures(); // OK, validated
+  
+  // Store image filenames into AFDB 
+  writeImage<MaskImageType>(m_ListOfStations["8"], "seg/Station8.mhd");
+  GetAFDB()->SetImageFilename("Station8", "seg/Station8.mhd");  
+  WriteAFDB();
+  StopSubStep();
+}
+//--------------------------------------------------------------------
+
+
 //--------------------------------------------------------------------
 template <class ImageType>
 void 
@@ -25,72 +67,42 @@ ExtractStation_8_SI_Limits()
     left; the lower border of the bronchus intermedius on the right"
 
   */
-  StartNewStep("[Station8] Inf/Sup mediastinum limits with carina/diaphragm junction");
-
-  // Get Carina Z position
-  MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after read Carina"); 
-  std::vector<MaskImagePointType> centroids;
-  clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
-  DD(centroids.size());
-  DD(centroids[0]); // BG
-  DD(centroids[1]);
-  m_CarinaZ = centroids[1][2];
-  DD(m_CarinaZ);
-  // add one slice to include carina ?
-  m_CarinaZ += m_Mediastinum->GetSpacing()[2];
-  DD(m_CarinaZ);
-  DD(Carina->GetReferenceCount());
-  // We dont need Carina structure from now
-  Carina->Delete();
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after delete Carina"); 
-  
-  // Get left lower lobe bronchus (L), take the upper border
-  // Get right bronchus intermedius (RML), take the lower border
+  StartNewStep("[Station8] Sup/Inf limits with LeftLower/RightMiddle Lobe and diaphragm");
 
-  /*
-    double m_CarinaZ = GetAFDB()->GetPoint3D("Carina", 2) + 
-    m_Mediastinum->GetSpacing()[2]; // add one slice to include carina
-    // double GOjunctionZ = GetAFDB()->GetPoint3D("GOjunction", 2);
-    */
+  /* -----------------------------------------------
+     NEW SUPERIOR LIMIT = LeftLowerLobeBronchus /
+     RightMiddleLobeBronchus See FindLineForS7S8Separation
+     -----------------------------------------------
+  */
+  ImagePointType A;
+  ImagePointType B;
+  FindLineForS7S8Separation(A, B);
+
+  // add one slice to be adjacent to Station7
+  B[2] += m_Working_Support->GetSpacing()[2];
+  A[2] += m_Working_Support->GetSpacing()[2];
+
+  // Use the line to remove the inferior part
+  m_Working_Support =
+    SliceBySliceSetBackgroundFromSingleLine<MaskImageType>(m_Working_Support, 
+                                                           GetBackgroundValue(), A, B, 2, 0, true);
+  
+  /* -----------------------------------------------
+     INFERIOR LIMIT = Diaphragm
+     -----------------------------------------------
+  */  
 
   // Found most inferior part of the lung
   MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
   // It should be already croped, so I took the origin and add 10mm above 
   m_DiaphragmInferiorLimit = Lungs->GetOrigin()[2]+10;
-  DD(m_DiaphragmInferiorLimit);
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after read Lung");  
-  Lungs->Delete(); // we don't need it, release memory
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after release Lung");  
-
-  /* Crop support :
-     Superior limit = carina
-     Inferior limit = DiaphragmInferiorLimit (old=Gastroesphogeal junction) */
+  GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
+
   m_Working_Support = 
     clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2, 
-                                                m_DiaphragmInferiorLimit, //GOjunctionZ, 
-                                                m_CarinaZ, true,
+                                                m_DiaphragmInferiorLimit,
+                                                B[2], true,
                                                 GetBackgroundValue());
-  
-  // Remove some structures that we know are excluded 
-  MaskImagePointer VertebralBody = 
-    GetAFDB()->template GetImage <MaskImageType>("VertebralBody");  
-  MaskImagePointer Aorta = 
-    GetAFDB()->template GetImage <MaskImageType>("Aorta");  
-
-  typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
-  typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); 
-  boolFilter->InPlaceOn();
-  boolFilter->SetInput1(m_Working_Support);
-  boolFilter->SetInput2(VertebralBody);    
-  boolFilter->SetOperationType(BoolFilterType::AndNot);
-  boolFilter->Update();    
-  boolFilter->SetInput1(boolFilter->GetOutput());
-  boolFilter->SetInput2(Aorta);    
-  boolFilter->SetOperationType(BoolFilterType::AndNot);
-  boolFilter->Update();    
-  m_Working_Support = boolFilter->GetOutput();
-
   // Done.
   StopCurrentStep<MaskImageType>(m_Working_Support);
   m_ListOfStations["8"] = m_Working_Support;
@@ -102,597 +114,346 @@ ExtractStation_8_SI_Limits()
 template <class ImageType>
 void 
 clitk::ExtractLymphStationsFilter<ImageType>::
-ExtractStation_8_AP_Limits() 
+ExtractStation_8_Ant_Limits() 
 {
+  //--------------------------------------------------------------------
+  StartNewStep("[Station8] Ant part: not post to Esophagus");
   /*
-    Station 8: paraeosphageal nodes
+    Consider Esophagus, dilate it and remove ant part. It remains part
+    on L & R, than can be partly removed by cutting what remains at
+    right of vertebral body.
+  */
+  
+  // Get Esophagus
+  m_Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
+
+  // In images from the original article, Atlas – UM, the oesophagus
+  //was included in nodal stations 3p and 8.  Having said that, in the
+  //description for station 8, it indicates that “the delineation of
+  //station 8 is limited to the soft tissues surrounding the
+  //oesophagus”.  In the recent article, The IASLC Lung Cancer Staging
+  //Project (J Thorac Oncol 4:5, 568-77), the images drawn by
+  //Dr. Aletta Frasier exclude this structure.  From an oncological
+  //prospective, the oesophagus should be excluded from these nodal
+  //stations.
+
+  // Resize Esophagus like current support
+  m_Esophagus = 
+    clitk::ResizeImageLike<MaskImageType>(m_Esophagus, m_Working_Support, GetBackgroundValue()); // Needed ?
 
-    Anteriorly, it is in contact with Station 7 and the
-    left main stem bronchus in its superior aspect (Fig. 3C–H) and
-    with the heart more inferiorly. Posteriorly, Station 8 abuts the
-    descending aorta and the anterior aspect of the vertebral body
-    until an imaginary horizontal line running 1 cm posterior to the
-    anterior border of the vertebral body (Fig. 3C). 
+  // Dilate to keep only not Anterior positions
+  MaskImagePointType radiusInMM = GetEsophagusDiltationForAnt();
+  m_Esophagus = clitk::Dilate<MaskImageType>(m_Esophagus, 
+                                             radiusInMM, 
+                                             GetBackgroundValue(), 
+                                             GetForegroundValue(), true);
+  // Keep what is Posterior to Esophagus
+  typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
+  typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
+  relPosFilter->VerboseStepFlagOff();
+  relPosFilter->WriteStepFlagOff();
+  relPosFilter->SetInput(m_Working_Support);
+  relPosFilter->SetInputObject(m_Esophagus);
+  relPosFilter->AddOrientationTypeString("PostTo");
+  //  relPosFilter->InverseOrientationFlagOff();
+  relPosFilter->SetDirection(2); // Z axis
+  relPosFilter->UniqueConnectedComponentBySliceFlagOff();
+  relPosFilter->SetIntermediateSpacing(3);
+  relPosFilter->IntermediateSpacingFlagOn();
+  relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("8", "Esophagus"));
+  relPosFilter->RemoveObjectFlagOff(); // Do not exclude here because it is dilated
+  relPosFilter->CombineWithOrFlagOff(); // NO !
+  relPosFilter->IgnoreEmptySliceObjectFlagOn();
+  relPosFilter->Update();
+  m_Working_Support = relPosFilter->GetOutput();
 
-    New classification IASCL 2009 :
+  // AutoCrop (OR SbS ?)
+  m_Working_Support = clitk::AutoCrop<MaskImageType>(m_Working_Support, GetBackgroundValue()); 
 
-    "Nodes lying adjacent to the wall of the esophagus and to the
-    right or left of the midline, excluding subcarinal nodes (S7)
-    Upper"
+  StopCurrentStep<MaskImageType>(m_Working_Support);
+  m_ListOfStations["8"] = m_Working_Support;
+}
+//--------------------------------------------------------------------
 
-  */
 
-  StartNewStep("[Station8] Post limits with vertebral body");
-  MaskImagePointer VertebralBody = 
-    GetAFDB()->template GetImage <MaskImageType>("VertebralBody");  
+//--------------------------------------------------------------------
+template <class ImageType>
+void 
+clitk::ExtractLymphStationsFilter<ImageType>::
+ExtractStation_8_Left_Sup_Limits() 
+{
+  //--------------------------------------------------------------------
+  StartNewStep("[Station8] Left limits: remove Left to line from Aorta to PulmonaryTrunk");
 
-  // Consider vertebral body slice by slice
-  typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
-  typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
-  typedef typename ExtractSliceFilterType::SliceType SliceType;
-  std::vector<typename SliceType::Pointer> vertebralSlices;
-  extractSliceFilter->SetInput(VertebralBody);
-  extractSliceFilter->SetDirection(2);
-  extractSliceFilter->Update();
-  extractSliceFilter->GetOutputSlices(vertebralSlices);
-
-  // For each slice, compute the most anterior point
-  std::map<int, typename SliceType::PointType> vertebralAntPositionBySlice;
-  for(uint i=0; i<vertebralSlices.size(); i++) {
-    typename SliceType::PointType p;
-    bool found = clitk::FindExtremaPointInAGivenDirection<SliceType>(vertebralSlices[i], 
-                                                                     GetBackgroundValue(), 
-                                                                     1, true, p);
+  /*
+    We consider a line from Left part of Aorta to left part of
+    PulmonaryTrunk and remove what is at Left.
+  */
+    
+  // Load the structures
+  MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
+  MaskImagePointer PulmonaryTrunk = GetAFDB()->template GetImage<MaskImageType>("PulmonaryTrunk");
+  
+  // Resize like the PT and define the slices
+  MaskImagePointType min, max;
+  clitk::GetMinMaxPointPosition<MaskImageType>(PulmonaryTrunk, min, max);
+  Aorta = clitk::CropImageAlongOneAxis<MaskImageType>(Aorta, 2, min[2], max[2], false, GetBackgroundValue());
+  std::vector<MaskSlicePointer> slices_aorta;
+  clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_aorta);
+  std::vector<MaskSlicePointer> slices_PT;
+  clitk::ExtractSlices<MaskImageType>(PulmonaryTrunk, 2, slices_PT);
+  
+  // Find the points at left
+  std::vector<MaskImagePointType> p_A;
+  std::vector<MaskImagePointType> p_B;
+  MaskImagePointType pA;
+  MaskImagePointType pB;
+  for(uint i=0; i<slices_PT.size(); i++) {    
+    typename MaskSliceType::PointType ps;
+    // In Aorta (assume a single CCL)
+    bool found = 
+      clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices_aorta[i], GetBackgroundValue(), 0, false, ps);
+    clitk::PointsUtils<MaskImageType>::Convert2DTo3D(ps, Aorta, i, pA);
+    
     if (found) {
-      vertebralAntPositionBySlice[i] = p;
+      // In PT : generally 2 CCL, we keep the most at left
+      found = 
+        clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices_PT[i], GetBackgroundValue(), 0, false, ps);
+      clitk::PointsUtils<MaskImageType>::Convert2DTo3D(ps, PulmonaryTrunk, i, pB);
     }
-    else { 
-      // no Foreground in this slice (it appends generally at the
-      // beginning and the end of the volume). Do nothing in this
-      // case.
+    
+    if (found) {
+      p_A.push_back(pA);
+      p_B.push_back(pB);
     }
   }
+  clitk::WriteListOfLandmarks<MaskImageType>(p_A, "S8-Aorta-Left-points.txt");
+  clitk::WriteListOfLandmarks<MaskImageType>(p_B, "S8-PT-Left-points.txt");
 
-  // Convert 2D points in slice into 3D points
-  std::vector<MaskImagePointType> vertebralAntPositions;
-  clitk::PointsUtils<MaskImageType>::Convert2DTo3DList(vertebralAntPositionBySlice, 
-                                                       VertebralBody, 
-                                                       vertebralAntPositions);
-
-  // DEBUG : write list of points
-  clitk::WriteListOfLandmarks<MaskImageType>(vertebralAntPositions, 
-                                             "vertebralMostAntPositions.txt");
-
-  // Cut support posteriorly 1cm the most anterior point of the
-  // VertebralBody. Do nothing below and above the VertebralBody. To
-  // do that compute several region, slice by slice and fill. 
-  MaskImageRegionType region;
-  MaskImageSizeType size;
-  MaskImageIndexType start;
-  size[2] = 1; // a single slice
-  start[0] = m_Working_Support->GetLargestPossibleRegion().GetIndex()[0];
-  size[0] = m_Working_Support->GetLargestPossibleRegion().GetSize()[0];
-  for(uint i=0; i<vertebralAntPositions.size(); i++) {
-    typedef typename itk::ImageRegionIterator<MaskImageType> IteratorType;
-    IteratorType iter = 
-      IteratorType(m_Working_Support, m_Working_Support->GetLargestPossibleRegion());
-    MaskImageIndexType index;
-    // Consider some cm posterior to most anterior positions (usually
-    // 1 cm).
-    vertebralAntPositions[i][1] += GetDistanceMaxToAnteriorPartOfTheSpine();
-    // Get index of this point
-    m_Working_Support->TransformPhysicalPointToIndex(vertebralAntPositions[i], index);
-    // Compute region (a single slice)
-    start[2] = index[2];
-    start[1] = m_Working_Support->GetLargestPossibleRegion().GetIndex()[1]+index[1];
-    size[1] = m_Working_Support->GetLargestPossibleRegion().GetSize()[1]-start[1];
-    region.SetSize(size);
-    region.SetIndex(start);
-    // Fill region
-    if (m_Working_Support->GetLargestPossibleRegion().IsInside(start))  {
-      itk::ImageRegionIterator<MaskImageType> it(m_Working_Support, region);
-      it.GoToBegin();
-      while (!it.IsAtEnd()) {
-        it.Set(GetBackgroundValue());
-        ++it;
-      }
-    }
-  }
+  // Remove part at Left  
+  clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support, 
+                                                                    p_A, p_B,
+                                                                    GetBackgroundValue(), 
+                                                                    0, -10);
 
   StopCurrentStep<MaskImageType>(m_Working_Support);
   m_ListOfStations["8"] = m_Working_Support;
+}
+//--------------------------------------------------------------------
 
-  //--------------------------------------------------------------------
-  StartNewStep("[Station8] Ant limits with S7 above Carina");
-  /*
-    Anteriorly, it is in contact with Station 7 and the
-    left main stem bronchus in its superior aspect (Fig. 3C–H) and
-    with the heart more inferiorly. 
-
-    1) line post wall bronchi (S7), till originRMLB
-    - LUL bronchus : to detect in trachea. But not needed here ??
-    2) heart ! -> not delineated.
-    ==> below S7, crop CT not to far away (ant), then try with threshold
-    
-    1) ==> how to share with S7 ? Prepare both support at the same time !
-    NEED vector of initial support for each station ? No -> map if it exist before, take it !!
-
-  */
-
-  // Ant limit from carina (start) to end of S7 = originRMLB
-  // Get Trachea
-  MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after read Trachea");
-  // Crop above Carina
-  //double m_CarinaZ = GetAFDB()->GetPoint3D("Carina", 2) + 
-  //    Trachea->GetSpacing()[2]; // add one slice to include carina;
-
-  MaskImagePointer m_Working_Trachea = 
-    clitk::CropImageAbove<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
-                                         GetBackgroundValue());
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after crop");
-
-  // Seprate into two main bronchi
-  MaskImagePointer LeftBronchus;
-  MaskImagePointer RightBronchus;
-
-  // Labelize and consider the two first (main) labels
-  m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after labelize");
-
-  // Carina position must at the first slice that separate the two
-  // main bronchus (not superiorly). We thus first check that the
-  // upper slice is composed of at least two labels
-  typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
-  SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
-  iter.SetFirstDirection(0);
-  iter.SetSecondDirection(1);
-  iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
-  int maxLabel=0;
-  while (!iter.IsAtReverseEndOfSlice()) {
-    while (!iter.IsAtReverseEndOfLine()) {    
-      if (iter.Get() > maxLabel) maxLabel = iter.Get();
-      --iter;
-    }
-    iter.PreviousLine();
-  }
-  if (maxLabel < 2) {
-    clitkExceptionMacro("First slice form Carina does not seems to seperate the two main bronchus. Abort");
-  }
 
-  // Compute 3D centroids of both parts to identify the left from the
-  // right bronchus
-  std::vector<ImagePointType> c;
-  clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
-  ImagePointType C1 = c[1];
-  ImagePointType C2 = c[2];
-
-  ImagePixelType leftLabel;
-  ImagePixelType rightLabel;  
-  if (C1[0] < C2[0]) { leftLabel = 1; rightLabel = 2; }
-  else { leftLabel = 2; rightLabel = 1; }
-
-  // Select LeftLabel (set one label to Backgroundvalue)
-  LeftBronchus = 
-    SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea, 
-                                                rightLabel, GetBackgroundValue(), false);
-  RightBronchus  = 
-    SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea, 
-                                                leftLabel, GetBackgroundValue(), false);
-  // Crop images
-  LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue()); 
-  RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue()); 
-
-  // Insert int AFDB if need after 
-  GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd", 
-                                               LeftBronchus, true);
-  GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd", 
-                                               RightBronchus, true);
-
-
-  // Now crop below OriginOfRightMiddleLobeBronchusZ
-  // It is not done before to keep entire bronchi.
-  
-  /*
-    double OriginOfRightMiddleLobeBronchusZ = 
-    GetAFDB()->GetPoint3D("OriginOfRightMiddleLobeBronchus", 2)+
-    LeftBronchus->GetSpacing()[2]; 
-    // ^--> Add one slice because the origin is the first slice without S7
-    DD(OriginOfRightMiddleLobeBronchusZ);
-    DD(OriginOfRightMiddleLobeBronchusZ-LeftBronchus->GetSpacing()[2]);  
-  */
+//--------------------------------------------------------------------
+template <class ImageType>
+void 
+clitk::ExtractLymphStationsFilter<ImageType>::
+ExtractStation_8_Single_CCL_Limits() 
+{
+  //--------------------------------------------------------------------
+  StartNewStep("[Station8 or 3P] Slice by slice, keep a single CCL (the closest to VertebralBody)");
 
-  MaskImagePointer OriginOfRightMiddleLobeBronchus = 
-    GetAFDB()->template GetImage<MaskImageType>("OriginOfRightMiddleLobeBronchus");
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after read OriginOfRightMiddleLobeBronchus"); 
-  std::vector<MaskImagePointType> centroids;
-  clitk::ComputeCentroids<MaskImageType>(OriginOfRightMiddleLobeBronchus, GetBackgroundValue(), centroids);
-  DD(centroids.size());
-  DD(centroids[0]); // BG
-  DD(centroids[1]);
-  m_OriginOfRightMiddleLobeBronchusZ = centroids[1][2];
-  DD(m_OriginOfRightMiddleLobeBronchusZ);
-  // add one slice to include carina ?
-  m_OriginOfRightMiddleLobeBronchusZ += LeftBronchus->GetSpacing()[2];
-  DD(m_OriginOfRightMiddleLobeBronchusZ);
-  DD(OriginOfRightMiddleLobeBronchus->GetReferenceCount());
-  // We dont need Carina structure from now
-  OriginOfRightMiddleLobeBronchus->Delete();
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after delete OriginOfRightMiddleLobeBronchus"); 
-
-
-  LeftBronchus = 
-    clitk::CropImageBelow<MaskImageType>(LeftBronchus, 2, 
-                                         m_OriginOfRightMiddleLobeBronchusZ, 
-                                         true, // AutoCrop
-                                         GetBackgroundValue());
-  RightBronchus = 
-    clitk::CropImageBelow<MaskImageType>(RightBronchus, 2, 
-                                         m_OriginOfRightMiddleLobeBronchusZ, 
-                                         true, // AutoCrop
-                                         GetBackgroundValue());
-
-  // Search for points that are the most left/post/ant and most
-  // right/post/ant of the left and right bronchus
-  // 15  = not 15 mm more distance than the middle point.
-  FindExtremaPointsInBronchus(LeftBronchus, 0, 10, m_RightMostInLeftBronchus, 
-                             m_AntMostInLeftBronchus, m_PostMostInLeftBronchus);
-  FindExtremaPointsInBronchus(RightBronchus, 1, 10, m_LeftMostInRightBronchus, 
-                             m_AntMostInRightBronchus, m_PostMostInRightBronchus);
-
-  // DEBUG : write the list of points
-  ListOfPointsType v;
-  v.reserve(m_RightMostInLeftBronchus.size()+m_AntMostInLeftBronchus.size()+
-            m_PostMostInLeftBronchus.size());
-  v.insert(v.end(), m_RightMostInLeftBronchus.begin(), m_RightMostInLeftBronchus.end() );
-  v.insert(v.end(), m_AntMostInLeftBronchus.begin(), m_AntMostInLeftBronchus.end() );
-  v.insert(v.end(), m_PostMostInLeftBronchus.begin(), m_PostMostInLeftBronchus.end() );
-  clitk::WriteListOfLandmarks<MaskImageType>(v, "LeftBronchusPoints.txt");
-
-  v.clear();
-  v.reserve(m_LeftMostInRightBronchus.size()+m_AntMostInRightBronchus.size()+
-            m_PostMostInRightBronchus.size());
-  v.insert(v.end(), m_LeftMostInRightBronchus.begin(), m_LeftMostInRightBronchus.end() );
-  v.insert(v.end(), m_AntMostInRightBronchus.begin(), m_AntMostInRightBronchus.end() );
-  v.insert(v.end(), m_PostMostInRightBronchus.begin(), m_PostMostInRightBronchus.end() );
-  clitk::WriteListOfLandmarks<MaskImageType>(v, "RightBronchusPoints.txt");
-
-
-  // Now uses these points to limit, slice by slice 
-  // http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
-  /*
-    Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
-    (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
-    This will equal zero if the point C is on the line formed by
-    points A and B, and will have a different sign depending on the
-    side. Which side this is depends on the orientation of your (x,y)
-    coordinates, but you can plug test values for A,B and C into this
-    formula to determine whether negative values are to the left or to
-    the right.
-    => to accelerate, start with formula, when change sign -> stop and fill
-  */
-  typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
-  SliceIteratorType siter = SliceIteratorType(m_Working_Support, 
-                                              m_Working_Support->GetLargestPossibleRegion());
-  siter.SetFirstDirection(0);
-  siter.SetSecondDirection(1);
-  siter.GoToBegin();
-  int i=0;
-  MaskImageType::PointType A;
-  MaskImageType::PointType B;
-  MaskImageType::PointType C;
-  while (!siter.IsAtEnd()) {
-    // Check that the current slice correspond to the current point
-    m_Working_Support->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
-    if (C[2] != m_PostMostInLeftBronchus[i][2]) {
-      // m_Working_Support start from GOjunction while list of point
-      // start at OriginOfRightMiddleLobeBronchusZ, so we must skip some slices.
-    }
-    else {
-      // Define A,B,C points
-      A = m_PostMostInLeftBronchus[i];
-      B = m_PostMostInRightBronchus[i];
-      C = A;
-      C[1] += 10; // I know I must keep this point
-      double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
-      bool isPositive = s<0;
-      while (!siter.IsAtEndOfSlice()) {
-        while (!siter.IsAtEndOfLine()) {
-          // Very slow, I know ... but image should be very small
-          m_Working_Support->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
-          double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
-          if (s == 0) siter.Set(2);
-          if (isPositive) {
-            if (s > 0) siter.Set(GetBackgroundValue());
-          }
-          else {
-            if (s < 0) siter.Set(GetBackgroundValue());
-          }
-          ++siter;
+  // Consider slices
+  std::vector<typename MaskSliceType::Pointer> slices;
+  std::vector<typename MaskSliceType::Pointer> slices_vb;
+  clitk::ExtractSlices<MaskImageType>(m_Working_Support, 2, slices);
+  MaskImagePointer VertebralBody = 
+    GetAFDB()->template GetImage <MaskImageType>("VertebralBody");  
+  VertebralBody = clitk::ResizeImageLike<MaskImageType>(VertebralBody, m_Working_Support, GetBackgroundValue());
+  clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, slices_vb);
+
+  for(uint i=0; i<slices.size(); i++) {
+    // Decompose in labels
+    slices[i] = Labelize<MaskSliceType>(slices[i], 0, true, 100);
+    // Compute centroids coordinate
+    std::vector<typename MaskSliceType::PointType> centroids;
+    std::vector<typename MaskSliceType::PointType> c;
+    clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), centroids);
+    clitk::ComputeCentroids<MaskSliceType>(slices_vb[i], GetBackgroundValue(), c);
+    if ((c.size() > 1) && (centroids.size() > 1)) {
+      // keep the one which is the closer to vertebralbody centroid
+      double distance=1000000;
+      int index=0;
+      for(uint j=1; j<centroids.size(); j++) {
+        double d = centroids[j].EuclideanDistanceTo(c[1]);
+        if (d<distance) {
+          distance = d;
+          index = j;
         }
-        siter.NextLine();
       }
-      ++i;
+      // remove all others label
+      slices[i] = KeepLabels<MaskSliceType>(slices[i], GetBackgroundValue(), 
+                                            GetForegroundValue(), index, index, true);
     }
-    siter.NextSlice();
   }
-
-  // Keep main 3D CCL :
-  m_Working_Support = Labelize<MaskImageType>(m_Working_Support, 0, false, 10);
-  m_Working_Support = KeepLabels<MaskImageType>(m_Working_Support, 
-                                                GetBackgroundValue(), 
-                                                GetForegroundValue(), 1, 1, true);
+  m_Working_Support = clitk::JoinSlices<MaskImageType>(slices, m_Working_Support, 2);
   
-  // Autocrop
-  m_Working_Support = clitk::AutoCrop<MaskImageType>(m_Working_Support, GetBackgroundValue()); 
-
-  // End of step
   StopCurrentStep<MaskImageType>(m_Working_Support);
-  //  m_ListOfStations["8"] = m_Working_Support;
-
-  //--------------------------------------------------------------------
-  StartNewStep("[Station8] Ant limits arround esophagus below Carina");
-  /*
-    Consider Esophagus, dilate it and remove ant part. It remains part
-    on L & R, than can be partly removed by cutting what remains at
-    right of vertebral body.
-  */
-  
-  // Get Esophagus
-  DD("Esophagus");
-  MaskImagePointer Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after read Esophagus");
-
-  // Crop Esophagus : keep only below the OriginOfRightMiddleLobeBronchusZ
-  Esophagus = 
-    clitk::CropImageAbove<MaskImageType>(Esophagus, 2, 
-                                         m_OriginOfRightMiddleLobeBronchusZ, 
-                                         true, // AutoCrop
-                                         GetBackgroundValue());
-
-  // Dilate to keep only not Anterior positions
-  MaskImagePointType radiusInMM = GetEsophagusDiltationForAnt();
-  Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(Esophagus);
-  Esophagus = clitk::Dilate<MaskImageType>(Esophagus, 
-                                           radiusInMM, 
-                                           GetBackgroundValue(), 
-                                           GetForegroundValue(), true);
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after dilate Esophagus");
-  writeImage<MaskImageType>(Esophagus, "enlarged-eso.mhd");
-
-  // Remove Anterior part according to this dilatated esophagus
-  typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
-  typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
-  relPosFilter->VerboseStepFlagOff();
-  relPosFilter->SetInput(m_Working_Support);
-  relPosFilter->SetInputObject(Esophagus);
-  relPosFilter->AddOrientationTypeString("P");
-  relPosFilter->InverseOrientationFlagOff();
-  relPosFilter->SetDirection(2); // Z axis
-  relPosFilter->UniqueConnectedComponentBySliceOff();
-  relPosFilter->SetIntermediateSpacing(3);
-  relPosFilter->ResampleBeforeRelativePositionFilterOn();
-  relPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS8());
-  relPosFilter->RemoveObjectFlagOff();
-  relPosFilter->Update();
-  m_Working_Support = relPosFilter->GetOutput();
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after SbS Rel pos Post");
-
-  // AutoCrop (OR SbS ?)
-  m_Working_Support = clitk::AutoCrop<MaskImageType>(m_Working_Support, GetBackgroundValue()); 
+  m_ListOfStations["8"] = m_Working_Support;  
+}
+//--------------------------------------------------------------------
 
-  writeImage<MaskImageType>(m_Working_Support, "step1.4.1.mhd");
 
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after autocrop");
+//--------------------------------------------------------------------
+template <class ImageType>
+void 
+clitk::ExtractLymphStationsFilter<ImageType>::
+ExtractStation_8_Left_Inf_Limits() 
+{
+  //--------------------------------------------------------------------
+  StartNewStep("[Station8] Left limits around esophagus with lines from VertebralBody-Aorta-Esophagus");
 
   // Estract slices for current support for slice by slice processing
   std::vector<typename MaskSliceType::Pointer> slices;
   clitk::ExtractSlices<MaskImageType>(m_Working_Support, 2, slices);
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after support slices");
-  
+
+  // Remove what is outside the mediastinum in this enlarged Esophagus -> it allows to select
+  // 'better' extrema points (not too post).
+  MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
+  clitk::AndNot<MaskImageType>(m_Esophagus, Lungs, GetBackgroundValue());
+  GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
+
   // Estract slices of Esophagus (resize like support before to have the same set of slices)
-  MaskImagePointer EsophagusForSlice = clitk::ResizeImageLike<MaskImageType>(Esophagus, m_Working_Support, GetBackgroundValue());
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after Eso resize");
+  MaskImagePointer EsophagusForSlice = 
+    clitk::ResizeImageLike<MaskImageType>(m_Esophagus, m_Working_Support, GetBackgroundValue());
   std::vector<typename MaskSliceType::Pointer> eso_slices;
   clitk::ExtractSlices<MaskImageType>(EsophagusForSlice, 2, eso_slices);
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after Eso slices");
 
   // Estract slices of Vertebral (resize like support before to have the same set of slices)
-  VertebralBody = GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after Read VertebralBody");
+  MaskImagePointer VertebralBody = GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
   VertebralBody = clitk::ResizeImageLike<MaskImageType>(VertebralBody, m_Working_Support, GetBackgroundValue());
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after VertebralBody Resize");
+  // Remove what is outside the support to not consider what is to
+  // posterior in the VertebralBody (post the horizontal line)
+  clitk::And<MaskImageType>(VertebralBody, m_Working_Support, GetBackgroundValue());
+  // writeImage<MaskImageType>(VertebralBody, "vb.mhd");
   std::vector<typename MaskSliceType::Pointer> vert_slices;
   clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, vert_slices);
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after VertebralBody slices");
+
+  // Estract slices of Aorta (resize like support before to have the same set of slices)
+  MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
+  Aorta = clitk::ResizeImageLike<MaskImageType>(Aorta, m_Working_Support, GetBackgroundValue());
+  std::vector<typename MaskSliceType::Pointer> aorta_slices;
+  clitk::ExtractSlices<MaskImageType>(Aorta, 2, aorta_slices);
 
   // Extract slices of Mediastinum (resize like support before to have the same set of slices)
   m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after read Mediastinum");
   m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Working_Support, GetBackgroundValue());
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after resize Mediastinum");
   std::vector<typename MaskSliceType::Pointer> mediast_slices;
   clitk::ExtractSlices<MaskImageType>(m_Mediastinum, 2, mediast_slices);
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after Mediastinum slices");
-
-  writeImage<MaskImageType>(EsophagusForSlice, "slices_eso.mhd");
-  writeImage<MaskImageType>(VertebralBody, "slices_vert.mhd");
-  writeImage<MaskImageType>(m_Mediastinum, "slices_medias.mhd");
-  writeImage<MaskImageType>(m_Working_Support, "slices_support.mhd");
-
-  // Find common slices between Eso and m_Working_Support
-  // int s=0;
-  // MaskImageIndexType z_Eso = Esophagus->GetLargestPossibleRegion().GetIndex();
-  // MaskImagePointType p_Eso; 
-  // Esophagus->TransformIndexToPhysicalPoint(z_Eso, p_Eso);
-  // MaskImageIndexType z_Support;  
-  // z_Support = m_Working_Support->GetLargestPossibleRegion().GetIndex();
-  // MaskImagePointType p_Support; 
-  // m_Working_Support->TransformIndexToPhysicalPoint(z_Support, p_Support);
-  // while (p_Eso[2] < p_Support[2]) {
-  //   z_Eso[2] ++;
-  //   Esophagus->TransformIndexToPhysicalPoint(z_Eso, p_Eso);    
-  // }
-  // s = z_Eso[2] - Esophagus->GetLargestPossibleRegion().GetIndex()[2];
-  // DD(s);
-
-  // Find common slices between m_Working_Support and Mediastinum
-  // int sm=0;
-  // MaskImageIndexType z_Mediast = m_Mediastinum->GetLargestPossibleRegion().GetIndex();
-  // MaskImagePointType p_Mediast; 
-  // m_Mediastinum->TransformIndexToPhysicalPoint(z_Mediast, p_Mediast);
-  // z_Support = m_Working_Support->GetLargestPossibleRegion().GetIndex();
-  // m_Working_Support->TransformIndexToPhysicalPoint(z_Support, p_Support);
-  // while (p_Mediast[2] < p_Support[2]) {
-  //   z_Mediast[2] ++;
-  //   m_Mediastinum->TransformIndexToPhysicalPoint(z_Mediast, p_Mediast);    
-  // }
-  // sm = z_Mediast[2] - m_Mediastinum->GetLargestPossibleRegion().GetIndex()[2];
-  // DD(sm);
-  
-  DD(EsophagusForSlice->GetLargestPossibleRegion().GetSize()[2]);
-  DD(m_Mediastinum->GetLargestPossibleRegion().GetSize()[2]);
-  DD(slices.size());
-  DD(vert_slices.size());
-  DD(eso_slices.size());
-  DD(mediast_slices.size());
 
-  // uint max = std::min(slices.size(), std::min(eso_slices.size()-s, mediast_slices.size()-sm));
-  // DD(max);
+  // List of points
+  std::vector<MaskImagePointType> p_MostLeftVertebralBody;
+  std::vector<MaskImagePointType> p_MostRightVertebralBody;
+  std::vector<MaskImagePointType> p_MostLeftAorta;
+  std::vector<MaskImagePointType> p_MostLeftEsophagus;
 
-  // DEBUG POINTS
-  std::ofstream osp;
-  openFileForWriting(osp, "point_s8.txt");
-  osp << "LANDMARKS1" << std::endl;
+  /*
+    In the following, we search for the LeftRight limits.  We search
+    for the most Right points in Esophagus and in VertebralBody and
+    consider a line between those to most right points. All points in
+    the support which are most right to this line are discarded. Same
+    for the left part. The underlying assumption is that the support
+    is concave between Eso/VertebralBody. Esophagus is a bit
+    dilatated. On VertebralBody we go right (or left) until we reach
+    the lung (but no more 20 mm).
+  */
 
-  // Loop slices
+  // Temporary 3D point
   MaskImagePointType p;
-  int pi=0;
-  for(uint i=0; i<slices.size() ; i++) {
-    DD(i);
-
-    // --------------------------------------------------------------------------
-    // Find the limit on the Right: most right point between Eso and
-    // Vertebra. (Right is left on screen, coordinate decrease)
-    typename MaskSliceType::PointType p_maxRight;
-    
-    // Find right limit of Esophagus
-    typename MaskSliceType::PointType p_maxRight_Eso;    
-    clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(eso_slices[i], GetBackgroundValue(), 
-                                                            1, false, p_maxRight_Eso);
-    // Debug point 
-    clitk::PointsUtils<MaskImageType>::Convert2DTo3D(p_maxRight_Eso, EsophagusForSlice, i, p);
-    osp << pi << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
-    ++pi;
-
-    // Find right limit of Vertebra
-    typename MaskSliceType::PointType p_maxRight_Vertebra;
-    clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vert_slices[i], GetBackgroundValue(), 
-                                                            1, false, p_maxRight_Vertebra);
-    // Debug point 
-    clitk::PointsUtils<MaskImageType>::Convert2DTo3D(p_maxRight_Vertebra, VertebralBody, i, p);
-    osp << pi << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
-    ++pi;
-    
-    // Find last point out of the mediastinum on this line :
-    typename MaskSliceType::IndexType index;
-    mediast_slices[i]->TransformPhysicalPointToIndex(p_maxRight_Vertebra, index);
-    double distance = 0.0;
-    while (mediast_slices[i]->GetPixel(index) != GetBackgroundValue()) {
-      index[0] --;
-      distance += mediast_slices[i]->GetSpacing()[0];
-    }
-    DD(distance);
-    if (distance < 20) { // Ok in this case, we found limit with lung
-      mediast_slices[i]->TransformIndexToPhysicalPoint(index, p_maxRight_Vertebra);
-      clitk::PointsUtils<MaskImageType>::Convert2DTo3D(p_maxRight_Vertebra, m_Mediastinum, i, p);
-      osp << pi << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
-      ++pi;
-    }
-    
-    // Choose the most extrema one
-    if (p_maxRight_Vertebra[0] < p_maxRight_Eso[0]) {
-      p_maxRight = p_maxRight_Vertebra;
-    }
-    else p_maxRight = p_maxRight_Eso;
 
-    // --------------------------------------------------------------------------
+  typename MaskSliceType::PointType minSlicePoint;
+  typename MaskSliceType::PointType maxSlicePoint;
+  clitk::GetMinMaxPointPosition<MaskSliceType>(vert_slices[0], minSlicePoint, maxSlicePoint);
 
-
-    /*
-    // Get most post of dilated Esophagus
-    typename MaskSliceType::PointType p_post;
-    bool f = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(eso_slices[i], GetBackgroundValue(), 1, false, p_post);
-    clitk::PointsUtils<MaskImageType>::Convert2DTo3D(p_post, EsophagusForSlice, i, p);
-    osp << pi << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
-    ++pi;
-    // DD(f);
+  // Loop slices
+  for(uint i=0; i<slices.size() ; i++) {
+    // Declare all needed 2D points (sp = slice point)
+    //    typename MaskSliceType::PointType sp_MostRightVertebralBody;
+    typename MaskSliceType::PointType sp_MostLeftVertebralBody;
+    typename MaskSliceType::PointType sp_MostLeftAorta;
+    typename MaskSliceType::PointType sp_temp;
+    typename MaskSliceType::PointType sp_MostLeftEsophagus;
+
+    /* -------------------------------------------------------------------------------------
+        Find first point not in mediastinum at LEFT
+    */
+    clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vert_slices[i], GetBackgroundValue(), 
+                                                            0, false, sp_MostLeftVertebralBody);
+    // clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostLeftVertebralBody, VertebralBody, i, p);
     // DD(p);
 
-    // Get most left of the vertebral body
-    typename MaskSliceType::PointType s_left;
-    f = f && clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vert_slices[i], GetBackgroundValue(), 0, false, s_left);
-    clitk::PointsUtils<MaskImageType>::Convert2DTo3D(s_left, VertebralBody, i, p);
-    osp << pi << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
-    ++pi;
-    // DD(f);
+    sp_MostLeftVertebralBody = 
+      clitk::FindExtremaPointInAGivenLine<MaskSliceType>(mediast_slices[i], 0, false, 
+                                                         sp_MostLeftVertebralBody, GetBackgroundValue(), 30);
+    // clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostLeftVertebralBody, VertebralBody, i, p);
     // DD(p);
 
-    // Find last point out of the mediastinum on this line :
-    typename MaskSliceType::IndexType index_left;
-    mediast_slices[i]->TransformPhysicalPointToIndex(s_left, index_left);
-    index_left[0] ++; // on more left to be inside the support
-    while (mediast_slices[i]->GetPixel(index_left) != GetBackgroundValue()) {
-    index_left[0] ++;
-    }
-    mediast_slices[i]->TransformIndexToPhysicalPoint(index_left, s_left);
-    clitk::PointsUtils<MaskImageType>::Convert2DTo3D(s_left, m_Mediastinum, i, p);
-    osp << pi << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
-    ++pi;
-    // DD(f);
+    sp_MostLeftVertebralBody[0] += 2; // 2mm margin
+    // clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostLeftVertebralBody, VertebralBody, i, p);
     // DD(p);
-    */
-
 
-    // Loop to suppress
-    // if (f) {
-    typename MaskSliceType::PointType p;
-    typedef itk::ImageRegionIteratorWithIndex<MaskSliceType> IteratorType;
-    IteratorType iter(slices[i], slices[i]->GetLargestPossibleRegion());
-    iter.GoToBegin();
-    while (!iter.IsAtEnd()) {
-      if (iter.Get() != GetBackgroundValue()) {
-        slices[i]->TransformIndexToPhysicalPoint(iter.GetIndex(), p);
+    // Convert 2D points into 3D
+    clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostLeftVertebralBody, VertebralBody, i, p);
+    p_MostLeftVertebralBody.push_back(p);
 
-        // Remove point too at RIGHT
-        if (p[0] < p_maxRight[0]) {
-          iter.Set(GetBackgroundValue());
-        }
+    /* -------------------------------------------------------------------------------------
+       Find first point not in mediastinum at RIGHT. Not used yet.
+    clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vert_slices[i], GetBackgroundValue(), 
+                                                            0, true, sp_MostRightVertebralBody);
+    sp_MostRightVertebralBody = 
+      clitk::FindExtremaPointInAGivenLine<MaskSliceType>(mediast_slices[i], 0, true, 
+                                                         sp_MostRightVertebralBody, GetBackgroundValue(),30);
+    sp_MostRightVertebralBody[0] -= 2; // 2 mm margin
+    
+    // Convert 2D points into 3D
+    clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostRightVertebralBody, VertebralBody, i, p);
+    p_MostRightVertebralBody.push_back(p);
+    */
 
-        /*
-        // Remove point from foreground if too right or to high
-        if ((p[1] > p_post[1]) && (p[0] > s_left[0])) {
-        iter.Set(GetBackgroundValue());
-        }
-        */
 
+    /* -------------------------------------------------------------------------------------
+       Find most Left point in Esophagus
+     */
+    bool found = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(eso_slices[i], GetBackgroundValue(), 0, false, sp_MostLeftEsophagus);
+    if (!found) { // No more Esophagus, I remove the previous point
+      //DD("no eso pop back");
+      p_MostLeftVertebralBody.pop_back();
+    }
+    else {
+      clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostLeftEsophagus, EsophagusForSlice, i, p);
+      p_MostLeftEsophagus.push_back(p);
+    }
+      
+    /* -------------------------------------------------------------------------------------
+       Find most Left point in Aorta 
+    */
+    if (found) {
+      clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(aorta_slices[i], GetBackgroundValue(), 0, false, sp_MostLeftAorta);
+      sp_MostLeftAorta = clitk::FindExtremaPointInAGivenLine<MaskSliceType>(mediast_slices[i], 0, false, sp_MostLeftAorta, GetBackgroundValue(), 10);
+      typename MaskSliceType::PointType temp=sp_MostLeftEsophagus;
+      temp[0] -= 10;
+      if (clitk::IsOnTheSameLineSide(sp_MostLeftAorta, sp_MostLeftVertebralBody, sp_MostLeftEsophagus, temp)) { 
+        // sp_MostLeftAorta is on the same side than temp (at Right) -> ignore Aorta
+        sp_MostLeftAorta = sp_MostLeftEsophagus;
       }
-      ++iter;
+      clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostLeftAorta, Aorta, i, p);
+      p_MostLeftAorta.push_back(p);
     }
-    // } // end if f
-    // s++;
-    // sm++;
-  }
-  osp.close();
-  
-
-  // Joint slices
-  DD("after loop");
-  m_Working_Support = clitk::JoinSlices<MaskImageType>(slices, m_Working_Support, 2);
-  clitk::PrintMemory(GetVerboseMemoryFlag(), "after JoinSlices");
 
-  // DEBUG
+  } // End of slice loop
+  
+  clitk::WriteListOfLandmarks<MaskImageType>(p_MostLeftVertebralBody, "S8-MostLeft-VB-points.txt");
+  //  clitk::WriteListOfLandmarks<MaskImageType>(p_MostRightVertebralBody, "S8-MostRight-VB-points.txt");
+  clitk::WriteListOfLandmarks<MaskImageType>(p_MostLeftAorta, "S8-MostLeft-Aorta-points.txt");
+  clitk::WriteListOfLandmarks<MaskImageType>(p_MostLeftEsophagus, "S8-MostLeft-eso-points.txt");
+
+  clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support, 
+                                                                    p_MostLeftVertebralBody, p_MostLeftAorta,
+                                                                    GetBackgroundValue(), 0, -10);
+
+  clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support, 
+                                                                    p_MostLeftAorta, p_MostLeftEsophagus,
+                                                                    GetBackgroundValue(), 0, -10);
+  // END
+  StopCurrentStep<MaskImageType>(m_Working_Support);
   m_ListOfStations["8"] = m_Working_Support;
   return;
 }
@@ -701,123 +462,88 @@ ExtractStation_8_AP_Limits()
 
 //--------------------------------------------------------------------
 template <class ImageType>
-typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
+void 
 clitk::ExtractLymphStationsFilter<ImageType>::
-EnlargeEsophagusDilatationRadiusInferiorly(MaskImagePointer & Esophagus)
+ExtractStation_8_Remove_Structures()
 {
-  // Check if Esophagus is delineated at least until GOjunction. Now,
-  // because we use AutoCrop, Origin[2] gives this max inferior
-  // position.
-  //  double GOjunctionZ = GetAFDB()->GetPoint3D("GOjunction", 2);
-
-  if (Esophagus->GetOrigin()[2] > m_DiaphragmInferiorLimit) {
-    std::cout << "Warning Esophagus is not delineated until Diaphragm. I mirror-pad it." 
-              << std::endl;
-    double extraSize = Esophagus->GetOrigin()[2]-m_DiaphragmInferiorLimit; 
-    
-    // Pad with few more slices
-    typedef itk::MirrorPadImageFilter<MaskImageType, MaskImageType> PadFilterType;
-    typename PadFilterType::Pointer padFilter = PadFilterType::New();
-    padFilter->SetInput(Esophagus);
-    MaskImageSizeType b;
-    b[0] = 0; b[1] = 0; b[2] = (uint)ceil(extraSize/Esophagus->GetSpacing()[2])+1;
-    padFilter->SetPadLowerBound(b);
-    padFilter->Update();
-    Esophagus = padFilter->GetOutput();
-  }
-  return Esophagus;
+  //--------------------------------------------------------------------
+  Remove_Structures("8", "Aorta");
+  Remove_Structures("8", "Esophagus");
+  Remove_Structures("8", "VertebralBody");  
+
+  // END
+  StopCurrentStep<MaskImageType>(m_Working_Support);
+  m_ListOfStations["8"] = m_Working_Support;
+  return;
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
-template <class ImageType>
+template <class TImageType>
 void 
-clitk::ExtractLymphStationsFilter<ImageType>::
-ExtractStation_8_LR_Limits() 
+clitk::ExtractLymphStationsFilter<TImageType>::
+FindExtremaPointsInBronchus(MaskImagePointer input, 
+                           int direction,
+                           double distance_max_from_center_point, 
+                           ListOfPointsType & LR, 
+                           ListOfPointsType & Ant, 
+                           ListOfPointsType & Post)
 {
-  /*
-    Station 8: paraeosphageal nodes
 
-    Laterally, it is within the pleural envelope and again abuts the
-    descending aorta on the left. Reasonably, the delineation of
-    Station 8 is limited to the soft tissue surrounding the esophagus
-    (Fig.  3C–H). 
-  */
+  // Other solution ==> with auto bounding box ! (but pb to prevent to
+  // be too distant from the center point
 
-  StartNewStep("[Station8] Right limits (around esophagus)");
-  // Get Esophagus
-  MaskImagePointer Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
-
-  // Autocrop to get first slice with starting Esophagus
-  Esophagus = clitk::AutoCrop<MaskImageType>(Esophagus, GetBackgroundValue()); 
-
-  // Dilate 
-  // LR dilatation -> large to keep point inside
-  // AP dilatation -> few mm 
-  // SI dilatation -> enough to cover Diaphragm (old=GOjunctionZ)
-  MaskImagePointType radiusInMM = GetEsophagusDiltationForRight();
-  Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(Esophagus);
-  Esophagus = clitk::Dilate<MaskImageType>(Esophagus, 
-                                           radiusInMM, 
-                                           GetBackgroundValue(), 
-                                           GetForegroundValue(), true);
-  writeImage<MaskImageType>(Esophagus, "dilateEso2.mhd");
-
-  writeImage<MaskImageType>(m_Working_Support, "before-relpos2.mhd");
-
-  // Remove Right (Left on the screen) part according to this
-  // dilatated esophagus
-  typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
-  typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
-  relPosFilter->VerboseStepFlagOff();
-  relPosFilter->WriteStepFlagOff();
-  relPosFilter->SetInput(m_Working_Support);
-  relPosFilter->SetInputObject(Esophagus);
-  relPosFilter->AddOrientationTypeString("L");
-  relPosFilter->InverseOrientationFlagOn(); // Not Left to
-  relPosFilter->SetDirection(2); // Z axis
-  relPosFilter->UniqueConnectedComponentBySliceOff(); // important
-  relPosFilter->SetIntermediateSpacing(4);
-  relPosFilter->ResampleBeforeRelativePositionFilterOn();
-  relPosFilter->SetFuzzyThreshold(0.9); // remove few part only
-  relPosFilter->RemoveObjectFlagOff();
-  relPosFilter->Update();
-  m_Working_Support = relPosFilter->GetOutput();
-
-  // Get a single 3D CCL
-  m_Working_Support = Labelize<MaskImageType>(m_Working_Support, 0, false, 10);
-  m_Working_Support = KeepLabels<MaskImageType>(m_Working_Support, 
-                                                GetBackgroundValue(), 
-                                                GetForegroundValue(), 1, 1, true);
-
-
-  /*
-  // Re-Add post to Esophagus -> sometimes previous relpos remove some
-  // correct part below esophagus.
-  MaskImagePointer Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
-  EnlargeEsophagusDilatationRadiusInferiorly(Esophagus);
-  writeImage<MaskImageType>(Esophagus, "e-again.mhd");
-  typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
-  typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
-  relPosFilter->VerboseStepOff();
-  relPosFilter->WriteStepOff();
-  relPosFilter->SetInput(m_Working_Support);
-  relPosFilter->SetInputObject(Esophagus);
-  relPosFilter->SetOrientationTypeString("P");
-  relPosFilter->InverseOrientationFlagOff(); 
-  relPosFilter->SetDirection(2); // Z axis
-  relPosFilter->UniqueConnectedComponentBySliceOff(); // important
-  relPosFilter->SetIntermediateSpacing(4);
-  relPosFilter->ResampleBeforeRelativePositionFilterOn();
-  relPosFilter->CombineWithOrFlagOn();
-  relPosFilter->SetFuzzyThreshold(0.9); // remove few part only
-  relPosFilter->RemoveObjectFlagOff();
-  relPosFilter->Update();
-  m_Working_Support = relPosFilter->GetOutput();
-  */
+  // Extract slices
+  std::vector<typename MaskSliceType::Pointer> slices;
+  clitk::ExtractSlices<MaskImageType>(input, 2, slices);
+  
+  // Loop on slices
+  bool found;
+  for(uint i=0; i<slices.size(); i++) {
+    /*
+    // Keep main CCL
+    slices[i] = Labelize<MaskSliceType>(slices[i], 0, true, 10);
+    slices[i] = KeepLabels<MaskSliceType>(slices[i], 
+    GetBackgroundValue(), 
+    GetForegroundValue(), 1, 1, true);
+    */
 
-  StopCurrentStep<MaskImageType>(m_Working_Support);
-  m_ListOfStations["8"] = m_Working_Support;
-}
+    // ------- Find rightmost or leftmost point  ------- 
+    MaskSliceType::PointType LRMost;
+    found = 
+      clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices[i], 
+                                                              GetBackgroundValue(), 
+                                                              0, // axis XY
+                                                              (direction==0?false:true),  // right or left according to direction
+                                                              LRMost);
+    // ------- Find postmost point  ------- 
+    MaskSliceType::PointType postMost;
+    found = 
+      clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices[i], 
+                                                              GetBackgroundValue(), 
+                                                              1, false, LRMost, 
+                                                              distance_max_from_center_point, 
+                                                              postMost);
+    // ------- Find antmost point  ------- 
+    MaskSliceType::PointType antMost;
+    found = 
+      clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices[i], 
+                                                              GetBackgroundValue(), 
+                                                              1, true, LRMost, 
+                                                              distance_max_from_center_point, 
+                                                              antMost);
+    // Only add point if found
+    if (found)  {
+      // ------- Convert 2D to 3D points --------
+      MaskImageType::PointType p;
+      clitk::PointsUtils<MaskImageType>::Convert2DTo3D(LRMost, input, i, p);
+      LR.push_back(p); 
+      clitk::PointsUtils<MaskImageType>::Convert2DTo3D(antMost, input, i, p);
+      Ant.push_back(p);
+      clitk::PointsUtils<MaskImageType>::Convert2DTo3D(postMost, input, i, p);
+      Post.push_back(p);
+    }
+  }
+} 
 //--------------------------------------------------------------------