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

// ----------------------
// --- ecc_features.c ---
// ----------------------

/*
 * Copyright (c) 2012 - 2014, Lionel Lacassagne, All rights reserved
 * University of Paris Sud, Laboratoire de Recherche en Informatique 
 */

// Caracteristiques d'une region / label

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <malloc.h>
#include <string.h>


#include "nrc_os_config.h"
#include "config.h"
#include "nrc.h"


#if TARGET_OS == LINUX
    #include <sys/types.h>
    #include <sys/stat.h>
    #include <fcntl.h>
    #include <unistd.h>
#endif


#include "ecc_features.h"

#define BUFF_SIZE 1024 // for snprintf


// -----------------------------------------------------
void RegionStats_Init(RegionStats * Stats, uint32 nemax)
// -----------------------------------------------------
{
    for (int i = 0; i < (int) nemax; i++) {
#if PARMERGE
        pthread_spin_init(&Stats[i].lock, PTHREAD_PROCESS_PRIVATE);
#endif
    }
}


// -------------------------------------------------------------
void RegionStats_Constructor(RegionStats ** Stats, uint32 nemax)
// -------------------------------------------------------------
{
    *Stats = RegionStats_pConstructor(nemax);
}


// -------------------------------------------------
RegionStats * RegionStats_pConstructor(uint32 nemax)
// -------------------------------------------------
{
    RegionStats * Stats;
    
    Stats = (RegionStats *) malloc((nemax) * sizeof(RegionStats));
    if (Stats == NULL) {
        nrerror("allocation failed in RegionStats_pConstructor");
    }
    
    RegionStats_Init(Stats, nemax);
    
    return Stats;
}


// ------------------------------------------------------------
void RegionStats_Destructor(RegionStats ** Stats, uint32 nemax)
// ------------------------------------------------------------
{
    RegionStats_pDestructor(*Stats, nemax);
}


// ------------------------------------------------------------
void RegionStats_pDestructor(RegionStats * Stats, uint32 nemax)
// ------------------------------------------------------------
{
    //RegionStats_Clear(Stats, nemax);
    free(Stats);
}

// ------------------------------------------------------
void RegionStats_Clear(RegionStats * Stats, uint32 nemax)
// ------------------------------------------------------
{
    for (int i = 0; i < (int) nemax; i++) {
#if FEATURES
        Stats[i].xmin = 65535;
        Stats[i].xmax = 0;
        Stats[i].ymin = 65535;
        Stats[i].ymax = 0;
        
        Stats[i].S = 0;
        
        Stats[i].Sx = 0;
        Stats[i].Sy = 0;
#endif
    }
}

#if FEATURES
// -----------------------------------------
void RegionStats_Clear1(RegionStats * stats)
// -----------------------------------------
{
    stats->xmin = 0;
    stats->xmax = 0;
    stats->ymin = 0;
    stats->ymax = 0;
    
    stats->S = 0;
    
    stats->Sx = 0;
    stats->Sy = 0;
}


// -----------------------------------------
int RegionStats_Create_File(char * filename)
// -----------------------------------------
{
    int fd;
    
    fd = open(filename, O_CREAT | O_TRUNC);
    if (fd < 0) {
        printf("RegionStats_Open_File : can't create file %s\n", filename);
        exit(1);
    }
    return fd;
}


// ---------------------------------------
int RegionStats_Open_File(char * filename)
// ---------------------------------------
{
    int fd;
    
    fd = open(filename, O_RDONLY);
    if (fd < 0) {
        printf("RegionStats_Open_File : can't open file %s\n", filename);
        exit(1);
    }
    return fd;
}


// --------------------------------
void RegionStats_Close_File(int fd)
// --------------------------------
{
    close(fd);
}


#if 0 // pb : fscanf requires manipulating FILE *
// -------------------------------- 
int RegionStats_Read_Header(int fd)
// -------------------------------- 
{
    int ne = 0;
    fscanf(fd, "%d", &ne);
    return ne;
}
#endif


// ------------------------------------------ 
void RegionStats_Write_Header(int ne, int fd)
// ------------------------------------------ 
{
    char buff[BUFF_SIZE];
    snprintf(buff, BUFF_SIZE, "%d\n", ne);
    write(fd, buff, strlen(buff) + 1);
}


#if 0
// -------------------------------------------------------------- 
void RegionStats_Read_Stats1(int fd, int ne, RegionStats * Stats)
// -------------------------------------------------------------- 
{
    int i;
    
    for (i = 1; i <= ne; i++) {
        fscanf(fd, "%d%d%d%d%d%d%d%d\n",
               &t,
               &(Stats[i].xmin),
               &(Stats[i].xmax),
               &(Stats[i].ymin),
               &(Stats[i].ymax),
               &(Stats[i].S),
               &(Stats[i].Sx),
               &(Stats[i].Sy));
    }
}
#endif


