/*
 * Medial Axis for Chamfer Distances in 2D:
 * computing Look-Up Tables and Neighbourhoods
 *
 * Copyright (C) Eric Remy and Edouard Thiel - March 2001
 *
 * http://www.lif-sud.univ-mrs.fr/~thiel/PRL2001/LutChamfer2D.c
 *
 * This program is free software under the terms of the 
 * GNU Lesser General Public License (LGPL) version 2.1.
 *
 * To compile:  cc LutChamfer2D.c -o LutChamfer2D
 * Usage:       ./LutChamfer2D
*/


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


/*
 * Types and macros. Modify macros to enlarge masks and look-up tables.
*/

typedef struct {
    int x, y, w;
} Weighting;

#define MAXWEIGHTING 100

typedef struct {
    Weighting vg[MAXWEIGHTING];
    int ng;
} MaskG;

#define INFINITY 65535

typedef unsigned short ImagePoint;
typedef ImagePoint *Image;

#define MAXLUTENTRY 20000

typedef int LookUpTable[MAXWEIGHTING][MAXLUTENTRY];


/*
 * Create a new image L*L and initialize the pixels to 0.
 *
 * Input:  L the side length.
 * Output: the pointer to the image, or NULL if memory allocation failed.
*/

Image NewImage (int L)
{
    int n = L*L;
    Image res;

    printf ("NewImage: allocating %dk ... ", n*sizeof(ImagePoint)/1024);
    res = calloc (n, sizeof(ImagePoint));
    if (res == NULL)
         printf ("malloc error\n");
    else printf ("ok\n");
    return res;
}


/*
 * Add a weigthing in the mask M.
 *
 * Input:  M the generator of the mask, (x,y,w) the weighting.
 * Output: the number of the weighting in the mask.
*/

int AddWeighting (MaskG *M, int x, int y, int w)
{
    int i = M->ng;

    if (i >= MAXWEIGHTING) {
        printf ("AddWeighting: weighting number limited to MAXWEIGHTING = %d\n",
            MAXWEIGHTING);
        return 0;
    }
    if (! (0 <= y && y <= x && 0 < x && 0 < w)) {
        printf ("AddWeighting: (x = %d, y = %d, w = %d)\n",
            x, y, w);
        printf ("  does not respect 0 <= y <= x, 0 < x, 0 < w\n");
        return 0;
    }
    M->vg[i].x = x;
    M->vg[i].y = y;
    M->vg[i].w = w;
    M->ng++;
    return i;
}


/*
 * Fast Cone Distance Transform Algorithm.
 *
 * Input:  L the side length, MgC the generator of the dC mask.
 * Output: CTg the L^2 distance image to the origin.
*/

void ComputeCTg (int L, MaskG *MgC, Image CTg)
{
    int x, y, xx, yy, i, min, val;

    CTg[0] = 0;

    for (x = 1; x < L; x++)
    for (y = 0; y <= x; y++)
    {
        min = INFINITY;

        for (i = 0; i < MgC->ng; i++)
        {
            xx = x - MgC->vg[i].x;
            yy = y - MgC->vg[i].y;

            if (0 <= yy && yy <= xx)
            {
                val = CTg[ yy*L + xx ] + MgC->vg[i].w;

                if (val < min) min = val;
            }
        }

        CTg[ y*L + x ] = min;
    }
}


/*
 * Fast Distance Transform in Z^2/8.
 *
 * Input:  MgC the generator of the dC mask, L the side length,
 *         DTg the shape (limited to Z^2/8).
 * Output: DTg the distance image to the border of the shape.
*/

void ComputeDTg (MaskG *MgC, int L, Image DTg)
{
    int x, y, xx, yy, min, val, i;

    for (y = L-1; y >= 0; y--)
    for (x = L-1; x >= y; x--)
    {
        if (DTg[ y*L + x ] != 0)
        {
            min = INFINITY;

            for (i = 0; i < MgC->ng; i++)
            {
                xx = x + MgC->vg[i].x;
                yy = y + MgC->vg[i].y;

                if (xx < L)
                {
                    val = DTg[ yy*L + xx ] + MgC->vg[i].w;
                    if (val < min) min = val;
                }
            }

            DTg[ y*L + x ] = min;
        }
    }
}


