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

/*
 * mlut_remy.c - 08/01/2009
 *
 * Computation of Medial Axis, LUT and Mlut for Distance Transforms
 * using Remy-Thiel algorithms.
 *
 * References:
 *
 * [1] E. Remy and E. Thiel, 2005. Exact Medial Axis with Euclidean Distance. 
 *     Image and Vision Computing, 23(2):167-175.
 *
 * [2] 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.
 *
 * [3] E. Thiel, 1994. Les distances de chanfrein en analyse d'images : 
 *     fondements et applications. These de Doctorat, Universite Joseph Fourier,
 *     Grenoble 1. http://pageperso.lif.univ-mrs.fr/~edouard.thiel/these/
*/


#include <npic.h>


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

/*
 * Lut column Computation Algorithm. See [3, chap.5 p.79], [1], [2].
 *
 * Input:  nCTg the distance cone to origin, vg the vector, 
 *         Rmax the greatest radius to be verified in Lut.
 * Output: The Lut column LutCol is filled with the correct values.
 * Assume that nCTg is a L^n image and that LutCol is [0..Rmax].
 *
 * Do nothing if nCTg is not ok. set not ok on error. Verbose.
 * Return computation time in second.microsecond
*/

double NpicMlutCompLutCol (Npic_image *nCTg, Npic_vec *vg, 
        Npic_l Rmax, Npic_l *LutCol)
{
    double time1, time2;
    Npic_l r1, r2, ra, rb;

    if (NpicImageIsOK (nCTg, __func__) != NPIC_SUCCESS) 
        return 0;
    if (NpicImageIsSquare (nCTg, __func__) != NPIC_SUCCESS) 
      { nCTg->gen.ok = NPIC_ERROR; return 0; }

    time1 = NpicGetTime();

    /* Initialize LutCol to 0 */
    for (ra = 0; ra <= Rmax; ra++)
        LutCol[ra] = 0;

    /* Scan CTg */
    switch (nCTg->type) {

        case NPIC_IMAGE_2L : {   /* - - - - - - - - - - - - - - - - - - - - - */

            Npic_image_2l *CTg = NpicCastImage(nCTg);
            int y, x, L = CTg->xmax;

            for (x = 0; x < L - vg->x; x++)
            for (y = 0; y <= x; y++)
            {
                r1 = CTg->pix[y][x] +1;
                r2 = CTg->pix[y+vg->y][x+vg->x] +1;

                if (r1 <= Rmax && r2 > LutCol[r1]) LutCol[r1] = r2;
            }
          } break;

        case NPIC_IMAGE_3L : {   /* - - - - - - - - - - - - - - - - - - - - - */

            Npic_image_3l *CTg = NpicCastImage(nCTg);
            int z, y, x, L = CTg->xmax;

            for (x = 0; x < L - vg->x; x++)
            for (y = 0; y <= x; y++)
            for (z = 0; z <= y; z++)
            {
                r1 = CTg->pix[z][y][x] +1;
                r2 = CTg->pix[z+vg->z][y+vg->y][x+vg->x] +1;

                if (r1 <= Rmax && r2 > LutCol[r1]) LutCol[r1] = r2;
            }
          } break;

        case NPIC_IMAGE_4L : {   /* - - - - - - - - - - - - - - - - - - - - - */

            Npic_image_4l *CTg = NpicCastImage(nCTg);
            int t, z, y, x, L = CTg->xmax;

            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++)
            {
                r1 = CTg->pix[t][z][y][x] +1;
                r2 = CTg->pix[t+vg->t][z+vg->z][y+vg->y][x+vg->x] +1;

                if (r1 <= Rmax && r2 > LutCol[r1]) LutCol[r1] = r2;
            }
          } break;

        case NPIC_IMAGE_5L : {   /* - - - - - - - - - - - - - - - - - - - - - */

            Npic_image_5l *CTg = NpicCastImage(nCTg);
            int s, t, z, y, x, L = CTg->xmax;

            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++)
            for (s = 0; s <= t; s++)
            {
                r1 = CTg->pix[s][t][z][y][x] +1;
                r2 = CTg->pix[s+vg->s][t+vg->t][z+vg->z][y+vg->y][x+vg->x] +1;

                if (r1 <= Rmax && r2 > LutCol[r1]) LutCol[r1] = r2;
            }
          } break;

        case NPIC_IMAGE_6L : {   /* - - - - - - - - - - - - - - - - - - - - - */

            Npic_image_6l *CTg = NpicCastImage(nCTg);
            int r, s, t, z, y, x, L = CTg->xmax;

            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++)
            for (s = 0; s <= t; s++)
            for (r = 0; r <= s; r++)
            {
                r1 = CTg->pix[r][s][t][z][y][x] +1;
                r2 = CTg->pix[r+vg->r][s+vg->s][t+vg->t][z+vg->z][y+vg->y][x+vg->x] +1;

                if (r1 <= Rmax && r2 > LutCol[r1]) LutCol[r1] = r2;
            }
          } break;

        default : NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
    }

    /* Finalize inclusions */
    rb = 0;
    for (ra = 0; ra <= Rmax; ra++)
    {
        if (LutCol[ra] > rb)
             rb = LutCol[ra];
        else LutCol[ra] = rb;
    }

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


