/*
 * 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_ma.ct" */

/*
 * mlut_ma.c - 07/01/2008
 *
 * Computation of Medial Axis on Distance Transforms using Look-Up Tables.
 *
 * The general LUT method for extracting the Medial Axis from a DT is presented
 * in [1,2]. The Mlut and LUT computation is shown in [3,4].
 *
 * References:
 *
 * [1] G. Borgefors, I. Ragnemalm, and G. Sanniti di Baja, 1991. 
 *     "The Euclidean Distance Transform : finding the local maxima and 
 *      reconstructing the shape". In 7th Scandinavian Conf. on Image Analysis, 
 *      volume 2, pages 974-981, Aalborg, Denmark.
 *
 * [2] G. Borgefors, 1993. "Centres of maximal disks in the 5-7-11 distance 
 *     transform". In 8th Scandinavian Conf. on Image Analysis, pages 105-111, 
 *     Tromso, Norway.
 * 
 * [3] E. Remy and E. Thiel, 2005. Exact Medial Axis with Euclidean Distance. 
 *     Image and Vision Computing, 23(2):167-175.
 *
 * [4] 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 ----------------------*/

/*
 * Compute a DT for a binary object, then extract Medial Axis from this DT.
 *
 * Input: np1 a binary image, nMdist a distance mask (half weighted or Euclidean),
 *        nMlut a MA test neighboorhood or NULL.
 * Output: np1 contains the DT, nMA contains the medial axis,
 *         nMlut (if it wasn't NULL) is computed or updated.
 *
 * Do nothing if images are not ok. Set not ok for np1 on error. Verbose.
*/

void NpicExtractMA_bin (Npic_image *np1, Npic_image *nMA,
        Npic_mask *nMdist, Npic_mask *nMlut)
{
    double t1;

    if (NpicImageIsOK (np1, __func__) != NPIC_SUCCESS) return;

    printf ("Computing DT ... \n");
    t1 = NpicDT (np1, nMdist);
    if (np1->gen.ok != NPIC_SUCCESS) return;
    printf ("DT computed in %.3lf s\n", t1);

    NpicExtractMA_DT (np1, nMA, nMdist, nMlut);
}


/*
 * Extract the Medial Axis from a DT.
 *
 * Input: nDT a DT, nMdist a distance mask (half weighted or Euclidean),
 *        nMlut a MA test neighboorhood or NULL.
 * Output: nMA contains the medial axis,
 *         nMlut (if it wasn't NULL) is computed or updated.
 *
 * Do nothing if images are not ok. Set not ok for nMA on error. Verbose.
*/

