/*
    filtres-cor.c - Edouard Thiel - 12/01/2001

    Compilation :
        cc filtres-cor.c bsutil.c -o filtres-cor \
            `~/helium/helium-cfg --cflags --libs` -lm

    Usage :
        filtres-cor [image-in.pgm]
*/

#include <helium.h>
#include <math.h>
#include "bsutil.h"


/* Variables globales */
He_node *princ, *panel1, *text_load, *canvas, *panel2, *mess;
int c_rouge;

BsMap *gbm1 = NULL, *gbm2 = NULL, *gbm3 = NULL, *gbm4 = NULL;
XImage *gxi = NULL;
int sobel_ok = FALSE;


/*-------------------------- A L G O R I T H M E S --------------------------*/

/*
 * Fixe l'épaisseur du bord de bm à e, puis affecte les pixels qui sont dans
 * le bord de l'image, à la valeur de leur pixel miroir dans l'image.
 *
 * Exemple :  e = 2                            k j i j k l k j
 *                                             g f e f g h g f
 *   (les pixels de           A B C D          c b A B C D c b
 *    l'intérieur de          E F G H   --->   g f E F G H g f
 *    l'image sont            I J K L          k j I J K L k j
 *    en majuscule)                            g f e f g h g f
 *                                             c b a b c d c b
*/

void BordMiroir (BsMap *bm, int e)
{
    int x, y;
    
    SetBorderBsMap (bm, e);
    
    /* Haut et bas */
    for (y = -bm->bord; y < 0; y++)
    for (x = 0; x < bm->xsiz; x++) {
	bm->tab[y][x] = bm->tab[-y][x];
	bm->tab[bm->ysiz-1-y][x] = bm->tab[bm->ysiz-1+y][x];
    }
    
    /* Gauche, droite et coins */
    for (x = -bm->bord; x < 0; x++)
    for (y = -bm->bord; y < bm->ymax; y++) {
	bm->tab[y][x] = bm->tab[y][-x];
	bm->tab[y][bm->xsiz-1-x] = bm->tab[y][bm->xsiz-1+x];
    }
}


/*
 * Filtre moyen de bm1 dans bm2 sur un voisinage L*L
*/

void FiltreMoyen (BsMap *bm1, BsMap *bm2, int L)
{
    int x, y, ll = L/2, Lcarre = L*L;
    
    BordMiroir (bm1, ll);
    
    for (y = 0; y < bm1->ysiz; y++)
    for (x = 0; x < bm1->xsiz; x++)
    {
	int s = 0, xx, yy;
        
        for (yy = -ll; yy <= ll; yy++)
        for (xx = -ll; xx <= ll; xx++)
            s += bm1->tab[y + yy][x + xx];
        
        bm2->tab[y][x] = s / Lcarre;
    }
}


/*
 * Tri à bulle (optimisé par dernière permutation).
*/

void tri_bulle (int *v, int n)
{
    int i, j, tmp, dp = 0;  /* dp est l'indice de la dernière permutation */

    while (dp < n) {
        j = n;
        for (i = n-2; i >= dp; i--)
            if (v[i] > v[i+1]) {
                j = i+1;
                tmp = v[j]; v[j] = v[i]; v[i] = tmp;
            }
        dp = j;
    }
}


/*
 * Filtre médian de bm1 dans bm2 sur un voisinage L*L
*/

void FiltreMedian (BsMap *bm1, BsMap *bm2, int L)
{
    int x, y, ll = L/2, Lcarre = L*L;
    
    BordMiroir (bm1, ll);
    
    for (y = 0; y < bm1->ysiz; y++)
    for (x = 0; x < bm1->xsiz; x++)
    {
	int s = 0, v[121], xx, yy;
        
        for (yy = -ll; yy <= ll; yy++)
        for (xx = -ll; xx <= ll; xx++)
            v[s++] = bm1->tab[y + yy][x + xx];
        
        tri_bulle (v, Lcarre);
        
        bm2->tab[y][x] = v[Lcarre/2];
    }
}


/*
 * Filtre de Nagao de bm1 dans bm2
*/

