#include <itkBinaryDilateImageFilter.h>
#include <itkMirrorPadImageFilter.h>

//--------------------------------------------------------------------
template <class ImageType>
void 
clitk::ExtractLymphStationsFilter<ImageType>::
ExtractStation_8_SI_Limits() 
{
  /*
    Station 8: paraeosphageal nodes
    
    "Superiorly, Station 8 begins at the level of the carina and,
    therefore, abuts Station 3P (Fig. 3C). Inferiorly, the lower limit
    of Station 8 was not specified in the Mountain and Dresler
    classification (1). We delineate this volume until it reaches the
    gastroesphogeal junction. "
    
    => new classification IASCL 2009: 

    "Lower border: the diaphragm"
    
    "Upper border: the upper border of the lower lobe bronchus on the
    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

  /*
    double m_CarinaZ = GetAFDB()->GetPoint3D("Carina", 2) + 
    m_Mediastinum->GetSpacing()[2]; // add one slice to include carina
    // double GOjunctionZ = GetAFDB()->GetPoint3D("GOjunction", 2);
    */

  // 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) */
  m_Working_Support = 
    clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2, 
                                                m_DiaphragmInferiorLimit, //GOjunctionZ, 
                                                m_CarinaZ, 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;
}
//--------------------------------------------------------------------


//--------------------------------------------------------------------
template <class ImageType>
void 
clitk::ExtractLymphStationsFilter<ImageType>::
ExtractStation_8_AP_Limits() 
{
  /*
    Station 8: paraeosphageal nodes

    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). 

    New classification IASCL 2009 :

    "Nodes lying adjacent to the wall of the esophagus and to the
    right or left of the midline, excluding subcarinal nodes (S7)
    Upper"

  */

  StartNewStep("[Station8] Post limits with vertebral body");
  MaskImagePointer VertebralBody = 
    GetAFDB()->template GetImage <MaskImageType>("VertebralBody");  

  // 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);
    if (found) {
      vertebralAntPositionBySlice[i] = p;
    }
    else { 
      // no Foreground in this slice (it appends generally at the
      // beginning and the end of the volume). Do nothing in this
      // case.
    }
  }

  // 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;
      }
    }
  }

  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]);  
  */

  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;
        }
        siter.NextLine();
      }
      ++i;
    }
    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);
  
  // 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()); 

  writeImage<MaskImageType>(m_Working_Support, "step1.4.1.mhd");

  clitk::PrintMemory(GetVerboseMemoryFlag(), "after autocrop");

  // 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");
  
  // 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");
  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");
  VertebralBody = clitk::ResizeImageLike<MaskImageType>(VertebralBody, m_Working_Support, GetBackgroundValue());
  clitk::PrintMemory(GetVerboseMemoryFlag(), "after VertebralBody Resize");
  std::vector<typename MaskSliceType::Pointer> vert_slices;
  clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, vert_slices);
  clitk::PrintMemory(GetVerboseMemoryFlag(), "after VertebralBody 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);

  // DEBUG POINTS
  std::ofstream osp;
  openFileForWriting(osp, "point_s8.txt");
  osp << "LANDMARKS1" << std::endl;

  // Loop slices
  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;

    // --------------------------------------------------------------------------


    /*
    // 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);
    // 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);
    // 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);
    // 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);

        // Remove point too at RIGHT
        if (p[0] < p_maxRight[0]) {
          iter.Set(GetBackgroundValue());
        }

        /*
        // Remove point from foreground if too right or to high
        if ((p[1] > p_post[1]) && (p[0] > s_left[0])) {
        iter.Set(GetBackgroundValue());
        }
        */

      }
      ++iter;
    }
    // } // 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
  m_ListOfStations["8"] = m_Working_Support;
  return;
}
//--------------------------------------------------------------------


//--------------------------------------------------------------------
template <class ImageType>
typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
clitk::ExtractLymphStationsFilter<ImageType>::
EnlargeEsophagusDilatationRadiusInferiorly(MaskImagePointer & Esophagus)
{
  // 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;
}
//--------------------------------------------------------------------


//--------------------------------------------------------------------
template <class ImageType>
void 
clitk::ExtractLymphStationsFilter<ImageType>::
ExtractStation_8_LR_Limits() 
{
  /*
    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). 
  */

  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();
  */

  StopCurrentStep<MaskImageType>(m_Working_Support);
  m_ListOfStations["8"] = m_Working_Support;
}
//--------------------------------------------------------------------