void NpicExtractMA_DT (Npic_image *nDT, Npic_image *nMA,
        Npic_mask *nMdist, Npic_mask *nMlut)
{
    int Mlut_tmp = 0, dtype, L = -1;
    Npic_l Rmax = -1, Rverif = 0;
    Npic_mask *nMg = NULL;
    Npic_image *nCTg = NULL;
    Npic_Lut Lut;
    double time1, time2;

    if (NpicImageIsOK_DS1 (nMA, nDT, __func__) != NPIC_SUCCESS)
        return;
    if (NpicMaskIsOK (nMdist, __func__) != NPIC_SUCCESS) 
        { nMA->gen.ok = NPIC_ERROR; return; }
    if (NpicMaskCompat (nMdist, nDT) != NPIC_SUCCESS ||
        NpicSameImage (nDT, nMA, NPIC_TYPE | NPIC_SIZE) != NPIC_SUCCESS) {
        nMA->gen.ok = NpicError (__func__, NPIC_ERR_INCOMPAT, ""); 
        return; 
    }

    printf ("Searching max value in DT ... "); fflush (stdout);
    Rmax = NpicMaxValue (nDT);
    printf ("%"NPIC_PL"\n", Rmax);
    NpicMlutInitLut (&Lut, Rmax);
    /* Now from this point, jump to lb_final on error */
  
    dtype = NpicDTFindDistType (nMdist);
    if (dtype == NPIC_DIST_WEIGHTED) {
        printf ("Number of vectors in distance half mask: %d\n", nMdist->gen.nb);
        nMg = NpicCompGeneMask (nMdist);
        if (nMg == NULL)  { nMA->gen.ok = NPIC_ERROR; goto lb_final; }
        printf ("Number of vectors in generator: %d\n", nMg->gen.nb); 
    }
 
    if (nMlut == NULL) {  
        Mlut_tmp = 1; /* we need a temporary nMlut */
        nMlut = NpicCreateMaskDP (nDT->gen.dim, NPIC_L);
    }
    if (NpicMaskIsOK (nMlut, __func__) != NPIC_SUCCESS) 
        { nMA->gen.ok = NPIC_ERROR; goto lb_final; }

    printf ("Number of vectors in Mlut: %d\n", nMlut->gen.nb);
    Rverif = atoi (NpicMaskGetPropD (nMlut, "RVerified", "0"));
    printf ("RVerified = %" NPIC_PL "\n", Rverif);

    switch (dtype) {
        case NPIC_DIST_WEIGHTED       : L = NpicMlutComp_L_WD (nMg, Rmax) + 5; break;
        case NPIC_DIST_EUCLID_SQUARED : L = sqrt (Rmax) + 5; break;
        default : 
            nMA->gen.ok = NpicError (__func__, NPIC_ERR_UNAV_DIST, ""); 
            goto lb_final;
    }

    printf ("Creating a %d^%d CTg image ...\n", L, nDT->gen.dim);
    switch (nDT->gen.dim) {
        case 2 : nCTg = NpicCreateImage_2l (L, L, 0, 0); break;
        case 3 : nCTg = NpicCreateImage_3l (L, L, L, 0, 0, 0); break;
        case 4 : nCTg = NpicCreateImage_4l (L, L, L, L, 0, 0, 0, 0); break;
        case 5 : nCTg = NpicCreateImage_5l (L, L, L, L, L, 0, 0, 0, 0, 0); break;
        case 6 : nCTg = NpicCreateImage_6l (L, L, L, L, L, L, 0, 0, 0, 0, 0, 0); break;
    }
    if (nCTg == NULL) { nMA->gen.ok = NPIC_ERROR; goto lb_final; }

    printf ("Computing CTg ...\n");
    time1 = NpicGetTime();
    switch (dtype) {
        case NPIC_DIST_WEIGHTED       : NpicMlutCTg_WD  (nCTg, nMg); break;
        case NPIC_DIST_EUCLID_SQUARED : NpicMlutCTg_SED (nCTg);      break;
        default : 
            nMA->gen.ok = NpicError (__func__, NPIC_ERR_UNAV_DIST, ""); 
            goto lb_final;
    }
    time2 = NpicGetTime();
    printf ("CTg computed in %.3lf s\n", time2 - time1);

    /* Enlarge Mlut if necessary */
    if (Rmax > Rverif) {
        printf ("Computing MLut for radii ]%"NPIC_PL" .. %"NPIC_PL"] ...\n", Rverif, Rmax);
        if (NpicMlutCompLutMask (nCTg, nMlut, &Lut, Rverif, Rmax, nMg) != NPIC_SUCCESS) {
            nMA->gen.ok = NPIC_ERROR; goto lb_final;
        } else {
            char s[200];
            sprintf (s, "%" NPIC_PL, Rmax);
            NpicMaskAddProp (nMlut, "RVerified", s);
       }
    }

    /* Compute Medial Axis */
    time1 = NpicGetTime();
    if (Rmax > Rverif) {
        printf ("Computing Medial Axis, with Lut already in memory ... \n"); 
        NpicExtractMA_Lut (nDT, nMA, nMlut, &Lut, Rmax);
    } else {
        printf ("Computing Medial Axis, one Lut col at a time in memory ... \n"); 
        NpicExtractMA_Mlut (nDT, nMA, nMlut, nCTg, Rmax);
    }
    time2 = NpicGetTime();
    printf ("MA computed in %.3lf s\n", time2 - time1);

lb_final:
    NpicMlutFreeLut (&Lut);  /* don't jump to lb_final before Lut has been initialized! */
    if (nCTg != NULL) NpicDestroyImage (nCTg);
    if (nMg != NULL) NpicDestroyMask (nMg);
    if (Mlut_tmp) NpicDestroyMask (nMlut);
}


/*
 * Extract Medial Axis from a DT having Lut completely pre-computed.
 *
 * Do nothing if nDT, nMA or nMlut is not ok. set not ok for nMA on error. Verbose.
*/

