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

/*
 * dist_sedt_saito.c - 11/12/2008
 *
 * Squared Euclidean Distance Transforms.
 *
 * See comments in dist_sedt_hirata.c
*/


#include <npic.h>


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

/*
 * Compute the Squared Euclidean Distance Transform (SEDT) on image np,
 * with the separable algorithm from Saito et al [1].
 *
 * Do nothing if np is not ok. Set ok on error. Verbose.
 * Return computation time in second.microsecond
*/

double NpicSEDT_Saito (Npic_image *np)
{
    double time1, time2;
    int L;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return 0;

    /* Change the border width of np to 1 if it was smaller */
    NpicSetBorderWidthMin (np, 1);
    
    /* Set border pixels to 0 */
    switch (np->type)
    {
        case NPIC_IMAGE_2L : 
        case NPIC_IMAGE_3L : 
        case NPIC_IMAGE_4L : 
        case NPIC_IMAGE_5L : 
        case NPIC_IMAGE_6L : 
            NpicFillBorder_i (np, 0);
            break;
    }

    L = np->gen.ymax;
    if (np->gen.dim >= 3 && np->gen.zmax > L) L = np->gen.zmax;
    if (np->gen.dim >= 4 && np->gen.tmax > L) L = np->gen.tmax;
    if (np->gen.dim >= 5 && np->gen.smax > L) L = np->gen.smax;
    if (np->gen.dim >= 6 && np->gen.rmax > L) L = np->gen.rmax;
    L += 2; /* for border pixels */

    time1 = NpicGetTime();

    switch (np->type) {

        case NPIC_IMAGE_2L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_2l *p = NpicCastImage(np);
            Npic_l k, j, Cw[L], *Cv=Cw+1;
            int y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(y, x, k, j, Cw, Cv)
            { Cv=Cw+1;  /* again because Cv is private */
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (y = 0; y < p->ymax; y++)
            {
                k = 0;
                for (x = 0; x <= p->xmax; x++)  /* up to border, 0 filled */
                {
                    if (p->pix[y][x] != 0)
                      { k++; p->pix[y][x] = k*k; }  /* forward */
                    else if (k != 0) {
                        for (j = k/2; j >= 1; j--) 
                            p->pix[y][x-j] = j*j;   /* backward */
                        k = 0;
                    }
                }
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (x = 0; x < p->xmax; x++)
            {
                for (y = -1; y <= p->ymax; y++)
                    Cv[y] = p->pix[y][x];

                for (y = 0; y < p->ymax; y++)
                    p->pix[y][x] = Npic_Saito_min (Cv, y, p->ymax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

        case NPIC_IMAGE_3L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_3l *p = NpicCastImage(np);
            Npic_l k, j, Cw[L], *Cv=Cw+1;
            int z, y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(z, y, x, k, j, Cw, Cv)
            { Cv=Cw+1;  /* again because Cv is private */
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            {
                k = 0;
                for (x = 0; x <= p->xmax; x++)  /* up to border, 0 filled */
                {
                    if (p->pix[z][y][x] != 0)
                      { k++; p->pix[z][y][x] = k*k; }  /* forward */
                    else if (k != 0) {
                        for (j = k/2; j >= 1; j--) 
                            p->pix[z][y][x-j] = j*j;   /* backward */
                        k = 0;
                    }
                }
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (z = 0; z < p->zmax; z++)
            for (x = 0; x < p->xmax; x++)
            {
                for (y = -1; y <= p->ymax; y++)
                    Cv[y] = p->pix[z][y][x];

                for (y = 0; y < p->ymax; y++)
                    p->pix[z][y][x] = Npic_Saito_min (Cv, y, p->ymax);
            }

            /*
             * Scan in z
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                for (z = -1; z <= p->zmax; z++)
                    Cv[z] = p->pix[z][y][x];

                for (z = 0; z < p->zmax; z++)
                    p->pix[z][y][x] = Npic_Saito_min (Cv, z, p->zmax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

        case NPIC_IMAGE_4L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_4l *p = NpicCastImage(np);
            Npic_l k, j, Cw[L], *Cv=Cw+1;
            int t, z, y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(t, z, y, x, k, j, Cw, Cv)
            { Cv=Cw+1;  /* again because Cv is private */
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            {
                k = 0;
                for (x = 0; x <= p->xmax; x++)  /* up to border, 0 filled */
                {
                    if (p->pix[t][z][y][x] != 0)
                      { k++; p->pix[t][z][y][x] = k*k; }  /* forward */
                    else if (k != 0) {
                        for (j = k/2; j >= 1; j--) 
                            p->pix[t][z][y][x-j] = j*j;   /* backward */
                        k = 0;
                    }
                }
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (x = 0; x < p->xmax; x++)
            {
                for (y = -1; y <= p->ymax; y++)
                    Cv[y] = p->pix[t][z][y][x];

                for (y = 0; y < p->ymax; y++)
                    p->pix[t][z][y][x] = Npic_Saito_min (Cv, y, p->ymax);
            }

            /*
             * Scan in z
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (t = 0; t < p->tmax; t++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                for (z = -1; z <= p->zmax; z++)
                    Cv[z] = p->pix[t][z][y][x];

                for (z = 0; z < p->zmax; z++)
                    p->pix[t][z][y][x] = Npic_Saito_min (Cv, z, p->zmax);
            }

            /*
             * Scan in t
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                for (t = -1; t <= p->tmax; t++)
                    Cv[t] = p->pix[t][z][y][x];

                for (t = 0; t < p->tmax; t++)
                    p->pix[t][z][y][x] = Npic_Saito_min (Cv, t, p->tmax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

        case NPIC_IMAGE_5L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_5l *p = NpicCastImage(np);
            Npic_l k, j, Cw[L], *Cv=Cw+1;
            int s, t, z, y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(s, t, z, y, x, k, j, Cw, Cv)
            { Cv=Cw+1;  /* again because Cv is private */
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            {
                k = 0;
                for (x = 0; x <= p->xmax; x++)  /* up to border, 0 filled */
                {
                    if (p->pix[s][t][z][y][x] != 0)
                      { k++; p->pix[s][t][z][y][x] = k*k; }  /* forward */
                    else if (k != 0) {
                        for (j = k/2; j >= 1; j--) 
                            p->pix[s][t][z][y][x-j] = j*j;   /* backward */
                        k = 0;
                    }
                }
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (x = 0; x < p->xmax; x++)
            {
                for (y = -1; y <= p->ymax; y++)
                    Cv[y] = p->pix[s][t][z][y][x];

                for (y = 0; y < p->ymax; y++)
                    p->pix[s][t][z][y][x] = Npic_Saito_min (Cv, y, p->ymax);
            }

            /*
             * Scan in z
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                for (z = -1; z <= p->zmax; z++)
                    Cv[z] = p->pix[s][t][z][y][x];

                for (z = 0; z < p->zmax; z++)
                    p->pix[s][t][z][y][x] = Npic_Saito_min (Cv, z, p->zmax);
            }

            /*
             * Scan in t
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                for (t = -1; t <= p->tmax; t++)
                    Cv[t] = p->pix[s][t][z][y][x];

                for (t = 0; t < p->tmax; t++)
                    p->pix[s][t][z][y][x] = Npic_Saito_min (Cv, t, p->tmax);
            }

            /*
             * Scan in s
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                for (s = -1; s <= p->smax; s++)
                    Cv[s] = p->pix[s][t][z][y][x];

                for (s = 0; s < p->smax; s++)
                    p->pix[s][t][z][y][x] = Npic_Saito_min (Cv, s, p->smax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

        case NPIC_IMAGE_6L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_6l *p = NpicCastImage(np);
            Npic_l k, j, Cw[L], *Cv=Cw+1;
            int r, s, t, z, y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(r, s, t, z, y, x, k, j, Cw, Cv)
            { Cv=Cw+1;  /* again because Cv is private */
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            {
                k = 0;
                for (x = 0; x <= p->xmax; x++)  /* up to border, 0 filled */
                {
                    if (p->pix[r][s][t][z][y][x] != 0)
                      { k++; p->pix[r][s][t][z][y][x] = k*k; }  /* forward */
                    else if (k != 0) {
                        for (j = k/2; j >= 1; j--) 
                            p->pix[r][s][t][z][y][x-j] = j*j;   /* backward */
                        k = 0;
                    }
                }
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (x = 0; x < p->xmax; x++)
            {
                for (y = -1; y <= p->ymax; y++)
                    Cv[y] = p->pix[r][s][t][z][y][x];

                for (y = 0; y < p->ymax; y++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_min (Cv, y, p->ymax);
            }

            /*
             * Scan in z
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                for (z = -1; z <= p->zmax; z++)
                    Cv[z] = p->pix[r][s][t][z][y][x];

                for (z = 0; z < p->zmax; z++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_min (Cv, z, p->zmax);
            }

            /*
             * Scan in t
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (s = 0; s < p->smax; s++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                for (t = -1; t <= p->tmax; t++)
                    Cv[t] = p->pix[r][s][t][z][y][x];

                for (t = 0; t < p->tmax; t++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_min (Cv, t, p->tmax);
            }

            /*
             * Scan in s
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                for (s = -1; s <= p->smax; s++)
                    Cv[s] = p->pix[r][s][t][z][y][x];

                for (s = 0; s < p->smax; s++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_min (Cv, s, p->smax);
            }

            /*
             * Scan in r
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                for (r = -1; r <= p->rmax; r++)
                    Cv[r] = p->pix[r][s][t][z][y][x];

                for (r = 0; r < p->rmax; r++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_min (Cv, r, p->rmax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

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

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


/*****************************************************************************/

/*
 * Compute the Squared Euclidean Distance Transform (SEDT_inf) on image np
 * with "infinite border", i.e border pixels propagating no distance value.
 * Derived from the separable algorithm from Saito et al [1].
 *
 * Do nothing if np is not ok. Set ok on error. Verbose.
 * Return computation time in second.microsecond
*/

double NpicSEDT_Saito_inf (Npic_image *np)
{
    double time1, time2;
    int L;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return 0;

    L = np->gen.ymax;
    if (np->gen.dim >= 3 && np->gen.zmax > L) L = np->gen.zmax;
    if (np->gen.dim >= 4 && np->gen.tmax > L) L = np->gen.tmax;
    if (np->gen.dim >= 5 && np->gen.smax > L) L = np->gen.smax;
    if (np->gen.dim >= 6 && np->gen.rmax > L) L = np->gen.rmax;

    time1 = NpicGetTime();

    switch (np->type) {

        case NPIC_IMAGE_2L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_2l *p = NpicCastImage(np);
            Npic_l k, k2, j, Cv[L];
            int y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(y, x, k, k2, j, Cv)
            {
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (y = 0; y < p->ymax; y++)
            {
                /* Forward */
                k = -1;
                for (x = 0; x < p->xmax; x++)
                {
                    if (p->pix[y][x] == 0) k = 0;
                    else if (k == -1) p->pix[y][x] = -1;    /* no distance */
                    else { k++; p->pix[y][x] = k*k; }
                }
                
                /* Backward */
                k = -1;
                for (x = p->xmax-1; x >= 0; x--)
                {
                    if (p->pix[y][x] == 0) k = 0;
                    else if (k != -1) {
                        k++; k2 = k*k; j = p->pix[y][x];
                        if (k2 < j || j == -1) 
                            p->pix[y][x] = k2; 
                    }
                }
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (y = 0; y < p->ymax; y++)
                    Cv[y] = p->pix[y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (y = 0; y < p->ymax; y++)
                    p->pix[y][x] = Npic_Saito_min_inf (Cv, y, p->ymax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

        case NPIC_IMAGE_3L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_3l *p = NpicCastImage(np);
            Npic_l k, k2, j, Cv[L];
            int z, y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(z, y, x, k, k2, j, Cv)
            {
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            {
                /* Forward */
                k = -1;
                for (x = 0; x < p->xmax; x++)
                {
                    if (p->pix[z][y][x] == 0) k = 0;
                    else if (k == -1) p->pix[z][y][x] = -1;    /* no distance */
                    else { k++; p->pix[z][y][x] = k*k; }
                }
                
                /* Backward */
                k = -1;
                for (x = p->xmax-1; x >= 0; x--)
                {
                    if (p->pix[z][y][x] == 0) k = 0;
                    else if (k != -1) {
                        k++; k2 = k*k; j = p->pix[z][y][x];
                        if (k2 < j || j == -1) 
                            p->pix[z][y][x] = k2; 
                    }
                }
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (z = 0; z < p->zmax; z++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (y = 0; y < p->ymax; y++)
                    Cv[y] = p->pix[z][y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (y = 0; y < p->ymax; y++)
                    p->pix[z][y][x] = Npic_Saito_min_inf (Cv, y, p->ymax);
            }

            /*
             * Scan in z
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (z = 0; z < p->zmax; z++)
                    Cv[z] = p->pix[z][y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (z = 0; z < p->zmax; z++)
                    p->pix[z][y][x] = Npic_Saito_min_inf (Cv, z, p->zmax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

        case NPIC_IMAGE_4L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_4l *p = NpicCastImage(np);
            Npic_l k, k2, j, Cv[L];
            int t, z, y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(t, z, y, x, k, k2, j, Cv)
            {
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            {
                /* Forward */
                k = -1;
                for (x = 0; x < p->xmax; x++)
                {
                    if (p->pix[t][z][y][x] == 0) k = 0;
                    else if (k == -1) p->pix[t][z][y][x] = -1;    /* no distance */
                    else { k++; p->pix[t][z][y][x] = k*k; }
                }
                
                /* Backward */
                k = -1;
                for (x = p->xmax-1; x >= 0; x--)
                {
                    if (p->pix[t][z][y][x] == 0) k = 0;
                    else if (k != -1) {
                        k++; k2 = k*k; j = p->pix[t][z][y][x];
                        if (k2 < j || j == -1) 
                            p->pix[t][z][y][x] = k2; 
                    }
                }
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (y = 0; y < p->ymax; y++)
                    Cv[y] = p->pix[t][z][y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (y = 0; y < p->ymax; y++)
                    p->pix[t][z][y][x] = Npic_Saito_min_inf (Cv, y, p->ymax);
            }

            /*
             * Scan in z
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (t = 0; t < p->tmax; t++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (z = 0; z < p->zmax; z++)
                    Cv[z] = p->pix[t][z][y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (z = 0; z < p->zmax; z++)
                    p->pix[t][z][y][x] = Npic_Saito_min_inf (Cv, z, p->zmax);
            }

            /*
             * Scan in t
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (t = 0; t < p->tmax; t++)
                    Cv[t] = p->pix[t][z][y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (t = 0; t < p->tmax; t++)
                    p->pix[t][z][y][x] = Npic_Saito_min_inf (Cv, t, p->tmax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

        case NPIC_IMAGE_5L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_5l *p = NpicCastImage(np);
            Npic_l k, k2, j, Cv[L];
            int s, t, z, y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(s, t, z, y, x, k, k2, j, Cv)
            {
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            {
                /* Forward */
                k = -1;
                for (x = 0; x < p->xmax; x++)
                {
                    if (p->pix[s][t][z][y][x] == 0) k = 0;
                    else if (k == -1) p->pix[s][t][z][y][x] = -1;    /* no distance */
                    else { k++; p->pix[s][t][z][y][x] = k*k; }
                }
                
                /* Backward */
                k = -1;
                for (x = p->xmax-1; x >= 0; x--)
                {
                    if (p->pix[s][t][z][y][x] == 0) k = 0;
                    else if (k != -1) {
                        k++; k2 = k*k; j = p->pix[s][t][z][y][x];
                        if (k2 < j || j == -1) 
                            p->pix[s][t][z][y][x] = k2; 
                    }
                }
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (y = 0; y < p->ymax; y++)
                    Cv[y] = p->pix[s][t][z][y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (y = 0; y < p->ymax; y++)
                    p->pix[s][t][z][y][x] = Npic_Saito_min_inf (Cv, y, p->ymax);
            }

            /*
             * Scan in z
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (z = 0; z < p->zmax; z++)
                    Cv[z] = p->pix[s][t][z][y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (z = 0; z < p->zmax; z++)
                    p->pix[s][t][z][y][x] = Npic_Saito_min_inf (Cv, z, p->zmax);
            }

            /*
             * Scan in t
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (t = 0; t < p->tmax; t++)
                    Cv[t] = p->pix[s][t][z][y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (t = 0; t < p->tmax; t++)
                    p->pix[s][t][z][y][x] = Npic_Saito_min_inf (Cv, t, p->tmax);
            }

            /*
             * Scan in s
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (s = 0; s < p->smax; s++)
                    Cv[s] = p->pix[s][t][z][y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (s = 0; s < p->smax; s++)
                    p->pix[s][t][z][y][x] = Npic_Saito_min_inf (Cv, s, p->smax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

        case NPIC_IMAGE_6L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_6l *p = NpicCastImage(np);
            Npic_l k, k2, j, Cv[L];
            int r, s, t, z, y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(r, s, t, z, y, x, k, k2, j, Cv)
            {
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            {
                /* Forward */
                k = -1;
                for (x = 0; x < p->xmax; x++)
                {
                    if (p->pix[r][s][t][z][y][x] == 0) k = 0;
                    else if (k == -1) p->pix[r][s][t][z][y][x] = -1;    /* no distance */
                    else { k++; p->pix[r][s][t][z][y][x] = k*k; }
                }
                
                /* Backward */
                k = -1;
                for (x = p->xmax-1; x >= 0; x--)
                {
                    if (p->pix[r][s][t][z][y][x] == 0) k = 0;
                    else if (k != -1) {
                        k++; k2 = k*k; j = p->pix[r][s][t][z][y][x];
                        if (k2 < j || j == -1) 
                            p->pix[r][s][t][z][y][x] = k2; 
                    }
                }
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (y = 0; y < p->ymax; y++)
                    Cv[y] = p->pix[r][s][t][z][y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (y = 0; y < p->ymax; y++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_min_inf (Cv, y, p->ymax);
            }

            /*
             * Scan in z
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (z = 0; z < p->zmax; z++)
                    Cv[z] = p->pix[r][s][t][z][y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (z = 0; z < p->zmax; z++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_min_inf (Cv, z, p->zmax);
            }

            /*
             * Scan in t
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (s = 0; s < p->smax; s++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (t = 0; t < p->tmax; t++)
                    Cv[t] = p->pix[r][s][t][z][y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (t = 0; t < p->tmax; t++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_min_inf (Cv, t, p->tmax);
            }

            /*
             * Scan in s
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (s = 0; s < p->smax; s++)
                    Cv[s] = p->pix[r][s][t][z][y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (s = 0; s < p->smax; s++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_min_inf (Cv, s, p->smax);
            }

            /*
             * Scan in r
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (r = 0; r < p->rmax; r++)
                    Cv[r] = p->pix[r][s][t][z][y][x];

                /* Compute min on column Cv and store result in p->pix */
                for (r = 0; r < p->rmax; r++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_min_inf (Cv, r, p->rmax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

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

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


/*****************************************************************************/

/*
 * Compute the Squared Euclidean Reverse Distance Transform (SEDT_rev) on image np.
 * Derived from the separable algorithm from Saito et al [1].
 *
 * Do nothing if np is not ok. Set ok on error. Verbose.
 * Return computation time in second.microsecond
*/

double NpicSEDT_Saito_rev (Npic_image *np)
{
    double time1, time2;
    int L;

    if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return 0;

    L = np->gen.xmax;
    if (np->gen.ymax > L) L = np->gen.ymax;
    if (np->gen.dim >= 3 && np->gen.zmax > L) L = np->gen.zmax;
    if (np->gen.dim >= 4 && np->gen.tmax > L) L = np->gen.tmax;
    if (np->gen.dim >= 5 && np->gen.smax > L) L = np->gen.smax;
    if (np->gen.dim >= 6 && np->gen.rmax > L) L = np->gen.rmax;

    time1 = NpicGetTime();

    switch (np->type) {

        case NPIC_IMAGE_2L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_2l *p = NpicCastImage(np);
            Npic_l Cv[L];
            int y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(y, x, Cv)
            {
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (y = 0; y < p->ymax; y++)
            {
                /* Copy column in Cv */
                for (x = 0; x < p->xmax; x++)
                    Cv[x] = p->pix[y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (x = 0; x < p->xmax; x++)
                    p->pix[y][x] = Npic_Saito_max_rev (Cv, x, p->xmax);
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (y = 0; y < p->ymax; y++)
                    Cv[y] = p->pix[y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (y = 0; y < p->ymax; y++)
                    p->pix[y][x] = Npic_Saito_max_rev (Cv, y, p->ymax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

        case NPIC_IMAGE_3L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_3l *p = NpicCastImage(np);
            Npic_l Cv[L];
            int z, y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(z, y, x, Cv)
            {
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            {
                /* Copy column in Cv */
                for (x = 0; x < p->xmax; x++)
                    Cv[x] = p->pix[z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (x = 0; x < p->xmax; x++)
                    p->pix[z][y][x] = Npic_Saito_max_rev (Cv, x, p->xmax);
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (z = 0; z < p->zmax; z++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (y = 0; y < p->ymax; y++)
                    Cv[y] = p->pix[z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (y = 0; y < p->ymax; y++)
                    p->pix[z][y][x] = Npic_Saito_max_rev (Cv, y, p->ymax);
            }

            /*
             * Scan in z
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (z = 0; z < p->zmax; z++)
                    Cv[z] = p->pix[z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (z = 0; z < p->zmax; z++)
                    p->pix[z][y][x] = Npic_Saito_max_rev (Cv, z, p->zmax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

        case NPIC_IMAGE_4L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_4l *p = NpicCastImage(np);
            Npic_l Cv[L];
            int t, z, y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(t, z, y, x, Cv)
            {
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            {
                /* Copy column in Cv */
                for (x = 0; x < p->xmax; x++)
                    Cv[x] = p->pix[t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (x = 0; x < p->xmax; x++)
                    p->pix[t][z][y][x] = Npic_Saito_max_rev (Cv, x, p->xmax);
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (y = 0; y < p->ymax; y++)
                    Cv[y] = p->pix[t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (y = 0; y < p->ymax; y++)
                    p->pix[t][z][y][x] = Npic_Saito_max_rev (Cv, y, p->ymax);
            }

            /*
             * Scan in z
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (t = 0; t < p->tmax; t++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (z = 0; z < p->zmax; z++)
                    Cv[z] = p->pix[t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (z = 0; z < p->zmax; z++)
                    p->pix[t][z][y][x] = Npic_Saito_max_rev (Cv, z, p->zmax);
            }

            /*
             * Scan in t
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (t = 0; t < p->tmax; t++)
                    Cv[t] = p->pix[t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (t = 0; t < p->tmax; t++)
                    p->pix[t][z][y][x] = Npic_Saito_max_rev (Cv, t, p->tmax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

        case NPIC_IMAGE_5L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_5l *p = NpicCastImage(np);
            Npic_l Cv[L];
            int s, t, z, y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(s, t, z, y, x, Cv)
            {
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            {
                /* Copy column in Cv */
                for (x = 0; x < p->xmax; x++)
                    Cv[x] = p->pix[s][t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (x = 0; x < p->xmax; x++)
                    p->pix[s][t][z][y][x] = Npic_Saito_max_rev (Cv, x, p->xmax);
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (y = 0; y < p->ymax; y++)
                    Cv[y] = p->pix[s][t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (y = 0; y < p->ymax; y++)
                    p->pix[s][t][z][y][x] = Npic_Saito_max_rev (Cv, y, p->ymax);
            }

            /*
             * Scan in z
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (z = 0; z < p->zmax; z++)
                    Cv[z] = p->pix[s][t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (z = 0; z < p->zmax; z++)
                    p->pix[s][t][z][y][x] = Npic_Saito_max_rev (Cv, z, p->zmax);
            }

            /*
             * Scan in t
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (t = 0; t < p->tmax; t++)
                    Cv[t] = p->pix[s][t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (t = 0; t < p->tmax; t++)
                    p->pix[s][t][z][y][x] = Npic_Saito_max_rev (Cv, t, p->tmax);
            }

            /*
             * Scan in s
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (s = 0; s < p->smax; s++)
                    Cv[s] = p->pix[s][t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (s = 0; s < p->smax; s++)
                    p->pix[s][t][z][y][x] = Npic_Saito_max_rev (Cv, s, p->smax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

        case NPIC_IMAGE_6L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_6l *p = NpicCastImage(np);
            Npic_l Cv[L];
            int r, s, t, z, y, x;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel private(r, s, t, z, y, x, Cv)
            {
#           endif

            /*
             * First scan in x
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            {
                /* Copy column in Cv */
                for (x = 0; x < p->xmax; x++)
                    Cv[x] = p->pix[r][s][t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (x = 0; x < p->xmax; x++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_max_rev (Cv, x, p->xmax);
            }

            /*
             * Scan in y
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (y = 0; y < p->ymax; y++)
                    Cv[y] = p->pix[r][s][t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (y = 0; y < p->ymax; y++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_max_rev (Cv, y, p->ymax);
            }

            /*
             * Scan in z
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (z = 0; z < p->zmax; z++)
                    Cv[z] = p->pix[r][s][t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (z = 0; z < p->zmax; z++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_max_rev (Cv, z, p->zmax);
            }

            /*
             * Scan in t
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (s = 0; s < p->smax; s++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (t = 0; t < p->tmax; t++)
                    Cv[t] = p->pix[r][s][t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (t = 0; t < p->tmax; t++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_max_rev (Cv, t, p->tmax);
            }

            /*
             * Scan in s
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (r = 0; r < p->rmax; r++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (s = 0; s < p->smax; s++)
                    Cv[s] = p->pix[r][s][t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (s = 0; s < p->smax; s++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_max_rev (Cv, s, p->smax);
            }

            /*
             * Scan in r
            */
#           ifdef NPIC_USE_OMP
#           pragma omp for schedule(dynamic)
#           endif
            for (s = 0; s < p->smax; s++)
            for (t = 0; t < p->tmax; t++)
            for (z = 0; z < p->zmax; z++)
            for (y = 0; y < p->ymax; y++)
            for (x = 0; x < p->xmax; x++)
            {
                /* Copy column in Cv */
                for (r = 0; r < p->rmax; r++)
                    Cv[r] = p->pix[r][s][t][z][y][x];

                /* Compute max on column Cv and store result in p->pix */
                for (r = 0; r < p->rmax; r++)
                    p->pix[r][s][t][z][y][x] = Npic_Saito_max_rev (Cv, r, p->rmax);
            }

#           ifdef NPIC_USE_OMP
            }  /* omp parallel */
#           endif
            break;
        }

        default :
            np->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 ---------------------*/

/*
 * Return min dE from pos on column Cv[-1..cmax]
 * Used in NpicSEDT_Saito.
*/

Npic_l Npic_Saito_min (Npic_l *Cv, int pos, int cmax)
{
    Npic_l k, k2, h, j = Cv[pos];
            
    for (k = 1; pos-k >= -1; k++) {
        k2 = k*k; if (k2 >= j) break;
        h = Cv[pos-k] + k2; if (h < j) j = h;
    }
            
    for (k = 1; pos+k <= cmax; k++) {
        k2 = k*k; if (k2 >= j) break;
        h = Cv[pos+k] + k2; if (h < j) j = h;
    }
            
    return j;
}


/*
 * Return min dE from pos on column Cv[0..cmax-1] ; -1 means no distance.
 * Used in NpicSEDT_Saito_inf.
*/

Npic_l Npic_Saito_min_inf (Npic_l *Cv, int pos, int cmax)
{
    Npic_l k, k2, h, j = Cv[pos];
    
    for (k = 1; pos-k >= 0; k++) {
        k2 = k*k; if (k2 >= j && j != -1) break;
        h = Cv[pos-k]; if (h == -1) continue;
        h += k2; if (h < j || j == -1) j = h;
    }
            
    for (k = 1; pos+k < cmax; k++) {
        k2 = k*k; if (k2 >= j && j != -1) break;
        h = Cv[pos+k]; if (h == -1) continue;
        h += k2; if (h < j || j == -1) j = h;
    }
    
    return j;
}


/*
 * Return max dE from pos on column Cv[0..cmax-1].
 * Used in NpicSEDT_Saito_rev.
*/

Npic_l Npic_Saito_max_rev (Npic_l *Cv, int pos, int cmax)
{
    Npic_l k, h, j = Cv[pos];
    
    for (k = 1; pos-k >= 0; k++) {
        h = Cv[pos-k]; if (h == 0) continue;
        h -= k*k; if (h > j) j = h;
    }
            
    for (k = 1; pos+k < cmax; k++) {
        h = Cv[pos+k]; if (h == 0) continue;
        h -= k*k; if (h > j) j = h;
    }
    
    return j;
}


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