// ---------------------------------------------------------------
void RegionStats_Write_Stats1(RegionStats * Stats, int ne, int fd)
// ---------------------------------------------------------------
{
    char buff[BUFF_SIZE];
    
    for (int i = 1; i <= ne; i++) {
        snprintf(buff, BUFF_SIZE, "%4d %5d %5d %5d %5d %7d %8d %8d\n",
                i,
                Stats[i].xmin,
                Stats[i].xmax,
                Stats[i].ymin,
                Stats[i].ymax,
                
                Stats[i].S,
                Stats[i].Sx,
                Stats[i].Sy);
        write(fd, buff, strlen(buff) + 1);
    }
}


// ---------------------------------------------------------------------------------------------------
void RegionStats_Write_Stats1_Sparse(RegionStats * Stats, uint32 * EQ, uint32 ne0, uint32 ne1, int fd)
// ---------------------------------------------------------------------------------------------------
{
    uint32 e;
    char buff[BUFF_SIZE];
    
    for (e = ne0; e <= ne1; e++) {
        if ((e == EQ[e]) && (Stats[e].S > 0)) {
            snprintf(buff, BUFF_SIZE, "%4d %5d %5d %5d %5d %7d %8d %8d\n",
                    e,
                    Stats[e].xmin,
                    Stats[e].xmax,
                    Stats[e].ymin,
                    Stats[e].ymax,
                    
                    Stats[e].S,
                    Stats[e].Sx,
                    Stats[e].Sy);
            write(fd, buff, strlen(buff) + 1);
        }
    }
}


// -----------------------------------------------------------------
void RegionStats_Write_pStats1(RegionStats ** Stats, int ne, int fd)
// -----------------------------------------------------------------
{
    char buff[BUFF_SIZE];
    
    for (int i = 1; i <= ne; i++) {
        snprintf(buff, BUFF_SIZE, "%4d %5d %5d %5d %5d %7d %8d %8d\n",
                i,
                Stats[i]->xmin,
                Stats[i]->xmax,
                Stats[i]->ymin,
                Stats[i]->ymax,
                Stats[i]->S,
                Stats[i]->Sx,
                Stats[i]->Sy);
        write(fd, buff, strlen(buff) + 1);
    }
}


#if 0
// --------------------------------------------------------------------------
void RegionStats_Load_Stats1(char * filename, int * ne, RegionStats ** Stats)
// --------------------------------------------------------------------------
{
    int fd;
    
    fd = RegionStats_Open_File(filename);
    
    *ne = RegionStats_Read_Header(fd);
    
    RegionStats_Constructor(Stats, *ne);
    RegionStats_Read_Stats1(fd, *ne, *Stats);
    RegionStats_Close_File(fd);
}


// --------------------------------------------------------------------------
void RegionStats_MLoad_Stats1(char * filename, int * ne, RegionStats * Stats)
// --------------------------------------------------------------------------
{
    int fd;
    
    fd = RegionStats_Open_File(filename);
    
    *ne = RegionStats_Read_Header(fd);
    RegionStats_Read_Stats1(fd, *ne, Stats);
    RegionStats_Close_File(fd);
}
#endif

// -----------------------------------------------------------------------
void RegionStats_Save_Stats1(RegionStats * Stats, int ne, char * filename)
// -----------------------------------------------------------------------
{
    int fd;
    
    fd = RegionStats_Create_File(filename);
    
    RegionStats_Write_Header(ne, fd);
    RegionStats_Write_Stats1(Stats, ne, fd);
    RegionStats_Close_File(fd);
}


// -------------------------------------------------------------------------
void RegionStats_Save_pStats1(RegionStats ** Stats, int ne, char * filename)
// -------------------------------------------------------------------------
{
    int fd;
    
    fd = RegionStats_Create_File(filename);
    
    RegionStats_Write_Header(ne, fd);
    RegionStats_Write_pStats1(Stats, ne, fd);
    RegionStats_Close_File(fd);
}


// ----------------------------------------------------------------------
void RegionStats_Display_Stats1(RegionStats * Stats, int ne, char * name)
// ----------------------------------------------------------------------
{
    if (name != NULL) {
        printf("%s : %d\n", name, ne);
    }
    else {
        printf("RegionStats : %d\n", ne);
    }
    for (int i = 1; i <= ne; i++) {
        printf("#%3d: %4d %4d %4d %4d %6d %8d %8d\n",
               i,
               Stats[i].xmin,
               Stats[i].xmax,
               Stats[i].ymin,
               Stats[i].ymax,
               Stats[i].S,
               Stats[i].Sx,
               Stats[i].Sy);
    }
}


// ------------------------------------------------------------------------
void RegionStats_Display_pStats1(RegionStats ** Stats, int ne, char * name)
// ------------------------------------------------------------------------
{
    if (name != NULL) {
        printf("%s : %d\n", name, ne);
    }
    else {
        printf("RegionStats : %d\n", ne);
    }
    for (int i = 1; i <= ne; i++) {
        printf("#%3d: %4d %4d %4d %4d %6d %8d %8d\n",
               i,
               Stats[i]->xmin,
               Stats[i]->xmax,
               Stats[i]->ymin,
               Stats[i]->ymax,
               Stats[i]->S,
               Stats[i]->Sx,
               Stats[i]->Sy);
    }
}


