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

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