void FiltreNagao (BsMap *bm1, BsMap *bm2)
{
    /* [Numéro fenêtre][Numéro point de la fenêtre] */

    short px[9][9] = { { -2, -1,  0, -2, -1,  0, -2, -1,  0 },
                       {  0,  1,  2,  0,  1,  2,  0,  1,  2 },
                       {  0,  1,  2,  0,  1,  2,  0,  1,  2 },
                       { -2, -1,  0, -2, -1,  0, -2, -1,  0 },
                       { -1,  0,  1, -1,  0,  1, -1,  0,  1 },
                       { -2, -1,  0,  1,  2, -1,  0,  1,  0 },
                       {  2,  1,  2,  0,  1,  2,  1,  2,  2 },
                       {  0, -1,  0,  1, -2, -1,  0,  1,  2 },
                       { -2, -2, -1, -2, -1,  0, -2, -1, -2 }  };
 
    short py[9][9] = { { -2, -2, -2, -1, -1, -1,  0,  0,  0 },
                       { -2, -2, -2, -1, -1, -1,  0,  0,  0 },
                       {  0,  0,  0,  1,  1,  1,  2,  2,  2 },
                       {  0,  0,  0,  1,  1,  1,  2,  2,  2 },
                       { -1, -1, -1,  0,  0,  0,  1,  1,  1 },
                       { -2, -2, -2, -2, -2, -1, -1, -1,  0 },
                       { -2, -1, -1,  0,  0,  0,  1,  1,  2 },
                       {  0,  1,  1,  1,  2,  2,  2,  2,  2 },
                       { -2, -1, -1,  0,  0,  0,  1,  1,  2 }  };

    int x, y, i, j, k;
    double s, c, variance[9], moyenne[9];
    
    BordMiroir (bm1, 2);
    
    for (y = 0; y < bm1->ysiz; y++)
    for (x = 0; x < bm1->xsiz; x++)
    {
        /* Calcul des variances sur les 9 fenêtres */
        for (i = 0; i < 9; i++)
        {
            /* Calcul moyenne sur la fenêtre i */
            s = 0;
            for (j = 0; j < 9; j++)
                s += bm1->tab[ y+py[i][j] ][ x+px[i][j] ];
            moyenne[i] = s/9;
            
            /* Calcul variance sur la fenêtre i */
            s = 0;
            for (j = 0; j < 9; j++) {
                c = bm1->tab[ y+py[i][j] ][ x+px[i][j] ] - moyenne[i];
                s += c*c;
            }
            variance[i] = s/9;
        }
        
        /* Recherche de la fenêtre k de variance minimale */
        k = 0;
	for (i = 1; i < 9; i++)
	    if (variance[k] > variance[i]) k = i;
            
        /* On prend la valeur moyenne de la fenêtre k */ 
        bm2->tab[y][x] = moyenne[k];
    }
}



/*
 * Filtre linéaire de bm1 dans bm2
*/

#define FILTREMAX 11

typedef struct {
    int taille;				/* voisinage centré taille*taille */
    double coef[FILTREMAX][FILTREMAX];  /* accès : coef[y][x] */
} Filtre;

void FiltreLineaire (BsMap *bm1, BsMap *bm2, Filtre *f)
{
    int x, y, ll = f->taille/2;

    BordMiroir (bm1, ll);
    
    for (y = 0; y < bm1->ysiz; y++)
    for (x = 0; x < bm1->xsiz; x++)
    {
	double s = 0; 
        int xx, yy;
        
        for (yy = -ll; yy <= ll; yy++)
        for (xx = -ll; xx <= ll; xx++)
            s += bm1->tab[y + yy][x + xx] * f->coef[yy+ll][xx+ll];
        
        bm2->tab[y][x] = s;
    }
}


/*
 * Filtre de Sobel en X de bm1 dans bm2
*/

void FiltreSobelX (BsMap *bm1, BsMap *bm2)
{
    Filtre fx = { 3, {{ -1,  0,  1 },
                      { -2,  0,  2 },
                      { -1,  0,  1 }}};
                      
    FiltreLineaire (bm1, bm2, &fx);
}


/*
 * Filtre de Sobel en Y de bm1 dans bm2
*/

void FiltreSobelY (BsMap *bm1, BsMap *bm2)
{
    Filtre fy = { 3, {{ -1, -2, -1 },
                      {  0,  0,  0 },
                      {  1,  2,  1 }}};
                      
    FiltreLineaire (bm1, bm2, &fy);
}


/*
 * Calcul du module du gradient à partir des images de Sobel en X et en Y
 * dans sobM
*/

void ModuleSobel (BsMap *sobX, BsMap *sobY, BsMap *sobM)
{
    int x, y;
    
    for (y = 0; y < sobM->ysiz; y++)
    for (x = 0; x < sobM->xsiz; x++) {
        double gx = sobX->tab[y][x], gy = sobY->tab[y][x];
        sobM->tab[y][x] = sqrt (gx*gx + gy*gy);
    }
}


