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

/*
 * npic-metbases.c - 07/10/2001
 *
*/

#include <npic.h>
#include <unistd.h>
#include <signal.h>


void ShowUsage ()
{
    printf (
        "npic-metbases - Compute metric bases and dimension.\n"
        "Usage:\n"
        "  npic-metbases  -h | -help | --help  : print help\n"
        "  npic-metbases options action : check or find bases and dimension\n"
        "\n"
        "Options:\n"
        "  -2|..|-6     : dimension of space\n"
        "  -size coord  : size of the rectangle\n"
        "  -m mask      : weighted distance mask in " NPIC_KNOWN_MASK_EXT " format,\n"
        "                 or \"sed\" for squared Euclidean distance (default).\n"
        "  -v coord     : add one vertex to the set S (empty by default)\n"
        "  -out out1    : save result in out1.\n"
        "  -list        : print additional results\n"
        "  -trace       : show the progression in stderr\n"
        "\n"
        "Action:\n"
        "  -test-rset      : test if the set S is a resolving set\n"
        "                    Labels in out1: -1 undef, 0 vertex, 1 resolved, >1 unresolved\n"
        "  -find-rsets d   : find resolving sets in metric dimension d, starting from S\n"
        "                    Labels in out1: occurrence number\n"
        "  -find-dim d1 d2 : find metric dimension in [d1..d2]\n"
        "                    Labels in out1: 0 background, 1 vertex of the first base\n"
       "\n"
        "out1: image in " NPIC_KNOWN_IMAGE_EXT " format\n"
        "coord: y:x or z:y:x or t:z:y:x or s:t:z:y:x or r:s:t:z:y:x\n"
        "\n"
        "Example 1: test a resolving set\n"
        "   ./npic-metbases -size 10:20 -m ../masks/d-5-7-11_2l.nmask -out tmp2.npz \\\n"
        "       -v 0:0 -v 9:17 -test-rset\n"
        "   ./npic-print tmp2.npz | sed -e 's/ 1/ ./g'\n"
         "\n"
        "Example 2: find bases of dimension 2\n"
        "   ./npic-metbases -size 8:11 -m ../masks/d-3-4_2l.nmask -out tmp2.npz \\\n"
        "       -find-rsets 2\n"
        "   ./npic-print tmp2.npz | sed -e 's/ 0/ ./g'\n"
         "\n"
        "Example 3: find bases given a forcing subset\n"
        "   ./npic-metbases -size 8:11 -m ../masks/d-3-4_2l.nmask -out tmp2.npz \\\n"
        "       -find-rsets 2 -v 2:1 -list\n"
        "   ./npic-print tmp2.npz | sed -e 's/ 0/ ./g'\n"
         "\n"
         "Example 4: searching metric dimension\n"
         "   for f in ../masks/d*_2l.nmask ; do \\\n"
         "       ./npic-metbases -size 10:15 -m \"$f\" -find-dim 1 4 -list ; done \n"
        "\n"
    );
}


#define NVER_MAX 64

typedef struct {
    Npic_vec rsize;                /* Size of the rectangle       */
    int rdim;                      /* Dim of the rectangle        */
    Npic_vec ver[NVER_MAX];        /* Coordinates of vertices     */
    int nver,                      /* Number of vertices in ver[] */
        fver;                      /* Number of forcing vertices  */
} Mba_datas;

typedef struct {
    Npic_vec v;               /* Coordinate in Cartesian space    */
    int n, b[NVER_MAX];       /* Metric dimension and coordinates */
} Multi_coord;


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


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


int bor_signal (int sig, void (*h)(int), int options)
{
    int r;
    struct sigaction s;
    s.sa_handler = h;
    sigemptyset (&s.sa_mask);
    s.sa_flags = options;
    r = sigaction (sig, &s, NULL);
    if (r < 0) perror ("bor_signal");
    return r;
}


int glo_print_trace = 0;

void catch_alarm (int sig)
{
    glo_print_trace = 1;
}


Npic_mask *read_mask (char *name)
{
    Npic_mask *nm = NULL;

    if (strcmp (name, "sed") != 0) {   
        printf ("Loading distance mask \"%s\"\n", name);
        nm = NpicReadDistanceMask (name);
        if (nm == NULL) exit (1);
        printf ("Number of mask vectors: %d\n", nm->gen.nb);
    }
    return nm;
}


int get_rdim (Npic_mask *mask, int rdim)
{
    int d;
    if (mask == NULL) return rdim;

    d = NpicMaskDim (mask->gen.type);
    if (rdim == -1 || rdim == d) return d;

    fprintf (stderr, "ERROR: bad mask dimension\n"); 
    exit (1);
}


void scan_vec (int rdim, char *coords, Npic_vec *v)
{
    int ok = 0;

    switch (rdim) {
        case 2: ok = sscanf (coords, "%d:%d", 
            &v->y, &v->x) == 2; break;
        case 3: ok = sscanf (coords, "%d:%d:%d", 
            &v->z, &v->y, &v->x) == 3; break;
        case 4: ok = sscanf (coords, "%d:%d:%d:%d", 
            &v->t, &v->z, &v->y, &v->x) == 4; break;
        case 5: ok = sscanf (coords, "%d:%d:%d:%d:%d", 
            &v->s, &v->t, &v->z, &v->y, &v->x) == 5; break;
        case 6: ok = sscanf (coords, "%d:%d:%d:%d:%d:%d", 
            &v->r, &v->s, &v->t, &v->z, &v->y, &v->x) == 6; break;
    }
    if (!ok) { fprintf (stderr, "ERROR: bad size or coordinate\n"); exit (1); }
}


