/*
 * The Npic library and tools
 *
 * Copyright (C) 2003 Edouard Thiel <Edouard.Thiel@lif.univ-mrs.fr>
 *
 * This program 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 "npic-wd-norm.ct" */

/*
 * npic-wd-norm.c - 09/01/2009
 *
*/

#include <npic.h>


void ShowUsage ()
{
    printf (
        "npic-wd-norm - Test if a weighted (= chamfer) distance induces a norm.\n"
        "Usage:\n"
        "  npic-wd-norm  -h | -help | --help   : print help\n"
        "  npic-wd-norm  mask r [-o out1] [-p j] : \n"
        "      compute a DT_inf on a (2r+1)^dim centred image, then scan for\n"
        "      homogeneity and positive definiteness with a sieve on visible vectors.\n"
        "\n"
        "  -o out1 : save DT_inf  (default is don't save). Wrong vectors\n"
        "            are marked as negative.\n"
        "  -p j    : print the j first series of wrong vectors (default is 10).\n"
        "  -p -1   : silent mode, only print \"NORM: YES\" or \"NORM: NO\".\n"
        "\n"
        "mask : mask file in " NPIC_KNOWN_MASK_EXT " format (see npic-mask -h)\n"
        "out1 : image file in " NPIC_KNOWN_IMAGE_EXT " format\n"
        "\n"
        "Example: 5-7-d16 is not a 2D norm\n"
        "  ./npic-wd-norm ../masks/d-5-7-d16_2l.nmask 10 -o tmp1.npz\n"
        "  ./npic-print tmp1.npz\n"
        "\n"
        "Example: bash script that checks a,b,c 3D masks\n"
        "  for ((a = 1; a < 10; a++)); do\n"
        "  for ((b = 1; b < 10; b++)); do\n" 
        "  for ((c = 1; c < 10; c++)); do\n"
        "      echo -n \"$a $b $c  \"\n"
        "      ./npic-mask tmp1.nmask -create -3l -dist -ghalf \\\n"
        "          -vec 0:0:1:$a -vec 0:1:1:$b -vec 1:1:1:$c > /dev/null\n"
        "      ./npic-wd-norm tmp1.nmask 20 -p -1\n"
        "  done ; done ; done | more\n"
    );
}


void sieve_draw_image (Npic_mask *m1, Npic_image *np1, int sieve_r)
{
    switch (np1->type) {
        case NPIC_IMAGE_2L : 
            NpicFillVal_i (np1, 1); 
            np1->im_2l.pix[sieve_r][sieve_r] = 0;
            break;
        case NPIC_IMAGE_3L : 
            NpicFillVal_i (np1, 1); 
            np1->im_3l.pix[sieve_r][sieve_r][sieve_r] = 0;
            break;
        case NPIC_IMAGE_4L : 
            NpicFillVal_i (np1, 1); 
            np1->im_4l.pix[sieve_r][sieve_r][sieve_r][sieve_r] = 0;
            break;
        case NPIC_IMAGE_5L : 
            NpicFillVal_i (np1, 1); 
            np1->im_5l.pix[sieve_r][sieve_r][sieve_r][sieve_r][sieve_r] = 0;
            break;
        case NPIC_IMAGE_6L : 
            NpicFillVal_i (np1, 1); 
            np1->im_6l.pix[sieve_r][sieve_r][sieve_r][sieve_r][sieve_r][sieve_r] = 0;
            break;
        case NPIC_IMAGE_2D : 
            NpicFillVal_d (np1, 1.0); 
            np1->im_2d.pix[sieve_r][sieve_r] = 0.0;
            break;
        case NPIC_IMAGE_3D : 
            NpicFillVal_d (np1, 1.0); 
            np1->im_3d.pix[sieve_r][sieve_r][sieve_r] = 0.0;
            break;
        case NPIC_IMAGE_4D : 
            NpicFillVal_d (np1, 1.0); 
            np1->im_4d.pix[sieve_r][sieve_r][sieve_r][sieve_r] = 0.0;
            break;
        case NPIC_IMAGE_5D : 
            NpicFillVal_d (np1, 1.0); 
            np1->im_5d.pix[sieve_r][sieve_r][sieve_r][sieve_r][sieve_r] = 0.0;
            break;
        case NPIC_IMAGE_6D : 
            NpicFillVal_d (np1, 1.0); 
            np1->im_6d.pix[sieve_r][sieve_r][sieve_r][sieve_r][sieve_r][sieve_r] = 0.0;
            break;
    }
}


