source: soft/giet_vm/applications/rosenfeld/src/ecc_features.c @ 817

Last change on this file since 817 was 798, checked in by meunier, 10 years ago
  • Bug fixes in Rosenfeld
File size: 55.9 KB
Line 
1// ----------------------
2// --- ecc_features.c ---
3// ----------------------
4
5/*
6 * Copyright (c) 2012 - 2014, Lionel Lacassagne, All rights reserved
7 * University of Paris Sud, Laboratoire de Recherche en Informatique
8 */
9
10// Caracteristiques d'une region / label
11
12//#pragma message("------------------")
13//#pragma message("--- Features.c ---")
14//#pragma message("------------------")
15
16#include <stdio.h>
17#include <stddef.h>
18#include <stdlib.h>
19#include <malloc.h>
20
21
22#ifdef CLI
23#include "nrc_os_config.h"
24#include "nrc.h"
25#endif
26
27
28
29#if TARGET_OS == LINUX
30    #include <sys/types.h>
31    #include <sys/stat.h>
32    #include <fcntl.h>
33#endif
34
35
36#ifdef OPENMP
37#include <omp.h>
38#endif
39
40#include "ecc_features.h"
41//#include "label.h"
42
43// ------------------------------------------------------------
44void RegionStats_Constructor(RegionStats **Stats, uint32 nemax)
45// ------------------------------------------------------------
46{
47    *Stats = RegionStats_pConstructor(nemax);
48}
49// ------------------------------------------------
50RegionStats* RegionStats_pConstructor(uint32 nemax)
51// ------------------------------------------------
52{
53    RegionStats *Stats;
54   
55    Stats = (RegionStats*) malloc((nemax)*sizeof(RegionStats));
56    if(Stats == NULL) nrerror("allocation failed in RegionStats_pConstructor");
57   
58    RegionStats_Clear(Stats, nemax);
59   
60    return Stats;
61}
62// -----------------------------------------------------------
63void RegionStats_Destructor(RegionStats **Stats, uint32 nemax)
64// -----------------------------------------------------------
65{
66    RegionStats_pDestructor(*Stats, nemax);
67}
68// -----------------------------------------------------------
69void RegionStats_pDestructor(RegionStats *Stats, uint32 nemax)
70// -----------------------------------------------------------
71{
72    RegionStats_Clear(Stats, nemax);
73    free(Stats);
74}
75// -----------------------------------------------------
76void RegionStats_Clear(RegionStats *Stats, uint32 nemax)
77// -----------------------------------------------------
78{
79    int i;
80    for (i = 0; i < (int) nemax; i++) {
81        Stats[i].xmin = 65535;
82        Stats[i].xmax = 0;
83        Stats[i].ymin = 65535;
84        Stats[i].ymax = 0;
85       
86        Stats[i].S = 0;
87       
88        Stats[i].x = 0;
89        Stats[i].y = 0;
90       
91        Stats[i].Sx = 0;
92        Stats[i].Sy = 0;
93       
94#ifdef REGION_STATS2
95        Stats[i].Sx2 = 0;
96        Stats[i].Sxy = 0;
97        Stats[i].Sy2 = 0;
98#endif
99    }
100}
101// ----------------------------------------
102void RegionStats_Clear1(RegionStats *stats)
103// ----------------------------------------
104{
105    stats->xmin = 0;
106    stats->xmax = 0;
107    stats->ymin = 0;
108    stats->ymax = 0;
109   
110    stats->S = 0;
111   
112    stats->x = 0;
113    stats->y = 0;
114   
115    stats->Sx = 0;
116    stats->Sy = 0;
117#ifdef REGION_STATS2
118    stats->Sx2 = 0;
119    stats->Sxy = 0;
120    stats->Sy2 = 0;
121#endif
122}
123// ------------------------------------------
124int RegionStats_Create_File(char *filename)
125// ------------------------------------------
126{
127    int fd;
128   
129    fd = open(filename, O_CREAT | O_TRUNC);
130    if (fd < 0) {
131        printf("RegionStats_Open_File : can't create file %s\n", filename);
132        giet_pthread_exit("");
133    }
134    return fd;
135}
136// ----------------------------------------
137int RegionStats_Open_File(char *filename)
138// ----------------------------------------
139{
140    int fd;
141   
142    fd = open(filename, O_RDONLY);
143    if (fd < 0) {
144        printf("RegionStats_Open_File : can't open file %s\n", filename);
145        giet_pthread_exit("");
146    }
147    return fd;
148}
149// ---------------------------------
150void RegionStats_Close_File(int fd)
151// ---------------------------------
152{
153    close(fd);
154}
155// ---------------------------------
156int RegionStats_Read_Header(int fd)
157// ---------------------------------
158{
159    int ne = 0;
160    // @QM giet
161    fscanf(fd, "%d", &ne);
162    return ne;
163}
164// -------------------------------------------
165void RegionStats_Write_Header(int ne, int fd)
166// -------------------------------------------
167{
168    fprintf(fd, "%d\n", ne);
169}
170// --------------------------------------------------------------
171void RegionStats_Read_Stats1(int fd, int ne, RegionStats *Stats)
172// --------------------------------------------------------------
173{
174    int i;
175   
176    for (i = 1; i <= ne; i++) {
177        fscanf(fd, "%d%d%d%d%d%d%d%d\n",
178               &t,
179               &(Stats[i].xmin),
180               &(Stats[i].xmax),
181               &(Stats[i].ymin),
182               &(Stats[i].ymax),
183               &(Stats[i].S),
184               &(Stats[i].Sx),
185               &(Stats[i].Sy));
186       
187        Stats[i].x = Stats[i].Sx / Stats[i].S;
188        Stats[i].y = Stats[i].Sy / Stats[i].S;
189       
190    }
191}
192// --------------------------------------------------------------
193void RegionStats_Read_Stats2(int fd, int ne, RegionStats *Stats)
194// --------------------------------------------------------------
195{
196#ifdef REGION_STATS2
197    int i, t;
198   
199    int32 Sx2, Sxy, Sy2;
200   
201    for (i = 1; i <= ne; i++) {
202        fscanf(f, "%d%d%d%d%d%d%d%d%d%d%d\n",
203               &t,
204               &(Stats[i].xmin),
205               &(Stats[i].xmax),
206               &(Stats[i].ymin),
207               &(Stats[i].ymax),
208               &(Stats[i].S),
209               &(Stats[i].Sx),
210               &(Stats[i].Sy),
211               &(Sx2),
212               &(Sxy),
213               &(Sy2));
214       
215        Stats[i].Sx2 = Sx2;
216        Stats[i].Sxy = Sxy;
217        Stats[i].Sy2 = Sy2;
218       
219        Stats[i].x = Stats[i].Sx / Stats[i].S;
220        Stats[i].y = Stats[i].Sy / Stats[i].S;
221    }
222#else
223    nrerror("RegionStats_Read_Stats2 REGION_STAT2 not defined");
224#endif
225}
226// ---------------------------------------------------------------
227void RegionStats_Write_Stats1(RegionStats *Stats, int ne, int fd)
228// ---------------------------------------------------------------
229{
230    int i;
231   
232    for (i = 1; i <= ne; i++) {
233        fprintf(fd, "%4d %5d %5d %5d %5d %7d %8d %8d\n",
234                i,
235                Stats[i].xmin,
236                Stats[i].xmax,
237                Stats[i].ymin,
238                Stats[i].ymax,
239               
240                Stats[i].S,
241                Stats[i].Sx,
242                Stats[i].Sy);
243    }
244}
245// --------------------------------------------------------------------------------------------------
246void RegionStats_Write_Stats1_Sparse(RegionStats *Stats, uint32 *EQ, uint32 ne0, uint32 ne1, int fd)
247// --------------------------------------------------------------------------------------------------
248{
249    uint32 e;
250   
251    for (e = ne0; e <= ne1; e++) {
252        if ((e == EQ[e]) && (Stats[e].S > 0)) {
253            fprintf(fd, "%4d %5d %5d %5d %5d %7d %8d %8d\n",
254                    e,
255                    Stats[e].xmin,
256                    Stats[e].xmax,
257                    Stats[e].ymin,
258                    Stats[e].ymax,
259                   
260                    Stats[e].S,
261                    Stats[e].Sx,
262                    Stats[e].Sy);
263        }
264    }
265}
266// ---------------------------------------------------------------
267void RegionStats_Write_Stats2(RegionStats *Stats, int ne, int fd)
268// ---------------------------------------------------------------
269{
270#ifdef REGION_STATS2
271    int i;
272    for (i = 1; i <= ne; i++) {
273        fprintf(fd, "%4d %4d %4d %4d %6d %8d %8d %8d %8d %8d\n",
274                i,
275                Stats[i].xmin,
276                Stats[i].xmax,
277                Stats[i].ymin,
278                Stats[i].ymax,
279               
280                Stats[i].S,
281                Stats[i].Sx,
282                Stats[i].Sy,
283               
284                (int32) Stats[i].Sx2,
285                (int32) Stats[i].Sxy,
286                (int32) Stats[i].Sy2);
287    }
288#else
289    nrerror("RegionStats_Write_Stats2: REGION_STATS2 not defined");
290#endif
291}
292// -----------------------------------------------------------------
293void RegionStats_Write_pStats1(RegionStats **Stats, int ne, int fd)
294// -----------------------------------------------------------------
295{
296    int i;
297   
298    for (i = 1; i <= ne; i++) {
299        fprintf(fd, "%4d %5d %5d %5d %5d %7d %8d %8d\n",
300                i,
301                Stats[i]->xmin,
302                Stats[i]->xmax,
303                Stats[i]->ymin,
304                Stats[i]->ymax,
305                Stats[i]->S,
306                Stats[i]->Sx,
307                Stats[i]->Sy);
308    }
309}
310// -----------------------------------------------------------------
311void RegionStats_Write_pStats2(RegionStats **Stats, int ne, int fd)
312// -----------------------------------------------------------------
313{
314#ifdef REGION_STATS2
315    int i;
316    for (i = 1; i <= ne; i++) {
317        fprintf(fd, "%3d %4d %4d %4d %4d %6d %8d %8d %8d %8d %8d\n",
318                i,
319                Stats[i]->xmin,
320                Stats[i]->xmax,
321                Stats[i]->ymin,
322                Stats[i]->ymax,
323                Stats[i]->S,
324                Stats[i]->Sx,
325                Stats[i]->Sy,
326               
327                (int32) Stats[i]->Sx2,
328                (int32) Stats[i]->Sxy,
329                (int32) Stats[i]->Sy2);
330    }
331#else
332    nrerror("RegionStats_Write_Stats2: REGION_STATS2 not defined");
333#endif
334}
335// -----------------------------------------------------------------------
336void RegionStats_Load_Stats1(char *filename, int *ne, RegionStats **Stats)
337// -----------------------------------------------------------------------
338{
339    int fd;
340   
341    fd = RegionStats_Open_File(filename);
342   
343    *ne = RegionStats_Read_Header(fd);
344   
345    RegionStats_Constructor(Stats, *ne);
346    RegionStats_Read_Stats1(fd, *ne, *Stats);
347    RegionStats_Close_File(fd);
348}
349// -----------------------------------------------------------------------
350void RegionStats_Load_Stats2(char *filename, int *ne, RegionStats **Stats)
351// -----------------------------------------------------------------------
352{
353#ifdef REGION_STATS2
354    int fd;
355   
356    fd = RegionStats_Open_File(filename);
357   
358    *ne = RegionStats_Read_Header(fd);
359   
360    RegionStats_Constructor(Stats, *ne);
361    RegionStats_Read_Stats2(fd, *ne, *Stats);
362    RegionStats_Close_File(fd);
363#else
364    nrerror("RegionStats_Load_Stats2 : REGION_STATS2 not defined");
365#endif
366}
367// -----------------------------------------------------------------------
368void RegionStats_MLoad_Stats1(char *filename, int *ne, RegionStats *Stats)
369// -----------------------------------------------------------------------
370{
371    int fd;
372   
373    fd = RegionStats_Open_File(filename);
374   
375    *ne = RegionStats_Read_Header(fd);
376    RegionStats_Read_Stats1(fd, *ne, Stats);
377    RegionStats_Close_File(fd);
378}
379// -----------------------------------------------------------------------
380void RegionStats_MLoad_Stats2(char *filename, int *ne, RegionStats *Stats)
381// -----------------------------------------------------------------------
382{
383#ifdef REGION_STATS2
384    int fd
385   
386    fd = RegionStats_Open_File(filename);
387   
388    *ne = RegionStats_Read_Header(fd);
389    RegionStats_Read_Stats2(fd, *ne, Stats);
390    RegionStats_Close_File(fd);
391#else
392    nrerror("RegionStats_MLoad_Stats2 : REGION_STATS2 not defined");
393#endif
394}
395// ---------------------------------------------------------------------
396void RegionStats_Save_Stats1(RegionStats *Stats, int ne, char *filename)
397// ---------------------------------------------------------------------
398{
399    int fd;
400   
401    fd = RegionStats_Create_File(filename);
402   
403    RegionStats_Write_Header(ne, fd);
404    RegionStats_Write_Stats1(Stats, ne, fd);
405    RegionStats_Close_File(fd);
406}
407// ---------------------------------------------------------------------
408void RegionStats_Save_Stats2(RegionStats *Stats, int ne, char *filename)
409// ---------------------------------------------------------------------
410{
411#ifdef REGION_STATS2
412    int fd;
413   
414    fd = RegionStats_Create_File(filename);
415   
416    RegionStats_Write_Header(ne, fd);
417    RegionStats_Write_Stats2(Stats, ne, fd);
418    RegionStats_Close_File(fd);
419#else
420    nrerror("RegionStats_Save_Stats2 : REGION_STATS2 not defined");
421#endif
422}
423// -----------------------------------------------------------------------
424void RegionStats_Save_pStats1(RegionStats **Stats, int ne, char *filename)
425// -----------------------------------------------------------------------
426{
427    int fd;
428   
429    fd = RegionStats_Create_File(filename);
430   
431    RegionStats_Write_Header(ne, fd);
432    RegionStats_Write_pStats1(Stats, ne, fd);
433    RegionStats_Close_File(fd);
434}
435// -----------------------------------------------------------------------
436void RegionStats_Save_pStats2(RegionStats **Stats, int ne, char *filename)
437// -----------------------------------------------------------------------
438{
439#ifdef REGION_STATS2
440    int fd;
441   
442    fd = RegionStats_Create_File(filename);
443   
444    RegionStats_Write_Header(ne, fd);
445    RegionStats_Write_pStats2(Stats, ne, fd);
446    RegionStats_Close_File(fd);
447#else
448    nrerror("RegionStats_Save_Stats2 : REGION_STATS2 not defined");
449#endif
450}
451// --------------------------------------------------------------------
452void RegionStats_Display_Stats1(RegionStats *Stats, int ne, char *name)
453// --------------------------------------------------------------------
454{
455    int i;
456   
457    if (name != NULL) {
458        printf("%s : %d\n", name, ne);
459    }
460    else {
461        printf("RegionStats : %d\n", ne);
462    }
463    for (i = 1; i <= ne; i++) {
464        printf("#%3d: %4d %4d %4d %4d %6d %8d %8d\n",
465               i,
466               Stats[i].xmin,
467               Stats[i].xmax,
468               Stats[i].ymin,
469               Stats[i].ymax,
470               Stats[i].S,
471               Stats[i].Sx,
472               Stats[i].Sy);
473    }
474}
475// --------------------------------------------------------------------
476void RegionStats_Display_Stats2(RegionStats *Stats, int ne, char *name)
477// --------------------------------------------------------------------
478{
479#ifdef REGION_STATS2
480    int i;
481   
482    if (name != NULL) {
483        printf("%s : %d\n", name, ne);
484    }
485    else {
486        printf("RegionStats : %d\n", ne);
487    }
488   
489    for (i = 1; i <= ne; i++) {
490        printf("#%3d: %4d %4d %4d %4d %6d %8d %8d %8d %8d %8d\n",
491               i,
492               Stats[i].xmin,
493               Stats[i].xmax,
494               Stats[i].ymin,
495               Stats[i].ymax,
496               Stats[i].S,
497               Stats[i].Sx,
498               Stats[i].Sy,
499               
500               (int32) Stats[i].Sx2,
501               (int32) Stats[i].Sxy,
502               (int32) Stats[i].Sy2);
503       
504    }
505#else
506    nrerror("RegionStats_Display_Stats2 : REGION_STATS2 not defined");
507#endif
508}
509// ----------------------------------------------------------------------
510void RegionStats_Display_pStats1(RegionStats **Stats, int ne, char *name)
511// ----------------------------------------------------------------------
512{
513    int i;
514   
515    if (name != NULL) {
516        printf("%s : %d\n", name, ne);
517    }
518    else {
519        printf("RegionStats : %d\n", ne);
520    }
521    for (i = 1; i <= ne; i++) {
522        printf("#%3d: %4d %4d %4d %4d %6d %8d %8d\n",
523               i,
524               Stats[i]->xmin,
525               Stats[i]->xmax,
526               Stats[i]->ymin,
527               Stats[i]->ymax,
528               Stats[i]->S,
529               Stats[i]->Sx,
530               Stats[i]->Sy);
531    }
532}
533// ----------------------------------------------------------------------
534void RegionStats_Display_pStats2(RegionStats **Stats, int ne, char *name)
535// ----------------------------------------------------------------------
536{
537#ifdef REGION_STATS2
538    int i;
539   
540    if (name != NULL) {
541        printf("%s : %d\n", name, ne);
542    }
543    else {
544        printf("RegionStats : %d\n", ne);
545    }
546    for (i = 1; i <= ne; i++) {
547        printf("#%3d: %4d %4d %4d %4d %6d %8d %8d %8d %8d %8d\n",
548               i,
549               Stats[i]->xmin,
550               Stats[i]->xmax,
551               Stats[i]->ymin,
552               Stats[i]->ymax,
553               Stats[i]->S,
554               Stats[i]->Sx,
555               Stats[i]->Sy,
556               
557               (int32) Stats[i]->Sx2,
558               (int32) Stats[i]->Sxy,
559               (int32) Stats[i]->Sy2);
560       
561    }
562#else
563    nrerror("RegionStats_Display_Stats2 : REGION_STATS2 not defined");
564#endif
565}
566// ---------------------------------------------------------------------------------------------
567void RegionStats_SetRectangle(RegionStats *Stats, int e, int ymin, int ymax, int xmin, int xmax)
568// ---------------------------------------------------------------------------------------------
569{
570    Stats[e].ymin = ymin;
571    Stats[e].xmin = xmin;
572    Stats[e].ymax = ymax;
573    Stats[e].xmax = xmax;
574}
575// -------------------------------------------------------
576void RegionStats_Copy1(RegionStats *src, RegionStats *dst)
577// -------------------------------------------------------
578{
579    dst->xmin = src->xmin;
580    dst->xmax = src->xmax;
581    dst->ymin = src->ymin;
582    dst->ymax = src->ymax;
583   
584    dst->S    = src->S;
585   
586    dst->x    = src->x;
587    dst->y    = src->y;
588   
589    dst->Sx   = src->Sx;
590    dst->Sy   = src->Sy;
591   
592#ifdef REGION_STATS2
593    dst->Sx2  = src->Sx2;
594    dst->Sxy  = src->Sxy;
595    dst->Sy2  = src->Sy2;
596    dst->teta = src->teta;
597#endif
598   
599#ifdef REGION_STATS3
600    dst->Sx3  = src->Sx3;
601    dst->Sx2y = src->Sx2y;
602    dst->Sxy2 = src->Sxy2;
603    dst->Sy3  = src->Sy3;
604   
605    dst->Mx2  = src->Mx2;
606    dst->Mxy  = src->Mxy;
607    dst->My2  = src->My2;
608   
609    dst->Mx3  = src->Mx3;
610    dst->Mx2y = src->Mx2y;
611    dst->Mxy2 = src->Mxy2;
612    dst->My3  = src->My3;
613#endif
614}
615// ===============================
616// === nouvelles versions 2009 ===
617// ===============================
618
619// -------------------------------------------
620RegionStats* RegionStatsVector(int i0, int i1)
621// -------------------------------------------
622// allocate a float vector with subscript range v[i0..i1]
623{
624    RegionStats * v;
625   
626    v = (RegionStats *) malloc((size_t) ((i1 - i0 + 1 + NR_END) * sizeof(RegionStats)));
627    if (!v) nrerror("allocation failure in RegionStatsVector()");
628    if (!v) return NULL;
629    return v - i0 + NR_END;
630}
631// ----------------------------------------------------------
632RegionStats* RegionStatsVector0(int i0, int i1)
633// ----------------------------------------------------------
634// allocate a float vector with subscript range v[i0..i1]
635{
636    RegionStats *v;
637   
638    v = (RegionStats *) calloc((size_t) (i1 - i0 + 1 + NR_END), sizeof(RegionStats));
639    if (!v) nrerror("allocation failure in RegionStatsVector0()");
640    if (!v) return NULL;
641    return v - i0 + NR_END;
642}
643// ----------------------------------------------------------------------
644void free_RegionStatsVector(RegionStats *v, int i0, int i1)
645// ----------------------------------------------------------
646// free a RegionStats vector allocated with vector()
647{
648    free((FREE_ARG) (v + i0 - NR_END));
649}
650// ------------------------------------------------------------
651RegionStats** RegionStatsMatrix(int i0, int i1, int j0, int j1)
652// ------------------------------------------------------------
653
654// allocate a RegionStats matrix with subscript range m[nrl..nrh][ncl..nch]
655{
656    long i;
657    long nrow = i1 - i0 + 1;
658    long ncol = j1 - j0 + 1;
659    RegionStats **m;
660   
661    // allocate pointers to rows
662    m = (RegionStats **) malloc((size_t) ((nrow + NR_END) * sizeof(RegionStats *)));
663    if (!m) nrerror("allocation failure 1 in RegionStatsMatrix()");
664    m += NR_END;
665    m -= i0;
666   
667    // allocate rows and set pointers to them
668    m[i0] = (RegionStats *) malloc((size_t) ((nrow * ncol + NR_END) * sizeof(RegionStats)));
669    if (!m[i0]) nrerror("allocation failure 2 in RegionStatsMatrix()");
670    m[i0] += NR_END;
671    m[i0] -= j0;
672   
673    for(i = i0 + 1; i <= i1; i++) {
674        m[i] = m[i - 1] + ncol;
675    }
676   
677    // return pointer to array of pointers to rows
678    return m;
679}
680// -------------------------------------------------------------
681RegionStats** RegionStatsMatrix0(int i0, int i1, int j0, int j1)
682// -------------------------------------------------------------
683
684// allocate a float matrix with subscript range m[nrl..nrh][ncl..nch]
685{
686    long i, nrow=i1-i0+1,ncol=j1-j0+1;
687    RegionStats **m;
688   
689    // allocate pointers to rows
690    m=(RegionStats**) malloc((size_t)((nrow+NR_END)*sizeof(RegionStats*)));
691    if (!m) nrerror("allocation failure 1 in RegionStatsMatrix()");
692    m += NR_END;
693    m -= i0;
694   
695    // allocate rows and set pointers to them
696    m[i0]=(RegionStats*) calloc((size_t)(nrow*ncol+NR_END), sizeof(RegionStats));
697    if (!m[i0]) nrerror("allocation failure 2 in RegionStatsMatrix()");
698    m[i0] += NR_END;
699    m[i0] -= j0;
700   
701    for(i=i0+1;i<=i1;i++) m[i]=m[i-1]+ncol;
702   
703    // return pointer to array of pointers to rows
704    return m;
705}
706// -------------------------------------------------------------------------
707void free_RegionStatsMatrix(RegionStats **m, int i0, int i1, int j0, int j1)
708// -------------------------------------------------------------------------
709{
710    free((FREE_ARG) (m[i0]+j0-NR_END));
711    free((FREE_ARG) (m+i0-NR_END));
712}
713// ----------------------------------
714void zero_RegionStats(RegionStats *x)
715// ----------------------------------
716{
717    x->xmin = 32767;
718    x->xmax = 0;
719    x->ymin = 32767;
720    x->ymax = 0;
721   
722    x->S  = 0;
723    x->Sx = 0;
724    x->Sy = 0;
725   
726#ifdef REGION_STATS2
727    x->Sx2 = 0;
728    x->Sxy = 0;
729    x->Sy2 = 0;
730#endif
731}
732// --------------------------------------------------------
733void zero_RegionStatsVector(RegionStats *v, int i0, int i1)
734// --------------------------------------------------------
735{
736    int i;
737    for(i=i0; i<=i1; i++) {
738        zero_RegionStats(&v[i]);
739    }
740}
741// -------------------------------------------------------------------------
742void zero_RegionStatsMatrix(RegionStats **m, int i0, int i1, int j0, int j1)
743// -------------------------------------------------------------------------
744{
745    int i, j;
746    for(i=i0; i<=i1; i++) {
747        for(j=j0; j<=j1; j++) {         
748            zero_RegionStats(&(m[i][j]));
749        }
750    }
751}
752// -------------------------------------------------
753void display_RegionStats(RegionStats *x, char *name)
754// -------------------------------------------------
755{
756    if(name != NULL) printf("%s : \n", name);
757   
758#ifndef REGION_STATS2
759    printf("%4d %4d %4d %4d %6d %8d %8d\n",
760           x->xmin,
761           x->xmax,
762           x->ymin,
763           x->ymax,
764           
765           x->S,
766           x->Sx,
767           x->Sy);
768#else
769    printf("%4d %4d %4d %4d %6d %8d %8d %8d %8d %8d\n",
770           x->xmin,
771           x->xmax,
772           x->ymin,
773           x->ymax,
774           
775           x->S,
776           x->Sx,
777           x->Sy,
778           
779           x->Sx2,
780           x->Sxy,
781           x->Sy2);
782#endif 
783}
784// -----------------------------------------------------------------------
785void display_RegionStatsVector(RegionStats *v, int i0, int i1, char *name)
786// -----------------------------------------------------------------------
787{
788    int i;
789   
790    if(name != NULL) printf("%s : [%d..%d]\n", name, i0, i1); else printf("RegionStats : [%d..%d]\n", i0, i1);
791    for(i=i0; i<=i1; i++) {
792        printf("#%3d: ", i);
793        display_RegionStats(&(v[i]), NULL);
794        //puts("");
795    }
796}
797// ----------------------------------------------------------------------------------------
798void display_RegionStatsMatrix(RegionStats **m, int i0, int i1, int j0, int j1, char *name)
799// ----------------------------------------------------------------------------------------
800{
801    int i, j;
802   
803    if (name != NULL) printf("%s : [%d..%d][%d..%d]\n", name, i0, i1, j0, j1); else printf("RegionStats : [%d..%d][%d..%d]\n", i0, i1, j0, j1);
804    for (i = i0; i <= i1; i++) {
805        for (j = j0; j <= j1; j++) {
806            printf("#%3d: ", i);
807            display_RegionStats(&(m[i][j]), NULL);
808        }
809    }
810}
811// ----------------------------------------------
812void save_RegionStats(RegionStats *x, char *name)
813// ----------------------------------------------
814{
815    int fd = -1;
816   
817    if (name == NULL) {
818        return;
819    }
820    // assume name != NULL if single element
821    // assume name == NULL if vector or matrix
822    fd = RegionStats_Create_File(name);
823    if (fd <= 0) {
824        printf("*** Erreur : ouverture du fichier %s dans %s\n", name, __func__);
825    }
826    fprintf(fd, "%s: ", name);
827   
828#ifndef REGION_STATS2
829    fprintf(fd, "%4d %4d %4d %4d %6d %8d %8d\n",
830           
831            x->xmin,
832            x->xmax,
833            x->ymin,
834            x->ymax,
835           
836            x->S,
837            x->Sx,
838            x->Sy);
839#else
840    fprintf(fd, "%4d %4d %4d %4d %6d %8d %8d %8d %8d %8d\n",
841           
842            x->xmin,
843            x->xmax,
844            x->ymin,
845            x->ymax,
846           
847            x->S,
848            x->Sx,
849            x->Sy,
850           
851            x->Sx2,
852            x->Sxy,
853            x->Sy2);
854#endif 
855   
856    if (name) {
857        RegionStats_Close_File(fd);
858    }   
859}
860// --------------------------------------------------------------------
861void save_RegionStatsVector(RegionStats *v, int i0, int i1, char *name)
862// --------------------------------------------------------------------
863{
864    int i;
865    int fd;
866   
867    if (name == NULL) {
868        name = "RegionStatsVector";
869    }
870    fd = RegionStats_Create_File(name);
871   
872    fprintf(fd, "%s : [%d..%d]\n", name, i0, i1);
873   
874    for (i = i0; i <= i1; i++) {
875        printf("#%3d: ", i);
876        save_RegionStats(&v[i], NULL);
877        printf("");
878    }
879    RegionStats_Close_File(fd);
880}
881// -------------------------------------------------------------------------------------
882void save_RegionStatsMatrix(RegionStats **m, int i0, int i1, int j0, int j1, char *name)
883// -------------------------------------------------------------------------------------
884{
885    int i, j;
886    int fd;
887   
888    if (name == NULL) {
889        name = "RegionStatsMatrix";
890    }
891    fd = RegionStats_Create_File(name);
892   
893    fprintf(fd, "%s : [%d..%d]\n", name, i0, i1);
894   
895    for (i = i0; i <= i1; i++) {
896        for (j = j0; j <= j1; j++) {
897            fprintf(fd, "#%3d: ", i);
898            save_RegionStats(&m[i][j], NULL);
899        }
900        printf("");
901    }
902    RegionStats_Close_File(fd);
903}
904// ------------------------------------------------------------------------------
905void RegionStats_Calc1_Features_1Pass(RegionStats *Stats, uint32 e, int i, int j)
906// ------------------------------------------------------------------------------
907{
908    // calcul sur 1 point et non sur toute l'image
909        // Rectangle
910       
911    if (i < Stats[e].ymin) Stats[e].ymin = i;
912        if (i > Stats[e].ymax) Stats[e].ymax = i;
913   
914    if (j < Stats[e].xmin) Stats[e].xmin = j;
915        if (j > Stats[e].xmax) Stats[e].xmax = j;       
916   
917        // Moment1
918        Stats[e].S  += 1;
919        Stats[e].Sx += j;
920        Stats[e].Sy += i;
921       
922        return;
923}
924// --------------------------------
925// --- fonctions de 2013 ----------
926// --------------------------------
927// -------------------------------------------------------------------------------------------
928void RegionStats_Calc_Rectangle_Moment1(uint32 **E, int height, int width, RegionStats *Stats)
929// -------------------------------------------------------------------------------------------
930{
931    int i, j;
932    uint32 x, y;
933    uint32 e;
934   
935    for (i = 0; i < height; i++) {
936        for (j = 0; j < width; j++) {
937           
938            e = E[i][j];
939            if (e) {
940               
941                x = j;
942                y = i;
943               
944                if (i<Stats[e].ymin) Stats[e].ymin = y;
945                if (i>Stats[e].ymax) Stats[e].ymax = y;
946               
947                if (j<Stats[e].xmin) Stats[e].xmin = x;
948                if (j>Stats[e].xmax) Stats[e].xmax = x;
949               
950                Stats[e].S  += 1;
951                Stats[e].Sx += x;
952                Stats[e].Sy += y;
953            }
954        }
955    }
956}
957// -----------------------------------------------------------------------------------------------------------------------------
958void RegionStats_calc_Status(RegionStats *Stats, uint32 ne, uint32 min_height, uint32 min_width, uint32 min_area, uint8 *status)
959// -----------------------------------------------------------------------------------------------------------------------------
960{
961    uint16 xmin, xmax, ymin, ymax, xsize, ysize;
962    uint32 size;
963    uint32 e;
964   
965    for (e = 1; e < ne; e++) {
966       
967        ymin = Stats[e].ymin;
968        ymax = Stats[e].ymax;
969       
970        xmin = Stats[e].xmin;
971        xmax = Stats[e].xmax;
972       
973        ysize = ymax - ymin + 1;
974        xsize = xmax - xmin + 1;
975        size = xsize * ysize;
976       
977        if ((size > min_area) && (xsize >= min_width) && (ysize >= min_height)) {
978            status[e] = 1;
979        }
980        else {
981            status[e] = 0;
982        }
983    }
984}
985// --------------------------------------------------------------------------
986uint32 RegionStats_UpdateEQ_with_Status(uint8 *status, uint32 ne, uint32 *EQ)
987// --------------------------------------------------------------------------
988{
989    uint32 e;
990    uint32 na = 0;
991    for (e = 1; e < ne; e++) {
992        if (status[e]) {
993            EQ[e] = ++na;
994        }
995        else {
996            EQ[e] = 0;
997        }
998    }
999    return na;
1000}
1001// ----------------------------------------------------------------------------
1002void RegionStats_UpdateStats_with_EQ(uint32 *EQ, uint32 ne, RegionStats *Stats)
1003// ----------------------------------------------------------------------------
1004{
1005    uint32 e, a;
1006    for (e = 1; e < ne; e++) {
1007        a = EQ[e];
1008        if (a != e) {
1009            // copy
1010            RegionStats_Copy1(&Stats[e], &Stats[a]);
1011        }
1012        else {
1013            // do nothing
1014        }
1015    }
1016}
1017// ---------------------------------------------------------------------------
1018void featuresComputation(uint32 **E, int height,int width, RegionStats *Stats)
1019// ---------------------------------------------------------------------------
1020{
1021    //uint32 nemax = height * width /2;   
1022    RegionStats_Calc_Rectangle_Moment1(E, height, width, Stats);   
1023}
1024// ------------------------------------------------------------------------------
1025void pointFeaturesComputation_Dummy(uint32 **E, int i, int j, RegionStats *Stats)
1026// ------------------------------------------------------------------------------
1027{
1028    // pour pointeur de fonction
1029    return;
1030}
1031// -------------------------------------------------------------------------
1032void pointFeaturesComputation( uint32 **E, int i, int j, RegionStats *Stats)
1033// -------------------------------------------------------------------------
1034{
1035    uint32 x, y;
1036    uint32 e;
1037   
1038    e = E[i][j];
1039    if (e) {
1040       
1041        x = j;
1042        y = i;
1043       
1044        if (i<Stats[e].ymin) Stats[e].ymin = y;
1045        if (i>Stats[e].ymax) Stats[e].ymax = y;
1046       
1047        if (j<Stats[e].xmin) Stats[e].xmin = x;
1048        if (j>Stats[e].xmax) Stats[e].xmax = x;
1049       
1050        Stats[e].S  += 1;
1051        Stats[e].Sx += x;
1052        Stats[e].Sy += y;
1053    }
1054}
1055// ----------------------------------------------------------------------------------
1056void lineFeaturesComputation_Dummy( uint32 **E, int i, int width, RegionStats *Stats)
1057// ----------------------------------------------------------------------------------
1058{
1059    // pour pointeur de fonction
1060    return;
1061}
1062// ----------------------------------------------------------------------------
1063void lineFeaturesComputation( uint32 **E, int i, int width, RegionStats *Stats)
1064// ----------------------------------------------------------------------------
1065{
1066    // line RegionStats_Calc_Rectangle_Moment1
1067    int j;
1068   
1069    uint32 x, y;
1070    uint32 e;
1071   
1072    for (j = 0; j < width; j++) {
1073       
1074        e = E[i][j];
1075        if (e) {
1076           
1077            x = j;
1078            y = i;
1079           
1080            if (i<Stats[e].ymin) Stats[e].ymin = y;
1081            if (i>Stats[e].ymax) Stats[e].ymax = y;
1082           
1083            if (j<Stats[e].xmin) Stats[e].xmin = x;
1084            if (j>Stats[e].xmax) Stats[e].xmax = x;
1085           
1086            Stats[e].S  += 1;
1087            Stats[e].Sx += x;
1088            Stats[e].Sy += y;
1089        }
1090    }
1091}
1092// ------------------------------------------------------------------------------------------
1093void bandFeaturesComputation_Dummy(uint32 **E, int i0, int i1, int width, RegionStats *Stats)
1094// ------------------------------------------------------------------------------------------
1095{
1096    return;
1097}
1098// ------------------------------------------------------------------------------------
1099void bandFeaturesComputation(uint32 **E, int i0, int i1, int width, RegionStats *Stats)
1100// ------------------------------------------------------------------------------------
1101{
1102    int i;
1103    for (i = i0; i <= i1; i++) {
1104        lineFeaturesComputation(E, i, width, Stats);
1105    }
1106}
1107// ---------------------------------------------------------------------------------------
1108void imageFeaturesComputation_Dummy(uint32 **E, int height, int width, RegionStats *Stats)
1109// ---------------------------------------------------------------------------------------
1110{
1111    // pour pointeur de fonction
1112    return;
1113}
1114// ---------------------------------------------------------------------------------
1115void imageFeaturesComputation(uint32 **E, int height, int width, RegionStats *Stats)
1116// ---------------------------------------------------------------------------------
1117{
1118    // image RegionStats_Calc_Rectangle_Moment1
1119    int i;
1120    for (i = 0; i < height; i++) {
1121        lineFeaturesComputation(E, i, width, Stats);
1122    }
1123}
1124// ---------------------------------------
1125// --- Fonctions 2014 --------------------
1126// ---------------------------------------
1127
1128// --------------------------------------------------------------------------------------
1129void RegionStats_Copy_Stats1_From_Index(RegionStats *Stats, int dst_index, int src_index)
1130// --------------------------------------------------------------------------------------                                 
1131{
1132    // R[dst] = R[src]
1133    RegionStats_Copy1(&Stats[src_index], &Stats[dst_index]);
1134}
1135// --------------------------------------------------------------------------------------------
1136void RegionStats_Accumulate_Stats1_From_Index(RegionStats *Stats, int dst_index, int src_index)
1137// --------------------------------------------------------------------------------------------                                 
1138{
1139    // R[dst] += R[src]
1140    Stats[dst_index].xmin = ui16min2(Stats[dst_index].xmin, Stats[src_index].xmin);
1141    Stats[dst_index].xmax = ui16max2(Stats[dst_index].xmax, Stats[src_index].xmax);
1142    Stats[dst_index].ymin = ui16min2(Stats[dst_index].ymin, Stats[src_index].ymin);
1143    Stats[dst_index].ymax = ui16max2(Stats[dst_index].ymax, Stats[src_index].ymax);
1144   
1145    Stats[dst_index].S  += Stats[src_index].S;
1146    Stats[dst_index].Sx += Stats[src_index].Sx;
1147    Stats[dst_index].Sy += Stats[src_index].Sy;   
1148}
1149// -----------------------------------------------------------------------------------------------------
1150void RegionStats_DisplayStats_Sparse(uint32 *EQ, uint32 ne0, uint32 ne1, RegionStats *Stats, char *name)
1151// -----------------------------------------------------------------------------------------------------
1152{
1153    // n'affiche que les racines.
1154    // ne pas utiliser apres un pack, car le test n'a plus de sens
1155    uint32 e;
1156    uint32 na; // compteur
1157   
1158    if (name) printf(name);
1159   
1160    na = RegionStats_Count_Roots_Sparse(Stats, EQ, ne0, ne1);
1161    printf("%d\n", na);
1162   
1163    for (e = ne0; e <= ne1; e++) {
1164        if ((e == EQ[e]) && (Stats[e].S > 0)) {
1165            printf("%5d ", e);
1166            display_RegionStats(&Stats[e], NULL);
1167        }
1168    }
1169}
1170// ----------------------------------------------------------------------------------------------------
1171void RegionStats_DisplayStats_Range(uint32 *EQ, uint32 ne0, uint32 ne1, RegionStats *Stats, char *name)
1172// ----------------------------------------------------------------------------------------------------
1173{
1174    // affichage dense (apres un pack)
1175
1176    uint32 e;
1177   
1178    if (name) printf(name);
1179   
1180    for (e = ne0; e <= ne1; e++) {
1181        printf("%5d ", e);
1182        display_RegionStats(&Stats[e], NULL);
1183    }
1184}
1185// --------------------------------------------------------------------------------------------------------
1186void RegionStats_Save_Stats1_Sparse(RegionStats *Stats, uint32 *EQ, uint32 ne0, uint32 ne1, char *filename)
1187// --------------------------------------------------------------------------------------------------------
1188{
1189    int fd;
1190    uint32 na = 0;
1191   
1192    fd = RegionStats_Create_File(filename);
1193    na = RegionStats_Count_Roots_Sparse(Stats, EQ, ne0, ne1);
1194   
1195    RegionStats_Write_Header(na, fd);
1196    RegionStats_Write_Stats1_Sparse(Stats, EQ, ne0, ne1, fd);
1197    RegionStats_Close_File(fd);
1198}
1199// ------------------------------------------------------------------------------------------
1200uint32 RegionStats_Count_Roots_Sparse(RegionStats *Stats, uint32 *EQ, uint32 ne0, uint32 ne1)
1201// ------------------------------------------------------------------------------------------
1202{
1203    uint32 e, c = 0; // compteur
1204   
1205    for (e = ne0; e <= ne1; e++) {
1206        if ((e == EQ[e]) && (Stats[e].S > 0)) {
1207            c++;
1208        }
1209    }
1210    return c;
1211}
1212// ---------------------------------------------------------------------------------
1213uint32 RegionStats_Count_Roots_Sparse1(RegionStats *Stats, uint32 *EQ, uint32 nemax)
1214// ---------------------------------------------------------------------------------
1215{
1216    return RegionStats_Count_Roots_Sparse(Stats, EQ, 1, nemax);
1217}
1218// -------------------------------------------------------------------------------------------
1219uint32 RegionStats_Count_Labels_Sparse(RegionStats *Stats, uint32 *EQ, uint32 ne0, uint32 ne1)
1220// -------------------------------------------------------------------------------------------
1221{
1222    uint32 e, c = 0; // compteur
1223   
1224    for (e = ne0; e <= ne1; e++) {
1225        if (Stats[e].S > 0) {
1226            c++;
1227        }
1228    }
1229    return c;
1230}
1231// ----------------------------*-----------------------------------------------------
1232uint32 RegionStats_Count_Labels_Sparse1(RegionStats *Stats, uint32 *EQ, uint32 nemax)
1233// ---------------------------*------------------------------------------------------
1234{
1235    return RegionStats_Count_Labels_Sparse(Stats, EQ, 1, nemax);
1236}
1237// ---------------------------------------------------------------------
1238void copy_features_ui32matrix(RegionStats *Stats, uint32 ne, uint32 **m)
1239// ---------------------------------------------------------------------
1240{
1241    int i;
1242    for (i = 0; i <= (int) ne; i++) {
1243       
1244        m[i][0] = i;
1245       
1246        // 16 bits: clean but requires 2 sorts
1247        //m[i][1] = Stats[i].xmin;
1248        //m[i][2] = Stats[i].xmax;
1249        //m[i][3] = Stats[i].ymin;
1250        //m[i][4] = Stats[i].ymax;
1251       
1252        // 32 bits: dirty, but requires only 1 sort
1253        m[i][1] = (Stats[i].ymin << 16) | Stats[i].ymax;
1254        m[i][2] = (Stats[i].xmin << 16) | Stats[i].xmax;
1255
1256        // 32 bits
1257        m[i][3] = Stats[i].S;
1258        m[i][4] = Stats[i].Sx;
1259        m[i][5] = Stats[i].Sy;
1260    }
1261}
1262// ---------------------------------------------------------------------
1263void copy_ui32matrix_features(uint32 **m, uint32 ne, RegionStats *Stats)
1264// ---------------------------------------------------------------------
1265{
1266    int i;
1267    for (i = 0; i <= (int) ne; i++) {
1268       
1269        Stats[i].xmin = m[i][2] >> 16;
1270        Stats[i].xmax = m[i][2] & 0xffff;
1271        Stats[i].ymin = m[i][1] >> 16;
1272        Stats[i].ymax = m[i][1] & 0xffff;
1273       
1274        Stats[i].S  = m[i][3];
1275        Stats[i].Sx = m[i][4];
1276        Stats[i].Sy = m[i][5];
1277    }
1278}
1279// ---------------------------------------------------------------------------
1280void sortv_ui32matrix_col(uint32 **m, int i0, int i1, int j0, int j1, int col)
1281// ---------------------------------------------------------------------------
1282{
1283    // nrsort2 for NRC2
1284    //sortv_ui32matrix_selection_min(m, i0, i1, j0, j1,  col);
1285   
1286    long nrl = i0;
1287    long nrh = i1;
1288    long nc  = col;
1289   
1290    /*
1291     * sort an matrix of int, with the selection algorithm.
1292     * the key is in column nc
1293     * the sort is performed, by doing a purmutation on the lines,
1294     * instead of copying the lines.
1295     */
1296        int i, j;
1297       
1298    uint32 x, min, pos;
1299        uint32 * ptr;
1300       
1301        for (i = nrl; i < nrh; i++) {
1302                min = m[i][nc];
1303                pos = i;
1304                for (j = i + 1; j <= nrh; j++) {
1305                        x = m[j][nc];
1306                        if (x < min) {
1307                                min = x;
1308                                pos = j;
1309                        }
1310                } // j
1311               
1312                // permutation des pointeurs de ligne de la matrice
1313                ptr    = m[i];
1314                m[i]   = m[pos];
1315                m[pos] = ptr;
1316               
1317        } // i
1318   
1319}
1320// ------------------------------------------------------------
1321void RegionStats_SortFeatures(RegionStats *Stats, uint32 nemax)
1322// ------------------------------------------------------------
1323{
1324    uint32 ** m = NULL;
1325   
1326    m = ui32matrix(0, nemax, 0, 7);
1327   
1328    copy_features_ui32matrix(Stats, nemax, m);
1329    sortv_ui32matrix_col(m, 0, nemax, 0, 5, 1);
1330    copy_ui32matrix_features(m, nemax, Stats);
1331}
1332// --------------------------------------------------------------------------------------
1333void imageFeaturesComputation_omp0(uint32 **E, int height, int width, RegionStats *Stats)
1334// --------------------------------------------------------------------------------------
1335{
1336    // version OpenMP 2.0 fausse (sans serialisation de la section critique)
1337    // pour evaluer l'impact de la synchro
1338    int i, j;
1339   
1340    uint32 x, y;
1341    uint32 e;
1342   
1343#ifdef OPENMP
1344#pragma omp parallel for private(height, width, i, j, x, y, e) shared (E, Stats)
1345#endif
1346    for (i = 0; i < height; i++) {
1347        for (j = 0; j < width; j++) {
1348           
1349            e = E[i][j];
1350            if (e) {
1351               
1352                x = j;
1353                y = i;
1354               
1355                // min max reduction
1356                if (y < Stats[e].ymin) Stats[e].ymin = y;
1357                if (y > Stats[e].ymax) Stats[e].ymax = y;
1358               
1359                if (x < Stats[e].xmin) Stats[e].xmin = x;
1360                if (x > Stats[e].xmax) Stats[e].xmax = x;
1361               
1362                // + reduction
1363                Stats[e].S  += 1;
1364                Stats[e].Sx += x;
1365                Stats[e].Sy += y;
1366            }
1367        }
1368    }
1369}
1370// --------------------------------------------------------------------------------------
1371void imageFeaturesComputation_omp2(uint32** E, int height, int width, RegionStats* Stats)
1372// --------------------------------------------------------------------------------------
1373{
1374    // version OpenMP 2.0 classique avec "critical"
1375   
1376    int i, j;
1377   
1378    uint32 x, y;
1379    uint32 e;
1380   
1381   
1382    #ifdef OPENMP
1383    //#pragma omp parallel for private(height, width, i, j, x, y, e) shared(E, Stats)
1384    #pragma omp parallel for shared(E, Stats) private(height, width, i, j, x, y, e) schedule(dynamic)
1385    //#pragma omp for private (j)
1386    #endif   
1387    for (i = 0; i < height; i++) {
1388        //printf("i = %d\n", i);
1389        //printf("omp_get_num_threads = %d\n", omp_get_num_threads());
1390       
1391        for (j = 0; j < width; j++) {
1392           
1393            e = E[i][j];
1394            if (e) {
1395               
1396                x = j;
1397                y = i;
1398               
1399                #ifdef OPENMP
1400                #pragma omp critical
1401                #endif
1402                {
1403                    // min max reduction
1404                    if (y < Stats[e].ymin) Stats[e].ymin = y;
1405                    if (y > Stats[e].ymax) Stats[e].ymax = y;
1406                   
1407                    if (x < Stats[e].xmin) Stats[e].xmin = x;
1408                    if (x > Stats[e].xmax) Stats[e].xmax = x;
1409                   
1410                    // + reduction
1411                    Stats[e].S  += 1;
1412                    Stats[e].Sx += x;
1413                    Stats[e].Sy += y;
1414                } // omp critical
1415            } // if e
1416        } // j
1417    } // i
1418}
1419// --------------------------------------------------------------------------------------
1420void imageFeaturesComputation_omp3(uint32** E, int height, int width, RegionStats* Stats)
1421// --------------------------------------------------------------------------------------
1422{
1423    // version OpenMP 2.0 classique avec "critical" (from Laurent Cabaret with optimal use of critical and atomic)
1424
1425  int i, j;
1426   
1427    uint32 x, y;
1428    uint32 e;
1429
1430
1431    #ifdef OPENMP
1432    //#pragma omp parallel for private(height, width, i, j, x, y, e) shared(E, Stats)
1433    #pragma omp parallel for shared(E, Stats) private(height, width, i, j, x, y, e) schedule(dynamic)
1434    #endif   
1435    for (i = 0; i < height; i++) {
1436
1437        for (j = 0; j < width; j++) {
1438           
1439            e = E[i][j];
1440            if (e) {
1441               
1442                x = j;
1443                y = i;
1444               
1445                #ifdef OPENMP
1446                #pragma omp critical
1447                #endif
1448                {
1449                    // min max reduction
1450                    if (y < Stats[e].ymin) Stats[e].ymin = y;
1451                    if (y > Stats[e].ymax) Stats[e].ymax = y;
1452                }
1453                #ifdef OPENMP
1454                #pragma omp critical
1455                #endif
1456                {
1457                    if (x < Stats[e].xmin) Stats[e].xmin = x;
1458                    if (x > Stats[e].xmax) Stats[e].xmax = x;
1459                }
1460                // + reduction
1461                #ifdef OPENMP
1462                #pragma omp atomic
1463                #endif
1464                Stats[e].S += 1;
1465               
1466                #ifdef OPENMP
1467                #pragma omp atomic
1468                #endif
1469                Stats[e].Sx += x;
1470               
1471                #ifdef OPENMP
1472                #pragma omp atomic
1473                Stats[e].Sy += y;
1474                #endif
1475            } // if e
1476        } // j
1477    } // i
1478}
1479// --------------------------------------------------------------------------------------------
1480void imageFeaturesComputation_range_omp2(uint32** E, int height, int width, RegionStats* Stats)
1481// --------------------------------------------------------------------------------------------
1482{
1483    // version OpenMP 2.0
1484   
1485    int i, j;
1486   
1487    uint32 x, y;
1488    uint32 e;
1489   
1490   
1491#ifdef OPENMP
1492    //#pragma omp parallel for private(height, width, i, j, x, y, e) shared(E, Stats)
1493#pragma omp parallel for shared(E, Stats) private(height, width, i, j, x, y, e) schedule(dynamic)
1494//#pragma omp for private (j)
1495#endif   
1496    for (i = 0; i < height; i++) {
1497        //printf("i = %d\n", i);
1498        //printf("omp_get_num_threads = %d\n", omp_get_num_threads());
1499       
1500        for (j = 0; j < width; j++) {
1501           
1502            e = E[i][j];
1503            if (e) {
1504               
1505                x = j;
1506                y = i;
1507               
1508#ifdef OPENMP
1509#pragma omp critical
1510#endif
1511                {
1512                    // min max reduction
1513                    if (y < Stats[e].ymin)  Stats[e].ymin = y;
1514                    if (y > Stats[e].ymax)  Stats[e].ymax = y;
1515                   
1516                    if (x < Stats[e].xmin)  Stats[e].xmin = x;
1517                    if (x > Stats[e].xmax)  Stats[e].xmax = x;
1518                   
1519                    // + reduction
1520                    Stats[e].S  += 1;
1521                    Stats[e].Sx += x;
1522                    Stats[e].Sy += y;
1523                } // omp critical
1524            } // if e
1525        } // j
1526    } // i
1527}
1528// --------------------------------------------------------------------------------------------------------
1529void imageFeaturesComputation_omp4(uint32** restrict E, int height, int width, RegionStats* restrict Stats)
1530// --------------------------------------------------------------------------------------------------------
1531{
1532    // version avec "task" (OpenMP 3.0) et "depend" (OpenMP 4.0)
1533#ifdef OPENMP
1534#pragma omp parallel private(height,width) shared(E,Stats)
1535    {
1536#endif // OPENMP
1537        int i, j;
1538       
1539        uint32 x, y;
1540        uint32 e;
1541       
1542#ifdef OPENMP4
1543        //#pragma omp task depend ( in:E[0:height-1][0:width-1]) depend( inout: E[1:height*width/4])
1544#pragma omp task depend( inout: E[1:height*width/4])   
1545#endif
1546        for (i = 0; i < height; i++) {
1547            for (j = 0; j < width; j++) {
1548               
1549                e = E[i][j];
1550                if (e) {
1551                   
1552                    x = j;
1553                    y = i;
1554                   
1555                   
1556                    // min max reduction
1557                    if (y < Stats[e].ymin)  Stats[e].ymin = y;
1558                    if (y > Stats[e].ymax)  Stats[e].ymax = y;
1559                   
1560                    if (x < Stats[e].xmin)  Stats[e].xmin = x;
1561                    if (x > Stats[e].xmax)  Stats[e].xmax = x;
1562                   
1563                    // + reduction
1564                    Stats[e].S  += 1;
1565                    Stats[e].Sx += x;
1566                    Stats[e].Sy += y;
1567                }
1568            }
1569        }
1570#ifdef OPENMP
1571    }
1572#endif // OPENMP
1573}
1574// ------------------------------------------------------------------------------
1575void calc_xmin(uint32** restrict E, int height, int width, uint16* restrict Xmin)
1576// ------------------------------------------------------------------------------
1577{
1578    int i, j;
1579    uint32 x;
1580    uint32 e;
1581/*   
1582#ifdef OPENMP
1583#pragma omp critical
1584#endif
1585    { printf("calc xmin %d x %d\n", width, height); }
1586*/
1587    for (i = 0; i < height; i++) {
1588        for (j = 0; j < width; j++) {
1589            e = E[i][j];
1590            if (e) {
1591                x = j;
1592                if (x < Xmin[e]) Xmin[e] = x;
1593            }
1594        }
1595    }
1596}
1597// ------------------------------------------------------------------------------
1598void calc_xmax(uint32** restrict E, int height, int width, uint16* restrict Xmax)
1599// ------------------------------------------------------------------------------
1600{
1601    int i, j;
1602    uint32 x;
1603    uint32 e;
1604/*   
1605#ifdef OPENMP
1606#pragma omp critical
1607#endif
1608    { printf("calc xmax %d x %d\n", width, height); }
1609 */
1610    for (i = 0; i < height; i++) {
1611        for (j = 0; j < width; j++) {
1612            e = E[i][j];
1613            if (e) {
1614                x = j;
1615                if (x > Xmax[e]) Xmax[e] = x;
1616            }
1617        }
1618    }
1619}
1620// ------------------------------------------------------------------------------
1621void calc_ymin(uint32** restrict E, int height, int width, uint16* restrict Ymin)
1622// ------------------------------------------------------------------------------
1623{
1624    int i, j;
1625    uint32 y;
1626    uint32 e;
1627/*   
1628#ifdef OPENMP
1629#pragma omp critical
1630#endif
1631    { printf("calc ymin %d x %d\n", width, height); }
1632*/
1633    for(i=0; i<height; i++) {
1634        for(j=0; j<width; j++) {
1635            e = E[i][j];
1636            if(e) {
1637                y = i;
1638                if(y < Ymin[e]) Ymin[e] = y;
1639            }
1640        }
1641    }
1642}
1643// ------------------------------------------------------------------------------
1644void calc_ymax(uint32** restrict E, int height, int width, uint16* restrict Ymax)
1645// ------------------------------------------------------------------------------
1646{
1647    int i, j;
1648    uint32 y;
1649    uint32 e;
1650/*
1651#ifdef OPENMP
1652#pragma omp critical
1653#endif
1654    { printf("calc ymax %d x %d\n", width, height); }
1655*/
1656    for(i=0; i<height; i++) {
1657        for(j=0; j<width; j++) {
1658            e = E[i][j];
1659            if(e) {
1660                y = i;
1661                if(y > Ymax[e]) Ymax[e] = y;
1662            }
1663        }
1664    }
1665}
1666// ------------------------------------------------------------------------
1667void calc_s(uint32** restrict E, int height, int width, uint32* restrict S)
1668// ------------------------------------------------------------------------
1669{
1670    int i, j;
1671    uint32 e;
1672/*   
1673#ifdef OPENMP
1674#pragma omp critical
1675#endif 
1676    { printf("calc s %d x %d\n", width, height); }
1677*/
1678    for(i=0; i<height; i++) {
1679        for(j=0; j<width; j++) {
1680            e = E[i][j];
1681            if(e) {
1682                S[e] += 1;
1683            }
1684        }
1685    }
1686}
1687// --------------------------------------------------------------------------
1688void calc_sx(uint32** restrict E, int height, int width, uint32* restrict Sx)
1689// --------------------------------------------------------------------------
1690{
1691    int i, j;
1692    uint32 e;
1693/*   
1694#ifdef OPENMP
1695#pragma omp critical
1696#endif 
1697    { printf("calc sx %d x %d\n", width, height); }
1698*/
1699    for(i=0; i<height; i++) {
1700        for(j=0; j<width; j++) {
1701            e = E[i][j];
1702            if(e) {
1703                Sx[e] += j;
1704            }
1705        }
1706    }
1707}
1708// --------------------------------------------------------------------------
1709void calc_sy(uint32** restrict E, int height, int width, uint32* restrict Sy)
1710// --------------------------------------------------------------------------
1711{
1712    int i, j;
1713    uint32 e;
1714/*
1715#ifdef OPENMP
1716#pragma omp critical
1717#endif
1718    { printf("calc sy %d x %d\n", width, height); }
1719 */
1720    for (i = 0; i < height; i++) {
1721        for (j = 0; j < width; j++) {
1722            e = E[i][j];
1723            if (e) {
1724                Sy[e] += i;
1725            }
1726        }
1727    }
1728}
1729// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
1730void imageFeaturesComputation_omp5(uint32** E, int height, int width, uint16* Xmin, uint16* Xmax, uint16* Ymin, uint16* Ymax, uint32* S, uint32* Sx, uint32* Sy)
1731// ---------------------------------------------------------------------------------------------------------------------------------------------------------------
1732{
1733    // version avec "task" (OpenMP 3.0) et "depend" (OpenMP 4.0)
1734#ifdef OPENMP
1735#pragma omp parallel shared(E,Xmin,Xmax,Ymin,Ymax,S,Sx,Sy)
1736    // ne sourtout pas mettre height et width en private
1737    {
1738#endif // OPENMP
1739
1740        int id; // thread number
1741       
1742#ifdef OPENMP
1743        id = omp_get_thread_num();
1744#else
1745        id = 1;
1746#endif
1747       
1748        //printf("thread id = %d h = %d w = %d\n", id, height, width);
1749       
1750        if (id == 0) { calc_xmin(E, height, width, Xmin); }
1751        if (id == 1) { calc_xmax(E, height, width, Xmax); }
1752        if (id == 2) { calc_ymin(E, height, width, Ymin); }
1753        if (id == 3) { calc_ymax(E, height, width, Ymax); }
1754       
1755        if (id == 4) { calc_s (E, height, width, S);  }
1756        if (id == 5) { calc_sx(E, height, width, Sx); }
1757        if (id == 6) { calc_sy(E, height, width, Sy); }
1758#ifdef OPENMP
1759    }
1760#endif // OPENMP
1761}
1762// ------------------------------------------------------
1763int RegionStats_Compare(RegionStats *S1, RegionStats *S2)
1764// ------------------------------------------------------
1765{
1766    //puts("----------------------------------------");
1767    //display_RegionStats(S1, "S1");
1768    //display_RegionStats(S2, "S2");
1769    if((S1->xmin == S2->xmin) &&
1770       (S1->xmax == S2->xmax) &&
1771       (S1->ymin == S2->ymin) &&
1772       (S1->ymax == S2->ymax) &&
1773       (S1->S    == S2->S   ) &&
1774       (S1->Sx   == S2->Sx  ) &&
1775       (S1->Sy   == S2->Sy  ))
1776        return 1;
1777    else 
1778        return 0;
1779}
1780// ----------------------------------------------------------------------------
1781int RegionStatsVector_Compare(RegionStats *S1, int i0, int i1, RegionStats *S2)
1782// ----------------------------------------------------------------------------
1783{
1784    int i;
1785    int c; // resultat de la comparaison 0 = identique, 1 = different
1786    int s = 0; // somme
1787   
1788    for(i=i0; i<=i1; i++) {
1789        c = RegionStats_Compare(&S1[i], &S2[i]);
1790        s += c;
1791       
1792        /*if(c) {
1793            puts("---------------------------------------------------");
1794            printf("e = %d\n", i);
1795            display_RegionStats(&S1[i], NULL);
1796            display_RegionStats(&S2[i], NULL);
1797        }*/
1798           
1799    }
1800    return s;
1801}
1802// ------------------------------------------------------------------------------------------
1803int RegionStatsVector_Match(RegionStats *S1, int i0, int i1, RegionStats *S2, int j0, int j1)
1804// ------------------------------------------------------------------------------------------
1805{
1806    int i, j, pos;
1807    int c; // resultat de la comparaison 1 = identique, 0 = different
1808    int a; // accumulateur de c
1809    int s = 0; // somme
1810    int perm = 0; // permutation de numero de features
1811    int n1 = i1-i0+1;
1812    int n2 = j1-j0+1;
1813   
1814    //printf("[RegionStatsVector_Match]: [%d..%d]=%d vs [%d..%d]=%d\n", i0, i1, n1, j0,j1,n2);
1815    if (n1 != n2) {
1816        printf("card(S1) = %d neq card(S2) = %d  ", n1, n2);
1817        return 1;
1818    }
1819   
1820    for (i = i0; i <= i1; i++) {
1821        a   =  0;
1822        pos = -1;
1823       
1824        for (j = j0; j <= j1; j++) {
1825            c = RegionStats_Compare(&S1[i], &S2[j]);
1826            a = a + c;
1827            if (c) pos = j;
1828        }
1829        s += a;
1830       
1831        if (a > 1) {
1832            printf("erreur: il y a plusieurs fois la composante S1[%d] dans S2\n", i);
1833            for (j = j0; j <= j1; j++) {
1834                c = RegionStats_Compare(&S1[i], &S2[j]);
1835                if (c) printf("S2[%d] ", j);
1836            }
1837            printf("");
1838            giet_pthread_exit("");
1839        }
1840       
1841        if (i != pos) {
1842            //printf("perm(%d,%d)", i, pos);
1843            perm++;
1844        }
1845        if (a) {
1846            //printf("S1[%d] = S2[%d]\n", i, pos);
1847        }
1848        else {
1849            //printf("S1[%d] not matched\n", i);
1850        }
1851       
1852    }
1853    //printf("%4d", n1 - s); // le nb d'erreur
1854    printf("perm = %d  ", perm);
1855    return n1 - s;
1856}
Note: See TracBrowser for help on using the repository browser.