add options ; "last chance" etc

diff --git a/segmentation/clitkExtractMediastinalVesselsFilter.txx b/segmentation/clitkExtractMediastinalVesselsFilter.txx
index 02576f498289c8ef0b3e04e9dadb6fd94fee64ca..2198529b8c60cd667e1362ae0ab9767ceaa8f8dd 100644 (file)
--- a/segmentation/clitkExtractMediastinalVesselsFilter.txx
+++ b/segmentation/clitkExtractMediastinalVesselsFilter.txx
@@ -22,13 +22,11 @@
 // clitk
 #include "clitkCommon.h"
 #include "clitkExtractMediastinalVesselsFilter.h"
-#include "clitkAutoCropFilter.h"
 #include "clitkSegmentationUtils.h"
-#include "clitkMorphoMathFilter.h"
+#include "clitkReconstructWithConditionalGrayscaleDilateImageFilter.h"
 
 // itk
 #include <itkBinaryThresholdImageFilter.h>
-#include <itkGrayscaleDilateImageFilter.h>
 #include <itkMinimumMaximumImageCalculator.h>
 
 //--------------------------------------------------------------------
@@ -42,8 +40,15 @@ ExtractMediastinalVesselsFilter():
   this->SetNumberOfRequiredInputs(1);
   SetBackgroundValue(0);
   SetForegroundValue(1);
-  SetThreshold(140);
-  SetTemporaryForegroundValue(1);
+  SetThresholdHigh(140);
+  SetThresholdLow(55);
+  SetErosionRadius(2);
+  SetDilatationRadius(9);
+  SetMaxDistancePostToCarina(10);
+  SetMaxDistanceAntToCarina(40);
+  SetMaxDistanceLeftToCarina(35);
+  SetMaxDistanceRightToCarina(35);
+  SetSoughtVesselSeedName("NoSeedNameGiven");
 }
 //--------------------------------------------------------------------
 
