Changeset 821 for soft/giet_vm/applications/rosenfeld/src-par/mca_rosenfeld.c

Timestamp:

May 6, 2016, 3:06:29 PM (9 years ago)

Author:

meunier

Message:

• Added several versions of rosenfeld: { SLOW, FAST } x { FEATURES, NO_FEATURES }
• Added native linux compilation support
• Added a script to check results natively
• Started to refactor nrc code

File:

• soft/giet_vm/applications/rosenfeld/src-par/mca_rosenfeld.c (modified) (28 diffs)

soft/giet_vm/applications/rosenfeld/src-par/mca_rosenfeld.c

@@ -12 +12 @@
 #include <string.h>
 #include <math.h>
-#include <user_barrier.h>
-#include <user_lock.h>
-
-#ifdef CLI
+#include <assert.h>
+#if PARMERGE
+#include <pthread.h>
+#endif
+
 #include "nrc_os_config.h"
+#include "config.h"
 #include "nrc.h"
+
+#if TARGET_OS == GIETVM
+    #include <user_barrier.h>
+    #include <user_lock.h>
+    #include <giet_config.h>
+#else
+    #include <stdbool.h>
 #endif
+
 
 #include "util.h"
 #include "ecc_common.h"
-
 #include "palette.h"
 #include "bmpNR.h"
-
+#include "clock.h"
 #include "str_ext.h"
+#include "ecc_features.h"
 
 // -----------
@@ -34 +44 @@
 #include "mca.h"
 
-extern giet_barrier_t main_barrier;
-
-// ----------------------------------
-uint32 FindRoot(uint32 * T, uint32 e)
-// ----------------------------------
+extern pthread_barrier_t main_barrier;
+extern int display_features;
+
+CLOCK_DEC;
+
+
+// -----------------------------------------
+static uint32 FindRoot(uint32 * T, uint32 e)
+// -----------------------------------------
 {
     uint32 r;
     
+    assert(e != 0);
     r = e;
     while (T[r] < r) {
         r = T[r];
     }
+    if (r == 0) {
+        printf("e = %d\n",e);
+        assert(0);
+    }
     return r;
 }
 
 
-// ---------------------------------------------------
-uint32 FindRoot_Dist(uint32 ** D, uint32 r, int shift)
-// ---------------------------------------------------
+// ----------------------------------------------------------
+static uint32 FindRoot_Dist(uint32 ** D, uint32 r, int shift)
+// ----------------------------------------------------------
 {
     uint32 e;
     uint32 e1;
     uint32 e0;
+
+    assert(r != 0);
     
     int mask = (1 << shift) - 1;
     
-    MCA_VERBOSE2(printf("FindRoot_Dist(%d) (alpha = %d) \n", r, shift));
+    MCA_VERBOSE2(printf("%s(%d, %d) \n", __func__, r, shift));
     do {
         e  = r;
@@ -66 +87 @@
         e0 = r & mask;
         r = D[e1][e0];
-        MCA_VERBOSE2(printf("FindRoot: D(%d) = D[%d,%d] = %d (alpha = %d)\n", e, e1, e0, r, shift));
+        MCA_VERBOSE2(printf("%s: D(%d) = D[%d,%d] = %d (alpha = %d)\n", __func__, e, e1, e0, r, shift));
     } while (r < e);
-    MCA_VERBOSE2(printf("FindRoot_Dist = %d \n\n", r));
+    MCA_VERBOSE2(printf("%s = %d \n\n", __func__, r));
+    assert(r != 0);
     return r;
 }
 
 
-// -----------------------------------------
-void SetRoot(uint32 * T, uint32 e, uint32 r)
-// -----------------------------------------
-{
-    while (T[e] < e) {
-        e = T[e];
-    }
-    T[e] = r;
-}
-
-
-// ----------------------------------------------------------
-void SetRoot_Dist(uint32 ** D, uint32 e, uint32 r, int shift)
-// ----------------------------------------------------------
+#if !FEATURES
+// --------------------------------------------------------------------------------
+static void SetRoot_Rosenfeld_Dist(uint32 ** D, uint32 root, uint32 eps, int shift)
+// --------------------------------------------------------------------------------
 {
     int mask = (1 << shift) - 1;
-    
-    uint32 e1 = e >> shift;
-    uint32 e0 = e & mask;
-    
-    D[e1][e0] = r;
-}
-
-
-// --------------------------------------------
-uint32 Union0(uint32 * T, uint32 ei, uint32 ej)
-// --------------------------------------------
-{
-    // version de la publication
-    // @QM : faut-il tester le cas == 0 ici aussi ?
-    uint32 ri, rj;
-    ri = (ei == 0) ? 0 : FindRoot(T, ei);
-    if (ei != ej) {
-        rj = (ej == 0) ? 0 : FindRoot(T, ej);
-        if (ri > rj) {
-            ri = rj;
-        }
-        SetRoot(T, ej, ri);
-    }
-    SetRoot(T, ei, ri);
-    return ri;
-}
-
-
-// -------------------------------------------------
-uint32 QuickUnion2(uint32 * T, uint32 e1, uint32 e2)
-// -------------------------------------------------
+    assert(root != 0 && eps != 0);
+    
+    uint32 r1 = root >> shift;
+    uint32 r0 = root & mask;
+    
+    D[r1][r0] = eps;
+}
+#endif // !FEATURES
+
+
+#if FEATURES && !PARMERGE
+// ----------------------------------------------------------------------------------------------------
+void SetRoot_Features_Rosenfeld_Dist(uint32 ** D, uint32 root, uint32 eps, int shift, RegionStats ** F)
+// ----------------------------------------------------------------------------------------------------
+{
+    assert(root != 0 && eps != 0);
+
+    MCA_VERBOSE2(printf("F(%d) += F(%d)\n", eps, root));
+    
+    int mask = (1 << shift) - 1;
+
+    // SetRoot_Rosenfeld_Dist
+    uint32 r1 = root >> shift;
+    uint32 r0 = root & mask;
+    
+    D[r1][r0] = eps;
+    
+    uint32 e1 = eps >> shift;
+    uint32 e0 = eps & mask;
+    
+    // version Dist de "RegionStats_Accumulate_Stats1_From_Index"
+    
+    // F(eps) = F(eps) U F(root)
+    
+    F[e1][e0].xmin = ui16min2(F[e1][e0].xmin, F[r1][r0].xmin);
+    F[e1][e0].xmax = ui16max2(F[e1][e0].xmax, F[r1][r0].xmax);
+    F[e1][e0].ymin = ui16min2(F[e1][e0].ymin, F[r1][r0].ymin);
+    F[e1][e0].ymax = ui16max2(F[e1][e0].ymax, F[r1][r0].ymax);
+    
+    F[e1][e0].S  += F[r1][r0].S;
+    F[e1][e0].Sx += F[r1][r0].Sx;
+    F[e1][e0].Sy += F[r1][r0].Sy;
+}
+#endif // FEATURES && !PARMERGE
+
+
+#if FEATURES && PARMERGE
+// -------------------------------------------------------------------------------------------------------------
+bool SetRoot_Parallel_Features_Rosenfeld_Dist(uint32 ** D, uint32 root, uint32 eps, int shift, RegionStats ** F)
+// -------------------------------------------------------------------------------------------------------------
+{
+    assert(root != 0 && eps != 0);
+
+    MCA_VERBOSE2(printf("F(%d) += F(%d)\n", eps, root));
+    
+    int mask = (1 << shift) - 1;
+
+    // SetRoot_Rosenfeld_Dist
+    uint32 r1 = root >> shift;
+    uint32 r0 = root & mask;
+    
+    uint32 e1 = eps >> shift;
+    uint32 e0 = eps & mask;
+
+    // Locking towards the root (first root, then eps)
+    pthread_spin_lock(&F[r1][r0].lock);
+    pthread_spin_lock(&F[e1][e0].lock);
+    // FIXME: merge these conditions later, when they both appear
+    if (D[e1][e0] != eps) {
+        // Someone change the root of epsilon, need to find the new root
+        printf("race cond 1\n");
+        pthread_spin_unlock(&F[e1][e0].lock);
+        pthread_spin_unlock(&F[r1][r0].lock);
+        return false;
+    }
+    if (D[r1][r0] != root) {
+        // Someone change the root of epsilon, need to find the new root
+        printf("race cond 2\n");
+        pthread_spin_unlock(&F[e1][e0].lock);
+        pthread_spin_unlock(&F[r1][r0].lock);
+        return false;
+    }
+
+    D[r1][r0] = eps;
+    
+    // F(eps) = F(eps) U F(root)
+    F[e1][e0].xmin = ui16min2(F[e1][e0].xmin, F[r1][r0].xmin);
+    F[e1][e0].xmax = ui16max2(F[e1][e0].xmax, F[r1][r0].xmax);
+    F[e1][e0].ymin = ui16min2(F[e1][e0].ymin, F[r1][r0].ymin);
+    F[e1][e0].ymax = ui16max2(F[e1][e0].ymax, F[r1][r0].ymax);
+    
+    F[e1][e0].S  += F[r1][r0].S;
+    F[e1][e0].Sx += F[r1][r0].Sx;
+    F[e1][e0].Sy += F[r1][r0].Sy;
+
+    pthread_spin_unlock(&F[e1][e0].lock);
+    pthread_spin_unlock(&F[r1][r0].lock);
+    return true;
+}
+#endif // FEATURES && PARMERGE
+
+
+
+#if FAST
+// --------------------------------------------------------
+static uint32 QuickUnion2(uint32 * T, uint32 e1, uint32 e2)
+// --------------------------------------------------------
 {
     // version QU de Union2
-    uint32 r1, r2, r;
-    
-    r1 = (e1 == 0) ? 0 : FindRoot(T, e1);
-    r2 = (e2 == 0) ? 0 : FindRoot(T, e2);
-    
-    r = ui32Min2(r1, r2);
-    
-    if (r1 != r) {
-        T[r1] = r; // SetRoot
-    }
-    if (r2 != r) {
-        T[r2] = r; // SetRoot
-    }
-    
-    return r;
-}
-
-
-// --------------------------------------------
-uint32 use1_QU_Rosenfeld(uint32 e1, uint32 * T)
-// --------------------------------------------
-{
-    return T[e1];
-}
-
-
-// -------------------------------------------------------
-uint32 use2_QU_Rosenfeld(uint32 e1, uint32 e2, uint32 * T)
-// -------------------------------------------------------
+    uint32 r1 = FindRoot(T, e1);
+    uint32 r2 = FindRoot(T, e2);
+    
+    assert(e1 != 0 && e2 != 0 && r1 != 0 && r2 != 0);
+    uint32 eps = ui32Min2(r1, r2);
+
+    if (r1 > eps) {
+        T[r1] = eps; // SetRoot sans besoin de remonter
+    }
+    if (r2 > eps) {
+        T[r2] = eps; // SetRoot sans besoin de remonter
+    }
+    assert(e1 != 0 && e2 != 0 && r1 != 0 && r2 != 0);
+    
+    return eps;
+}
+#endif // FAST
+
+
+#if FAST
+// ---------------------------------------------------
+static uint32 use1_QU_Rosenfeld(uint32 e1, uint32 * T)
+// ---------------------------------------------------
+{
+    return FindRoot(T, e1);
+}
+#endif // FAST
+
+
+#if FAST
+// --------------------------------------------------------------
+static uint32 use2_QU_Rosenfeld(uint32 e1, uint32 e2, uint32 * T)
+// --------------------------------------------------------------
 {
     return QuickUnion2(T, e1, e2);
 }
