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

/*
 * npic-mlut.c - 22/11/2008
 *
*/

#include <npic.h>


void ShowUsage ()
{
    printf (
        "npic-mlut - Compute Medial Axis, LUT and Mlut for Distance Transforms.\n"
        "Usage:\n"
        "  npic-mlut  -h | -help | --help   : print help\n"
        "  npic-mlut options ...\n"
        "\n"
        "Options:\n"
        "  -wd   d-mask        : use weighted distance mask\n"
        "  -sed                : or use square Euclidean distance\n"
        "  -2l|-3l|-4l|-5l|-6l : specify dimension (if it cannot be guessed)\n"
        "  -mlut t-mask        : medial axis test mask (a.k.a Mlut), to use or to expand\n"
        "  -create             : create or overwrite t-mask\n"
        "  -rad  Rtarget       : the radius to reach (if it cannot be guessed)\n"
        "  -axis in1 out1      : compute medial axis from distance transform in1\n"
        "  -col  n R           : show the Lut column [n][1..R] for vector number n in t-mask\n"
        "  -ctg  out2          : save CTg\n"
        "\n"
        "d-mask, t-mask : masks in " NPIC_KNOWN_MASK_EXT " format\n"
        "in1, out1, out2 : images in " NPIC_KNOWN_IMAGE_EXT " format\n"
        "\n"
        "Example: compute and show MA from WDT for d4\n"
        "  ./npic-wdt -dt ../masks/d4_2l.nmask ../images/K1b_2c.pgm.gz tmp1.npz\n"
        "  ./npic-mlut -wd ../masks/d4_2l.nmask -axis tmp1.npz tmp2.npz\n"
        "  ./npic-val tmp2.npz tmp3.npz -thr-gt 0 255\n"
        "  ./npic-conv tmp3.npz tmp3.pgm -c\n"
        "  gimp tmp3.pgm&\n"
        " \n"
        "Example: compute and show MA from SEDT\n"
        "  ./npic-sedt -dt -hirata ../images/klette_3c.pan.gz tmp1.npz\n"
        "  ./npic-mlut -sed -axis tmp1.npz tmp2.npz -mlut ../masks/t-sed_3l.nmask.gz\n"
        "  # or: ./npic-mlut -sed -axis tmp1.npz tmp2.npz\n"
        "  ./npic-geomv tmp2.npz tmp3.geomv -show\n"
        "\n"
        "Example: compute Mlut for SED in dim 4, then resume computation\n"
        "  ./npic-mlut -sed -4l -rad 30 -mlut tmp1.nmask -create\n"
        "  ./npic-mlut -sed -4l -rad 50 -mlut tmp1.nmask\n"
        "\n"
        "Example: show a Lut column\n"
        "  ./npic-mlut -wd ../masks/d-5-7-11_2l.nmask -col 2 65\n"
        "\n"
    );
}

enum { D_NONE, D_SE, D_W };


/* Try to find L such that the ball of radius Rtarget will fit in image L^n. */