void add_vertex (char **coords, char *coord, int *nver)
{
    if (*nver >= NVER_MAX) {
        fprintf (stderr, "ERROR: too many vertices\n"); 
        exit (1);
    }

    /* At this point, dimension is unknown */
    coords[*nver] = coord;
    *nver += 1;
}


Npic_image *new_image (int rdim, Npic_mask *nm, Npic_vec *rsize)
{
    int pixty = NPIC_L;
    if (nm != NULL) pixty = NpicMaskPixelType (nm->type);

    /* create a L or D image */
    return NpicCreateImageDPS (rdim, pixty, rsize);
}


Npic_image *create_dt (Npic_image *np1, Npic_mask *nm1)
{
    Npic_image *np2 = NULL;

    switch (np1->type) {
        case NPIC_IMAGE_2L : {
            Npic_image_2l *p2;
            np2 = NpicCreateImage_2l (np1->gen.ymax*2+1, np1->gen.xmax*2+1, 0, 0);
            NpicFillWhole_i (np2, 1);
            p2 = NpicCastImage (np2);
            p2->pix[p2->ymax/2][p2->xmax/2] = 0;
        } break;

        case NPIC_IMAGE_3L : {
            Npic_image_3l *p2;
            np2 = NpicCreateImage_3l (np1->gen.zmax*2+1, np1->gen.ymax*2+1, np1->gen.xmax*2+1, 0, 0, 0);
            NpicFillWhole_i (np2, 1);
            p2 = NpicCastImage (np2);
            p2->pix[p2->zmax/2][p2->ymax/2][p2->xmax/2] = 0;
        } break;

        case NPIC_IMAGE_4L : {
            Npic_image_4l *p2;
            np2 = NpicCreateImage_4l (np1->gen.tmax*2+1, np1->gen.zmax*2+1, np1->gen.ymax*2+1, np1->gen.xmax*2+1, 0, 0, 0, 0);
            NpicFillWhole_i (np2, 1);
            p2 = NpicCastImage (np2);
            p2->pix[p2->tmax/2][p2->zmax/2][p2->ymax/2][p2->xmax/2] = 0;
        } break;

        case NPIC_IMAGE_5L : {
            Npic_image_5l *p2;
            np2 = NpicCreateImage_5l (np1->gen.smax*2+1, np1->gen.tmax*2+1, np1->gen.zmax*2+1, np1->gen.ymax*2+1, np1->gen.xmax*2+1, 0, 0, 0, 0, 0);
            NpicFillWhole_i (np2, 1);
            p2 = NpicCastImage (np2);
            p2->pix[p2->smax/2][p2->tmax/2][p2->zmax/2][p2->ymax/2][p2->xmax/2] = 0;
        } break;

        case NPIC_IMAGE_6L : {
            Npic_image_6l *p2;
            np2 = NpicCreateImage_6l (np1->gen.rmax*2+1, np1->gen.smax*2+1, np1->gen.tmax*2+1, np1->gen.zmax*2+1, np1->gen.ymax*2+1, np1->gen.xmax*2+1, 0, 0, 0, 0, 0, 0);
            NpicFillWhole_i (np2, 1);
            p2 = NpicCastImage (np2);
            p2->pix[p2->rmax/2][p2->smax/2][p2->tmax/2][p2->zmax/2][p2->ymax/2][p2->xmax/2] = 0;
        } break;

    }

    printf ("Computing DT ...\n");
    if (nm1 == NULL)
         NpicSEDT_Hirata_inf (np2);
    else NpicWDT_inf (np2, nm1);

    return np2;
}


void print_vertices (FILE *f, Mba_datas *ba)
{
    int i;

    for (i = 0; i < ba->nver; i++)
    switch (ba->rdim) {
	case 2 : fprintf (f, "( %d , %d ) ", ba->ver[i].y, ba->ver[i].x); break;
	case 3 : fprintf (f, "( %d , %d , %d ) ", ba->ver[i].z, ba->ver[i].y, ba->ver[i].x); break;
	case 4 : fprintf (f, "( %d , %d , %d , %d ) ", ba->ver[i].t, ba->ver[i].z, ba->ver[i].y, ba->ver[i].x); break;
	case 5 : fprintf (f, "( %d , %d , %d , %d , %d ) ", ba->ver[i].s, ba->ver[i].t, ba->ver[i].z, ba->ver[i].y, ba->ver[i].x); break;
	case 6 : fprintf (f, "( %d , %d , %d , %d , %d , %d ) ", ba->ver[i].r, ba->ver[i].s, ba->ver[i].t, ba->ver[i].z, ba->ver[i].y, ba->ver[i].x); break;
    }
}


/* Label all vertices to the value c in np1 */