// ----------------------------------------------------------------------------------------------
void RegionStats_SetRectangle(RegionStats * Stats, int e, int ymin, int ymax, int xmin, int xmax)
// ----------------------------------------------------------------------------------------------
{
    Stats[e].ymin = ymin;
    Stats[e].xmin = xmin;
    Stats[e].ymax = ymax;
    Stats[e].xmax = xmax;
}


// ---------------------------------------------------------
void RegionStats_Copy1(RegionStats * src, RegionStats * dst)
// ---------------------------------------------------------
{
    dst->xmin = src->xmin;
    dst->xmax = src->xmax;
    dst->ymin = src->ymin;
    dst->ymax = src->ymax;
    
    dst->S    = src->S;
    
    dst->Sx   = src->Sx;
    dst->Sy   = src->Sy;
}
#endif

// ===============================
// === nouvelles versions 2009 ===
// ===============================

// --------------------------------------------
RegionStats * RegionStatsVector(int i0, int i1)
// --------------------------------------------
// allocate a float vector with subscript range v[i0..i1]
{
    RegionStats * v;
    
    v = (RegionStats *) malloc((size_t) ((i1 - i0 + 1 + NR_END) * sizeof(RegionStats)));
    if (!v) {
        nrerror("allocation failure in %s()", __func__);
        return NULL;
    }
    RegionStats_Init(v, i1 - i0 + 1 + NR_END);
    return v - i0 + NR_END;
}


#if TARGET_OS == GIETVM
// -----------------------------------------------------------------
RegionStats * remote_RegionStatsVector(int i0, int i1, int x, int y)
// -----------------------------------------------------------------
// allocate a float vector with subscript range v[i0..i1]
{
    RegionStats * v;
    
    v = (RegionStats *) remote_malloc((size_t) ((i1 - i0 + 1 + NR_END) * sizeof(RegionStats)), x, y);
    if (!v) {
        nrerror("allocation failure in %s()", __func__);
        return NULL;
    }
    RegionStats_Init(v, i1 - i0 + 1 + NR_END);
    return v - i0 + NR_END;
}
#endif


// ---------------------------------------------
RegionStats * RegionStatsVector0(int i0, int i1)
// ---------------------------------------------
// allocate a float vector with subscript range v[i0..i1]
{
    RegionStats * v;
    
    v = (RegionStats *) calloc((size_t) (i1 - i0 + 1 + NR_END), sizeof(RegionStats));
    if (!v) {
        nrerror("allocation failure in RegionStatsVector0()");
        return NULL;
    }
    return v - i0 + NR_END;
}


// ---------------------------------------------------------
void free_RegionStatsVector(RegionStats * v, int i0, int i1)
// ---------------------------------------------------------
// free a RegionStats vector allocated with vector()
{
    free(v + i0 - NR_END);
}


// -------------------------------------------------------------
RegionStats ** RegionStatsMatrix(int i0, int i1, int j0, int j1)
// -------------------------------------------------------------

// allocate a RegionStats matrix with subscript range m[nrl..nrh][ncl..nch]
{
    long nrow = i1 - i0 + 1;
    long ncol = j1 - j0 + 1;
    RegionStats ** m;
    
    // allocate pointers to rows
    m = (RegionStats **) malloc((size_t) ((nrow + NR_END) * sizeof(RegionStats *)));
    if (!m) {
        nrerror("allocation failure 1 in RegionStatsMatrix()");
    }
    m += NR_END;
    m -= i0;
    
    // allocate rows and set pointers to them
    m[i0] = (RegionStats *) malloc((size_t) ((nrow * ncol + NR_END) * sizeof(RegionStats)));
    if (!m[i0]) {
        nrerror("allocation failure 2 in RegionStatsMatrix()");
    }
    m[i0] += NR_END;
    m[i0] -= j0;
    
    for (int i = i0 + 1; i <= i1; i++) {
        m[i] = m[i - 1] + ncol;
    }
    
    // return pointer to array of pointers to rows
    return m;
}


// --------------------------------------------------------------
RegionStats ** RegionStatsMatrix0(int i0, int i1, int j0, int j1)
// --------------------------------------------------------------
// allocate a float matrix with subscript range m[nrl..nrh][ncl..nch]
{
    long i;
    long nrow = i1 - i0 + 1;
    long ncol = j1 - j0 + 1;
    RegionStats ** m;
    
    // allocate pointers to rows
    m= (RegionStats **) malloc((size_t) ((nrow + NR_END) * sizeof(RegionStats*)));
    if (!m) {
        nrerror("allocation failure 1 in RegionStatsMatrix()");
    }
    m += NR_END;
    m -= i0;
    
    // allocate rows and set pointers to them
    m[i0] = (RegionStats *) calloc((size_t) (nrow * ncol + NR_END), sizeof(RegionStats));
    if (!m[i0]) {
        nrerror("allocation failure 2 in RegionStatsMatrix()");
    }
    m[i0] += NR_END;
    m[i0] -= j0;
    
    for (i = i0 + 1; i <= i1; i++) {
        m[i] = m[i - 1] + ncol;
    }
    
    // return pointer to array of pointers to rows
    return m;
}


