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

/*
 * npic-bissec.c - 14/01/2009
 *
*/

#include <npic.h>

void ShowUsage ()
{
    printf (
        "npic-bissec - Compute bissector of two points for two distances.\n"
        "Usage:\n"
        "  npic-bissec  -h | -help | --help  : print help\n"
        "  npic-bissec out1 options : create image of bissector between points\n"
        "\n"
        "Options:\n"
        "  -2|..|-6        : dimension\n"
        "  -size   coord   : size of image\n"
        "  -p1|-p2 coord   : coordinates of centers\n"
        "  -m1|-m2 mask    : weighted distance mask in " NPIC_KNOWN_MASK_EXT " format,\n"
        "                    or \"sed\" for squared Euclidean distance (default).\n"
        "  -a1|-a2 alpha   : normalisation factor, default is 1.0\n"
        "  -col p r1 r2 f  : color of points, regions and frontier, in 0..255, \n"
        "                    default is 0 99 200 255\n"
        "  -thr h          : threshold for comparing normalized distances, default 1E-5\n"
        "\n"
        "The result out1 will contain the bissector of points p1 and p2, computed for\n"
        "distance m1 from p1 and distance m2 from p2. Distances are normalized.\n"
        "\n"
        "coord: y:x or z:y:x or t:z:y:x s:t:z:y:x or r:s:t:z:y:x\n"
        "out1: image in " NPIC_KNOWN_IMAGE_EXT " format\n"
        "\n"
        "Example 1: using d8 and sed\n"
        "  ./npic-bissec tmp1.npz -2 -size 40:30 -p1 20:2 -p2 28:26 \\\n"
        "      -m1 ../masks/d8_2l.nmask -m2 sed\n"
        "  ./npic-print tmp1.npz\n"
        "\n"
        "Exemple 2: using d<3,4> and d<5,7,11>\n"
        "   ./npic-bissec tmp1.pgm -2 -size 700:500 -p1 300:10 -p2 351:490 \\\n"
        "       -m1 ../masks/d-3-4_2l.nmask -a1 3 -m2 ../masks/d-5-7-11_2l.nmask -a2 5\n"
        "   gimp tmp1.pgm &\n"
        "\n"
        "Exemple 3 in dim 3: using sed and d6\n"
        "   ./npic-bissec tmp1.npz -3 -size 80:80:80 -p1 40:10:40 -p2 37:70:45\\\n"
        "       -m1 sed -m2 ../masks/d6_3l.nmask -col 1 0 0 2 -thr 1\n"
        "   ./npic-geomv tmp1.npz tmp1.geomv -show &\n"
        "\n"
    );
}


typedef struct {
        int c_centers, c_region1, c_region2, c_frontier;
        double frontier_threshold;
} Bissec_datas;


void normalise_wd (Npic_image *np1, Npic_mask *nm1, double alpha)
{
    switch (np1->type) {
        case NPIC_IMAGE_2D : {
            Npic_image_2d *p1 = NpicCastImage (np1);
            int y, x;

            for (y = 0; y < p1->ymax; y++)
            for (x = 0; x < p1->xmax; x++)
                p1->pix[y][x] = p1->pix[y][x] / alpha;
        } break;

        case NPIC_IMAGE_3D : {
            Npic_image_3d *p1 = NpicCastImage (np1);
            int z, y, x;

            for (z = 0; z < p1->zmax; z++)
            for (y = 0; y < p1->ymax; y++)
            for (x = 0; x < p1->xmax; x++)
                p1->pix[z][y][x] = p1->pix[z][y][x] / alpha;
        } break;

        case NPIC_IMAGE_4D : {
            Npic_image_4d *p1 = NpicCastImage (np1);
            int t, z, y, x;

            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++)
                p1->pix[t][z][y][x] = p1->pix[t][z][y][x] / alpha;
        } break;

        case NPIC_IMAGE_5D : {
            Npic_image_5d *p1 = NpicCastImage (np1);
            int s, t, z, y, x;

            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++)
                p1->pix[s][t][z][y][x] = p1->pix[s][t][z][y][x] / alpha;
        } break;

        case NPIC_IMAGE_6D : {
            Npic_image_6d *p1 = NpicCastImage (np1);
            int r, s, t, z, y, x;

            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++)
                p1->pix[r][s][t][z][y][x] = p1->pix[r][s][t][z][y][x] / alpha;
        } break;

    }    
}