void mark_vertices (Mba_datas *ba, Npic_image *np1, int c)
{
    switch (np1->type) {
        case NPIC_IMAGE_2L : {
            Npic_image_2l *p1 = NpicCastImage (np1);
            int i;

            for (i = 0; i < ba->nver; i++)
                if ( ba->ver[i].x >= 0 && ba->ver[i].x < p1->xmax 
                  && ba->ver[i].y >= 0 && ba->ver[i].y < p1->ymax
                   )
                p1->pix[ba->ver[i].y][ba->ver[i].x] = c;
        } break;

        case NPIC_IMAGE_3L : {
            Npic_image_3l *p1 = NpicCastImage (np1);
            int i;

            for (i = 0; i < ba->nver; i++)
                if ( ba->ver[i].x >= 0 && ba->ver[i].x < p1->xmax 
                  && ba->ver[i].y >= 0 && ba->ver[i].y < p1->ymax
                  && ba->ver[i].z >= 0 && ba->ver[i].z < p1->zmax
                   )
                p1->pix[ba->ver[i].z][ba->ver[i].y][ba->ver[i].x] = c;
        } break;

        case NPIC_IMAGE_4L : {
            Npic_image_4l *p1 = NpicCastImage (np1);
            int i;

            for (i = 0; i < ba->nver; i++)
                if ( ba->ver[i].x >= 0 && ba->ver[i].x < p1->xmax 
                  && ba->ver[i].y >= 0 && ba->ver[i].y < p1->ymax
                  && ba->ver[i].z >= 0 && ba->ver[i].z < p1->zmax
                  && ba->ver[i].t >= 0 && ba->ver[i].t < p1->tmax
                   )
                p1->pix[ba->ver[i].t][ba->ver[i].z][ba->ver[i].y][ba->ver[i].x] = c;
        } break;

        case NPIC_IMAGE_5L : {
            Npic_image_5l *p1 = NpicCastImage (np1);
            int i;

            for (i = 0; i < ba->nver; i++)
                if ( ba->ver[i].x >= 0 && ba->ver[i].x < p1->xmax 
                  && ba->ver[i].y >= 0 && ba->ver[i].y < p1->ymax
                  && ba->ver[i].z >= 0 && ba->ver[i].z < p1->zmax
                  && ba->ver[i].t >= 0 && ba->ver[i].t < p1->tmax
                  && ba->ver[i].s >= 0 && ba->ver[i].s < p1->smax
                   )
                p1->pix[ba->ver[i].s][ba->ver[i].t][ba->ver[i].z][ba->ver[i].y][ba->ver[i].x] = c;
        } break;

        case NPIC_IMAGE_6L : {
            Npic_image_6l *p1 = NpicCastImage (np1);
            int i;

            for (i = 0; i < ba->nver; i++)
                if ( ba->ver[i].x >= 0 && ba->ver[i].x < p1->xmax 
                  && ba->ver[i].y >= 0 && ba->ver[i].y < p1->ymax
                  && ba->ver[i].z >= 0 && ba->ver[i].z < p1->zmax
                  && ba->ver[i].t >= 0 && ba->ver[i].t < p1->tmax
                  && ba->ver[i].s >= 0 && ba->ver[i].s < p1->smax
                  && ba->ver[i].r >= 0 && ba->ver[i].r < p1->rmax
                   )
                p1->pix[ba->ver[i].r][ba->ver[i].s][ba->ver[i].t][ba->ver[i].z][ba->ver[i].y][ba->ver[i].x] = c;
        } break;

    }
}


/* Increment the value in np1 of each vertex */

void incr_vertices (Mba_datas *ba, Npic_image *np1)
{
    switch (np1->type) {
        case NPIC_IMAGE_2L : {
            Npic_image_2l *p1 = NpicCastImage (np1);
            int i;

            for (i = 0; i < ba->nver; i++)
                if ( ba->ver[i].x >= 0 && ba->ver[i].x < p1->xmax 
                  && ba->ver[i].y >= 0 && ba->ver[i].y < p1->ymax
                   )
                p1->pix[ba->ver[i].y][ba->ver[i].x] ++;
        } break;

        case NPIC_IMAGE_3L : {
            Npic_image_3l *p1 = NpicCastImage (np1);
            int i;

            for (i = 0; i < ba->nver; i++)
                if ( ba->ver[i].x >= 0 && ba->ver[i].x < p1->xmax 
                  && ba->ver[i].y >= 0 && ba->ver[i].y < p1->ymax
                  && ba->ver[i].z >= 0 && ba->ver[i].z < p1->zmax
                   )
                p1->pix[ba->ver[i].z][ba->ver[i].y][ba->ver[i].x] ++;
        } break;

        case NPIC_IMAGE_4L : {
            Npic_image_4l *p1 = NpicCastImage (np1);
            int i;

            for (i = 0; i < ba->nver; i++)
                if ( ba->ver[i].x >= 0 && ba->ver[i].x < p1->xmax 
                  && ba->ver[i].y >= 0 && ba->ver[i].y < p1->ymax
                  && ba->ver[i].z >= 0 && ba->ver[i].z < p1->zmax
                  && ba->ver[i].t >= 0 && ba->ver[i].t < p1->tmax
                   )
                p1->pix[ba->ver[i].t][ba->ver[i].z][ba->ver[i].y][ba->ver[i].x] ++;
        } break;

        case NPIC_IMAGE_5L : {
            Npic_image_5l *p1 = NpicCastImage (np1);
            int i;

            for (i = 0; i < ba->nver; i++)
                if ( ba->ver[i].x >= 0 && ba->ver[i].x < p1->xmax 
                  && ba->ver[i].y >= 0 && ba->ver[i].y < p1->ymax
                  && ba->ver[i].z >= 0 && ba->ver[i].z < p1->zmax
                  && ba->ver[i].t >= 0 && ba->ver[i].t < p1->tmax
                  && ba->ver[i].s >= 0 && ba->ver[i].s < p1->smax
                   )
                p1->pix[ba->ver[i].s][ba->ver[i].t][ba->ver[i].z][ba->ver[i].y][ba->ver[i].x] ++;
        } break;

        case NPIC_IMAGE_6L : {
            Npic_image_6l *p1 = NpicCastImage (np1);
            int i;

            for (i = 0; i < ba->nver; i++)
                if ( ba->ver[i].x >= 0 && ba->ver[i].x < p1->xmax 
                  && ba->ver[i].y >= 0 && ba->ver[i].y < p1->ymax
                  && ba->ver[i].z >= 0 && ba->ver[i].z < p1->zmax
                  && ba->ver[i].t >= 0 && ba->ver[i].t < p1->tmax
                  && ba->ver[i].s >= 0 && ba->ver[i].s < p1->smax
                  && ba->ver[i].r >= 0 && ba->ver[i].r < p1->rmax
                   )
                p1->pix[ba->ver[i].r][ba->ver[i].s][ba->ver[i].t][ba->ver[i].z][ba->ver[i].y][ba->ver[i].x] ++;
        } break;

    }
}