// --------------------------------------------------------------------------
void free_RegionStatsMatrix(RegionStats ** m, int i0, int i1, int j0, int j1)
// --------------------------------------------------------------------------
{
    free(m[i0] + j0 - NR_END);
    free(m + i0 - NR_END);
}


#if FEATURES
// -----------------------------------
void zero_RegionStats(RegionStats * x)
// -----------------------------------
{
    x->xmin = 32767;
    x->xmax = 0;
    x->ymin = 32767;
    x->ymax = 0;
    
    x->S  = 0;
    x->Sx = 0;
    x->Sy = 0;
}


// ---------------------------------------------------------
void zero_RegionStatsVector(RegionStats * v, int i0, int i1)
// ---------------------------------------------------------
{
    for (int i = i0; i <= i1; i++) {
        zero_RegionStats(&v[i]);
    }
}


// --------------------------------------------------------------------------
void zero_RegionStatsMatrix(RegionStats ** m, int i0, int i1, int j0, int j1)
// --------------------------------------------------------------------------
{
    for (int i = i0; i <= i1; i++) {
        for (int j = j0; j <= j1; j++) {         
            zero_RegionStats(&(m[i][j]));
        }
    }
}


// ---------------------------------------------------
void display_RegionStats(RegionStats * x, char * name)
// ---------------------------------------------------
{
    if (name != NULL) {
        printf("%s : \n", name);
    }
 
    printf("%4d %4d %4d %4d %6d %8d %8d\n",
           x->xmin,
           x->xmax,
           x->ymin,
           x->ymax,
           
           x->S,
           x->Sx,
           x->Sy);
}


// ------------------------------------------------------------------------
void display_RegionStatsVector(RegionStats * v, int i0, int i1, char *name)
// ------------------------------------------------------------------------
{
    if (name != NULL) {
        printf("%s : [%d..%d]\n", name, i0, i1);
    }
    else {
        printf("RegionStats : [%d..%d]\n", i0, i1);
    }
    for (int i = i0; i <= i1; i++) {
        printf("#%3d: ", i);
        display_RegionStats(&(v[i]), NULL);
    }
}


// ------------------------------------------------------------------------------------------
void display_RegionStatsMatrix(RegionStats ** m, int i0, int i1, int j0, int j1, char * name)
// ------------------------------------------------------------------------------------------
{
    
    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);
    }
    for (int i = i0; i <= i1; i++) {
        for (int j = j0; j <= j1; j++) {
            printf("#%3d: ", i);
            display_RegionStats(&(m[i][j]), NULL);
        }
    }
}


// ------------------------------------------------
void save_RegionStats(RegionStats * x, char * name)
// ------------------------------------------------
{
    int fd = -1;
    char buff[BUFF_SIZE];
    
    if (name == NULL) {
        return;
    }
    // assume name != NULL if single element
    // assume name == NULL if vector or matrix
    fd = RegionStats_Create_File(name);
    if (fd <= 0) {
        printf("*** Erreur : ouverture du fichier %s dans %s\n", name, __func__);
    }
    snprintf(buff, BUFF_SIZE, "%s: %4d %4d %4d %4d %6d %8d %8d\n",
            name,
            x->xmin,
            x->xmax,
            x->ymin,
            x->ymax,
            
            x->S,
            x->Sx,
            x->Sy);
    write(fd, buff, strlen(buff) + 1);
    
    if (name) {
        RegionStats_Close_File(fd);
    }    
}


// ----------------------------------------------------------------------
void save_RegionStatsVector(RegionStats * v, int i0, int i1, char * name)
// ----------------------------------------------------------------------
{
    int fd;
    char buff[BUFF_SIZE];
    
    if (name == NULL) {
        name = "RegionStatsVector";
    }
    fd = RegionStats_Create_File(name);
    
    snprintf(buff, BUFF_SIZE, "%s : [%d..%d]\n", name, i0, i1);
    write(fd, buff, strlen(buff) + 1);
    
    for (int i = i0; i <= i1; i++) {
        printf("#%3d: ", i);
        save_RegionStats(&v[i], NULL);
    }
    RegionStats_Close_File(fd);
}


// ---------------------------------------------------------------------------------------
void save_RegionStatsMatrix(RegionStats ** m, int i0, int i1, int j0, int j1, char * name)
// ---------------------------------------------------------------------------------------
{
    int fd;
    char buff[BUFF_SIZE];
    
    if (name == NULL) {
        name = "RegionStatsMatrix";
    }
    fd = RegionStats_Create_File(name);
    
    snprintf(buff, BUFF_SIZE, "%s : [%d..%d]\n", name, i0, i1);
    write(fd, buff, strlen(buff) + 1);
    
    for (int i = i0; i <= i1; i++) {
        for (int j = j0; j <= j1; j++) {
            snprintf(buff, BUFF_SIZE, "#%3d: ", i);
            write(fd, buff, strlen(buff) + 1);
            save_RegionStats(&m[i][j], NULL);
        }
    }
    RegionStats_Close_File(fd);
}