void normalise_sed (Npic_image *np1, double alpha)
{
    switch (np1->type) {
        case NPIC_IMAGE_2D : {
            Npic_image_2d *p1 = NpicCastImage (np1);
            int y, x;

            for (y = 0; y < p1->ymax; y++)
            for (x = 0; x < p1->xmax; x++)
                p1->pix[y][x] = sqrt (p1->pix[y][x]) / alpha;
        } break;

        case NPIC_IMAGE_3D : {
            Npic_image_3d *p1 = NpicCastImage (np1);
            int z, y, x;

            for (z = 0; z < p1->zmax; z++)
            for (y = 0; y < p1->ymax; y++)
            for (x = 0; x < p1->xmax; x++)
                p1->pix[z][y][x] = sqrt (p1->pix[z][y][x]) / alpha;
        } break;

        case NPIC_IMAGE_4D : {
            Npic_image_4d *p1 = NpicCastImage (np1);
            int t, z, y, x;

            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++)
                p1->pix[t][z][y][x] = sqrt (p1->pix[t][z][y][x]) / alpha;
        } break;

        case NPIC_IMAGE_5D : {
            Npic_image_5d *p1 = NpicCastImage (np1);
            int s, t, z, y, x;

            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++)
                p1->pix[s][t][z][y][x] = sqrt (p1->pix[s][t][z][y][x]) / alpha;
        } break;

        case NPIC_IMAGE_6D : {
            Npic_image_6d *p1 = NpicCastImage (np1);
            int r, s, t, z, y, x;

            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++)
                p1->pix[r][s][t][z][y][x] = sqrt (p1->pix[r][s][t][z][y][x]) / alpha;
        } break;

    }    
}