/*
 * Cherche le min et le max dans bm, puis fait une expansion de
 * dynamique des niveaux de gris de [min..max] à [0..255] dans bm.
*/

void Normaliser (BsMap *bm)
{
    int x, y, min = bm->tab[0][0], max = bm->tab[0][0];
    
    for (y = 0; y < bm->ysiz; y++)
    for (x = 0; x < bm->xsiz; x++)
        if (min > bm->tab[y][x]) min = bm->tab[y][x];
        else if (max < bm->tab[y][x]) max = bm->tab[y][x];
        
    for (y = 0; y < bm->ysiz; y++)
    for (x = 0; x < bm->xsiz; x++)
        bm->tab[y][x] = (bm->tab[y][x] - min) * 255 / (max - min);
}


/*---------------------------- I N T E R F A C E ----------------------------*/

void init_couleurs ()
{
    c_rouge   = HeAllocRgb (255,   0,   0, he_black);
}

void voir_image (BsMap *bm)
{
    int l;
    
    /* On adapte la taille à l'image */
    HeSetX (canvas, 0);
    HeSetWidth  (canvas, bm ? bm->xsiz : 50);
    HeSetHeight (canvas, bm ? bm->ysiz : 50);
    HeJustify (panel2, canvas, HE_TOP);
    HeSetWidth (panel2, HeGetExtWidth(canvas));
    HeFit (princ);
    HeExpand (panel2, NULL, HE_RIGHT);
    
    l = HeGetWidth (princ) - HeGetExtWidth(canvas); 
    if (l > 0) HeSetX (canvas, l/2);
    
    HeSetMessageLabel (mess, bm ? bm->title : "");
    
    /* On met à jour gxi et on provoque un réaffichage */
    FreeXImage (gxi); 
    gxi = NewXImageFromBsMap (bm, LUT_VGA_GREY);
    HePostRepaint (canvas);
}

void canvas_repaint_proc (He_node *hn, Window win)
{
    if (gxi == NULL) HeDrawBg (canvas, he_black);
    else PaintXImage (gxi, win);
}


/*
 * Lors d'un clic de la souris dans le canvas, affiche le vecteur gradient du
 * point cliqué dans la couleur c_rouge. Les images du gradient en X et Y,
 * calculées avec Sobel, sont mémorisées dans les images globales gbm3 et gbm4.
*/

void canvas_event_proc (He_node *hn, He_event *hev)
{
    if (!sobel_ok) return;
    
    switch (hev->type) {
    
        case ButtonPress : {
            int gx, gy;
            
            gx = gbm3->tab[hev->sy][hev->sx],
            gy = gbm4->tab[hev->sy][hev->sx];
            XSetForeground (he_display, he_default_gc, c_rouge);
            XDrawLine (he_display, hev->win, he_default_gc,
                hev->sx, hev->sy, hev->sx + gx, hev->sy + gy);
        } break;
    }
}

void algos_proc (He_node *hn)
{
    char *nom = HeGetButtonLabel (hn);
    
    if (gbm1 == NULL) return;
    SetTitleBsMap (gbm2, nom);
    
    if (strcmp (nom, "Moy3") == 0) {
        FiltreMoyen (gbm1, gbm2, 3);
    } else if (strcmp (nom, "Moy5") == 0) {
        FiltreMoyen (gbm1, gbm2, 5);
    } else if (strcmp (nom, "Moy7") == 0) {
        FiltreMoyen (gbm1, gbm2, 7);
    } else if (strcmp (nom, "Med3") == 0) {
        FiltreMedian (gbm1, gbm2, 3);
    } else if (strcmp (nom, "Med5") == 0) {
        FiltreMedian (gbm1, gbm2, 5);
    } else if (strcmp (nom, "Med7") == 0) {
        FiltreMedian (gbm1, gbm2, 7);
    } else if (strcmp (nom, "Nagao") == 0) {
        FiltreNagao (gbm1, gbm2);
    } else if (strcmp (nom, "Sobel X") == 0) {
        FiltreSobelX (gbm1, gbm2);
        Normaliser (gbm2);
    } else if (strcmp (nom, "Sobel Y") == 0) {
        FiltreSobelY (gbm1, gbm2);
        Normaliser (gbm2);
    } else if (strcmp (nom, "Sobel Mod") == 0) {
        FiltreSobelX (gbm1, gbm3);
        FiltreSobelY (gbm1, gbm4);
        ModuleSobel (gbm3, gbm4, gbm2);
        Normaliser (gbm2);
        SetTitleBsMap (gbm2, "Cliquer pour voir le gradient");
        sobel_ok = TRUE;
    }
    voir_image (gbm2);
}