void fill_mc (Multi_coord *mc, Mba_datas *ba, Npic_image *np1, Npic_image *np2)
{
    switch (np1->type) {
        case NPIC_IMAGE_2L : {
            Npic_image_2l *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int i, k, y, x, uy, ux;
            
            /* Store the multi-coordinates in mc */
            k = 0;
            for (y = 0; y < p1->ymax; y++)
            for (x = 0; x < p1->xmax; x++)
            {
                /* Cartesian coordinates */
                mc[k].v.y = y, mc[k].v.x = x;

                /* Metric dimension */
                mc[k].n = ba->nver;

                /* Compute metric coordinates */
                for (i = 0 ; i < ba->nver; i++) {
                    uy = ba->ver[i].y - y + p2->ymax/2,
                    ux = ba->ver[i].x - x + p2->xmax/2;

                    if ( 0 <= ux && ux < p2->xmax 
                      && 0 <= uy && uy < p2->ymax
                       )
                         mc[k].b[i] = p2->pix[uy][ux];
                    else mc[k].b[i] = -1;
                }

                k++;
            }
        } break;

        case NPIC_IMAGE_3L : {
            Npic_image_3l *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int i, k, z, y, x, uz, uy, ux;
            
            /* Store the multi-coordinates in mc */
            k = 0;
            for (z = 0; z < p1->zmax; z++)
            for (y = 0; y < p1->ymax; y++)
            for (x = 0; x < p1->xmax; x++)
            {
                /* Cartesian coordinates */
                mc[k].v.z = z, mc[k].v.y = y, mc[k].v.x = x;

                /* Metric dimension */
                mc[k].n = ba->nver;

                /* Compute metric coordinates */
                for (i = 0 ; i < ba->nver; i++) {
                    uz = ba->ver[i].z - z + p2->zmax/2,
                    uy = ba->ver[i].y - y + p2->ymax/2,
                    ux = ba->ver[i].x - x + p2->xmax/2;

                    if ( 0 <= ux && ux < p2->xmax 
                      && 0 <= uy && uy < p2->ymax
                      && 0 <= uz && uz < p2->zmax
                       )
                         mc[k].b[i] = p2->pix[uz][uy][ux];
                    else mc[k].b[i] = -1;
                }

                k++;
            }
        } break;

        case NPIC_IMAGE_4L : {
            Npic_image_4l *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int i, k, t, z, y, x, ut, uz, uy, ux;
            
            /* Store the multi-coordinates in mc */
            k = 0;
            for (t = 0; t < p1->tmax; t++)
            for (z = 0; z < p1->zmax; z++)
            for (y = 0; y < p1->ymax; y++)
            for (x = 0; x < p1->xmax; x++)
            {
                /* Cartesian coordinates */
                mc[k].v.t = t, mc[k].v.z = z, mc[k].v.y = y, mc[k].v.x = x;

                /* Metric dimension */
                mc[k].n = ba->nver;

                /* Compute metric coordinates */
                for (i = 0 ; i < ba->nver; i++) {
                    ut = ba->ver[i].t - t + p2->tmax/2,
                    uz = ba->ver[i].z - z + p2->zmax/2,
                    uy = ba->ver[i].y - y + p2->ymax/2,
                    ux = ba->ver[i].x - x + p2->xmax/2;

                    if ( 0 <= ux && ux < p2->xmax 
                      && 0 <= uy && uy < p2->ymax
                      && 0 <= uz && uz < p2->zmax
                      && 0 <= ut && ut < p2->tmax
                       )
                         mc[k].b[i] = p2->pix[ut][uz][uy][ux];
                    else mc[k].b[i] = -1;
                }

                k++;
            }
        } break;

        case NPIC_IMAGE_5L : {
            Npic_image_5l *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int i, k, s, t, z, y, x, us, ut, uz, uy, ux;
            
            /* Store the multi-coordinates in mc */
            k = 0;
            for (s = 0; s < p1->smax; s++)
            for (t = 0; t < p1->tmax; t++)
            for (z = 0; z < p1->zmax; z++)
            for (y = 0; y < p1->ymax; y++)
            for (x = 0; x < p1->xmax; x++)
            {
                /* Cartesian coordinates */
                mc[k].v.s = s, mc[k].v.t = t, mc[k].v.z = z, mc[k].v.y = y, mc[k].v.x = x;

                /* Metric dimension */
                mc[k].n = ba->nver;

                /* Compute metric coordinates */
                for (i = 0 ; i < ba->nver; i++) {
                    us = ba->ver[i].s - s + p2->smax/2,
                    ut = ba->ver[i].t - t + p2->tmax/2,
                    uz = ba->ver[i].z - z + p2->zmax/2,
                    uy = ba->ver[i].y - y + p2->ymax/2,
                    ux = ba->ver[i].x - x + p2->xmax/2;

                    if ( 0 <= ux && ux < p2->xmax 
                      && 0 <= uy && uy < p2->ymax
                      && 0 <= uz && uz < p2->zmax
                      && 0 <= ut && ut < p2->tmax
                      && 0 <= us && us < p2->smax
                       )
                         mc[k].b[i] = p2->pix[us][ut][uz][uy][ux];
                    else mc[k].b[i] = -1;
                }

                k++;
            }
        } break;

        case NPIC_IMAGE_6L : {
            Npic_image_6l *p1 = NpicCastImage (np1),
                          *p2 = NpicCastImage (np2);
            int i, k, r, s, t, z, y, x, ur, us, ut, uz, uy, ux;
            
            /* Store the multi-coordinates in mc */
            k = 0;
            for (r = 0; r < p1->rmax; r++)
            for (s = 0; s < p1->smax; s++)
            for (t = 0; t < p1->tmax; t++)
            for (z = 0; z < p1->zmax; z++)
            for (y = 0; y < p1->ymax; y++)
            for (x = 0; x < p1->xmax; x++)
            {
                /* Cartesian coordinates */
                mc[k].v.r = r, mc[k].v.s = s, mc[k].v.t = t, mc[k].v.z = z, mc[k].v.y = y, mc[k].v.x = x;

                /* Metric dimension */
                mc[k].n = ba->nver;

                /* Compute metric coordinates */
                for (i = 0 ; i < ba->nver; i++) {
                    ur = ba->ver[i].r - r + p2->rmax/2,
                    us = ba->ver[i].s - s + p2->smax/2,
                    ut = ba->ver[i].t - t + p2->tmax/2,
                    uz = ba->ver[i].z - z + p2->zmax/2,
                    uy = ba->ver[i].y - y + p2->ymax/2,
                    ux = ba->ver[i].x - x + p2->xmax/2;

                    if ( 0 <= ux && ux < p2->xmax 
                      && 0 <= uy && uy < p2->ymax
                      && 0 <= uz && uz < p2->zmax
                      && 0 <= ut && ut < p2->tmax
                      && 0 <= us && us < p2->smax
                      && 0 <= ur && ur < p2->rmax
                       )
                         mc[k].b[i] = p2->pix[ur][us][ut][uz][uy][ux];
                    else mc[k].b[i] = -1;
                }

                k++;
            }
        } break;

    }
}