void compare_distances (Npic_image *np1, Npic_image *np2, Npic_image *np3,
        Bissec_datas *bis)
{
    switch (np1->type) {
        case NPIC_IMAGE_2D : {
            Npic_image_2d *p1 = NpicCastImage (np1);
            Npic_image_2d *p2 = NpicCastImage (np2);
            Npic_image_2c *p3 = NpicCastImage (np3);
            Npic_d a1, a2;
            int y, x;

            for (y = 0; y < p1->ymax; y++)
            for (x = 0; x < p1->xmax; x++)
            {    
                a1 = p1->pix[y][x];
                a2 = p2->pix[y][x];
                if (a1 == 0 || a2 == 0)
                    p3->pix[y][x] = bis->c_centers;
                else if (fabs (a1-a2) < bis->frontier_threshold)
                    p3->pix[y][x] = bis->c_frontier;
                else if (a1 < a2)
                    p3->pix[y][x] = bis->c_region1;
                else
                    p3->pix[y][x] = bis->c_region2;
            }
        } break;

        case NPIC_IMAGE_3D : {
            Npic_image_3d *p1 = NpicCastImage (np1);
            Npic_image_3d *p2 = NpicCastImage (np2);
            Npic_image_3c *p3 = NpicCastImage (np3);
            Npic_d a1, a2;
            int z, y, x;

            for (z = 0; z < p1->zmax; z++)
            for (y = 0; y < p1->ymax; y++)
            for (x = 0; x < p1->xmax; x++)
            {    
                a1 = p1->pix[z][y][x];
                a2 = p2->pix[z][y][x];
                if (a1 == 0 || a2 == 0)
                    p3->pix[z][y][x] = bis->c_centers;
                else if (fabs (a1-a2) < bis->frontier_threshold)
                    p3->pix[z][y][x] = bis->c_frontier;
                else if (a1 < a2)
                    p3->pix[z][y][x] = bis->c_region1;
                else
                    p3->pix[z][y][x] = bis->c_region2;
            }
        } break;

        case NPIC_IMAGE_4D : {
            Npic_image_4d *p1 = NpicCastImage (np1);
            Npic_image_4d *p2 = NpicCastImage (np2);
            Npic_image_4c *p3 = NpicCastImage (np3);
            Npic_d a1, a2;
            int t, z, y, x;

            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++)
            {    
                a1 = p1->pix[t][z][y][x];
                a2 = p2->pix[t][z][y][x];
                if (a1 == 0 || a2 == 0)
                    p3->pix[t][z][y][x] = bis->c_centers;
                else if (fabs (a1-a2) < bis->frontier_threshold)
                    p3->pix[t][z][y][x] = bis->c_frontier;
                else if (a1 < a2)
                    p3->pix[t][z][y][x] = bis->c_region1;
                else
                    p3->pix[t][z][y][x] = bis->c_region2;
            }
        } break;

        case NPIC_IMAGE_5D : {
            Npic_image_5d *p1 = NpicCastImage (np1);
            Npic_image_5d *p2 = NpicCastImage (np2);
            Npic_image_5c *p3 = NpicCastImage (np3);
            Npic_d a1, a2;
            int s, t, z, y, x;

            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++)
            {    
                a1 = p1->pix[s][t][z][y][x];
                a2 = p2->pix[s][t][z][y][x];
                if (a1 == 0 || a2 == 0)
                    p3->pix[s][t][z][y][x] = bis->c_centers;
                else if (fabs (a1-a2) < bis->frontier_threshold)
                    p3->pix[s][t][z][y][x] = bis->c_frontier;
                else if (a1 < a2)
                    p3->pix[s][t][z][y][x] = bis->c_region1;
                else
                    p3->pix[s][t][z][y][x] = bis->c_region2;
            }
        } break;

        case NPIC_IMAGE_6D : {
            Npic_image_6d *p1 = NpicCastImage (np1);
            Npic_image_6d *p2 = NpicCastImage (np2);
            Npic_image_6c *p3 = NpicCastImage (np3);
            Npic_d a1, a2;
            int r, s, t, z, y, x;

            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++)
            {    
                a1 = p1->pix[r][s][t][z][y][x];
                a2 = p2->pix[r][s][t][z][y][x];
                if (a1 == 0 || a2 == 0)
                    p3->pix[r][s][t][z][y][x] = bis->c_centers;
                else if (fabs (a1-a2) < bis->frontier_threshold)
                    p3->pix[r][s][t][z][y][x] = bis->c_frontier;
                else if (a1 < a2)
                    p3->pix[r][s][t][z][y][x] = bis->c_region1;
                else
                    p3->pix[r][s][t][z][y][x] = bis->c_region2;
            }
        } break;

    }    
}


void comp_dist (Npic_image *np1, Npic_vec *v1, Npic_mask *nm1, double alpha)
{
    double time1;

    /* Set all points to 1, except point v1 to 0 */
    switch (np1->type) {
        case NPIC_IMAGE_2L : {
            Npic_image_2l *p1 = NpicCastImage (np1);
            NpicFillWhole_i (np1, 1);
            p1->pix[v1->y][v1->x] = 0;
        } break;
        case NPIC_IMAGE_3L : {
            Npic_image_3l *p1 = NpicCastImage (np1);
            NpicFillWhole_i (np1, 1);
            p1->pix[v1->z][v1->y][v1->x] = 0;
        } break;
        case NPIC_IMAGE_4L : {
            Npic_image_4l *p1 = NpicCastImage (np1);
            NpicFillWhole_i (np1, 1);
            p1->pix[v1->t][v1->z][v1->y][v1->x] = 0;
        } break;
        case NPIC_IMAGE_5L : {
            Npic_image_5l *p1 = NpicCastImage (np1);
            NpicFillWhole_i (np1, 1);
            p1->pix[v1->s][v1->t][v1->z][v1->y][v1->x] = 0;
        } break;
        case NPIC_IMAGE_6L : {
            Npic_image_6l *p1 = NpicCastImage (np1);
            NpicFillWhole_i (np1, 1);
            p1->pix[v1->r][v1->s][v1->t][v1->z][v1->y][v1->x] = 0;
        } break;
        case NPIC_IMAGE_2D : {
            Npic_image_2d *p1 = NpicCastImage (np1);
            NpicFillWhole_d (np1, 1);
            p1->pix[v1->y][v1->x] = 0;
        } break;
        case NPIC_IMAGE_3D : {
            Npic_image_3d *p1 = NpicCastImage (np1);
            NpicFillWhole_d (np1, 1);
            p1->pix[v1->z][v1->y][v1->x] = 0;
        } break;
        case NPIC_IMAGE_4D : {
            Npic_image_4d *p1 = NpicCastImage (np1);
            NpicFillWhole_d (np1, 1);
            p1->pix[v1->t][v1->z][v1->y][v1->x] = 0;
        } break;
        case NPIC_IMAGE_5D : {
            Npic_image_5d *p1 = NpicCastImage (np1);
            NpicFillWhole_d (np1, 1);
            p1->pix[v1->s][v1->t][v1->z][v1->y][v1->x] = 0;
        } break;
        case NPIC_IMAGE_6D : {
            Npic_image_6d *p1 = NpicCastImage (np1);
            NpicFillWhole_d (np1, 1);
            p1->pix[v1->r][v1->s][v1->t][v1->z][v1->y][v1->x] = 0;
        } break;
    }

    /* Set each point to its distance to point v1 */
    if (nm1 != NULL) {
        printf ("Computing WDT ... "); fflush (stdin);
        time1 = NpicWDT_inf (np1, nm1);
        printf ("%.6lf s\n", time1);
    } else {
        printf ("Computing SEDT ... "); fflush (stdin);
        time1 = NpicSEDT_Hirata_inf (np1);
        printf ("%.6lf s\n", time1);
    }

    /* Normalise */
    NpicConvertImage_d (np1);
    if (nm1 != NULL)
          normalise_wd (np1, nm1, alpha);
    else normalise_sed (np1, alpha);
}


void comp_bissec (Npic_image *np1, Npic_image *np2, Npic_image *np3,
    Npic_mask *nm1, Npic_mask *nm2, double alpha1, double alpha2, 
    Npic_vec *v1, Npic_vec *v2, Bissec_datas *bis)
{
    comp_dist (np1, v1, nm1, alpha1);
    comp_dist (np2, v2, nm2, alpha2);
    compare_distances (np1, np2, np3, bis);
}


void check_point (int dim, Npic_vec *vp, Npic_vec *vs)
{
    if ( vp->x < 0 || vp->x >= vs->x ||
         vp->y < 0 || vp->y >= vs->y ||
         (dim >= 3 && (vp->z < 0 || vp->z >= vs->z)) ||
         (dim >= 4 && (vp->t < 0 || vp->t >= vs->t)) ||
         (dim >= 5 && (vp->s < 0 || vp->s >= vs->s)) ||
         (dim >= 6 && (vp->r < 0 || vp->r >= vs->r)) )
    {
        fprintf (stderr, "ERROR: point coordinates out of image size\n"); exit (1);
    }
}


void scan_vec (int dim, char *s, Npic_vec *v)
{
    int ok = 0;

    switch (dim) {
        case 2: ok = sscanf (s, "%d:%d", 
            &v->y, &v->x) == 2; break;
        case 3: ok = sscanf (s, "%d:%d:%d", 
            &v->z, &v->y, &v->x) == 3; break;
        case 4: ok = sscanf (s, "%d:%d:%d:%d", 
            &v->t, &v->z, &v->y, &v->x) == 4; break;
        case 5: ok = sscanf (s, "%d:%d:%d:%d:%d", 
            &v->s, &v->t, &v->z, &v->y, &v->x) == 5; break;
        case 6: ok = sscanf (s, "%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); }
}


Npic_mask *read_mask (char *name, int dim)
{
    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);
        if (dim != NpicMaskDim (nm->gen.type))
          { fprintf (stderr, "ERROR: bad mask dimension\n"); exit (1); }
    }
    return nm;
}

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

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


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

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