// -------------------------------------------------------------------------------
void RegionStats_Calc1_Features_1Pass(RegionStats * Stats, uint32 e, int i, int j)
// -------------------------------------------------------------------------------
{
    // calcul sur 1 point et non sur toute l'image
        // Rectangle
        
    if (i < Stats[e].ymin) {
        Stats[e].ymin = i;
    }
        if (i > Stats[e].ymax) {
        Stats[e].ymax = i;
    }
    
    if (j < Stats[e].xmin) {
        Stats[e].xmin = j;
    }
        if (j > Stats[e].xmax) {
        Stats[e].xmax = j;      
    }
    
        // Moment1
        Stats[e].S  += 1;
        Stats[e].Sx += j;
        Stats[e].Sy += i;
        
        return;
}


// --------------------------------
// --- fonctions de 2013 ----------
// --------------------------------
// ---------------------------------------------------------------------------------------------
void RegionStats_Calc_Rectangle_Moment1(uint32 ** E, int height, int width, RegionStats * Stats)
// ---------------------------------------------------------------------------------------------
{
    uint32 x, y;
    uint32 e;
    
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width; j++) {
            
            e = E[i][j];
            if (e) {
                x = j;
                y = i;
                
                if (i < Stats[e].ymin) {
                    Stats[e].ymin = y;
                }
                if (i > Stats[e].ymax) {
                    Stats[e].ymax = y;
                }
                
                if (j < Stats[e].xmin) {
                    Stats[e].xmin = x;
                }
                if (j > Stats[e].xmax) {
                    Stats[e].xmax = x;
                }
                
                Stats[e].S  += 1;
                Stats[e].Sx += x;
                Stats[e].Sy += y;
            }
        }
    }
}


// -------------------------------------------------------------------------------------------------------------------------------
void RegionStats_calc_Status(RegionStats * Stats, uint32 ne, uint32 min_height, uint32 min_width, uint32 min_area, uint8 * status)
// -------------------------------------------------------------------------------------------------------------------------------
{
    uint16 xmin, xmax, ymin, ymax, xsize, ysize;
    uint32 size;
    uint32 e;
    
    for (e = 1; e < ne; e++) {
        
        ymin = Stats[e].ymin;
        ymax = Stats[e].ymax;
        
        xmin = Stats[e].xmin;
        xmax = Stats[e].xmax;
        
        ysize = ymax - ymin + 1;
        xsize = xmax - xmin + 1;
        size = xsize * ysize;
        
        if ((size > min_area) && (xsize >= min_width) && (ysize >= min_height)) {
            status[e] = 1;
        }
        else {
            status[e] = 0;
        }
    }
}


// ----------------------------------------------------------------------------
uint32 RegionStats_UpdateEQ_with_Status(uint8 * status, uint32 ne, uint32 * EQ)
// ----------------------------------------------------------------------------
{
    uint32 e;
    uint32 na = 0;
    for (e = 1; e < ne; e++) {
        if (status[e]) {
            EQ[e] = ++na;
        }
        else {
            EQ[e] = 0;
        }
    }
    return na;
}


// ------------------------------------------------------------------------------
void RegionStats_UpdateStats_with_EQ(uint32 * EQ, uint32 ne, RegionStats * Stats)
// ------------------------------------------------------------------------------
{
    uint32 e, a;
    for (e = 1; e < ne; e++) {
        a = EQ[e];
        if (a != e) {
            // copy
            RegionStats_Copy1(&Stats[e], &Stats[a]);
        }
        else {
            // do nothing
        }
    }
}


// -----------------------------------------------------------------------------
void featuresComputation(uint32 ** E, int height,int width, RegionStats * Stats)
// -----------------------------------------------------------------------------
{
    //uint32 nemax = height * width /2;   
    RegionStats_Calc_Rectangle_Moment1(E, height, width, Stats);    
}


// ---------------------------------------------------------------------------
void pointFeaturesComputation( uint32 ** E, int i, int j, RegionStats * Stats)
// ---------------------------------------------------------------------------
{
    uint32 x, y;
    uint32 e;
    
    e = E[i][j];
    if (e) {
        
        x = j;
        y = i;
        
        if (i < Stats[e].ymin) {
            Stats[e].ymin = y;
        }
        if (i > Stats[e].ymax) {
            Stats[e].ymax = y;
        }
        
        if (j < Stats[e].xmin) {
            Stats[e].xmin = x;
        }
        if (j > Stats[e].xmax) {
            Stats[e].xmax = x;
        }
        
        Stats[e].S  += 1;
        Stats[e].Sx += x;
        Stats[e].Sy += y;
    }
}