int compare_mc (const void *mc1, const void *mc2)
{
    int i, n = ((Multi_coord*)mc1)->n;
    for (i = 0; i < n; i++)
    {
        if (((Multi_coord*)mc1)->b[i] < ((Multi_coord*)mc2)->b[i]) return -1;
        if (((Multi_coord*)mc1)->b[i] > ((Multi_coord*)mc2)->b[i]) return  1;
    }
    return 0;
}


/* Search the unresolved points, i.e. points sharing same metric representation.
   The input datas mc is already sorted ; this function performs a "uniq" on mc.
   Return the total number of unresolved points.
   If the result is 0, then the set ver[0..nver-1] is a resolving set.
*/

int search_multiples (Multi_coord *mc, Mba_datas *ba, Npic_image *np1, int nb, 
    int mark_multi, int print_multi, int find_all)
{
    int res = 0;

    switch (np1->type) {
        case NPIC_IMAGE_2L : {
            Npic_image_2l *p1 = NpicCastImage (np1);
            int c, i, j, k, l; 

            for (k = 0; k < nb; k = j) {

                /* Check if [k] is defined */
                c = 0;
                for (i = 0; i < 2; i++)
                    if (mc[k].b[i] == -1) { c = -1; break; }

                /* Check for multiple occurences in mc */
                j = k+1;
                if (c != -1) {
                    for ( ; j < nb; j++)
                        if (compare_mc (mc+k, mc+j) != 0) break;
                    c = j-k;
                    if (c > 1) {
                        res += c; 
                        if (! find_all) return 1;
                    }
                }

                /* Mark the points [k..j[ */
                if (mark_multi) {
                    for (l = k; l < j; l++)
                        p1->pix[mc[l].v.y][mc[l].v.x] = c;
                }

                /* Print coords */
                if (print_multi && c > 1) {
                    for (l = k; l < j; l++) {
                        printf (" y = %3d ,  x = %3d ,  c = %3d ,  "
                                "s = (", mc[l].v.y, mc[l].v.x, c);
                        for (i = 0; i < ba->nver; i++) 
                            printf ("%c %4d ", i == 0 ? ' ':',', mc[l].b[i]);
                        printf (" )\n");
                    }
                }
            }
        } break;

        case NPIC_IMAGE_3L : {
            Npic_image_3l *p1 = NpicCastImage (np1);
            int c, i, j, k, l; 

            for (k = 0; k < nb; k = j) {

                /* Check if [k] is defined */
                c = 0;
                for (i = 0; i < 3; i++)
                    if (mc[k].b[i] == -1) { c = -1; break; }

                /* Check for multiple occurences in mc */
                j = k+1;
                if (c != -1) {
                    for ( ; j < nb; j++)
                        if (compare_mc (mc+k, mc+j) != 0) break;
                    c = j-k;
                    if (c > 1) {
                        res += c; 
                        if (! find_all) return 1;
                    }
                }

                /* Mark the points [k..j[ */
                if (mark_multi) {
                    for (l = k; l < j; l++)
                        p1->pix[mc[l].v.z][mc[l].v.y][mc[l].v.x] = c;
                }

                /* Print coords */
                if (print_multi && c > 1) {
                    for (l = k; l < j; l++) {
                        printf (" z = %3d ,  y = %3d ,  x = %3d ,  c = %3d ,  "
                                "s = (", mc[l].v.z, mc[l].v.y, mc[l].v.x, c);
                        for (i = 0; i < ba->nver; i++) 
                            printf ("%c %4d ", i == 0 ? ' ':',', mc[l].b[i]);
                        printf (" )\n");
                    }
                }
            }
        } break;

        case NPIC_IMAGE_4L : {
            Npic_image_4l *p1 = NpicCastImage (np1);
            int c, i, j, k, l; 

            for (k = 0; k < nb; k = j) {

                /* Check if [k] is defined */
                c = 0;
                for (i = 0; i < 4; i++)
                    if (mc[k].b[i] == -1) { c = -1; break; }

                /* Check for multiple occurences in mc */
                j = k+1;
                if (c != -1) {
                    for ( ; j < nb; j++)
                        if (compare_mc (mc+k, mc+j) != 0) break;
                    c = j-k;
                    if (c > 1) {
                        res += c; 
                        if (! find_all) return 1;
                    }
                }

                /* Mark the points [k..j[ */
                if (mark_multi) {
                    for (l = k; l < j; l++)
                        p1->pix[mc[l].v.t][mc[l].v.z][mc[l].v.y][mc[l].v.x] = c;
                }

                /* Print coords */
                if (print_multi && c > 1) {
                    for (l = k; l < j; l++) {
                        printf (" t = %3d ,  z = %3d ,  y = %3d ,  x = %3d ,  c = %3d ,  "
                                "s = (", mc[l].v.t, mc[l].v.z, mc[l].v.y, mc[l].v.x, c);
                        for (i = 0; i < ba->nver; i++) 
                            printf ("%c %4d ", i == 0 ? ' ':',', mc[l].b[i]);
                        printf (" )\n");
                    }
                }
            }
        } break;

        case NPIC_IMAGE_5L : {
            Npic_image_5l *p1 = NpicCastImage (np1);
            int c, i, j, k, l; 

            for (k = 0; k < nb; k = j) {

                /* Check if [k] is defined */
                c = 0;
                for (i = 0; i < 5; i++)
                    if (mc[k].b[i] == -1) { c = -1; break; }

                /* Check for multiple occurences in mc */
                j = k+1;
                if (c != -1) {
                    for ( ; j < nb; j++)
                        if (compare_mc (mc+k, mc+j) != 0) break;
                    c = j-k;
                    if (c > 1) {
                        res += c; 
                        if (! find_all) return 1;
                    }
                }

                /* Mark the points [k..j[ */
                if (mark_multi) {
                    for (l = k; l < j; l++)
                        p1->pix[mc[l].v.s][mc[l].v.t][mc[l].v.z][mc[l].v.y][mc[l].v.x] = c;
                }

                /* Print coords */
                if (print_multi && c > 1) {
                    for (l = k; l < j; l++) {
                        printf (" s = %3d ,  t = %3d ,  z = %3d ,  y = %3d ,  x = %3d ,  c = %3d ,  "
                                "s = (", mc[l].v.s, mc[l].v.t, mc[l].v.z, mc[l].v.y, mc[l].v.x, c);
                        for (i = 0; i < ba->nver; i++) 
                            printf ("%c %4d ", i == 0 ? ' ':',', mc[l].b[i]);
                        printf (" )\n");
                    }
                }
            }
        } break;

        case NPIC_IMAGE_6L : {
            Npic_image_6l *p1 = NpicCastImage (np1);
            int c, i, j, k, l; 

            for (k = 0; k < nb; k = j) {

                /* Check if [k] is defined */
                c = 0;
                for (i = 0; i < 6; i++)
                    if (mc[k].b[i] == -1) { c = -1; break; }

                /* Check for multiple occurences in mc */
                j = k+1;
                if (c != -1) {
                    for ( ; j < nb; j++)
                        if (compare_mc (mc+k, mc+j) != 0) break;
                    c = j-k;
                    if (c > 1) {
                        res += c; 
                        if (! find_all) return 1;
                    }
                }

                /* Mark the points [k..j[ */
                if (mark_multi) {
                    for (l = k; l < j; l++)
                        p1->pix[mc[l].v.r][mc[l].v.s][mc[l].v.t][mc[l].v.z][mc[l].v.y][mc[l].v.x] = c;
                }

                /* Print coords */
                if (print_multi && c > 1) {
                    for (l = k; l < j; l++) {
                        printf (" r = %3d ,  s = %3d ,  t = %3d ,  z = %3d ,  y = %3d ,  x = %3d ,  c = %3d ,  "
                                "s = (", mc[l].v.r, mc[l].v.s, mc[l].v.t, mc[l].v.z, mc[l].v.y, mc[l].v.x, c);
                        for (i = 0; i < ba->nver; i++) 
                            printf ("%c %4d ", i == 0 ? ' ':',', mc[l].b[i]);
                        printf (" )\n");
                    }
                }
            }
        } break;

    }

    return res;
}