@@ -58,145 +63,149 @@ GenerateOutputInformation() {
   m_Input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
   
   // ------------------------------------------------------------------
-  // Sup-Inf crop -> Carina
-  CropSupInf();  
+  // Crop the initial image superiorly and inferiorly. 
+  // TODO : add options for x cm above/below
+  CropInputImage();  
 
   // ------------------------------------------------------------------
-  // Binarize
-  StartNewStep("Binarize with treshold = "+toString(GetThreshold()));
+  // Binarize the image. Need two thresholds, one high to select
+  // structures (CCL) that are almost not connected (after erosion),
+  // and one low thresholds to select the real contours. Will be
+  // reconstructed later. 
+  StartNewStep("Binarize with high threshold = "+toString(GetThresholdHigh()));
   typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> BinarizeFilterType; 
   typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
   binarizeFilter->SetInput(m_Input);
-  binarizeFilter->SetLowerThreshold(GetThreshold());
-  binarizeFilter->SetInsideValue(GetTemporaryForegroundValue());
+  binarizeFilter->SetLowerThreshold(GetThresholdHigh());
+  binarizeFilter->SetInsideValue(GetForegroundValue());
   binarizeFilter->SetOutsideValue(GetBackgroundValue());
   binarizeFilter->Update();
   m_Mask = binarizeFilter->GetOutput();
-  clitk::writeImage<MaskImageType>(m_Mask, "m.mhd");
   StopCurrentStep<MaskImageType>(m_Mask);
 
-  // Keep main CCL ? 
-  m_Mask = clitk::Labelize<MaskImageType>(m_Mask, GetBackgroundValue(), false, 10);
-  m_Mask = KeepLabels<MaskImageType>(m_Mask, GetBackgroundValue(), GetTemporaryForegroundValue(), 1, 1, true);
-  clitk::writeImage<MaskImageType>(m_Mask, "m2.mhd");
-  
+  // ------------------------------------------------------------------
+  StartNewStep("Binarize with low threshold = "+toString(GetThresholdLow()));
+  binarizeFilter = BinarizeFilterType::New();
+  binarizeFilter->SetInput(m_Input);
+  binarizeFilter->SetLowerThreshold(GetThresholdLow());
+  binarizeFilter->SetInsideValue(GetForegroundValue());
+  binarizeFilter->SetOutsideValue(GetBackgroundValue());
+  binarizeFilter->Update();
+  MaskImagePointer m_Mask2 = binarizeFilter->GetOutput();
+  StopCurrentStep<MaskImageType>(m_Mask2);
+
   // ------------------------------------------------------------------
   // Extract slices
-  StartNewStep("Detect vessels (slice by slice reconstruction)");
+  StartNewStep("Detect objects : erosion, then slice by slice reconstruction");
   std::vector<MaskSlicePointer> slices_mask;
   clitk::ExtractSlices<MaskImageType>(m_Mask, 2, slices_mask);
-  
-  DD(slices_mask.size());
-  
+  std::vector<MaskSlicePointer> slices_mask2;
+  clitk::ExtractSlices<MaskImageType>(m_Mask2, 2, slices_mask2);
+  int radius = GetErosionRadius();
+
+  // List of working slices (debug only)
   std::vector<MaskSlicePointer> debug_eroded;
-  std::vector<MaskSlicePointer> debug_labeled;
-  std::vector<MaskSlicePointer> debug_slabeled;
+  // std::vector<MaskSlicePointer> debug_labeled;
   
-  int radius = 3;
-  DD(radius); // TO PUT IN OPTION
-
   // ------------------------------------------------------------------
-  // Loop Slice by Slice -> erode find CCL and reconstruct
-  clitk::MorphoMathFilter<MaskSliceType>::Pointer f= clitk::MorphoMathFilter<MaskSliceType>::New();
+  // Loop Slice by Slice in order to break connectivity between
+  // CCL. Erode and reconstruct all labels at the same time without
+  // merging them.
   for(uint i=0; i<slices_mask.size(); i++) {
-    // Erosion
-    f->SetInput(slices_mask[i]);
-    f->SetBackgroundValue(GetBackgroundValue());
-    f->SetForegroundValue(GetTemporaryForegroundValue());
-    f->SetRadius(radius);
-    f->SetOperationType(0); // Erode
-    f->VerboseFlagOff();
-    f->Update();
-    MaskSlicePointer eroded = f->GetOutput();
-    //    writeImage<MaskSliceType>(eroded, "er-"+toString(i)+".mhd");
+    // Erosion kernel
+    typedef itk::BinaryBallStructuringElement<MaskSliceType::PixelType,2> KernelType;
+    KernelType structuringElement;
+    structuringElement.SetRadius(radius);
+    structuringElement.CreateStructuringElement();
+    
+    // Erosion -> we break the connectivity between structure
+    typedef itk::BinaryErodeImageFilter<MaskSliceType, MaskSliceType, KernelType> ErodeFilterType;
+    typename ErodeFilterType::Pointer eroder = ErodeFilterType::New();
+    eroder->SetInput(slices_mask[i]);
+    eroder->SetBackgroundValue(GetBackgroundValue());
+    eroder->SetForegroundValue(GetForegroundValue());
+    eroder->SetBoundaryToForeground(true); // ??
+    eroder->SetKernel(structuringElement);
+    eroder->Update();
+    MaskSlicePointer eroded = eroder->GetOutput();
     debug_eroded.push_back(eroded);
       
-    // CCL
-    int nb;
+    // Labelize (CCL)
     MaskSlicePointer labeled = 
-      clitk::LabelizeAndCountNumberOfObjects<MaskSliceType>(eroded, GetBackgroundValue(), false, 1, nb);
-
-    // Relabel, large CCL with large label number
-    for(int n=nb; n>0; n--) {
-      //        DD(n);
-      int li = n;
-      int lo = 2*(nb+1)-li;
-      labeled = clitk::SetBackground<MaskSliceType, MaskSliceType>(labeled, labeled, li, lo, true);
-    }
-    debug_labeled.push_back(labeled);
-
-    // Create kernel for GrayscaleDilateImageFilter
-    typedef itk::BinaryBallStructuringElement<MaskSliceType::PixelType,MaskSliceType::ImageDimension > KernelType;
-    KernelType k;
-    k.SetRadius(radius+1);
-    k.CreateStructuringElement();
-    
-    // Keep the MAX -> we prefer the opposite su change the label
-    typedef itk::GrayscaleDilateImageFilter<MaskSliceType, MaskSliceType, KernelType> FilterType;
-    FilterType::Pointer m = FilterType::New();
-    m->SetKernel(k);
-    m->SetInput(labeled);
-    // DD(m->GetAlgorithm());
-    // m->SetAlgorithm(3);
-    m->Update();
-    MaskSlicePointer s = m->GetOutput();
-
-
-    // Remove Initial BG
-    s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, slices_mask[i], 
-                                                   GetBackgroundValue(), GetBackgroundValue(), true);
-
+      clitk::Labelize<MaskSliceType>(eroded, GetBackgroundValue(), true, 1); // Fully connected !
+    // debug_labeled.push_back(labeled);
+
+    // Make Reconstruction filter : dilation all labels at the same
+    // time, prevent to merge them.
+    typedef clitk::ReconstructWithConditionalGrayscaleDilateImageFilter<MaskSliceType> ReconstructFilterType;
+    typename ReconstructFilterType::Pointer reconstructor = ReconstructFilterType::New();
+    reconstructor->SetInput(labeled);
+    reconstructor->SetIterationNumber(radius+GetDilatationRadius());
+    reconstructor->Update();
+    MaskSlicePointer s = reconstructor->GetOutput();
+
+    // Remove Initial BG of the second tresholded image
+    s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, slices_mask2[i], 
+                                                           GetBackgroundValue(), GetBackgroundValue(), true);
     m_slice_recon.push_back(s);
+
   } // end loop
-  DD("end loop");
   
+  // Build 3D images from the slice by slice processing
   MaskImageType::Pointer eroded = clitk::JoinSlices<MaskImageType>(debug_eroded, m_Mask, 2);
-  clitk::writeImage<MaskImageType>(eroded, "eroded.mhd");
-
-  DD("l");
-  MaskImageType::Pointer l = clitk::JoinSlices<MaskImageType>(debug_labeled, m_Mask, 2);
-  clitk::writeImage<MaskImageType>(l, "labeled.mhd");
-  
-  DD("r");
+  writeImage<MaskImageType>(eroded, "erode.mhd");
+  //MaskImageType::Pointer l = clitk::JoinSlices<MaskImageType>(debug_labeled, m_Mask, 2);  
   MaskImageType::Pointer r = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
-  clitk::writeImage<MaskImageType>(r, "recon.mhd");
+  writeImage<MaskImageType>(r, "recon1.mhd");
   
   // ------------------------------------------------------------------
-  // Loop Slice by Slice -> BCA not found yet
-  /*  MaskImagePointType BCA_p;
-  GetAFDB()->GetPoint3D("BrachioCephalicArteryFirstInferiorPoint", BCA_p);
-  DD(BCA_p);
-  MaskImagePointType bif1;
-  MaskImagePointType bif2;
-  TrackBifurcationFromPoint(r, BCA_p, bif1, bif2);
-  DD(bif1);
-  DD(bif2);
-  */
-  // Find max label
+  // Track the SoughtVessel from the given first point
+  // superiorly. This is done by TrackBifurcationFromPoint
+  MaskImagePointType SoughtVesselSeedPoint;
+  GetAFDB()->GetPoint3D(m_SoughtVesselSeedName, SoughtVesselSeedPoint);
+
+  // Find the label with the maximum value to set the result
   typedef itk::MinimumMaximumImageCalculator<MaskImageType> MinMaxFilterType;
   MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
   ff->SetImage(r);
   ff->ComputeMaximum();
   LabelType newLabel = ff->GetMaximum()+1; 
-  DD(newLabel);
 