// -----------------------------------------------------------------------------
void lineFeaturesComputation(uint32 ** E, int i, int width, RegionStats * Stats)
// -----------------------------------------------------------------------------
{
    // line RegionStats_Calc_Rectangle_Moment1
    
    uint32 x, y;
    uint32 e;
    
    for (int j = 0; j < width; j++) {
        
        e = E[i][j];
        if (e) {
            
            x = j;
            y = i;
            
            if (i < Stats[e].ymin) {
                Stats[e].ymin = y;
            }
            if (i > Stats[e].ymax) {
                Stats[e].ymax = y;
            }
            
            if (j < Stats[e].xmin) {
                Stats[e].xmin = x;
            }
            if (j > Stats[e].xmax) {
                Stats[e].xmax = x;
            }
            
            Stats[e].S  += 1;
            Stats[e].Sx += x;
            Stats[e].Sy += y;
        }
    }
}


// --------------------------------------------------------------------------------------
void bandFeaturesComputation(uint32 ** E, int i0, int i1, int width, RegionStats * Stats)
// --------------------------------------------------------------------------------------
{
    for (int i = i0; i <= i1; i++) {
        lineFeaturesComputation(E, i, width, Stats);
    }
}


// -----------------------------------------------------------------------------------
void imageFeaturesComputation(uint32 ** E, int height, int width, RegionStats * Stats)
// -----------------------------------------------------------------------------------
{
    // image RegionStats_Calc_Rectangle_Moment1
    for (int i = 0; i < height; i++) {
        lineFeaturesComputation(E, i, width, Stats);
    }
}


// ---------------------------------------
// --- Fonctions 2014 --------------------
// ---------------------------------------

// ---------------------------------------------------------------------------------------
void RegionStats_Copy_Stats1_From_Index(RegionStats * Stats, int dst_index, int src_index)
// ---------------------------------------------------------------------------------------                                 
{
    // R[dst] = R[src]
    RegionStats_Copy1(&Stats[src_index], &Stats[dst_index]);
}


// ---------------------------------------------------------------------------------------------
void RegionStats_Accumulate_Stats1_From_Index(RegionStats * Stats, int dst_index, int src_index)
// ---------------------------------------------------------------------------------------------                                 
{
    // R[dst] += R[src]
    Stats[dst_index].xmin = ui16min2(Stats[dst_index].xmin, Stats[src_index].xmin);
    Stats[dst_index].xmax = ui16max2(Stats[dst_index].xmax, Stats[src_index].xmax);
    Stats[dst_index].ymin = ui16min2(Stats[dst_index].ymin, Stats[src_index].ymin);
    Stats[dst_index].ymax = ui16max2(Stats[dst_index].ymax, Stats[src_index].ymax);
    
    Stats[dst_index].S  += Stats[src_index].S;
    Stats[dst_index].Sx += Stats[src_index].Sx;
    Stats[dst_index].Sy += Stats[src_index].Sy;   
}


// ---------------------------------------------------------------------------------------------------------------------------
void RegionStats_DisplayStats_Sparse(uint32 * EQ, uint32 ne0, uint32 ne1, RegionStats * Stats, char * name, int * start_index)
// ---------------------------------------------------------------------------------------------------------------------------
{
    // n'affiche que les racines.
    // ne pas utiliser apres un pack, car le test n'a plus de sens
    uint32 e;
    uint32 na; // compteur
    
    if (name) {
        printf(name);
    }
    
    //na = RegionStats_Count_Roots_Sparse(Stats, EQ, ne0, ne1);
    //printf("%d\n", na);
    
    for (e = ne0; e <= ne1; e++) {
        if (e == EQ[e] && Stats[e].S > 0) {
            if (start_index != NULL) {
                printf("%5d ", *start_index);
                *start_index = *start_index + 1;
            }
            else {
                printf("%5d ", e);
            }
            display_RegionStats(&Stats[e], NULL);
        }
    }
}


// -------------------------------------------------------------------------------------------------------
void RegionStats_DisplayStats_Range(uint32 * EQ, uint32 ne0, uint32 ne1, RegionStats * Stats, char * name)
// -------------------------------------------------------------------------------------------------------
{
    // affichage dense (apres un pack)

    uint32 e;
    
    if (name) {
        printf(name);
    }
    
    for (e = ne0; e <= ne1; e++) {
        printf("%5d ", e);
        display_RegionStats(&Stats[e], NULL);
    }
}


// -----------------------------------------------------------------------------------------------------------
void RegionStats_Save_Stats1_Sparse(RegionStats * Stats, uint32 * EQ, uint32 ne0, uint32 ne1, char * filename)
// -----------------------------------------------------------------------------------------------------------
{
    int fd;
    uint32 na = 0;
    
    fd = RegionStats_Create_File(filename);
    na = RegionStats_Count_Roots_Sparse(Stats, EQ, ne0, ne1);
    
    RegionStats_Write_Header(na, fd);
    RegionStats_Write_Stats1_Sparse(Stats, EQ, ne0, ne1, fd);
    RegionStats_Close_File(fd);
}


// --------------------------------------------------------------------------------------------
uint32 RegionStats_Count_Roots_Sparse(RegionStats * Stats, uint32 * EQ, uint32 ne0, uint32 ne1)
// --------------------------------------------------------------------------------------------
{
    uint32 e, c = 0; // compteur
    
    for (e = ne0; e <= ne1; e++) {
        if ((e == EQ[e]) && (Stats[e].S > 0)) {
            c++;
        }
    }
    return c;
}