void quit_proc (He_node *hn)
{
    HeQuit(0);
}

void load_proc (He_node *hn)
{
    char *nom = HeGetTextValue (text_load);
    
    FreeBsMap(gbm1); gbm1 = NULL;
    FreeBsMap(gbm2); gbm2 = NULL;
    FreeBsMap(gbm3); gbm3 = NULL;
    FreeBsMap(gbm4); gbm4 = NULL;
    sobel_ok = FALSE;
    
    gbm1 = ReadBsMapFromPGM (nom, 0);
    if (gbm1 == NULL) {
	HeSimpleDialog (
	    HE_DIALOG_BELL,
	    HE_DIALOG_TITLE,   "Attention",
	    HE_DIALOG_MESSAGE, "Erreur de lecture du fichier",
	    HE_DIALOG_QUOTED,  nom,
	    HE_DIALOG_BUTTOK,  "Ok",
	    0);
    } else {
        gbm2 = NewBsMap (gbm1->ysiz, gbm1->xsiz, 0);
        gbm3 = NewBsMap (gbm1->ysiz, gbm1->xsiz, 0);
        gbm4 = NewBsMap (gbm1->ysiz, gbm1->xsiz, 0);
        if (gbm2 == NULL || gbm3 == NULL || gbm4 == NULL) HeQuit(1);
    }
    voir_image (gbm1);
}

void reaffi_proc (He_node *hn)
{
    voir_image (gbm1);
}


/*---------------------------- P R I N C I P A L ----------------------------*/

int main (int argc, char **argv)
{
    int code_sortie;

    HeInit (&argc, &argv);
    
    init_couleurs ();
    
    princ = HeCreateFrame();
    HeSetFrameLabel (princ, "Filtres (cor)");

    panel1 = HeCreatePanel (princ);
    
    text_load = HeCreateText (panel1);
    HeSetTextVisibleLen (text_load, 35);
    HeSetTextCompletion (text_load, TRUE);
    HeSetTextNotifyProc (text_load, load_proc);
    HeSetPanelLayout (panel1, HE_LINE_FEED);
    
    HeCreateButtonP (panel1, "Charger", load_proc, NULL);
    HeCreateButtonP (panel1, "Reafficher", reaffi_proc, NULL);
    HeCreateButtonP (panel1, "Quit", quit_proc, NULL);
    HeSetPanelLayout (panel1, HE_LINE_FEED);
    
    HeCreateButtonP (panel1, "Moy3", algos_proc, NULL);
    HeCreateButtonP (panel1, "Moy5", algos_proc, NULL);
    HeCreateButtonP (panel1, "Moy7", algos_proc, NULL);
    HeCreateButtonP (panel1, "Med3", algos_proc, NULL);
    HeCreateButtonP (panel1, "Med5", algos_proc, NULL);
    HeCreateButtonP (panel1, "Med7", algos_proc, NULL);
    HeSetPanelLayout (panel1, HE_LINE_FEED);
    
    HeCreateButtonP (panel1, "Nagao", algos_proc, NULL);
    HeCreateButtonP (panel1, "Sobel X", algos_proc, NULL);
    HeCreateButtonP (panel1, "Sobel Y", algos_proc, NULL);
    HeCreateButtonP (panel1, "Sobel Mod", algos_proc, NULL);
    HeSetPanelLayout (panel1, HE_LINE_FEED);
    
    HeFit(panel1); 
    
    canvas = HeCreateCanvas (princ);
    HeSetCanvasRepaintProc (canvas, canvas_repaint_proc);
    HeSetCanvasEventProc (canvas, canvas_event_proc);
    HeJustify (canvas, panel1, HE_TOP); HeMoveY (canvas, 5);
    
    panel2 = HeCreatePanel (princ);
    mess = HeCreateMessage (panel2);
    HeFit (panel2);
    HeJustify (panel2, canvas, HE_TOP);
    HeFit (princ);
   
    if (argc > 1) { 
        HeSetTextValue (text_load, argv[1]); 
        load_proc (text_load);
    } else voir_image (NULL);
    
    code_sortie = HeMainLoop (princ);

    FreeBsMap (gbm1);
    FreeBsMap (gbm2);
    FreeBsMap (gbm3);
    FreeBsMap (gbm4);
    FreeXImage (gxi); 

    return code_sortie;
}