/*
 * The Npic library
 *
 * Copyright (C) 2003 Edouard Thiel <Edouard.Thiel@lif.univ-mrs.fr>
 *
 * This library is free software under the terms of the 
 * GNU Lesser General Public License (LGPL) version 2.1.
*/

/* DO NOT EDIT !!!  Generated by npic-templa from "mlut_wdtg.ct" */

/*
 * mlut_wdtg.c - 08/01/2009
 *
 * WDT in G(Z^n) for Medial Axis extraction.
 *
 * References:
 *
 * [1] E. Remy and E. Thiel, 2002. Medial axis for chamfer distances: 
 *     computing look-up tables and neighbourhoods in 2D or 3D.
 *     Pattern Recognition Letters, 23(6):649-661.
*/


#include <npic.h>


/*--------------------- P U B L I C - I N T E R F A C E ----------------------*/

/*
 * Fast Cone Distance Transform Algorithm for WD.
 *
 * Input:  nCTg is a L^n image,
 *         nMg is the generator of a weighted distance mask.
 * Output: nCTg is the L^n distance image to the origin.
 *
 * Do nothing if nCTg is not ok. set not ok on error. Verbose.
*/

void NpicMlutCTg_WD (Npic_image *nCTg, Npic_mask *nMg)
{
    if (NpicImageIsOK (nCTg, __func__) != NPIC_SUCCESS)
        return;
    if (NpicImageIsSquare (nCTg, __func__) != NPIC_SUCCESS ||
        NpicMaskIsOK (nMg, __func__) != NPIC_SUCCESS) 
      { nCTg->gen.ok = NPIC_ERROR; return; }
    if (NpicMaskCompat (nMg, nCTg) != NPIC_SUCCESS) {
        nCTg->gen.ok = NpicError (__func__, NPIC_ERR_INCOMPAT, ""); 
        return; 
    }

    NpicFillWhole_i (nCTg, -1);

    switch (nCTg->type) {

        case NPIC_IMAGE_2L : {   /* - - - - - - - - - - - - - - - - - - - - - */
        
            Npic_image_2l *CTg = NpicCastImage(nCTg);
            Npic_mask_2l *Mg = NpicCastMask (nMg);
            int y, x, ya, xa, L = CTg->xmax, i;
            Npic_l min, val;

            CTg->pix[0][0] = 0;

            for (x = 1; x < L; x++)
            for (y = 0; y <= x; y++)
            {
                for (i = 0, min = -1; i < Mg->nb; i++)
                {
                    xa = x - Mg->v[i].x; ya = y - Mg->v[i].y;
                    
                    /* testing xa < L is useless, but who knows ? */
                    if (0 <= ya && ya <= xa && xa < L)
                    {
                        val = CTg->pix[ya][xa] + Mg->v[i].h;
                        if (val < min || min == -1) min = val;
                    }
                }
                CTg->pix[y][x] = min;

                /* Consistency check */
                if (min == -1) {
                    nCTg->gen.ok = NpicError (__func__, NPIC_ERROR, 
                        ": nMg insufficient to compute CTg"); 
                    return; 
                }
            }
          } break;

        case NPIC_IMAGE_3L : {   /* - - - - - - - - - - - - - - - - - - - - - */
        
            Npic_image_3l *CTg = NpicCastImage(nCTg);
            Npic_mask_3l *Mg = NpicCastMask (nMg);
            int z, y, x, za, ya, xa, L = CTg->xmax, i;
            Npic_l min, val;

            CTg->pix[0][0][0] = 0;

            for (x = 1; x < L; x++)
            for (y = 0; y <= x; y++)
            for (z = 0; z <= y; z++)
            {
                for (i = 0, min = -1; i < Mg->nb; i++)
                {
                    xa = x - Mg->v[i].x; ya = y - Mg->v[i].y;
                    za = z - Mg->v[i].z;
                    
                    /* testing xa < L is useless, but who knows ? */
                    if (0 <= za && za <= ya && ya <= xa && xa < L)
                    {
                        val = CTg->pix[za][ya][xa] + Mg->v[i].h;
                        if (val < min || min == -1) min = val;
                    }
                }
                CTg->pix[z][y][x] = min;

                /* Consistency check */
                if (min == -1) {
                    nCTg->gen.ok = NpicError (__func__, NPIC_ERROR, 
                        ": nMg insufficient to compute CTg"); 
                    return; 
                }
            }
          } break;

        case NPIC_IMAGE_4L : {   /* - - - - - - - - - - - - - - - - - - - - - */
        
            Npic_image_4l *CTg = NpicCastImage(nCTg);
            Npic_mask_4l *Mg = NpicCastMask (nMg);
            int t, z, y, x, ta, za, ya, xa, L = CTg->xmax, i;
            Npic_l min, val;

            CTg->pix[0][0][0][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++)
            {
                for (i = 0, min = -1; i < Mg->nb; i++)
                {
                    xa = x - Mg->v[i].x; ya = y - Mg->v[i].y;
                    za = z - Mg->v[i].z; ta = t - Mg->v[i].t;
                    
                    /* testing xa < L is useless, but who knows ? */
                    if (0 <= ta && ta <= za && za <= ya && ya <= xa && xa < L)
                    {
                        val = CTg->pix[ta][za][ya][xa] + Mg->v[i].h;
                        if (val < min || min == -1) min = val;
                    }
                }
                CTg->pix[t][z][y][x] = min;

                /* Consistency check */
                if (min == -1) {
                    nCTg->gen.ok = NpicError (__func__, NPIC_ERROR, 
                        ": nMg insufficient to compute CTg"); 
                    return; 
                }
            }
          } break;

        case NPIC_IMAGE_5L : {   /* - - - - - - - - - - - - - - - - - - - - - */
        
            Npic_image_5l *CTg = NpicCastImage(nCTg);
            Npic_mask_5l *Mg = NpicCastMask (nMg);
            int s, t, z, y, x, sa, ta, za, ya, xa, L = CTg->xmax, i;
            Npic_l min, val;

            CTg->pix[0][0][0][0][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++)
            for (s = 0; s <= t; s++)
            {
                for (i = 0, min = -1; i < Mg->nb; i++)
                {
                    xa = x - Mg->v[i].x; ya = y - Mg->v[i].y;
                    za = z - Mg->v[i].z; ta = t - Mg->v[i].t;
                    sa = s - Mg->v[i].s;
                    
                    /* testing xa < L is useless, but who knows ? */
                    if (0 <= sa && sa <= ta && ta <= za && za <= ya && ya <= xa && xa < L)
                    {
                        val = CTg->pix[sa][ta][za][ya][xa] + Mg->v[i].h;
                        if (val < min || min == -1) min = val;
                    }
                }
                CTg->pix[s][t][z][y][x] = min;

                /* Consistency check */
                if (min == -1) {
                    nCTg->gen.ok = NpicError (__func__, NPIC_ERROR, 
                        ": nMg insufficient to compute CTg"); 
                    return; 
                }
            }
          } break;

        case NPIC_IMAGE_6L : {   /* - - - - - - - - - - - - - - - - - - - - - */
        
            Npic_image_6l *CTg = NpicCastImage(nCTg);
            Npic_mask_6l *Mg = NpicCastMask (nMg);
            int r, s, t, z, y, x, ra, sa, ta, za, ya, xa, L = CTg->xmax, i;
            Npic_l min, val;

            CTg->pix[0][0][0][0][0][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++)
            for (s = 0; s <= t; s++)
            for (r = 0; r <= s; r++)
            {
                for (i = 0, min = -1; i < Mg->nb; i++)
                {
                    xa = x - Mg->v[i].x; ya = y - Mg->v[i].y;
                    za = z - Mg->v[i].z; ta = t - Mg->v[i].t;
                    sa = s - Mg->v[i].s; ra = r - Mg->v[i].r;
                    
                    /* testing xa < L is useless, but who knows ? */
                    if (0 <= ra && ra <= sa && sa <= ta && ta <= za && za <= ya && ya <= xa && xa < L)
                    {
                        val = CTg->pix[ra][sa][ta][za][ya][xa] + Mg->v[i].h;
                        if (val < min || min == -1) min = val;
                    }
                }
                CTg->pix[r][s][t][z][y][x] = min;

                /* Consistency check */
                if (min == -1) {
                    nCTg->gen.ok = NpicError (__func__, NPIC_ERROR, 
                        ": nMg insufficient to compute CTg"); 
                    return; 
                }
            }
          } break;

        default : nCTg->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
    }
}