// -----------------------------------------------------------------------------------
uint32 RegionStats_Count_Roots_Sparse1(RegionStats * Stats, uint32 * EQ, uint32 nemax)
// -----------------------------------------------------------------------------------
{
    return RegionStats_Count_Roots_Sparse(Stats, EQ, 1, nemax);
}


// ---------------------------------------------------------------------------------------------
uint32 RegionStats_Count_Labels_Sparse(RegionStats * Stats, uint32 * EQ, uint32 ne0, uint32 ne1)
// ---------------------------------------------------------------------------------------------
{
    uint32 e, c = 0; // compteur
    
    for (e = ne0; e <= ne1; e++) {
        if (Stats[e].S > 0) {
            c++;
        }
    }
    return c;
}


// ------------------------------------------------------------------------------------
uint32 RegionStats_Count_Labels_Sparse1(RegionStats * Stats, uint32 * EQ, uint32 nemax)
// ------------------------------------------------------------------------------------
{
    return RegionStats_Count_Labels_Sparse(Stats, EQ, 1, nemax);
}


// -----------------------------------------------------------------------
void copy_features_ui32matrix(RegionStats * Stats, uint32 ne, uint32 ** m)
// -----------------------------------------------------------------------
{
    for (int i = 0; i <= (int) ne; i++) {
        
        m[i][0] = i;
        
        // 16 bits: clean but requires 2 sorts
        //m[i][1] = Stats[i].xmin;
        //m[i][2] = Stats[i].xmax;
        //m[i][3] = Stats[i].ymin;
        //m[i][4] = Stats[i].ymax;
        
        // 32 bits: dirty, but requires only 1 sort
        m[i][1] = (Stats[i].ymin << 16) | Stats[i].ymax;
        m[i][2] = (Stats[i].xmin << 16) | Stats[i].xmax;

        // 32 bits
        m[i][3] = Stats[i].S;
        m[i][4] = Stats[i].Sx;
        m[i][5] = Stats[i].Sy;
    }
}


// -----------------------------------------------------------------------
void copy_ui32matrix_features(uint32 ** m, uint32 ne, RegionStats * Stats)
// -----------------------------------------------------------------------
{
    for (int i = 0; i <= (int) ne; i++) {
        
        Stats[i].xmin = m[i][2] >> 16;
        Stats[i].xmax = m[i][2] & 0xffff;
        Stats[i].ymin = m[i][1] >> 16;
        Stats[i].ymax = m[i][1] & 0xffff;
        
        Stats[i].S  = m[i][3];
        Stats[i].Sx = m[i][4];
        Stats[i].Sy = m[i][5];
    }
}


// ----------------------------------------------------------------------------
void sortv_ui32matrix_col(uint32 ** m, int i0, int i1, int j0, int j1, int col)
// ----------------------------------------------------------------------------
{
    // nrsort2 for NRC2
    //sortv_ui32matrix_selection_min(m, i0, i1, j0, j1,  col);
    
    long nrl = i0;
    long nrh = i1;
    long nc  = col;
    
    /*
     * sort an matrix of int, with the selection algorithm.
     * the key is in column nc
     * the sort is performed, by doing a purmutation on the lines,
     * instead of copying the lines.
     */
        
    uint32 x, min, pos;
        uint32 * ptr;
        
        for (int i = nrl; i < nrh; i++) {
                min = m[i][nc];
                pos = i;
                for (int j = i + 1; j <= nrh; j++) {
                        x = m[j][nc];
                        if (x < min) {
                                min = x;
                                pos = j;
                        }
                } // j
                
                // permutation des pointeurs de ligne de la matrice
                ptr    = m[i];
                m[i]   = m[pos];
                m[pos] = ptr;
        } // i
}


// -------------------------------------------------------------
void RegionStats_SortFeatures(RegionStats * Stats, uint32 nemax)
// -------------------------------------------------------------
{
    uint32 ** m = NULL;
    
    m = ui32matrix(0, nemax, 0, 7);
    
    copy_features_ui32matrix(Stats, nemax, m);
    sortv_ui32matrix_col(m, 0, nemax, 0, 5, 1);
    copy_ui32matrix_features(m, nemax, Stats);
}


// --------------------------------------------------------------------------------
void calc_xmin(uint32 ** restrict E, int height, int width, uint16 * restrict Xmin)
// --------------------------------------------------------------------------------
{
    uint32 x;
    uint32 e;
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width; j++) {
            e = E[i][j];
            if (e) {
                x = j;
                if (x < Xmin[e]) {
                    Xmin[e] = x;
                }
            }
        }
    }
}


// --------------------------------------------------------------------------------
void calc_xmax(uint32 ** restrict E, int height, int width, uint16 * restrict Xmax)
// --------------------------------------------------------------------------------
{
    uint32 x;
    uint32 e;
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width; j++) {
            e = E[i][j];
            if (e) {
                x = j;
                if (x > Xmax[e]) Xmax[e] = x;
            }
        }
    }
}


