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


#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, z, t, 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];

#define MAXCOORDS 1000

typedef struct {
    int X[MAXCOORDS], Y[MAXCOORDS], Z[MAXCOORDS];
} HyperCoords;

#define COORD(H,x,y,z,t) ((H)->X[x] + (H)->Y[y] + (H)->Z[z] + t)


/*
 * Pre-compute coordinates in hypercone to save memory usage :
 * need 24 times less memory than a hypervolume.
 *
 * Input:  L the side length of the hypercone.
 * Output: H the pre-computed coordinates in the hypercone.
 *
 * To acces a (x,y,z,t) point in a hypercone Cone :
 *      Cone[ H->X[x] + H->Y[y] + H->Z[z] + t ]
 *  or  Cone[ COORD(H,x,y,z,t) ]
*/

void ComputeHyperCoords (int L, HyperCoords *H)
{
    int i;

    for (i = 0; i < L; i++)
    {
        H->X[i] = i*(i+1)*(i+2)*(i+3)/24;
        H->Y[i] = i*(i+1)*(i+2)/6;
        H->Z[i] = i*(i+1)/2;
    }
}


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

Image NewImage (int L)
{
    int n = L*(L+1)*(L+2)*(L+3)/24;
    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,z,t,w) the weighting.
 * Output: the number of the weighting in the mask.
*/

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

    if (i >= MAXWEIGHTING) {
        printf ("AddWeighting: weighting number limited to MAXWEIGHTING = %d\n",
            MAXWEIGHTING);
        return 0;
    }
    if (! (0 <= t && t <= z && z <= y && y <= x && 0 < x && 0 < w)) {
        printf ("AddWeighting: (x = %d, y = %d, z = %d, t = %d, w = %d)\n",
            x, y, z, t, w);
        printf ("  does not respect 0 <= t <= z <= y <= x, 0 < x, 0 < w\n");
        return 0;
    }
    M->vg[i].x = x;
    M->vg[i].y = y;
    M->vg[i].z = z;
    M->vg[i].t = t;
    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^4 distance image to the origin.
*/

void ComputeCTg (int L, MaskG *MgC, Image CTg, HyperCoords *H)
{
    int x, y, z, t, xx, yy, zz, tt, i, min, val;

    CTg[0] = 0;

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

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

            if (0 <= tt && tt <= zz && zz <= yy && yy <= xx)
            {
                val = CTg[ COORD(H,xx,yy,zz,tt) ] + MgC->vg[i].w;

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

        CTg[ COORD(H,x,y,z,t) ] = min;
    }
}


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

void ComputeDTg (MaskG *MgC, int L, Image DTg, HyperCoords *H)
{
    int x, y, z, t, xx, yy, zz, tt, min, val, i;

    for (t = L-1; t >= 0; t--)
    for (z = L-1; z >= t; z--)
    for (y = L-1; y >= z; y--)
    for (x = L-1; x >= y; x--)
    {
        if (DTg[ COORD(H,x,y,z,t) ] != 0)
        {
            min = INFINITY;

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

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

            DTg[ COORD(H,x,y,z,t) ] = 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, HyperCoords *H)
{
    int x, y, z, t, 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++)
    for (z = 0; z <= y; z++)
    for (t = 0; t <= z; t++)
    {
        /* Radius of the ball where p1 is located */
        r1 = CTg[ COORD(H,x,y,z,t) ] +1;

        /* Same for p2 */
        r2 = CTg[ COORD(H,x+vg->x,y+vg->y,z+vg->z,t+vg->t) ] +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^4/384.
 *
 * Input:  x,y,z,t 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,z,t is detected as MA in the DTg.
*/

int IsMAg (int x, int y, int z, int t, MaskG *MgL, LookUpTable Lut,
           Image DTg, int L, HyperCoords *H)
{
    int xx, yy, zz, tt, val, i;

    val = DTg[ COORD(H,x,y,z,t) ];

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

        if (0 <= tt && tt <= zz && zz <= yy && yy <= xx)
        {
            if ( DTg[ COORD(H,xx,yy,zz,tt) ] >= 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, HyperCoords *H)
{
    int x, y, z, t, p, R, i;

    ComputeCTg (L, MgC, CTg, H);

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

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

        /* Copy Bd^(R) inter Z^4/384 in DTg */
        for (x = 0; x < L; x++)
        for (y = 0; y <= x; y++)
        for (z = 0; z <= y; z++)
        for (t = 0; t <= z; t++)
        {
            p = COORD(H,x,y,z,t);
            if (CTg[p] <= R) DTg[p] = 1; else DTg[p] = 0;
        }

        ComputeDTg (MgC, L, DTg, H);

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

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

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

                if (IsMAg (x, y, z, t, MgL, Lut, DTg, L, H))
                {
                    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 &&
            MgC->vg[i].z == 0 && MgC->vg[i].t == 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, %2d, %2d, %3d)\n",
            M->vg[i].x, M->vg[i].y, M->vg[i].z, M->vg[i].t, 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, HyperCoords *H)
{
    int x, y, z, t, 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++)
    for (z = 0; z <= y; z++)
    for (t = 0; t <= z; t++)
    {
        val = CTg[ COORD(H,x,y,z,t) ];
        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, "lut4D");
    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;
    HyperCoords H;

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

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

    ComputeHyperCoords (L, &H);
    CTg = NewImage (L);
    DTg = NewImage (L);

    printf ("\nList of masks:\n");
    printf ("  1: 3,4,5,6   2: 7,10,12,13   0: custom\n");
    printf ("Choose a mask number (1 for instance): ");
    scanf ("%d", &choice);

    MgC.ng = 0;

    switch (choice) {
        case 1:
            AddWeighting (&MgC, 1, 0, 0, 0, 3);
            AddWeighting (&MgC, 1, 1, 0, 0, 4);
            AddWeighting (&MgC, 1, 1, 1, 0, 5);
            AddWeighting (&MgC, 1, 1, 1, 1, 6);
            break;
        case 2:
            AddWeighting (&MgC, 1, 0, 0, 0, 7);
            AddWeighting (&MgC, 1, 1, 0, 0, 10);
            AddWeighting (&MgC, 1, 1, 1, 0, 12);
            AddWeighting (&MgC, 1, 1, 1, 1, 13);
            break;
        case 0: {
            int x, y, z, t, w;
            printf ("Type '0 0 0 0 0' to stop.\n");
            do {
                printf ("Enter 'x y z t w' : ");
                scanf ("%d %d %d %d %d", &x, &y, &z, &t, &w);
                if (x != 0) AddWeighting (&MgC, x, y, z, t, 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, &H);
    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, &H);
        printf ("ok\n\n");
        fclose (f);
    }

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