/* Search the unresolved points, i.e. points sharing same metric representation.
   Return the total number of unresolved points.
   If the result is 0, then the set ver[0..nver-1] is a resolving set.
*/

int test_resolving_set (Mba_datas *ba, Npic_image *np1, Npic_image *np2,
    int mark_multi, int print_multi)
{
    Multi_coord *mc;
    int nb = np1->gen.pixtot, res = 0;

    /* Create an array of multi-coords */
    mc = malloc (sizeof (Multi_coord) * nb);
    if (mc == NULL) { 
        fprintf (stderr, "ERROR malloc %d bytes\n", nb);
        exit (1); 
    }

    /* Fill in the array mc */
    fill_mc (mc, ba, np1, np2);

    /* Sort the array mc */
    qsort (mc, nb, sizeof (Multi_coord), compare_mc);

    /* Search unresolved points. 
       Labels in np1 : -1 undefined, 1 resolved, >1 multiple */
    res = search_multiples (mc, ba, np1, nb, mark_multi, print_multi, 1);

    /* Mark the vertices to 0 in np1 */
    mark_vertices (ba, np1, 0);

    /* Destroy array */
    free (mc);

    return res;  /* Number or unresolved points */
}


/* Return 1 if v1 == v2, else 0 */

int points_are_equal (int rdim, Npic_vec *v1, Npic_vec *v2)
{
    return (              v1->x == v2->x
          &&              v1->y == v2->y
          && (rdim < 3 || v1->z == v2->z)
          && (rdim < 4 || v1->t == v2->t)
          && (rdim < 5 || v1->s == v2->s)
          && (rdim < 6 || v1->r == v2->r)
           ) ? 1 : 0;
}


