Ignore:
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
Location:
soft/giet_vm/applications/rosenfeld/src-par
Files:
2 added
1 deleted
4 edited

Legend:

Unmodified
Added
Removed

  • soft/giet_vm/applications/rosenfeld/src-par/mca.c

    r805 r821  
    1515#include <malloc.h>
    1616
    17 #ifdef CLI
    1817#include "nrc_os_config.h"
     18#include "config.h"
    1919#include "nrc.h"
    20 #endif
    21 
     20
     21#if TARGET_OS == GIETVM
     22    #include <giet_config.h>
     23#endif
    2224
    2325#include "util.h"
    2426#include "ecc_common.h"
     27#include "ecc_features.h"
    2528#include "mca_matrix_dist.h"
    2629
    … …  
    3437
    3538
    36 // ----------------------
     39// -----------------------
    3740void MCA_Error(char * msg)
    38 // ----------------------
     41// -----------------------
    3942{
    4043    printf("MCA ERROR: %s\n", msg);
    4144    printf("now exiting to system...\n");
    42     exit(1);   
    43 }
    44 
    45 
    46 // ---------------------
    47 void MCA_Zero(MCA * mca)
    48 // ---------------------
    49 {
     45    exit(1);
    5046}
    5147
    … …  
    5854    mca = (MCA *) malloc(sizeof(MCA));
    5955
    60     if (mca) {
    61         MCA_Zero(mca);
    62     }
    63     else {
     56    if (mca == NULL) {
    6457        MCA_Error("allocation failed in MCA_pConstructor_Empty");
    6558    }
    … …  
    9184
    9285    mca->j0 = 0;
    93     mca->j1 = width-1;
     86    mca->j1 = width - 1;
    9487}
    9588
    … …  
    10295 
    10396    mca->i0 = 0;
    104     mca->i1 = height-1;
     97    mca->i1 = height - 1;
    10598}
    10699
    … …  
    161154    int j0_par, j1_par;
    162155    int height_par, height_mod;
    163    
     156
    164157    int pw2;
    165158    int32 ne_par;          // quantite par bande
    166159    uint32 nemax_par;      // la puissance de 2 >=
    167160    uint32 e0_par, e1_par; // indice par bande [start..end]
     161    int nb_level;
    168162   
    169163    MCA ** mcas;
    170164    MCA *  mca_par;
    171165   
    172     MCA_VERBOSE1(printf("====================\n"));
    173     MCA_VERBOSE1(printf("== MCA_Initialize ==\n"));
    174     MCA_VERBOSE1(printf("====================\n"));
    175    
    176     MCA_VERBOSE1(printf("height = %d\n", height));
    177     MCA_VERBOSE1(printf("width  = %d\n", width));
     166    printf("*** %s ***\n", __func__);
     167    MCA_VERBOSE1(printf("   height = %d\n", height));
     168    MCA_VERBOSE1(printf("   width  = %d\n", width));
    178169   
    179170    // array of pointers to mca workers
    … …  
    187178    height_par = height / np;
    188179    height_mod = height % np;
    189     //if(height % np) height_par++;
    190     MCA_VERBOSE1(printf("height_par = %d x %d + %d\n", height_par, np, height_mod));
    191    
    192     MCA_VERBOSE1(printf("========================\n"));
     180
     181    MCA_VERBOSE1(printf("   height_par = %d x %d + %d\n", height_par, np, height_mod));
     182    MCA_VERBOSE1(printf("   ========================\n"));
    193183   
    194184    i1_par_previous = 0;
    195    
    196     for (int p = 0; p < np; p++) {
     185
     186    // puissance de 2 de chaque bande
     187    ne_par = height_par * width + 1;
     188    MCA_VERBOSE1(printf("   ne_par    = %d\n", ne_par));
     189    pw2 = i32log2(ne_par);
     190    if (ne_par > (1 << pw2)) {
     191        pw2++;
     192    }
     193    nemax_par = 1 << pw2;
     194
     195    MCA_VERBOSE1(printf("   nemax_par = %d\n", nemax_par));
     196
     197    nb_level = i32log2(np);
     198    if ((1 << nb_level) < np) {
     199        nb_level++;
     200    }
     201
     202#if PYR_BARRIERS
     203    // ------------------------------------------
     204    // -- Allocation des barriÚres pyramidales --
     205    // ------------------------------------------
     206
     207    pthread_barrier_t * barriers = NULL;
     208    if (nb_level > 0) {
     209        barriers = malloc(sizeof(pthread_barrier_t) * nb_level);
     210
     211        // Initially all threads are active except thread 0
     212        int nb_active = np - 1;
     213        pthread_barrier_init(&barriers[0], NULL, nb_active);
     214        for (int i = 1; i < nb_level; i++) {
     215            // thread 0 never does any merge
     216            for (int p = 1; p < np; p++) {
     217                if ((p + (1 << (i - 1))) % (1 << i) == 0) {
     218                    // thread inactive at level i
     219                    nb_active -= 1;
     220                }
     221            }
     222            pthread_barrier_init(&barriers[i], NULL, nb_active);
     223        }
     224    }
     225#endif
     226
     227    for (int p = 0; p < np; p++) {
     228
    197229        // ----------------- //
    198230        // -- constructor -- //
    199231        // ----------------- //
    200         MCA_VERBOSE1(printf("-- p = %d ----------------\n", p));
     232        MCA_VERBOSE2(printf("-- p = %d ----------------\n", p));
    201233   
    202234        // alloc of mca workers into array of pointers
    … …  
    208240        mca_par->p   = p;
    209241        mca_par->mca = mca; // pointer to master
     242#if TARGET_OS == GIETVM
     243        int x, y; // cluster coordinates
     244        // We have p == 4 => x = 0; y = 1
     245        x = (p / NB_PROCS_MAX) / Y_SIZE;
     246        y = (p / NB_PROCS_MAX) % Y_SIZE;
     247        MCA_VERBOSE2(printf("p = %d (x = %d, y = %d)\n", p, x, y));
     248#endif
    210249       
    211250        // ------------------------------------- //
    … …  
    214253       
    215254        // hauteur de chaque bande
    216        
    217         //printf("i1_par_previous = %d\n", i1_par_previous);
    218        
    219255        if (p == 0) {
    220256            i0_par = 0;
    … …  
    232268        i1_par_previous = i1_par;
    233269       
    234         MCA_VERBOSE1(printf("i0_par = %d\n", i0_par));
    235         MCA_VERBOSE1(printf("i1_par = %d\n", i1_par));
    236        
    237         // puissance de 2 de chaque bande
    238         ne_par = height_par * width + 1;
    239         //if (p == 0) {
    240         //    ne_par++;
    241         //}
    242         MCA_VERBOSE1(printf("ne_par    = %d\n", ne_par));
    243         pw2 = i32log2(ne_par);
    244         if (ne_par > (1 << pw2)) {
    245             pw2++;
    246         }
    247         nemax_par = 1 << pw2;
    248         MCA_VERBOSE1(printf("nemax_par = %d\n", nemax_par));
     270        MCA_VERBOSE2(printf("i0_par = %d\n", i0_par));
     271        MCA_VERBOSE2(printf("i1_par = %d\n", i1_par));
    249272       
    250273        // etiquettes
    … …  
    260283        MCA_VERBOSE2(printf("e0_par = %d\n", e0_par));
    261284        MCA_VERBOSE2(printf("e1_par = %d\n", e1_par));
    262        
     285
    263286        mca_par->width  = width;
    264287        mca_par->height = height_par;
    265    
    266288        mca_par->i0 = i0_par;
    267289        mca_par->i1 = i1_par;
    268290        mca_par->j0 = 0;
    269291        mca_par->j1 = width - 1;
    270    
    271292        mca_par->e0 = e0_par;
    272293        mca_par->e1 = e1_par;
    273        
    274294        mca_par->alpha = pw2;
     295        mca_par->np = np;
     296        // Pour les barriÚres pyramidales
     297        mca_par->nb_level = nb_level;
     298#if PYR_BARRIERS
     299        mca_par->barriers = barriers;
     300#else
     301        mca_par->barriers = NULL;
     302#endif
     303        mca_par->F = NULL; // default init
     304        mca_par->stats = NULL; // default init
    275305 
    276306        // ---------------- //
    277307        // -- allocation -- //
    278308        // ---------------- //
    279    
     309#if TARGET_OS == GIETVM
     310        mca_par->X = remote_ui8matrix(i0_par, i1_par, 0, width - 1, x, y);
     311        mca_par->E = remote_dist_ui32matrix(i0_par, i1_par, 0, width - 1, x, y); // distributed matrix with border
     312       
     313        if (p == 0) {
     314            mca_par->T = remote_ui32vector(e0_par - 1, e1_par, x, y); // car e0 = 1, on a besoin que T[0] = 0 pour FindRoot
     315#if FEATURES
     316            mca_par->stats = remote_RegionStatsVector(e0_par - 1, e1_par, x, y);
     317#endif
     318        }
     319        else {
     320            mca_par->T = remote_ui32vector(e0_par, e1_par, x, y);
     321#if FEATURES
     322            mca_par->stats = remote_RegionStatsVector(e0_par, e1_par, x, y);
     323#endif
     324        }
     325       
     326        mca_par->D = (uint32 **) remote_vvector(0, np - 1, x, y);
     327#if FEATURES
     328        mca_par->F = (RegionStats **) remote_vvector(0, np - 1, x, y);
     329#endif
     330#else // !GIETVM
    280331        mca_par->X = ui8matrix (i0_par, i1_par, 0, width - 1);
    281332        mca_par->E = dist_ui32matrix(i0_par, i1_par, 0, width - 1); // distributed matrix with border
    … …  
    283334        if (p == 0) {
    284335            mca_par->T = ui32vector(e0_par - 1, e1_par); // car e0 = 1, on a besoin que T[0] = 0 pour FindRoot
     336#if FEATURES
     337            mca_par->stats = RegionStatsVector(e0_par - 1, e1_par);
     338
     339#endif
    285340        }
    286341        else {
    287342            mca_par->T = ui32vector(e0_par, e1_par);
     343#if FEATURES
     344            mca_par->stats = RegionStatsVector(e0_par, e1_par);
     345#endif
    288346        }
    289347       
    290348        mca_par->D = (uint32 **) vvector(0, np - 1);
    291        
     349#if FEATURES
     350        mca_par->F = (RegionStats **) vvector(0, np - 1);
     351#endif
     352#endif   
    292353        MCA_VERBOSE2(printf("X = %p\n", mca_par->X));
    293354        MCA_VERBOSE2(printf("E = %p\n", mca_par->E));
    294355        MCA_VERBOSE2(printf("T = %p\n", mca_par->T));
    295356        MCA_VERBOSE2(printf("D = %p\n", mca_par->D));
    296    
    297357    } // p
    298358   
    299     // pour debug
    300     MCA_VERBOSE2(printf("init des tables d'EQ a l'identite\n"));
    301     for (int p = 0; p < np; p++) {
    302    
     359
     360    for (int p = 0; p < np; p++) {
    303361        MCA * mca_par = mcas[p];
    304362       
    … …  
    319377    MCA_VERBOSE2(printf("display des tables d'EQ\n"));
    320378    for (int p = 0; p < np; p++) {
    321        
    322379        MCA * mca_par = mcas[p];
    323380       
    … …  
    326383        uint32 e1 = mca_par->e1;
    327384       
    328         MCA_VERBOSE1(printf("p = %d T[%d..%d]\n", p, e0, e1));
     385        MCA_VERBOSE2(printf("p = %d T[%d..%d]\n", p, e0, e1));
    329386        if (p == 0) {
    330             MCA_VERBOSE1(display_ui32vector_number(T, e0 - 1, e0 + 10, "%5d", "T"));
    331         }
    332         else {
    333             MCA_VERBOSE1(display_ui32vector_number(T, e0, e0 + 10, "%5d", "T"));
     387            MCA_VERBOSE2(display_ui32vector_number(T, e0 - 1, e0 + 10, "%5d", "T"));
     388        }
     389        else {
     390            MCA_VERBOSE2(display_ui32vector_number(T, e0, e0 + 10, "%5d", "T"));
    334391        }
    335392        MCA_VERBOSE2(printf("\n"));
    … …  
    343400    // table d'indirection distribuee D
    344401    MCA_VERBOSE2(printf("nemax_par = %d\n", nemax_par));
    345    
    346     for (int p = 0; p < np; p++) {
    347    
     402    for (int p = 0; p < np; p++) {
    348403        MCA * mca_p = mcas[p];
    349404        uint32 ** D = mca_p->D;
     405        RegionStats ** F  = mca_p->F;
    350406       
    351407        for (int k = 0; k < np; k++) {
    352             //mcas[p]->D[k] = (void*) (&(mcas[k]->T[mcas[k]->e0]))  - mcas[k]->e0; //k * nemax_par;
    353            
    354408            MCA * mca_k = mcas[k];
    355409            uint32 * T = mca_k->T;
    356            
    357410            D[k] = T + k * nemax_par; // il faut soustraire le "MSB"
     411#if FEATURES
     412            RegionStats * stat = mca_k->stats;
     413            F[k] = stat + k * nemax_par; // il faut soustraire le "MSB"
     414#endif
    358415        } // k
    359416    } // p
    … …  
    376433            MCA_VERBOSE2(display_ui32vector(D[k], 0, 9, "%5d", "D\n"));
    377434        }
    378         printf("\n");
     435        MCA_VERBOSE2(printf("\n"));
    379436    }
    380437   
    … …  
    413470// -----------------------------------
    414471{
    415     int p, np = mca->np;
     472    int np = mca->np;
    416473   
    417474    MCA ** mcas = mca->mcas;
    418475    MCA *  mca_par;
    419    
    420     printf("============================\n");
    421     printf("== MCA_Display_Parameters ==\n");
    422     printf("============================\n");
    423    
    424     printf("height = %d\n", mca->height);
    425     printf("width  = %d\n", mca->width);
    426     printf("np     = %d\n", mca->np);
    427    
    428     for (p = 0; p < np; p++) {
     476    (void) mca_par;
     477   
     478    printf("*** MCA_Display_Parameters ***\n");
     479   
     480    MCA_VERBOSE1(printf("   height = %d\n", mca->height));
     481    MCA_VERBOSE1(printf("   width  = %d\n", mca->width));
     482    MCA_VERBOSE1(printf("   np     = %d\n", mca->np));
     483   
     484    for (int p = 0; p < np; p++) {
    429485        mca_par = mcas[p];
    430486       
    431         printf("Display MCA[%d]\n", p);
    432         printf("p = %d\n", mca_par->p);
    433         printf("i0 = %8d  i1 = %8d\n", mca_par->i0, mca_par->i1);
    434         printf("j0 = %8d  j1 = %8d\n", mca_par->j0, mca_par->j1);
    435         printf("e0 = %8d  e1 = %8d\n", mca_par->e0, mca_par->e1);
     487        MCA_VERBOSE2(printf("Display MCA[%d]\n", p));
     488        MCA_VERBOSE2(printf("p = %d\n", mca_par->p));
     489        MCA_VERBOSE2(printf("i0 = %8d  i1 = %8d\n", mca_par->i0, mca_par->i1));
     490        MCA_VERBOSE2(printf("j0 = %8d  j1 = %8d\n", mca_par->j0, mca_par->j1));
     491        MCA_VERBOSE2(printf("e0 = %8d  e1 = %8d\n", mca_par->e0, mca_par->e1));
    436492    }
    437493}
    … …  
    442498// -------------------------
    443499{
    444     int p, np = mca->np;
     500    int np = mca->np;
    445501   
    446502    MCA ** mcas = mca->mcas;
    … …  
    451507    uint32 e0, e1;
    452508   
    453     printf("==================\n");
    454     printf("== MCA_Finalize ==\n");
    455     printf("==================\n");
    456    
    457     for (p = 0; p < np; p++) {
    458    
     509    printf("*** MCA_Finalize ***\n");
     510   
     511#if PYR_BARRIERS
     512    free(mcas[0]->barriers);
     513#endif
     514
     515    for (int p = 0; p < np; p++) {
    459516        mca_par = mcas[p];
    460517   
    … …  
    475532        if (p == 0) {
    476533            free_ui32vector(mca_par->T, e0 - 1, e1); // car e0 = 1, on a besoin que T[0] = 0 pour FindRoot
     534#if FEATURES
     535            free_RegionStatsVector(mca_par->stats, e0 - 1, e1);
     536#endif
    477537        }
    478538        else {
    479539            free_ui32vector(mca_par->T, e0, e1);
     540#if FEATURES
     541            free_RegionStatsVector(mca_par->stats, e0, e1);
     542#endif
    480543        }
    481544       
    482545        free_vvector((void **) mca_par->D, 0, np - 1);
    483        
    484     }
    485     printf("[MCA_Finalize]: fin boucle\n");
     546#if FEATURES
     547        free_vvector((void **) mca_par->F, 0, np - 1);
     548#endif
     549        free(mca_par);
     550    }
    486551    free(mcas);
    487     printf("[MCA_Finalize]: mcas freed\n");
    488     free(mca); // plante si free XET, ne plante pas si pas de free ...
    489     printf("[MCA_Finalize]: mca freed\n");
    490 }
    491 
    492 
    493 // @QM check my modifs...
     552    free(mca);
     553}
     554
     555
    494556// -------------------------------
    495557void MCA_Scatter_ImageX(MCA * mca)
    … …  
    498560    // diffusion de l'image binaire source
    499561   
    500     int      np = mca->np;
     562    int np = mca->np;
    501563    uint8 ** X  = mca->mca->X;
    502564   
    503565    if (mca->p == 0) {
    504         MCA_VERBOSE1(printf("------------------------\n"));
    505         MCA_VERBOSE1(printf("-- MCA_Scatter_ImageX --\n"));
    506         MCA_VERBOSE1(printf("------------------------\n"));
     566        printf("*** MCA_Scatter_ImageX ***\n");
    507567    }
    508568   
    … …  
    524584
    525585
    526 // @QM check my modifs...
    527586// ------------------------------
    528587void MCA_Gather_ImageL(MCA * mca)
    … …  
    532591    int np = mca->np;
    533592    uint32 ** E = mca->mca->E;
     593
     594    if (mca->p == 0) {
     595        printf("*** MCA_Gather_ImageL ***\n");
     596    }
    534597
    535598    int i0 = mca->i0;

  • soft/giet_vm/applications/rosenfeld/src-par/mca_main.c

    r805 r821  
    11/* ------------------ */
    2 /* --- mca.c --- */
     2/* --- mca_main.c --- */
    33/* ------------------ */
    44
    … …  
    1212#include <string.h>
    1313#include <math.h>
    14 
    15 #include <user_lock.h>
    16 
    17 #ifdef CLI
     14#include <malloc.h>
     15
    1816#include "nrc_os_config.h"
     17#include "config.h"
    1918#include "nrc.h"
    20 #endif
    21 
     19
     20#if TARGET_OS == GIETVM
     21    #include <user_lock.h>
     22    #include <malloc.h>
     23    #include <giet_config.h>
     24    #include <user_barrier.h>
     25#else
     26    #include <unistd.h>
     27#endif
    2228
    2329#include "util.h"
    2430#include "ecc_common.h"
    2531#include "ecc_features.h"
    26 
    2732#include "palette.h"
    2833#include "bmpNR.h"
    29 
     34#include "mca_matrix_dist.h"
     35#include "mca_rosenfeld.h"
     36#include "clock.h"
    3037#include "str_ext.h"
     38
    3139
    3240/* -- local -- */
    3341#include "mca.h"
    34 #include "mca_test.h"
    35 
    36 
     42
     43#define MAX_THREADS 256
     44#define DEFAULT_NTHREADS 1
     45#define DEFAULT_IN_FILENAME "/misc/cadastre.pgm"
     46#define DEFAULT_OUT_FILENAME "out.bmp"
     47
     48pthread_t thread_table[MAX_THREADS];
     49pthread_barrier_t main_barrier;
     50int display_features = 0;
     51int generate_output_image = 0;
     52
     53CLOCK_DEC;
     54
     55static void usage(char * name) {
     56    printf("Usage: %s <options>\n", name);
     57    printf("options:\n");
     58    printf("  -i <input_file>  : Input file (default = %s)\n", DEFAULT_IN_FILENAME);
     59    printf("  -o <output_file> : Output file (default = %s)\n", DEFAULT_OUT_FILENAME);
     60    printf("  -nN              : N = number of threads (default = %d).\n", DEFAULT_NTHREADS);
     61    printf("  -d               : Display features (default = false, requires features computation).\n");
     62    printf("  -g               : Generate output image (default = false).\n");
     63    printf("  -h               : Print out command line options.\n\n");
     64}
     65
     66
     67
     68// --------------------------------------------------------------------------
     69void init_forme_boulon1(uint8 *** X0, int * i0, int * i1, int * j0, int * j1)
     70// --------------------------------------------------------------------------
     71{
     72    uint8 ** X;
     73    int i =  0;
     74    int h =  28;
     75    int w =  30;
     76   
     77    X = ui8matrix(0, h - 1, 0, w - 1);
     78    zero_ui8matrix(X, 0, h - 1, 0, w - 1);
     79   
     80    *X0 = X;
     81    *i0 = 0;
     82    *i1 = h - 1;
     83    *j0 = 0;
     84    *j1 = w - 1;
     85   
     86    //                                 0000000001111111111122222222223
     87    //                                 0123456789012345678901234567890
     88    set_ui8vector_str(X[i++], 0, w - 1, "                         111  "); // 00
     89    set_ui8vector_str(X[i++], 0, w - 1, "                        11111 "); // 01
     90    set_ui8vector_str(X[i++], 0, w - 1, "                      1111111 "); // 02
     91    set_ui8vector_str(X[i++], 0, w - 1, "                     11111111 "); // 03
     92    set_ui8vector_str(X[i++], 0, w - 1, "                    1111111111"); // 04
     93    set_ui8vector_str(X[i++], 0, w - 1, "                   11111111111"); // 05
     94    set_ui8vector_str(X[i++], 0, w - 1, "                 1111111111111"); // 06
     95    set_ui8vector_str(X[i++], 0, w - 1, "               11111111111111 "); // 07
     96    set_ui8vector_str(X[i++], 0, w - 1, "              11111111111111  "); // 08
     97    set_ui8vector_str(X[i++], 0, w - 1, "             11111111111111   "); // 09
     98    set_ui8vector_str(X[i++], 0, w - 1, "     11    11111111111111     "); // 10
     99    set_ui8vector_str(X[i++], 0, w - 1, "    111   11111111111111      "); // 11
     100    set_ui8vector_str(X[i++], 0, w - 1, "   11111111111111111111       "); // 12
     101    set_ui8vector_str(X[i++], 0, w - 1, " 11111111111111111111         "); // 13
     102    set_ui8vector_str(X[i++], 0, w - 1, "1111111111111111111           "); // 14
     103    set_ui8vector_str(X[i++], 0, w - 1, " 11111111111111111            "); // 15
     104    set_ui8vector_str(X[i++], 0, w - 1, " 1111111111111111             "); // 16
     105    set_ui8vector_str(X[i++], 0, w - 1, " 111111111111111              "); // 17
     106    set_ui8vector_str(X[i++], 0, w - 1, "  111111111111                "); // 18
     107    set_ui8vector_str(X[i++], 0, w - 1, "  1111111111                  "); // 29
     108    set_ui8vector_str(X[i++], 0, w - 1, "  1111111111                  "); // 20
     109    set_ui8vector_str(X[i++], 0, w - 1, "   111111111                  "); // 21
     110    set_ui8vector_str(X[i++], 0, w - 1, "   111111111                  "); // 22
     111    set_ui8vector_str(X[i++], 0, w - 1, "    11111111                  "); // 23
     112    set_ui8vector_str(X[i++], 0, w - 1, "    1111111                   "); // 24
     113    set_ui8vector_str(X[i++], 0, w - 1, "     11111                    "); // 25
     114    set_ui8vector_str(X[i++], 0, w - 1, "     111                      "); // 26
     115    set_ui8vector_str(X[i++], 0, w - 1, "                              "); // 27
     116   
     117    //printf("[init_forme_boulon1]: h = %d i = %d\n", h, i);
     118    if (i != h) {
     119        MCA_Error("init_forme_boulon1 i != h");
     120    }
     121
     122   
     123    //display_ui8matrix_positive(X, 0, h-1, 0, w-1, 4, "forme_boulon1"); printf("");
     124    //write_ui8matrix_positive(  X, 0, h-1, 0, w-1, 4, "forme_boulon1.txt");
     125}
     126
     127
     128// QM : The cost of this function is horrible
     129// but it is only for testing purpose
     130// Renumbers object in a contiguous way, for an image which has already
     131// been processed with several threads
     132// --------------------------------------------------------------------
     133static void renumber_image(uint32 ** E, int i0, int i1, int j0, int j1)
     134// --------------------------------------------------------------------
     135{
     136    int size = 10;
     137    int idx = 1; // next label to give, first invalid index in the equiv table
     138    uint32 * equiv = malloc(sizeof(uint32) * size);
     139    equiv[0] = 0; // unused
     140    int found;
     141
     142    for (int i = i0; i <= i1; i++) {
     143        for (int j = j0; j <= j1; j++) {
     144            if (E[i][j] != 0) {
     145                found = 0;
     146                for (int k = 1; k < idx; k++) {
     147                    if (equiv[k] == E[i][j]) {
     148                        E[i][j] = k;
     149                        found = 1;
     150                        break;
     151                    }
     152                }
     153                if (found == 0) {
     154                    equiv[idx] = E[i][j];
     155                    E[i][j] = idx;
     156                    idx += 1;
     157                    if (idx == size) {
     158                        size = size * 2;
     159                        equiv = realloc(equiv, sizeof(uint32) * size);
     160                    }
     161                }
     162            }
     163        }
     164    }
     165    free(equiv);
     166}
     167
     168
     169// ----------------------------
     170void mca_test1(int num_threads)
     171// ----------------------------
     172{
     173    int i0, i1, j0, j1;
     174    int height, width;
     175   
     176    uint8 ** X0;
     177    uint32 ** E;
     178    MCA * mca;
     179
     180    pthread_barrier_init(&main_barrier, NULL, num_threads);
     181
     182    // -- Allocation --
     183    init_forme_boulon1(&X0, &i0, &i1, &j0, &j1);
     184   
     185    height = i1 - i0 + 1;
     186    width  = j1 - j0 + 1;
     187   
     188    E = ui32matrix(i0, i1, j0, j1);
     189   
     190    zero_ui32matrix(E, i0, i1, j0, j1);
     191   
     192    mca = MCA_pConstructor_Empty();
     193   
     194    // -- set param
     195    MCA_Set_Size(mca, width, height);
     196    MCA_Set_ImageX(mca, X0);
     197    MCA_Set_ImageL(mca, E);
     198    MCA_Set_NP(mca, num_threads);
     199   
     200    // -- MCA init
     201    MCA_Initialize(mca);
     202    MCA_Display_Parameters(mca);
     203   
     204    display_ui8matrix_positive(mca->X, i0, i1, j0, j1, 5, "X0");
     205#if FEATURES
     206    for (int i = 1; i < num_threads; i++) {
     207        pthread_create(&thread_table[i], NULL, MCA_Label_Features_Rosenfeld, (void *) mca->mcas[i]);
     208    }
     209    MCA_Label_Features_Rosenfeld(mca->mcas[0]);
     210#else
     211    for (int i = 1; i < num_threads; i++) {
     212        pthread_create(&thread_table[i], NULL, MCA_Label_Rosenfeld, (void *) mca->mcas[i]);
     213    }
     214    MCA_Label_Rosenfeld(mca->mcas[0]);
     215#endif
     216    for (int i = 1; i < num_threads; i++) {
     217        pthread_join(thread_table[i], NULL);
     218    }
     219    display_ui32matrix_positive(mca->E, i0, i1, j0, j1, 5, "Efinal");
     220
     221   
     222    // -- free --
     223    printf("Finalize\n");
     224    MCA_Finalize(mca);
     225   
     226    printf("Free_matrix\n");
     227    free_ui8matrix (X0, i0, i1, j0, j1);
     228    free_ui32matrix(E,  i0, i1, j0, j1);
     229}
     230
     231
     232
     233// -----------------------------------------------------------
     234void mca_test2(int num_threads, char * infile, char * outfile)
     235// -----------------------------------------------------------
     236{
     237    int i0, i1, j0, j1;
     238    int height, width;
     239   
     240    uint8 ** X;
     241    uint8 ** E8;
     242    uint32 ** E;
     243    MCA * mca;
     244
     245    RGBQuad palette[256];
     246
     247    pthread_barrier_init(&main_barrier, NULL, num_threads);
     248
     249    Palette_18ColorsBW(palette);
     250   
     251    printf("Loading file %s... ", infile);
     252    X = LoadPGM_ui8matrix(infile, &i0, &i1, &j0, &j1);
     253    printf("done.\n");
     254
     255    printf("Allocating memory... ");
     256    height = i1 - i0 + 1;
     257    width  = j1 - j0 + 1;
     258   
     259    E8 = ui8matrix (i0, i1, j0, j1);
     260    E  = ui32matrix(i0, i1, j0, j1);
     261   
     262    zero_ui8matrix(E8, i0, i1, j0, j1);
     263    zero_ui32matrix(E, i0, i1, j0, j1);
     264
     265    // pre-traitements
     266    binarisation_ui8matrix(X, i0, i1, j0, j1, 20, 1, X); // pour le traitement
     267    printf("done.\n");
     268
     269    printf("Allocating and initializing MCA... \n");
     270    mca = MCA_pConstructor_Empty();
     271   
     272    // -- set param
     273    MCA_Set_Size(mca, width, height);
     274    MCA_Set_ImageX(mca, X);
     275    MCA_Set_ImageL(mca, E);
     276    MCA_Set_NP(mca, num_threads);
     277   
     278    // -- MCA init
     279    MCA_Initialize(mca);
     280    MCA_Display_Parameters(mca);
     281    printf("End of MCA allocation and initialization.\n");
     282   
     283    CLOCK_APP_CREATE;
     284#if FEATURES
     285    for (int i = 1; i < num_threads; i++) {
     286        pthread_create(&thread_table[i], NULL, MCA_Label_Features_Rosenfeld, (void *) mca->mcas[i]);
     287    }
     288    MCA_Label_Features_Rosenfeld(mca->mcas[0]);
     289#else
     290    for (int i = 1; i < num_threads; i++) {
     291        pthread_create(&thread_table[i], NULL, MCA_Label_Rosenfeld, (void *) mca->mcas[i]);
     292    }
     293    MCA_Label_Rosenfeld(mca->mcas[0]);
     294#endif
     295    for (int i = 1; i < num_threads; i++) {
     296        pthread_join(thread_table[i], NULL);
     297    }
     298    CLOCK_APP_JOIN;
     299
     300    if (generate_output_image) {
     301#if TARGET_OS != GIETVM
     302        renumber_image(mca->E, i0, i1, j0, j1);
     303#else
     304        printf("Warning: the output image has not been renumbered, it cannot be used as a comparison with the reference\n");
     305#endif
     306        mod_ui32matrix_ui8matrix(mca->E, i0, i1, j0, j1, E8);
     307        printf("Saving file %s for verification... ", outfile);
     308        SaveBMP2_ui8matrix(E8, width, height, palette, outfile);
     309        printf("done.\n");
     310    }
     311
     312    MCA_Finalize(mca);
     313    printf("Deallocating memory...");
     314    free_ui8matrix (X,  i0, i1, j0, j1);
     315    free_ui8matrix (E8, i0, i1, j0, j1);
     316    free_ui32matrix(E,  i0, i1, j0, j1);
     317    printf("done.\n");
     318}
     319
     320
     321// --------------------------------------------------------------
     322int main_test_mca(int num_threads, char * infile, char * outfile)
     323// --------------------------------------------------------------
     324{
     325    CLOCK_INIT(num_threads, 4); // 4 = Number of steps in body
     326    CLOCK_APP_START;
     327
     328    mca_test2(num_threads, infile, outfile);
     329
     330    CLOCK_APP_END;
     331    CLOCK_FINALIZE;
     332    PRINT_CLOCK;
     333    CLOCK_FREE;
     334   
     335    return 0;
     336}
     337
     338
     339#if TARGET_OS == GIETVM
     340// ------------------------------------
     341__attribute__((constructor)) int main()
     342// ------------------------------------
     343#else
    37344// -----------------------------
    38 __attribute__((constructor)) void main()
     345int main(int argc, char ** argv)
    39346// -----------------------------
    40 {
     347#endif
     348{
     349    char * infile = DEFAULT_IN_FILENAME;
     350    char * outfile = DEFAULT_OUT_FILENAME;
     351
     352    int ch;
     353    int num_threads = DEFAULT_NTHREADS;
     354
     355    printf("*** Starting application Rosenfeld ***\n");
     356
     357#if TARGET_OS != GIETVM // @QM I think the giet has some random (uninitialized) values for argc and argv
     358    while ((ch = getopt(argc, argv, "i:o:n:hdg")) != EOF) {
     359        switch (ch) {
     360        case 'i':
     361            infile = optarg;
     362            break;
     363        case 'o':
     364            outfile = optarg;
     365            break;
     366        case 'n':
     367            num_threads = atoi(optarg);
     368            break;
     369        case 'h':
     370            usage(argv[0]);
     371            return 0;
     372            break;
     373        case 'd':
     374#if !FEATURES
     375            fprintf(stderr, "*** Error: Features display requires features computation\n");
     376            return 1;
     377#endif
     378            display_features = 1;
     379            break;
     380        case 'g':
     381            generate_output_image = 1;
     382            break;
     383        default:
     384            usage(argv[0]);
     385            return 1;
     386            break;
     387        }
     388    }
     389
     390    // Check arguments
     391    if (num_threads < 1) {
     392        fprintf(stderr, "*** Error: The number of threads must at least be 1\n");
     393        usage(argv[0]);
     394        return -1;
     395    }
     396#endif
     397
     398#if TARGET_OS == GIETVM
     399    {
     400        unsigned int xsize, ysize, nprocs;
     401        giet_procs_number(&xsize, &ysize, &nprocs);
     402        num_threads = xsize * ysize * nprocs;
     403    }
     404#endif
     405
     406    if (num_threads > MAX_THREADS) {
     407        printf("*** Error: The maximum number of threads is %d, i.e. less than the current number of threads.\n", MAX_THREADS);
     408        printf("Please recompile with a bigger MAX_THREADS value.\n");
     409        exit(1);
     410    }
     411
     412    printf("Parameters:\n");
     413    printf("- Number of threads: %d\n", num_threads);
     414    printf("- Input file: %s\n", infile);
     415    printf("- Output file: %s\n", outfile);
     416#if FAST
     417    printf("- Using decision trees (fast): yes\n");
     418#elif SLOW
     419    printf("- Using decision trees (fast): no\n");
     420#endif
     421#if FEATURES
     422    printf("- Computing features: yes\n");
     423#else
     424    printf("- Computing features: no\n");
     425#endif
     426
     427
    41428#if TARGET_OS == GIETVM
    42429    giet_tty_alloc(1);
    43     heap_init(0, 0);
    44 #endif
    45     lock_init(&print_lock);
    46     main_test_mca();
    47 
    48     exit(0);
    49 }
    50 
     430    printf("Initializing heaps... ");
     431    for (int i = 0; i < X_SIZE; i++) {
     432        for (int j = 0; j < X_SIZE; j++) {
     433            heap_init(i, j);
     434        }
     435    }
     436    printf("done.\n");
     437#endif
     438
     439    pthread_mutex_init(&print_lock, PTHREAD_PROCESS_PRIVATE);
     440    main_test_mca(num_threads, infile, outfile);
     441
     442    return 0;
     443}
     444

  • soft/giet_vm/applications/rosenfeld/src-par/mca_matrix_dist.c

    r805 r821  
    1414#include <malloc.h>
    1515
    16 #ifdef CLI
    1716#include "nrc_os_config.h"
    1817#include "nrc.h"
    19 #endif
    2018
    2119
    … …  
    6159}
    6260
     61#if TARGET_OS == GIETVM
     62// ---------------------------------------------------------------------------
     63uint32 ** remote_dist_ui32matrix(int i0, int i1, int j0, int j1, int x, int y)
     64// ---------------------------------------------------------------------------
     65{
     66    int i;
     67    int nrow = i1 - i0 + 1;
     68    int ncol = j1 - j0 + 1;
     69    uint32 ** m;
     70   
     71    // allocate pointers to rows
     72    m = (uint32 **) remote_malloc((nrow + 2) * sizeof(uint32 *), x, y);
     73    if (!m) {
     74        nrerror("allocation failure 1 in dist_ui32matrix()");
     75    }
     76    m -= i0;
     77    m += 1;
     78   
     79    // allocate rows and set pointers to them
     80    m[i0] = (uint32 *) remote_malloc((nrow * ncol + 1) * sizeof(uint32), x, y);
     81    if (!m[i0]) {
     82        nrerror("allocation failure 2 in dist_ui32matrix()");
     83    }
     84    m[i0] -= j0;
     85   
     86    for (i = i0 + 1; i <= i1; i++) {
     87        m[i] = m[i - 1] + ncol;
     88    }
     89   
     90    // make borders to point to first and last lines
     91    m[i0 - 1] = m[i0];
     92    m[i1 + 1] = m[i1];
     93   
     94    return m;
     95}
     96#endif
     97
     98
    6399
    64100// -------------------------------------------------------------------
    … …  
    66102// -------------------------------------------------------------------
    67103{
    68     free((FREE_ARG) (m[i0] + j0));
    69     free((FREE_ARG) (m + i0 - 1));
     104    free(m[i0] + j0);
     105    free(m + i0 - 1);
    70106}
    71107

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

    r805 r821  
    1212#include <string.h>
    1313#include <math.h>
    14 #include <user_barrier.h>
    15 #include <user_lock.h>
    16 
    17 #ifdef CLI
     14#include <assert.h>
     15#if PARMERGE
     16#include <pthread.h>
     17#endif
     18
    1819#include "nrc_os_config.h"
     20#include "config.h"
    1921#include "nrc.h"
     22
     23#if TARGET_OS == GIETVM
     24    #include <user_barrier.h>
     25    #include <user_lock.h>
     26    #include <giet_config.h>
     27#else
     28    #include <stdbool.h>
    2029#endif
     30
    2131
    2232#include "util.h"
    2333#include "ecc_common.h"
    24 
    2534#include "palette.h"
    2635#include "bmpNR.h"
    27 
     36#include "clock.h"
    2837#include "str_ext.h"
     38#include "ecc_features.h"
    2939
    3040// -----------
    … …  
    3444#include "mca.h"
    3545
    36 extern giet_barrier_t main_barrier;
    37 
    38 // ----------------------------------
    39 uint32 FindRoot(uint32 * T, uint32 e)
    40 // ----------------------------------
     46extern pthread_barrier_t main_barrier;
     47extern int display_features;
     48
     49CLOCK_DEC;
     50
     51
     52// -----------------------------------------
     53static uint32 FindRoot(uint32 * T, uint32 e)
     54// -----------------------------------------
    4155{
    4256    uint32 r;
    4357   
     58    assert(e != 0);
    4459    r = e;
    4560    while (T[r] < r) {
    4661        r = T[r];
    4762    }
     63    if (r == 0) {
     64        printf("e = %d\n",e);
     65        assert(0);
     66    }
    4867    return r;
    4968}
    5069
    5170
    52 // ---------------------------------------------------
    53 uint32 FindRoot_Dist(uint32 ** D, uint32 r, int shift)
    54 // ---------------------------------------------------
     71// ----------------------------------------------------------
     72static uint32 FindRoot_Dist(uint32 ** D, uint32 r, int shift)
     73// ----------------------------------------------------------
    5574{
    5675    uint32 e;
    5776    uint32 e1;
    5877    uint32 e0;
     78
     79    assert(r != 0);
    5980   
    6081    int mask = (1 << shift) - 1;
    6182   
    62     MCA_VERBOSE2(printf("FindRoot_Dist(%d) (alpha = %d) \n", r, shift));
     83    MCA_VERBOSE2(printf("%s(%d, %d) \n", __func__, r, shift));
    6384    do {
    6485        e  = r;
    … …  
    6687        e0 = r & mask;
    6788        r = D[e1][e0];
    68         MCA_VERBOSE2(printf("FindRoot: D(%d) = D[%d,%d] = %d (alpha = %d)\n", e, e1, e0, r, shift));
     89        MCA_VERBOSE2(printf("%s: D(%d) = D[%d,%d] = %d (alpha = %d)\n", __func__, e, e1, e0, r, shift));
    6990    } while (r < e);
    70     MCA_VERBOSE2(printf("FindRoot_Dist = %d \n\n", r));
     91    MCA_VERBOSE2(printf("%s = %d \n\n", __func__, r));
     92    assert(r != 0);
    7193    return r;
    7294}
    7395
    7496
    75 // -----------------------------------------
    76 void SetRoot(uint32 * T, uint32 e, uint32 r)
    77 // -----------------------------------------
    78 {
    79     while (T[e] < e) {
    80         e = T[e];
    81     }
    82     T[e] = r;
    83 }
    84 
    85 
    86 // ----------------------------------------------------------
    87 void SetRoot_Dist(uint32 ** D, uint32 e, uint32 r, int shift)
    88 // ----------------------------------------------------------
     97#if !FEATURES
     98// --------------------------------------------------------------------------------
     99static void SetRoot_Rosenfeld_Dist(uint32 ** D, uint32 root, uint32 eps, int shift)
     100// --------------------------------------------------------------------------------
    89101{
    90102    int mask = (1 << shift) - 1;
    91    
    92     uint32 e1 = e >> shift;
    93     uint32 e0 = e & mask;
    94    
    95     D[e1][e0] = r;
    96 }
    97 
    98 
    99 // --------------------------------------------
    100 uint32 Union0(uint32 * T, uint32 ei, uint32 ej)
    101 // --------------------------------------------
    102 {
    103     // version de la publication
    104     // @QM : faut-il tester le cas == 0 ici aussi ?
    105     uint32 ri, rj;
    106     ri = (ei == 0) ? 0 : FindRoot(T, ei);
    107     if (ei != ej) {
    108         rj = (ej == 0) ? 0 : FindRoot(T, ej);
    109         if (ri > rj) {
    110             ri = rj;
    111         }
    112         SetRoot(T, ej, ri);
    113     }
    114     SetRoot(T, ei, ri);
    115     return ri;
    116 }
    117 
    118 
    119 // -------------------------------------------------
    120 uint32 QuickUnion2(uint32 * T, uint32 e1, uint32 e2)
    121 // -------------------------------------------------
     103    assert(root != 0 && eps != 0);
     104   
     105    uint32 r1 = root >> shift;
     106    uint32 r0 = root & mask;
     107   
     108    D[r1][r0] = eps;
     109}
     110#endif // !FEATURES
     111
     112
     113#if FEATURES && !PARMERGE
     114// ----------------------------------------------------------------------------------------------------
     115void SetRoot_Features_Rosenfeld_Dist(uint32 ** D, uint32 root, uint32 eps, int shift, RegionStats ** F)
     116// ----------------------------------------------------------------------------------------------------
     117{
     118    assert(root != 0 && eps != 0);
     119
     120    MCA_VERBOSE2(printf("F(%d) += F(%d)\n", eps, root));
     121   
     122    int mask = (1 << shift) - 1;
     123
     124    // SetRoot_Rosenfeld_Dist
     125    uint32 r1 = root >> shift;
     126    uint32 r0 = root & mask;
     127   
     128    D[r1][r0] = eps;
     129   
     130    uint32 e1 = eps >> shift;
     131    uint32 e0 = eps & mask;
     132   
     133    // version Dist de "RegionStats_Accumulate_Stats1_From_Index"
     134   
     135    // F(eps) = F(eps) U F(root)
     136   
     137    F[e1][e0].xmin = ui16min2(F[e1][e0].xmin, F[r1][r0].xmin);
     138    F[e1][e0].xmax = ui16max2(F[e1][e0].xmax, F[r1][r0].xmax);
     139    F[e1][e0].ymin = ui16min2(F[e1][e0].ymin, F[r1][r0].ymin);
     140    F[e1][e0].ymax = ui16max2(F[e1][e0].ymax, F[r1][r0].ymax);
     141   
     142    F[e1][e0].S  += F[r1][r0].S;
     143    F[e1][e0].Sx += F[r1][r0].Sx;
     144    F[e1][e0].Sy += F[r1][r0].Sy;
     145}
     146#endif // FEATURES && !PARMERGE
     147
     148
     149#if FEATURES && PARMERGE
     150// -------------------------------------------------------------------------------------------------------------
     151bool SetRoot_Parallel_Features_Rosenfeld_Dist(uint32 ** D, uint32 root, uint32 eps, int shift, RegionStats ** F)
     152// -------------------------------------------------------------------------------------------------------------
     153{
     154    assert(root != 0 && eps != 0);
     155
     156    MCA_VERBOSE2(printf("F(%d) += F(%d)\n", eps, root));
     157   
     158    int mask = (1 << shift) - 1;
     159
     160    // SetRoot_Rosenfeld_Dist
     161    uint32 r1 = root >> shift;
     162    uint32 r0 = root & mask;
     163   
     164    uint32 e1 = eps >> shift;
     165    uint32 e0 = eps & mask;
     166
     167    // Locking towards the root (first root, then eps)
     168    pthread_spin_lock(&F[r1][r0].lock);
     169    pthread_spin_lock(&F[e1][e0].lock);
     170    // FIXME: merge these conditions later, when they both appear
     171    if (D[e1][e0] != eps) {
     172        // Someone change the root of epsilon, need to find the new root
     173        printf("race cond 1\n");
     174        pthread_spin_unlock(&F[e1][e0].lock);
     175        pthread_spin_unlock(&F[r1][r0].lock);
     176        return false;
     177    }
     178    if (D[r1][r0] != root) {
     179        // Someone change the root of epsilon, need to find the new root
     180        printf("race cond 2\n");
     181        pthread_spin_unlock(&F[e1][e0].lock);
     182        pthread_spin_unlock(&F[r1][r0].lock);
     183        return false;
     184    }
     185
     186    D[r1][r0] = eps;
     187   
     188    // F(eps) = F(eps) U F(root)
     189    F[e1][e0].xmin = ui16min2(F[e1][e0].xmin, F[r1][r0].xmin);
     190    F[e1][e0].xmax = ui16max2(F[e1][e0].xmax, F[r1][r0].xmax);
     191    F[e1][e0].ymin = ui16min2(F[e1][e0].ymin, F[r1][r0].ymin);
     192    F[e1][e0].ymax = ui16max2(F[e1][e0].ymax, F[r1][r0].ymax);
     193   
     194    F[e1][e0].S  += F[r1][r0].S;
     195    F[e1][e0].Sx += F[r1][r0].Sx;
     196    F[e1][e0].Sy += F[r1][r0].Sy;
     197
     198    pthread_spin_unlock(&F[e1][e0].lock);
     199    pthread_spin_unlock(&F[r1][r0].lock);
     200    return true;
     201}
     202#endif // FEATURES && PARMERGE
     203
     204
     205
     206#if FAST
     207// --------------------------------------------------------
     208static uint32 QuickUnion2(uint32 * T, uint32 e1, uint32 e2)
     209// --------------------------------------------------------
    122210{
    123211    // version QU de Union2
    124     uint32 r1, r2, r;
    125    
    126     r1 = (e1 == 0) ? 0 : FindRoot(T, e1);
    127     r2 = (e2 == 0) ? 0 : FindRoot(T, e2);
    128    
    129     r = ui32Min2(r1, r2);
    130    
    131     if (r1 != r) {
    132         T[r1] = r; // SetRoot
    133     }
    134     if (r2 != r) {
    135         T[r2] = r; // SetRoot
    136     }
    137    
    138     return r;
    139 }
    140 
    141 
    142 // --------------------------------------------
    143 uint32 use1_QU_Rosenfeld(uint32 e1, uint32 * T)
    144 // --------------------------------------------
    145 {
    146     return T[e1];
    147 }
    148 
    149 
    150 // -------------------------------------------------------
    151 uint32 use2_QU_Rosenfeld(uint32 e1, uint32 e2, uint32 * T)
    152 // -------------------------------------------------------
     212    uint32 r1 = FindRoot(T, e1);
     213    uint32 r2 = FindRoot(T, e2);
     214   
     215    assert(e1 != 0 && e2 != 0 && r1 != 0 && r2 != 0);
     216    uint32 eps = ui32Min2(r1, r2);
     217
     218    if (r1 > eps) {
     219        T[r1] = eps; // SetRoot sans besoin de remonter
     220    }
     221    if (r2 > eps) {
     222        T[r2] = eps; // SetRoot sans besoin de remonter
     223    }
     224    assert(e1 != 0 && e2 != 0 && r1 != 0 && r2 != 0);
     225   
     226    return eps;
     227}
     228#endif // FAST
     229
     230
     231#if FAST
     232// ---------------------------------------------------
     233static uint32 use1_QU_Rosenfeld(uint32 e1, uint32 * T)
     234// ---------------------------------------------------
     235{
     236    return FindRoot(T, e1);
     237}
     238#endif // FAST
     239
     240
     241#if FAST
     242// --------------------------------------------------------------
     243static uint32 use2_QU_Rosenfeld(uint32 e1, uint32 e2, uint32 * T)
     244// --------------------------------------------------------------
    153245{
    154246    return QuickUnion2(T, e1, e2);
    155247}
    156 
    157 
    158 // ----------------------------------------------------------------
    159 uint32 updateTable_Rosenfeld(uint32 * T, uint32 e, uint32 epsilon)
    160 // ----------------------------------------------------------------
    161 {
    162     // notations e == v, epsilon == u avec v > u (v forcement different de u)
    163     return Union0(T, e, epsilon); // original
    164 }
    165 
    166 
    167 // ----------------------------------------------------------------
    168 void vuse2_Rosenfeld(uint32 e1, uint32 e2, uint32 * T, uint32 ** D)
    169 // ----------------------------------------------------------------
    170 {
    171     uint32 e;
    172     uint32 a1;
    173     uint32 a2;
    174    
    175     a1 = (e1 == 0) ? 0 : FindRoot(T, e1);
    176     a2 = (e2 == 0) ? 0 : FindRoot(T, e2);
    177    
    178     if (a1 == a2) {
     248#endif // FAST
     249
     250
     251#if FAST && !FEATURES
     252// ---------------------------------------------------------------------------------------
     253static void vuse2_Rosenfeld_Dist(uint32 ed, uint32 el, uint32 * T, uint32 ** D, int alpha)
     254// ---------------------------------------------------------------------------------------
     255{
     256    uint32 rd = FindRoot_Dist(D, ed, alpha);
     257   
     258    uint32 rl = T[el]; // car le premier acces est local
     259    rl = FindRoot_Dist(D, rl, alpha);
     260   
     261    assert(ed != 0 && el != 0 && rd != 0 && rl != 0);
     262    if (rd == rl) {
    179263        return; // evite la backdoor
    180264    }
    181265   
    182     // forcement positifs car appel depuis optimizedBorder qui a fait un test
    183     if (a1 < a2) {
    184         e = a1;
    185         updateTable_Rosenfeld(T, a2, e);
     266    // forcement positifs car appel depuis optimizedBorder
     267    // qui a fait un test
     268    if (rd < rl) {
     269        SetRoot_Rosenfeld_Dist(D, rl, rd, alpha);
    186270    }
    187271    else {
    188         e = a2;
    189         updateTable_Rosenfeld(T, a1, e);
    190     }
    191 }
    192 
    193 
    194 // ---------------------------------------------------------------------------
    195 void vuse3_Rosenfeld(uint32 e1, uint32 e2, uint32 e3, uint32 * T, uint32 ** D)
    196 // ---------------------------------------------------------------------------
    197 {
    198     uint32 e;
    199     uint32 a1;
    200     uint32 a2;
    201     uint32 a3;
    202    
    203     a1 = (e1 == 0) ? 0 : FindRoot(T, e1);
    204     a2 = (e2 == 0) ? 0 : FindRoot(T, e2);
    205     a3 = (e3 == 0) ? 0 : FindRoot(T, e3);
    206    
    207     if (a1 == a2 && a2 == a3) {
     272        SetRoot_Rosenfeld_Dist(D, rd, rl, alpha);
     273    }
     274}
     275#endif // FAST && !FEATURES
     276
     277
     278#if FAST && !FEATURES
     279// -----------------------------------------------------------------------------------------------------
     280static void vuse3_Rosenfeld_Dist(uint32 ed1, uint32 ed2, uint32 el3, uint32 * T, uint32 ** D, int alpha)
     281// -----------------------------------------------------------------------------------------------------
     282{
     283    uint32 r1 = FindRoot_Dist(D, ed1, alpha);
     284    uint32 r2 = FindRoot_Dist(D, ed2, alpha);
     285   
     286    // QM
     287    //uint32 r3 = FindRoot(T, el3); // local - distant
     288    uint32 r3 = T[el3]; // local - distant
     289    r3 = FindRoot_Dist(D, r3, alpha);
     290
     291    assert(ed1 != 0 && ed2 != 0 && el3 != 0 && r1 != 0 && r2 != 0 && r3 != 0);
     292   
     293    if (r1 == r2 && r2 == r3) {
    208294        return;
    209295    }
    210296   
    211     e = ui32Min3(a1, a2, a3);  // forcement positifs car appel depuis optimizedBorder qui a fait un test
    212    
    213     if (a1 > e) {
    214         updateTable_Rosenfeld(T, a1, e);
    215     }
    216     a2 = T[a2];
    217     if (a2 > e) {
    218         updateTable_Rosenfeld(T, a2, e);
    219     }
    220     a3 = T[a3];
    221     if (a3 > e) {
    222         updateTable_Rosenfeld(T, a3, e);
    223     }
    224 }
    225 
    226 
    227 // ----------------------------------------------
    228 uint32 solveTable_Rosenfeld(uint32 * T, uint32 ne)
    229 // ----------------------------------------------
    230 {
    231     // equivalent a Flatten
    232     // fermeture transitive sans pack
    233     // (presence de trous dans les numeros d'etiquettes)
    234    
    235     uint32 e;
    236    
    237     for (e = 1; e <= ne; e++) {   
    238         T[e] = T[T[e]];
    239     }
    240     return ne;
    241 }
    242 
    243 
    244 // ----------------------------------------------------------------------------------
    245 uint32 optimizedAccess_DT_Rosenfeld(uint32 ** E, int i, int j, uint32 * T, uint32 ne)
    246 // ----------------------------------------------------------------------------------
    247 {
    248     // Decision Tree 8-connexe avec Quick-Union
    249     uint32 a, b, c, d, e;
    250    
    251     b = E[i - 1][j];
    252     if (b) {
    253         e = use1_QU_Rosenfeld(b, T);
     297    uint32 eps = ui32Min3(r1, r2, r3);  // forcement positifs car appel depuis optimizedBorder qui a fait un test
     298   
     299    if (r1 > eps) {
     300        SetRoot_Rosenfeld_Dist(D, r1, eps, alpha);
     301    }
     302    //r2 = T[r2]; // @QM est-ce indispensable s'il n'y a pas de features ? (cf. slow no features)
     303    // comment est-on sur que r2 (ou r3) est local ???
     304    if (r2 > eps) {
     305        SetRoot_Rosenfeld_Dist(D, r2, eps, alpha);
     306    }
     307    //r3 = T[r3];
     308    if (r3 > eps) {
     309        SetRoot_Rosenfeld_Dist(D, r3, eps, alpha);
     310    }
     311}
     312#endif // FAST && !FEATURES
     313
     314
     315#if FAST && FEATURES && !PARMERGE
     316// -----------------------------------------------------------------------------------------------------------
     317void vuse2_Features_Rosenfeld_Dist(uint32 ed, uint32 el, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
     318// -----------------------------------------------------------------------------------------------------------
     319{
     320    assert(ed != 0 && el != 0);
     321
     322    uint32 rd = FindRoot_Dist(D, ed, alpha);
     323   
     324    uint32 rl = T[el]; // car le premier acces est local
     325    assert(rl != 0);
     326    rl = FindRoot_Dist(D, rl, alpha);
     327   
     328    assert(rd != 0 && rl != 0);
     329
     330    if (rd == rl) {
     331        return; // evite la backdoor
     332    }
     333   
     334    // forcement positifs car appel depuis optimizedBorder
     335    // qui a fait un test
     336    if (rd < rl) {
     337        SetRoot_Features_Rosenfeld_Dist(D, rl, rd, alpha, F);
    254338    }
    255339    else {
    256         c = E[i - 1][j + 1];
    257         if (c) {
    258             a = E[i - 1][j - 1];
     340        SetRoot_Features_Rosenfeld_Dist(D, rd, rl, alpha, F);
     341    }
     342}
     343#endif // FAST && FEATURES && !PARMERGE
     344
     345
     346#if FAST && FEATURES && !PARMERGE
     347// -------------------------------------------------------------------------------------------------------------------------
     348void vuse3_Features_Rosenfeld_Dist(uint32 ed1, uint32 ed2, uint32 el3, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
     349// -------------------------------------------------------------------------------------------------------------------------
     350{
     351    assert(ed1 != 0 && ed2 != 0 && el3 != 0);
     352
     353    uint32 r1 = FindRoot_Dist(D, ed1, alpha);
     354    uint32 r2 = FindRoot_Dist(D, ed2, alpha);
     355   
     356    //uint32 r3 = FindRoot(T, el3); // local - distant
     357    uint32 r3 = T[el3]; // local - distant
     358    assert(r3 != 0);
     359    r3 = FindRoot_Dist(D, r3, alpha);
     360   
     361    assert(r1 != 0 && r2 != 0 && r3 != 0);
     362
     363    if (r1 == r2 && r2 == r3) {
     364        return;
     365    }
     366   
     367    uint32 eps = ui32Min3(r1, r2, r3);  // forcement positifs car appel depuis optimizedBorder qui a fait un test
     368   
     369    if (r1 > eps) {
     370        SetRoot_Features_Rosenfeld_Dist(D, r1, eps, alpha, F);
     371    }
     372    //r2 = T[r2];
     373    if (r2 > eps && r2 != r1) {
     374        SetRoot_Features_Rosenfeld_Dist(D, r2, eps, alpha, F);
     375    }
     376    //r3 = T[r3];
     377    if (r3 > eps && r3 != r2 && r3 != r1) {
     378        SetRoot_Features_Rosenfeld_Dist(D, r3, eps, alpha, F);
     379    }
     380}
     381#endif // FAST && FEATURES && !PARMERGE
     382
     383
     384#if FAST && FEATURES && PARMERGE
     385// --------------------------------------------------------------------------------------------------------------------
     386void vuse2_Parallel_Features_Rosenfeld_Dist(uint32 ed, uint32 el, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
     387// --------------------------------------------------------------------------------------------------------------------
     388{
     389    bool ok;
     390    assert(ed != 0 && el != 0);
     391    uint32 rl = T[el]; // car le premier acces est local
     392    assert(rl != 0);
     393
     394    uint32 rd;
     395   
     396    do {
     397        rd = FindRoot_Dist(D, ed, alpha); // no lock
     398        rl = FindRoot_Dist(D, rl, alpha);
     399
     400        assert(rd != 0 && rl != 0);
     401
     402        if (rd == rl) {
     403            return; // evite la backdoor
     404        }
     405
     406        // forcement positifs car appel depuis optimizedBorder
     407        // qui a fait un test
     408        if (rd < rl) {
     409            ok = SetRoot_Parallel_Features_Rosenfeld_Dist(D, rl, rd, alpha, F);
     410        }
     411        else {
     412            ok = SetRoot_Parallel_Features_Rosenfeld_Dist(D, rd, rl, alpha, F);
     413        }
     414    } while (!ok);
     415}
     416#endif // FAST && FEATURES && PARMERGE
     417
     418
     419#if FAST && FEATURES && PARMERGE
     420// ----------------------------------------------------------------------------------------------------------------------------------
     421void vuse3_Parallel_Features_Rosenfeld_Dist(uint32 ed1, uint32 ed2, uint32 el3, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
     422// ----------------------------------------------------------------------------------------------------------------------------------
     423{
     424    bool ok1, ok2, ok3;
     425    assert(ed1 != 0 && ed2 != 0 && el3 != 0);
     426
     427    uint32 r1;
     428    uint32 r2;
     429    uint32 r3 = T[el3]; // local - distant
     430    assert(r3 != 0);
     431
     432    do {
     433        r1 = FindRoot_Dist(D, ed1, alpha);
     434        r2 = FindRoot_Dist(D, ed2, alpha);
     435        r3 = FindRoot_Dist(D, r3, alpha);
     436   
     437        assert(r1 != 0 && r2 != 0 && r3 != 0);
     438
     439        if (r1 == r2 && r2 == r3) {
     440            return;
     441        }
     442   
     443        uint32 eps = ui32Min3(r1, r2, r3);  // forcement positifs car appel depuis optimizedBorder qui a fait un test
     444   
     445        ok1 = true;
     446        ok2 = true;
     447        ok3 = true;
     448        if (r1 > eps) {
     449            ok1 = SetRoot_Parallel_Features_Rosenfeld_Dist(D, r1, eps, alpha, F);
     450        }
     451        if (r2 > eps && r2 != r1) {
     452            ok2 = SetRoot_Parallel_Features_Rosenfeld_Dist(D, r2, eps, alpha, F);
     453        }
     454        if (r3 > eps && r3 != r2 && r3 != r1) {
     455            ok3 = SetRoot_Parallel_Features_Rosenfeld_Dist(D, r3, eps, alpha, F);
     456        }
     457    } while (!(ok1 && ok2 && ok3));
     458}
     459#endif // FAST && FEATURES && PARMERGE
     460
     461
     462
     463
     464#if FAST && !FEATURES
     465// ------------------------------------------------------------------------------------------------------
     466static void optimizedBorder_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha)
     467// ------------------------------------------------------------------------------------------------------
     468{
     469    uint32 a, b, c, x;
     470   
     471    x = E[i][j];
     472    if (x) {
     473        b = E[i - 1][j];
     474        if (b) {
     475            vuse2_Rosenfeld_Dist(b, x, T, D, alpha); // dist, local
     476        }
     477        else {
     478            c = E[i - 1][j + 1];
     479            if (c) {
     480                a = E[i - 1][j - 1];
     481                if (a) {
     482                    vuse3_Rosenfeld_Dist(a, c, x, T, D, alpha); // dist, local
     483                }
     484                else {
     485                    vuse2_Rosenfeld_Dist(c, x, T, D, alpha); // dist, local
     486                }
     487            }
     488            else {
     489                a = E[i - 1][j - 1];
     490                if (a) {
     491                    vuse2_Rosenfeld_Dist(a, x, T, D, alpha); // dist, local
     492                }
     493            }
     494        }
     495    }
     496}
     497#endif // FAST && !FEATURES
     498
     499
     500#if FAST && !FEATURES
     501// ---------------------------------------------------------------------------------------------------
     502static void optimizedBorderLeft_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha)
     503// ---------------------------------------------------------------------------------------------------
     504{
     505    uint32 x = E[i][j];
     506    if (x) {
     507        uint32 b = E[i - 1][j];
     508        if (b) {
     509            vuse2_Rosenfeld_Dist(b, x, T, D, alpha); // dist, local
     510        }
     511        else {
     512            uint32 c = E[i - 1][j + 1];
     513            if (c) {
     514                vuse2_Rosenfeld_Dist(c, x, T, D, alpha); // dist, local
     515            }
     516        }
     517    }
     518}
     519#endif // FAST && !FEATURES
     520
     521
     522#if FAST && !FEATURES
     523// -----------------------------------------------------------------------------------------------------------
     524static void optimizedBorderRight_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha)
     525// -----------------------------------------------------------------------------------------------------------
     526{
     527    // copie de optimizedBorder_Rosenfeld
     528    // test d'existance de ex en local local
     529
     530    uint32 b = E[i - 1][j];
     531    uint32 x = E[i][j];
     532   
     533    if (x) {
     534        if (b) {
     535            vuse2_Rosenfeld_Dist(b, x, T, D, alpha); // dist, local
     536        }
     537        else {
     538            uint32 a = E[i - 1][j - 1];
    259539            if (a) {
    260                 e = use2_QU_Rosenfeld(a, c, T);
    261             }
    262             else {
    263                 d = E[i][j - 1];
    264                 if (d) {
    265                     e = use2_QU_Rosenfeld(c, d, T);
    266                 }
    267                 else {
    268                     e = use1_QU_Rosenfeld(c, T);
    269                 }
    270             }
    271         }
    272         else {
    273             a = E[i - 1][j - 1];
    274             if (a) {
    275                 e = use1_QU_Rosenfeld(a, T);
    276             }
    277             else {
    278                 d = E[i][j - 1];
    279                 if (d) {
    280                     e = use1_QU_Rosenfeld(d, T);
    281                 }
    282                 else {
    283                     e = ++ne;
    284                 }
    285             }
    286         }
    287     }
    288     E[i][j] = e;
    289     return ne;
    290 }
    291 
    292 
    293 // ------------------------------------------------------------------------------------------
    294 void optimizedBorder_Rosenfeld(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha)
    295 // ------------------------------------------------------------------------------------------
    296 {
    297     // copie de optimizedBorder_Rosenfeld
    298     uint32 a, b, c, x;
    299    
    300     b = E[i - 1][j];
    301     x = E[i][j];
    302    
    303     if (b) {
    304         //printf("%d = %d\n", b, x);
    305         vuse2_Rosenfeld(b, x, T, D);
    306     }
    307     else {
    308         c = E[i - 1][j + 1];
    309         if (c) {
    310             a = E[i - 1][j - 1];
    311             if (a) {
    312                 //printf("%d = %d = %d\n", a, c, x);
    313                 vuse3_Rosenfeld(a, c, x, T, D);
    314             }
    315             else {
    316                 //printf("%d = %d\n", c, x);
    317                 vuse2_Rosenfeld(c, x, T, D);
    318             }
    319         }
    320         else {
    321             a = E[i - 1][j - 1];
    322             if (a) {
    323                 //printf("%d = %d\n", a, x);
    324                 vuse2_Rosenfeld(a, x, T, D);
    325             }
    326         }
    327     }
    328 }
    329 
    330 
    331 // -----------------------------------------------------------------------------------------------------------------
    332 void borderMerging_Fast_Rosenfeld_Dist(uint8 **X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha)
    333 // -----------------------------------------------------------------------------------------------------------------
    334 {
    335     for (int j = 0; j < width; j++) {
    336         if (X[i][j])  {
    337             optimizedBorder_Rosenfeld(E, i, j, T, D, alpha);
    338         }
    339     }
    340     return;
    341 }
    342 
    343 
    344 // ------------------------------------------------------------------------------------------------------------------
    345 void borderMerging_Slow_Rosenfeld_Dist(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha)
    346 // ------------------------------------------------------------------------------------------------------------------
    347 {
    348     // copie de borderMerging_Rosenfeld_UF_Fast2_8C
    349    
     540                vuse2_Rosenfeld_Dist(a, x, T, D, alpha); // dist, local
     541            }
     542        }
     543    }
     544}
     545#endif // FAST && !FEATURES
     546
     547
     548#if FAST && !FEATURES
     549// ------------------------------------------------------------------------------------------------------------------------
     550static void borderMerging_Fast_Rosenfeld_Dist(uint8 **X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha)
     551// ------------------------------------------------------------------------------------------------------------------------
     552{
     553    // Prologue
     554    optimizedBorderLeft_Rosenfeld_Dist(E, i, 0, T, D, alpha);
     555    // Boucle principale
     556    for (int j = 1; j < width - 1; j++) {
     557        optimizedBorder_Rosenfeld_Dist(E, i, j, T, D, alpha);
     558    }
     559    // Epilogue
     560    optimizedBorderRight_Rosenfeld_Dist(E, i, width - 1, T, D, alpha);
     561}
     562#endif // FAST && !FEATURES
     563
     564
     565#if SLOW && !FEATURES
     566// -------------------------------------------------------------------------------------------------------------------------
     567static void borderMerging_Slow_Rosenfeld_Dist(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha)
     568// -------------------------------------------------------------------------------------------------------------------------
     569{
    350570    int j;
    351571   
    352     uint32 e;
    353    
     572    uint32 eps;
    354573    uint32 e1, e2, e3, ex;
    355574    uint32 r1, r2, r3, rx;
    … …  
    358577    // -- prologue --
    359578    // --------------
    360     MCA_VERBOSE2(printf("[borderMerging_Slow_Rosenfeld_Dist] i = %d\n", i));
     579    MCA_VERBOSE2(printf("[%s] i = %d\n", __func__, i));
    361580   
    362581    j = 0;
    … …  
    365584    if (ex) {
    366585       
    367         MCA_VERBOSE2(printf("[borderMerging_Slow_Rosenfeld_Dist] j = %d\n", j));
     586        MCA_VERBOSE2(printf("[%s] j = %d\n", __func__, j));
    368587       
    369588        e2 = E[i - 1][j];
    370589        e3 = E[i - 1][j + 1];
    371        
    372         r2 = FindRoot_Dist(D, e2, alpha);
    373         r3 = FindRoot_Dist(D, e3, alpha);
    374         rx = FindRoot(T, ex); // we already tested that ex != 0
    375        
     590
     591        // test pour eviter acces distant
     592        r2 = e2 ? FindRoot_Dist(D, e2, alpha) : 0;
     593        r3 = e3 ? FindRoot_Dist(D, e3, alpha) : 0;
     594
     595        rx = T[ex];
     596        rx = FindRoot_Dist(D, rx, alpha);
     597 
    376598        MCA_VERBOSE2(printf("\n"));
    377599        MCA_VERBOSE2(printf("e2 = %4d -> %4d\n", e2, r2));
    … …  
    379601        MCA_VERBOSE2(printf("ex = %4d -> %4d\n", ex, rx));
    380602       
    381         e = ui32MinNonNul3(r2, r3, rx);
     603        eps = ui32MinNonNul3(r2, r3, rx);
    382604       
    383605        // Quick-Union
    384         if (r2 > e) {
    385             SetRoot_Dist(D, r2, e, alpha);
    386             MCA_VERBOSE2(printf("D[%4d] <- %d\n", r2, e));
    387         }
    388         if (r3 > e) {
    389             SetRoot_Dist(D, r3, e, alpha);
    390             MCA_VERBOSE2(printf("D[%4d] <- %d\n", r3, e));
    391         }
    392         if (rx > e) {
    393             SetRoot(T, rx, e);
    394             MCA_VERBOSE2(printf("D[%4d] <- %d\n", rx, e));
     606        if (r2 > eps) {
     607            SetRoot_Rosenfeld_Dist(D, r2, eps, alpha);
     608            MCA_VERBOSE2(printf("D[%4d] <- %d\n", r2, eps));
     609        }
     610        if (r3 > eps) {
     611            SetRoot_Rosenfeld_Dist(D, r3, eps, alpha);
     612            MCA_VERBOSE2(printf("D[%4d] <- %d\n", r3, eps));
     613        }
     614        if (rx > eps) {
     615            SetRoot_Rosenfeld_Dist(D, rx, eps, alpha);
     616            MCA_VERBOSE2(printf("D[%4d] <- %d\n", rx, eps));
    395617        }
    396618        MCA_VERBOSE2(printf("\n"));
    397         // attention SetRoot fait un while inutile
    398619    }
    399620   
    … …  
    408629        // que le cas general (pour faire un code simple)
    409630        if (ex) {
    410             MCA_VERBOSE2(printf("[borderMerging_Slow_Rosenfeld_Dist] j = %d\n", j));
     631            MCA_VERBOSE2(printf("[%s] j = %d\n", __func__, j));
    411632           
    412633            e1 = E[i - 1][j - 1];
    … …  
    414635            e3 = E[i - 1][j + 1];
    415636       
    416             r1 = FindRoot_Dist(D, e1, alpha);
    417             r2 = FindRoot_Dist(D, e2, alpha);
    418             r3 = FindRoot_Dist(D, e3, alpha);
    419             rx = FindRoot(T, ex); // we already tested that ex != 0
    420        
     637            // test pour eviter acces distant
     638            r1 = e1 ? FindRoot_Dist(D, e1, alpha) : 0;
     639            r2 = e2 ? FindRoot_Dist(D, e2, alpha) : 0;
     640            r3 = e3 ? FindRoot_Dist(D, e3, alpha) : 0;
     641
     642            rx = T[ex];
     643            rx = FindRoot_Dist(D, rx, alpha);
     644
    421645            MCA_VERBOSE2(printf("\n"));
    422646            MCA_VERBOSE2(printf("e1 = %4d -> %4d\n", e1, r1));
    … …  
    425649            MCA_VERBOSE2(printf("ex = %4d -> %4d\n", ex, rx));
    426650           
    427             e = ui32MinNonNul4(r1, r2, r3, rx);
     651            eps = ui32MinNonNul4(r1, r2, r3, rx);
    428652           
    429653            // Quick-Union
    430             if (r1 > e) {
    431                 SetRoot_Dist(D, r1, e, alpha);
    432                 MCA_VERBOSE2(printf("D[%4d] <- %d\n", r1, e));
    433             }
    434             if (r2 > e) {
    435                 SetRoot_Dist(D, r2, e, alpha);
    436                 MCA_VERBOSE2(printf("D[%4d] <- %d\n", r2, e));
    437             }
    438             if (r3 > e) {
    439                 SetRoot_Dist(D, r3, e, alpha);
    440                 MCA_VERBOSE2(printf("D[%4d] <- %d\n", r3, e));
    441             }
    442             if (rx > e) {
    443                 // @QM pourquoi pas T[e] = rx; ?
    444                 //SetRoot(T, rx, e);
    445                 T[e] = rx;
    446                 MCA_VERBOSE2(printf("D[%4d] <- %d\n", rx, e));
     654            if (r1 > eps) {
     655                SetRoot_Rosenfeld_Dist(D, r1, eps, alpha);
     656                MCA_VERBOSE2(printf("D[%4d] <- %d\n", r1, eps));
     657            }
     658            if (r2 > eps) {
     659                SetRoot_Rosenfeld_Dist(D, r2, eps, alpha);
     660                MCA_VERBOSE2(printf("D[%4d] <- %d\n", r2, eps));
     661            }
     662            if (r3 > eps) {
     663                SetRoot_Rosenfeld_Dist(D, r3, eps, alpha);
     664                MCA_VERBOSE2(printf("D[%4d] <- %d\n", r3, eps));
     665            }
     666            if (rx > eps) {
     667                SetRoot_Rosenfeld_Dist(D, rx, eps, alpha);
     668                MCA_VERBOSE2(printf("D[%4d] <- %d\n", rx, eps));
    447669            }
    448670            MCA_VERBOSE2(printf("\n"));
    … …  
    460682    if (ex) {
    461683       
    462         MCA_VERBOSE2(printf("[borderMerging_Slow_Rosenfeld_Dist] j = %d\n", j));
     684        MCA_VERBOSE2(printf("[%s] j = %d\n", __func__, j));
    463685       
    464686        e1 = E[i - 1][j - 1];
    465687        e2 = E[i - 1][j];
    466        
    467         r1 = FindRoot_Dist(D, e1, alpha);
    468         r2 = FindRoot_Dist(D, e2, alpha);
    469         rx = FindRoot(T, ex); // we already tested that ex != 0
    470        
     688
     689        // test pour eviter acces distant
     690        r1 = e1 ? FindRoot_Dist(D, e1, alpha) : 0;
     691        r2 = e2 ? FindRoot_Dist(D, e2, alpha) : 0;
     692
     693        rx = T[ex];
     694        rx = FindRoot_Dist(D, rx, alpha);
     695
    471696        MCA_VERBOSE2(printf("\n"));
    472697        MCA_VERBOSE2(printf("e1 = %4d -> %4d\n", e1, r1));
    … …  
    474699        MCA_VERBOSE2(printf("ex = %4d -> %4d\n", ex, rx));
    475700       
    476         e = ui32MinNonNul3(r1, r2, rx);
     701        eps = ui32MinNonNul3(r1, r2, rx);
    477702       
    478703        // Quick-Union
    479         if (r1 > e) {
    480             SetRoot_Dist(D, r1, e, alpha);
    481             MCA_VERBOSE2(printf("D[%4d] <- %d\n", r1, e));
    482         }
    483         if (r2 > e) {
    484             SetRoot_Dist(D, r2, e, alpha);
    485             MCA_VERBOSE2(printf("D[%4d] <- %d\n", r2, e));
    486         }
    487         if (rx > e) {
    488             SetRoot(T, rx, e);
    489             MCA_VERBOSE2(printf("D[%4d] <- %d\n", rx, e));
     704        if (r1 > eps) {
     705            SetRoot_Rosenfeld_Dist(D, r1, eps, alpha);
     706            MCA_VERBOSE2(printf("D[%4d] <- %d\n", r1, eps));
     707        }
     708        if (r2 > eps) {
     709            SetRoot_Rosenfeld_Dist(D, r2, eps, alpha);
     710            MCA_VERBOSE2(printf("D[%4d] <- %d\n", r2, eps));
     711        }
     712        if (rx > eps) {
     713            SetRoot_Rosenfeld_Dist(D, rx, eps, alpha);
     714            MCA_VERBOSE2(printf("D[%4d] <- %d\n", rx, eps));
    490715        }
    491716        MCA_VERBOSE2(printf("\n"));
    … …  
    493718    return;
    494719}
    495 
    496 
    497 // -------------------------------------------------------------------------------------------------------------
    498 void borderMerging_Rosenfeld_Dist(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha)
    499 // -------------------------------------------------------------------------------------------------------------
    500 {
     720#endif // SLOW && !FEATURES
     721
     722
     723#if SLOW && FEATURES
     724// ----------------------------------------------------------------------------------------------------------------------------------------------------
     725static void borderMerging_Slow_Features_Rosenfeld_Dist(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
     726// ----------------------------------------------------------------------------------------------------------------------------------------------------
     727{
     728    int j = 0;
     729   
     730    uint32 eps;
     731   
     732    uint32 e1, e2, e3, ex;
     733    uint32 r1, r2, r3, rx;
     734   
     735    // --------------
     736    // -- prologue --
     737    // --------------
     738    MCA_VERBOSE2(printf("[%s] i = %d\n", __func__, i));
     739   
     740    ex = E[i][j];
     741   
     742    if (ex) {
     743       
     744        MCA_VERBOSE2(printf("[%s] j = %d\n", __func__, j));
     745       
     746        e2 = E[i - 1][j];
     747        e3 = E[i - 1][j + 1];
     748       
     749        if (e2 || e3) {
     750       
     751            // test pour eviter acces distant
     752            r2 = e2 ? FindRoot_Dist(D, e2, alpha) : 0;
     753            r3 = e3 ? FindRoot_Dist(D, e3, alpha) : 0;
     754
     755            rx = T[ex];
     756            rx = FindRoot_Dist(D, rx, alpha);
     757           
     758            eps = ui32MinNonNul3(r2, r3, rx);
     759           
     760            MCA_VERBOSE2(printf("\n"));
     761            MCA_VERBOSE2(printf("e2  = %5d -> r2 = %5d\n", e2, r2));
     762            MCA_VERBOSE2(printf("e3  = %5d -> r3 = %5d\n", e3, r3));
     763            MCA_VERBOSE2(printf("ex  = %5d -> rx = %5d\n", ex, rx));
     764            MCA_VERBOSE2(printf("eps = %5d\n", eps));
     765           
     766            // Quick-Union
     767            // @QM
     768            if (r2 > eps) {
     769                SetRoot_Features_Rosenfeld_Dist(D, r2, eps, alpha, F);
     770                MCA_VERBOSE2(printf("D[%5d] <- %d\n", r2, eps));
     771            }
     772            if (r3 > 0) {
     773                r3 = FindRoot_Dist(D, r3, alpha);
     774            }
     775            // Pour le cas où r2 == r3, il ne faut pas ajouter deux fois les features
     776            //if (r3 > eps && r3 != r2) {
     777            if (r3 > eps) {
     778                SetRoot_Features_Rosenfeld_Dist(D, r3, eps, alpha, F);
     779                MCA_VERBOSE2(printf("D[%5d] <- %d\n", r3, eps));
     780            }
     781            rx = FindRoot_Dist(D, rx, alpha);
     782            //if (rx > eps && rx != r3 && rx != r2) {
     783            if (rx > eps) {
     784                SetRoot_Features_Rosenfeld_Dist(D, rx, eps, alpha, F);
     785                MCA_VERBOSE2(printf("D[%5d] <- %d\n", rx, eps));
     786            }
     787            MCA_VERBOSE2(printf("---------------------------\n"));
     788        }
     789    }
     790   
     791    // -----------------------
     792    // -- boucle principale --
     793    // -----------------------
     794   
     795    for (j = 0 + 1; j < width - 1; j++) {
     796       
     797        ex = E[i][j];
     798       
     799        if (ex) {
     800           
     801            MCA_VERBOSE2(printf("[%s] j = %d\n", __func__, j));
     802           
     803            e1 = E[i - 1][j - 1];
     804            e2 = E[i - 1][j];
     805            e3 = E[i - 1][j + 1];
     806           
     807            if (e1 || e2 || e3) {
     808                // test pour eviter un acces distant
     809                r1 = e1 ? FindRoot_Dist(D, e1, alpha) : 0;
     810                r2 = e2 ? FindRoot_Dist(D, e2, alpha) : 0;
     811                r3 = e3 ? FindRoot_Dist(D, e3, alpha) : 0;
     812
     813                rx = T[ex];
     814                rx = FindRoot_Dist(D, rx, alpha);
     815               
     816                eps = ui32MinNonNul4(r1, r2, r3, rx);
     817
     818                MCA_VERBOSE2(printf("\n"));
     819                MCA_VERBOSE2(printf("e1  = %5d -> r1 = %5d\n", e1, r1));
     820                MCA_VERBOSE2(printf("e2  = %5d -> r2 = %5d\n", e2, r2));
     821                MCA_VERBOSE2(printf("e3  = %5d -> r3 = %5d\n", e3, r3));
     822                MCA_VERBOSE2(printf("ex  = %5d -> rx = %5d\n", ex, rx));
     823                MCA_VERBOSE2(printf("eps = %5d\n", eps));
     824               
     825                // Quick-Union
     826                // @QM
     827                if (r1 > eps) {
     828                    SetRoot_Features_Rosenfeld_Dist(D, r1, eps, alpha, F);
     829                    MCA_VERBOSE2(printf("D[%5d] <- %d\n", r1, eps));
     830                }
     831                if (r2 > 0) {
     832                    r2 = FindRoot_Dist(D, r2, alpha);
     833                }
     834                //if (r2 > eps && r2 != r1) {
     835                if (r2 > eps) {
     836                    SetRoot_Features_Rosenfeld_Dist(D, r2, eps, alpha, F);
     837                    MCA_VERBOSE2(printf("D[%5d] <- %d\n", r2, eps));
     838                }
     839                if (r3 > 0) {
     840                    r3 = FindRoot_Dist(D, r3, alpha);
     841                }
     842                //if (r3 > eps && r3 != r2 && r3 != r1) {
     843                if (r3 > eps) {
     844                    SetRoot_Features_Rosenfeld_Dist(D, r3, eps, alpha, F);
     845                    MCA_VERBOSE2(printf("D[%5d] <- %d\n", r3, eps));
     846                }
     847                rx = FindRoot_Dist(D, rx, alpha);
     848                //if (rx > eps && rx != r3 && rx != r2 && rx != r1) {
     849                if (rx > eps) {
     850                    SetRoot_Features_Rosenfeld_Dist(D, rx, eps, alpha, F);
     851                    MCA_VERBOSE2(printf("D[%5d] <- %d\n", rx, eps));
     852                }
     853                MCA_VERBOSE2(puts("---------------------------\n"));
     854               
     855                // attention SetRoot fait un while inutile
     856            }
     857        }
     858    }
     859   
     860    // --------------
     861    // -- epilogue --
     862    // --------------
     863   
     864    j = width - 1;
     865    ex = E[i][j];
     866   
     867    if (ex) {
     868       
     869        MCA_VERBOSE2(printf("[%s] j = %d\n", __func__, j));
     870       
     871        e1 = E[i - 1][j - 1];
     872        e2 = E[i - 1][j];
     873       
     874        if (e1 || e2) {
     875       
     876            // test pour eviter acces distant
     877            r1 = e1 ? FindRoot_Dist(D, e1, alpha) : 0;
     878            r2 = e2 ? FindRoot_Dist(D, e2, alpha) : 0;
     879
     880            rx = T[ex];
     881            rx = FindRoot_Dist(D, rx, alpha);
     882           
     883            eps = ui32MinNonNul3(r1, r2, rx);
     884           
     885            MCA_VERBOSE2(printf("\n"));
     886            MCA_VERBOSE2(printf("e1  = %5d -> r1 = %5d\n", e1, r1));
     887            MCA_VERBOSE2(printf("e2  = %5d -> r2 = %5d\n", e2, r2));
     888            MCA_VERBOSE2(printf("ex  = %5d -> rx = %5d\n", ex, rx));
     889            MCA_VERBOSE2(printf("eps = %5d\n", eps));
     890           
     891            // Quick-Union
     892            if (r1 > eps) {
     893                SetRoot_Features_Rosenfeld_Dist(D, r1, eps, alpha, F);
     894                MCA_VERBOSE2(printf("D[%5d] <- %d\n", r1, eps));
     895            }
     896            if (r2 > 0) {
     897                r2 = FindRoot_Dist(D, r2, alpha);
     898            }
     899            //if (r2 > eps && r2 != r1) {
     900            if (r2 > eps) {
     901                SetRoot_Features_Rosenfeld_Dist(D, r2, eps, alpha, F);
     902                MCA_VERBOSE2(printf("D[%5d] <- %d\n", r2, eps));
     903            }
     904            rx = FindRoot_Dist(D, rx, alpha);
     905            //if (rx > eps && rx != r2 && rx != r1) {
     906            if (rx > eps) {
     907                SetRoot_Features_Rosenfeld_Dist(D, rx, eps, alpha, F);
     908                MCA_VERBOSE2(printf("D[%5d] <- %d\n", rx, eps));
     909            }
     910            MCA_VERBOSE2(printf("---------------------------\n"));
     911        }
     912    }
     913    return;
     914}
     915#endif // SLOW && FEATURES
     916
     917
     918#if FAST && FEATURES && !PARMERGE
     919// --------------------------------------------------------------------------------------------------------------------------
     920void optimizedBorder_Features_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
     921// --------------------------------------------------------------------------------------------------------------------------
     922{
     923    // copie de optimizedBorder_Rosenfeld
     924    uint32 a, b, c, x;
     925   
     926    x = E[i][j];
     927   
     928    if (x) {
     929        b = E[i - 1][j];
     930        if (b) {
     931            vuse2_Features_Rosenfeld_Dist(b, x, T, D, alpha, F); // dist, local
     932        }
     933        else {
     934            c = E[i - 1][j + 1];
     935            if (c) {
     936                a = E[i - 1][j - 1];
     937                if (a) {
     938                    vuse3_Features_Rosenfeld_Dist(a, c, x, T, D, alpha, F); // dist, local
     939                }
     940                else {
     941                    vuse2_Features_Rosenfeld_Dist(c, x, T, D, alpha, F); // dist, local
     942                }
     943            }
     944            else {
     945                a = E[i - 1][j - 1];
     946                if (a) {
     947                    vuse2_Features_Rosenfeld_Dist(a, x, T, D, alpha, F); // dist, local
     948                }
     949            }
     950        }
     951    }
     952}
     953#endif // FAST && FEATURES && !PARMERGE
     954
     955
     956#if FAST && FEATURES && !PARMERGE
     957// ------------------------------------------------------------------------------------------------------------------------------
     958void optimizedBorderLeft_Features_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
     959// ------------------------------------------------------------------------------------------------------------------------------
     960{
     961    uint32 x = E[i][j];
     962   
     963    if (x) {
     964        uint32 b = E[i - 1][j];
     965        if (b) {
     966            vuse2_Features_Rosenfeld_Dist(b, x, T, D, alpha, F); // dist, local
     967        }
     968        else {
     969            uint32 c = E[i - 1][j + 1];
     970            if (c) {
     971                vuse2_Features_Rosenfeld_Dist(c, x, T, D, alpha, F); // dist, local
     972            }
     973        }
     974    }
     975}
     976#endif // FAST && FEATURES && !PARMERGE
     977
     978
     979#if FAST && FEATURES && !PARMERGE
     980// -------------------------------------------------------------------------------------------------------------------------------
     981void optimizedBorderRight_Features_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
     982// -------------------------------------------------------------------------------------------------------------------------------
     983{
     984    // copie de optimizedBorder_Rosenfeld
     985    // test d'existance de ex en local local
     986   
     987    uint32 x = E[i][j];
     988   
     989    if (x) {
     990        uint32 b = E[i - 1][j];
     991        if (b) {
     992            vuse2_Features_Rosenfeld_Dist(b, x, T, D, alpha, F); // dist, local
     993        }
     994        else {
     995            uint32 a = E[i - 1][j - 1];
     996            if (a) {
     997                vuse2_Features_Rosenfeld_Dist(a, x, T, D, alpha, F); // dist, local
     998            }
     999        }
     1000    }
     1001}
     1002#endif // FAST && FEATURES && !PARMERGE
     1003
     1004
     1005#if FAST && FEATURES && PARMERGE
     1006// -----------------------------------------------------------------------------------------------------------------------------------
     1007void optimizedBorder_Parallel_Features_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
     1008// -----------------------------------------------------------------------------------------------------------------------------------
     1009{
     1010    // copie de optimizedBorder_Rosenfeld
     1011    uint32 a, b, c, x;
     1012   
     1013    x = E[i][j];
     1014   
     1015    if (x) {
     1016        b = E[i - 1][j];
     1017        if (b) {
     1018            vuse2_Parallel_Features_Rosenfeld_Dist(b, x, T, D, alpha, F); // dist, local
     1019        }
     1020        else {
     1021            c = E[i - 1][j + 1];
     1022            if (c) {
     1023                a = E[i - 1][j - 1];
     1024                if (a) {
     1025                    vuse3_Parallel_Features_Rosenfeld_Dist(a, c, x, T, D, alpha, F); // dist, local
     1026                }
     1027                else {
     1028                    vuse2_Parallel_Features_Rosenfeld_Dist(c, x, T, D, alpha, F); // dist, local
     1029                }
     1030            }
     1031            else {
     1032                a = E[i - 1][j - 1];
     1033                if (a) {
     1034                    vuse2_Parallel_Features_Rosenfeld_Dist(a, x, T, D, alpha, F); // dist, local
     1035                }
     1036            }
     1037        }
     1038    }
     1039}
     1040#endif // FAST && FEATURES && PARMERGE
     1041
     1042
     1043#if FAST && FEATURES && PARMERGE
     1044// ---------------------------------------------------------------------------------------------------------------------------------------
     1045void optimizedBorderLeft_Parallel_Features_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
     1046// ---------------------------------------------------------------------------------------------------------------------------------------
     1047{
     1048    uint32 x = E[i][j];
     1049   
     1050    if (x) {
     1051        uint32 b = E[i - 1][j];
     1052        if (b) {
     1053            vuse2_Parallel_Features_Rosenfeld_Dist(b, x, T, D, alpha, F); // dist, local
     1054        }
     1055        else {
     1056            uint32 c = E[i - 1][j + 1];
     1057            if (c) {
     1058                vuse2_Parallel_Features_Rosenfeld_Dist(c, x, T, D, alpha, F); // dist, local
     1059            }
     1060        }
     1061    }
     1062}
     1063#endif // FAST && FEATURES && PARMERGE
     1064
     1065
     1066#if FAST && FEATURES && PARMERGE
     1067// ----------------------------------------------------------------------------------------------------------------------------------------
     1068void optimizedBorderRight_Parallel_Features_Rosenfeld_Dist(uint32 ** E, int i, int j, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
     1069// ----------------------------------------------------------------------------------------------------------------------------------------
     1070{
     1071    // copie de optimizedBorder_Rosenfeld
     1072    // test d'existance de ex en local local
     1073   
     1074    uint32 x = E[i][j];
     1075   
     1076    if (x) {
     1077        uint32 b = E[i - 1][j];
     1078        if (b) {
     1079            vuse2_Parallel_Features_Rosenfeld_Dist(b, x, T, D, alpha, F); // dist, local
     1080        }
     1081        else {
     1082            uint32 a = E[i - 1][j - 1];
     1083            if (a) {
     1084                vuse2_Parallel_Features_Rosenfeld_Dist(a, x, T, D, alpha, F); // dist, local
     1085            }
     1086        }
     1087    }
     1088}
     1089#endif // FAST && FEATURES && PARMERGE
     1090
     1091
     1092#if FAST && FEATURES
     1093// ---------------------------------------------------------------------------------------------------------------------------------------------
     1094void borderMerging_Fast_Features_Rosenfeld_Dist(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
     1095// ---------------------------------------------------------------------------------------------------------------------------------------------
     1096{
     1097    MCA_VERBOSE2(printf("[%s]", __func__));
     1098   
     1099#if PARMERGE
     1100    optimizedBorderLeft_Parallel_Features_Rosenfeld_Dist(E, i, 0, T, D, alpha, F);
     1101#else
     1102    optimizedBorderLeft_Features_Rosenfeld_Dist(E, i, 0, T, D, alpha, F);
     1103#endif
     1104   
     1105    for (int j = 1; j < width - 1; j++) {
     1106#if PARMERGE
     1107        optimizedBorder_Parallel_Features_Rosenfeld_Dist(E, i, j, T, D, alpha, F);
     1108#else
     1109        optimizedBorder_Features_Rosenfeld_Dist(E, i, j, T, D, alpha, F);
     1110#endif
     1111    }
     1112   
     1113#if PARMERGE
     1114    optimizedBorderRight_Parallel_Features_Rosenfeld_Dist(E, i, width - 1, T, D, alpha, F);
     1115#else
     1116    optimizedBorderRight_Features_Rosenfeld_Dist(E, i, width - 1, T, D, alpha, F);
     1117#endif
     1118}
     1119#endif // FAST && FEATURES
     1120
     1121
     1122#if !FEATURES
     1123// --------------------------------------------------------------------------------------------------------------------
     1124static void borderMerging_Rosenfeld_Dist(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha)
     1125// --------------------------------------------------------------------------------------------------------------------
     1126{
     1127#if SLOW
    5011128    borderMerging_Slow_Rosenfeld_Dist(X, i, width, E, T, D, alpha);
    502 }
    503 
    504 
    505 // ---------------------------------------------------------------------------------------------
    506 uint32 line0Labeling_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
    507 // ---------------------------------------------------------------------------------------------
     1129#elif FAST
     1130    borderMerging_Fast_Rosenfeld_Dist(X, i, width, E, T, D, alpha);
     1131#else
     1132#error "Please define SLOW or FAST for the Rosenfeld version"
     1133#endif
     1134}
     1135#endif // !FEATURES
     1136
     1137
     1138#if FEATURES
     1139// -----------------------------------------------------------------------------------------------------------------------------------------------
     1140static void borderMerging_Features_Rosenfeld_Dist(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ** D, int alpha, RegionStats ** F)
     1141// -----------------------------------------------------------------------------------------------------------------------------------------------
     1142{
     1143#if SLOW
     1144    borderMerging_Slow_Features_Rosenfeld_Dist(X, i, width, E, T, D, alpha, F);
     1145#elif FAST
     1146    borderMerging_Fast_Features_Rosenfeld_Dist(X, i, width, E, T, D, alpha, F);
     1147#else
     1148#error "Please define SLOW or FAST for the Rosenfeld version"
     1149#endif
     1150}
     1151#endif // FEATURES
     1152
     1153
     1154// ----------------------------------------------------------------------------------------------------
     1155static uint32 line0Labeling_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
     1156// ----------------------------------------------------------------------------------------------------
    5081157{
    5091158    int j;
    … …  
    5421191
    5431192
    544 // -------------------------------------------------------------------------------------------------
    545 uint32 lineLabeling_Slow_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
    546 // -------------------------------------------------------------------------------------------------
     1193#if SLOW
     1194// --------------------------------------------------------------------------------------------------------
     1195static uint32 lineLabeling_Slow_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
     1196// --------------------------------------------------------------------------------------------------------
    5471197{
    5481198    // version lineLabeling_Rosenfeld_UF_QU_8C avec Quick-Union
    … …  
    6331283                   
    6341284                    e = ui32MinNonNul4(r1, r2, r3, r4);
    635                     giet_pthread_assert(e != 0, "e = 0\n");
    6361285                   
    6371286                    // Quick-Union
    … …  
    7071356    return ne;
    7081357}
    709 
    710 
    711 // -------------------------------------------------------------------------------------------------
    712 uint32 lineLabeling_Fast_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
    713 // -------------------------------------------------------------------------------------------------
    714 {
     1358#endif // SLOW
     1359
     1360
     1361#if FAST
     1362// ---------------------------------------------------------------------------------------------
     1363static uint32 optimizedAccessLeft_DT_Rosenfeld(uint32 ** E, int i, int j, uint32 * T, uint32 ne)
     1364// ---------------------------------------------------------------------------------------------
     1365{
     1366    // Decision Tree 8-connexe avec Quick-Union
     1367    uint32 b, c, e;
     1368   
     1369    b = E[i - 1][j];
     1370    if (b) {
     1371        e = use1_QU_Rosenfeld(b, T);
     1372    }
     1373    else {
     1374        c = E[i - 1][j + 1];
     1375        if (c) {
     1376            e = use1_QU_Rosenfeld(c, T);
     1377        }
     1378        else {
     1379            e = ++ne;
     1380        }
     1381    }
     1382    E[i][j] = e;
     1383    return ne;
     1384}
     1385#endif // FAST
     1386
     1387
     1388#if FAST
     1389// ----------------------------------------------------------------------------------------------
     1390static uint32 optimizedAccessRight_DT_Rosenfeld(uint32 ** E, int i, int j, uint32 * T, uint32 ne)
     1391// ----------------------------------------------------------------------------------------------
     1392{
     1393    // Decision Tree 8-connexe avec Quick-Union
     1394    uint32 a, b, d, e;
     1395   
     1396    b = E[i - 1][j];
     1397    if (b) {
     1398        e = use1_QU_Rosenfeld(b, T);
     1399    }
     1400    else {
     1401        a = E[i - 1][j - 1];
     1402        if (a) {
     1403            e = use1_QU_Rosenfeld(a, T);
     1404        }
     1405        else {
     1406            d = E[i][j - 1];
     1407            if (d) {
     1408                e = use1_QU_Rosenfeld(d, T);
     1409            }
     1410            else {
     1411                e = ++ne;
     1412            }
     1413        }
     1414    }
     1415    E[i][j] = e;
     1416    return ne;
     1417}
     1418#endif // FAST
     1419
     1420
     1421#if FAST
     1422// -----------------------------------------------------------------------------------------
     1423static uint32 optimizedAccess_DT_Rosenfeld(uint32 ** E, int i, int j, uint32 * T, uint32 ne)
     1424// -----------------------------------------------------------------------------------------
     1425{
     1426    // Decision Tree 8-connexe avec Quick-Union
     1427    uint32 a, b, c, d, e;
     1428   
     1429    b = E[i - 1][j];
     1430    if (b) {
     1431        e = use1_QU_Rosenfeld(b, T);
     1432    }
     1433    else {
     1434        c = E[i - 1][j + 1];
     1435        if (c) {
     1436            a = E[i - 1][j - 1];
     1437            if (a) {
     1438                e = use2_QU_Rosenfeld(a, c, T);
     1439            }
     1440            else {
     1441                d = E[i][j - 1];
     1442                if (d) {
     1443                    e = use2_QU_Rosenfeld(c, d, T);
     1444                }
     1445                else {
     1446                    e = use1_QU_Rosenfeld(c, T);
     1447                }
     1448            }
     1449        }
     1450        else {
     1451            a = E[i - 1][j - 1];
     1452            if (a) {
     1453                e = use1_QU_Rosenfeld(a, T);
     1454            }
     1455            else {
     1456                d = E[i][j - 1];
     1457                if (d) {
     1458                    e = use1_QU_Rosenfeld(d, T);
     1459                }
     1460                else {
     1461                    e = ++ne;
     1462                }
     1463            }
     1464        }
     1465    }
     1466    E[i][j] = e;
     1467    return ne;
     1468}
     1469#endif // FAST
     1470
     1471
     1472
     1473#if FAST
     1474// --------------------------------------------------------------------------------------------------------
     1475static uint32 lineLabeling_Fast_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
     1476// --------------------------------------------------------------------------------------------------------
     1477{
     1478    uint8 x;
    7151479    // avec DT et QU
    716     int j;
    717     uint8 x;
    718    
    719     for (j = 0; j < width; j++) {
    720         x = X[i][j];
     1480    // Left Border
     1481    x = X[i][0];
     1482    if (x) {
     1483        ne = optimizedAccessLeft_DT_Rosenfeld(E, i, 0, T, ne);
     1484    }
     1485    else {
     1486        E[i][0] = 0;
     1487    }
     1488    // Middle
     1489    for (int j = 1; j < width - 1; j++) {
     1490        uint8 x = X[i][j];
    7211491        if (x) {
    7221492            ne = optimizedAccess_DT_Rosenfeld(E, i, j, T, ne);
    … …  
    7261496        }
    7271497    }
     1498    // Right Border
     1499    x = X[i][width - 1];
     1500    if (x) {
     1501        ne = optimizedAccessRight_DT_Rosenfeld(E, i, width - 1, T, ne);
     1502    }
     1503    else {
     1504        E[i][width - 1] = 0;
     1505    }
    7281506    return ne;
    7291507}
    730 
    731 
    732 // --------------------------------------------------------------------------------------------
    733 uint32 lineLabeling_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
    734 // --------------------------------------------------------------------------------------------
    735 {
     1508#endif // FAST
     1509
     1510
     1511// ---------------------------------------------------------------------------------------------------
     1512static uint32 lineLabeling_Rosenfeld(uint8 ** X, int i, int width, uint32 ** E, uint32 * T, uint32 ne)
     1513// ---------------------------------------------------------------------------------------------------
     1514{
     1515#if SLOW
    7361516    return lineLabeling_Slow_Rosenfeld(X, i, width, E, T, ne);
    737     //return lineLabeling_Fast_Rosenfeld(X, i, width, E, T, ne);
    738 }
    739 
    740 
    741 // ----------------------------------------------------------------
    742 uint32 countTable_Range_Rosenfeld(uint32 * T, uint32 e0, uint32 e1)
    743 // ----------------------------------------------------------------
     1517#elif FAST
     1518    return lineLabeling_Fast_Rosenfeld(X, i, width, E, T, ne);
     1519#else
     1520#error "Please define SLOW or FAST for the Rosenfeld version"
     1521#endif
     1522}
     1523
     1524
     1525// -----------------------------------------------------------------------
     1526static uint32 countTable_Range_Rosenfeld(uint32 * T, uint32 e0, uint32 e1)
     1527// -----------------------------------------------------------------------
    7441528{
    7451529    uint32 e;
    … …  
    7551539
    7561540
    757 // --------------------------------------------------------------
    758 void solveTable_Range_Rosenfeld(uint32 * T, uint32 e0, uint32 e1)
    759 // --------------------------------------------------------------
     1541#if !FEATURES
     1542// ---------------------------------------------------------------------
     1543static void solveTable_Range_Rosenfeld(uint32 * T, uint32 e0, uint32 e1)
     1544// ---------------------------------------------------------------------
    7601545{
    7611546    uint32 e, r;
    … …  
    7661551            T[e] = r; // racine de la classe d'equivalence
    7671552        }
    768     }   
    769 }
    770 
    771 
     1553    }
     1554}
     1555#endif // !FEATURES
     1556
     1557
     1558#if FEATURES
     1559// ----------------------------------------------------------------------------------------------------------
     1560static void solveTable_solveFeatures_Range_Rosenfeld(uint32 * T, uint32 e0, uint32 e1, RegionStats * Stats)
     1561// ----------------------------------------------------------------------------------------------------------
     1562{
     1563    uint32 e, r;
     1564   
     1565    for (e = e0; e <= e1; e++) {
     1566        r = T[T[e]];
     1567        assert(r != 0);
     1568        if (r < e) {
     1569            T[e] = r; // racine de la classe d'equivalence
     1570            RegionStats_Accumulate_Stats1_From_Index(Stats, r, e);
     1571        }
     1572    }
     1573}
     1574#endif // FEATURES
     1575
     1576
     1577#if !FEATURES
    7721578// -------------------------------------
    7731579void MCA_Label_Rosenfeld_PAR1(MCA * mca)
    … …  
    7751581{
    7761582    if (mca->p == 0) {
    777         MCA_VERBOSE1(printf("------------------------------\n"));
    778         MCA_VERBOSE1(printf("-- MCA_Label_Rosenfeld_PAR1 --\n"));
    779         MCA_VERBOSE1(printf("------------------------------\n"));
    780     }
    781    
     1583        printf("*** %s ***\n", __func__);
     1584    }
     1585   
     1586    CLOCK_THREAD_START_STEP(mca->p, 0);
     1587
    7821588    int i0 = mca->i0;
    7831589    int i1 = mca->i1;
    … …  
    7951601    if (mca->p == 0) {
    7961602        set_ui32vector_j(T, e0 - 1, e1); // car e0 = 1, on a besoin que T[0] = 0 pour FindRoot
    797         // @QM : maintenant que c'est testé partout, en a-t-on encore besoin ? A priori non (a tester)
    7981603    }
    7991604    else {
    … …  
    8091614
    8101615    MCA_VERBOSE2(display_ui32matrix_positive(E, i0, i1, 0, width - 1, 5, "Ep"); printf("\n"));
    811     if (mca->p == 0) { 
     1616    if (mca->p == 0) {
    8121617        MCA_VERBOSE2(display_ui32vector_number(T, e0, ne, "%5d", "Tp_avant"));
    8131618    }
    … …  
    8151620    // fermeture transitive sans pack
    8161621    solveTable_Range_Rosenfeld(T, e0, ne);
    817     nr = countTable_Range_Rosenfeld(T, e0, ne);
    8181622    mca->ne = ne; // Plus grande etiquette de l'intervalle [e0..e1]
    8191623
     1624    MCA_VERBOSE2(nr = countTable_Range_Rosenfeld(T, e0, ne));
    8201625    MCA_VERBOSE2(printf("p = %d : e = [%d..%d] -> ne = %d -> nr = %d\n", mca->p, e0, ne, (ne - e0 + 1), nr));
    821     if (mca->p == 0) { 
     1626    if (mca->p == 0) {
    8221627        MCA_VERBOSE2(display_ui32vector_number(T, e0, ne, "%5d", "Tp_apres"));
    8231628    }
    824 }
    825 
    826 
     1629   
     1630    CLOCK_THREAD_END_STEP(mca->p, 0);
     1631}
     1632#endif // !FEATURES
     1633
     1634
     1635#if !FEATURES
    8271636// -------------------------------------
    8281637void MCA_Label_Rosenfeld_PYR2(MCA * mca)
    … …  
    8301639{
    8311640    // input
    832     int np = mca->mca->np;
    833    
    834     // variables
    835     int n = np;
    836     int nb_level = i32log2(np);
    837     if ((1 << nb_level) < np) {
    838         nb_level++; // correction pour traiter n non puissance de 2
    839     }
     1641    int p = mca->p;
     1642    int nb_level = mca->nb_level;
    8401643
    8411644    if (mca->p == 0) {
    842         MCA_VERBOSE1(printf("------------------------------\n"));
    843         MCA_VERBOSE1(printf("-- MCA_Label_Rosenfeld_PYR2 --\n"));
    844         MCA_VERBOSE1(printf("------------------------------\n"));
     1645        printf("*** %s ***\n", __func__);
    8451646    }
    8461647   
    … …  
    8621663    uint32 ** D = mca->D;
    8631664
    864     // @QM
    865     // en fait, c'est compliqué.
    866     // On pourrait optimiser en faisant faire un "break" aux procs qui n'ont plus jamais
    867     // à faire d'itération, mais le problÚme est alors qu'il faut utiliser des barriÚres avec
    868     // un nombre de procs à attendre différent à chaque fois, et qu'il faut les
    869     // initialiser => il faut précalculer toutes ces valeurs et avoir une alloc dynamique
    870     // du nombre de barriÚres.
    871     // De plus, le problÚme est décuplé si le nombre de lignes n'est pas une puissance de 2, car
    872     // dans ce cas certains threads ne doivent rien faire à une itération courante i,
    873     // mais doivent être actifs à i + 1 => encore plus dur de calculer le nombre
    874     // de threads à attendre à chaque barriÚre + surtout savoir s'il faut break ou continue
     1665    CLOCK_THREAD_START_STEP(p, 1);
     1666#if PYR_BARRIERS
     1667    // Version optimisée qui fait faire un break aux processeurs qui n'ont plus
     1668    // à faire de merge.
     1669    // Implique de pré-calculer le nombre de threads à chaque barriÚre
     1670    if (p != 0) { // thread 0 never has any merge to do
     1671        int been_active = 0;
     1672        for (int level = 0; level < nb_level; level++) {
     1673            if ((p + (1 << level)) % (1 << (level + 1)) == 0) {
     1674                borderMerging_Rosenfeld_Dist(X, i, width, E, T, D, alpha);  // en (i) et (i-1)
     1675                been_active = 1;
     1676            }
     1677            else if (been_active) {
     1678                break;
     1679            }
     1680            pthread_barrier_wait(&mca->barriers[level]);
     1681        }
     1682    }
     1683    pthread_barrier_wait(&main_barrier);
     1684#else
    8751685    for (int level = 1; level <= nb_level; level++) {
    876         if ((mca->p + (1 << (level - 1))) % (1 << level) == 0) {
     1686        if ((p + (1 << (level - 1))) % (1 << level) == 0) {
    8771687            // thread actif
    878             //MCA_VERBOSE1(printf("### level = %d - p = %d\n", level, mca->p));
    8791688            borderMerging_Rosenfeld_Dist(X, i, width, E, T, D, alpha);  // en (i) et (i-1)
    8801689        }
    881         barrier_wait(&main_barrier);
    882     }
     1690        pthread_barrier_wait(&main_barrier);
     1691    }
     1692#endif
     1693    CLOCK_THREAD_END_STEP(p, 1);
    8831694   
    8841695
    … …  
    8871698    // ---------------------------------
    8881699   
     1700    CLOCK_THREAD_START_STEP(p, 2);
    8891701    for (uint32 e = e0; e <= e1; e++) {
    8901702        uint32 r = T[e]; // acces local
    8911703        if (r < e) {
    8921704            r = FindRoot_Dist(D, e, alpha); // acces distant
    893         }
    894         T[e] = r;
    895         MCA_VERBOSE2(printf("p%d : T[%d] <- %d\n", mca->p, e, r));
    896     }
    897 }
     1705            T[e] = r; // @QM était en dehors du "if" (je pense que déjà demandé)
     1706        }
     1707        MCA_VERBOSE2(printf("p%d : T[%d] <- %d\n", p, e, r));
     1708    }
     1709    CLOCK_THREAD_END_STEP(p, 2);
     1710}
     1711#endif // !FEATURES
    8981712
    8991713
    … …  
    9041718    // input
    9051719    if (mca->p == 0) {
    906         MCA_VERBOSE1(printf("------------------------------\n"));
    907         MCA_VERBOSE1(printf("-- MCA_Label_Rosenfeld_PAR3 --\n"));
    908         MCA_VERBOSE1(printf("------------------------------\n"));
     1720        printf("*** %s ***\n", __func__);
    9091721    }
    9101722   
    … …  
    9171729    uint32 * T = mca->T;
    9181730
     1731    CLOCK_THREAD_START_STEP(mca->p, 3);
    9191732    for (int i = i0; i <= i1; i++) {
    9201733        for (int j = j0; j <= j1; j++) {
    … …  
    9251738        }
    9261739    }
    927 }
    928 
    929 
     1740    CLOCK_THREAD_END_STEP(mca->p, 3);
     1741}
     1742
     1743
     1744#if FEATURES
     1745// -----------------------------------------------------
     1746static void MCA_Label_Features_Rosenfeld_PAR1(MCA * mca)
     1747// -----------------------------------------------------
     1748{
     1749    if (mca->p == 0) {
     1750        printf("*** %s ***\n", __func__);
     1751    }
     1752   
     1753    CLOCK_THREAD_START_STEP(mca->p, 0);
     1754
     1755    int i0 = mca->i0;
     1756    int i1 = mca->i1;
     1757    int width = mca->width;
     1758
     1759    uint32 e0 = mca->e0;
     1760    uint32 e1 = mca->e1;
     1761    uint32 ne = e0 - 1;
     1762    uint32 nr = 0;
     1763
     1764    // local memory zones
     1765    uint8 **  X = mca->X;
     1766    uint32 ** E = mca->E;
     1767    uint32 *  T = mca->T;
     1768
     1769    RegionStats * stats = mca->stats;
     1770
     1771    // reset sous optimal (pour le moment = voir region32)
     1772    if (mca->p == 0) {
     1773        set_ui32vector_j(T, e0 - 1, e1); // car e0 = 1, on a besoin que T[0] = 0 pour FindRoot
     1774        zero_RegionStatsVector(stats, e0 - 1, e1);
     1775    }
     1776    else {
     1777        set_ui32vector_j(T, e0, e1);
     1778        zero_RegionStatsVector(stats, e0, e1);
     1779    }
     1780
     1781    if (mca->p == 0) {
     1782        MCA_DISPLAY2(display_ui8matrix_positive(X, i0, i1, 0, width - 1, 5, "Xp"); printf("\n"));
     1783    }
     1784
     1785    // ---------------------------- //
     1786    // -- Etiquetage d'une bande -- //
     1787    // ---------------------------- //
     1788
     1789    ne = line0Labeling_Rosenfeld(X, i0, width, E, T, ne);
     1790    lineFeaturesComputation(E, i0, width, stats);
     1791
     1792    for (int i = i0 + 1; i <= i1; i++) {
     1793        ne = lineLabeling_Rosenfeld(X, i, width, E, T, ne); // Slow or Fast
     1794        lineFeaturesComputation(E, i, width, stats);
     1795    }
     1796    mca->ne = ne; //plus grande etiquette de l'intervalle [e0..e1]
     1797
     1798    if (mca->p == 0) {
     1799        MCA_VERBOSE2(printf("ne = %d\n", ne));
     1800        MCA_DISPLAY2(display_ui32matrix_positive(E, i0, i1, 0, width - 1, 5, "Ep"); printf("\n"));
     1801        MCA_DISPLAY2(display_ui32vector_number(T, e0, ne, "%5d", "Tp_avant"));
     1802    }
     1803
     1804    // ------------------------------------------------------ //
     1805    // -- Fermeture transitive sans pack de chaque table T -- //
     1806    // ------------------------------------------------------ //
     1807
     1808    solveTable_solveFeatures_Range_Rosenfeld(T, e0, ne, stats);
     1809
     1810    if (mca->p == 0) {
     1811        MCA_VERBOSE2(nr = countTable_Range_Rosenfeld(T, e0, ne);
     1812                printf("p = %d : e = [%d..%d] -> ne = %d -> nr = %d\n", mca->p, e0, ne, (ne - e0 + 1), nr));
     1813        MCA_DISPLAY2(display_ui32vector_number(T, e0, ne, "%5d", "Tp_apres"));
     1814    }
     1815    CLOCK_THREAD_END_STEP(mca->p, 0);
     1816}
     1817#endif // FEATURES
     1818
     1819
     1820#if FEATURES && !PARMERGE
     1821// -----------------------------------------------------
     1822static void MCA_Label_Features_Rosenfeld_PYR2(MCA * mca)
     1823// -----------------------------------------------------
     1824{
     1825    int p = mca->p;
     1826    int nb_level = mca->nb_level;
     1827
     1828    if (mca->p == 0) {
     1829        printf("*** %s ***\n", __func__);
     1830    }
     1831   
     1832    // ------------------------------
     1833    // -- pyramidal border merging --
     1834    // ------------------------------
     1835   
     1836    // local variables
     1837    int i = mca->i0;
     1838    int width = mca->width;
     1839    int alpha = mca->alpha;
     1840    uint32 e0 = mca->e0;
     1841    uint32 e1 = mca->ne;
     1842
     1843    // local memory zones
     1844    uint8 **  X = mca->X;
     1845    uint32 ** E = mca->E;
     1846    uint32 *  T = mca->T;
     1847    uint32 ** D = mca->D;
     1848    RegionStats ** F = mca->F;
     1849
     1850    CLOCK_THREAD_START_STEP(p, 1);
     1851#if PYR_BARRIERS
     1852    // Version optimisée qui fait faire un break aux processeurs qui n'ont plus
     1853    // à faire de merge.
     1854    // Implique de pré-calculer le nombre de threads à chaque barriÚre
     1855    if (p != 0) { // thread 0 never has any merge to do
     1856        int been_active = 0;
     1857        for (int level = 0; level < nb_level; level++) {
     1858            if ((p + (1 << level)) % (1 << (level + 1)) == 0) {
     1859                borderMerging_Features_Rosenfeld_Dist(X, i, width, E, T, D, alpha, F);  // (i) et (i-1)
     1860                been_active = 1;
     1861            }
     1862            else if (been_active) {
     1863                break;
     1864            }
     1865            pthread_barrier_wait(&mca->barriers[level]);
     1866        }
     1867    }
     1868    pthread_barrier_wait(&main_barrier);
     1869#else
     1870    for (int level = 1; level <= nb_level; level++) {
     1871        if ((p + (1 << (level - 1))) % (1 << level) == 0) {
     1872            // thread actif
     1873            borderMerging_Features_Rosenfeld_Dist(X, i, width, E, T, D, alpha, F);  // (i) et (i-1)
     1874        }
     1875        pthread_barrier_wait(&main_barrier);
     1876    }
     1877#endif
     1878    CLOCK_THREAD_END_STEP(p, 1);
     1879
     1880
     1881    /**
     1882     * To remove?
     1883    // -- Affichage de debug
     1884    if (mca->p == 0) {
     1885        MCA_VERBOSE1(puts("-----------------------------"));
     1886        MCA_VERBOSE1(puts("[PYR2]: avant pack sequentiel"));
     1887        MCA_VERBOSE1(puts("-----------------------------"));
     1888   
     1889        for (int p = 0; p < mca->np; p++) {
     1890   
     1891            MCA* mca_par = mcas[p];
     1892            uint32 e0 = mca_par->e0;
     1893            uint32 e1 = mca_par->ne;
     1894           
     1895            uint32*  T = mca_par->T;
     1896            RegionStats* Stats = mca_par->Stats;
     1897       
     1898            RegionStats_DisplayStats_Sparse(T, e0, e1, Stats, NULL);
     1899            puts("");
     1900        }
     1901    }
     1902    */
     1903
     1904    // ---------------------------------
     1905    // -- parallel transitive closure --
     1906    // ---------------------------------
     1907    // identique a la version sans Features
     1908     
     1909    CLOCK_THREAD_START_STEP(p, 2);
     1910    for (uint32 e = e0; e <= e1; e++) {
     1911        uint32 r = T[e]; // acces local
     1912        if (r < e) {
     1913            r = FindRoot_Dist(D, e, alpha); // acces distant
     1914            T[e] = r;
     1915        }
     1916        MCA_VERBOSE2(printf("p%d : T[%d] <- %d\n", p, e, r));
     1917    }
     1918    CLOCK_THREAD_END_STEP(p, 2);
     1919
     1920    // To avoid uninitialized accesses
     1921    CLOCK_THREAD_START_STEP(p, 3);
     1922    CLOCK_THREAD_END_STEP(p, 3);
     1923}
     1924#endif // FEATURES && !PARMERGE
     1925
     1926
     1927#if FEATURES && PARMERGE
     1928// -----------------------------------------------------
     1929static void MCA_Label_Features_Rosenfeld_PAR2(MCA * mca)
     1930// -----------------------------------------------------
     1931{
     1932    int p = mca->p;
     1933    int nb_level = mca->nb_level;
     1934
     1935    if (mca->p == 0) {
     1936        printf("*** %s ***\n", __func__);
     1937    }
     1938   
     1939    // ------------------------------
     1940    // -- parallel border merging --
     1941    // ------------------------------
     1942   
     1943    // local variables
     1944    int i = mca->i0;
     1945    int width = mca->width;
     1946    int alpha = mca->alpha;
     1947    uint32 e0 = mca->e0;
     1948    uint32 e1 = mca->ne;
     1949
     1950    // local memory zones
     1951    uint8 **  X = mca->X;
     1952    uint32 ** E = mca->E;
     1953    uint32 *  T = mca->T;
     1954    uint32 ** D = mca->D;
     1955    RegionStats ** F = mca->F;
     1956
     1957    CLOCK_THREAD_START_STEP(p, 1);
     1958    if (p != 0) { // thread 0 never has any merge to do
     1959        borderMerging_Features_Rosenfeld_Dist(X, i, width, E, T, D, alpha, F);  // (i) et (i-1)
     1960    }
     1961    pthread_barrier_wait(&main_barrier);
     1962    CLOCK_THREAD_END_STEP(p, 1);
     1963
     1964
     1965    // ---------------------------------
     1966    // -- parallel transitive closure --
     1967    // ---------------------------------
     1968    // identique a la version sans Features
     1969     
     1970    CLOCK_THREAD_START_STEP(p, 2);
     1971    for (uint32 e = e0; e <= e1; e++) {
     1972        uint32 r = T[e]; // acces local
     1973        if (r < e) {
     1974            r = FindRoot_Dist(D, e, alpha); // acces distant
     1975            T[e] = r;
     1976        }
     1977        MCA_VERBOSE2(printf("p%d : T[%d] <- %d\n", p, e, r));
     1978    }
     1979    CLOCK_THREAD_END_STEP(p, 2);
     1980
     1981    // To avoid uninitialized accesses
     1982    CLOCK_THREAD_START_STEP(p, 3);
     1983    CLOCK_THREAD_END_STEP(p, 3);
     1984}
     1985#endif // FEATURES
     1986
     1987
     1988
     1989
     1990#if !FEATURES
    9301991// =============================================================
     1992#if TARGET_OS == GIETVM
    9311993__attribute__((constructor)) void MCA_Label_Rosenfeld(MCA * mca)
     1994#else
     1995void MCA_Label_Rosenfeld(MCA * mca)
     1996#endif
    9321997// =============================================================
    9331998{
     1999#if TARGET_OS == GIETVM
     2000    unsigned int x, y, lpid;
     2001    giet_proc_xyp(&x, &y, &lpid);
     2002    // Mettre à jour mca->p en fonction de x, y, lpid
     2003    // pour que les allocations faites par le main soient locales,
     2004    // i.e.
     2005    mca->p = (x * Y_SIZE + y) * NB_PROCS_MAX + lpid;
     2006    // We have :
     2007    // mca->p = 4 pour (x = 0, y = 1, lpid = 0)
     2008    // mca->p = 5 pour (x = 0, y = 1, lpid = 1)
     2009    MCA_VERBOSE2(printf("mca->p = %d pour (x = %d, y = %d, lpid = %d)\n", mca->p, x, y, lpid));
     2010#endif
     2011
     2012    CLOCK_THREAD_START(mca->p);
     2013    CLOCK_THREAD_COMPUTE_START(mca->p);
     2014
    9342015    MCA_Scatter_ImageX(mca);
    935     barrier_wait(&main_barrier);
     2016    pthread_barrier_wait(&main_barrier);
    9362017
    9372018    MCA_Label_Rosenfeld_PAR1(mca);
    938     barrier_wait(&main_barrier);
    939    
    940     //MCA_Gather_ImageL(mca);
    941     //barrier_wait(&main_barrier);
    942     //MCA_VERBOSE2(display_ui32matrix_positive(mca->E, mca->i0, mca->i1, 0, mca->width - 1, 5, "E2"));
    943     //barrier_wait(&main_barrier);
    944    
    945     //MCA_Label_Rosenfeld_SEQ2(mca);
     2019    pthread_barrier_wait(&main_barrier);
     2020   
    9462021    MCA_Label_Rosenfeld_PYR2(mca);
    947     barrier_wait(&main_barrier);
    948     //MCA_VERBOSE2(display_ui32matrix_positive(mca->E, mca->i0, mca->i1, 0, mca->width - 1, 5, "EPYR"));
    949     //barrier_wait(&main_barrier);
     2022    pthread_barrier_wait(&main_barrier);
    9502023   
    9512024    MCA_Label_Rosenfeld_PAR3(mca);
    952     barrier_wait(&main_barrier);
     2025    pthread_barrier_wait(&main_barrier);
    9532026
    9542027    MCA_Gather_ImageL(mca);
    955     barrier_wait(&main_barrier);
    956     //MCA_VERBOSE2(display_ui32matrix_positive(mca->E, mca->i0, mca->i1, 0, mca->width - 1, 5, "E3"));
    957     //barrier_wait(&main_barrier);
    958    
     2028    pthread_barrier_wait(&main_barrier);
     2029
     2030    CLOCK_THREAD_COMPUTE_END(mca->p);
     2031    CLOCK_THREAD_END(mca->p);
     2032
     2033#if TARGET_OS == GIETVM
    9592034    if (mca->p != 0) {
    960         giet_pthread_exit(NULL);
    961     }
    962 }
     2035        exit(0);
     2036    }
     2037#endif
     2038}
     2039#endif // !FEATURES
     2040
     2041
     2042#if FEATURES
     2043// ======================================================================
     2044#if TARGET_OS == GIETVM
     2045__attribute__((constructor)) void * MCA_Label_Features_Rosenfeld(void * arg)
     2046#else
     2047void * MCA_Label_Features_Rosenfeld(void * arg)
     2048#endif
     2049// ======================================================================
     2050{
     2051    MCA * mca = (MCA *) arg;
     2052#if TARGET_OS == GIETVM
     2053    unsigned int x, y, lpid;
     2054    giet_proc_xyp(&x, &y, &lpid);
     2055    // Mettre à jour mca->p en fonction de x, y, lpid
     2056    // pour que les allocations faites par le main soient locales,
     2057    // i.e.
     2058    mca->p = (x * Y_SIZE + y) * NB_PROCS_MAX + lpid;
     2059    // We have :
     2060    // mca->p = 4 pour (x = 0, y = 1, lpid = 0)
     2061    // mca->p = 5 pour (x = 0, y = 1, lpid = 1)
     2062    MCA_VERBOSE2(printf("mca->p = %d pour (x = %d, y = %d, lpid = %d)\n", mca->p, x, y, lpid));
     2063#endif
     2064
     2065    CLOCK_THREAD_START(mca->p);
     2066    CLOCK_THREAD_COMPUTE_START(mca->p);
     2067
     2068    MCA_Scatter_ImageX(mca);
     2069    pthread_barrier_wait(&main_barrier);
     2070
     2071    MCA_Label_Features_Rosenfeld_PAR1(mca);
     2072    pthread_barrier_wait(&main_barrier);
     2073   
     2074#if PARMERGE
     2075    MCA_Label_Features_Rosenfeld_PAR2(mca);
     2076#else
     2077    MCA_Label_Features_Rosenfeld_PYR2(mca);
     2078#endif
     2079    pthread_barrier_wait(&main_barrier);
     2080   
     2081    MCA_Label_Rosenfeld_PAR3(mca);
     2082    pthread_barrier_wait(&main_barrier);
     2083
     2084    MCA_Gather_ImageL(mca);
     2085    pthread_barrier_wait(&main_barrier);
     2086
     2087    CLOCK_THREAD_COMPUTE_END(mca->p);
     2088 
     2089    if (display_features) {
     2090        if (mca->p == 0) {
     2091            int i = 1;
     2092            printf("[STATS]\n");
     2093            for (int p = 0; p < mca->np; p++) {
     2094                MCA * mca_par = mca->mca->mcas[p];
     2095                uint32 e0 = mca_par->e0;
     2096                uint32 ne = mca_par->ne - mca_par->e0; // number of elements
     2097                uint32 * T = mca_par->T;
     2098                RegionStats * stats = mca_par->stats;
     2099                RegionStats_DisplayStats_Sparse(T, e0, e0 + ne, stats, NULL, &i);
     2100            }
     2101            printf("[/STATS]\n");
     2102        }
     2103    }
     2104
     2105    CLOCK_THREAD_END(mca->p);
     2106
     2107#if TARGET_OS == GIETVM
     2108    if (mca->p != 0) {
     2109        exit(0);
     2110    }
     2111#endif
     2112
     2113    return NULL;
     2114}
     2115#endif // FEATURES
     2116
    9632117
    9642118// Local Variables:
Note: See TracChangeset for help on using the changeset viewer.