/*
 * Manage Lut
*/

int NpicMlutInitLut (Npic_Lut *Lut, Npic_l Rtarget)
{
    if (Lut == NULL)
        return NpicError (__func__, NPIC_ERR_NULL_PTR, "");
    if (Rtarget <= 0)
        return NpicError (__func__, NPIC_ERROR, ": Rtarget <= 0");

    Lut->col_nb = 0;                /* Number of allocated columns */
    Lut->col_tot = 0;               /* size of vector of column adresses */
    Lut->col_len = Rtarget+1;       /* size of one column */
    Lut->v = NULL;                  /* Lut->v[0..col_nb-1][0..Rtarget] */

    return NPIC_SUCCESS;
}

#define NPIC_COL_STEP        1024

int NpicMlutNewCol (Npic_Lut *Lut)
{
    Npic_l **old_v = NULL;

    if (Lut == NULL)
        return NpicError (__func__, NPIC_ERR_NULL_PTR, "");
    
    if (Lut->col_tot == Lut->col_nb) {
        Lut->col_tot += NPIC_COL_STEP;
        old_v = Lut->v;
        Lut->v = realloc (old_v, Lut->col_tot * sizeof (Npic_l*));
    } 

    if (Lut->v == NULL) {
        Lut->v = old_v;
        return NpicError (__func__, NPIC_ERR_MALLOC, "");
    }

    Lut->v[Lut->col_nb] = malloc (Lut->col_len * sizeof (Npic_l));
    if (Lut->v[Lut->col_nb] == NULL)
        return NpicError (__func__, NPIC_ERR_MALLOC, "");

    Lut->col_nb++;

    return NPIC_SUCCESS;
}


void NpicMlutFreeLut (Npic_Lut *Lut)
{
    int i;

    if (Lut == NULL)
      { NpicError (__func__, NPIC_ERR_NULL_PTR, ""); return ; }

    for (i = 0; i < Lut->col_nb; i++)
        free (Lut->v[i]);
    free (Lut->v);

    NpicMlutInitLut (Lut, 1);
}


/*
 * Full Computation Algorithm of Lut and Mlut for SEDT [1] or WDT [2].
 *
 * Input:  nCTg a L^n image, nMlut the generator of the
 *         Lut mask, Lut the look-up table, Rknown the last verified radius, 
 *         Rtarget the maximum radius to be verified,
 *         nMg == NULL for SEDT, else the generator of the WD mask.
 * Output: nMlut and Lut are filled with the correct values.
 *
 * Do nothing if nCTg is not ok. Verbose. Return error code.
 *
 * Usage: to compute nMlut and Lut from beginning to Rtarget 
 *        with L such that L*L > Rtarget for SED :
 *
 *   Npic_mask  *nMlut = NpicCreateMask (NPIC_MASK_..);
 *   Npic_image *nCTg  = NpicCreateImage_.. (L, L, ..., 0, 0, ...);
 *   Npic_Lut Lut;
 *   char s[200]; 
 *
 *   NpicMaskAddProp (nMlut, "Nature", "MedialAxis");
 *   NpicMaskAddProp (nMlut, "RVerified", "0");
 *   NpicMlutCTg_SED (nCTg); 
 *   NpicMlutInitLut (&Lut, Rtarget);
 *
 *   NpicMlutCompLutMask (nCTg, nMlut, &Lut, 0, Rtarget, NULL);
 *
 *   sprintf (s, "%d", Rtarget);
 *   NpicMaskAddProp (nMlut, "RVerified", s);
 *   NpicWriteMask (nMlut, "name.nmask");
 *
 *   NpicDestroyMask (nMlut);
 *   NpicDestroyImage (nCTg);
 *   NpicMlutFreeLut (&Lut);
*/

int NpicMlutCompLutMask (Npic_image *nCTg, Npic_mask *nMlut, Npic_Lut *Lut, 
        Npic_l Rknown, Npic_l Rtarget, Npic_mask *nMg)
{
    Npic_image *nDTg;
    char *Possible;
    Npic_l R;

    if (NpicImageIsOK (nCTg, __func__) != NPIC_SUCCESS)
        return NPIC_ERROR;
    if (NpicImageIsSquare (nCTg, __func__) != NPIC_SUCCESS ||
        NpicMaskIsOK (nMlut, __func__) != NPIC_SUCCESS)
        return nCTg->gen.ok = NPIC_ERROR;
    if (NpicMaskCompat (nMlut, nCTg) != NPIC_SUCCESS)
        return nCTg->gen.ok = NpicError (__func__, NPIC_ERR_INCOMPAT, ": Mlut and CTg");
    if (nMg != NULL && NpicMaskIsOK (nMg, __func__) != NPIC_SUCCESS) 
        return nCTg->gen.ok = NPIC_ERROR;
    if (nMg != NULL && NpicMaskCompat (nMg, nCTg) != NPIC_SUCCESS)
        return nCTg->gen.ok = NpicError (__func__, NPIC_ERR_INCOMPAT, ": Mg and CTg");

    /* Init DTg to 0 once. During the next loop, we can avoid to re-init DTg 
       because the ball grows with R. */ 
    nDTg = NpicDupImage (nCTg);
    if (nDTg == NULL) return nCTg->gen.ok = NPIC_ERROR;
    NpicFillWholeZero (nDTg);
    
    /* Mark possible values in Possible[] */
    Possible = malloc ((Rtarget+1) * sizeof(char));
    if (Possible == NULL) {
        NpicDestroyImage (nDTg); nCTg->gen.ok = NPIC_ERROR;
        return NpicError (__func__, NPIC_ERR_NULL_PTR, ""); 
    }
    NpicMlutPossibleValues (nCTg, Possible, Rtarget);

    /* Compute Lut columns for current Lut mask */
    if (NpicMlutRecompLutCols (nCTg, nMlut, Lut, Rtarget) != NPIC_SUCCESS) {
        free (Possible); NpicDestroyImage (nDTg);
        return nCTg->gen.ok = NPIC_ERROR;
    }

    /* Main loop */
    for (R = Rknown+1; R <= Rtarget; R++)
    if (Possible[R])
    {
        printf ("R = %"NPIC_PL" / %"NPIC_PL"\r", R, Rtarget); fflush (stdout);

        if (nMg == NULL)
             NpicMlutDTg_Hirata (nCTg, nDTg, R);   /* SEDT */
        else NpicMlutDTg_WD (nCTg, nDTg, nMg, R);  /* WDT  */
        if (nDTg->gen.ok != NPIC_SUCCESS) break;

        if (NpicMlutScanAndInsert (nCTg, nDTg, nMlut, Lut, R, Rtarget) != NPIC_SUCCESS) {
            free (Possible); NpicDestroyImage (nDTg);
            return nCTg->gen.ok = NPIC_ERROR;
        }
    }
    printf ("\n");

    free (Possible);
    NpicDestroyImage (nDTg);
    return NPIC_SUCCESS;
}


/*
 * Create and compute one Lut col. Return the column or NULL.
 * The result must be freed after use.
 *
 * Do nothing if nDT, nCTg or nMlut is not ok. Verbose.
*/

Npic_l *NpicAllocAndCompLutCol (Npic_mask *nMlut, Npic_image *nCTg, 
    Npic_l Rmax, int n)
{
    Npic_l *col;
    Npic_vec vg;

    if (NpicImageIsOK (nCTg, __func__) != NPIC_SUCCESS ||
        NpicMaskIsOK (nMlut, __func__) != NPIC_SUCCESS)
        return NULL;
    if (NpicMaskCompat (nMlut, nCTg) != NPIC_SUCCESS) {
        NpicError (__func__, NPIC_ERR_INCOMPAT, ""); 
        return NULL; 
    }
    if (n < 0 || n >= nMlut->gen.nb) {
        NpicError (__func__, NPIC_ERROR, ": n = %d is not in [0 .. %d[", 
            n, nMlut->gen.nb); 
        return NULL; 
    }

    col = malloc (sizeof(Npic_l) * (Rmax+1));
    if (col == NULL) {
        NpicError (__func__, NPIC_ERR_MALLOC, ""); 
        return NULL; 
    }

    switch (nCTg->type) {

        case NPIC_IMAGE_2L : {
            Npic_mask_2l *Mlut = NpicCastMask (nMlut);
            vg.x = Mlut->v[n].x; vg.y = Mlut->v[n].y;
            NpicMlutCompLutCol (nCTg, &vg, Rmax, col);
          } break;

        case NPIC_IMAGE_3L : {
            Npic_mask_3l *Mlut = NpicCastMask (nMlut);
            vg.x = Mlut->v[n].x; vg.y = Mlut->v[n].y;
            vg.z = Mlut->v[n].z; 
            NpicMlutCompLutCol (nCTg, &vg, Rmax, col);
          } break;

        case NPIC_IMAGE_4L : {
            Npic_mask_4l *Mlut = NpicCastMask (nMlut);
            vg.x = Mlut->v[n].x; vg.y = Mlut->v[n].y;
            vg.z = Mlut->v[n].z; vg.t = Mlut->v[n].t;
            NpicMlutCompLutCol (nCTg, &vg, Rmax, col);
          } break;

        case NPIC_IMAGE_5L : {
            Npic_mask_5l *Mlut = NpicCastMask (nMlut);
            vg.x = Mlut->v[n].x; vg.y = Mlut->v[n].y;
            vg.z = Mlut->v[n].z; vg.t = Mlut->v[n].t;
            vg.s = Mlut->v[n].s; 
            NpicMlutCompLutCol (nCTg, &vg, Rmax, col);
          } break;

        case NPIC_IMAGE_6L : {
            Npic_mask_6l *Mlut = NpicCastMask (nMlut);
            vg.x = Mlut->v[n].x; vg.y = Mlut->v[n].y;
            vg.z = Mlut->v[n].z; vg.t = Mlut->v[n].t;
            vg.s = Mlut->v[n].s; vg.r = Mlut->v[n].r;
            NpicMlutCompLutCol (nCTg, &vg, Rmax, col);
          } break;

        default : NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");

    }

    return col;
}


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

/*
 * Fast extraction of MA points from Bd inter G(Z^n).
 *
 * Input:  nDTg the distance transform of the section of the ball,
 *         nMlut = the generator of the Lut mask,
 *         Lut the look-up table, P(t,z,y,x) the point to test;
 * Output: returns 1 if P is detected as MA in the DTg.
 * Assume nDTg is ok and is a L^n image.
*/

int NpicMlutIsMAg_2l (Npic_image_2l *DTg, Npic_mask_2l *Mlut, Npic_Lut *Lut, 
        int y, int x)
{
    int ya, xa, i;
    Npic_l val;

    val = DTg->pix[y][x];
    if (val <= 0 || val >= Lut->col_len) return 0;  /* cannot test */

    for (i = 0; i < Mlut->nb; i++)
    {
        xa = x - Mlut->v[i].x;
        ya = y - Mlut->v[i].y;

        if ( 0 <= ya && ya <= xa &&
             DTg->pix[ya][xa] >= Lut->v[i][val] )
            return 0;  /* P+v[i] forbids P */
    }

    return 1;  /* P is a MA point ! */
}

int NpicMlutIsMAg_3l (Npic_image_3l *DTg, Npic_mask_3l *Mlut, Npic_Lut *Lut, 
        int z, int y, int x)
{
    int za, ya, xa, i;
    Npic_l val;

    val = DTg->pix[z][y][x];
    if (val <= 0 || val >= Lut->col_len) return 0;  /* cannot test */

    for (i = 0; i < Mlut->nb; i++)
    {
        xa = x - Mlut->v[i].x;
        ya = y - Mlut->v[i].y;
        za = z - Mlut->v[i].z;

        if ( 0 <= za && za <= ya && ya <= xa &&
             DTg->pix[za][ya][xa] >= Lut->v[i][val] )
            return 0;  /* P+v[i] forbids P */
    }

    return 1;  /* P is a MA point ! */
}

int NpicMlutIsMAg_4l (Npic_image_4l *DTg, Npic_mask_4l *Mlut, Npic_Lut *Lut, 
        int t, int z, int y, int x)
{
    int ta, za, ya, xa, i;
    Npic_l val;

    val = DTg->pix[t][z][y][x];
    if (val <= 0 || val >= Lut->col_len) return 0;  /* cannot test */

    for (i = 0; i < Mlut->nb; i++)
    {
        xa = x - Mlut->v[i].x;
        ya = y - Mlut->v[i].y;
        za = z - Mlut->v[i].z;
        ta = t - Mlut->v[i].t;

        if ( 0 <= ta && ta <= za && za <= ya && ya <= xa &&
             DTg->pix[ta][za][ya][xa] >= Lut->v[i][val] )
            return 0;  /* P+v[i] forbids P */
    }

    return 1;  /* P is a MA point ! */
}

int NpicMlutIsMAg_5l (Npic_image_5l *DTg, Npic_mask_5l *Mlut, Npic_Lut *Lut, 
        int s, int t, int z, int y, int x)
{
    int sa, ta, za, ya, xa, i;
    Npic_l val;

    val = DTg->pix[s][t][z][y][x];
    if (val <= 0 || val >= Lut->col_len) return 0;  /* cannot test */

    for (i = 0; i < Mlut->nb; i++)
    {
        xa = x - Mlut->v[i].x;
        ya = y - Mlut->v[i].y;
        za = z - Mlut->v[i].z;
        ta = t - Mlut->v[i].t;
        sa = s - Mlut->v[i].s;

        if ( 0 <= sa && sa <= ta && ta <= za && za <= ya && ya <= xa &&
             DTg->pix[sa][ta][za][ya][xa] >= Lut->v[i][val] )
            return 0;  /* P+v[i] forbids P */
    }

    return 1;  /* P is a MA point ! */
}

int NpicMlutIsMAg_6l (Npic_image_6l *DTg, Npic_mask_6l *Mlut, Npic_Lut *Lut, 
        int r, int s, int t, int z, int y, int x)
{
    int ra, sa, ta, za, ya, xa, i;
    Npic_l val;

    val = DTg->pix[r][s][t][z][y][x];
    if (val <= 0 || val >= Lut->col_len) return 0;  /* cannot test */

    for (i = 0; i < Mlut->nb; i++)
    {
        xa = x - Mlut->v[i].x;
        ya = y - Mlut->v[i].y;
        za = z - Mlut->v[i].z;
        ta = t - Mlut->v[i].t;
        sa = s - Mlut->v[i].s;
        ra = r - Mlut->v[i].r;

        if ( 0 <= ra && ra <= sa && sa <= ta && ta <= za && za <= ya && ya <= xa &&
             DTg->pix[ra][sa][ta][za][ya][xa] >= Lut->v[i][val] )
            return 0;  /* P+v[i] forbids P */
    }

    return 1;  /* P is a MA point ! */
}


/* Compute Lut columns for current Lut mask */

int NpicMlutRecompLutCols (Npic_image *nCTg, Npic_mask *nMlut, Npic_Lut *Lut, int Rtarget)
{
    int i, k;
    Npic_vec vg;

    switch (nMlut->type) {

        case NPIC_MASK_2L : {
            Npic_mask_2l *Mlut = NpicCastMask (nMlut);

            printf ("Allocating %d columns ...\n", Mlut->nb);
            for (i = 0; i < Mlut->nb; i++) {
                k = NpicMlutNewCol (Lut);
                if (k != NPIC_SUCCESS) return k;
            }

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for private(vg) schedule(dynamic)
#           endif
            for (i = 0; i < Mlut->nb; i++) {
                vg.x = Mlut->v[i].x; 
                vg.y = Mlut->v[i].y;

#               ifdef NPIC_USE_OMP
                if (omp_get_thread_num() == 0)
#               endif
                printf ("Computing column %d / %d ...\r", i+1, Mlut->nb); 
                fflush (stdout);

                NpicMlutCompLutCol (nCTg, &vg, Rtarget, Lut->v[i]);
            }
          } break;

        case NPIC_MASK_3L : {
            Npic_mask_3l *Mlut = NpicCastMask (nMlut);

            printf ("Allocating %d columns ...\n", Mlut->nb);
            for (i = 0; i < Mlut->nb; i++) {
                k = NpicMlutNewCol (Lut);
                if (k != NPIC_SUCCESS) return k;
            }

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for private(vg) schedule(dynamic)
#           endif
            for (i = 0; i < Mlut->nb; i++) {
                vg.x = Mlut->v[i].x; 
                vg.y = Mlut->v[i].y;
                vg.z = Mlut->v[i].z;

#               ifdef NPIC_USE_OMP
                if (omp_get_thread_num() == 0)
#               endif
                printf ("Computing column %d / %d ...\r", i+1, Mlut->nb); 
                fflush (stdout);

                NpicMlutCompLutCol (nCTg, &vg, Rtarget, Lut->v[i]);
            }
          } break;

        case NPIC_MASK_4L : {
            Npic_mask_4l *Mlut = NpicCastMask (nMlut);

            printf ("Allocating %d columns ...\n", Mlut->nb);
            for (i = 0; i < Mlut->nb; i++) {
                k = NpicMlutNewCol (Lut);
                if (k != NPIC_SUCCESS) return k;
            }

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for private(vg) schedule(dynamic)
#           endif
            for (i = 0; i < Mlut->nb; i++) {
                vg.x = Mlut->v[i].x; 
                vg.y = Mlut->v[i].y;
                vg.z = Mlut->v[i].z;
                vg.t = Mlut->v[i].t;

#               ifdef NPIC_USE_OMP
                if (omp_get_thread_num() == 0)
#               endif
                printf ("Computing column %d / %d ...\r", i+1, Mlut->nb); 
                fflush (stdout);

                NpicMlutCompLutCol (nCTg, &vg, Rtarget, Lut->v[i]);
            }
          } break;

        case NPIC_MASK_5L : {
            Npic_mask_5l *Mlut = NpicCastMask (nMlut);

            printf ("Allocating %d columns ...\n", Mlut->nb);
            for (i = 0; i < Mlut->nb; i++) {
                k = NpicMlutNewCol (Lut);
                if (k != NPIC_SUCCESS) return k;
            }

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for private(vg) schedule(dynamic)
#           endif
            for (i = 0; i < Mlut->nb; i++) {
                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;

#               ifdef NPIC_USE_OMP
                if (omp_get_thread_num() == 0)
#               endif
                printf ("Computing column %d / %d ...\r", i+1, Mlut->nb); 
                fflush (stdout);

                NpicMlutCompLutCol (nCTg, &vg, Rtarget, Lut->v[i]);
            }
          } break;

        case NPIC_MASK_6L : {
            Npic_mask_6l *Mlut = NpicCastMask (nMlut);

            printf ("Allocating %d columns ...\n", Mlut->nb);
            for (i = 0; i < Mlut->nb; i++) {
                k = NpicMlutNewCol (Lut);
                if (k != NPIC_SUCCESS) return k;
            }

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for private(vg) schedule(dynamic)
#           endif
            for (i = 0; i < Mlut->nb; i++) {
                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;

#               ifdef NPIC_USE_OMP
                if (omp_get_thread_num() == 0)
#               endif
                printf ("Computing column %d / %d ...\r", i+1, Mlut->nb); 
                fflush (stdout);

                NpicMlutCompLutCol (nCTg, &vg, Rtarget, Lut->v[i]);
            }
          } break;

        default : return NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
    }
    printf ("\n");

    return NPIC_SUCCESS;
}