-  // Get all centroids of the first slice
-  std::vector<MaskSlicePointType> centroids2D;
-  clitk::ComputeCentroids<MaskSliceType>(m_slice_recon[0], GetBackgroundValue(), centroids2D);
-  DD(centroids2D.size());
+  // the following bifurcations point will the centroids of the
+  // components obtain when (hopefully!) the SoughtVessel
+  // split into CommonArtery and SubclavianArtery.
   std::vector<MaskImagePointType> bifurcations;
-  clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(centroids2D, 0, r, bifurcations);  
-  DD(bifurcations.size());
-  for(uint i=1; i<bifurcations.size()+1; i++) {
-    DD(i);
-    DD(bifurcations.size());
-    TrackBifurcationFromPoint(r, m_slice_recon, bifurcations[i], newLabel+i, bifurcations);
-    DD("end track");
-    DD(bifurcations.size());
-    MaskImageType::Pointer rr = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
-    clitk::writeImage<MaskImageType>(rr, "recon"+toString(i)+".mhd");
-  }
+  TrackBifurcationFromPoint(r, m_slice_recon, SoughtVesselSeedPoint, newLabel, bifurcations);
+
+  // Build the final 3D image from the previous slice by slice processing
+  m_SoughtVessel = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
+  writeImage<MaskImageType>(m_SoughtVessel, "recon2.mhd");
+  
+  // Set binary image, (remove other labels).  
+  // TODO: keep labeled image to track SubclavianArtery and CommonArtery
+  m_SoughtVessel = 
+    clitk::Binarize<MaskImageType>(m_SoughtVessel, newLabel, newLabel, 
+                                   GetBackgroundValue(), GetForegroundValue());
+
+  writeImage<MaskImageType>(m_SoughtVessel, "afterbinarize.mhd");
+
+  m_SoughtVessel = clitk::AutoCrop<MaskImageType>(m_SoughtVessel, GetBackgroundValue());
 
+  writeImage<MaskImageType>(m_SoughtVessel, "afterautocrop.mhd");
+
+  // Clean the image : Opening (not in Z direction)
+  typename MaskImageType::SizeType rad;
+  rad[0] = rad[1] = 2;
+  rad[2] = 0;
+  m_SoughtVessel = clitk::Opening<MaskImageType>(m_SoughtVessel, rad, 
+                                                          GetBackgroundValue(), GetForegroundValue());
+
+  writeImage<MaskImageType>(m_SoughtVessel, "afteropen.mhd");
+
+  // Clean the image : keep main CCL slice by slice
+  m_SoughtVessel = clitk::SliceBySliceKeepMainCCL<MaskImageType>(m_SoughtVessel, 
+                                                                          GetBackgroundValue(), 
+                                                                          GetForegroundValue());  
 }
 //--------------------------------------------------------------------
 
@@ -216,11 +225,11 @@ template <class TImageType>
 void 
 clitk::ExtractMediastinalVesselsFilter<TImageType>::
 GenerateData() {
-  DD("GenerateData");
-  // Final Step -> graft output (if SetNthOutput => redo)
-  MaskImagePointer BrachioCephalicArtery = 
-    GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
-  this->GraftNthOutput(0, BrachioCephalicArtery);
+  // Save in the AFDB (not write on the disk here)
+  GetAFDB()->SetImageFilename(GetSoughtVesselName(), GetOutputFilename());  
+  WriteAFDB();
+  // Final Step -> graft output
+  this->GraftNthOutput(0, m_SoughtVessel);
 }
 //--------------------------------------------------------------------
   