int main (int argc, char *argv[])
{
    char *out1, *size = NULL, *mask1 = "sed", *mask2 = "sed", 
        *point1 = NULL, *point2 = NULL;
    Npic_image *np1 = NULL, *np2 = NULL, *np3 = NULL;
    Npic_mask *nm1 = NULL, *nm2 = NULL;
    double alpha1 = 1, alpha2 = 1;
    int dim = -1;
    Npic_vec v1, v2, vs;
    Bissec_datas bis;

    bis.c_centers = 0;
    bis.c_region1 = 99;
    bis.c_region2 = 200;
    bis.c_frontier = 255;
    bis.frontier_threshold = 1E-5;

    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);
    out1 = argv[1];
    ArgShift (&argc, &argv, 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], "-m1") == 0) {
            ArgcExit (argc, 2+1); mask1 = argv[2]; ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-m2") == 0) {
            ArgcExit (argc, 2+1); mask2 = argv[2]; ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-a1") == 0) {
            ArgcExit (argc, 2+1); alpha1 = atof(argv[2]); ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-a2") == 0) {
            ArgcExit (argc, 2+1); alpha2 = atof(argv[2]); ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-p1") == 0) {
            ArgcExit (argc, 2+1); point1 = argv[2]; ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-p2") == 0) {
            ArgcExit (argc, 2+1); point2 = argv[2]; ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-col") == 0) {
            ArgcExit (argc, 5+1); 
            bis.c_centers  = atoi(argv[2]); 
            bis.c_region1  = atoi(argv[3]); 
            bis.c_region2  = atoi(argv[4]); 
            bis.c_frontier = atoi(argv[5]); 
            ArgShift (&argc, &argv, 5);
        } else if (strcmp (argv[1], "-thr") == 0) {
            ArgcExit (argc, 2+1); 
            bis.frontier_threshold = atof(argv[2]); 
            ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-2") == 0) {
            dim = 2; ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-3") == 0) {
            dim = 3; ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-4") == 0) {
            dim = 4; ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-5") == 0) {
            dim = 5; ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-6") == 0) {
            dim = 6; ArgShift (&argc, &argv, 1);
        } else {
            fprintf (stderr, "ERROR: unknown argument \"%s\", "
                     "type \"npic-bissec -h\" to get help.\n", argv[1]);
            exit (1);
        }
    }

    if (size == NULL || dim == -1) 
      { fprintf (stderr, "ERROR: size or dimension missing for out1\n"); exit (1); }

    scan_vec (dim, size, &vs);

    if (point1 == NULL || point2 == NULL) 
      { fprintf (stderr, "ERROR: point1 or point2 missing\n"); exit (1); }

    scan_vec (dim, point1, &v1);
    scan_vec (dim, point2, &v2);

    check_point (dim, &v1, &vs);
    check_point (dim, &v2, &vs);

    if (alpha1 == 0 || alpha2 == 0) 
      { fprintf (stderr, "ERROR: normalisation factor = 0\n"); exit (1); }

    nm1 = read_mask (mask1, dim);
    nm2 = read_mask (mask2, dim);

    np1 = new_image (dim, nm1, &vs);
    np2 = new_image (dim, nm2, &vs);
    np3 = NpicCreateImageDPS (dim, NPIC_C, &vs);

    if (np1 == NULL || np2 == NULL || np3 == NULL) exit (1);

    comp_bissec (np1, np2, np3, nm1, nm2, alpha1, alpha2, &v1, &v2, &bis);

    printf ("Saving image \"%s\"\n", out1);
    if (NpicWriteImage (np3, out1) != NPIC_SUCCESS) exit(1);
    
    NpicDestroyImage (np1);
    NpicDestroyImage (np2);
    NpicDestroyImage (np3);
    NpicDestroyMask (nm1);
    NpicDestroyMask (nm2);
    exit (0);
}