#define D_THRESH 1E-8

int sieve_scan (Npic_image *np1, Npic_image *np2, int sieve_r, int sieve_j)
{
    int nj = 0, nk, k;

    switch (np1->type)
    {

        case NPIC_IMAGE_2L : {
            Npic_image_2l *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int y, x, yd, xd,
                yb, xb, yk, xk;

            /* Orthant by orthant */
            for (yd = 1; yd >= -1; yd -=2)
            for (xd = 1; xd >= -1; xd -=2)
            for (y = 0; y <= sieve_r; y++)
            for (x = 0; x <= sieve_r; x++)
            {
                yb = sieve_r + y*yd;
                xb = sieve_r + x*xd;

                /* Visible point ? */
                if (p2->pix[yb][xb] != 1) continue;

                /* Periods */
                nk = 0;
                for (k = 2; 
                            y*k <= sieve_r && 
                            x*k <= sieve_r ; k++) 
                {
                    yk = sieve_r + y*yd*k;
                    xk = sieve_r + x*xd*k;
                    p2->pix[yk][xk] = 0;  /* mark as non visible */

                    /* Test positive definiteness and homogeneity */
                    if (p1->pix[yb][xb] <= 0 || p1->pix[yb][xb] * k != p1->pix[yk][xk])
                    {
                        nk++; if (nk == 1) nj++;
                        if (nj <= sieve_j) { 
                            if (nk == 1) printf ("  [%d,%d] = %" NPIC_PL, 
                                yb-sieve_r, xb-sieve_r, 
                                p1->pix[yb][xb] );
                            printf ("  [%d,%d] = %" NPIC_PL,
                                yk-sieve_r, xk-sieve_r, 
                                p1->pix[yk][xk] );
                        }
                        /* Mark wrong value as negative */
                        p1->pix[yk][xk] = -abs(p1->pix[yk][xk]);
                    }
                }
                if (nk > 0 && nj <= sieve_j) printf ("\n");
            }
             break;
        }

        case NPIC_IMAGE_3L : {
            Npic_image_3l *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int z, y, x, zd, yd, xd,
                zb, yb, xb, zk, yk, xk;

            /* Orthant by orthant */
            for (zd = 1; zd >= -1; zd -=2)
            for (yd = 1; yd >= -1; yd -=2)
            for (xd = 1; xd >= -1; xd -=2)
            for (z = 0; z <= sieve_r; z++)
            for (y = 0; y <= sieve_r; y++)
            for (x = 0; x <= sieve_r; x++)
            {
                zb = sieve_r + z*zd;
                yb = sieve_r + y*yd;
                xb = sieve_r + x*xd;

                /* Visible point ? */
                if (p2->pix[zb][yb][xb] != 1) continue;

                /* Periods */
                nk = 0;
                for (k = 2; 
                            z*k <= sieve_r &&
                            y*k <= sieve_r && 
                            x*k <= sieve_r ; k++) 
                {
                    zk = sieve_r + z*zd*k;
                    yk = sieve_r + y*yd*k;
                    xk = sieve_r + x*xd*k;
                    p2->pix[zk][yk][xk] = 0;  /* mark as non visible */

                    /* Test positive definiteness and homogeneity */
                    if (p1->pix[zb][yb][xb] <= 0 || p1->pix[zb][yb][xb] * k != p1->pix[zk][yk][xk])
                    {
                        nk++; if (nk == 1) nj++;
                        if (nj <= sieve_j) { 
                            if (nk == 1) printf ("  [%d,%d,%d] = %" NPIC_PL, 
                                zb-sieve_r, yb-sieve_r, xb-sieve_r, 
                                p1->pix[zb][yb][xb] );
                            printf ("  [%d,%d,%d] = %" NPIC_PL,
                                zk-sieve_r, yk-sieve_r, xk-sieve_r, 
                                p1->pix[zk][yk][xk] );
                        }
                        /* Mark wrong value as negative */
                        p1->pix[zk][yk][xk] = -abs(p1->pix[zk][yk][xk]);
                    }
                }
                if (nk > 0 && nj <= sieve_j) printf ("\n");
            }
             break;
        }

        case NPIC_IMAGE_4L : {
            Npic_image_4l *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int t, z, y, x, td, zd, yd, xd,
                tb, zb, yb, xb, tk, zk, yk, xk;

            /* Orthant by orthant */
            for (td = 1; td >= -1; td -=2)
            for (zd = 1; zd >= -1; zd -=2)
            for (yd = 1; yd >= -1; yd -=2)
            for (xd = 1; xd >= -1; xd -=2)
            for (t = 0; t <= sieve_r; t++)
            for (z = 0; z <= sieve_r; z++)
            for (y = 0; y <= sieve_r; y++)
            for (x = 0; x <= sieve_r; x++)
            {
                tb = sieve_r + t*td;
                zb = sieve_r + z*zd;
                yb = sieve_r + y*yd;
                xb = sieve_r + x*xd;

                /* Visible point ? */
                if (p2->pix[tb][zb][yb][xb] != 1) continue;

                /* Periods */
                nk = 0;
                for (k = 2; 
                            t*k <= sieve_r &&
                            z*k <= sieve_r &&
                            y*k <= sieve_r && 
                            x*k <= sieve_r ; k++) 
                {
                    tk = sieve_r + t*td*k;
                    zk = sieve_r + z*zd*k;
                    yk = sieve_r + y*yd*k;
                    xk = sieve_r + x*xd*k;
                    p2->pix[tk][zk][yk][xk] = 0;  /* mark as non visible */

                    /* Test positive definiteness and homogeneity */
                    if (p1->pix[tb][zb][yb][xb] <= 0 || p1->pix[tb][zb][yb][xb] * k != p1->pix[tk][zk][yk][xk])
                    {
                        nk++; if (nk == 1) nj++;
                        if (nj <= sieve_j) { 
                            if (nk == 1) printf ("  [%d,%d,%d,%d] = %" NPIC_PL, 
                                tb-sieve_r, zb-sieve_r, yb-sieve_r, xb-sieve_r, 
                                p1->pix[tb][zb][yb][xb] );
                            printf ("  [%d,%d,%d,%d] = %" NPIC_PL,
                                tk-sieve_r, zk-sieve_r, yk-sieve_r, xk-sieve_r, 
                                p1->pix[tk][zk][yk][xk] );
                        }
                        /* Mark wrong value as negative */
                        p1->pix[tk][zk][yk][xk] = -abs(p1->pix[tk][zk][yk][xk]);
                    }
                }
                if (nk > 0 && nj <= sieve_j) printf ("\n");
            }
             break;
        }

        case NPIC_IMAGE_5L : {
            Npic_image_5l *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int s, t, z, y, x, sd, td, zd, yd, xd,
                sb, tb, zb, yb, xb, sk, tk, zk, yk, xk;

            /* Orthant by orthant */
            for (sd = 1; sd >= -1; sd -=2)
            for (td = 1; td >= -1; td -=2)
            for (zd = 1; zd >= -1; zd -=2)
            for (yd = 1; yd >= -1; yd -=2)
            for (xd = 1; xd >= -1; xd -=2)
            for (s = 0; s <= sieve_r; s++)
            for (t = 0; t <= sieve_r; t++)
            for (z = 0; z <= sieve_r; z++)
            for (y = 0; y <= sieve_r; y++)
            for (x = 0; x <= sieve_r; x++)
            {
                sb = sieve_r + s*sd;
                tb = sieve_r + t*td;
                zb = sieve_r + z*zd;
                yb = sieve_r + y*yd;
                xb = sieve_r + x*xd;

                /* Visible point ? */
                if (p2->pix[sb][tb][zb][yb][xb] != 1) continue;

                /* Periods */
                nk = 0;
                for (k = 2; 
                            s*k <= sieve_r &&
                            t*k <= sieve_r &&
                            z*k <= sieve_r &&
                            y*k <= sieve_r && 
                            x*k <= sieve_r ; k++) 
                {
                    sk = sieve_r + s*sd*k;
                    tk = sieve_r + t*td*k;
                    zk = sieve_r + z*zd*k;
                    yk = sieve_r + y*yd*k;
                    xk = sieve_r + x*xd*k;
                    p2->pix[sk][tk][zk][yk][xk] = 0;  /* mark as non visible */

                    /* Test positive definiteness and homogeneity */
                    if (p1->pix[sb][tb][zb][yb][xb] <= 0 || p1->pix[sb][tb][zb][yb][xb] * k != p1->pix[sk][tk][zk][yk][xk])
                    {
                        nk++; if (nk == 1) nj++;
                        if (nj <= sieve_j) { 
                            if (nk == 1) printf ("  [%d,%d,%d,%d,%d] = %" NPIC_PL, 
                                sb-sieve_r, tb-sieve_r, zb-sieve_r, yb-sieve_r, xb-sieve_r, 
                                p1->pix[sb][tb][zb][yb][xb] );
                            printf ("  [%d,%d,%d,%d,%d] = %" NPIC_PL,
                                sk-sieve_r, tk-sieve_r, zk-sieve_r, yk-sieve_r, xk-sieve_r, 
                                p1->pix[sk][tk][zk][yk][xk] );
                        }
                        /* Mark wrong value as negative */
                        p1->pix[sk][tk][zk][yk][xk] = -abs(p1->pix[sk][tk][zk][yk][xk]);
                    }
                }
                if (nk > 0 && nj <= sieve_j) printf ("\n");
            }
             break;
        }

        case NPIC_IMAGE_6L : {
            Npic_image_6l *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int r, s, t, z, y, x, rd, sd, td, zd, yd, xd,
                rb, sb, tb, zb, yb, xb, rk, sk, tk, zk, yk, xk;

            /* Orthant by orthant */
            for (rd = 1; rd >= -1; rd -=2)
            for (sd = 1; sd >= -1; sd -=2)
            for (td = 1; td >= -1; td -=2)
            for (zd = 1; zd >= -1; zd -=2)
            for (yd = 1; yd >= -1; yd -=2)
            for (xd = 1; xd >= -1; xd -=2)
            for (r = 0; r <= sieve_r; r++)
            for (s = 0; s <= sieve_r; s++)
            for (t = 0; t <= sieve_r; t++)
            for (z = 0; z <= sieve_r; z++)
            for (y = 0; y <= sieve_r; y++)
            for (x = 0; x <= sieve_r; x++)
            {
                rb = sieve_r + r*rd;
                sb = sieve_r + s*sd;
                tb = sieve_r + t*td;
                zb = sieve_r + z*zd;
                yb = sieve_r + y*yd;
                xb = sieve_r + x*xd;

                /* Visible point ? */
                if (p2->pix[rb][sb][tb][zb][yb][xb] != 1) continue;

                /* Periods */
                nk = 0;
                for (k = 2; r*k <= sieve_r &&
                            s*k <= sieve_r &&
                            t*k <= sieve_r &&
                            z*k <= sieve_r &&
                            y*k <= sieve_r && 
                            x*k <= sieve_r ; k++) 
                {
                    rk = sieve_r + r*rd*k;
                    sk = sieve_r + s*sd*k;
                    tk = sieve_r + t*td*k;
                    zk = sieve_r + z*zd*k;
                    yk = sieve_r + y*yd*k;
                    xk = sieve_r + x*xd*k;
                    p2->pix[rk][sk][tk][zk][yk][xk] = 0;  /* mark as non visible */

                    /* Test positive definiteness and homogeneity */
                    if (p1->pix[rb][sb][tb][zb][yb][xb] <= 0 || p1->pix[rb][sb][tb][zb][yb][xb] * k != p1->pix[rk][sk][tk][zk][yk][xk])
                    {
                        nk++; if (nk == 1) nj++;
                        if (nj <= sieve_j) { 
                            if (nk == 1) printf ("  [%d,%d,%d,%d,%d,%d] = %" NPIC_PL, 
                                rb-sieve_r, sb-sieve_r, tb-sieve_r, zb-sieve_r, yb-sieve_r, xb-sieve_r, 
                                p1->pix[rb][sb][tb][zb][yb][xb] );
                            printf ("  [%d,%d,%d,%d,%d,%d] = %" NPIC_PL,
                                rk-sieve_r, sk-sieve_r, tk-sieve_r, zk-sieve_r, yk-sieve_r, xk-sieve_r, 
                                p1->pix[rk][sk][tk][zk][yk][xk] );
                        }
                        /* Mark wrong value as negative */
                        p1->pix[rk][sk][tk][zk][yk][xk] = -abs(p1->pix[rk][sk][tk][zk][yk][xk]);
                    }
                }
                if (nk > 0 && nj <= sieve_j) printf ("\n");
            }
             break;
        }

        case NPIC_IMAGE_2D : {
            Npic_image_2d *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int y, x, yd, xd,
                yb, xb, yk, xk;

            /* Orthant by orthant */
            for (yd = 1; yd >= -1; yd -=2)
            for (xd = 1; xd >= -1; xd -=2)
            for (y = 0; y <= sieve_r; y++)
            for (x = 0; x <= sieve_r; x++)
            {
                yb = sieve_r + y*yd;
                xb = sieve_r + x*xd;

                /* Visible point ? */
                if (p2->pix[yb][xb] != 1) continue;

                /* Periods */
                nk = 0;
                for (k = 2; 
                            y*k <= sieve_r && 
                            x*k <= sieve_r ; k++) 
                {
                    yk = sieve_r + y*yd*k;
                    xk = sieve_r + x*xd*k;
                    p2->pix[yk][xk] = 0;  /* mark as non visible */

                    /* Test positive definiteness and homogeneity */
                    if (fabs (p1->pix[yb][xb]) <= 0 || fabs (p1->pix[yb][xb] * k - p1->pix[yk][xk]) > D_THRESH)
                    {
                        nk++; if (nk == 1) nj++;
                        if (nj <= sieve_j) { 
                            if (nk == 1) printf ("  [%d,%d] = %lf", 
                                yb-sieve_r, xb-sieve_r, 
                                p1->pix[yb][xb] );
                            printf ("  [%d,%d] = %lf",
                                yk-sieve_r, xk-sieve_r, 
                                p1->pix[yk][xk] );
                        }
                        /* Mark wrong value as negative */
                        p1->pix[yk][xk] = -fabs(p1->pix[yk][xk]);
                    }
                }
                if (nk > 0 && nj <= sieve_j) printf ("\n");
            }
             break;
        }

        case NPIC_IMAGE_3D : {
            Npic_image_3d *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int z, y, x, zd, yd, xd,
                zb, yb, xb, zk, yk, xk;

            /* Orthant by orthant */
            for (zd = 1; zd >= -1; zd -=2)
            for (yd = 1; yd >= -1; yd -=2)
            for (xd = 1; xd >= -1; xd -=2)
            for (z = 0; z <= sieve_r; z++)
            for (y = 0; y <= sieve_r; y++)
            for (x = 0; x <= sieve_r; x++)
            {
                zb = sieve_r + z*zd;
                yb = sieve_r + y*yd;
                xb = sieve_r + x*xd;

                /* Visible point ? */
                if (p2->pix[zb][yb][xb] != 1) continue;

                /* Periods */
                nk = 0;
                for (k = 2; 
                            z*k <= sieve_r &&
                            y*k <= sieve_r && 
                            x*k <= sieve_r ; k++) 
                {
                    zk = sieve_r + z*zd*k;
                    yk = sieve_r + y*yd*k;
                    xk = sieve_r + x*xd*k;
                    p2->pix[zk][yk][xk] = 0;  /* mark as non visible */

                    /* Test positive definiteness and homogeneity */
                    if (fabs (p1->pix[zb][yb][xb]) <= 0 || fabs (p1->pix[zb][yb][xb] * k - p1->pix[zk][yk][xk]) > D_THRESH)
                    {
                        nk++; if (nk == 1) nj++;
                        if (nj <= sieve_j) { 
                            if (nk == 1) printf ("  [%d,%d,%d] = %lf", 
                                zb-sieve_r, yb-sieve_r, xb-sieve_r, 
                                p1->pix[zb][yb][xb] );
                            printf ("  [%d,%d,%d] = %lf",
                                zk-sieve_r, yk-sieve_r, xk-sieve_r, 
                                p1->pix[zk][yk][xk] );
                        }
                        /* Mark wrong value as negative */
                        p1->pix[zk][yk][xk] = -fabs(p1->pix[zk][yk][xk]);
                    }
                }
                if (nk > 0 && nj <= sieve_j) printf ("\n");
            }
             break;
        }

        case NPIC_IMAGE_4D : {
            Npic_image_4d *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int t, z, y, x, td, zd, yd, xd,
                tb, zb, yb, xb, tk, zk, yk, xk;

            /* Orthant by orthant */
            for (td = 1; td >= -1; td -=2)
            for (zd = 1; zd >= -1; zd -=2)
            for (yd = 1; yd >= -1; yd -=2)
            for (xd = 1; xd >= -1; xd -=2)
            for (t = 0; t <= sieve_r; t++)
            for (z = 0; z <= sieve_r; z++)
            for (y = 0; y <= sieve_r; y++)
            for (x = 0; x <= sieve_r; x++)
            {
                tb = sieve_r + t*td;
                zb = sieve_r + z*zd;
                yb = sieve_r + y*yd;
                xb = sieve_r + x*xd;

                /* Visible point ? */
                if (p2->pix[tb][zb][yb][xb] != 1) continue;

                /* Periods */
                nk = 0;
                for (k = 2; 
                            t*k <= sieve_r &&
                            z*k <= sieve_r &&
                            y*k <= sieve_r && 
                            x*k <= sieve_r ; k++) 
                {
                    tk = sieve_r + t*td*k;
                    zk = sieve_r + z*zd*k;
                    yk = sieve_r + y*yd*k;
                    xk = sieve_r + x*xd*k;
                    p2->pix[tk][zk][yk][xk] = 0;  /* mark as non visible */

                    /* Test positive definiteness and homogeneity */
                    if (fabs (p1->pix[tb][zb][yb][xb]) <= 0 || fabs (p1->pix[tb][zb][yb][xb] * k - p1->pix[tk][zk][yk][xk]) > D_THRESH)
                    {
                        nk++; if (nk == 1) nj++;
                        if (nj <= sieve_j) { 
                            if (nk == 1) printf ("  [%d,%d,%d,%d] = %lf", 
                                tb-sieve_r, zb-sieve_r, yb-sieve_r, xb-sieve_r, 
                                p1->pix[tb][zb][yb][xb] );
                            printf ("  [%d,%d,%d,%d] = %lf",
                                tk-sieve_r, zk-sieve_r, yk-sieve_r, xk-sieve_r, 
                                p1->pix[tk][zk][yk][xk] );
                        }
                        /* Mark wrong value as negative */
                        p1->pix[tk][zk][yk][xk] = -fabs(p1->pix[tk][zk][yk][xk]);
                    }
                }
                if (nk > 0 && nj <= sieve_j) printf ("\n");
            }
             break;
        }

        case NPIC_IMAGE_5D : {
            Npic_image_5d *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int s, t, z, y, x, sd, td, zd, yd, xd,
                sb, tb, zb, yb, xb, sk, tk, zk, yk, xk;

            /* Orthant by orthant */
            for (sd = 1; sd >= -1; sd -=2)
            for (td = 1; td >= -1; td -=2)
            for (zd = 1; zd >= -1; zd -=2)
            for (yd = 1; yd >= -1; yd -=2)
            for (xd = 1; xd >= -1; xd -=2)
            for (s = 0; s <= sieve_r; s++)
            for (t = 0; t <= sieve_r; t++)
            for (z = 0; z <= sieve_r; z++)
            for (y = 0; y <= sieve_r; y++)
            for (x = 0; x <= sieve_r; x++)
            {
                sb = sieve_r + s*sd;
                tb = sieve_r + t*td;
                zb = sieve_r + z*zd;
                yb = sieve_r + y*yd;
                xb = sieve_r + x*xd;

                /* Visible point ? */
                if (p2->pix[sb][tb][zb][yb][xb] != 1) continue;

                /* Periods */
                nk = 0;
                for (k = 2; 
                            s*k <= sieve_r &&
                            t*k <= sieve_r &&
                            z*k <= sieve_r &&
                            y*k <= sieve_r && 
                            x*k <= sieve_r ; k++) 
                {
                    sk = sieve_r + s*sd*k;
                    tk = sieve_r + t*td*k;
                    zk = sieve_r + z*zd*k;
                    yk = sieve_r + y*yd*k;
                    xk = sieve_r + x*xd*k;
                    p2->pix[sk][tk][zk][yk][xk] = 0;  /* mark as non visible */

                    /* Test positive definiteness and homogeneity */
                    if (fabs (p1->pix[sb][tb][zb][yb][xb]) <= 0 || fabs (p1->pix[sb][tb][zb][yb][xb] * k - p1->pix[sk][tk][zk][yk][xk]) > D_THRESH)
                    {
                        nk++; if (nk == 1) nj++;
                        if (nj <= sieve_j) { 
                            if (nk == 1) printf ("  [%d,%d,%d,%d,%d] = %lf", 
                                sb-sieve_r, tb-sieve_r, zb-sieve_r, yb-sieve_r, xb-sieve_r, 
                                p1->pix[sb][tb][zb][yb][xb] );
                            printf ("  [%d,%d,%d,%d,%d] = %lf",
                                sk-sieve_r, tk-sieve_r, zk-sieve_r, yk-sieve_r, xk-sieve_r, 
                                p1->pix[sk][tk][zk][yk][xk] );
                        }
                        /* Mark wrong value as negative */
                        p1->pix[sk][tk][zk][yk][xk] = -fabs(p1->pix[sk][tk][zk][yk][xk]);
                    }
                }
                if (nk > 0 && nj <= sieve_j) printf ("\n");
            }
             break;
        }

        case NPIC_IMAGE_6D : {
            Npic_image_6d *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int r, s, t, z, y, x, rd, sd, td, zd, yd, xd,
                rb, sb, tb, zb, yb, xb, rk, sk, tk, zk, yk, xk;

            /* Orthant by orthant */
            for (rd = 1; rd >= -1; rd -=2)
            for (sd = 1; sd >= -1; sd -=2)
            for (td = 1; td >= -1; td -=2)
            for (zd = 1; zd >= -1; zd -=2)
            for (yd = 1; yd >= -1; yd -=2)
            for (xd = 1; xd >= -1; xd -=2)
            for (r = 0; r <= sieve_r; r++)
            for (s = 0; s <= sieve_r; s++)
            for (t = 0; t <= sieve_r; t++)
            for (z = 0; z <= sieve_r; z++)
            for (y = 0; y <= sieve_r; y++)
            for (x = 0; x <= sieve_r; x++)
            {
                rb = sieve_r + r*rd;
                sb = sieve_r + s*sd;
                tb = sieve_r + t*td;
                zb = sieve_r + z*zd;
                yb = sieve_r + y*yd;
                xb = sieve_r + x*xd;

                /* Visible point ? */
                if (p2->pix[rb][sb][tb][zb][yb][xb] != 1) continue;

                /* Periods */
                nk = 0;
                for (k = 2; r*k <= sieve_r &&
                            s*k <= sieve_r &&
                            t*k <= sieve_r &&
                            z*k <= sieve_r &&
                            y*k <= sieve_r && 
                            x*k <= sieve_r ; k++) 
                {
                    rk = sieve_r + r*rd*k;
                    sk = sieve_r + s*sd*k;
                    tk = sieve_r + t*td*k;
                    zk = sieve_r + z*zd*k;
                    yk = sieve_r + y*yd*k;
                    xk = sieve_r + x*xd*k;
                    p2->pix[rk][sk][tk][zk][yk][xk] = 0;  /* mark as non visible */

                    /* Test positive definiteness and homogeneity */
                    if (fabs (p1->pix[rb][sb][tb][zb][yb][xb]) <= 0 || fabs (p1->pix[rb][sb][tb][zb][yb][xb] * k - p1->pix[rk][sk][tk][zk][yk][xk]) > D_THRESH)
                    {
                        nk++; if (nk == 1) nj++;
                        if (nj <= sieve_j) { 
                            if (nk == 1) printf ("  [%d,%d,%d,%d,%d,%d] = %lf", 
                                rb-sieve_r, sb-sieve_r, tb-sieve_r, zb-sieve_r, yb-sieve_r, xb-sieve_r, 
                                p1->pix[rb][sb][tb][zb][yb][xb] );
                            printf ("  [%d,%d,%d,%d,%d,%d] = %lf",
                                rk-sieve_r, sk-sieve_r, tk-sieve_r, zk-sieve_r, yk-sieve_r, xk-sieve_r, 
                                p1->pix[rk][sk][tk][zk][yk][xk] );
                        }
                        /* Mark wrong value as negative */
                        p1->pix[rk][sk][tk][zk][yk][xk] = -fabs(p1->pix[rk][sk][tk][zk][yk][xk]);
                    }
                }
                if (nk > 0 && nj <= sieve_j) printf ("\n");
            }
             break;
        }

    }  /* switch */

    return nj;
}


int sieve_test_norm (Npic_mask *m1, Npic_image **nq1, Npic_image **nq2, 
        int sieve_r, int sieve_j)
{
    int L = sieve_r*2+1, b;

    b = NpicMaskLargestCoord (m1);

    /* Create images */
    switch (m1->type) {
        case NPIC_MASK_2L : *nq1 = NpicCreateImage_2l (L, L, b, b); break;
        case NPIC_MASK_2D : *nq1 = NpicCreateImage_2d (L, L, b, b); break;
        case NPIC_MASK_3L : *nq1 = NpicCreateImage_3l (L, L, L, b, b, b); break;
        case NPIC_MASK_3D : *nq1 = NpicCreateImage_3d (L, L, L, b, b, b); break;
        case NPIC_MASK_4L : *nq1 = NpicCreateImage_4l (L, L, L, L, b, b, b, b); break;
        case NPIC_MASK_4D : *nq1 = NpicCreateImage_4d (L, L, L, L, b, b, b, b); break;
        case NPIC_MASK_5L : *nq1 = NpicCreateImage_5l (L, L, L, L, L, b, b, b, b, b); break;
        case NPIC_MASK_5D : *nq1 = NpicCreateImage_5d (L, L, L, L, L, b, b, b, b, b); break;
        case NPIC_MASK_6L : *nq1 = NpicCreateImage_6l (L, L, L, L, L, L, b, b, b, b, b, b); break;
        case NPIC_MASK_6D : *nq1 = NpicCreateImage_6d (L, L, L, L, L, L, b, b, b, b, b, b); break;
    }
    *nq2 = NpicDupImage (*nq1);

    if (NpicImageIsOK (*nq1, "npic-wd-norm") != NPIC_SUCCESS ||
        NpicImageIsOK (*nq2, "npic-wd-norm") != NPIC_SUCCESS) return -1;

    /* Draw image : all pixels set to 1 except the center set to 0 */
    sieve_draw_image (m1, *nq1, sieve_r);
    NpicCopyImageI (*nq2, *nq1);

    /* Compute DT_inf */
    NpicWDT_inf (*nq1, m1);

    /* Scan DT_inf orthant by orthant */
    return sieve_scan (*nq1, *nq2, sieve_r, sieve_j);
}


void ArgcExit (int argc, int n)
{
    if (argc < n) {
        fprintf (stderr, "ERROR: %d argument(s) missing, "
                         "type \"npic-wd-norm -h\" to get help.\n", n-argc);
        exit (1);
    }
}

void ArgShift (int *argc, char **argv[], int n)
{
    *argc -= n; *argv += n;
}

void ArgUnknown (char *arg)
{
    fprintf (stderr, "ERROR: unknown argument \"%s\", "
             "type \"npic-wd-norm -h\" to get help.\n", arg);

}

int main (int argc, char *argv[])
{
    char *out1, *mask;
    Npic_image *np1 = NULL, *np2 = NULL;
    Npic_mask *m1 = NULL;
    int sieve_r = 0, sieve_j = 10, nj;

    ArgcExit (argc, 2);
    
    if (strcmp (argv[1], "-h") == 0 || strcmp (argv[1], "-help") == 0 ||
        strcmp (argv[1], "--help") == 0) { ShowUsage (); exit (0); }
    
    ArgcExit (argc, 3);
    mask = argv[1];
    sieve_r = atoi(argv[2]);
    out1 = "";
    ArgShift (&argc, &argv, 2);

    while (argc > 1) {
        if (strcmp (argv[1], "-o") == 0) {
            ArgcExit (argc, 3);
            out1 = argv[2];
            ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-p") == 0) {
            ArgcExit (argc, 3);
            sieve_j = atoi(argv[2]);
            ArgShift (&argc, &argv, 2);
        } else { ArgUnknown (argv[1]); exit (1); }
    } 

    if (sieve_j >= 0) printf ("Loading distance mask \"%s\"\n", mask);
    m1 = NpicReadDistanceMask (mask);
    if (m1 == NULL) exit (1);
    if (sieve_j >= 0) printf ("Number of mask vectors: %d\n", m1->gen.nb);

    nj = sieve_test_norm (m1, &np1, &np2, sieve_r, sieve_j);
    if (sieve_j >= 0 && nj > 0)    
        printf ("Number of series of wrong vectors detected: %d\n", nj);
    printf ("NORM: %s\n", nj == 0 ? "YES" : "NO");

    if (out1 != NULL && out1[0] != 0 && np1 != NULL) {
        if (sieve_j >= 0) printf ("Saving image \"%s\"\n", out1);
        if (NpicWriteImage (np1, out1) != NPIC_SUCCESS) exit(1);
    }

    NpicDestroyImage (np2);
    NpicDestroyImage (np1);
    NpicDestroyMask (m1);
    exit (0);
}