@@ -229,40 +238,52 @@ GenerateData() {
 template <class TImageType>
 void 
 clitk::ExtractMediastinalVesselsFilter<TImageType>::
-CropSupInf() { 
-  StartNewStep("Inf/Sup limits (carina) and crop with mediastinum");
+CropInputImage() { 
+  StartNewStep("Crop the input image: SI,AP limits with carina and crop with mediastinum");
+  /*
+    Need : Trachea, Carina (roi not point), 
+   */
   // Get Trachea and Carina
   MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");  
   
-  // Get or compute Carina
+  // Compute Carina position
   double m_CarinaZ;
-  // Get Carina Z position
   MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
-  
   std::vector<MaskImagePointType> centroids;
   clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
   m_CarinaZ = centroids[1][2];
-  // DD(m_CarinaZ);
-  // add one slice to include carina ?
+  // add one slice to include carina 
   m_CarinaZ += Carina->GetSpacing()[2];
   // We dont need Carina structure from now
   Carina->Delete();
-  GetAFDB()->SetDouble("CarinaZ", m_CarinaZ);
+  GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
   
-  // Crop
-  m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 2, 
-                                                       m_CarinaZ, false, GetBackgroundValue());  
+  // Crop Inf, remove below Carina
+  m_Input = 
+    clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 2, m_CarinaZ, false, GetBackgroundValue());  
 
-  // Crop not post to centroid
+  // Crop post
   double m_CarinaY = centroids[1][1];
-  DD(m_CarinaY);
-  m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 1, // OLD ABOVE
-                                                         m_CarinaY, false, GetBackgroundValue());  
-  // Crop not ant to centroid
+  m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 1,
+                                                         m_CarinaY+GetMaxDistancePostToCarina(), 
+                                                         false, GetBackgroundValue());  
+  // Crop ant 
   m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 1, 
-                                                         m_CarinaY-80, false, GetBackgroundValue());  
-
-  // AutoCrop with Mediastinum
+                                                       m_CarinaY-GetMaxDistanceAntToCarina(), 
+                                                       false, GetBackgroundValue());  
+
+  // Crop Right
+  double m_CarinaX = centroids[1][0];
+  m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 0, 
+                                                       m_CarinaX-GetMaxDistanceRightToCarina(), 
+                                                       false, GetBackgroundValue());  
+  // Crop Left
+  m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 0, 
+                                                         m_CarinaX+GetMaxDistanceLeftToCarina(), 
+                                                         false, GetBackgroundValue());  
+
+  /*
+  // AutoCrop with Mediastinum, generaly only allow to remove few slices (superiorly)
   m_Mediastinum  = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
   // Resize like input (sup to carina)
   m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Input, GetBackgroundValue());
@@ -270,6 +291,7 @@ CropSupInf() {
   m_Mediastinum = clitk::AutoCrop<MaskImageType>(m_Mediastinum, GetBackgroundValue());
   // Resize input
   m_Input = clitk::ResizeImageLike<ImageType>(m_Input, m_Mediastinum, GetBackgroundValue());
+  */
 
   // End
   StopCurrentStep<ImageType>(m_Input);