/*
 * Lut column Computation Algorithm.
 *
 * Input:  CTg the distance cone to origin, L the side length, vg the
 *         weighting and ivg its number, Rmax the greatest radius
 *         to be verified in Lut.
 * Output: The Lut column Lut[ivg] is filled with the correct values.
*/

void ComputeLutCol (Image CTg, int L, Weighting *vg, int ivg,
                    int Rmax, LookUpTable Lut)
{
    int x, y, r, r1, r2, ra, rb;

    /* Initializes Lut[ivg] to 0 */
    for (r = 0; r <= Rmax; r++)
        Lut[ivg][r] = 0;

    for (x = 0; x < L - vg->x; x++)
    for (y = 0; y <= x; y++)
    {
        /* Radius of the ball where p1 is located */
        r1 = CTg[ y*L + x ] +1;

        /* Same for p2 */
        r2 = CTg[ (y+vg->y)*L + x+vg->x ] +1;

        if (r1 <= Rmax && r2 > Lut[ivg][r1]) Lut[ivg][r1] = r2;
    }

    rb = 0;
    for (ra = 0; ra <= Rmax; ra++)
    {
        if (Lut[ivg][ra] > rb)
             rb = Lut[ivg][ra];
        else Lut[ivg][ra] = rb;
    }
}


/*
 * Fast extraction of MA points from Bd inter Z^2/8.
 *
 * Input:  x,y the point to test, MgL the generator of the Lut mask,
 *         Lut the look-up table, DTg the distance transform of the section of
 *         the ball, L the side length.
 * Output: returns 1 if point x,y is detected as MA in the DTg.
*/

int IsMAg (int x, int y, MaskG *MgL, LookUpTable Lut, Image DTg, int L)
{
    int xx, yy, val, i;

    val = DTg[ y*L + x ];

    for (i = 0; i < MgL->ng; i++)
    {
        xx = x - MgL->vg[i].x;
        yy = y - MgL->vg[i].y;

        if (0 <= yy && yy <= xx)
        {
            if ( DTg[ yy*L + xx ] >= Lut[i][val] )
                return 0;
        }
    }

    return 1;
}


/*
 * Full Lut Computation Algorithm with determination of MgLut.
 *
 * Input:  CTg and DTg two images, L the side length, MgC the generator of
 *         the dC mask, MgL the generator of the Lut mask, Lut the look-up
 *         table, Rknown the last verified radius, Rtarget the maximum radius
 *         to be verified.
 * Output: MgL and Lut are filled with the correct values.
 *
 * To compute MgL and Lut from beginning to the maximum radius allowed with L:
 *    MgL.ng = 0;
 *    Rknown = 0;
 *    Rtarget = GreatestRadius (&MgC, L);
 *    ComputeAndVerifyLut (CTg, DTg, L, &MgC, &MgL, Lut, Rknown, Rtarget);
 *    Rknown = Rtarget;
*/

void ComputeAndVerifyLut (Image CTg, Image DTg, int L,
                          MaskG *MgC, MaskG *MgL, LookUpTable Lut,
                          int Rknown, int Rtarget)
{
    int x, y, p, R, i;

    ComputeCTg (L, MgC, CTg);

    for (i = 0; i < MgL->ng; i++)
        ComputeLutCol (CTg, L, MgL->vg + i, i, Rtarget, Lut);

    for (R = Rknown+1; R <= Rtarget; R++)
    {
        printf ("R = %d / %d\r", R, Rtarget); fflush (stdout);

        /* Copy Bd^(R) inter Z^2/8 in DTg */
        for (x = 0; x < L; x++)
        for (y = 0; y <= x; y++)
        {
            p = y*L + x;
            if (CTg[p] <= R) DTg[p] = 1; else DTg[p] = 0;
        }

        ComputeDTg (MgC, L, DTg);

        for (x = 1; x < L; x++)
        for (y = 0; y <= x; y++)
        {
            p = y*L + x;
            if ( DTg[p] > 0 && IsMAg (x, y, MgL, Lut, DTg, L) )
            {
                printf ("(%2d, %2d, %3d) added for R = %d\n",
                    x, y, CTg[p], R);

                /* Add a new weighting to MgL */
                i = AddWeighting (MgL, x, y, CTg[p]);

                /* New column in Lut */
                ComputeLutCol (CTg, L, MgL->vg + i, i, Rtarget, Lut);

                if (IsMAg (x, y, MgL, Lut, DTg, L))
                {
                    printf ("\nComputeAndVerifyLut: ERROR for R = %d\n", R);
                    return;
                }
            }
        }
    }
}