/* Set v to 0 */

void point_init (Npic_vec *v)
{
    v->x = 0;
    v->y = 0;
    v->z = 0;
    v->t = 0;
    v->s = 0;
    v->r = 0;
}


/* Set v as the next point after src in the rectangle of size rsize.
   Return 0 success, -1 none */

int point_next (int rdim, Npic_vec *rsize, Npic_vec *src, Npic_vec *v)
{
    *v = *src;  /* Init */

    v->x++;
    if (v->x < rsize->x) return 0; else v->x = 0;

    v->y++;
    if (v->y < rsize->y) return 0; else v->y = 0;

    if (rdim < 3) return -1;
    v->z++;
    if (v->z < rsize->z) return 0; else v->z = 0;

    if (rdim < 4) return -1;
    v->t++;
    if (v->t < rsize->t) return 0; else v->t = 0;

    if (rdim < 5) return -1;
    v->s++;
    if (v->s < rsize->s) return 0; else v->s = 0;

    if (rdim < 6) return -1;
    v->r++;
    if (v->r < rsize->r) return 0; else v->r = 0;

   return -1;    
}


/* Set [k] as the next available point after [j].
   If j == -1, init [k] as the first available point.
   A point is available if it's not in the forcing set [0..fver-1].
   Return 0 if success, else -1.
*/

int candidate_next (Mba_datas *ba, int j, int k)
{
    int i;

    /* First or next point */
    if (j < 0) point_init (&ba->ver[k]);
    else if (point_next (ba->rdim, &ba->rsize, &ba->ver[j], &ba->ver[k]) < 0) return -1;

    /* Check if the point is not in the forcing set [0..fver-1] */
    for (i = 0; i < ba->fver; i++) {
        if (points_are_equal (ba->rdim, &ba->ver[i], &ba->ver[k])) {
            /* Try next point */
            if (point_next (ba->rdim, &ba->rsize, &ba->ver[i], &ba->ver[k]) < 0) return -1;
            /* Verify from beginning */
            i = -1;
        }
    }

    return 0;
}


/* Init the scan of all possible forced sets of dimension nver.
   Return 0 for success, else -1.
*/

int scan_sets_init (Mba_datas *ba)
{
    int i;

    for (i = ba->fver; i < ba->nver; i++)
        if (candidate_next (ba, (i == ba->fver) ? -1 : i-1, i) < 0) return -1;
    return 0;
}


/* Next position in the scan of all possible forced sets of dimension nver.
   Return 0 for success, else -1 (meaning that the scan is finished).
*/

int scan_sets_next (Mba_datas *ba)
{
    int i, j, ok;

    /* Search the first [i] which can be moved */
    for (i = ba->nver-1; i >= ba->fver; i--) {

        if (candidate_next (ba, i, i) < 0) continue;
    
        /* Init the points [i+1..nver-1] */
        ok = 1;
        for (j = i+1; j < ba->nver; j++)
            if (candidate_next (ba, j-1, j) < 0) { ok = 0; break; }
        if (ok) return 0;
    }

    return -1;
}


/* Find all resolving (forced) sets of dimension dim1.
   Return the number of resolving sets, -1 on error.
*/

int find_resolving_sets (Mba_datas *ba, Npic_image *np1, Npic_image *np2, 
    int fver, int dim1, int mark_rsets, int print_list, int find_all, int print_trace)
{
    Multi_coord *mc;
    int nb = np1->gen.pixtot, res = 0;

    if (dim1 <= ba->nver || dim1 >= NVER_MAX)
      { fprintf (stderr, "ERROR, bad metric dimension"); exit (1); }

    ba->fver = fver;
    ba->nver = dim1;

    /* Create an array of multi-coords */
    mc = malloc (sizeof (Multi_coord) * nb);
    if (mc == NULL) { 
        fprintf (stderr, "ERROR malloc %d bytes\n", nb);
        exit (1); 
    }

    /* Init */
    if (scan_sets_init (ba) < 0) return -1;

    /* Scan all possible sets for vertices [fver..nver-1]*/
    do {
        /* Print a trace */
        if (glo_print_trace) {
            glo_print_trace = 0; alarm (1);
            fprintf (stderr, "testing ");
            print_vertices (stderr, ba); 
            fprintf (stderr, "            \r");
        }

        /* Fill in the array mc */
        fill_mc (mc, ba, np1, np2);

        /* Sort the array mc */
        qsort (mc, nb, sizeof (Multi_coord), compare_mc);

        /* Search unresolved points */
        if (search_multiples (mc, ba, np1, nb, 0, 0, 0) == 0) {
            if (mark_rsets) incr_vertices (ba, np1);
            if (print_list) { 
                if (print_trace) fprintf (stderr, "%80c\r", ' ');
                printf ("find: ");
                print_vertices (stdout, ba); 
                printf ("\n");
            }
            res++;
            if (!find_all) break;
         }

    } while (scan_sets_next (ba) == 0);

    /* Destroy array */
    free (mc);
    if (print_trace) fprintf (stderr, "%80c\r", ' ');

    return res;
}


/* Find the metric dimension in the range [dim1..dim2].
 * The dimension is searched starting from the forced set [0..fver-1],
 * thus fver has to be 0 if the actual metric dimension is wanted.
 * Return the metric dimension, -1 on error or dimension not found.
*/