@@ -281,39 +303,36 @@ CropSupInf() {
 template <class TImageType>
 void 
 clitk::ExtractMediastinalVesselsFilter<TImageType>::
-//SearchBrachioCephalicArtery(int & BCA_first_slice, LabelType & BCA_first_label) { 
 TrackBifurcationFromPoint(MaskImagePointer & recon, 
                           std::vector<MaskSlicePointer> & slices_recon, 
                           MaskImagePointType point3D, 
                           LabelType newLabel,
                           std::vector<MaskImagePointType> & bifurcations) {
-  StartNewStep("Search for BCA first slice and label");
-  DD(newLabel);
+  StartNewStep("Track the SoughtVessel from the seed point");
 
-  // Extract slices
-  //  std::vector<MaskSlicePointer> slices_recon;
-  //clitk::ExtractSlices<MaskImageType>(recon, 2, slices_recon);
-
-  // Find first slice
+  // Find first slice index
   MaskImageIndexType index;
   recon->TransformPhysicalPointToIndex(point3D, index);
-  DD(point3D);
-  DD(index);
+  int numberOfBifurcation = 0;
+  typedef typename MaskSliceType::PointType SlicePointType;
+  SlicePointType previousCenter;
 
-  uint i=index[2]; 
+  // Get current label at the point3D of interest
+  uint currentSlice=index[2]; 
   bool found = false;
-  LabelType previous_largest_label=recon->GetPixel(index);
-  DD(previous_largest_label);
+  LabelType previous_slice_label=recon->GetPixel(index);
+  //  DD(slices_recon.size());
   do {
-    DD(i);
+    //    DD(currentSlice);
     // Consider current reconstructed slice
-    MaskSlicePointer s = slices_recon[i];
+    MaskSlicePointer s = slices_recon[currentSlice];
     MaskSlicePointer previous;
-    if (i==index[2]) previous = s;
-    else previous = slices_recon[i-1];
+    if (currentSlice == index[2]) previous = s;
+    else {
+      previous = slices_recon[currentSlice-1];
+    }
     
     // Get centroids of the labels in the current slice
-    typedef typename MaskSliceType::PointType SlicePointType;
     static const unsigned int Dim = MaskSliceType::ImageDimension;
     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
     typedef itk::LabelMap< LabelObjectType > LabelMapType;
@@ -330,73 +349,128 @@ TrackBifurcationFromPoint(MaskImagePointer & recon,
     // Look what centroid inside the previous largest one
     std::vector<SlicePointType> centroids;
     std::vector<LabelType> centroids_label;
+    std::vector<double> labels_size;
     for(uint c=0; c<labelMap->GetNumberOfLabelObjects(); c++) {
       int label = labelMap->GetLabels()[c];
-      DD(label);
+      //      DD(label);
       SlicePointType center = labelMap->GetLabelObject(label)->GetCentroid();
-      SlicePointType center_previous = center;
-      center_previous[2] -= m_Input->GetSpacing()[2];      
+      //      DD(center);
       // Get label into previous slice
-      typename MaskSliceType::IndexType index;
-      previous->TransformPhysicalPointToIndex(center_previous, index);
-      LabelType l = previous->GetPixel(index);
-      DD(l);
-      if (l == previous_largest_label) {
+      typename MaskSliceType::IndexType centerIndex;
+      previous->TransformPhysicalPointToIndex(center, centerIndex);
+      LabelType labelInPreviousSlice = previous->GetPixel(centerIndex);
+      // if this current centroid was in the current label, add it to bifurcations
+      if (labelInPreviousSlice == previous_slice_label) {
         centroids.push_back(center);
         centroids_label.push_back(label);
+        labels_size.push_back(labelMap->GetLabelObject(label)->GetPhysicalSize());
+      }
+    }
+
+    if (centroids.size() == 0) {
+      // Last attempt to find -> check if previous centroid is inside a CCL
+      //      if in s -> get value, getcentroid add.
+      DD(currentSlice);
+      DD("Last change to find");
+      typename MaskSliceType::IndexType previousCenterIndex;
+      s->TransformPhysicalPointToIndex(previousCenter, previousCenterIndex);
+      DD(previousCenter);
+      LabelType labelInSlice = s->GetPixel(previousCenterIndex);
+      DD(labelInSlice);
+      if (labelInSlice != GetBackgroundValue()) {
+        centroids.push_back(labelMap->GetLabelObject(labelInSlice)->GetCentroid());
+        centroids_label.push_back(labelInSlice);
+        labels_size.push_back(labelMap->GetLabelObject(labelInSlice)->GetPhysicalSize());
       }
     }
-    DD(centroids.size());
+
+    
+    //    DD(centroids.size());
     
     // If several centroids, we found a bifurcation
     if (centroids.size() > 1) {
-      found = true;
-      for(uint c=0; c<centroids.size(); c++) {
-        ImagePointType bif;
-        clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[c], m_Mask, i, bif);
-        bifurcations.push_back(bif);
-        s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, centroids_label[c], newLabel+c+1, true);
-        //        slices_recon[i] = s; // (useful ?)
+      numberOfBifurcation++;
+      // If the number of bifurcation is greater than the required one, we stop
+      if (numberOfBifurcation > GetMaxNumberOfFoundBifurcation()) {
+        found = true;
+        DD("max bif reach");
+        for(uint c=0; c<centroids.size(); c++) {
+          ImagePointType bif;
+          clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[c], m_Mask, currentSlice, bif);
+          bifurcations.push_back(bif);
+        }
+      }
+      // Else we continue along the main (largest) connected component
+      else {
+        int indexOfLargest = 0;
+        for(uint b=0; b<centroids.size(); b++) {
+          if (labels_size[b] > labels_size[indexOfLargest]) {
+            indexOfLargest = b;
+          }
+        }
+        SlicePointType c = centroids[indexOfLargest];
+        LabelType l = centroids_label[indexOfLargest];
+        centroids.clear();
+        centroids.push_back(c);
+        centroids_label.push_back(l);
       }
-      DD("FOUND");
     }
     
     // if only one centroids, we change the current image with the current label 
     if (centroids.size() == 1) {
       s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, centroids_label[0], newLabel, true);
-      previous_largest_label = newLabel;
-      /*typedef itk::BinaryThresholdImageFilter<MaskSliceType, MaskSliceType> BinarizeFilterType; 
-      typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
-      binarizeFilter->SetInput(s);
-      binarizeFilter->SetLowerThreshold(centroids_label[0]);
-      binarizeFilter->SetUpperThreshold(centroids_label[0]+1);
-      binarizeFilter->SetInsideValue(previous_largest_label);
-      binarizeFilter->SetOutsideValue(GetBackgroundValue());
-      binarizeFilter->Update();
-      s = binarizeFilter->GetOutput();*/
-      slices_recon[i] = s; // (not useful ?)
+      slices_recon[currentSlice] = s;
+      previous_slice_label = newLabel;
+      // It can happend that several CCL share this same label. To
+      // prevent this case, we only consider the one that contains
+      // the centroid. 
+      MaskSlicePointer temp = clitk::Binarize<MaskSliceType>(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue());
+      //      writeImage<MaskSliceType>(temp, "bin-"+toString(currentSlice)+".mhd");
+      temp = clitk::Labelize<MaskSliceType>(temp, GetBackgroundValue(), true, 1);
+      //writeImage<MaskSliceType>(temp, "label-"+toString(currentSlice)+".mhd");
+      typename MaskSliceType::IndexType centroids_index;
+      temp->TransformPhysicalPointToIndex(centroids[0], centroids_index);
+      typename MaskSliceType::PixelType v = temp->GetPixel(centroids_index); 
+
+      // It can happend that the centroid is inside the BG, so we keep
+      // the largest CCL (the first);
+      if (v == GetBackgroundValue()) {
+        DD(currentSlice);
+        DD("inside BG");
+        DD(centroids[0]);
+        v = 1; // largest one
+      }
+
+      //DD(v);
+      temp = clitk::Binarize<MaskSliceType>(temp, v, v, GetBackgroundValue(), newLabel);      
+      //writeImage<MaskSliceType>(temp, "relabel-"+toString(currentSlice)+".mhd");      
+      s = temp;
+      slices_recon[currentSlice] = s;
+
+      // I need to recompute the centroid if we have removed some
+      // connected component.
+      clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), centroids);
+      previousCenter = centroids[1];
     }
 
     if (centroids.size() == 0) {
-      DD("no centroid, I stop");
+      DD("ZERO");
       found = true;
     }
 
-    if (i == slices_recon.size()-1) found = true;
+    if (currentSlice == slices_recon.size()-1) {
+      DD("end of slices");
+      found = true;
+    }
 
     // iterate
-    ++i;
+    ++currentSlice;
   } while (!found);
 
-
-  //MaskImageType::Pointer rr = clitk::JoinSlices<MaskImageType>(slices_recon, m_Mask, 2);
-  //clitk::writeImage<MaskImageType>(rr, "recon2.mhd");
-
   // End
   StopCurrentStep();
 }
 //--------------------------------------------------------------------
 
 
-
 #endif //#define CLITKEXTRACTMEDIASTINALVESSELSFILTER_TXX