/*
 * Compute the greatest verifiable radius of balls.
 *
 * Input:  MgC the generator of a mask, L the side length.
 * Output: returns the greatest verifiable radius in the image.
*/

int GreatestRadius (MaskG *MgC, int L)
{
    int i, res, a = 0, xmax = 0;

    for (i = 0; i < MgC->ng; i++)
    {
        if (MgC->vg[i].x > xmax) xmax = MgC->vg[i].x;

        if (MgC->vg[i].x == 1 && MgC->vg[i].y == 0)
            a = MgC->vg[i].w;
    }
    res = a * (L - xmax) -1;

    if (res >= MAXLUTENTRY) {
        printf ("GreatestRadius: maximum radius limited to MAXLUTENTRY = %d\n",
            MAXLUTENTRY);
        res = MAXLUTENTRY-1;
    }
    return res;
}


/*
 * Print a mask M to a file.
 *
 * Input:  f a file, M the generator of a mask.
 * Output: to file f.
*/

void PrintMask (FILE *f, MaskG *M)
{
    int i;

    for (i = 0; i < M->ng; i++)
        fprintf (f, "(%2d, %2d, %3d)\n",
            M->vg[i].x, M->vg[i].y, M->vg[i].w);
    fprintf (f, "\n");
}


/*
 * Print the LUT to a file.
 *
 * Input:  f a file, MgL the lut mask, Lut the look-up table, CTg the distance
 *         cone to origin, L the side length, Rknown the maximum verified radius.
 * Output: scan the lut entries for possible distance values, test if a line
 *         contains irregularities for the local maximum criteria LMC, and
 *         print them to file f.
*/

void PrintLut (FILE *f, MaskG *MgL, LookUpTable Lut, Image CTg,
               int L, int Rknown)
{
    int x, y, i, j, val, Possible[MAXLUTENTRY];

    /* Mark possible values in Possible[] */
    for (i = 1; i < MAXLUTENTRY; i++)
        Possible[i] = 0;

    for (x = 1; x < L; x++)
    for (y = 0; y <= x; y++)
    {
        val = CTg[ y*L + x ];
        if (val < MAXLUTENTRY) Possible[val] = 1;
    }

    fprintf (f, "Look-Up Table: values which differ from LMC"
                " (L = %d, R <= %d)\n", L, Rknown);

    for (j = 1; j < Rknown; j++)
    if (Possible[j])
    {
        int irreg = 0;

        for (i = 0; i < MgL->ng; i++)
            if (Lut[i][j] != j + MgL->vg[i].w) irreg = 1;

	if (irreg == 1)  /* Show only lines with irregularities */
        {
            fprintf (f, "Lut[][%5d ] = ", j);

            for (i = 0; i < MgL->ng; i++)
            {
                if (Lut[i][j] != j + MgL->vg[i].w)
                     fprintf (f, "%5d ", Lut[i][j]);
                else fprintf (f, "      ");
            }

            fprintf (f, "\n");
        }
    }
    fprintf (f, "\n");
}


/*
 * Compute a file name for the Lut of a mask.
 *
 * Input:  MgC the generator of a mask, name a string, n its size.
 * Output: name contains a file name identifying the mask.
*/

void ComputeLutFileName (MaskG *MgC, char *name, int n)
{
    char word[256];
    int i;

    sprintf (name, "lut2D");
    for (i = 0; i < MgC->ng; i++)
    {
        sprintf (word, "-%d", MgC->vg[i].w);
        if (strlen(name) + strlen(word) < n-8) strcat (name, word);
        else { strcat (name, "-x"); break; }
    }
    strcat (name, ".txt");
}