int find_dim (Mba_datas *ba, Npic_image *np1, Npic_image *np2, int dim1, int dim2,
    int mark_rsets, int print_list, int print_trace)
{
    int res, dim, fver = ba->nver;

    if (dim1 < 1 || dim1 <= fver || dim1 >= NVER_MAX)
      { fprintf (stderr, "ERROR, bad metric dimension"); exit (1); }

    for (dim = dim1; dim <= dim2; dim++) {
        printf ("Trying dim %d ...\n", dim);
        res = find_resolving_sets (ba, np1, np2, fver, dim, mark_rsets, 
            print_list, 0, print_trace);
        if (res > 0) return dim;
    }

    return -1;
}


int main (int argc, char *argv[])
{
    char *mask1 = "sed", *size = NULL, *out1 = NULL, *coords[NVER_MAX];
    Npic_mask *nm1 = NULL;
    Mba_datas ba;
    Npic_image *np1, *np2;
    enum { A_NONE, A_TEST_RSET, A_FIND_RSETS, A_FIND_DIM};
    int i, r, action = A_NONE, dim1 = 0, dim2 = 0, print_list = 0,
        print_trace = 0;

    ba.nver = 0;
    ba.rdim = -1;

    ArgcExit (argc, 1+1);
    
    if (strcmp (argv[1], "-h") == 0 || strcmp (argv[1], "-help") == 0 ||
        strcmp (argv[1], "--help") == 0) { ShowUsage (); exit (0); }
    
    ArgcExit (argc, 1+1);

    while (argc > 1) {
        if (strcmp (argv[1], "-size") == 0) {
            ArgcExit (argc, 2+1); 
            size = argv[2]; 
            ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-m") == 0) {
            ArgcExit (argc, 2+1); 
            mask1 = argv[2]; 
            ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-v") == 0) {
            ArgcExit (argc, 2+1); 
            add_vertex (coords, argv[2], &ba.nver); 
            ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-out") == 0) {
            ArgcExit (argc, 2+1); 
            out1 = argv[2]; 
            ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-list") == 0) {
            ArgcExit (argc, 1+1);
            print_list = 1;
            ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-trace") == 0) {
            ArgcExit (argc, 1+1);
            print_trace = 1;
            ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-test-rset") == 0) {
            ArgcExit (argc, 1+1);
            action = A_TEST_RSET;
            ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-find-rsets") == 0) {
            ArgcExit (argc, 2+1); 
            action = A_FIND_RSETS;
            dim1 = atoi (argv[2]); 
            ArgShift (&argc, &argv, 2);
         } else if (strcmp (argv[1], "-find-dim") == 0) {
            ArgcExit (argc, 3+1); 
            action = A_FIND_DIM;
            dim1 = atoi (argv[2]); 
            dim2 = atoi (argv[3]); 
            ArgShift (&argc, &argv, 3);
        } else if (strcmp (argv[1], "-2") == 0) {
            ba.rdim = 2; 
            ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-3") == 0) {
            ba.rdim = 3; 
            ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-4") == 0) {
            ba.rdim = 4; 
            ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-5") == 0) {
            ba.rdim = 5; 
            ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-6") == 0) {
            ba.rdim = 6; 
            ArgShift (&argc, &argv, 1);
        } else {
            fprintf (stderr, "ERROR: unknown argument \"%s\", "
                     "type \"npic-metbases -h\" to get help.\n", argv[1]);
            exit (1);
        }
    }

    if (action == A_NONE) 
      { fprintf (stderr, "ERROR: action missing\n"); exit (1); }

    nm1 = read_mask (mask1);

    ba.rdim = get_rdim (nm1, ba.rdim);
    if (ba.rdim == -1)
      { fprintf (stderr, "ERROR: dimension missing\n"); exit (1); }

    if (size == NULL)
      { fprintf (stderr, "ERROR: size missing\n"); exit (1); }
    scan_vec (ba.rdim, size, &ba.rsize);

    for (i = 0; i < ba.nver; i++)
        scan_vec (ba.rdim, coords[i], &ba.ver[i]);

    np1 = new_image (ba.rdim, nm1, &ba.rsize);
    if (np1 == NULL) exit (1);

    np2 = create_dt (np1, nm1);

    if (print_trace) {
        bor_signal (SIGALRM, catch_alarm, SA_RESTART);
        alarm (1);
    }

    switch (action) {
        case A_TEST_RSET : {
            printf ("Testing resolving set ...\n");
            r = test_resolving_set (&ba, np1, np2, out1 != NULL, print_list);
            printf ("RESOLVING-SET:     %s\n", (r == 0) ? "YES" : "NO");
            printf ("UNRESOLVED-POINTS: %d\n", r);
        } break;
        case A_FIND_RSETS : {
            printf ("Searching resolving sets ...\n");
             r = find_resolving_sets (&ba, np1, np2, ba.nver, dim1, 
                out1 != NULL, print_list, 1, print_trace);
            printf ("RESOLVING-DIM:  %s\n", (r > 0) ? "YES" : "NO");
            printf ("RESOLVING-SETS: %d\n", r);
        } break;
        case A_FIND_DIM : {
            printf ("Searching dimension ...\n");
            r = find_dim (&ba, np1, np2, dim1, dim2, out1 != NULL, print_list, print_trace);
            printf ("DIM-FOUND:  %s\n", (r != -1) ? "YES" : "NO");
            printf ("METRIC-DIM: %d\n", r);
         } break;
        default : ;
     }

    if (print_trace) alarm (0);

    if (out1 != NULL) {
        printf ("Saving \"%s\" ... \n", out1);
        NpicWriteImage (np1, out1);
    }

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