int comp_L_WD (Npic_mask *nMg, Npic_l Rtarget)
{
    int i, j, L = 1;

    switch (nMg->type) {

        case NPIC_MASK_2L : {
            Npic_mask_2l *Mg = NpicCastMask (nMg);
            for (i = 1, j = 0; i < Mg->nb; i++)
                if (Mg->v[i].h < Mg->v[j].h) j = i;
            if (Mg->nb > 0) L = abs(Mg->v[j].x) * Rtarget / Mg->v[j].h + 1;
        } break;

        case NPIC_MASK_3L : {
            Npic_mask_3l *Mg = NpicCastMask (nMg);
            for (i = 1, j = 0; i < Mg->nb; i++)
                if (Mg->v[i].h < Mg->v[j].h) j = i;
            if (Mg->nb > 0) L = abs(Mg->v[j].x) * Rtarget / Mg->v[j].h + 1;
        } break;

        case NPIC_MASK_4L : {
            Npic_mask_4l *Mg = NpicCastMask (nMg);
            for (i = 1, j = 0; i < Mg->nb; i++)
                if (Mg->v[i].h < Mg->v[j].h) j = i;
            if (Mg->nb > 0) L = abs(Mg->v[j].x) * Rtarget / Mg->v[j].h + 1;
        } break;

        case NPIC_MASK_5L : {
            Npic_mask_5l *Mg = NpicCastMask (nMg);
            for (i = 1, j = 0; i < Mg->nb; i++)
                if (Mg->v[i].h < Mg->v[j].h) j = i;
            if (Mg->nb > 0) L = abs(Mg->v[j].x) * Rtarget / Mg->v[j].h + 1;
        } break;

        case NPIC_MASK_6L : {
            Npic_mask_6l *Mg = NpicCastMask (nMg);
            for (i = 1, j = 0; i < Mg->nb; i++)
                if (Mg->v[i].h < Mg->v[j].h) j = i;
            if (Mg->nb > 0) L = abs(Mg->v[j].x) * Rtarget / Mg->v[j].h + 1;
        } break;
    }
    return L;
}


/* Print a Lut column */

void affi_lutcol (Npic_l *lutcol, Npic_l col_R, int n, Npic_mask *nMlut,
    Npic_image *nCTg, int dist)
{
    char *Possible;
    int R;
    Npic_l w = 1;

    /* Mark possible values in Possible[] */
    Possible = malloc ((col_R+1) * sizeof(char));
    if (Possible == NULL) {
        fprintf (stderr, "ERROR: malloc error\n");
        return;
    }
    NpicMlutPossibleValues (nCTg, Possible, col_R);

    printf ("\n# Column : %d\n", n);
    switch (nCTg->type) {
        case NPIC_IMAGE_2L : {
            Npic_image_2l *CTg = NpicCastImage (nCTg);
            Npic_mask_2l *m = NpicCastMask (nMlut);
            w = CTg->pix[m->v[n].y][m->v[n].x];
            printf ("# Vector : ( %d, %d )\n", 
                m->v[n].y, m->v[n].x);
            printf ("# Weight : %"NPIC_PL"\n", w);
            printf ("# Appear : %"NPIC_PL"\n", m->v[n].h);
        } break;
        case NPIC_IMAGE_3L : {
            Npic_image_3l *CTg = NpicCastImage (nCTg);
            Npic_mask_3l *m = NpicCastMask (nMlut);
            w = CTg->pix[m->v[n].z][m->v[n].y][m->v[n].x];
            printf ("# Vector : ( %d, %d, %d )\n", 
                m->v[n].z, m->v[n].y, m->v[n].x);
            printf ("# Weight : %"NPIC_PL"\n", w);
            printf ("# Appear : %"NPIC_PL"\n", m->v[n].h);
        } break;
        case NPIC_IMAGE_4L : {
            Npic_image_4l *CTg = NpicCastImage (nCTg);
            Npic_mask_4l *m = NpicCastMask (nMlut);
            w = CTg->pix[m->v[n].t][m->v[n].z][m->v[n].y][m->v[n].x];
            printf ("# Vector : ( %d, %d, %d, %d )\n", 
                m->v[n].t, m->v[n].z, m->v[n].y, m->v[n].x);
            printf ("# Weight : %"NPIC_PL"\n", w);
            printf ("# Appear : %"NPIC_PL"\n", m->v[n].h);
        } break;
        case NPIC_IMAGE_5L : {
            Npic_image_5l *CTg = NpicCastImage (nCTg);
            Npic_mask_5l *m = NpicCastMask (nMlut);
            w = CTg->pix[m->v[n].s][m->v[n].t][m->v[n].z][m->v[n].y][m->v[n].x];
            printf ("# Vector : ( %d, %d, %d, %d, %d )\n", 
                m->v[n].s, m->v[n].t, m->v[n].z, m->v[n].y, m->v[n].x);
            printf ("# Weight : %"NPIC_PL"\n", w);
            printf ("# Appear : %"NPIC_PL"\n", m->v[n].h);
        } break;
        case NPIC_IMAGE_6L : {
            Npic_image_6l *CTg = NpicCastImage (nCTg);
            Npic_mask_6l *m = NpicCastMask (nMlut);
            w = CTg->pix[m->v[n].r][m->v[n].s][m->v[n].t][m->v[n].z][m->v[n].y][m->v[n].x];
            printf ("# Vector : ( %d, %d, %d, %d, %d, %d )\n", 
                m->v[n].r, m->v[n].s, m->v[n].t, m->v[n].z, m->v[n].y, m->v[n].x);
            printf ("# Weight : %"NPIC_PL"\n", w);
            printf ("# Appear : %"NPIC_PL"\n", m->v[n].h);
        } break;
    }

    printf ("#    R    col[R]   %s\n#\n", dist == D_W ? "!LMC" : "");
    for (R = 1; R <= col_R; R++) 
    if (Possible[R]) {
        printf ("  %5d   %5"NPIC_PL"  ", R, lutcol[R]);
        if (dist == D_W && lutcol[R] != R+w) printf (" %5"NPIC_PL"  ", lutcol[R]);
        printf ("\n");
    }
    printf ("\n");

    free (Possible);
}