/*
 * Main program with some examples and custom mode.
*/

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

    MaskG MgC, MgL;
    LookUpTable Lut;
    Image CTg, DTg;
    int L, Rknown, Rtarget, choice;
    char name[256];
    FILE *f;

    printf ("LutChamfer2D - (C) Eric Remy and Edouard Thiel - March 2001\n\n");

    printf ("Enter side length L of image (50 for instance): ");
    scanf ("%d", &L);

    CTg = NewImage (L);
    DTg = NewImage (L);

    printf ("\nList of masks:\n");
    printf ("  1: 2,3   4: 5,7,11        7: 62,88,139,196,224\n");
    printf ("  2: 3,4   5: 14,20,31,44   8: 68,96,152,215,245,280,340,346,413\n");
    printf ("  3: 5,7   6: 72,102,161    0: custom\n");
    printf ("Choose a mask number (4 for instance): ");
    scanf ("%d", &choice);

    MgC.ng = 0;

    switch (choice) {
        case 1:
            AddWeighting (&MgC, 1, 0, 2);
            AddWeighting (&MgC, 1, 1, 3);
            break;
        case 2:
            AddWeighting (&MgC, 1, 0, 3);
            AddWeighting (&MgC, 1, 1, 4);
            break;
        case 3:
            AddWeighting (&MgC, 1, 0, 5);
            AddWeighting (&MgC, 1, 1, 7);
            break;
        case 4:
            AddWeighting (&MgC, 1, 0, 5);
            AddWeighting (&MgC, 1, 1, 7);
            AddWeighting (&MgC, 2, 1, 11);
            break;
        case 5:
            AddWeighting (&MgC, 1, 0, 14);
            AddWeighting (&MgC, 1, 1, 20);
            AddWeighting (&MgC, 2, 1, 31);
            AddWeighting (&MgC, 3, 1, 44);
            break;
        case 6:
            AddWeighting (&MgC, 1, 0, 72);
            AddWeighting (&MgC, 1, 1, 102);
            AddWeighting (&MgC, 2, 1, 161);
            break;
        case 7:
            AddWeighting (&MgC, 1, 0, 62);
            AddWeighting (&MgC, 1, 1, 88);
            AddWeighting (&MgC, 2, 1, 139);
            AddWeighting (&MgC, 3, 1, 196);
            AddWeighting (&MgC, 3, 2, 224);
            break;
        case 8:
            AddWeighting (&MgC, 1, 0, 68);
            AddWeighting (&MgC, 1, 1, 96);
            AddWeighting (&MgC, 2, 1, 152);
            AddWeighting (&MgC, 3, 1, 215);
            AddWeighting (&MgC, 3, 2, 245);
            AddWeighting (&MgC, 4, 1, 280);
            AddWeighting (&MgC, 4, 3, 340);
            AddWeighting (&MgC, 5, 1, 346);
            AddWeighting (&MgC, 6, 1, 413);
            break;
        case 0: {
            int x, y, w;
            printf ("Type '0 0 0' to stop.\n");
            do {
                printf ("Enter 'x y w' : ");
                scanf ("%d %d %d", &x, &y, &w);
                if (x != 0) AddWeighting (&MgC, x, y, w);
            } while (x != 0);
        } break;
    }

    printf ("\nChamfer Mask:\n"); PrintMask (stdout, &MgC);

    /*
     * Computation of Lut and MgL
    */
    printf ("Computed Lut Mask:\n");
    MgL.ng = 0;
    Rknown = 0;
    Rtarget = GreatestRadius (&MgC, L);
    ComputeAndVerifyLut (CTg, DTg, L, &MgC, &MgL, Lut, Rknown, Rtarget);
    Rknown = Rtarget;

    /*
     * Write results into a file
    */
    ComputeLutFileName (&MgC, name, 256);
    printf ("\n\nSaving results in %s ... ", name);
    f = fopen (name, "w");
    if (f == NULL) printf ("error while opening file\n\n");
    else {
        fprintf (f, "Chamfer Mask:\n"); PrintMask (f, &MgC);
        fprintf (f, "Computed Lut Mask:\n"); PrintMask (f, &MgL);
        PrintLut (f, &MgL, Lut, CTg, L, Rknown);
        printf ("ok\n\n");
        fclose (f);
    }

    free (CTg); free (DTg);
    return 0;
}