/* Compute possible values */

void NpicMlutPossibleValues (Npic_image *nCTg, char *Possible, int Rtarget)
{
    int i, L = nCTg->gen.xmax;
    Npic_l val;

    for (i = 1; i <= Rtarget; i++)
        Possible[i] = 0;

    switch (nCTg->type) {

        case NPIC_IMAGE_2L : {
            Npic_image_2l *CTg = NpicCastImage (nCTg);
            int y, x;

            for (x = 1; x <  L; x++)
            for (y = 0; y <= x; y++)
            {
                val = CTg->pix[y][x];
                if (0 <= val && val <= Rtarget) Possible[val] = 1;
            }
          } break;

        case NPIC_IMAGE_3L : {
            Npic_image_3l *CTg = NpicCastImage (nCTg);
            int z, y, x;

            for (x = 1; x <  L; x++)
            for (y = 0; y <= x; y++)
            for (z = 0; z <= y; z++)
            {
                val = CTg->pix[z][y][x];
                if (0 <= val && val <= Rtarget) Possible[val] = 1;
            }
          } break;

        case NPIC_IMAGE_4L : {
            Npic_image_4l *CTg = NpicCastImage (nCTg);
            int t, z, y, x;

            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->pix[t][z][y][x];
                if (0 <= val && val <= Rtarget) Possible[val] = 1;
            }
          } break;

        case NPIC_IMAGE_5L : {
            Npic_image_5l *CTg = NpicCastImage (nCTg);
            int s, t, z, y, x;

            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++)
            {
                val = CTg->pix[s][t][z][y][x];
                if (0 <= val && val <= Rtarget) Possible[val] = 1;
            }
          } break;

        case NPIC_IMAGE_6L : {
            Npic_image_6l *CTg = NpicCastImage (nCTg);
            int r, s, t, z, y, x;

            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++)
            {
                val = CTg->pix[r][s][t][z][y][x];
                if (0 <= val && val <= Rtarget) Possible[val] = 1;
            }
          } break;

        default : NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
    }
}


/* Scan DTg for MA points */