-
-
-// ----------------------------------------------------------------
-uint32 updateTable_Rosenfeld(uint32 * T, uint32 e, uint32 epsilon)
-// ----------------------------------------------------------------
-{
-    // notations e == v, epsilon == u avec v > u (v forcement different de u)
-    return Union0(T, e, epsilon); // original
-}
-
-
-// ----------------------------------------------------------------
-void vuse2_Rosenfeld(uint32 e1, uint32 e2, uint32 * T, uint32 ** D)
-// ----------------------------------------------------------------
-{
-    uint32 e;
-    uint32 a1;
-    uint32 a2;
-    
-    a1 = (e1 == 0) ? 0 : FindRoot(T, e1);
-    a2 = (e2 == 0) ? 0 : FindRoot(T, e2);
-    
-    if (a1 == a2) {
+#endif // FAST
+
+
+#if FAST && !FEATURES
+// ---------------------------------------------------------------------------------------
+static void vuse2_Rosenfeld_Dist(uint32 ed, uint32 el, uint32 * T, uint32 ** D, int alpha)
+// ---------------------------------------------------------------------------------------
+{
+    uint32 rd = FindRoot_Dist(D, ed, alpha);
+    
+    uint32 rl = T[el]; // car le premier acces est local
+    rl = FindRoot_Dist(D, rl, alpha);
+    
+    assert(ed != 0 && el != 0 && rd != 0 && rl != 0);
+    if (rd == rl) {
         return; // evite la backdoor
     }
     
-    // forcement positifs car appel depuis optimizedBorder qui a fait un test
-    if (a1 < a2) {
-        e = a1;
-        updateTable_Rosenfeld(T, a2, e);
+    // forcement positifs car appel depuis optimizedBorder
+    // qui a fait un test
+    if (rd < rl) {
+        SetRoot_Rosenfeld_Dist(D, rl, rd, alpha);
     }
     else {
-        e = a2;
-        updateTable_Rosenfeld(T, a1, e);
-    }
-}
-
-
-// ---------------------------------------------------------------------------
-void vuse3_Rosenfeld(uint32 e1, uint32 e2, uint32 e3, uint32 * T, uint32 ** D)
-// ---------------------------------------------------------------------------
-{
-    uint32 e;
-    uint32 a1;
-    uint32 a2;
-    uint32 a3;
-    
-    a1 = (e1 == 0) ? 0 : FindRoot(T, e1);
-    a2 = (e2 == 0) ? 0 : FindRoot(T, e2);
-    a3 = (e3 == 0) ? 0 : FindRoot(T, e3);
-    
-    if (a1 == a2 && a2 == a3) {
+        SetRoot_Rosenfeld_Dist(D, rd, rl, alpha);
+    }
+}
+#endif // FAST && !FEATURES
+
+
+#if FAST && !FEATURES
+// -----------------------------------------------------------------------------------------------------
+static void vuse3_Rosenfeld_Dist(uint32 ed1, uint32 ed2, uint32 el3, uint32 * T, uint32 ** D, int alpha)
+// -----------------------------------------------------------------------------------------------------
+{
+    uint32 r1 = FindRoot_Dist(D, ed1, alpha);
+    uint32 r2 = FindRoot_Dist(D, ed2, alpha);
+    
+    // QM
+    //uint32 r3 = FindRoot(T, el3); // local - distant
+    uint32 r3 = T[el3]; // local - distant
+    r3 = FindRoot_Dist(D, r3, alpha);
+
+    assert(ed1 != 0 && ed2 != 0 && el3 != 0 && r1 != 0 && r2 != 0 && r3 != 0);
+    
+    if (r1 == r2 && r2 == r3) {
         return;
     }
     
-    e = ui32Min3(a1, a2, a3);  // forcement positifs car appel depuis optimizedBorder qui a fait un test
-    
-    if (a1 > e) {
-        updateTable_Rosenfeld(T, a1, e);
-    }
-    a2 = T[a2];
-    if (a2 > e) {
-        updateTable_Rosenfeld(T, a2, e);
-    }
-    a3 = T[a3];
-    if (a3 > e) {
-        updateTable_Rosenfeld(T, a3, e);
-    }
-}
-
-
-// ----------------------------------------------
-uint32 solveTable_Rosenfeld(uint32 * T, uint32 ne)
-// ----------------------------------------------
-{
-    // equivalent a Flatten
-    // fermeture transitive sans pack
-    // (presence de trous dans les numeros d'etiquettes)
-    
-    uint32 e;
-    
-    for (e = 1; e <= ne; e++) {   
-        T[e] = T[T[e]];
-    }
-    return ne;
-}
-
-
-// ----------------------------------------------------------------------------------
-uint32 optimizedAccess_DT_Rosenfeld(uint32 ** E, int i, int j, uint32 * T, uint32 ne)
-// ----------------------------------------------------------------------------------
-{
-    // Decision Tree 8-connexe avec Quick-Union
-    uint32 a, b, c, d, e;
-    
-    b = E[i - 1][j];
-    if (b) {
-        e = use1_QU_Rosenfeld(b, T);
+    uint32 eps = ui32Min3(r1, r2, r3);  // forcement positifs car appel depuis optimizedBorder qui a fait un test
+    
+    if (r1 > eps) {
+        SetRoot_Rosenfeld_Dist(D, r1, eps, alpha);
+    }
+    //r2 = T[r2]; // @QM est-ce indispensable s'il n'y a pas de features ? (cf. slow no features)
+    // comment est-on sur que r2 (ou r3) est local ???
+    if (r2 > eps) {
+        SetRoot_Rosenfeld_Dist(D, r2, eps, alpha);
+    }
+    //r3 = T[r3];
+    if (r3 > eps) {
+        SetRoot_Rosenfeld_Dist(D, r3, eps, alpha);
+    }
+}
+#endif // FAST && !FEATURES
+
+
+#if FAST && FEATURES && !PARMERGE
+// -----------------------------------------------------------------------------------------------------------
+void vuse2_Features_Rosenfeld_Dist(uint32 ed, uint32 el, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
+// -----------------------------------------------------------------------------------------------------------
+{
+    assert(ed != 0 && el != 0);
+
+    uint32 rd = FindRoot_Dist(D, ed, alpha);
+    
+    uint32 rl = T[el]; // car le premier acces est local
+    assert(rl != 0);
+    rl = FindRoot_Dist(D, rl, alpha);
+    
+    assert(rd != 0 && rl != 0);
+
+    if (rd == rl) {
+        return; // evite la backdoor
+    }
+    
+    // forcement positifs car appel depuis optimizedBorder
+    // qui a fait un test
+    if (rd < rl) {
+        SetRoot_Features_Rosenfeld_Dist(D, rl, rd, alpha, F);
     }
     else {
-        c = E[i - 1][j + 1];
-        if (c) {
-            a = E[i - 1][j - 1];
+        SetRoot_Features_Rosenfeld_Dist(D, rd, rl, alpha, F);
+    }
+}
+#endif // FAST && FEATURES && !PARMERGE
+
+
+#if FAST && FEATURES && !PARMERGE
+// -------------------------------------------------------------------------------------------------------------------------
+void vuse3_Features_Rosenfeld_Dist(uint32 ed1, uint32 ed2, uint32 el3, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
+// -------------------------------------------------------------------------------------------------------------------------
+{
+    assert(ed1 != 0 && ed2 != 0 && el3 != 0);
+
+    uint32 r1 = FindRoot_Dist(D, ed1, alpha);
+    uint32 r2 = FindRoot_Dist(D, ed2, alpha);
+    
+    //uint32 r3 = FindRoot(T, el3); // local - distant
+    uint32 r3 = T[el3]; // local - distant
+    assert(r3 != 0);
+    r3 = FindRoot_Dist(D, r3, alpha);
+    
+    assert(r1 != 0 && r2 != 0 && r3 != 0);
+
+    if (r1 == r2 && r2 == r3) {
+        return;
+    }
+    
+    uint32 eps = ui32Min3(r1, r2, r3);  // forcement positifs car appel depuis optimizedBorder qui a fait un test
+    
+    if (r1 > eps) {
+        SetRoot_Features_Rosenfeld_Dist(D, r1, eps, alpha, F);
+    }
+    //r2 = T[r2];
+    if (r2 > eps && r2 != r1) {
+        SetRoot_Features_Rosenfeld_Dist(D, r2, eps, alpha, F);
+    }
+    //r3 = T[r3];
+    if (r3 > eps && r3 != r2 && r3 != r1) {
+        SetRoot_Features_Rosenfeld_Dist(D, r3, eps, alpha, F);
+    }
+}
+#endif // FAST && FEATURES && !PARMERGE
+
+
+#if FAST && FEATURES && PARMERGE
+// --------------------------------------------------------------------------------------------------------------------
+void vuse2_Parallel_Features_Rosenfeld_Dist(uint32 ed, uint32 el, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
+// --------------------------------------------------------------------------------------------------------------------
+{
+    bool ok;
+    assert(ed != 0 && el != 0);
+    uint32 rl = T[el]; // car le premier acces est local
+    assert(rl != 0);
+
+    uint32 rd;
+    
+    do {
+        rd = FindRoot_Dist(D, ed, alpha); // no lock
+        rl = FindRoot_Dist(D, rl, alpha);
+
+        assert(rd != 0 && rl != 0);
+
+        if (rd == rl) {
+            return; // evite la backdoor
+        }
+
+        // forcement positifs car appel depuis optimizedBorder
+        // qui a fait un test
+        if (rd < rl) {
+            ok = SetRoot_Parallel_Features_Rosenfeld_Dist(D, rl, rd, alpha, F);
+        }
+        else {
+            ok = SetRoot_Parallel_Features_Rosenfeld_Dist(D, rd, rl, alpha, F);
+        }
+    } while (!ok);
+}
+#endif // FAST && FEATURES && PARMERGE
+
+
+#if FAST && FEATURES && PARMERGE
+// ----------------------------------------------------------------------------------------------------------------------------------
+void vuse3_Parallel_Features_Rosenfeld_Dist(uint32 ed1, uint32 ed2, uint32 el3, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
+// ----------------------------------------------------------------------------------------------------------------------------------
+{
+    bool ok1, ok2, ok3;
+    assert(ed1 != 0 && ed2 != 0 && el3 != 0);
+
+    uint32 r1;
+    uint32 r2;
+    uint32 r3 = T[el3]; // local - distant
+    assert(r3 != 0);
+
+    do {
+        r1 = FindRoot_Dist(D, ed1, alpha);
+        r2 = FindRoot_Dist(D, ed2, alpha);
+        r3 = FindRoot_Dist(D, r3, alpha);
+    
+        assert(r1 != 0 && r2 != 0 && r3 != 0);
+
+        if (r1 == r2 && r2 == r3) {
+            return;
+        }
+    
+        uint32 eps = ui32Min3(r1, r2, r3);  // forcement positifs car appel depuis optimizedBorder qui a fait un test
+    
+        ok1 = true;
+        ok2 = true;
+        ok3 = true;
+        if (r1 > eps) {
+            ok1 = SetRoot_Parallel_Features_Rosenfeld_Dist(D, r1, eps, alpha, F);
+        }
+        if (r2 > eps && r2 != r1) {
+            ok2 = SetRoot_Parallel_Features_Rosenfeld_Dist(D, r2, eps, alpha, F);
+        }
+        if (r3 > eps && r3 != r2 && r3 != r1) {
+            ok3 = SetRoot_Parallel_Features_Rosenfeld_Dist(D, r3, eps, alpha, F);
+        }
+    } while (!(ok1 && ok2 && ok3));
+}
+#endif // FAST && FEATURES && PARMERGE
+
+
+
+
+#if FAST && !FEATURES
+// ------------------------------------------------------------------------------------------------------
+static void optimizedBorder_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha)
+// ------------------------------------------------------------------------------------------------------
+{
+    uint32 a, b, c, x;
+    
+    x = E[i][j];
+    if (x) {
+        b = E[i - 1][j];
+        if (b) {
+            vuse2_Rosenfeld_Dist(b, x, T, D, alpha); // dist, local
+        }
+        else {
+            c = E[i - 1][j + 1];
+            if (c) {
+                a = E[i - 1][j - 1];
+                if (a) {
+                    vuse3_Rosenfeld_Dist(a, c, x, T, D, alpha); // dist, local
+                }
+                else {
+                    vuse2_Rosenfeld_Dist(c, x, T, D, alpha); // dist, local
+                }
+            }
+            else {
+                a = E[i - 1][j - 1];
+                if (a) {
+                    vuse2_Rosenfeld_Dist(a, x, T, D, alpha); // dist, local
+                }
+            }
+        }
+    }
+}
+#endif // FAST && !FEATURES
+
+
+#if FAST && !FEATURES
+// ---------------------------------------------------------------------------------------------------
+static void optimizedBorderLeft_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha)
+// ---------------------------------------------------------------------------------------------------
+{
+    uint32 x = E[i][j];
+    if (x) {
+        uint32 b = E[i - 1][j];
+        if (b) {
+            vuse2_Rosenfeld_Dist(b, x, T, D, alpha); // dist, local
+        }
+        else {
+            uint32 c = E[i - 1][j + 1];
+            if (c) {
+                vuse2_Rosenfeld_Dist(c, x, T, D, alpha); // dist, local
+            }
+        }
+    }
+}
+#endif // FAST && !FEATURES
+
+
+#if FAST && !FEATURES
+// -----------------------------------------------------------------------------------------------------------
+static void optimizedBorderRight_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha)
+// -----------------------------------------------------------------------------------------------------------
+{
+    // copie de optimizedBorder_Rosenfeld
+    // test d'existance de ex en local local
+
+    uint32 b = E[i - 1][j];
+    uint32 x = E[i][j];
+    
+    if (x) {
+        if (b) {
+            vuse2_Rosenfeld_Dist(b, x, T, D, alpha); // dist, local
+        }
+        else {
+            uint32 a = E[i - 1][j - 1];
             if (a) {
-                e = use2_QU_Rosenfeld(a, c, T);
-            }
-            else {
-                d = E[i][j - 1];
-                if (d) {
-                    e = use2_QU_Rosenfeld(c, d, T);
-                }
-                else {
-                    e = use1_QU_Rosenfeld(c, T);
-                }
-            }
-        }
-        else {
-            a = E[i - 1][j - 1];
-            if (a) {
-                e = use1_QU_Rosenfeld(a, T);
-            }
-            else {
-                d = E[i][j - 1];
-                if (d) {
-                    e = use1_QU_Rosenfeld(d, T);
-                }
-                else {
-                    e = ++ne;
-                }
-            }
-        }
-    }
-    E[i][j] = e;
-    return ne;
-}
-
-
-// ------------------------------------------------------------------------------------------
-void optimizedBorder_Rosenfeld(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha)
-// ------------------------------------------------------------------------------------------
-{
-    // copie de optimizedBorder_Rosenfeld
-    uint32 a, b, c, x;
-    
-    b = E[i - 1][j];
-    x = E[i][j];
-    
-    if (b) {
-        //printf("%d = %d\n", b, x);
-        vuse2_Rosenfeld(b, x, T, D);
-    }
-    else {
-        c = E[i - 1][j + 1];
-        if (c) {
-            a = E[i - 1][j - 1];
-            if (a) {
-                //printf("%d = %d = %d\n", a, c, x);
-                vuse3_Rosenfeld(a, c, x, T, D);
-            }
-            else {
-                //printf("%d = %d\n", c, x);
-                vuse2_Rosenfeld(c, x, T, D);
-            }
-        }
-        else {
-            a = E[i - 1][j - 1];
-            if (a) {
-                //printf("%d = %d\n", a, x);
-                vuse2_Rosenfeld(a, x, T, D);
-            }
-        }
-    }
-}
-
-
-// -----------------------------------------------------------------------------------------------------------------
-void borderMerging_Fast_Rosenfeld_Dist(uint8 **X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha)
-// -----------------------------------------------------------------------------------------------------------------
-{
-    for (int j = 0; j < width; j++) {
-        if (X[i][j])  {
-            optimizedBorder_Rosenfeld(E, i, j, T, D, alpha);
-        }
-    }
-    return;
-}
-
-
-// ------------------------------------------------------------------------------------------------------------------
-void borderMerging_Slow_Rosenfeld_Dist(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha)
-// ------------------------------------------------------------------------------------------------------------------
-{
-    // copie de borderMerging_Rosenfeld_UF_Fast2_8C
-    
+                vuse2_Rosenfeld_Dist(a, x, T, D, alpha); // dist, local
+            }
+        }
+    }
+}
+#endif // FAST && !FEATURES
+
+
+#if FAST && !FEATURES
+// ------------------------------------------------------------------------------------------------------------------------
+static void borderMerging_Fast_Rosenfeld_Dist(uint8 **X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha)
+// ------------------------------------------------------------------------------------------------------------------------
+{
+    // Prologue
+    optimizedBorderLeft_Rosenfeld_Dist(E, i, 0, T, D, alpha);
+    // Boucle principale
+    for (int j = 1; j < width - 1; j++) {
+        optimizedBorder_Rosenfeld_Dist(E, i, j, T, D, alpha);
+    }
+    // Epilogue
+    optimizedBorderRight_Rosenfeld_Dist(E, i, width - 1, T, D, alpha);
+}
+#endif // FAST && !FEATURES
+
+
+#if SLOW && !FEATURES
+// -------------------------------------------------------------------------------------------------------------------------
+static void borderMerging_Slow_Rosenfeld_Dist(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha)
+// -------------------------------------------------------------------------------------------------------------------------
+{
     int j;
     
-    uint32 e;
-    
+    uint32 eps;
     uint32 e1, e2, e3, ex;
     uint32 r1, r2, r3, rx;
@@ -358 +577 @@
     // -- prologue --
     // --------------
-    MCA_VERBOSE2(printf("[borderMerging_Slow_Rosenfeld_Dist] i = %d\n", i));
+    MCA_VERBOSE2(printf("[%s] i = %d\n", __func__, i));
     
     j = 0;
@@ -365 +584 @@
     if (ex) {
         
-        MCA_VERBOSE2(printf("[borderMerging_Slow_Rosenfeld_Dist] j = %d\n", j));
+        MCA_VERBOSE2(printf("[%s] j = %d\n", __func__, j));
         
         e2 = E[i - 1][j];
         e3 = E[i - 1][j + 1];
-        
-        r2 = FindRoot_Dist(D, e2, alpha);
-        r3 = FindRoot_Dist(D, e3, alpha);
-        rx = FindRoot(T, ex); // we already tested that ex != 0
-        
+
+        // test pour eviter acces distant
+        r2 = e2 ? FindRoot_Dist(D, e2, alpha) : 0;
+        r3 = e3 ? FindRoot_Dist(D, e3, alpha) : 0;
+
+        rx = T[ex];
+        rx = FindRoot_Dist(D, rx, alpha);
+ 
         MCA_VERBOSE2(printf("\n"));
         MCA_VERBOSE2(printf("e2 = %4d -> %4d\n", e2, r2));
@@ -379 +601 @@
         MCA_VERBOSE2(printf("ex = %4d -> %4d\n", ex, rx));
         
-        e = ui32MinNonNul3(r2, r3, rx);
+        eps = ui32MinNonNul3(r2, r3, rx);
         
         // Quick-Union
-        if (r2 > e) {
-            SetRoot_Dist(D, r2, e, alpha);
-            MCA_VERBOSE2(printf("D[%4d] <- %d\n", r2, e));
-        }
-        if (r3 > e) {
-            SetRoot_Dist(D, r3, e, alpha);
-            MCA_VERBOSE2(printf("D[%4d] <- %d\n", r3, e));
-        }
-        if (rx > e) {
-            SetRoot(T, rx, e);
-            MCA_VERBOSE2(printf("D[%4d] <- %d\n", rx, e));
+        if (r2 > eps) {
+            SetRoot_Rosenfeld_Dist(D, r2, eps, alpha);
+            MCA_VERBOSE2(printf("D[%4d] <- %d\n", r2, eps));
+        }
+        if (r3 > eps) {
+            SetRoot_Rosenfeld_Dist(D, r3, eps, alpha);
+            MCA_VERBOSE2(printf("D[%4d] <- %d\n", r3, eps));
+        }
+        if (rx > eps) {
+            SetRoot_Rosenfeld_Dist(D, rx, eps, alpha);
+            MCA_VERBOSE2(printf("D[%4d] <- %d\n", rx, eps));
         }
         MCA_VERBOSE2(printf("\n"));
-        // attention SetRoot fait un while inutile
     }
     
@@ -408 +629 @@
         // que le cas general (pour faire un code simple)
         if (ex) {
-            MCA_VERBOSE2(printf("[borderMerging_Slow_Rosenfeld_Dist] j = %d\n", j));
+            MCA_VERBOSE2(printf("[%s] j = %d\n", __func__, j));
             
             e1 = E[i - 1][j - 1];
@@ -414 +635 @@
             e3 = E[i - 1][j + 1];
         
-            r1 = FindRoot_Dist(D, e1, alpha);
-            r2 = FindRoot_Dist(D, e2, alpha);
-            r3 = FindRoot_Dist(D, e3, alpha);
-            rx = FindRoot(T, ex); // we already tested that ex != 0
-        
+            // test pour eviter acces distant
+            r1 = e1 ? FindRoot_Dist(D, e1, alpha) : 0;
+            r2 = e2 ? FindRoot_Dist(D, e2, alpha) : 0;
+            r3 = e3 ? FindRoot_Dist(D, e3, alpha) : 0;
+
+            rx = T[ex];
+            rx = FindRoot_Dist(D, rx, alpha);
+
             MCA_VERBOSE2(printf("\n"));
             MCA_VERBOSE2(printf("e1 = %4d -> %4d\n", e1, r1));
@@ -425 +649 @@
             MCA_VERBOSE2(printf("ex = %4d -> %4d\n", ex, rx));
             
-            e = ui32MinNonNul4(r1, r2, r3, rx);
+            eps = ui32MinNonNul4(r1, r2, r3, rx);
             
             // Quick-Union
-            if (r1 > e) {
-                SetRoot_Dist(D, r1, e, alpha);
-                MCA_VERBOSE2(printf("D[%4d] <- %d\n", r1, e));
-            }
-            if (r2 > e) {
-                SetRoot_Dist(D, r2, e, alpha);
-                MCA_VERBOSE2(printf("D[%4d] <- %d\n", r2, e));
-            }
-            if (r3 > e) {
-                SetRoot_Dist(D, r3, e, alpha);
-                MCA_VERBOSE2(printf("D[%4d] <- %d\n", r3, e));
-            }
-            if (rx > e) {
-                // @QM pourquoi pas T[e] = rx; ?
-                //SetRoot(T, rx, e); 
-                T[e] = rx;
-                MCA_VERBOSE2(printf("D[%4d] <- %d\n", rx, e));
+            if (r1 > eps) {
+                SetRoot_Rosenfeld_Dist(D, r1, eps, alpha);
+                MCA_VERBOSE2(printf("D[%4d] <- %d\n", r1, eps));
+            }
+            if (r2 > eps) {
+                SetRoot_Rosenfeld_Dist(D, r2, eps, alpha);
+                MCA_VERBOSE2(printf("D[%4d] <- %d\n", r2, eps));
+            }
+            if (r3 > eps) {
+                SetRoot_Rosenfeld_Dist(D, r3, eps, alpha);
+                MCA_VERBOSE2(printf("D[%4d] <- %d\n", r3, eps));
+            }
+            if (rx > eps) {
+                SetRoot_Rosenfeld_Dist(D, rx, eps, alpha);
+                MCA_VERBOSE2(printf("D[%4d] <- %d\n", rx, eps));
             }
             MCA_VERBOSE2(printf("\n"));
@@ -460 +682 @@
     if (ex) {
         
-        MCA_VERBOSE2(printf("[borderMerging_Slow_Rosenfeld_Dist] j = %d\n", j));
+        MCA_VERBOSE2(printf("[%s] j = %d\n", __func__, j));
         
         e1 = E[i - 1][j - 1];
         e2 = E[i - 1][j];
-        
-        r1 = FindRoot_Dist(D, e1, alpha);
-        r2 = FindRoot_Dist(D, e2, alpha);
-        rx = FindRoot(T, ex); // we already tested that ex != 0
-        
+
+        // test pour eviter acces distant
+        r1 = e1 ? FindRoot_Dist(D, e1, alpha) : 0;
+        r2 = e2 ? FindRoot_Dist(D, e2, alpha) : 0;
+
+        rx = T[ex];
+        rx = FindRoot_Dist(D, rx, alpha);
+
         MCA_VERBOSE2(printf("\n"));
         MCA_VERBOSE2(printf("e1 = %4d -> %4d\n", e1, r1));
@@ -474 +699 @@
         MCA_VERBOSE2(printf("ex = %4d -> %4d\n", ex, rx));
         
-        e = ui32MinNonNul3(r1, r2, rx);
+        eps = ui32MinNonNul3(r1, r2, rx);
         
         // Quick-Union
-        if (r1 > e) {
-            SetRoot_Dist(D, r1, e, alpha);
-            MCA_VERBOSE2(printf("D[%4d] <- %d\n", r1, e));
-        }
-        if (r2 > e) {
-            SetRoot_Dist(D, r2, e, alpha);
-            MCA_VERBOSE2(printf("D[%4d] <- %d\n", r2, e));
-        }
-        if (rx > e) {
-            SetRoot(T, rx, e);
-            MCA_VERBOSE2(printf("D[%4d] <- %d\n", rx, e));
+        if (r1 > eps) {
+            SetRoot_Rosenfeld_Dist(D, r1, eps, alpha);
+            MCA_VERBOSE2(printf("D[%4d] <- %d\n", r1, eps));
+        }
+        if (r2 > eps) {
+            SetRoot_Rosenfeld_Dist(D, r2, eps, alpha);
+            MCA_VERBOSE2(printf("D[%4d] <- %d\n", r2, eps));
+        }
+        if (rx > eps) {
+            SetRoot_Rosenfeld_Dist(D, rx, eps, alpha);
+            MCA_VERBOSE2(printf("D[%4d] <- %d\n", rx, eps));
         }
         MCA_VERBOSE2(printf("\n"));
@@ -493 +718 @@
     return;
 }
-
-
-// -------------------------------------------------------------------------------------------------------------
-void borderMerging_Rosenfeld_Dist(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha)
-// -------------------------------------------------------------------------------------------------------------
-{
+#endif // SLOW && !FEATURES
+
+
+#if SLOW && FEATURES
+// ----------------------------------------------------------------------------------------------------------------------------------------------------
+static void borderMerging_Slow_Features_Rosenfeld_Dist(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
+// ----------------------------------------------------------------------------------------------------------------------------------------------------
+{
+    int j = 0;
+    
+    uint32 eps;
+    
+    uint32 e1, e2, e3, ex;
+    uint32 r1, r2, r3, rx;
+    
+    // --------------
+    // -- prologue --
+    // --------------
+    MCA_VERBOSE2(printf("[%s] i = %d\n", __func__, i));
+    
+    ex = E[i][j];
+    
+    if (ex) {
+        
+        MCA_VERBOSE2(printf("[%s] j = %d\n", __func__, j));
+        
+        e2 = E[i - 1][j];
+        e3 = E[i - 1][j + 1];
+        
+        if (e2 || e3) {
+        
+            // test pour eviter acces distant
+            r2 = e2 ? FindRoot_Dist(D, e2, alpha) : 0;
+            r3 = e3 ? FindRoot_Dist(D, e3, alpha) : 0;
+
+            rx = T[ex];
+            rx = FindRoot_Dist(D, rx, alpha);
+            
+            eps = ui32MinNonNul3(r2, r3, rx);
+            
+            MCA_VERBOSE2(printf("\n"));
+            MCA_VERBOSE2(printf("e2  = %5d -> r2 = %5d\n", e2, r2));
+            MCA_VERBOSE2(printf("e3  = %5d -> r3 = %5d\n", e3, r3));
+            MCA_VERBOSE2(printf("ex  = %5d -> rx = %5d\n", ex, rx));
+            MCA_VERBOSE2(printf("eps = %5d\n", eps));
+            
+            // Quick-Union
+            // @QM
+            if (r2 > eps) {
+                SetRoot_Features_Rosenfeld_Dist(D, r2, eps, alpha, F);
+                MCA_VERBOSE2(printf("D[%5d] <- %d\n", r2, eps));
+            }
+            if (r3 > 0) {
+                r3 = FindRoot_Dist(D, r3, alpha);
+            }
+            // Pour le cas où r2 == r3, il ne faut pas ajouter deux fois les features
+            //if (r3 > eps && r3 != r2) {
+            if (r3 > eps) {
+                SetRoot_Features_Rosenfeld_Dist(D, r3, eps, alpha, F);
+                MCA_VERBOSE2(printf("D[%5d] <- %d\n", r3, eps));
+            }
+            rx = FindRoot_Dist(D, rx, alpha);
+            //if (rx > eps && rx != r3 && rx != r2) {
+            if (rx > eps) {
+                SetRoot_Features_Rosenfeld_Dist(D, rx, eps, alpha, F);
+                MCA_VERBOSE2(printf("D[%5d] <- %d\n", rx, eps));
+            }
+            MCA_VERBOSE2(printf("---------------------------\n"));
+        }
+    }
+    
+    // -----------------------
+    // -- boucle principale --
+    // -----------------------
+    
+    for (j = 0 + 1; j < width - 1; j++) {
+        
+        ex = E[i][j];
+        
+        if (ex) {
+            
+            MCA_VERBOSE2(printf("[%s] j = %d\n", __func__, j));
+            
+            e1 = E[i - 1][j - 1];
+            e2 = E[i - 1][j];
+            e3 = E[i - 1][j + 1];
+            
+            if (e1 || e2 || e3) {
+                // test pour eviter un acces distant
+                r1 = e1 ? FindRoot_Dist(D, e1, alpha) : 0;
+                r2 = e2 ? FindRoot_Dist(D, e2, alpha) : 0;
+                r3 = e3 ? FindRoot_Dist(D, e3, alpha) : 0;
+
+                rx = T[ex];
+                rx = FindRoot_Dist(D, rx, alpha);
+                
+                eps = ui32MinNonNul4(r1, r2, r3, rx);
+
+                MCA_VERBOSE2(printf("\n"));
+                MCA_VERBOSE2(printf("e1  = %5d -> r1 = %5d\n", e1, r1));
+                MCA_VERBOSE2(printf("e2  = %5d -> r2 = %5d\n", e2, r2));
+                MCA_VERBOSE2(printf("e3  = %5d -> r3 = %5d\n", e3, r3));
+                MCA_VERBOSE2(printf("ex  = %5d -> rx = %5d\n", ex, rx));
+                MCA_VERBOSE2(printf("eps = %5d\n", eps));
+                
+                // Quick-Union
+                // @QM
+                if (r1 > eps) {
+                    SetRoot_Features_Rosenfeld_Dist(D, r1, eps, alpha, F);
+                    MCA_VERBOSE2(printf("D[%5d] <- %d\n", r1, eps));
+                }
+                if (r2 > 0) {
+                    r2 = FindRoot_Dist(D, r2, alpha);
+                }
+                //if (r2 > eps && r2 != r1) {
+                if (r2 > eps) {
+                    SetRoot_Features_Rosenfeld_Dist(D, r2, eps, alpha, F);
+                    MCA_VERBOSE2(printf("D[%5d] <- %d\n", r2, eps));
+                }
+                if (r3 > 0) {
+                    r3 = FindRoot_Dist(D, r3, alpha);
+                }
+                //if (r3 > eps && r3 != r2 && r3 != r1) {
+                if (r3 > eps) {
+                    SetRoot_Features_Rosenfeld_Dist(D, r3, eps, alpha, F);
+                    MCA_VERBOSE2(printf("D[%5d] <- %d\n", r3, eps));
+                }
+                rx = FindRoot_Dist(D, rx, alpha);
+                //if (rx > eps && rx != r3 && rx != r2 && rx != r1) {
+                if (rx > eps) {
+                    SetRoot_Features_Rosenfeld_Dist(D, rx, eps, alpha, F);
+                    MCA_VERBOSE2(printf("D[%5d] <- %d\n", rx, eps));
+                }
+                MCA_VERBOSE2(puts("---------------------------\n"));
+                
+                // attention SetRoot fait un while inutile
+            }
+        }
+    }
+    
+    // --------------
+    // -- epilogue --
+    // --------------
+    
+    j = width - 1;
+    ex = E[i][j];
+    
+    if (ex) {
+        
+        MCA_VERBOSE2(printf("[%s] j = %d\n", __func__, j));
+        
+        e1 = E[i - 1][j - 1];
+        e2 = E[i - 1][j];
+        
+        if (e1 || e2) {
+        
+            // test pour eviter acces distant
+            r1 = e1 ? FindRoot_Dist(D, e1, alpha) : 0;
+            r2 = e2 ? FindRoot_Dist(D, e2, alpha) : 0;
+
+            rx = T[ex];
+            rx = FindRoot_Dist(D, rx, alpha);
+            
+            eps = ui32MinNonNul3(r1, r2, rx);
+            
+            MCA_VERBOSE2(printf("\n"));
+            MCA_VERBOSE2(printf("e1  = %5d -> r1 = %5d\n", e1, r1));
+            MCA_VERBOSE2(printf("e2  = %5d -> r2 = %5d\n", e2, r2));
+            MCA_VERBOSE2(printf("ex  = %5d -> rx = %5d\n", ex, rx));
+            MCA_VERBOSE2(printf("eps = %5d\n", eps));
+            
+            // Quick-Union
+            if (r1 > eps) {
+                SetRoot_Features_Rosenfeld_Dist(D, r1, eps, alpha, F);
+                MCA_VERBOSE2(printf("D[%5d] <- %d\n", r1, eps));
+            }
+            if (r2 > 0) {
+                r2 = FindRoot_Dist(D, r2, alpha);
+            }
+            //if (r2 > eps && r2 != r1) {
+            if (r2 > eps) {
+                SetRoot_Features_Rosenfeld_Dist(D, r2, eps, alpha, F);
+                MCA_VERBOSE2(printf("D[%5d] <- %d\n", r2, eps));
+            }
+            rx = FindRoot_Dist(D, rx, alpha);
+            //if (rx > eps && rx != r2 && rx != r1) {
+            if (rx > eps) {
+                SetRoot_Features_Rosenfeld_Dist(D, rx, eps, alpha, F);
+                MCA_VERBOSE2(printf("D[%5d] <- %d\n", rx, eps));
+            }
+            MCA_VERBOSE2(printf("---------------------------\n"));
+        }
+    }
+    return;
+}
+#endif // SLOW && FEATURES
+
+
+#if FAST && FEATURES && !PARMERGE
+// --------------------------------------------------------------------------------------------------------------------------
+void optimizedBorder_Features_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
+// --------------------------------------------------------------------------------------------------------------------------
+{
+    // copie de optimizedBorder_Rosenfeld
+    uint32 a, b, c, x;
+    
+    x = E[i][j];
+    
+    if (x) {
+        b = E[i - 1][j];
+        if (b) {
+            vuse2_Features_Rosenfeld_Dist(b, x, T, D, alpha, F); // dist, local
+        }
+        else {
+            c = E[i - 1][j + 1];
+            if (c) {
+                a = E[i - 1][j - 1];
+                if (a) {
+                    vuse3_Features_Rosenfeld_Dist(a, c, x, T, D, alpha, F); // dist, local
+                }
+                else {
+                    vuse2_Features_Rosenfeld_Dist(c, x, T, D, alpha, F); // dist, local
+                }
+            }
+            else {
+                a = E[i - 1][j - 1];
+                if (a) {
+                    vuse2_Features_Rosenfeld_Dist(a, x, T, D, alpha, F); // dist, local
+                }
+            }
+        }
+    }
+}
+#endif // FAST && FEATURES && !PARMERGE
+
+
+#if FAST && FEATURES && !PARMERGE
+// ------------------------------------------------------------------------------------------------------------------------------
+void optimizedBorderLeft_Features_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
+// ------------------------------------------------------------------------------------------------------------------------------
+{
+    uint32 x = E[i][j];
+    
+    if (x) {
+        uint32 b = E[i - 1][j];
+        if (b) {
+            vuse2_Features_Rosenfeld_Dist(b, x, T, D, alpha, F); // dist, local
+        }
+        else {
+            uint32 c = E[i - 1][j + 1];
+            if (c) {
+                vuse2_Features_Rosenfeld_Dist(c, x, T, D, alpha, F); // dist, local
+            }
+        }
+    }
+}
+#endif // FAST && FEATURES && !PARMERGE
+
+
+#if FAST && FEATURES && !PARMERGE
+// -------------------------------------------------------------------------------------------------------------------------------
+void optimizedBorderRight_Features_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
+// -------------------------------------------------------------------------------------------------------------------------------
+{
+    // copie de optimizedBorder_Rosenfeld
+    // test d'existance de ex en local local
+    
+    uint32 x = E[i][j];
+    
+    if (x) {
+        uint32 b = E[i - 1][j];
+        if (b) {
+            vuse2_Features_Rosenfeld_Dist(b, x, T, D, alpha, F); // dist, local
+        }
+        else {
+            uint32 a = E[i - 1][j - 1];
+            if (a) {
+                vuse2_Features_Rosenfeld_Dist(a, x, T, D, alpha, F); // dist, local
+            }
+        }
+    }
+}
+#endif // FAST && FEATURES && !PARMERGE
+
+
+#if FAST && FEATURES && PARMERGE
+// -----------------------------------------------------------------------------------------------------------------------------------
+void optimizedBorder_Parallel_Features_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
+// -----------------------------------------------------------------------------------------------------------------------------------
+{
+    // copie de optimizedBorder_Rosenfeld
+    uint32 a, b, c, x;
+    
+    x = E[i][j];
+    
+    if (x) {
+        b = E[i - 1][j];
+        if (b) {
+            vuse2_Parallel_Features_Rosenfeld_Dist(b, x, T, D, alpha, F); // dist, local
+        }
+        else {
+            c = E[i - 1][j + 1];
+            if (c) {
+                a = E[i - 1][j - 1];
+                if (a) {
+                    vuse3_Parallel_Features_Rosenfeld_Dist(a, c, x, T, D, alpha, F); // dist, local
+                }
+                else {
+                    vuse2_Parallel_Features_Rosenfeld_Dist(c, x, T, D, alpha, F); // dist, local
+                }
+            }
+            else {
+                a = E[i - 1][j - 1];
+                if (a) {
+                    vuse2_Parallel_Features_Rosenfeld_Dist(a, x, T, D, alpha, F); // dist, local
+                }
+            }
+        }
+    }
+}
+#endif // FAST && FEATURES && PARMERGE
+
+
+#if FAST && FEATURES && PARMERGE
+// ---------------------------------------------------------------------------------------------------------------------------------------
+void optimizedBorderLeft_Parallel_Features_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
+// ---------------------------------------------------------------------------------------------------------------------------------------
+{
+    uint32 x = E[i][j];
+    
+    if (x) {
+        uint32 b = E[i - 1][j];
+        if (b) {
+            vuse2_Parallel_Features_Rosenfeld_Dist(b, x, T, D, alpha, F); // dist, local
+        }
+        else {
+            uint32 c = E[i - 1][j + 1];
+            if (c) {
+                vuse2_Parallel_Features_Rosenfeld_Dist(c, x, T, D, alpha, F); // dist, local
+            }
+        }
+    }
+}
+#endif // FAST && FEATURES && PARMERGE
+
+
+#if FAST && FEATURES && PARMERGE
+// ----------------------------------------------------------------------------------------------------------------------------------------
+void optimizedBorderRight_Parallel_Features_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
+// ----------------------------------------------------------------------------------------------------------------------------------------
+{
+    // copie de optimizedBorder_Rosenfeld
+    // test d'existance de ex en local local
+    
+    uint32 x = E[i][j];
+    
+    if (x) {
+        uint32 b = E[i - 1][j];
+        if (b) {
+            vuse2_Parallel_Features_Rosenfeld_Dist(b, x, T, D, alpha, F); // dist, local
+        }
+        else {
+            uint32 a = E[i - 1][j - 1];
+            if (a) {
+                vuse2_Parallel_Features_Rosenfeld_Dist(a, x, T, D, alpha, F); // dist, local
+            }
+        }
+    }
+}
+#endif // FAST && FEATURES && PARMERGE
+
+
+#if FAST && FEATURES
+// ---------------------------------------------------------------------------------------------------------------------------------------------
+void borderMerging_Fast_Features_Rosenfeld_Dist(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
+// ---------------------------------------------------------------------------------------------------------------------------------------------
+{
+    MCA_VERBOSE2(printf("[%s]", __func__));
+    
+#if PARMERGE
+    optimizedBorderLeft_Parallel_Features_Rosenfeld_Dist(E, i, 0, T, D, alpha, F);
+#else
+    optimizedBorderLeft_Features_Rosenfeld_Dist(E, i, 0, T, D, alpha, F);
+#endif
+    
+    for (int j = 1; j < width - 1; j++) {
+#if PARMERGE
+        optimizedBorder_Parallel_Features_Rosenfeld_Dist(E, i, j, T, D, alpha, F);
+#else
+        optimizedBorder_Features_Rosenfeld_Dist(E, i, j, T, D, alpha, F);
+#endif
+    }
+    
+#if PARMERGE
+    optimizedBorderRight_Parallel_Features_Rosenfeld_Dist(E, i, width - 1, T, D, alpha, F);
+#else
+    optimizedBorderRight_Features_Rosenfeld_Dist(E, i, width - 1, T, D, alpha, F);
+#endif
+}
+#endif // FAST && FEATURES
+
+
+#if !FEATURES
+// --------------------------------------------------------------------------------------------------------------------
+static void borderMerging_Rosenfeld_Dist(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha)
+// --------------------------------------------------------------------------------------------------------------------
+{
+#if SLOW
     borderMerging_Slow_Rosenfeld_Dist(X, i, width, E, T, D, alpha);
-}
-
-
-// ---------------------------------------------------------------------------------------------
-uint32 line0Labeling_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
-// ---------------------------------------------------------------------------------------------
+#elif FAST
+    borderMerging_Fast_Rosenfeld_Dist(X, i, width, E, T, D, alpha);
+#else
+#error "Please define SLOW or FAST for the Rosenfeld version"
+#endif
+}
+#endif // !FEATURES
+
+
+#if FEATURES
+// -----------------------------------------------------------------------------------------------------------------------------------------------
+static void borderMerging_Features_Rosenfeld_Dist(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
+// -----------------------------------------------------------------------------------------------------------------------------------------------
+{
+#if SLOW
+    borderMerging_Slow_Features_Rosenfeld_Dist(X, i, width, E, T, D, alpha, F);
+#elif FAST
+    borderMerging_Fast_Features_Rosenfeld_Dist(X, i, width, E, T, D, alpha, F);
+#else
+#error "Please define SLOW or FAST for the Rosenfeld version"
+#endif
+}
+#endif // FEATURES
+
+
+// ----------------------------------------------------------------------------------------------------
+static uint32 line0Labeling_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
+// ----------------------------------------------------------------------------------------------------
 {
     int j;
@@ -542 +1191 @@
 
 
-// -------------------------------------------------------------------------------------------------
-uint32 lineLabeling_Slow_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
-// -------------------------------------------------------------------------------------------------
+#if SLOW
+// --------------------------------------------------------------------------------------------------------
+static uint32 lineLabeling_Slow_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
+// --------------------------------------------------------------------------------------------------------
 {
     // version lineLabeling_Rosenfeld_UF_QU_8C avec Quick-Union
@@ -633 +1283 @@
                     
                     e = ui32MinNonNul4(r1, r2, r3, r4);
-                    giet_pthread_assert(e != 0, "e = 0\n");
                     
                     // Quick-Union
@@ -707 +1356 @@
     return ne;
 }
-
-
-// -------------------------------------------------------------------------------------------------
-uint32 lineLabeling_Fast_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
-// -------------------------------------------------------------------------------------------------
-{
+#endif // SLOW
+
+
+#if FAST
+// ---------------------------------------------------------------------------------------------
+static uint32 optimizedAccessLeft_DT_Rosenfeld(uint32 ** E, int i, int j, uint32 * T, uint32 ne)
+// ---------------------------------------------------------------------------------------------
+{
+    // Decision Tree 8-connexe avec Quick-Union
+    uint32 b, c, e;
+    
+    b = E[i - 1][j];
+    if (b) {
+        e = use1_QU_Rosenfeld(b, T);
+    }
+    else {
+        c = E[i - 1][j + 1];
+        if (c) {
+            e = use1_QU_Rosenfeld(c, T);
+        }
+        else {
+            e = ++ne;
+        }
+    }
+    E[i][j] = e;
+    return ne;
+}
+#endif // FAST
+
+
+#if FAST
+// ----------------------------------------------------------------------------------------------
+static uint32 optimizedAccessRight_DT_Rosenfeld(uint32 ** E, int i, int j, uint32 * T, uint32 ne)
+// ----------------------------------------------------------------------------------------------
+{
+    // Decision Tree 8-connexe avec Quick-Union
+    uint32 a, b, d, e;
+    
+    b = E[i - 1][j];
+    if (b) {
+        e = use1_QU_Rosenfeld(b, T);
+    }
+    else {
+        a = E[i - 1][j - 1];
+        if (a) {
+            e = use1_QU_Rosenfeld(a, T);
+        }
+        else {
+            d = E[i][j - 1];
+            if (d) {
+                e = use1_QU_Rosenfeld(d, T);
+            }
+            else {
+                e = ++ne;
+            }
+        }
+    }
+    E[i][j] = e;
+    return ne;
+}
+#endif // FAST
+
+
+#if FAST
+// -----------------------------------------------------------------------------------------
+static uint32 optimizedAccess_DT_Rosenfeld(uint32 ** E, int i, int j, uint32 * T, uint32 ne)
+// -----------------------------------------------------------------------------------------
+{
+    // Decision Tree 8-connexe avec Quick-Union
+    uint32 a, b, c, d, e;
+    
+    b = E[i - 1][j];
+    if (b) {
+        e = use1_QU_Rosenfeld(b, T);
+    }
+    else {
+        c = E[i - 1][j + 1];
+        if (c) {
+            a = E[i - 1][j - 1];
+            if (a) {
+                e = use2_QU_Rosenfeld(a, c, T);
+            }
+            else {
+                d = E[i][j - 1];
+                if (d) {
+                    e = use2_QU_Rosenfeld(c, d, T);
+                }
+                else {
+                    e = use1_QU_Rosenfeld(c, T);
+                }
+            }
+        }
+        else {
+            a = E[i - 1][j - 1];
+            if (a) {
+                e = use1_QU_Rosenfeld(a, T);
+            }
+            else {
+                d = E[i][j - 1];
+                if (d) {
+                    e = use1_QU_Rosenfeld(d, T);
+                }
+                else {
+                    e = ++ne;
+                }
+            }
+        }
+    }
+    E[i][j] = e;
+    return ne;
+}
+#endif // FAST
+
+
+
+#if FAST
+// --------------------------------------------------------------------------------------------------------
+static uint32 lineLabeling_Fast_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
+// --------------------------------------------------------------------------------------------------------
+{
+    uint8 x;
     // avec DT et QU
-    int j;
-    uint8 x;
-    
-    for (j = 0; j < width; j++) {
-        x = X[i][j];
+    // Left Border
+    x = X[i][0];
+    if (x) {
+        ne = optimizedAccessLeft_DT_Rosenfeld(E, i, 0, T, ne);
+    }
+    else {
+        E[i][0] = 0;
+    }
+    // Middle
+    for (int j = 1; j < width - 1; j++) {
+        uint8 x = X[i][j];
         if (x) {
             ne = optimizedAccess_DT_Rosenfeld(E, i, j, T, ne);
@@ -726 +1496 @@
         }
     }
+    // Right Border
+    x = X[i][width - 1];
+    if (x) {
+        ne = optimizedAccessRight_DT_Rosenfeld(E, i, width - 1, T, ne);
+    }
+    else {
+        E[i][width - 1] = 0;
+    }
     return ne;
 }
-
-
-// --------------------------------------------------------------------------------------------
-uint32 lineLabeling_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
-// --------------------------------------------------------------------------------------------
-{
+#endif // FAST
+
+
+// ---------------------------------------------------------------------------------------------------
+static uint32 lineLabeling_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
+// ---------------------------------------------------------------------------------------------------
+{
+#if SLOW
     return lineLabeling_Slow_Rosenfeld(X, i, width, E, T, ne);
-    //return lineLabeling_Fast_Rosenfeld(X, i, width, E, T, ne);
-}
-
-
-// ----------------------------------------------------------------
-uint32 countTable_Range_Rosenfeld(uint32 * T, uint32 e0, uint32 e1)
-// ----------------------------------------------------------------
+#elif FAST
+    return lineLabeling_Fast_Rosenfeld(X, i, width, E, T, ne);
+#else
+#error "Please define SLOW or FAST for the Rosenfeld version"
+#endif
+}
+
+
+// -----------------------------------------------------------------------
+static uint32 countTable_Range_Rosenfeld(uint32 * T, uint32 e0, uint32 e1)
+// -----------------------------------------------------------------------
 {
     uint32 e;
@@ -755 +1539 @@
 
 
-// --------------------------------------------------------------
-void solveTable_Range_Rosenfeld(uint32 * T, uint32 e0, uint32 e1)
-// --------------------------------------------------------------
+#if !FEATURES
+// ---------------------------------------------------------------------
+static void solveTable_Range_Rosenfeld(uint32 * T, uint32 e0, uint32 e1)
+// ---------------------------------------------------------------------
 {
     uint32 e, r;
@@ -766 +1551 @@
             T[e] = r; // racine de la classe d'equivalence
         }
-    }    
-}
-
-
+    }
+}
+#endif // !FEATURES
+
+
+#if FEATURES
+// ----------------------------------------------------------------------------------------------------------
+static void solveTable_solveFeatures_Range_Rosenfeld(uint32 * T, uint32 e0, uint32 e1, RegionStats * Stats)
+// ----------------------------------------------------------------------------------------------------------
+{
+    uint32 e, r;
+    
+    for (e = e0; e <= e1; e++) {
+        r = T[T[e]];
+        assert(r != 0);
+        if (r < e) {
+            T[e] = r; // racine de la classe d'equivalence
+            RegionStats_Accumulate_Stats1_From_Index(Stats, r, e);
+        }
+    }
+}
+#endif // FEATURES
+
+
+#if !FEATURES
 // -------------------------------------
 void MCA_Label_Rosenfeld_PAR1(MCA * mca)
@@ -775 +1581 @@
 {
     if (mca->p == 0) { 
-        MCA_VERBOSE1(printf("------------------------------\n"));
-        MCA_VERBOSE1(printf("-- MCA_Label_Rosenfeld_PAR1 --\n"));
-        MCA_VERBOSE1(printf("------------------------------\n"));
-    }
-    
+        printf("*** %s ***\n", __func__);
+    }
+    
+    CLOCK_THREAD_START_STEP(mca->p, 0);
+
     int i0 = mca->i0;
     int i1 = mca->i1;
@@ -795 +1601 @@
     if (mca->p == 0) {
         set_ui32vector_j(T, e0 - 1, e1); // car e0 = 1, on a besoin que T[0] = 0 pour FindRoot
-        // @QM : maintenant que c'est testé partout, en a-t-on encore besoin ? A priori non (a tester)
     }
     else {
@@ -809 +1614 @@
 
     MCA_VERBOSE2(display_ui32matrix_positive(E, i0, i1, 0, width - 1, 5, "Ep"); printf("\n"));
-    if (mca->p == 0) { 
+    if (mca->p == 0) {
         MCA_VERBOSE2(display_ui32vector_number(T, e0, ne, "%5d", "Tp_avant"));
     }
@@ -815 +1620 @@
     // fermeture transitive sans pack
     solveTable_Range_Rosenfeld(T, e0, ne);
-    nr = countTable_Range_Rosenfeld(T, e0, ne);
     mca->ne = ne; // Plus grande etiquette de l'intervalle [e0..e1]
 
+    MCA_VERBOSE2(nr = countTable_Range_Rosenfeld(T, e0, ne));
     MCA_VERBOSE2(printf("p = %d : e = [%d..%d] -> ne = %d -> nr = %d\n", mca->p, e0, ne, (ne - e0 + 1), nr));
-    if (mca->p == 0) { 
+    if (mca->p == 0) {
         MCA_VERBOSE2(display_ui32vector_number(T, e0, ne, "%5d", "Tp_apres"));
     }
-}
-
-
+    
+    CLOCK_THREAD_END_STEP(mca->p, 0);
+}
+#endif // !FEATURES
+
+
+#if !FEATURES
 // -------------------------------------
 void MCA_Label_Rosenfeld_PYR2(MCA * mca)
@@ -830 +1639 @@
 {
     // input
-    int np = mca->mca->np;
-    
-    // variables
-    int n = np;
-    int nb_level = i32log2(np);
-    if ((1 << nb_level) < np) {
-        nb_level++; // correction pour traiter n non puissance de 2
-    }
+    int p = mca->p;
+    int nb_level = mca->nb_level;
 
     if (mca->p == 0) {
-        MCA_VERBOSE1(printf("------------------------------\n"));
-        MCA_VERBOSE1(printf("-- MCA_Label_Rosenfeld_PYR2 --\n"));
-        MCA_VERBOSE1(printf("------------------------------\n"));
+        printf("*** %s ***\n", __func__);
     }
     
@@ -862 +1663 @@
     uint32 ** D = mca->D;
 
-    // @QM
-    // en fait, c'est compliqué.
-    // On pourrait optimiser en faisant faire un "break" aux procs qui n'ont plus jamais
-    // à faire d'itération, mais le problÚme est alors qu'il faut utiliser des barriÚres avec
-    // un nombre de procs à attendre différent à chaque fois, et qu'il faut les
-    // initialiser => il faut précalculer toutes ces valeurs et avoir une alloc dynamique
-    // du nombre de barriÚres.
-    // De plus, le problÚme est décuplé si le nombre de lignes n'est pas une puissance de 2, car
-    // dans ce cas certains threads ne doivent rien faire à une itération courante i,
-    // mais doivent être actifs à i + 1 => encore plus dur de calculer le nombre
-    // de threads à attendre à chaque barriÚre + surtout savoir s'il faut break ou continue
+    CLOCK_THREAD_START_STEP(p, 1);
+#if PYR_BARRIERS
+    // Version optimisée qui fait faire un break aux processeurs qui n'ont plus
+    // Ã  faire de merge.
+    // Implique de pré-calculer le nombre de threads à chaque barriÚre
+    if (p != 0) { // thread 0 never has any merge to do
+        int been_active = 0;
+        for (int level = 0; level < nb_level; level++) {
+            if ((p + (1 << level)) % (1 << (level + 1)) == 0) {
+                borderMerging_Rosenfeld_Dist(X, i, width, E, T, D, alpha);  // en (i) et (i-1)
+                been_active = 1;
+            }
+            else if (been_active) {
+                break;
+            }
+            pthread_barrier_wait(&mca->barriers[level]);
+        }
+    }
+    pthread_barrier_wait(&main_barrier);
+#else
     for (int level = 1; level <= nb_level; level++) {
-        if ((mca->p + (1 << (level - 1))) % (1 << level) == 0) {
+        if ((p + (1 << (level - 1))) % (1 << level) == 0) {
             // thread actif
-            //MCA_VERBOSE1(printf("### level = %d - p = %d\n", level, mca->p));
             borderMerging_Rosenfeld_Dist(X, i, width, E, T, D, alpha);  // en (i) et (i-1)
         }
-        barrier_wait(&main_barrier);
-    }
+        pthread_barrier_wait(&main_barrier);
+    }
+#endif
+    CLOCK_THREAD_END_STEP(p, 1);
     
 
@@ -887 +1698 @@
     // ---------------------------------
     
+    CLOCK_THREAD_START_STEP(p, 2);
     for (uint32 e = e0; e <= e1; e++) {
         uint32 r = T[e]; // acces local
         if (r < e) {
             r = FindRoot_Dist(D, e, alpha); // acces distant
-        }
-        T[e] = r;
-        MCA_VERBOSE2(printf("p%d : T[%d] <- %d\n", mca->p, e, r));
-    }
-}
+            T[e] = r; // @QM était en dehors du "if" (je pense que déjà demandé)
+        }
+        MCA_VERBOSE2(printf("p%d : T[%d] <- %d\n", p, e, r));
+    }
+    CLOCK_THREAD_END_STEP(p, 2);
+}
+#endif // !FEATURES
 
 
@@ -904 +1718 @@
     // input
     if (mca->p == 0) {
-        MCA_VERBOSE1(printf("------------------------------\n"));
-        MCA_VERBOSE1(printf("-- MCA_Label_Rosenfeld_PAR3 --\n"));
-        MCA_VERBOSE1(printf("------------------------------\n"));
+        printf("*** %s ***\n", __func__);
     }
     
@@ -917 +1729 @@
     uint32 * T = mca->T;
 
+    CLOCK_THREAD_START_STEP(mca->p, 3);
     for (int i = i0; i <= i1; i++) {
         for (int j = j0; j <= j1; j++) {
@@ -925 +1738 @@
         }
     }
-}
-
-
+    CLOCK_THREAD_END_STEP(mca->p, 3);
+}
+
+
+#if FEATURES
+// -----------------------------------------------------
+static void MCA_Label_Features_Rosenfeld_PAR1(MCA * mca)
+// -----------------------------------------------------
+{
+    if (mca->p == 0) { 
+        printf("*** %s ***\n", __func__);
+    }
+    
+    CLOCK_THREAD_START_STEP(mca->p, 0);
+
+    int i0 = mca->i0;
+    int i1 = mca->i1;
+    int width = mca->width;
+
+    uint32 e0 = mca->e0;
+    uint32 e1 = mca->e1;
+    uint32 ne = e0 - 1;
+    uint32 nr = 0;
+
+    // local memory zones
+    uint8 **  X = mca->X;
+    uint32 ** E = mca->E;
+    uint32 *  T = mca->T;
+
+    RegionStats * stats = mca->stats;
+
+    // reset sous optimal (pour le moment = voir region32)
+    if (mca->p == 0) {
+        set_ui32vector_j(T, e0 - 1, e1); // car e0 = 1, on a besoin que T[0] = 0 pour FindRoot
+        zero_RegionStatsVector(stats, e0 - 1, e1);
+    }
+    else {
+        set_ui32vector_j(T, e0, e1);
+        zero_RegionStatsVector(stats, e0, e1);
+    }
+
+    if (mca->p == 0) {
+        MCA_DISPLAY2(display_ui8matrix_positive(X, i0, i1, 0, width - 1, 5, "Xp"); printf("\n"));
+    }
+
+    // ---------------------------- //
+    // -- Etiquetage d'une bande -- //
+    // ---------------------------- //
+
+    ne = line0Labeling_Rosenfeld(X, i0, width, E, T, ne);
+    lineFeaturesComputation(E, i0, width, stats);
+
+    for (int i = i0 + 1; i <= i1; i++) {
+        ne = lineLabeling_Rosenfeld(X, i, width, E, T, ne); // Slow or Fast
+        lineFeaturesComputation(E, i, width, stats);
+    }
+    mca->ne = ne; //plus grande etiquette de l'intervalle [e0..e1]
+
+    if (mca->p == 0) {
+        MCA_VERBOSE2(printf("ne = %d\n", ne));
+        MCA_DISPLAY2(display_ui32matrix_positive(E, i0, i1, 0, width - 1, 5, "Ep"); printf("\n"));
+        MCA_DISPLAY2(display_ui32vector_number(T, e0, ne, "%5d", "Tp_avant"));
+    }
+
+    // ------------------------------------------------------ //
+    // -- Fermeture transitive sans pack de chaque table T -- //
+    // ------------------------------------------------------ //
+
+    solveTable_solveFeatures_Range_Rosenfeld(T, e0, ne, stats);
+
+    if (mca->p == 0) {
+        MCA_VERBOSE2(nr = countTable_Range_Rosenfeld(T, e0, ne);
+                printf("p = %d : e = [%d..%d] -> ne = %d -> nr = %d\n", mca->p, e0, ne, (ne - e0 + 1), nr));
+        MCA_DISPLAY2(display_ui32vector_number(T, e0, ne, "%5d", "Tp_apres"));
+    }
+    CLOCK_THREAD_END_STEP(mca->p, 0);
+}
+#endif // FEATURES
+
+
+#if FEATURES && !PARMERGE
+// -----------------------------------------------------
+static void MCA_Label_Features_Rosenfeld_PYR2(MCA * mca)
+// -----------------------------------------------------
+{
+    int p = mca->p;
+    int nb_level = mca->nb_level;
+
+    if (mca->p == 0) {
+        printf("*** %s ***\n", __func__);
+    }
+    
+    // ------------------------------
+    // -- pyramidal border merging --
+    // ------------------------------
+    
+    // local variables
+    int i = mca->i0;
+    int width = mca->width;
+    int alpha = mca->alpha;
+    uint32 e0 = mca->e0;
+    uint32 e1 = mca->ne;
+
+    // local memory zones
+    uint8 **  X = mca->X;
+    uint32 ** E = mca->E;
+    uint32 *  T = mca->T;
+    uint32 ** D = mca->D;
+    RegionStats ** F = mca->F;
+
+    CLOCK_THREAD_START_STEP(p, 1);
+#if PYR_BARRIERS
+    // Version optimisée qui fait faire un break aux processeurs qui n'ont plus
+    // Ã  faire de merge.
+    // Implique de pré-calculer le nombre de threads à chaque barriÚre
+    if (p != 0) { // thread 0 never has any merge to do
+        int been_active = 0;
+        for (int level = 0; level < nb_level; level++) {
+            if ((p + (1 << level)) % (1 << (level + 1)) == 0) {
+                borderMerging_Features_Rosenfeld_Dist(X, i, width, E, T, D, alpha, F);  // (i) et (i-1)
+                been_active = 1;
+            }
+            else if (been_active) {
+                break;
+            }
+            pthread_barrier_wait(&mca->barriers[level]);
+        }
+    }
+    pthread_barrier_wait(&main_barrier);
+#else
+    for (int level = 1; level <= nb_level; level++) {
+        if ((p + (1 << (level - 1))) % (1 << level) == 0) {
+            // thread actif
+            borderMerging_Features_Rosenfeld_Dist(X, i, width, E, T, D, alpha, F);  // (i) et (i-1)
+        }
+        pthread_barrier_wait(&main_barrier);
+    }
+#endif
+    CLOCK_THREAD_END_STEP(p, 1);
+
+
+    /**
+     * To remove?
+    // -- Affichage de debug
+    if (mca->p == 0) {
+        MCA_VERBOSE1(puts("-----------------------------"));
+        MCA_VERBOSE1(puts("[PYR2]: avant pack sequentiel"));
+        MCA_VERBOSE1(puts("-----------------------------"));
+    
+        for (int p = 0; p < mca->np; p++) {
+    
+            MCA* mca_par = mcas[p];
+            uint32 e0 = mca_par->e0;
+            uint32 e1 = mca_par->ne;
+            
+            uint32*  T = mca_par->T;
+            RegionStats* Stats = mca_par->Stats;
+        
+            RegionStats_DisplayStats_Sparse(T, e0, e1, Stats, NULL);
+            puts("");
+        }
+    }
+    */
+
+    // ---------------------------------
+    // -- parallel transitive closure --
+    // ---------------------------------
+    // identique a la version sans Features
+      
+    CLOCK_THREAD_START_STEP(p, 2);
+    for (uint32 e = e0; e <= e1; e++) {
+        uint32 r = T[e]; // acces local
+        if (r < e) {
+            r = FindRoot_Dist(D, e, alpha); // acces distant
+            T[e] = r;
+        }
+        MCA_VERBOSE2(printf("p%d : T[%d] <- %d\n", p, e, r));
+    }
+    CLOCK_THREAD_END_STEP(p, 2);
+
+    // To avoid uninitialized accesses
+    CLOCK_THREAD_START_STEP(p, 3);
+    CLOCK_THREAD_END_STEP(p, 3);
+}
+#endif // FEATURES && !PARMERGE
+
+
+#if FEATURES && PARMERGE
+// -----------------------------------------------------
+static void MCA_Label_Features_Rosenfeld_PAR2(MCA * mca)
+// -----------------------------------------------------
+{
+    int p = mca->p;
+    int nb_level = mca->nb_level;
+
+    if (mca->p == 0) {
+        printf("*** %s ***\n", __func__);
+    }
+    
+    // ------------------------------
+    // -- parallel border merging --
+    // ------------------------------
+    
+    // local variables
+    int i = mca->i0;
+    int width = mca->width;
+    int alpha = mca->alpha;
+    uint32 e0 = mca->e0;
+    uint32 e1 = mca->ne;
+
+    // local memory zones
+    uint8 **  X = mca->X;
+    uint32 ** E = mca->E;
+    uint32 *  T = mca->T;
+    uint32 ** D = mca->D;
+    RegionStats ** F = mca->F;
+
+    CLOCK_THREAD_START_STEP(p, 1);
+    if (p != 0) { // thread 0 never has any merge to do
+        borderMerging_Features_Rosenfeld_Dist(X, i, width, E, T, D, alpha, F);  // (i) et (i-1)
+    }
+    pthread_barrier_wait(&main_barrier);
+    CLOCK_THREAD_END_STEP(p, 1);
+
+
+    // ---------------------------------
+    // -- parallel transitive closure --
+    // ---------------------------------
+    // identique a la version sans Features
+     
+    CLOCK_THREAD_START_STEP(p, 2);
+    for (uint32 e = e0; e <= e1; e++) {
+        uint32 r = T[e]; // acces local
+        if (r < e) {
+            r = FindRoot_Dist(D, e, alpha); // acces distant
+            T[e] = r;
+        }
+        MCA_VERBOSE2(printf("p%d : T[%d] <- %d\n", p, e, r));
+    }
+    CLOCK_THREAD_END_STEP(p, 2);
+
+    // To avoid uninitialized accesses
+    CLOCK_THREAD_START_STEP(p, 3);
+    CLOCK_THREAD_END_STEP(p, 3);
+}
+#endif // FEATURES
+
+
+
+
+#if !FEATURES
 // =============================================================
+#if TARGET_OS == GIETVM
 __attribute__((constructor)) void MCA_Label_Rosenfeld(MCA * mca)
+#else
+void MCA_Label_Rosenfeld(MCA * mca)
+#endif
 // =============================================================
 {
+#if TARGET_OS == GIETVM
+    unsigned int x, y, lpid;
+    giet_proc_xyp(&x, &y, &lpid);
+    // Mettre à jour mca->p en fonction de x, y, lpid
+    // pour que les allocations faites par le main soient locales,
+    // i.e. 
+    mca->p = (x * Y_SIZE + y) * NB_PROCS_MAX + lpid;
+    // We have :
+    // mca->p = 4 pour (x = 0, y = 1, lpid = 0)
+    // mca->p = 5 pour (x = 0, y = 1, lpid = 1)
+    MCA_VERBOSE2(printf("mca->p = %d pour (x = %d, y = %d, lpid = %d)\n", mca->p, x, y, lpid));
+#endif
+
+    CLOCK_THREAD_START(mca->p);
+    CLOCK_THREAD_COMPUTE_START(mca->p);
+
     MCA_Scatter_ImageX(mca);
-    barrier_wait(&main_barrier);
+    pthread_barrier_wait(&main_barrier);
 
     MCA_Label_Rosenfeld_PAR1(mca);
-    barrier_wait(&main_barrier);
-    
-    //MCA_Gather_ImageL(mca);
-    //barrier_wait(&main_barrier);
-    //MCA_VERBOSE2(display_ui32matrix_positive(mca->E, mca->i0, mca->i1, 0, mca->width - 1, 5, "E2"));
-    //barrier_wait(&main_barrier);
-    
-    //MCA_Label_Rosenfeld_SEQ2(mca);
+    pthread_barrier_wait(&main_barrier);
+    
     MCA_Label_Rosenfeld_PYR2(mca);
-    barrier_wait(&main_barrier);
-    //MCA_VERBOSE2(display_ui32matrix_positive(mca->E, mca->i0, mca->i1, 0, mca->width - 1, 5, "EPYR"));
-    //barrier_wait(&main_barrier);
+    pthread_barrier_wait(&main_barrier);
     
     MCA_Label_Rosenfeld_PAR3(mca);
-    barrier_wait(&main_barrier);
+    pthread_barrier_wait(&main_barrier);
 
     MCA_Gather_ImageL(mca);
-    barrier_wait(&main_barrier);
-    //MCA_VERBOSE2(display_ui32matrix_positive(mca->E, mca->i0, mca->i1, 0, mca->width - 1, 5, "E3"));
-    //barrier_wait(&main_barrier);
-    
+    pthread_barrier_wait(&main_barrier);
+
+    CLOCK_THREAD_COMPUTE_END(mca->p);
+    CLOCK_THREAD_END(mca->p);
+
+#if TARGET_OS == GIETVM
     if (mca->p != 0) {
-        giet_pthread_exit(NULL);
-    }
-}
+        exit(0);
+    }
+#endif
+}
+#endif // !FEATURES
+
+
+#if FEATURES
+// ======================================================================
+#if TARGET_OS == GIETVM
+__attribute__((constructor)) void * MCA_Label_Features_Rosenfeld(void * arg)
+#else
+void * MCA_Label_Features_Rosenfeld(void * arg)
+#endif
+// ======================================================================
+{
+    MCA * mca = (MCA *) arg;
+#if TARGET_OS == GIETVM
+    unsigned int x, y, lpid;
+    giet_proc_xyp(&x, &y, &lpid);
+    // Mettre à jour mca->p en fonction de x, y, lpid
+    // pour que les allocations faites par le main soient locales,
+    // i.e. 
+    mca->p = (x * Y_SIZE + y) * NB_PROCS_MAX + lpid;
+    // We have :
+    // mca->p = 4 pour (x = 0, y = 1, lpid = 0)
+    // mca->p = 5 pour (x = 0, y = 1, lpid = 1)
+    MCA_VERBOSE2(printf("mca->p = %d pour (x = %d, y = %d, lpid = %d)\n", mca->p, x, y, lpid));
+#endif
+
+    CLOCK_THREAD_START(mca->p);
+    CLOCK_THREAD_COMPUTE_START(mca->p);
+
+    MCA_Scatter_ImageX(mca);
+    pthread_barrier_wait(&main_barrier);
+
+    MCA_Label_Features_Rosenfeld_PAR1(mca);
+    pthread_barrier_wait(&main_barrier);
+   
+#if PARMERGE 
+    MCA_Label_Features_Rosenfeld_PAR2(mca);
+#else
+    MCA_Label_Features_Rosenfeld_PYR2(mca);
+#endif
+    pthread_barrier_wait(&main_barrier);
+    
+    MCA_Label_Rosenfeld_PAR3(mca);
+    pthread_barrier_wait(&main_barrier);
+
+    MCA_Gather_ImageL(mca);
+    pthread_barrier_wait(&main_barrier);
+
+    CLOCK_THREAD_COMPUTE_END(mca->p);
+ 
+    if (display_features) {
+        if (mca->p == 0) {
+            int i = 1;
+            printf("[STATS]\n");
+            for (int p = 0; p < mca->np; p++) {
+                MCA * mca_par = mca->mca->mcas[p];
+                uint32 e0 = mca_par->e0;
+                uint32 ne = mca_par->ne - mca_par->e0; // number of elements
+                uint32 * T = mca_par->T;
+                RegionStats * stats = mca_par->stats;
+                RegionStats_DisplayStats_Sparse(T, e0, e0 + ne, stats, NULL, &i);
+            }
+            printf("[/STATS]\n");
+        }
+    }
+
+    CLOCK_THREAD_END(mca->p);
+
+#if TARGET_OS == GIETVM
+    if (mca->p != 0) {
+        exit(0);
+    }
+#endif
+
+    return NULL;
+}
+#endif // FEATURES
+
 
 // Local Variables: