source: soft/giet_vm/applications/rosenfeld/src-par/mca.c @ 821

Last change on this file since 821 was 821, checked in by meunier, 8 years ago
  • Added several versions of rosenfeld: { SLOW, FAST } x { FEATURES, NO_FEATURES }
  • Added native linux compilation support
  • Added a script to check results natively
  • Started to refactor nrc code
File size: 17.0 KB
Line 
1/* ------------- */
2/* --- mca.c --- */
3/* ------------- */
4
5/*
6 * Copyright (c) 2016 Lionel Lacassagne, LIP6, UPMC, CNRS
7 * Init  : 2016/03/03
8 * modif : 2016/03/07
9 */
10
11#include <stdio.h>
12#include <stdlib.h>
13#include <string.h>
14#include <math.h>
15#include <malloc.h>
16
17#include "nrc_os_config.h"
18#include "config.h"
19#include "nrc.h"
20
21#if TARGET_OS == GIETVM
22    #include <giet_config.h>
23#endif
24
25#include "util.h"
26#include "ecc_common.h"
27#include "ecc_features.h"
28#include "mca_matrix_dist.h"
29
30#include "palette.h"
31#include "bmpNR.h"
32#include "str_ext.h"
33
34
35// -- local --
36#include "mca.h"
37
38
39// -----------------------
40void MCA_Error(char * msg)
41// -----------------------
42{
43    printf("MCA ERROR: %s\n", msg);
44    printf("now exiting to system...\n");
45    exit(1);
46}
47
48
49// -------------------------------
50MCA * MCA_pConstructor_Empty(void)
51// -------------------------------
52{
53    MCA * mca;
54    mca = (MCA *) malloc(sizeof(MCA));
55
56    if (mca == NULL) {
57        MCA_Error("allocation failed in MCA_pConstructor_Empty");
58    }
59    return mca;
60}
61
62
63// ---------------------------------------
64void MCA_Set_ImageX(MCA * mca, uint8 ** X)
65// ---------------------------------------
66{
67    mca->X = X;
68}
69
70
71// --------------------------------------
72void MCA_Set_ImageL(MCA * mca, uint32 ** E)
73// --------------------------------------
74{
75    mca->E = E;
76}
77
78
79// ------------------------------------
80void MCA_Set_Width(MCA * mca, int width)
81// ------------------------------------
82{
83    mca->width = width;
84
85    mca->j0 = 0;
86    mca->j1 = width - 1;
87}
88
89
90// --------------------------------------
91void MCA_Set_Height(MCA * mca, int height)
92// --------------------------------------
93{
94    mca->height = height;
95 
96    mca->i0 = 0;
97    mca->i1 = height - 1;
98}
99
100
101// -----------------------------------------------
102void MCA_Set_Size(MCA * mca, int width, int height)
103// -----------------------------------------------
104{
105    MCA_Set_Width(mca, width);
106    MCA_Set_Height(mca, height);
107}
108
109
110// ----------------------------------------------------
111void MCA_Set_Dimension(MCA * mca, int width, int height)
112// ----------------------------------------------------
113{
114    MCA_Set_Width(mca, width);
115    MCA_Set_Height(mca, height);
116}
117
118
119// ------------------------------
120void MCA_Set_NP(MCA * mca, int np)
121// ------------------------------
122{
123    mca->np = np;
124}
125
126
127// ------------------------------------------------------------------
128uint32 MCA_CalcMaxLabels(int connection, uint32 height, uint32 width)
129// ------------------------------------------------------------------
130{
131    uint32 nemax = 0;
132   
133    if (connection == 4) {
134        nemax =  ((height + 1) / 2) * ((width + 1) / 2) + (height / 2) * (width / 2); // 4C
135    }
136    if (connection == 8) {
137        nemax = ((height + 1) / 2) * ((width + 1) / 2); // 8C
138    }
139    return nemax;
140}
141
142
143// ---------------------------
144void MCA_Initialize(MCA * mca)
145// ---------------------------
146{
147    // input
148    int np     = mca->np;
149    int width  = mca->width;
150    int height = mca->height;
151   
152    // variables
153    int i0_par, i1_par, i1_par_previous;
154    int j0_par, j1_par;
155    int height_par, height_mod;
156
157    int pw2;
158    int32 ne_par;          // quantite par bande
159    uint32 nemax_par;      // la puissance de 2 >=
160    uint32 e0_par, e1_par; // indice par bande [start..end]
161    int nb_level;
162   
163    MCA ** mcas;
164    MCA *  mca_par;
165   
166    printf("*** %s ***\n", __func__);
167    MCA_VERBOSE1(printf("   height = %d\n", height));
168    MCA_VERBOSE1(printf("   width  = %d\n", width));
169   
170    // array of pointers to mca workers
171    mcas = (MCA **) malloc(np * sizeof(MCA *));
172    if (mcas == NULL) {
173        MCA_Error("MCA_Initialize1");
174    }
175    mca->mcas = mcas;
176
177    // hauteur de chaque bande
178    height_par = height / np;
179    height_mod = height % np;
180
181    MCA_VERBOSE1(printf("   height_par = %d x %d + %d\n", height_par, np, height_mod));
182    MCA_VERBOSE1(printf("   ========================\n"));
183   
184    i1_par_previous = 0;
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
229        // ----------------- //
230        // -- constructor -- //
231        // ----------------- //
232        MCA_VERBOSE2(printf("-- p = %d ----------------\n", p));
233   
234        // alloc of mca workers into array of pointers
235        mca_par = MCA_pConstructor_Empty();
236        if (mca_par == NULL) {
237            MCA_Error("MCA_Initialize2\n");
238        }
239        mcas[p]      = mca_par;
240        mca_par->p   = p;
241        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
249       
250        // ------------------------------------- //
251        // -- calcul des parametres: passe #1 -- //
252        // ------------------------------------- //
253       
254        // hauteur de chaque bande
255        if (p == 0) {
256            i0_par = 0;
257        }
258        else {
259            i0_par = i1_par_previous + 1;
260        }
261        if (height_mod) {
262            i1_par = i0_par + height_par;
263            height_mod = height_mod - 1;
264        }
265        else {
266            i1_par = i0_par + height_par - 1;
267        }
268        i1_par_previous = i1_par;
269       
270        MCA_VERBOSE2(printf("i0_par = %d\n", i0_par));
271        MCA_VERBOSE2(printf("i1_par = %d\n", i1_par));
272       
273        // etiquettes
274        if (p == 0) {
275            e0_par = 1;
276            e1_par = nemax_par - 1;
277        }
278        else {
279            e0_par = p * nemax_par;
280            e1_par = e0_par + nemax_par - 1;
281        }
282   
283        MCA_VERBOSE2(printf("e0_par = %d\n", e0_par));
284        MCA_VERBOSE2(printf("e1_par = %d\n", e1_par));
285
286        mca_par->width  = width;
287        mca_par->height = height_par;
288        mca_par->i0 = i0_par;
289        mca_par->i1 = i1_par;
290        mca_par->j0 = 0;
291        mca_par->j1 = width - 1;
292        mca_par->e0 = e0_par;
293        mca_par->e1 = e1_par;
294        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
305 
306        // ---------------- //
307        // -- allocation -- //
308        // ---------------- //
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
331        mca_par->X = ui8matrix (i0_par, i1_par, 0, width - 1);
332        mca_par->E = dist_ui32matrix(i0_par, i1_par, 0, width - 1); // distributed matrix with border
333       
334        if (p == 0) {
335            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
340        }
341        else {
342            mca_par->T = ui32vector(e0_par, e1_par);
343#if FEATURES
344            mca_par->stats = RegionStatsVector(e0_par, e1_par);
345#endif
346        }
347       
348        mca_par->D = (uint32 **) vvector(0, np - 1);
349#if FEATURES
350        mca_par->F = (RegionStats **) vvector(0, np - 1);
351#endif
352#endif   
353        MCA_VERBOSE2(printf("X = %p\n", mca_par->X));
354        MCA_VERBOSE2(printf("E = %p\n", mca_par->E));
355        MCA_VERBOSE2(printf("T = %p\n", mca_par->T));
356        MCA_VERBOSE2(printf("D = %p\n", mca_par->D));
357    } // p
358   
359
360    for (int p = 0; p < np; p++) {
361        MCA * mca_par = mcas[p];
362       
363        uint32 * T = mca_par->T;
364        uint32 e0 = mca_par->e0;
365        uint32 e1 = mca_par->e1;
366   
367        MCA_VERBOSE2(printf("p = %d T[%d..%d]\n", p, e0, e1));
368        if (p == 0) {
369            set_ui32vector_j(T, e0 - 1, e1); // car e0 = 1, on a besoin que T[0] = 0 pour FindRoot
370        }
371        else {
372            set_ui32vector_j(T, e0, e1);
373        }
374        MCA_VERBOSE2(printf("\n"));
375    }
376   
377    MCA_VERBOSE2(printf("display des tables d'EQ\n"));
378    for (int p = 0; p < np; p++) {
379        MCA * mca_par = mcas[p];
380       
381        uint32 * T = mca_par->T;
382        uint32 e0 = mca_par->e0;
383        uint32 e1 = mca_par->e1;
384       
385        MCA_VERBOSE2(printf("p = %d T[%d..%d]\n", p, e0, e1));
386        if (p == 0) {
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"));
391        }
392        MCA_VERBOSE2(printf("\n"));
393    }
394    //exit(-1);
395   
396    // ------------------------------------------------------------- //
397    // -- calcul des parametres: passe #2 = parametres distribues -- //
398    // ------------------------------------------------------------- //
399   
400    // table d'indirection distribuee D
401    MCA_VERBOSE2(printf("nemax_par = %d\n", nemax_par));
402    for (int p = 0; p < np; p++) {
403        MCA * mca_p = mcas[p];
404        uint32 ** D = mca_p->D;
405        RegionStats ** F  = mca_p->F;
406       
407        for (int k = 0; k < np; k++) {
408            MCA * mca_k = mcas[k];
409            uint32 * T = mca_k->T;
410            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
415        } // k
416    } // p
417   
418    MCA_VERBOSE2(printf("table d'indirection distribuee D\n"));
419   
420    for (int p = 0; p < np; p++) {
421        MCA_VERBOSE2(printf("== p = %d ==========\n", p));
422       
423        MCA * mca_p = mcas[p];
424        uint32 ** D = mca_p->D;
425       
426        for (int k = 0; k < np; k++) {
427            MCA * mca_k = mcas[k];
428            uint32 * T = mca_k->T;
429           
430            uint32 e0 = mca_k->e0;
431            uint32 e1 = mca_k->e1;
432            MCA_VERBOSE2(display_ui32vector_number(T, e0, e0 + 9, "%5d", "T"));
433            MCA_VERBOSE2(display_ui32vector(D[k], 0, 9, "%5d", "D\n"));
434        }
435        MCA_VERBOSE2(printf("\n"));
436    }
437   
438    /**/
439    //exit(-1);
440   
441    //printf("[MCA_Initialize] positionnement des pointeurs de lignes (i0-1) et (i1+1)");
442   
443    for (int p = 0; p < np; p++) {
444   
445        //printf(" --p = %d ------------------\n", p);
446        // cas general
447       
448        if (p > 0) {
449            //printf("i0_(%d) = %d i1_{%d} = %d\n", p, mcas[p]->i0, p-1, mcas[p-1]->i1);
450            mcas[p]->E[mcas[p]->i0 - 1] = mcas[p - 1]->E[mcas[p - 1]->i1];
451           
452            /*printf("E[%2d] = E[%2d] = %p\n", mcas[p    ]->i0 - 1,
453                                             mcas[p - 1]->i1,
454                                             mcas[p - 1]->E[mcas[p - 1]->i1]);*/
455        }
456        if (p < np - 1) {
457            //printf("i1_(%d) = %d i0_{%d} = %d\n", p, mcas[p]->i1, p+1, mcas[p-1]->i0);
458            mcas[p]->E[mcas[p]->i1 + 1] = mcas[p + 1]->E[mcas[p + 1]->i0];
459           
460            /*printf("E[%2d] = E[%2d] = %p\n", mcas[p    ]->i1 + 1,
461                                             mcas[p + 1]->i0,
462                                             mcas[p + 1]->E[mcas[p + 1]->i1]);*/
463        }
464    }
465}
466
467
468// -----------------------------------
469void MCA_Display_Parameters(MCA * mca)
470// -----------------------------------
471{
472    int np = mca->np;
473   
474    MCA ** mcas = mca->mcas;
475    MCA *  mca_par;
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++) {
485        mca_par = mcas[p];
486       
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));
492    }
493}
494
495
496// -------------------------
497void MCA_Finalize(MCA * mca)
498// -------------------------
499{
500    int np = mca->np;
501   
502    MCA ** mcas = mca->mcas;
503    MCA *  mca_par;
504   
505    int i0, i1;
506    int j0, j1;
507    uint32 e0, e1;
508   
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++) {
516        mca_par = mcas[p];
517   
518        i0 = mca_par->i0;
519        i1 = mca_par->i1;
520        j0 = mca_par->j0;
521        j1 = mca_par->j1;
522        e0 = mca_par->e0;
523        e1 = mca_par->e1;
524   
525        // ---------- //
526        // -- free -- //
527        // ---------- //
528       
529        free_ui8matrix (mca_par->X, i0, i1, j0, j1);
530        free_dist_ui32matrix(mca_par->E, i0, i1, j0, j1);
531       
532        if (p == 0) {
533            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
537        }
538        else {
539            free_ui32vector(mca_par->T, e0, e1);
540#if FEATURES
541            free_RegionStatsVector(mca_par->stats, e0, e1);
542#endif
543        }
544       
545        free_vvector((void **) mca_par->D, 0, np - 1);
546#if FEATURES
547        free_vvector((void **) mca_par->F, 0, np - 1);
548#endif
549        free(mca_par);
550    }
551    free(mcas);
552    free(mca);
553}
554
555
556// -------------------------------
557void MCA_Scatter_ImageX(MCA * mca)
558// -------------------------------
559{
560    // diffusion de l'image binaire source
561   
562    int np = mca->np;
563    uint8 ** X  = mca->mca->X;
564   
565    if (mca->p == 0) { 
566        printf("*** MCA_Scatter_ImageX ***\n");
567    }
568   
569    int i0    = mca->i0;
570    int i1    = mca->i1;
571    int width = mca->width;
572
573    uint8 **  X_par = mca->X;
574    uint32 ** E_par = mca->E;
575
576    //printf("copie de [%d..%d]x[%d..%d]\n", i0, i1, 0, width - 1);
577    for (int i = i0; i <= i1; i++) {
578        for (int j = 0; j <= width - 1; j++) {
579            X_par[i][j] = X[i][j];
580            E_par[i][j] = 0; // inutile normalement car ecriture de 0
581        }
582    }
583}
584
585
586// ------------------------------
587void MCA_Gather_ImageL(MCA * mca)
588// ------------------------------
589{
590    // recuperation de l'image d'etiquettes
591    int np = mca->np;
592    uint32 ** E = mca->mca->E;
593
594    if (mca->p == 0) { 
595        printf("*** MCA_Gather_ImageL ***\n");
596    }
597
598    int i0 = mca->i0;
599    int i1 = mca->i1;
600    int width = mca->width;
601
602    uint32 ** E_par = mca->E;
603
604    for (int i = i0; i <= i1; i++) {
605        for (int j = 0; j <= width - 1; j++) {
606            E[i][j] = E_par[i][j];
607        }
608    }
609}
610
611
612// Local Variables:
613// tab-width: 4
614// c-basic-offset: 4
615// c-file-offsets:((innamespace . 0)(inline-open . 0))
616// indent-tabs-mode: nil
617// End:
618
619// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4
620
621
Note: See TracBrowser for help on using the repository browser.