// --------------------------------------------------------------------------------
void calc_ymin(uint32 ** restrict E, int height, int width, uint16 * restrict Ymin)
// --------------------------------------------------------------------------------
{
    uint32 y;
    uint32 e;
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width; j++) {
            e = E[i][j];
            if (e) {
                y = i;
                if (y < Ymin[e]) {
                    Ymin[e] = y;
                }
            }
        }
    }
}


// --------------------------------------------------------------------------------
void calc_ymax(uint32 ** restrict E, int height, int width, uint16 * restrict Ymax)
// --------------------------------------------------------------------------------
{
    uint32 y;
    uint32 e;
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width; j++) {
            e = E[i][j];
            if (e) {
                y = i;
                if (y > Ymax[e]) {
                    Ymax[e] = y;
                }
            }
        }
    }
}


// --------------------------------------------------------------------------
void calc_s(uint32 ** restrict E, int height, int width, uint32 * restrict S)
// --------------------------------------------------------------------------
{
    uint32 e;
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width; j++) {
            e = E[i][j];
            if (e) {
                S[e] += 1;
            }
        }
    }
}


// ----------------------------------------------------------------------------
void calc_sx(uint32 ** restrict E, int height, int width, uint32 * restrict Sx)
// ----------------------------------------------------------------------------
{
    uint32 e;
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width; j++) {
            e = E[i][j];
            if (e) {
                Sx[e] += j;
            }
        }
    }
}


// ----------------------------------------------------------------------------
void calc_sy(uint32 ** restrict E, int height, int width, uint32 * restrict Sy)
// ----------------------------------------------------------------------------
{
    uint32 e;
    for (int i = 0; i < height; i++) {
        for (int j = 0; j < width; j++) {
            e = E[i][j];
            if (e) {
                Sy[e] += i;
            }
        }
    }
}


// --------------------------------------------------------
int RegionStats_Compare(RegionStats * S1, RegionStats * S2)
// --------------------------------------------------------
{
    //puts("----------------------------------------");
    //display_RegionStats(S1, "S1");
    //display_RegionStats(S2, "S2");
    if ((S1->xmin == S2->xmin) &&
       (S1->xmax == S2->xmax) &&
       (S1->ymin == S2->ymin) &&
       (S1->ymax == S2->ymax) &&
       (S1->S  == S2->S) &&
       (S1->Sx == S2->Sx) &&
       (S1->Sy == S2->Sy)) {
        return 1;
    }
    else {
        return 0;
    }
}


// ------------------------------------------------------------------------------
int RegionStatsVector_Compare(RegionStats * S1, int i0, int i1, RegionStats * S2)
// ------------------------------------------------------------------------------
{
    int c; // resultat de la comparaison 0 = identique, 1 = different
    int s = 0; // somme
    
    for (int i = i0; i <= i1; i++) {
        c = RegionStats_Compare(&S1[i], &S2[i]);
        s += c;
    }
    return s;
}


// --------------------------------------------------------------------------------------------
int RegionStatsVector_Match(RegionStats * S1, int i0, int i1, RegionStats * S2, int j0, int j1)
// --------------------------------------------------------------------------------------------
{
    int j, pos;
    int c; // resultat de la comparaison 1 = identique, 0 = different
    int a; // accumulateur de c
    int s = 0; // somme
    int perm = 0; // permutation de numero de features
    int n1 = i1 - i0 + 1;
    int n2 = j1 - j0 + 1;
    
    //printf("[RegionStatsVector_Match]: [%d..%d]=%d vs [%d..%d]=%d\n", i0, i1, n1, j0,j1,n2);
    if (n1 != n2) {
        printf("card(S1) = %d neq card(S2) = %d  ", n1, n2);
        return 1;
    }
    
    for (int i = i0; i <= i1; i++) {
        a   =  0;
        pos = -1;
        
        for (j = j0; j <= j1; j++) {
            c = RegionStats_Compare(&S1[i], &S2[j]);
            a = a + c;
            if (c) {
                pos = j;
            }
        }
        s += a;
        
        if (a > 1) {
            printf("erreur: il y a plusieurs fois la composante S1[%d] dans S2\n", i);
            for (j = j0; j <= j1; j++) {
                c = RegionStats_Compare(&S1[i], &S2[j]);
                if (c) {
                    printf("S2[%d] ", j);
                }
            }
            printf("\n");
            exit(1);
        }
        
        if (i != pos) {
            //printf("perm(%d,%d)", i, pos);
            perm++;
        }
        if (a) {
            //printf("S1[%d] = S2[%d]\n", i, pos);
        }
        else {
            //printf("S1[%d] not matched\n", i);
        }
        
    }
    //printf("%4d", n1 - s); // le nb d'erreur
    printf("perm = %d  ", perm);
    return n1 - s;
}
#endif // FEATURES

// Local Variables:
// tab-width: 4
// c-basic-offset: 4
// c-file-offsets:((innamespace . 0)(inline-open . 0))
// indent-tabs-mode: nil
// End:

// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4