void NpicExtractMA_Lut (Npic_image *nDT, Npic_image *nMA,
        Npic_mask *nMlut, Npic_Lut *Lut, Npic_l Rmax)
{
    int ncol, k;

    if (NpicImageIsOK_DS1 (nMA, nDT, __func__) != NPIC_SUCCESS)
        return;
    if (NpicMaskIsOK (nMlut, __func__) != NPIC_SUCCESS) 
      { nMA->gen.ok = NPIC_ERROR; return; }
    if (NpicMaskCompat (nMlut, nDT) != NPIC_SUCCESS ||
        NpicSameImage (nMA, nDT, NPIC_TYPE | NPIC_SIZE) != NPIC_SUCCESS) {
        nMA->gen.ok = NpicError (__func__, NPIC_ERR_INCOMPAT, ""); 
        return; 
    }

    switch (nDT->type) {

        case NPIC_IMAGE_2L : {
            Npic_image_2l *DT = NpicCastImage (nDT),
                          *MA = NpicCastImage (nMA);
            Npic_mask_2l *Mlut = NpicCastMask (nMlut);
            int y, x;

            for (ncol = 0; ncol < Mlut->nb; ncol++)
                if (Mlut->v[ncol].h > Rmax) break;

            for (y = 0; y < DT->ymax; y++)
            for (x = 0; x < DT->xmax; x++)
            {
                k = NpicMlutIsMAlut_2l (DT, Mlut, Lut, ncol, y, x);
                MA->pix[y][x] = k ? DT->pix[y][x] : 0;
            }
          } break;

        case NPIC_IMAGE_3L : {
            Npic_image_3l *DT = NpicCastImage (nDT),
                          *MA = NpicCastImage (nMA);
            Npic_mask_3l *Mlut = NpicCastMask (nMlut);
            int z, y, x;

            for (ncol = 0; ncol < Mlut->nb; ncol++)
                if (Mlut->v[ncol].h > Rmax) break;

            for (z = 0; z < DT->zmax; z++)
            {
                printf ("z = %d /%d      \r", z, DT->zmax); fflush (stdout);
                for (y = 0; y < DT->ymax; y++)
                for (x = 0; x < DT->xmax; x++)
                {
                    k = NpicMlutIsMAlut_3l (DT, Mlut, Lut, ncol, z, y, x);
                    MA->pix[z][y][x] = k ? DT->pix[z][y][x] : 0;
                }
            }
            printf ("\n");
          } break;

        case NPIC_IMAGE_4L : {
            Npic_image_4l *DT = NpicCastImage (nDT),
                          *MA = NpicCastImage (nMA);
            Npic_mask_4l *Mlut = NpicCastMask (nMlut);
            int t, z, y, x;

            for (ncol = 0; ncol < Mlut->nb; ncol++)
                if (Mlut->v[ncol].h > Rmax) break;

            for (t = 0; t < DT->tmax; t++)
            for (z = 0; z < DT->zmax; z++)
            {
                printf ("t = %d /%d  ", t, DT->tmax);
                printf ("z = %d /%d      \r", z, DT->zmax); fflush (stdout);
                for (y = 0; y < DT->ymax; y++)
                for (x = 0; x < DT->xmax; x++)
                {
                    k = NpicMlutIsMAlut_4l (DT, Mlut, Lut, ncol, t, z, y, x);
                    MA->pix[t][z][y][x] = k ? DT->pix[t][z][y][x] : 0;
                }
            }
            printf ("\n");
          } break;

        case NPIC_IMAGE_5L : {
            Npic_image_5l *DT = NpicCastImage (nDT),
                          *MA = NpicCastImage (nMA);
            Npic_mask_5l *Mlut = NpicCastMask (nMlut);
            int s, t, z, y, x;

            for (ncol = 0; ncol < Mlut->nb; ncol++)
                if (Mlut->v[ncol].h > Rmax) break;

            for (s = 0; s < DT->smax; s++)
            for (t = 0; t < DT->tmax; t++)
            for (z = 0; z < DT->zmax; z++)
            {
                printf ("s = %d /%d  ", s, DT->smax);
                printf ("t = %d /%d  ", t, DT->tmax);
                printf ("z = %d /%d      \r", z, DT->zmax); fflush (stdout);
                for (y = 0; y < DT->ymax; y++)
                for (x = 0; x < DT->xmax; x++)
                {
                    k = NpicMlutIsMAlut_5l (DT, Mlut, Lut, ncol, s, t, z, y, x);
                    MA->pix[s][t][z][y][x] = k ? DT->pix[s][t][z][y][x] : 0;
                }
            }
            printf ("\n");
          } break;

        case NPIC_IMAGE_6L : {
            Npic_image_6l *DT = NpicCastImage (nDT),
                          *MA = NpicCastImage (nMA);
            Npic_mask_6l *Mlut = NpicCastMask (nMlut);
            int r, s, t, z, y, x;

            for (ncol = 0; ncol < Mlut->nb; ncol++)
                if (Mlut->v[ncol].h > Rmax) break;

            for (r = 0; r < DT->rmax; r++)
            for (s = 0; s < DT->smax; s++)
            for (t = 0; t < DT->tmax; t++)
            for (z = 0; z < DT->zmax; z++)
            {
                printf ("r = %d /%d  ", r, DT->rmax);
                printf ("s = %d /%d  ", s, DT->smax);
                printf ("t = %d /%d  ", t, DT->tmax);
                printf ("z = %d /%d      \r", z, DT->zmax); fflush (stdout);
                for (y = 0; y < DT->ymax; y++)
                for (x = 0; x < DT->xmax; x++)
                {
                    k = NpicMlutIsMAlut_6l (DT, Mlut, Lut, ncol, r, s, t, z, y, x);
                    MA->pix[r][s][t][z][y][x] = k ? DT->pix[r][s][t][z][y][x] : 0;
                }
            }
            printf ("\n");
          } break;

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


/*
 * Extract Medial Axis from a DT having Mlut pre-computed;
 * compute and store one Lut column at a time.
 *
 * Do nothing if nDT, nMA, nCTg or nMlut is not ok. set not ok for nMA on error. Verbose.
*/

void NpicExtractMA_Mlut (Npic_image *nDT, Npic_image *nMA, Npic_mask *nMlut, 
        Npic_image *nCTg, Npic_l Rmax)
{
    int ncol, k, i;
    Npic_l *col;
    Npic_vec vg;

    if (NpicImageIsOK_DS2 (nMA, nDT, nCTg, __func__) != NPIC_SUCCESS)
        return;
    if (NpicMaskIsOK (nMlut, __func__) != NPIC_SUCCESS)
      { nMA->gen.ok = NPIC_ERROR; return; }
    if (NpicMaskCompat (nMlut, nDT) != NPIC_SUCCESS ||
        NpicSameImage (nDT, nCTg, NPIC_TYPE) != NPIC_SUCCESS ||
        NpicSameImage (nMA, nDT, NPIC_TYPE | NPIC_SIZE) != NPIC_SUCCESS) {
        nMA->gen.ok = NpicError (__func__, NPIC_ERR_INCOMPAT, ""); 
        return; 
    }

    col = malloc (sizeof(Npic_l) * (Rmax+1));
    if (col == NULL) {
        nMA->gen.ok = NpicError (__func__, NPIC_ERR_MALLOC, ""); 
        return; 
    }

    NpicCopyImageI (nMA, nDT);

    switch (nDT->type) {

        case NPIC_IMAGE_2L : {
            Npic_image_2l *DT = NpicCastImage (nDT),
                          *MA = NpicCastImage (nMA);
            Npic_mask_2l *Mlut = NpicCastMask (nMlut);
            int y, x;

            for (ncol = 0; ncol < Mlut->nb; ncol++)
                    if (Mlut->v[ncol].h > Rmax) break;

            for (i = 0; i < ncol; i++)
            {
                printf ("Col %d /%d \r", i+1, ncol); fflush (stdout);

                vg.x = Mlut->v[i].x; vg.y = Mlut->v[i].y;
                NpicMlutCompLutCol (nCTg, &vg, Rmax, col);

                for (y = 0; y < DT->ymax; y++)
                for (x = 0; x < DT->xmax; x++)
                if (MA->pix[y][x] != 0)
                {
                    k = NpicMlutIsMAcol_2l (DT, Mlut, col, Rmax, &vg, y, x);
                    if (k == 0) MA->pix[y][x] = 0;
                }
            }
            printf ("\n");
          } break;

        case NPIC_IMAGE_3L : {
            Npic_image_3l *DT = NpicCastImage (nDT),
                          *MA = NpicCastImage (nMA);
            Npic_mask_3l *Mlut = NpicCastMask (nMlut);
            int z, y, x;

            for (ncol = 0; ncol < Mlut->nb; ncol++)
                    if (Mlut->v[ncol].h > Rmax) break;

            for (i = 0; i < ncol; i++)
            {
                printf ("Col %d /%d \r", i+1, ncol); fflush (stdout);

                vg.x = Mlut->v[i].x; vg.y = Mlut->v[i].y;
                vg.z = Mlut->v[i].z; 
                NpicMlutCompLutCol (nCTg, &vg, Rmax, col);

                for (z = 0; z < DT->zmax; z++)
                for (y = 0; y < DT->ymax; y++)
                for (x = 0; x < DT->xmax; x++)
                if (MA->pix[z][y][x] != 0)
                {
                    k = NpicMlutIsMAcol_3l (DT, Mlut, col, Rmax, &vg, z, y, x);
                    if (k == 0) MA->pix[z][y][x] = 0;
                }
            }
            printf ("\n");
          } break;

        case NPIC_IMAGE_4L : {
            Npic_image_4l *DT = NpicCastImage (nDT),
                          *MA = NpicCastImage (nMA);
            Npic_mask_4l *Mlut = NpicCastMask (nMlut);
            int t, z, y, x;

            for (ncol = 0; ncol < Mlut->nb; ncol++)
                    if (Mlut->v[ncol].h > Rmax) break;

            for (i = 0; i < ncol; i++)
            {
                printf ("Col %d /%d \r", i+1, ncol); fflush (stdout);

                vg.x = Mlut->v[i].x; vg.y = Mlut->v[i].y;
                vg.z = Mlut->v[i].z; vg.t = Mlut->v[i].t;
                NpicMlutCompLutCol (nCTg, &vg, Rmax, col);

                for (t = 0; t < DT->tmax; t++)
                for (z = 0; z < DT->zmax; z++)
                for (y = 0; y < DT->ymax; y++)
                for (x = 0; x < DT->xmax; x++)
                if (MA->pix[t][z][y][x] != 0)
                {
                    k = NpicMlutIsMAcol_4l (DT, Mlut, col, Rmax, &vg, t, z, y, x);
                    if (k == 0) MA->pix[t][z][y][x] = 0;
                }
            }
            printf ("\n");
          } break;

        case NPIC_IMAGE_5L : {
            Npic_image_5l *DT = NpicCastImage (nDT),
                          *MA = NpicCastImage (nMA);
            Npic_mask_5l *Mlut = NpicCastMask (nMlut);
            int s, t, z, y, x;

            for (ncol = 0; ncol < Mlut->nb; ncol++)
                    if (Mlut->v[ncol].h > Rmax) break;

            for (i = 0; i < ncol; i++)
            {
                printf ("Col %d /%d \r", i+1, ncol); fflush (stdout);

                vg.x = Mlut->v[i].x; vg.y = Mlut->v[i].y;
                vg.z = Mlut->v[i].z; vg.t = Mlut->v[i].t;
                vg.s = Mlut->v[i].s; 
                NpicMlutCompLutCol (nCTg, &vg, Rmax, col);

                for (s = 0; s < DT->smax; s++)
                for (t = 0; t < DT->tmax; t++)
                for (z = 0; z < DT->zmax; z++)
                for (y = 0; y < DT->ymax; y++)
                for (x = 0; x < DT->xmax; x++)
                if (MA->pix[s][t][z][y][x] != 0)
                {
                    k = NpicMlutIsMAcol_5l (DT, Mlut, col, Rmax, &vg, s, t, z, y, x);
                    if (k == 0) MA->pix[s][t][z][y][x] = 0;
                }
            }
            printf ("\n");
          } break;

        case NPIC_IMAGE_6L : {
            Npic_image_6l *DT = NpicCastImage (nDT),
                          *MA = NpicCastImage (nMA);
            Npic_mask_6l *Mlut = NpicCastMask (nMlut);
            int r, s, t, z, y, x;

            for (ncol = 0; ncol < Mlut->nb; ncol++)
                    if (Mlut->v[ncol].h > Rmax) break;

            for (i = 0; i < ncol; i++)
            {
                printf ("Col %d /%d \r", i+1, ncol); fflush (stdout);

                vg.x = Mlut->v[i].x; vg.y = Mlut->v[i].y;
                vg.z = Mlut->v[i].z; vg.t = Mlut->v[i].t;
                vg.s = Mlut->v[i].s; vg.r = Mlut->v[i].r;
                NpicMlutCompLutCol (nCTg, &vg, Rmax, col);

                for (r = 0; r < DT->rmax; r++)
                for (s = 0; s < DT->smax; s++)
                for (t = 0; t < DT->tmax; t++)
                for (z = 0; z < DT->zmax; z++)
                for (y = 0; y < DT->ymax; y++)
                for (x = 0; x < DT->xmax; x++)
                if (MA->pix[r][s][t][z][y][x] != 0)
                {
                    k = NpicMlutIsMAcol_6l (DT, Mlut, col, Rmax, &vg, r, s, t, z, y, x);
                    if (k == 0) MA->pix[r][s][t][z][y][x] = 0;
                }
            }
            printf ("\n");
          } break;

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

    free (col);
}


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


int NpicMlutIsMAlut_2l (Npic_image_2l *DT, Npic_mask_2l *Mlut, Npic_Lut *Lut,
        int ncol, int y, int x)
{
    int i, ya, xa;
    Npic_gsym gs;
    Npic_l val = DT->pix[y][x];

    if (val <= 0 || val >= Lut->col_len) return 0;

    for (i = 0; i < ncol; i++) {
        gs.x[1] = Mlut->v[i].x; gs.x[2] = Mlut->v[i].y; 
        gs.dim = 2; gs.half = 0; /* full mask */
                
        NpicGsymInit (&gs);
        while (NpicGsymNext (&gs)) {
            xa = x + gs.x[1]; ya = y + gs.x[2]; 
            if ( xa >= 0 && xa < DT->xmax &&
                 ya >= 0 && ya < DT->ymax &&
                 DT->pix[ya][xa] >= Lut->v[i][val] )
              return 0;
        }
    }

    return 1;
}


int NpicMlutIsMAlut_3l (Npic_image_3l *DT, Npic_mask_3l *Mlut, Npic_Lut *Lut,
        int ncol, int z, int y, int x)
{
    int i, za, ya, xa;
    Npic_gsym gs;
    Npic_l val = DT->pix[z][y][x];

    if (val <= 0 || val >= Lut->col_len) return 0;

    for (i = 0; i < ncol; i++) {
        gs.x[1] = Mlut->v[i].x; gs.x[2] = Mlut->v[i].y; 
        gs.x[3] = Mlut->v[i].z; 
        gs.dim = 3; gs.half = 0; /* full mask */
                
        NpicGsymInit (&gs);
        while (NpicGsymNext (&gs)) {
            xa = x + gs.x[1]; ya = y + gs.x[2]; 
            za = z + gs.x[3]; 
            if ( xa >= 0 && xa < DT->xmax &&
                 ya >= 0 && ya < DT->ymax &&
                 za >= 0 && za < DT->zmax &&
                 DT->pix[za][ya][xa] >= Lut->v[i][val] )
              return 0;
        }
    }

    return 1;
}


int NpicMlutIsMAlut_4l (Npic_image_4l *DT, Npic_mask_4l *Mlut, Npic_Lut *Lut,
        int ncol, int t, int z, int y, int x)
{
    int i, ta, za, ya, xa;
    Npic_gsym gs;
    Npic_l val = DT->pix[t][z][y][x];

    if (val <= 0 || val >= Lut->col_len) return 0;

    for (i = 0; i < ncol; i++) {
        gs.x[1] = Mlut->v[i].x; gs.x[2] = Mlut->v[i].y; 
        gs.x[3] = Mlut->v[i].z; gs.x[4] = Mlut->v[i].t;
        gs.dim = 4; gs.half = 0; /* full mask */
                
        NpicGsymInit (&gs);
        while (NpicGsymNext (&gs)) {
            xa = x + gs.x[1]; ya = y + gs.x[2]; 
            za = z + gs.x[3]; ta = t + gs.x[4];
            if ( xa >= 0 && xa < DT->xmax &&
                 ya >= 0 && ya < DT->ymax &&
                 za >= 0 && za < DT->zmax &&
                 ta >= 0 && ta < DT->tmax &&
                 DT->pix[ta][za][ya][xa] >= Lut->v[i][val] )
              return 0;
        }
    }

    return 1;
}


int NpicMlutIsMAlut_5l (Npic_image_5l *DT, Npic_mask_5l *Mlut, Npic_Lut *Lut,
        int ncol, int s, int t, int z, int y, int x)
{
    int i, sa, ta, za, ya, xa;
    Npic_gsym gs;
    Npic_l val = DT->pix[s][t][z][y][x];

    if (val <= 0 || val >= Lut->col_len) return 0;

    for (i = 0; i < ncol; i++) {
        gs.x[1] = Mlut->v[i].x; gs.x[2] = Mlut->v[i].y; 
        gs.x[3] = Mlut->v[i].z; gs.x[4] = Mlut->v[i].t;
        gs.x[5] = Mlut->v[i].s; 
        gs.dim = 5; gs.half = 0; /* full mask */
                
        NpicGsymInit (&gs);
        while (NpicGsymNext (&gs)) {
            xa = x + gs.x[1]; ya = y + gs.x[2]; 
            za = z + gs.x[3]; ta = t + gs.x[4];
            sa = s + gs.x[5]; 
            if ( xa >= 0 && xa < DT->xmax &&
                 ya >= 0 && ya < DT->ymax &&
                 za >= 0 && za < DT->zmax &&
                 ta >= 0 && ta < DT->tmax &&
                 sa >= 0 && sa < DT->smax &&
                 DT->pix[sa][ta][za][ya][xa] >= Lut->v[i][val] )
              return 0;
        }
    }

    return 1;
}


int NpicMlutIsMAlut_6l (Npic_image_6l *DT, Npic_mask_6l *Mlut, Npic_Lut *Lut,
        int ncol, int r, int s, int t, int z, int y, int x)
{
    int i, ra, sa, ta, za, ya, xa;
    Npic_gsym gs;
    Npic_l val = DT->pix[r][s][t][z][y][x];

    if (val <= 0 || val >= Lut->col_len) return 0;

    for (i = 0; i < ncol; i++) {
        gs.x[1] = Mlut->v[i].x; gs.x[2] = Mlut->v[i].y; 
        gs.x[3] = Mlut->v[i].z; gs.x[4] = Mlut->v[i].t;
        gs.x[5] = Mlut->v[i].s; gs.x[6] = Mlut->v[i].r;
        gs.dim = 6; gs.half = 0; /* full mask */
                
        NpicGsymInit (&gs);
        while (NpicGsymNext (&gs)) {
            xa = x + gs.x[1]; ya = y + gs.x[2]; 
            za = z + gs.x[3]; ta = t + gs.x[4];
            sa = s + gs.x[5]; ra = r + gs.x[6];
            if ( xa >= 0 && xa < DT->xmax &&
                 ya >= 0 && ya < DT->ymax &&
                 za >= 0 && za < DT->zmax &&
                 ta >= 0 && ta < DT->tmax &&
                 sa >= 0 && sa < DT->smax &&
                 ra >= 0 && ra < DT->rmax &&
                 DT->pix[ra][sa][ta][za][ya][xa] >= Lut->v[i][val] )
              return 0;
        }
    }

    return 1;
}


int NpicMlutIsMAcol_2l (Npic_image_2l *DT, Npic_mask_2l *Mlut, Npic_l *col,
        Npic_l Rmax, Npic_vec *vg, int y, int x)
{
    int ya, xa;
    Npic_gsym gs;
    Npic_l val = DT->pix[y][x];

    if (val <= 0 || val > Rmax) return 0;

    gs.x[1] = vg->x; gs.x[2] = vg->y; 
    gs.dim = 2; gs.half = 0; /* full mask */
                
    NpicGsymInit (&gs);
    while (NpicGsymNext (&gs)) {
        xa = x + gs.x[1]; ya = y + gs.x[2]; 
        if ( xa >= 0 && xa < DT->xmax &&
             ya >= 0 && ya < DT->ymax &&
             DT->pix[ya][xa] >= col[val] )
          return 0;
    }

    return 1;
}


int NpicMlutIsMAcol_3l (Npic_image_3l *DT, Npic_mask_3l *Mlut, Npic_l *col,
        Npic_l Rmax, Npic_vec *vg, int z, int y, int x)
{
    int za, ya, xa;
    Npic_gsym gs;
    Npic_l val = DT->pix[z][y][x];

    if (val <= 0 || val > Rmax) return 0;

    gs.x[1] = vg->x; gs.x[2] = vg->y; 
    gs.x[3] = vg->z; 
    gs.dim = 3; gs.half = 0; /* full mask */
                
    NpicGsymInit (&gs);
    while (NpicGsymNext (&gs)) {
        xa = x + gs.x[1]; ya = y + gs.x[2]; 
        za = z + gs.x[3]; 
        if ( xa >= 0 && xa < DT->xmax &&
             ya >= 0 && ya < DT->ymax &&
             za >= 0 && za < DT->zmax &&
             DT->pix[za][ya][xa] >= col[val] )
          return 0;
    }

    return 1;
}


int NpicMlutIsMAcol_4l (Npic_image_4l *DT, Npic_mask_4l *Mlut, Npic_l *col,
        Npic_l Rmax, Npic_vec *vg, int t, int z, int y, int x)
{
    int ta, za, ya, xa;
    Npic_gsym gs;
    Npic_l val = DT->pix[t][z][y][x];

    if (val <= 0 || val > Rmax) return 0;

    gs.x[1] = vg->x; gs.x[2] = vg->y; 
    gs.x[3] = vg->z; gs.x[4] = vg->t;
    gs.dim = 4; gs.half = 0; /* full mask */
                
    NpicGsymInit (&gs);
    while (NpicGsymNext (&gs)) {
        xa = x + gs.x[1]; ya = y + gs.x[2]; 
        za = z + gs.x[3]; ta = t + gs.x[4];
        if ( xa >= 0 && xa < DT->xmax &&
             ya >= 0 && ya < DT->ymax &&
             za >= 0 && za < DT->zmax &&
             ta >= 0 && ta < DT->tmax &&
             DT->pix[ta][za][ya][xa] >= col[val] )
          return 0;
    }

    return 1;
}


int NpicMlutIsMAcol_5l (Npic_image_5l *DT, Npic_mask_5l *Mlut, Npic_l *col,
        Npic_l Rmax, Npic_vec *vg, int s, int t, int z, int y, int x)
{
    int sa, ta, za, ya, xa;
    Npic_gsym gs;
    Npic_l val = DT->pix[s][t][z][y][x];

    if (val <= 0 || val > Rmax) return 0;

    gs.x[1] = vg->x; gs.x[2] = vg->y; 
    gs.x[3] = vg->z; gs.x[4] = vg->t;
    gs.x[5] = vg->s; 
    gs.dim = 5; gs.half = 0; /* full mask */
                
    NpicGsymInit (&gs);
    while (NpicGsymNext (&gs)) {
        xa = x + gs.x[1]; ya = y + gs.x[2]; 
        za = z + gs.x[3]; ta = t + gs.x[4];
        sa = s + gs.x[5]; 
        if ( xa >= 0 && xa < DT->xmax &&
             ya >= 0 && ya < DT->ymax &&
             za >= 0 && za < DT->zmax &&
             ta >= 0 && ta < DT->tmax &&
             sa >= 0 && sa < DT->smax &&
             DT->pix[sa][ta][za][ya][xa] >= col[val] )
          return 0;
    }

    return 1;
}


int NpicMlutIsMAcol_6l (Npic_image_6l *DT, Npic_mask_6l *Mlut, Npic_l *col,
        Npic_l Rmax, Npic_vec *vg, int r, int s, int t, int z, int y, int x)
{
    int ra, sa, ta, za, ya, xa;
    Npic_gsym gs;
    Npic_l val = DT->pix[r][s][t][z][y][x];

    if (val <= 0 || val > Rmax) return 0;

    gs.x[1] = vg->x; gs.x[2] = vg->y; 
    gs.x[3] = vg->z; gs.x[4] = vg->t;
    gs.x[5] = vg->s; gs.x[6] = vg->r;
    gs.dim = 6; gs.half = 0; /* full mask */
                
    NpicGsymInit (&gs);
    while (NpicGsymNext (&gs)) {
        xa = x + gs.x[1]; ya = y + gs.x[2]; 
        za = z + gs.x[3]; ta = t + gs.x[4];
        sa = s + gs.x[5]; ra = r + gs.x[6];
        if ( xa >= 0 && xa < DT->xmax &&
             ya >= 0 && ya < DT->ymax &&
             za >= 0 && za < DT->zmax &&
             ta >= 0 && ta < DT->tmax &&
             sa >= 0 && sa < DT->smax &&
             ra >= 0 && ra < DT->rmax &&
             DT->pix[ra][sa][ta][za][ya][xa] >= col[val] )
          return 0;
    }

    return 1;
}


/* 
 * Try to find L such that the weighted ball of radius Rtarget will fit in image L^n. 
*/

int NpicMlutComp_L_WD (Npic_mask *nMg, Npic_l Rtarget)
{
    int i, j, L = 1;

    switch (nMg->type) {

        case NPIC_MASK_2L : {
            Npic_mask_2l *Mg = NpicCastMask (nMg);
            for (i = 1, j = 0; i < Mg->nb; i++)
                if (Mg->v[i].h < Mg->v[j].h) j = i;
            if (Mg->nb > 0) L = abs(Mg->v[j].x) * Rtarget / Mg->v[j].h + 1;
        } break;

        case NPIC_MASK_3L : {
            Npic_mask_3l *Mg = NpicCastMask (nMg);
            for (i = 1, j = 0; i < Mg->nb; i++)
                if (Mg->v[i].h < Mg->v[j].h) j = i;
            if (Mg->nb > 0) L = abs(Mg->v[j].x) * Rtarget / Mg->v[j].h + 1;
        } break;

        case NPIC_MASK_4L : {
            Npic_mask_4l *Mg = NpicCastMask (nMg);
            for (i = 1, j = 0; i < Mg->nb; i++)
                if (Mg->v[i].h < Mg->v[j].h) j = i;
            if (Mg->nb > 0) L = abs(Mg->v[j].x) * Rtarget / Mg->v[j].h + 1;
        } break;

        case NPIC_MASK_5L : {
            Npic_mask_5l *Mg = NpicCastMask (nMg);
            for (i = 1, j = 0; i < Mg->nb; i++)
                if (Mg->v[i].h < Mg->v[j].h) j = i;
            if (Mg->nb > 0) L = abs(Mg->v[j].x) * Rtarget / Mg->v[j].h + 1;
        } break;

        case NPIC_MASK_6L : {
            Npic_mask_6l *Mg = NpicCastMask (nMg);
            for (i = 1, j = 0; i < Mg->nb; i++)
                if (Mg->v[i].h < Mg->v[j].h) j = i;
            if (Mg->nb > 0) L = abs(Mg->v[j].x) * Rtarget / Mg->v[j].h + 1;
        } break;
    }
    return L;
}


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