void ArgcExit (int argc, int n)
{
    if (argc < n) {
        fprintf (stderr, "ERROR: %d argument(s) missing, "
                         "type \"npic-mlut -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 *in1 = NULL, *out1 = NULL, *out2 = NULL, *tmask = NULL, *dmask = NULL;
    int col_n = -1, dim = -1, creat = 0, L = -1;
    Npic_l Rtarget = -1, Rverif = 0, Rmax = -1, col_R = -1;
    int dist = D_NONE;
    Npic_mask *nMlut = NULL, *nMh = NULL, *nMg = NULL;
    Npic_image *nCTg = NULL, *np1 = NULL, *np2 = NULL;
    Npic_Lut Lut;

    ArgcExit (argc, 1+1);
    
    if (strcmp (argv[1], "-h") == 0 || strcmp (argv[1], "-help") == 0 ||
        strcmp (argv[1], "--help") == 0) { ShowUsage (); exit (0); }
    
    while (argc > 1) {
        if (strcmp (argv[1], "-wd") == 0) {
            ArgcExit (argc, 2+1);
            dist = D_W; dmask = argv[2];
            ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-sed") == 0) {
            ArgcExit (argc, 1+1);
            dist = D_SE;
            ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-2l") == 0) {
            dim = 2;
            ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-3l") == 0) {
            dim = 3;
            ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-4l") == 0) {
            dim = 4;
            ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-5l") == 0) {
            dim = 5;
            ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-6l") == 0) {
            dim = 6;
            ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-mlut") == 0) {
            ArgcExit (argc, 2 +1);
            tmask = argv[2];
            ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-create") == 0) {
            creat = 1;
            ArgShift (&argc, &argv, 1);
        } else if (strcmp (argv[1], "-rad") == 0) {
            ArgcExit (argc, 2+1);
            Rtarget = atoi (argv[2]);
            ArgShift (&argc, &argv, 2);
        } else if (strcmp (argv[1], "-axis") == 0) {
            ArgcExit (argc, 3+1);
            in1 = argv[2]; out1 = argv[3];
            ArgShift (&argc, &argv, 3);
        } else if (strcmp (argv[1], "-col") == 0) {
            ArgcExit (argc, 3+1);
            col_n = atoi (argv[2]); col_R = atoi (argv[3]);
            ArgShift (&argc, &argv, 3);
        } else if (strcmp (argv[1], "-ctg") == 0) {
            ArgcExit (argc, 2+1);
            out2 = argv[2];
            ArgShift (&argc, &argv, 2);
        } else {
            fprintf (stderr, "ERROR: unknown argument \"%s\", "
                     "type \"npic-mlut -h\" to get help.\n", argv[1]);
            exit (1);
        }
    }

    if (dist == D_NONE) 
      { fprintf (stderr, "ERROR: need  -sed  or  -wd d-mask\n"); exit (1); }
    if (in1 == NULL && Rtarget < 0 && col_n == -1) 
      { fprintf (stderr, "ERROR: need  -axis in1 out1  or  -rad Rtarget  or  -col n R\n"); exit (1); }

    if (in1 != NULL) {
        printf ("Loading DT image in \"%s\"\n", in1);
        np1 = NpicReadImage (in1);
        if (np1 == NULL) exit (1);
        NpicConvertImage_l (np1);
        dim = np1->gen.dim;
    }

    if (dist == D_W) {
        printf ("Loading distance mask \"%s\"\n", dmask);
        nMh = NpicReadDistanceMask (dmask);
        if (nMh == NULL) exit (1);
        printf ("Number of vectors in half mask: %d\n", nMh->gen.nb);
        switch (nMh->type) {
            case NPIC_MASK_2L : dim = 2; break;
            case NPIC_MASK_3L : dim = 3; break;
            case NPIC_MASK_4L : dim = 4; break;
            case NPIC_MASK_5L : dim = 5; break;
            case NPIC_MASK_6L : dim = 6; break;
        }
        nMg = NpicCompGeneMask (nMh);
        if (nMg == NULL) exit (1);
        printf ("Number of vectors in generator: %d\n", nMg->gen.nb);
    }

    if (dim == -1)
      { fprintf (stderr, "ERROR: need  -2l|-3l|-4l|-5l|-6l\n"); exit (1); }

    if (in1 != NULL) {
        double t1 = 0;
        printf ("Computing distance transform ... "); fflush (stdout);
        if (dist == D_W)
             t1 = NpicWDT (np1, nMh);
        else t1 = NpicSEDT_Hirata (np1);
        printf ("%.3lf s\n", t1);

        printf ("Searching max value in DT ... "); fflush (stdout);
        Rmax = NpicMaxValue (np1);
        printf ("%"NPIC_PL"\n", Rmax);
    }

    if (creat || tmask == NULL) {
        char s[200];
        printf ("Creating new mask ...\n");
        nMlut = NpicCreateMaskDP (dim, NPIC_L);
        if (nMlut == NULL) exit (1);
        NpicMaskAddProp (nMlut, "Nature", "MedialAxis");
        NpicMaskAddProp (nMlut, "RVerified", "0");
        if (dist == D_W) strncpy (s, dmask, sizeof(s));
        else sprintf (s, "d-sed_%dl.nmask", dim); 
        NpicMaskAddProp (nMlut, "DistanceMask", s);
    } else {
        printf ("Loading medial axis test mask \"%s\"\n", tmask);
        nMlut = NpicReadMedialAxisMask (tmask);
        if (nMlut == NULL) exit (1);
        printf ("Number of mask vectors: %d\n", nMlut->gen.nb);
        Rverif = atoi (NpicMaskGetPropD (nMlut, "RVerified", "0"));
        printf ("RVerified = %" NPIC_PL "\n", Rverif);
    }

    if (Rmax > Rtarget) Rtarget = Rmax;
    if (col_R > Rtarget) Rtarget = col_R;
    NpicMlutInitLut (&Lut, Rtarget);
    if (dist == D_W) {
        L = comp_L_WD (nMg, Rtarget) + 5;
    } else {
        L = sqrt (Rtarget) + 5;
    }

    printf ("Creating a %d^%d CTg image ...\n", L, dim);
    switch (dim) {
        case 2 : nCTg = NpicCreateImage_2l (L, L, 0, 0); break;
        case 3 : nCTg = NpicCreateImage_3l (L, L, L, 0, 0, 0); break;
        case 4 : nCTg = NpicCreateImage_4l (L, L, L, L, 0, 0, 0, 0); break;
        case 5 : nCTg = NpicCreateImage_5l (L, L, L, L, L, 0, 0, 0, 0, 0); break;
        case 6 : nCTg = NpicCreateImage_6l (L, L, L, L, L, L, 0, 0, 0, 0, 0, 0); break;
    }
    if (nCTg == NULL) exit (1);
    printf ("Computing CTg ...\n");
    if (dist == D_W)
         NpicMlutCTg_WD  (nCTg, nMg);
    else NpicMlutCTg_SED (nCTg);  

    if (out2 != NULL && nCTg != NULL) {
        int k = nCTg->gen.ok; nCTg->gen.ok = NPIC_SUCCESS;
        printf ("Saving %s in \"%s\"\n", 
            (k == NPIC_SUCCESS) ? "CTg" : "wrong CTg", out2);
        if (NpicWriteImage (nCTg, out2) != NPIC_SUCCESS) exit(1);
        nCTg->gen.ok = k;
    }
    if (nCTg->gen.ok != NPIC_SUCCESS) exit (1);

    /* Enlarge Mlut if necessary */
    if (Rtarget > Rverif) {

        printf ("Computing MLut for radii ]%"NPIC_PL" .. %"NPIC_PL"] ...\n", Rverif, Rtarget);
        if (NpicMlutCompLutMask (nCTg, nMlut, &Lut, Rverif, Rtarget, nMg) != NPIC_SUCCESS) 
            exit (1);

        /* Save Mlut if updated */
        if (tmask != NULL) {
            char s[200];
            sprintf (s, "%" NPIC_PL, Rtarget);
            NpicMaskAddProp (nMlut, "RVerified", s);
            printf ("Saving medial axis test mask \"%s\" ...\n", tmask);
            if (NpicWriteMask (nMlut, tmask) != NPIC_SUCCESS) exit (1);
        }
    }

    /* Compute Medial Axis from DT, nMLut, Lut or nCTg */
    if (in1 != NULL) {
        np2 = NpicDupImage (np1);
        if (np2 == NULL) exit (1);

        if (Rtarget > Rverif) {
            printf ("Computing Medial Axis, with Lut already in memory ... \n"); 
            NpicExtractMA_Lut (np1, np2, nMlut, &Lut, Rmax);
        } else {
            printf ("Computing Medial Axis, one Lut col at a time in memory ... \n"); 
            NpicExtractMA_Mlut (np1, np2, nMlut, nCTg, Rmax);
        }
    }

    if (out1 != NULL && np2 != NULL && np2->gen.ok == NPIC_SUCCESS) {
        printf ("Saving medial axis in \"%s\"\n", out1);
        if (NpicWriteImage (np2, out1) != NPIC_SUCCESS) exit(1);
    }

    if (col_n != -1) {
        if (col_R < 0 || col_R > Rtarget) col_R = Rtarget;
        if (col_n < 0 || col_n >= nMlut->gen.nb) {
            fprintf (stderr, "ERROR: bad column number\n");
            exit (1);
        }
        if (Rtarget > Rverif) affi_lutcol (Lut.v[col_n], col_R, col_n, nMlut, nCTg, dist);
        else { 
            Npic_l *lutcol = NpicAllocAndCompLutCol (nMlut, nCTg, col_R, col_n);
            if (lutcol == NULL) exit (1);
            affi_lutcol (lutcol, col_R, col_n, nMlut, nCTg, dist);
            free (lutcol);
        }
    }

    NpicMlutFreeLut (&Lut);
    NpicDestroyImage (nCTg);
    NpicDestroyImage (np1);
    NpicDestroyImage (np2);
    NpicDestroyMask (nMh);
    NpicDestroyMask (nMg);
    NpicDestroyMask (nMlut);

    exit (0);
}