/*
 * Fast Distance Transform in G(Z^n).
 *
 * Input:  nCTg the distance cone to the origin, nDTg init to 0
 *         R the radius of the ball in CTg.
 * Output: nDTg the distance image to the background of the ball.
 * Assume that nCTg and nDTg are L^n images.
 *
 * Do nothing if nCTg or nDTg is not ok. set not ok for nDTg on error. Verbose.
 * Return computation time in second.microsecond
*/

double NpicMlutDTg_WD (Npic_image *nCTg, Npic_image *nDTg, Npic_mask *nMg, int R)
{
    double time1, time2;

    if (NpicImageIsOK_DS1 (nDTg, nCTg, __func__) != NPIC_SUCCESS)
        return 0;
    if (NpicImageIsSquare (nCTg, __func__) != NPIC_SUCCESS ||
        NpicMaskIsOK (nMg, __func__) != NPIC_SUCCESS) 
      { nDTg->gen.ok = NPIC_ERROR; return 0; }
    if (NpicMaskCompat (nMg, nCTg) != NPIC_SUCCESS ||
        NpicSameImage (nDTg, nCTg, NPIC_TYPE | NPIC_SIZE) != NPIC_SUCCESS) {
        nDTg->gen.ok = NpicError (__func__, NPIC_ERR_INCOMPAT, ""); 
        return 0; 
    }

    time1 = NpicGetTime();

    /*
     * In all the loops we must have:  0 <= r <= s <= t <= z <= y <= x <= xM < L 
    */
    switch (nCTg->type) {

        case NPIC_IMAGE_2L : {   /* - - - - - - - - - - - - - - - - - - - - - */
        
            Npic_image_2l *CTg = NpicCastImage(nCTg),
                          *DTg = NpicCastImage(nDTg);
            Npic_mask_2l *Mg = NpicCastMask (nMg);
            int y, x, ya, xa, L = CTg->xmax, xM, i;
            Npic_l min, val;

            /* compute bound xM[ and verify that xM < L */
            for (xM = 0; xM < L; xM++) 
                if (CTg->pix[0][xM] > R) break;
            if (xM == L) {
                DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large");
                return 0;
            }

            for (x = xM; x >= 0; x--)
            for (y = x ; y >= 0; y--)
            if (CTg->pix[y][x] <= R)
            {
                for (i = 0, min = -1; i < Mg->nb; i++)
                {
                    xa = x + Mg->v[i].x; ya = y + Mg->v[i].y;
        
                    /* Testing xx < L shoud be sufficient, but who knows ? */
                    if (0 <= ya && ya <= xa && xa < L)
                    {
                        val = DTg->pix[ya][xa] + Mg->v[i].h;
                        if (val < min || min == -1) min = val;
                    }
                }
                DTg->pix[y][x] = min;

                /* Consistency check */
                if (min == -1) {
                    nDTg->gen.ok = NpicError (__func__, NPIC_ERROR, 
                        ": nMg insufficient to compute DTg"); 
                    return 0; 
                }
            }
          } break;

        case NPIC_IMAGE_3L : {   /* - - - - - - - - - - - - - - - - - - - - - */
        
            Npic_image_3l *CTg = NpicCastImage(nCTg),
                          *DTg = NpicCastImage(nDTg);
            Npic_mask_3l *Mg = NpicCastMask (nMg);
            int z, y, x, za, ya, xa, L = CTg->xmax, xM, i;
            Npic_l min, val;

            /* compute bound xM[ and verify that xM < L */
            for (xM = 0; xM < L; xM++) 
                if (CTg->pix[0][0][xM] > R) break;
            if (xM == L) {
                DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large");
                return 0;
            }

            for (x = xM; x >= 0; x--)
            for (y = x ; y >= 0; y--)
            for (z = y ; z >= 0; z--)
            if (CTg->pix[z][y][x] <= R)
            {
                for (i = 0, min = -1; i < Mg->nb; i++)
                {
                    xa = x + Mg->v[i].x; ya = y + Mg->v[i].y;
                    za = z + Mg->v[i].z;
        
                    /* Testing xx < L shoud be sufficient, but who knows ? */
                    if (0 <= za && za <= ya && ya <= xa && xa < L)
                    {
                        val = DTg->pix[za][ya][xa] + Mg->v[i].h;
                        if (val < min || min == -1) min = val;
                    }
                }
                DTg->pix[z][y][x] = min;

                /* Consistency check */
                if (min == -1) {
                    nDTg->gen.ok = NpicError (__func__, NPIC_ERROR, 
                        ": nMg insufficient to compute DTg"); 
                    return 0; 
                }
            }
          } break;

        case NPIC_IMAGE_4L : {   /* - - - - - - - - - - - - - - - - - - - - - */
        
            Npic_image_4l *CTg = NpicCastImage(nCTg),
                          *DTg = NpicCastImage(nDTg);
            Npic_mask_4l *Mg = NpicCastMask (nMg);
            int t, z, y, x, ta, za, ya, xa, L = CTg->xmax, xM, i;
            Npic_l min, val;

            /* compute bound xM[ and verify that xM < L */
            for (xM = 0; xM < L; xM++) 
                if (CTg->pix[0][0][0][xM] > R) break;
            if (xM == L) {
                DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large");
                return 0;
            }

            for (x = xM; x >= 0; x--)
            for (y = x ; y >= 0; y--)
            for (z = y ; z >= 0; z--)
            for (t = z ; t >= 0; t--)
            if (CTg->pix[t][z][y][x] <= R)
            {
                for (i = 0, min = -1; i < Mg->nb; i++)
                {
                    xa = x + Mg->v[i].x; ya = y + Mg->v[i].y;
                    za = z + Mg->v[i].z; ta = t + Mg->v[i].t;
        
                    /* Testing xx < L shoud be sufficient, but who knows ? */
                    if (0 <= ta && ta <= za && za <= ya && ya <= xa && xa < L)
                    {
                        val = DTg->pix[ta][za][ya][xa] + Mg->v[i].h;
                        if (val < min || min == -1) min = val;
                    }
                }
                DTg->pix[t][z][y][x] = min;

                /* Consistency check */
                if (min == -1) {
                    nDTg->gen.ok = NpicError (__func__, NPIC_ERROR, 
                        ": nMg insufficient to compute DTg"); 
                    return 0; 
                }
            }
          } break;

        case NPIC_IMAGE_5L : {   /* - - - - - - - - - - - - - - - - - - - - - */
        
            Npic_image_5l *CTg = NpicCastImage(nCTg),
                          *DTg = NpicCastImage(nDTg);
            Npic_mask_5l *Mg = NpicCastMask (nMg);
            int s, t, z, y, x, sa, ta, za, ya, xa, L = CTg->xmax, xM, i;
            Npic_l min, val;

            /* compute bound xM[ and verify that xM < L */
            for (xM = 0; xM < L; xM++) 
                if (CTg->pix[0][0][0][0][xM] > R) break;
            if (xM == L) {
                DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large");
                return 0;
            }

            for (x = xM; x >= 0; x--)
            for (y = x ; y >= 0; y--)
            for (z = y ; z >= 0; z--)
            for (t = z ; t >= 0; t--)
            for (s = t ; s >= 0; s--)
            if (CTg->pix[s][t][z][y][x] <= R)
            {
                for (i = 0, min = -1; i < Mg->nb; i++)
                {
                    xa = x + Mg->v[i].x; ya = y + Mg->v[i].y;
                    za = z + Mg->v[i].z; ta = t + Mg->v[i].t;
                    sa = s + Mg->v[i].s;
        
                    /* Testing xx < L shoud be sufficient, but who knows ? */
                    if (0 <= sa && sa <= ta && ta <= za && za <= ya && ya <= xa && xa < L)
                    {
                        val = DTg->pix[sa][ta][za][ya][xa] + Mg->v[i].h;
                        if (val < min || min == -1) min = val;
                    }
                }
                DTg->pix[s][t][z][y][x] = min;

                /* Consistency check */
                if (min == -1) {
                    nDTg->gen.ok = NpicError (__func__, NPIC_ERROR, 
                        ": nMg insufficient to compute DTg"); 
                    return 0; 
                }
            }
          } break;

        case NPIC_IMAGE_6L : {   /* - - - - - - - - - - - - - - - - - - - - - */
        
            Npic_image_6l *CTg = NpicCastImage(nCTg),
                          *DTg = NpicCastImage(nDTg);
            Npic_mask_6l *Mg = NpicCastMask (nMg);
            int r, s, t, z, y, x, ra, sa, ta, za, ya, xa, L = CTg->xmax, xM, i;
            Npic_l min, val;

            /* compute bound xM[ and verify that xM < L */
            for (xM = 0; xM < L; xM++) 
                if (CTg->pix[0][0][0][0][0][xM] > R) break;
            if (xM == L) {
                DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large");
                return 0;
            }

            for (x = xM; x >= 0; x--)
            for (y = x ; y >= 0; y--)
            for (z = y ; z >= 0; z--)
            for (t = z ; t >= 0; t--)
            for (s = t ; s >= 0; s--)
            for (r = s ; r >= 0; r--)
            if (CTg->pix[r][s][t][z][y][x] <= R)
            {
                for (i = 0, min = -1; i < Mg->nb; i++)
                {
                    xa = x + Mg->v[i].x; ya = y + Mg->v[i].y;
                    za = z + Mg->v[i].z; ta = t + Mg->v[i].t;
                    sa = s + Mg->v[i].s; ra = r + Mg->v[i].r;
        
                    /* Testing xx < L shoud be sufficient, but who knows ? */
                    if (0 <= ra && ra <= sa && sa <= ta && ta <= za && za <= ya && ya <= xa && xa < L)
                    {
                        val = DTg->pix[ra][sa][ta][za][ya][xa] + Mg->v[i].h;
                        if (val < min || min == -1) min = val;
                    }
                }
                DTg->pix[r][s][t][z][y][x] = min;

                /* Consistency check */
                if (min == -1) {
                    nDTg->gen.ok = NpicError (__func__, NPIC_ERROR, 
                        ": nMg insufficient to compute DTg"); 
                    return 0; 
                }
            }
          } break;

        default : nDTg->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
    }

    time2 = NpicGetTime();
    return time2 - time1;
}


/*-------------------- P R I V A T E - F U N C T I O N S ---------------------*/


/*----------------------------------------------------------------------------*/