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

/*
 * dist_sedt_quadra.c - 25/09/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.
 * Quadratic algorithm, very slow, provided for checks.
 *
 * Do nothing if np is not ok. set not ok on error. Verbose.
 * Return computation time in second.microsecond
*/

double NpicSEDT_quadra (Npic_image *np)
{
    double time1, time2;

    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;
    }

    time1 = NpicGetTime();

    switch (np->type) {

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

            Npic_image_2l *p = NpicCastImage (np);
            int y1, x1, y2, x2;
            Npic_l j, k;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for schedule(dynamic) \
                       private(y1, x1, y2, x2, j, k)
#           endif
            for (y1 = 0; y1 < p->ymax; y1++)
            for (x1 = 0; x1 < p->xmax; x1++)
            {
                if (p->pix[y1][x1] == 0) continue;
                
                k =  p->ymax*p->ymax + p->xmax*p->xmax;
                
                /* Search also in border pixels */
                for (y2 = -1; y2 <= p->ymax; y2++)
                for (x2 = -1; x2 <= p->xmax; x2++)
                if (p->pix[y2][x2] == 0)
                {
                    j = (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (j < k) k = j;
                }
                
                p->pix[y1][x1] = k;
            }
        
            break;
        }

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

            Npic_image_3l *p = NpicCastImage (np);
            int z1, y1, x1, z2, y2, x2;
            Npic_l j, k;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for schedule(dynamic) \
                       private(z1, y1, x1, z2, y2, x2, j, k)
#           endif
            for (z1 = 0; z1 < p->zmax; z1++)
            for (y1 = 0; y1 < p->ymax; y1++)
            for (x1 = 0; x1 < p->xmax; x1++)
            {
                if (p->pix[z1][y1][x1] == 0) continue;
                
                k = p->zmax*p->zmax + p->ymax*p->ymax + p->xmax*p->xmax;
                
                /* Search also in border pixels */
                for (z2 = -1; z2 <= p->zmax; z2++)
                for (y2 = -1; y2 <= p->ymax; y2++)
                for (x2 = -1; x2 <= p->xmax; x2++)
                if (p->pix[z2][y2][x2] == 0)
                {
                    j = (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (j < k) k = j;
                }
                
                p->pix[z1][y1][x1] = k;
            }
        
            break;
        }

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

            Npic_image_4l *p = NpicCastImage (np);
            int t1, z1, y1, x1, t2, z2, y2, x2;
            Npic_l j, k;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for schedule(dynamic) \
                       private(t1, z1, y1, x1, t2, z2, y2, x2, j, k)
#           endif
            for (t1 = 0; t1 < p->tmax; t1++)
            for (z1 = 0; z1 < p->zmax; z1++)
            for (y1 = 0; y1 < p->ymax; y1++)
            for (x1 = 0; x1 < p->xmax; x1++)
            {
                if (p->pix[t1][z1][y1][x1] == 0) continue;
                
                k = p->tmax*p->tmax 
                  + p->zmax*p->zmax + p->ymax*p->ymax + p->xmax*p->xmax;
                
                /* Search also in border pixels */
                for (t2 = -1; t2 <= p->tmax; t2++)
                for (z2 = -1; z2 <= p->zmax; z2++)
                for (y2 = -1; y2 <= p->ymax; y2++)
                for (x2 = -1; x2 <= p->xmax; x2++)
                if (p->pix[t2][z2][y2][x2] == 0)
                {
                    j = (t2-t1)*(t2-t1) 
                      + (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (j < k) k = j;
                }
                
                p->pix[t1][z1][y1][x1] = k;
            }
        
            break;
        }

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

            Npic_image_5l *p = NpicCastImage (np);
            int s1, t1, z1, y1, x1, s2, t2, z2, y2, x2;
            Npic_l j, k;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for schedule(dynamic) \
                       private(s1, t1, z1, y1, x1, s2, t2, z2, y2, x2, j, k)
#           endif
            for (s1 = 0; s1 < p->smax; s1++)
            for (t1 = 0; t1 < p->tmax; t1++)
            for (z1 = 0; z1 < p->zmax; z1++)
            for (y1 = 0; y1 < p->ymax; y1++)
            for (x1 = 0; x1 < p->xmax; x1++)
            {
                if (p->pix[s1][t1][z1][y1][x1] == 0) continue;
                
                k = p->smax*p->smax + p->tmax*p->tmax 
                  + p->zmax*p->zmax + p->ymax*p->ymax + p->xmax*p->xmax;
                
                /* Search also in border pixels */
                for (s2 = -1; s2 <= p->smax; s2++)
                for (t2 = -1; t2 <= p->tmax; t2++)
                for (z2 = -1; z2 <= p->zmax; z2++)
                for (y2 = -1; y2 <= p->ymax; y2++)
                for (x2 = -1; x2 <= p->xmax; x2++)
                if (p->pix[s2][t2][z2][y2][x2] == 0)
                {
                    j = (s2-s1)*(s2-s1) + (t2-t1)*(t2-t1) 
                      + (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (j < k) k = j;
                }
                
                p->pix[s1][t1][z1][y1][x1] = k;
            }
        
            break;
        }

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

            Npic_image_6l *p = NpicCastImage (np);
            int r1, s1, t1, z1, y1, x1, r2, s2, t2, z2, y2, x2;
            Npic_l j, k;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for schedule(dynamic) \
                       private(r1, s1, t1, z1, y1, x1, r2, s2, t2, z2, y2, x2, j, k)
#           endif
            for (r1 = 0; r1 < p->rmax; r1++)
            for (s1 = 0; s1 < p->smax; s1++)
            for (t1 = 0; t1 < p->tmax; t1++)
            for (z1 = 0; z1 < p->zmax; z1++)
            for (y1 = 0; y1 < p->ymax; y1++)
            for (x1 = 0; x1 < p->xmax; x1++)
            {
                if (p->pix[r1][s1][t1][z1][y1][x1] == 0) continue;
                
                k = p->rmax*p->rmax + p->smax*p->smax + p->tmax*p->tmax 
                  + p->zmax*p->zmax + p->ymax*p->ymax + p->xmax*p->xmax;
                
                /* Search also in border pixels */
                for (r2 = -1; r2 <= p->rmax; r2++)
                for (s2 = -1; s2 <= p->smax; s2++)
                for (t2 = -1; t2 <= p->tmax; t2++)
                for (z2 = -1; z2 <= p->zmax; z2++)
                for (y2 = -1; y2 <= p->ymax; y2++)
                for (x2 = -1; x2 <= p->xmax; x2++)
                if (p->pix[r2][s2][t2][z2][y2][x2] == 0)
                {
                    j = (r2-r1)*(r2-r1) + (s2-s1)*(s2-s1) + (t2-t1)*(t2-t1) 
                      + (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (j < k) k = j;
                }
                
                p->pix[r1][s1][t1][z1][y1][x1] = k;
            }
        
            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.
 * Quadratic algorithm, very slow, provided for checks.
 *
 * Do nothing if np is not ok. set not ok on error. Verbose.
 * Return computation time in second.microsecond
*/

double NpicSEDT_quadra_inf (Npic_image *np)
{
    double time1, time2;

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

    time1 = NpicGetTime();

    switch (np->type) {

        case NPIC_IMAGE_2L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_2l *p = NpicCastImage (np);
            int y1, x1, y2, x2;
            Npic_l j, k;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for schedule(dynamic) \
                       private(y1, x1, y2, x2, j, k)
#           endif
            for (y1 = 0; y1 < p->ymax; y1++)
            for (x1 = 0; x1 < p->xmax; x1++)
            {
                if (p->pix[y1][x1] == 0) continue;
                
                k = -1;
                for (y2 = 0; y2 < p->ymax; y2++)
                for (x2 = 0; x2 < p->xmax; x2++)
                if (p->pix[y2][x2] == 0)
                {
                    j = (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (j < k || k == -1) k = j;
                }
                
                p->pix[y1][x1] = k;
            }

            break;
        }

        case NPIC_IMAGE_3L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_3l *p = NpicCastImage (np);
            int z1, y1, x1, z2, y2, x2;
            Npic_l j, k;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for schedule(dynamic) \
                       private(z1, y1, x1, z2, y2, x2, j, k)
#           endif
            for (z1 = 0; z1 < p->zmax; z1++)
            for (y1 = 0; y1 < p->ymax; y1++)
            for (x1 = 0; x1 < p->xmax; x1++)
            {
                if (p->pix[z1][y1][x1] == 0) continue;
                
                k = -1;
                for (z2 = 0; z2 < p->zmax; z2++)
                for (y2 = 0; y2 < p->ymax; y2++)
                for (x2 = 0; x2 < p->xmax; x2++)
                if (p->pix[z2][y2][x2] == 0)
                {
                    j = (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (j < k || k == -1) k = j;
                }
                
                p->pix[z1][y1][x1] = k;
            }

            break;
        }

        case NPIC_IMAGE_4L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_4l *p = NpicCastImage (np);
            int t1, z1, y1, x1, t2, z2, y2, x2;
            Npic_l j, k;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for schedule(dynamic) \
                       private(t1, z1, y1, x1, t2, z2, y2, x2, j, k)
#           endif
            for (t1 = 0; t1 < p->tmax; t1++)
            for (z1 = 0; z1 < p->zmax; z1++)
            for (y1 = 0; y1 < p->ymax; y1++)
            for (x1 = 0; x1 < p->xmax; x1++)
            {
                if (p->pix[t1][z1][y1][x1] == 0) continue;
                
                k = -1;
                for (t2 = 0; t2 < p->tmax; t2++)
                for (z2 = 0; z2 < p->zmax; z2++)
                for (y2 = 0; y2 < p->ymax; y2++)
                for (x2 = 0; x2 < p->xmax; x2++)
                if (p->pix[t2][z2][y2][x2] == 0)
                {
                    j = (t2-t1)*(t2-t1) 
                      + (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (j < k || k == -1) k = j;
                }
                
                p->pix[t1][z1][y1][x1] = k;
            }

            break;
        }

        case NPIC_IMAGE_5L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_5l *p = NpicCastImage (np);
            int s1, t1, z1, y1, x1, s2, t2, z2, y2, x2;
            Npic_l j, k;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for schedule(dynamic) \
                       private(s1, t1, z1, y1, x1, s2, t2, z2, y2, x2, j, k)
#           endif
            for (s1 = 0; s1 < p->smax; s1++)
            for (t1 = 0; t1 < p->tmax; t1++)
            for (z1 = 0; z1 < p->zmax; z1++)
            for (y1 = 0; y1 < p->ymax; y1++)
            for (x1 = 0; x1 < p->xmax; x1++)
            {
                if (p->pix[s1][t1][z1][y1][x1] == 0) continue;
                
                k = -1;
                for (s2 = 0; s2 < p->smax; s2++)
                for (t2 = 0; t2 < p->tmax; t2++)
                for (z2 = 0; z2 < p->zmax; z2++)
                for (y2 = 0; y2 < p->ymax; y2++)
                for (x2 = 0; x2 < p->xmax; x2++)
                if (p->pix[s2][t2][z2][y2][x2] == 0)
                {
                    j = (s2-s1)*(s2-s1) + (t2-t1)*(t2-t1) 
                      + (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (j < k || k == -1) k = j;
                }
                
                p->pix[s1][t1][z1][y1][x1] = k;
            }

            break;
        }

        case NPIC_IMAGE_6L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_6l *p = NpicCastImage (np);
            int r1, s1, t1, z1, y1, x1, r2, s2, t2, z2, y2, x2;
            Npic_l j, k;

#           ifdef NPIC_USE_OMP
#           pragma omp parallel for schedule(dynamic) \
                       private(r1, s1, t1, z1, y1, x1, r2, s2, t2, z2, y2, x2, j, k)
#           endif
            for (r1 = 0; r1 < p->rmax; r1++)
            for (s1 = 0; s1 < p->smax; s1++)
            for (t1 = 0; t1 < p->tmax; t1++)
            for (z1 = 0; z1 < p->zmax; z1++)
            for (y1 = 0; y1 < p->ymax; y1++)
            for (x1 = 0; x1 < p->xmax; x1++)
            {
                if (p->pix[r1][s1][t1][z1][y1][x1] == 0) continue;
                
                k = -1;
                for (r2 = 0; r2 < p->rmax; r2++)
                for (s2 = 0; s2 < p->smax; s2++)
                for (t2 = 0; t2 < p->tmax; t2++)
                for (z2 = 0; z2 < p->zmax; z2++)
                for (y2 = 0; y2 < p->ymax; y2++)
                for (x2 = 0; x2 < p->xmax; x2++)
                if (p->pix[r2][s2][t2][z2][y2][x2] == 0)
                {
                    j = (r2-r1)*(r2-r1) + (s2-s1)*(s2-s1) + (t2-t1)*(t2-t1) 
                      + (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (j < k || k == -1) k = j;
                }
                
                p->pix[r1][s1][t1][z1][y1][x1] = k;
            }

            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.
 * Quadratic algorithm, very slow, provided for checks.
 *
 * Do nothing if np is not ok. set not ok on error. Verbose.
 * Return computation time in second.microsecond
*/

double NpicSEDT_quadra_rev (Npic_image *np)
{
    Npic_image *np2;
    double time1, time2;

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

    np2 = NpicDupImage (np);
    NpicCopyImageI (np2, np);
 
    time1 = NpicGetTime();

    switch (np->type) {

        case NPIC_IMAGE_2L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_2l *p1 = NpicCastImage (np),
                          *p2 = NpicCastImage (np2);
            int y1, x1, y2, x2;
            Npic_l j;

#           ifdef NPIC_USE_OMP
            Npic_l k;

#           pragma omp parallel for schedule(dynamic) \
                       private(y1, x1, y2, x2, j, k) 
            for (y1 = 0; y1 < p2->ymax; y1++)
            for (x1 = 0; x1 < p2->xmax; x1++)
            {
                j = 0;

                for (y2 = 0; y2 < p2->ymax; y2++)
                for (x2 = 0; x2 < p2->xmax; x2++)
                {
                    k = p2->pix[y2][x2];
                    if (k <= j) continue;
                    k -= (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (k > j) j = k;
                }

                p1->pix[y1][x1] = j;
            }

#           else  /* NPIC_USE_OMP */

            /* This code is faster but not parallelisable */
            for (y1 = 0; y1 < p2->ymax; y1++)
            for (x1 = 0; x1 < p2->xmax; x1++)
            {
                if (p2->pix[y1][x1] == 0) continue;
                
                for (y2 = 0; y2 < p2->ymax; y2++)
                for (x2 = 0; x2 < p2->xmax; x2++)
                {
                    j = p2->pix[y1][x1] 
                      - (y2-y1)*(y2-y1) - (x2-x1)*(x2-x1);
                    if (j > p1->pix[y2][x2]) p1->pix[y2][x2] = j;
                }   
            }
#           endif  /* NPIC_USE_OMP */

            break;
        }

        case NPIC_IMAGE_3L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_3l *p1 = NpicCastImage (np),
                          *p2 = NpicCastImage (np2);
            int z1, y1, x1, z2, y2, x2;
            Npic_l j;

#           ifdef NPIC_USE_OMP
            Npic_l k;

#           pragma omp parallel for schedule(dynamic) \
                       private(z1, y1, x1, z2, y2, x2, j, k) 
            for (z1 = 0; z1 < p2->zmax; z1++)
            for (y1 = 0; y1 < p2->ymax; y1++)
            for (x1 = 0; x1 < p2->xmax; x1++)
            {
                j = 0;

                for (z2 = 0; z2 < p2->zmax; z2++)
                for (y2 = 0; y2 < p2->ymax; y2++)
                for (x2 = 0; x2 < p2->xmax; x2++)
                {
                    k = p2->pix[z2][y2][x2];
                    if (k <= j) continue;
                    k -= (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (k > j) j = k;
                }

                p1->pix[z1][y1][x1] = j;
            }

#           else  /* NPIC_USE_OMP */

            /* This code is faster but not parallelisable */
            for (z1 = 0; z1 < p2->zmax; z1++)
            for (y1 = 0; y1 < p2->ymax; y1++)
            for (x1 = 0; x1 < p2->xmax; x1++)
            {
                if (p2->pix[z1][y1][x1] == 0) continue;
                
                for (z2 = 0; z2 < p2->zmax; z2++)
                for (y2 = 0; y2 < p2->ymax; y2++)
                for (x2 = 0; x2 < p2->xmax; x2++)
                {
                    j = p2->pix[z1][y1][x1] 
                      - (z2-z1)*(z2-z1) - (y2-y1)*(y2-y1) - (x2-x1)*(x2-x1);
                    if (j > p1->pix[z2][y2][x2]) p1->pix[z2][y2][x2] = j;
                }   
            }
#           endif  /* NPIC_USE_OMP */

            break;
        }

        case NPIC_IMAGE_4L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_4l *p1 = NpicCastImage (np),
                          *p2 = NpicCastImage (np2);
            int t1, z1, y1, x1, t2, z2, y2, x2;
            Npic_l j;

#           ifdef NPIC_USE_OMP
            Npic_l k;

#           pragma omp parallel for schedule(dynamic) \
                       private(t1, z1, y1, x1, t2, z2, y2, x2, j, k) 
            for (t1 = 0; t1 < p2->tmax; t1++)
            for (z1 = 0; z1 < p2->zmax; z1++)
            for (y1 = 0; y1 < p2->ymax; y1++)
            for (x1 = 0; x1 < p2->xmax; x1++)
            {
                j = 0;

                for (t2 = 0; t2 < p2->tmax; t2++)
                for (z2 = 0; z2 < p2->zmax; z2++)
                for (y2 = 0; y2 < p2->ymax; y2++)
                for (x2 = 0; x2 < p2->xmax; x2++)
                {
                    k = p2->pix[t2][z2][y2][x2];
                    if (k <= j) continue;
                    k -= (t2-t1)*(t2-t1) 
                       + (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (k > j) j = k;
                }

                p1->pix[t1][z1][y1][x1] = j;
            }

#           else  /* NPIC_USE_OMP */

            /* This code is faster but not parallelisable */
            for (t1 = 0; t1 < p2->tmax; t1++)
            for (z1 = 0; z1 < p2->zmax; z1++)
            for (y1 = 0; y1 < p2->ymax; y1++)
            for (x1 = 0; x1 < p2->xmax; x1++)
            {
                if (p2->pix[t1][z1][y1][x1] == 0) continue;
                
                for (t2 = 0; t2 < p2->tmax; t2++)
                for (z2 = 0; z2 < p2->zmax; z2++)
                for (y2 = 0; y2 < p2->ymax; y2++)
                for (x2 = 0; x2 < p2->xmax; x2++)
                {
                    j = p2->pix[t1][z1][y1][x1] 
                      - (t2-t1)*(t2-t1)
                      - (z2-z1)*(z2-z1) - (y2-y1)*(y2-y1) - (x2-x1)*(x2-x1);
                    if (j > p1->pix[t2][z2][y2][x2]) p1->pix[t2][z2][y2][x2] = j;
                }   
            }
#           endif  /* NPIC_USE_OMP */

            break;
        }

        case NPIC_IMAGE_5L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_5l *p1 = NpicCastImage (np),
                          *p2 = NpicCastImage (np2);
            int s1, t1, z1, y1, x1, s2, t2, z2, y2, x2;
            Npic_l j;

#           ifdef NPIC_USE_OMP
            Npic_l k;

#           pragma omp parallel for schedule(dynamic) \
                       private(s1, t1, z1, y1, x1, s2, t2, z2, y2, x2, j, k) 
            for (s1 = 0; s1 < p2->smax; s1++)
            for (t1 = 0; t1 < p2->tmax; t1++)
            for (z1 = 0; z1 < p2->zmax; z1++)
            for (y1 = 0; y1 < p2->ymax; y1++)
            for (x1 = 0; x1 < p2->xmax; x1++)
            {
                j = 0;

                for (s2 = 0; s2 < p2->smax; s2++)
                for (t2 = 0; t2 < p2->tmax; t2++)
                for (z2 = 0; z2 < p2->zmax; z2++)
                for (y2 = 0; y2 < p2->ymax; y2++)
                for (x2 = 0; x2 < p2->xmax; x2++)
                {
                    k = p2->pix[s2][t2][z2][y2][x2];
                    if (k <= j) continue;
                    k -= (s2-s1)*(s2-s1) + (t2-t1)*(t2-t1) 
                       + (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (k > j) j = k;
                }

                p1->pix[s1][t1][z1][y1][x1] = j;
            }

#           else  /* NPIC_USE_OMP */

            /* This code is faster but not parallelisable */
            for (s1 = 0; s1 < p2->smax; s1++)
            for (t1 = 0; t1 < p2->tmax; t1++)
            for (z1 = 0; z1 < p2->zmax; z1++)
            for (y1 = 0; y1 < p2->ymax; y1++)
            for (x1 = 0; x1 < p2->xmax; x1++)
            {
                if (p2->pix[s1][t1][z1][y1][x1] == 0) continue;
                
                for (s2 = 0; s2 < p2->smax; s2++)
                for (t2 = 0; t2 < p2->tmax; t2++)
                for (z2 = 0; z2 < p2->zmax; z2++)
                for (y2 = 0; y2 < p2->ymax; y2++)
                for (x2 = 0; x2 < p2->xmax; x2++)
                {
                    j = p2->pix[s1][t1][z1][y1][x1] 
                      - (s2-s1)*(s2-s1) - (t2-t1)*(t2-t1)
                      - (z2-z1)*(z2-z1) - (y2-y1)*(y2-y1) - (x2-x1)*(x2-x1);
                    if (j > p1->pix[s2][t2][z2][y2][x2]) p1->pix[s2][t2][z2][y2][x2] = j;
                }   
            }
#           endif  /* NPIC_USE_OMP */

            break;
        }

        case NPIC_IMAGE_6L : {  /* - - - - - - - - - - - - - - - - - */
        
            Npic_image_6l *p1 = NpicCastImage (np),
                          *p2 = NpicCastImage (np2);
            int r1, s1, t1, z1, y1, x1, r2, s2, t2, z2, y2, x2;
            Npic_l j;

#           ifdef NPIC_USE_OMP
            Npic_l k;

#           pragma omp parallel for schedule(dynamic) \
                       private(r1, s1, t1, z1, y1, x1, r2, s2, t2, z2, y2, x2, j, k) 
            for (r1 = 0; r1 < p2->rmax; r1++)
            for (s1 = 0; s1 < p2->smax; s1++)
            for (t1 = 0; t1 < p2->tmax; t1++)
            for (z1 = 0; z1 < p2->zmax; z1++)
            for (y1 = 0; y1 < p2->ymax; y1++)
            for (x1 = 0; x1 < p2->xmax; x1++)
            {
                j = 0;

                for (r2 = 0; r2 < p2->rmax; r2++)
                for (s2 = 0; s2 < p2->smax; s2++)
                for (t2 = 0; t2 < p2->tmax; t2++)
                for (z2 = 0; z2 < p2->zmax; z2++)
                for (y2 = 0; y2 < p2->ymax; y2++)
                for (x2 = 0; x2 < p2->xmax; x2++)
                {
                    k = p2->pix[r2][s2][t2][z2][y2][x2];
                    if (k <= j) continue;
                    k -= (r2-r1)*(r2-r1) + (s2-s1)*(s2-s1) + (t2-t1)*(t2-t1) 
                       + (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1);
                    if (k > j) j = k;
                }

                p1->pix[r1][s1][t1][z1][y1][x1] = j;
            }

#           else  /* NPIC_USE_OMP */

            /* This code is faster but not parallelisable */
            for (r1 = 0; r1 < p2->rmax; r1++)
            for (s1 = 0; s1 < p2->smax; s1++)
            for (t1 = 0; t1 < p2->tmax; t1++)
            for (z1 = 0; z1 < p2->zmax; z1++)
            for (y1 = 0; y1 < p2->ymax; y1++)
            for (x1 = 0; x1 < p2->xmax; x1++)
            {
                if (p2->pix[r1][s1][t1][z1][y1][x1] == 0) continue;
                
                for (r2 = 0; r2 < p2->rmax; r2++)
                for (s2 = 0; s2 < p2->smax; s2++)
                for (t2 = 0; t2 < p2->tmax; t2++)
                for (z2 = 0; z2 < p2->zmax; z2++)
                for (y2 = 0; y2 < p2->ymax; y2++)
                for (x2 = 0; x2 < p2->xmax; x2++)
                {
                    j = p2->pix[r1][s1][t1][z1][y1][x1] 
                      - (r2-r1)*(r2-r1) - (s2-s1)*(s2-s1) - (t2-t1)*(t2-t1)
                      - (z2-z1)*(z2-z1) - (y2-y1)*(y2-y1) - (x2-x1)*(x2-x1);
                    if (j > p1->pix[r2][s2][t2][z2][y2][x2]) p1->pix[r2][s2][t2][z2][y2][x2] = j;
                }   
            }
#           endif  /* NPIC_USE_OMP */

            break;
        }

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

    time2 = NpicGetTime();

    NpicDestroyImage (np2);
    return time2-time1;
}


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


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