int NpicMlutScanAndInsert (Npic_image *nCTg, Npic_image *nDTg, Npic_mask *nMlut,
        Npic_Lut *Lut, Npic_l R, Npic_l Rtarget)
{
    int i, k;
    Npic_vec vg;

    switch (nDTg->type) {

            case NPIC_IMAGE_2L : {     /* - - - - - - - - - - - - - - - - - - */

                Npic_image_2l *DTg = NpicCastImage (nDTg);
                Npic_mask_2l  *Mlut = NpicCastMask (nMlut);
                int y, x, L = DTg->xmax;

                for (x = 1; x <  L && DTg->pix[0][x] != 0; x++)
                for (y = 0; y <= x && DTg->pix[y][x] != 0; y++)
                if (NpicMlutIsMAg_2l (DTg, Mlut, Lut, y, x))
                {
                    /* Add a new weighting to Mlut */
                    NpicMaskInsert_2l (nMlut, y, x, R);
                    i = nMlut->gen.nb-1;
 
                    printf ("i= %4d  (y,x)= ( %2d, %2d)"
                            "  added for R= %7"NPIC_PL"\n",
                        i+1, y, x, R);
 
                    /* Insert new column in Lut */
                    k = NpicMlutNewCol (Lut);
                    if (k != NPIC_SUCCESS) return k;

                    /* Compute Lut column on CTg */
                    vg.x = Mlut->v[i].x; 
                    vg.y = Mlut->v[i].y;
                    NpicMlutCompLutCol (nCTg, &vg, Rtarget, Lut->v[i]);

                    /* Consistency check */
                    if (NpicMlutIsMAg_2l (DTg, Mlut, Lut, y, x))
                    {
                        NpicError (__func__, NPIC_ERROR, 
                            ": CONSISTENCY ERROR for R = %"NPIC_PL"\n", R);
                        return NPIC_ERROR;
                    }
                }
              } break;

            case NPIC_IMAGE_3L : {     /* - - - - - - - - - - - - - - - - - - */

                Npic_image_3l *DTg = NpicCastImage (nDTg);
                Npic_mask_3l  *Mlut = NpicCastMask (nMlut);
                int z, y, x, L = DTg->xmax;

                for (x = 1; x <  L && DTg->pix[0][0][x] != 0; x++)
                for (y = 0; y <= x && DTg->pix[0][y][x] != 0; y++)
                for (z = 0; z <= y && DTg->pix[z][y][x] != 0; z++)
                if (NpicMlutIsMAg_3l (DTg, Mlut, Lut, z, y, x))
                {
                    /* Add a new weighting to Mlut */
                    NpicMaskInsert_3l (nMlut, z, y, x, R);
                    i = nMlut->gen.nb-1;
 
                    printf ("i= %4d  (z,y,x)= ( %2d, %2d, %2d)"
                            "  added for R= %7"NPIC_PL"\n",
                        i+1, z, y, x, R);
 
                    /* Insert new column in Lut */
                    k = NpicMlutNewCol (Lut);
                    if (k != NPIC_SUCCESS) return k;

                    /* Compute Lut column on CTg */
                    vg.x = Mlut->v[i].x; 
                    vg.y = Mlut->v[i].y;
                    vg.z = Mlut->v[i].z;
                    NpicMlutCompLutCol (nCTg, &vg, Rtarget, Lut->v[i]);

                    /* Consistency check */
                    if (NpicMlutIsMAg_3l (DTg, Mlut, Lut, z, y, x))
                    {
                        NpicError (__func__, NPIC_ERROR, 
                            ": CONSISTENCY ERROR for R = %"NPIC_PL"\n", R);
                        return NPIC_ERROR;
                    }
                }
              } break;

            case NPIC_IMAGE_4L : {     /* - - - - - - - - - - - - - - - - - - */

                Npic_image_4l *DTg = NpicCastImage (nDTg);
                Npic_mask_4l  *Mlut = NpicCastMask (nMlut);
                int t, z, y, x, L = DTg->xmax;

                for (x = 1; x <  L && DTg->pix[0][0][0][x] != 0; x++)
                for (y = 0; y <= x && DTg->pix[0][0][y][x] != 0; y++)
                for (z = 0; z <= y && DTg->pix[0][z][y][x] != 0; z++)
                for (t = 0; t <= z && DTg->pix[t][z][y][x] != 0; t++)
                if (NpicMlutIsMAg_4l (DTg, Mlut, Lut, t, z, y, x))
                {
                    /* Add a new weighting to Mlut */
                    NpicMaskInsert_4l (nMlut, t, z, y, x, R);
                    i = nMlut->gen.nb-1;
 
                    printf ("i= %4d  (t,z,y,x)= ( %2d, %2d, %2d, %2d)"
                            "  added for R= %7"NPIC_PL"\n",
                        i+1, t, z, y, x, R);
 
                    /* Insert new column in Lut */
                    k = NpicMlutNewCol (Lut);
                    if (k != NPIC_SUCCESS) return k;

                    /* Compute Lut column on CTg */
                    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, Rtarget, Lut->v[i]);

                    /* Consistency check */
                    if (NpicMlutIsMAg_4l (DTg, Mlut, Lut, t, z, y, x))
                    {
                        NpicError (__func__, NPIC_ERROR, 
                            ": CONSISTENCY ERROR for R = %"NPIC_PL"\n", R);
                        return NPIC_ERROR;
                    }
                }
              } break;

            case NPIC_IMAGE_5L : {     /* - - - - - - - - - - - - - - - - - - */

                Npic_image_5l *DTg = NpicCastImage (nDTg);
                Npic_mask_5l  *Mlut = NpicCastMask (nMlut);
                int s, t, z, y, x, L = DTg->xmax;

                for (x = 1; x <  L && DTg->pix[0][0][0][0][x] != 0; x++)
                for (y = 0; y <= x && DTg->pix[0][0][0][y][x] != 0; y++)
                for (z = 0; z <= y && DTg->pix[0][0][z][y][x] != 0; z++)
                for (t = 0; t <= z && DTg->pix[0][t][z][y][x] != 0; t++)
                for (s = 0; s <= t && DTg->pix[s][t][z][y][x] != 0; s++)
                if (NpicMlutIsMAg_5l (DTg, Mlut, Lut, s, t, z, y, x))
                {
                    /* Add a new weighting to Mlut */
                    NpicMaskInsert_5l (nMlut, s, t, z, y, x, R);
                    i = nMlut->gen.nb-1;
 
                    printf ("i= %4d  (s,t,z,y,x)= ( %2d, %2d, %2d, %2d, %2d)"
                            "  added for R= %7"NPIC_PL"\n",
                        i+1, s, t, z, y, x, R);
 
                    /* Insert new column in Lut */
                    k = NpicMlutNewCol (Lut);
                    if (k != NPIC_SUCCESS) return k;

                    /* Compute Lut column on CTg */
                    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, Rtarget, Lut->v[i]);

                    /* Consistency check */
                    if (NpicMlutIsMAg_5l (DTg, Mlut, Lut, s, t, z, y, x))
                    {
                        NpicError (__func__, NPIC_ERROR, 
                            ": CONSISTENCY ERROR for R = %"NPIC_PL"\n", R);
                        return NPIC_ERROR;
                    }
                }
              } break;

            case NPIC_IMAGE_6L : {     /* - - - - - - - - - - - - - - - - - - */

                Npic_image_6l *DTg = NpicCastImage (nDTg);
                Npic_mask_6l  *Mlut = NpicCastMask (nMlut);
                int r, s, t, z, y, x, L = DTg->xmax;

                for (x = 1; x <  L && DTg->pix[0][0][0][0][0][x] != 0; x++)
                for (y = 0; y <= x && DTg->pix[0][0][0][0][y][x] != 0; y++)
                for (z = 0; z <= y && DTg->pix[0][0][0][z][y][x] != 0; z++)
                for (t = 0; t <= z && DTg->pix[0][0][t][z][y][x] != 0; t++)
                for (s = 0; s <= t && DTg->pix[0][s][t][z][y][x] != 0; s++)
                for (r = 0; r <= s && DTg->pix[r][s][t][z][y][x] != 0; r++)
                if (NpicMlutIsMAg_6l (DTg, Mlut, Lut, r, s, t, z, y, x))
                {
                    /* Add a new weighting to Mlut */
                    NpicMaskInsert_6l (nMlut, r, s, t, z, y, x, R);
                    i = nMlut->gen.nb-1;
 
                    printf ("i= %4d  (r,s,t,z,y,x)= ( %2d, %2d, %2d, %2d, %2d, %2d)"
                            "  added for R= %7"NPIC_PL"\n",
                        i+1, r, s, t, z, y, x, R);
 
                    /* Insert new column in Lut */
                    k = NpicMlutNewCol (Lut);
                    if (k != NPIC_SUCCESS) return k;

                    /* Compute Lut column on CTg */
                    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, Rtarget, Lut->v[i]);

                    /* Consistency check */
                    if (NpicMlutIsMAg_6l (DTg, Mlut, Lut, r, s, t, z, y, x))
                    {
                        NpicError (__func__, NPIC_ERROR, 
                            ": CONSISTENCY ERROR for R = %"NPIC_PL"\n", R);
                        return NPIC_ERROR;
                    }
                }
              } break;

        default : return NpicError (__func__, NPIC_ERR_UNEX_NPIC, "");
    }

    return NPIC_SUCCESS;
}


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