gLib/gmt.c

Go to the documentation of this file.

#include "gmt.h"
#include "jbxl_state.h"


WSGraph  Laplacian(WSGraph vp, int mode)
{
    int  i, j;
    int  da, db, dc, dd, de, df, dg, dh;
    WSGraph lp;

    lp = make_WSGraph(vp.xs, vp.ys, 1);
    if (lp.gp==NULL) return lp;

    if (vp.gp==NULL) {
        lp.state = JBXL_GRAPH_NODATA_ERROR;
        return lp;
    }

    if (mode==4) {
        for (j=1; j<vp.ys-1; j++) {
            for (i=1; i<vp.xs-1; i++) {
                da = Px(vp, i+1, j) + Px(vp, i-1, j);
                db = Px(vp, i, j);
                dc = Px(vp, i, j+1) + Px(vp, i, j-1);
                Px(lp, i, j) = da - 4*db + dc;
            }
        }
    }

    else if (mode==8) {
        for (j=1; j<vp.ys-1; j++) {
            for (i=1; i<vp.xs-1; i++) {
                da = Px(vp, i+1, j) + Px(vp, i-1, j);
                db = Px(vp, i, j+1) + Px(vp, i, j-1);
                dc = Px(vp, i, j);
                dd = Px(vp, i+1, j+1) + Px(vp, i-1, j+1);
                de = Px(vp, i+1, j-1) + Px(vp, i-1, j-1);
                Px(lp, i, j) = da + db - 8*dc + dd + de;
            }
        }
    }

    else {
        for (j=2; j<vp.ys-2; j++) {
            for (i=2; i<vp.xs-2; i++) {
                da = Px(vp, i, j);
                db = Px(vp, i+1, j)   + Px(vp, i-1, j) + Px(vp, i, j+1) + Px(vp, i, j-1);
                dc = Px(vp, i-1, j-2) + Px(vp, i, j-2) + Px(vp, i+1, j-2);
                dd = Px(vp, i-1, j+2) + Px(vp, i, j+2) + Px(vp, i+1, j+2);
                de = Px(vp, i-2, j-1) + Px(vp, i-2, j) + Px(vp, i-2, j+1);
                df = Px(vp, i+2, j-1) + Px(vp, i+2, j) + Px(vp, i+2, j+1);
                dg = Px(vp, i-2, j-2) + Px(vp, i+2, j-2);
                dh = Px(vp, i-2, j+2) + Px(vp, i+2, j+2);
                Px(lp, i, j) = (sWord)((-24*da-8*db+4*(dc+dd+de+df)+2*(dg+dh))/64);
            }
        }
    }

    return lp;
}



WSGraph  xSobel(WSGraph vp) 
{
    int  i, j, k;
    int  da, db, dc, dd, de, nr;
    WSGraph xp;
    
    memset(&xp, 0, sizeof(WSGraph));
    if (vp.gp==NULL) {
        xp.state = JBXL_GRAPH_NODATA_ERROR;
        return xp;
    }
    xp.state = JBXL_NORMAL;

    if (vp.zs<=0) vp.zs = 1; 
    xp = make_WSGraph(vp.xs, vp.ys, vp.zs);
    if (xp.gp==NULL) return xp;

    for (k=0; k<vp.zs; k++) {
        for (j=1; j<vp.ys-1; j++) {
            for (i=1; i<vp.xs-1; i++) {
                da = Vx(vp, i+1, j-1, k) - Vx(vp, i-1, j-1, k);
                db = Vx(vp, i+1, j,   k) - Vx(vp, i-1, j,   k);
                dc = Vx(vp, i+1, j+1, k) - Vx(vp, i-1, j+1, k);
                if (k==0 || k==vp.zs-1) {
                    dd = de = 0;
                    nr = 8;
                }
                else {
                    dd = Vx(vp, i+1, j, k-1) - Vx(vp, i-1, j, k-1);
                    de = Vx(vp, i+1, j, k+1) - Vx(vp, i-1, j, k+1);
                    nr = 12;
                }
                Vx(xp, i, j, k) = (sWord)((da + 2*db + dc + dd + de)/nr);
            }
        }
    }

    return xp;
}



FSGraph  fxSobel(FSGraph vp) 
{
    int    i, j, k;
    double da, db, dc, dd, de, nr;
    FSGraph xp;
    
    memset(&xp, 0, sizeof(FSGraph));
    if (vp.gp==NULL) {
        xp.state = JBXL_GRAPH_NODATA_ERROR;
        return xp;
    }

    if (vp.zs<=0) vp.zs = 1; 
    xp = make_FSGraph(vp.xs, vp.ys, vp.zs);
    if (xp.gp==NULL) return xp;

    for (k=0; k<vp.zs; k++) {
        for (j=1; j<vp.ys-1; j++) {
            for (i=1; i<vp.xs-1; i++) {
                da = Vx(vp, i+1, j-1, k) - Vx(vp, i-1, j-1, k);
                db = Vx(vp, i+1, j,   k) - Vx(vp, i-1, j,   k);
                dc = Vx(vp, i+1, j+1, k) - Vx(vp, i-1, j+1, k);
                if (k==0 || k==vp.zs-1) {
                    dd = de = 0.;
                    nr = 8.;
                }
                else {
                    dd = Vx(vp, i+1, j, k-1) - Vx(vp, i-1, j, k-1);
                    de = Vx(vp, i+1, j, k+1) - Vx(vp, i-1, j, k+1);
                    nr = 12.;
                }
                Vx(xp, i, j, k) = (da + 2.*db + dc + dd + de)/nr;
            }
        }
    }

    return xp;
}



WSGraph  ySobel(WSGraph vp) 
{
    int  i, j, k;
    int  da, db, dc, dd, de, nr;
    WSGraph xp;
    
    memset(&xp, 0, sizeof(WSGraph));
    if (vp.gp==NULL) {
        xp.state = JBXL_GRAPH_NODATA_ERROR;
        return xp;
    }

    if (vp.zs<=0) vp.zs = 1; 
    xp = make_WSGraph(vp.xs, vp.ys, vp.zs);
    if (xp.gp==NULL) return xp;

    for (k=0; k<vp.zs; k++) {
        for (j=1; j<vp.ys-1; j++) {
            for (i=1; i<vp.xs-1; i++) {
                da = Vx(vp, i-1, j+1, k) - Vx(vp, i-1, j-1, k);
                db = Vx(vp, i,   j+1, k) - Vx(vp, i,   j-1, k);
                dc = Vx(vp, i+1, j+1, k) - Vx(vp, i+1, j-1, k);
                if (k==0 || k==vp.zs-1) {
                    dd = de = 0;
                    nr = 8;
                }
                else {
                    dd = Vx(vp, i, j+1, k-1) - Vx(vp, i, j-1, k-1);
                    de = Vx(vp, i, j+1, k+1) - Vx(vp, i, j-1, k+1);
                    nr = 12;
                }
                Vx(xp, i, j, k) = (sWord)((da + 2*db + dc + dd + de)/nr);
            }
        }
    }

    return xp;
}



FSGraph  fySobel(FSGraph vp) 
{
    int    i, j, k;
    double da, db, dc, dd, de, nr;
    FSGraph xp;
    
    memset(&xp, 0, sizeof(FSGraph));
    if (vp.gp==NULL) {
        xp.state = JBXL_GRAPH_NODATA_ERROR;
        return xp;
    }

    if (vp.zs<=0) vp.zs = 1; 
    xp = make_FSGraph(vp.xs, vp.ys, vp.zs);
    if (xp.gp==NULL) return xp;

    for (k=0; k<vp.zs; k++) {
        for (j=1; j<vp.ys-1; j++) {
            for (i=1; i<vp.xs-1; i++) {
                da = Vx(vp, i-1, j+1, k) - Vx(vp, i-1, j-1, k);
                db = Vx(vp, i,   j+1, k) - Vx(vp, i,   j-1, k);
                dc = Vx(vp, i+1, j+1, k) - Vx(vp, i+1, j-1, k);
                if (k==0 || k==vp.zs-1) {
                    dd = de = 0.;
                    nr = 8.;
                }
                else {
                    dd = Vx(vp, i, j+1, k-1) - Vx(vp, i, j-1, k-1);
                    de = Vx(vp, i, j+1, k+1) - Vx(vp, i, j-1, k+1);
                    nr = 12.;
                }
                Vx(xp, i, j, k) = (da + 2.*db + dc + dd + de)/nr;
            }
        }
    }

    return xp;
}



WSGraph  zSobel(WSGraph vp) 
{
    int  i, j, k;
    int  da, db, dc, dd, de;
    WSGraph xp;
    
    memset(&xp, 0, sizeof(WSGraph));
    if (vp.gp==NULL) {
        xp.state = JBXL_GRAPH_NODATA_ERROR;
        return xp;
    }

    if (vp.zs<=1) {
        //fprintf(stderr,"ZSOBEL: no 3D data inputed.\n");
        xp.state = JBXL_GRAPH_IVDARG_ERROR;
        return xp;
    }

    xp = make_WSGraph(vp.xs, vp.ys, vp.zs);
    if (xp.gp==NULL) return xp;

    for (k=1; k<vp.zs-1; k++) {
        for (j=1; j<vp.ys-1; j++) {
            for (i=1; i<vp.xs-1; i++) {
                da = Vx(vp, i-1, j, k+1) - Vx(vp, i-1, j, k-1);
                db = Vx(vp, i+1, j, k+1) - Vx(vp, i+1, j, k-1);
                dc = Vx(vp, i,   j, k+1) - Vx(vp, i,   j, k-1);
                dd = Vx(vp, i, j-1, k+1) - Vx(vp, i, j-1, k-1);
                de = Vx(vp, i, j+1, k+1) - Vx(vp, i, j+1, k-1);
                Vx(xp, i, j, k) = (sWord)((da + db + 2*dc + dd + de)/12);
            }
        }
    }

    for (j=1; j<vp.ys-1; j++) {
        for (i=1; i<vp.xs-1; i++) {
            da = Vx(vp, i-1, j, 1) - Vx(vp, i-1, j, 0);
            db = Vx(vp, i+1, j, 1) - Vx(vp, i+1, j, 0);
            dc = Vx(vp, i,   j, 1) - Vx(vp, i,   j, 0);
            dd = Vx(vp, i, j-1, 1) - Vx(vp, i, j-1, 0);
            de = Vx(vp, i, j+1, 1) - Vx(vp, i, j+1, 0);
            Vx(xp, i, j, 0) = (sWord)((da + db + 2*dc + dd + de)/12);

            da = Vx(vp, i-1, j, vp.zs-1) - Vx(vp, i-1, j, vp.zs-2);
            db = Vx(vp, i+1, j, vp.zs-1) - Vx(vp, i+1, j, vp.zs-2);
            dc = Vx(vp, i,   j, vp.zs-1) - Vx(vp, i,   j, vp.zs-2);
            dd = Vx(vp, i, j-1, vp.zs-1) - Vx(vp, i, j-1, vp.zs-2);
            de = Vx(vp, i, j+1, vp.zs-1) - Vx(vp, i, j+1, vp.zs-2);
            Vx(xp, i, j, vp.zs-1) = (sWord)((da + db + 2*dc + dd + de)/12);
        }
    }

    return xp;
}



FSGraph  fzSobel(FSGraph vp) 
{
    int    i, j, k;
    double da, db, dc, dd, de;
    FSGraph xp;
    
    memset(&xp, 0, sizeof(FSGraph));
    if (vp.gp==NULL) {
        xp.state = JBXL_GRAPH_NODATA_ERROR;
        return xp;
    }

    if (vp.zs<=1) {
        //fprintf(stderr,"FZSOBEL: no 3D data inputed.\n");
        xp.state = JBXL_GRAPH_IVDARG_ERROR;
        return xp;
    }

    xp = make_FSGraph(vp.xs, vp.ys, vp.zs);
    if (xp.gp==NULL) return xp;

    for (k=1; k<vp.zs-1; k++) {
        for (j=1; j<vp.ys-1; j++) {
            for (i=1; i<vp.xs-1; i++) {
                da = Vx(vp, i-1, j, k+1) - Vx(vp, i-1, j, k-1);
                db = Vx(vp, i+1, j, k+1) - Vx(vp, i+1, j, k-1);
                dc = Vx(vp, i,   j, k+1) - Vx(vp, i,   j, k-1);
                dd = Vx(vp, i, j-1, k+1) - Vx(vp, i, j-1, k-1);
                de = Vx(vp, i, j+1, k+1) - Vx(vp, i, j+1, k-1);
                Vx(xp, i, j, k) = (da + db + 2.*dc + dd + de)/12.;
            }
        }
    }

    for (j=1; j<vp.ys-1; j++) {
        for (i=1; i<vp.xs-1; i++) {
            da = Vx(vp, i-1, j, 1) - Vx(vp, i-1, j, 0);
            db = Vx(vp, i+1, j, 1) - Vx(vp, i+1, j, 0);
            dc = Vx(vp, i,   j, 1) - Vx(vp, i,   j, 0);
            dd = Vx(vp, i, j-1, 1) - Vx(vp, i, j-1, 0);
            de = Vx(vp, i, j+1, 1) - Vx(vp, i, j+1, 0);
            Vx(xp, i, j, 0) = (da + db + 2.*dc + dd + de)/12.;

            da = Vx(vp, i-1, j, vp.zs-1) - Vx(vp, i-1, j, vp.zs-2);
            db = Vx(vp, i+1, j, vp.zs-1) - Vx(vp, i+1, j, vp.zs-2);
            dc = Vx(vp, i,   j, vp.zs-1) - Vx(vp, i,   j, vp.zs-2);
            dd = Vx(vp, i, j-1, vp.zs-1) - Vx(vp, i, j-1, vp.zs-2);
            de = Vx(vp, i, j+1, vp.zs-1) - Vx(vp, i, j+1, vp.zs-2);
            Vx(xp, i, j, vp.zs-1) = (da + db + 2.*dc + dd + de)/12.;
        }
    }

    return xp;
}



WSGraph  xxSobel(WSGraph vp)
{
    int  x, y, z, xs, ys, zs, cx, cy, cz, ps;
    int  da, db, dc, dd, de;
    int  df, dg, dh, di, dj, dk, dl, dm;
    WSGraph px;

    memset(&px, 0, sizeof(WSGraph));
    if (vp.gp==NULL) {
        px.state = JBXL_GRAPH_NODATA_ERROR;
        return px;
    }

    xs = vp.xs;
    ys = vp.ys;
    zs = vp.zs;
    ps = xs*ys;

    if (zs<5) { 
        px = make_WSGraph(xs, ys, 1);
        if (px.gp==NULL) return px;

        for (y=2; y<ys-2; y++) {
            cy = y*xs;
            for (x=2; x<xs-2; x++) {
                cx = cy + x;
                da = vp.gp[cx-2*xs+2] - 2*vp.gp[cx-2*xs] + vp.gp[cx-2*xs-2];
                db = vp.gp[cx  -xs+2] - 2*vp.gp[cx  -xs] + vp.gp[cx  -xs-2];
                dc = vp.gp[cx     +2] - 2*vp.gp[cx]      + vp.gp[cx     -2];
                dd = vp.gp[cx  +xs+2] - 2*vp.gp[cx  +xs] + vp.gp[cx  +xs-2];
                de = vp.gp[cx+2*xs+2] - 2*vp.gp[cx+2*xs] + vp.gp[cx+2*xs-2];
                px.gp[cx] = (sWord)((da + 4*db + 6*dc + 4*dd + de)/64);
            }
        }
    }

    else {
        px = make_WSGraph(xs, ys, zs);
        if (px.gp==NULL) return px;

        for (z=2; z<zs-2; z++) {
            cz = z*ps;
            for (y=2; y<ys-2; y++) {
                cy = cz + y*xs;
                for (x=2; x<xs-2; x++) {
                    cx = cy + x;
                    da = vp.gp[cx   +2]    - 2*vp.gp[cx]       + vp.gp[cx-2];    
                    db = vp.gp[cx+xs+2]    - 2*vp.gp[cx+xs]    + vp.gp[cx+xs-2];     
                    dc = vp.gp[cx-xs+2]    - 2*vp.gp[cx-xs]    + vp.gp[cx-xs-2];    
                    dd = vp.gp[cx+ps+2]    - 2*vp.gp[cx+ps]    + vp.gp[cx+ps-2];     
                    de = vp.gp[cx-ps+2]    - 2*vp.gp[cx-ps]    + vp.gp[cx-ps-2];     
                    df = vp.gp[cx+xs+ps+2] - 2*vp.gp[cx+xs+ps] + vp.gp[cx+xs+ps-2];
                    dg = vp.gp[cx+xs-ps+2] - 2*vp.gp[cx+xs-ps] + vp.gp[cx+xs-ps-2];
                    dh = vp.gp[cx-xs+ps+2] - 2*vp.gp[cx-xs+ps] + vp.gp[cx-xs+ps-2];
                    di = vp.gp[cx-xs-ps+2] - 2*vp.gp[cx-xs-ps] + vp.gp[cx-xs-ps-2];
                    dj = vp.gp[cx+2*xs+2]  - 2*vp.gp[cx+2*xs]  + vp.gp[cx+2*xs-2];  
                    dk = vp.gp[cx-2*xs+2]  - 2*vp.gp[cx-2*xs]  + vp.gp[cx-2*xs-2];  
                    dl = vp.gp[cx+2*ps+2]  - 2*vp.gp[cx+2*ps]  + vp.gp[cx+2*ps-2]; 
                    dm = vp.gp[cx-2*ps+2]  - 2*vp.gp[cx-2*ps]  + vp.gp[cx-2*ps-2]; 
                    px.gp[cx] = (sWord)((8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144);
                }
            }
        }
    }

    return px;
}



FSGraph  fxxSobel(FSGraph vp)
{
    int  x, y, z, xs, ys, zs, cx, cy, cz, ps;
    double da, db, dc, dd, de;
    double df, dg, dh, di, dj, dk, dl, dm;
    FSGraph px;

    memset(&px, 0, sizeof(FSGraph));
    if (vp.gp==NULL) {
        px.state = JBXL_GRAPH_NODATA_ERROR;
        return px;
    }

    xs = vp.xs;
    ys = vp.ys;
    zs = vp.zs;
    ps = xs*ys;

    if (zs<5) { 
        px = make_FSGraph(xs, ys, 1);
        if (px.gp==NULL) return px;

        for (y=2; y<ys-2; y++) {
            cy = y*xs;
            for (x=2; x<xs-2; x++) {
                cx = cy + x;
                da = vp.gp[cx-2*xs+2] - 2*vp.gp[cx-2*xs] + vp.gp[cx-2*xs-2];
                db = vp.gp[cx  -xs+2] - 2*vp.gp[cx  -xs] + vp.gp[cx  -xs-2];
                dc = vp.gp[cx     +2] - 2*vp.gp[cx]      + vp.gp[cx     -2];
                dd = vp.gp[cx  +xs+2] - 2*vp.gp[cx  +xs] + vp.gp[cx  +xs-2];
                de = vp.gp[cx+2*xs+2] - 2*vp.gp[cx+2*xs] + vp.gp[cx+2*xs-2];
                px.gp[cx] = (da + 4*db + 6*dc + 4*dd + de)/64.;
            }
        }
    }

    else {
        px = make_FSGraph(xs, ys, zs);
        if (px.gp==NULL) return px;

        for (z=2; z<zs-2; z++) {
            cz = z*ps;
            for (y=2; y<ys-2; y++) {
                cy = cz + y*xs;
                for (x=2; x<xs-2; x++) {
                    cx = cy + x;
                    da = vp.gp[cx   +2]    - 2*vp.gp[cx]       + vp.gp[cx-2];    
                    db = vp.gp[cx+xs+2]    - 2*vp.gp[cx+xs]    + vp.gp[cx+xs-2];     
                    dc = vp.gp[cx-xs+2]    - 2*vp.gp[cx-xs]    + vp.gp[cx-xs-2];    
                    dd = vp.gp[cx+ps+2]    - 2*vp.gp[cx+ps]    + vp.gp[cx+ps-2];     
                    de = vp.gp[cx-ps+2]    - 2*vp.gp[cx-ps]    + vp.gp[cx-ps-2];     
                    df = vp.gp[cx+xs+ps+2] - 2*vp.gp[cx+xs+ps] + vp.gp[cx+xs+ps-2];
                    dg = vp.gp[cx+xs-ps+2] - 2*vp.gp[cx+xs-ps] + vp.gp[cx+xs-ps-2];
                    dh = vp.gp[cx-xs+ps+2] - 2*vp.gp[cx-xs+ps] + vp.gp[cx-xs+ps-2];
                    di = vp.gp[cx-xs-ps+2] - 2*vp.gp[cx-xs-ps] + vp.gp[cx-xs-ps-2];
                    dj = vp.gp[cx+2*xs+2]  - 2*vp.gp[cx+2*xs]  + vp.gp[cx+2*xs-2];  
                    dk = vp.gp[cx-2*xs+2]  - 2*vp.gp[cx-2*xs]  + vp.gp[cx-2*xs-2];  
                    dl = vp.gp[cx+2*ps+2]  - 2*vp.gp[cx+2*ps]  + vp.gp[cx+2*ps-2]; 
                    dm = vp.gp[cx-2*ps+2]  - 2*vp.gp[cx-2*ps]  + vp.gp[cx-2*ps-2]; 
                    px.gp[cx] = (8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144.;
                }
            }
        }
    }

    return px;
}



WSGraph  yySobel(WSGraph vp)
{
    int  x, y, z, xs, ys, zs, cx, cy, cz, ps;
    int  da, db, dc, dd, de;
    int  df, dg, dh, di, dj, dk, dl, dm;
    WSGraph py;
    
    memset(&py, 0, sizeof(WSGraph));
    if (vp.gp==NULL) {
        py.state = JBXL_GRAPH_NODATA_ERROR;
        return py;
    }

    xs = vp.xs;
    ys = vp.ys;
    zs = vp.zs;
    ps = xs*ys;

    if (zs<5) { 
        py = make_WSGraph(xs, ys, 1);
        if (py.gp==NULL) return py;

        for (y=2; y<ys-2; y++) {
            cy = y*xs;
            for (x=2; x<xs-2; x++) {
                cx = cy + x;
                da = vp.gp[cx+2*xs-2] - 2*vp.gp[cx-2] + vp.gp[cx-2*xs-2];
                db = vp.gp[cx+2*xs-1] - 2*vp.gp[cx-1] + vp.gp[cx-2*xs-1]; 
                dc = vp.gp[cx+2*xs]   - 2*vp.gp[cx]   + vp.gp[cx-2*xs]; 
                dd = vp.gp[cx+2*xs+1] - 2*vp.gp[cx+1] + vp.gp[cx-2*xs+1]; 
                de = vp.gp[cx+2*xs+2] - 2*vp.gp[cx+2] + vp.gp[cx-2*xs+2];
                py.gp[cx] = (sWord)((da + 4*db + 6*dc + 4*dd + de)/64);
            }
        }
    }
    else {
        py = make_WSGraph(xs, ys, zs);
        if (py.gp==NULL) return py;

        for (z=2; z<zs-2; z++) {
            cz = z*ps;
            for (y=2; y<ys-2; y++) {
                cy = cz + y*xs;
                for (x=2; x<xs-2; x++) {
                    cx = cy + x;
                    da = vp.gp[cx+2*xs]      - 2*vp.gp[cx]      + vp.gp[cx-2*xs];    
                    db = vp.gp[cx+1+2*xs]    - 2*vp.gp[cx+1]    + vp.gp[cx+1-2*xs]; 
                    dc = vp.gp[cx-1+2*xs]    - 2*vp.gp[cx-1]    + vp.gp[cx-1-2*xs];
                    dd = vp.gp[cx+ps+2*xs]   - 2*vp.gp[cx+ps]   + vp.gp[cx+ps-2*xs]; 
                    de = vp.gp[cx-ps+2*xs]   - 2*vp.gp[cx-ps]   + vp.gp[cx-ps-2*xs]; 
                    df = vp.gp[cx+1+ps+2*xs] - 2*vp.gp[cx+1+ps] + vp.gp[cx+1+ps-2*xs]; 
                    dg = vp.gp[cx+1-ps+2*xs] - 2*vp.gp[cx+1-ps] + vp.gp[cx+1-ps-2*xs]; 
                    dh = vp.gp[cx-1+ps+2*xs] - 2*vp.gp[cx-1+ps] + vp.gp[cx-1+ps-2*xs]; 
                    di = vp.gp[cx-1-ps+2*xs] - 2*vp.gp[cx-1-ps] + vp.gp[cx-1-ps-2*xs]; 
                    dj = vp.gp[cx+2+2*xs]    - 2*vp.gp[cx+2]    + vp.gp[cx+2-2*xs];
                    dk = vp.gp[cx-2+2*xs]    - 2*vp.gp[cx-2]    + vp.gp[cx-2-2*xs];
                    dl = vp.gp[cx+2*ps+2*xs] - 2*vp.gp[cx+2*ps] + vp.gp[cx+2*ps-2*xs];
                    dm = vp.gp[cx-2*ps+2*xs] - 2*vp.gp[cx-2*ps] + vp.gp[cx-2*ps-2*xs];
                    py.gp[cx] = (sWord)((8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144);
                }
            }
        }
    }

    return py;
}



FSGraph  fyySobel(FSGraph vp)
{
    int  x, y, z, xs, ys, zs, cx, cy, cz, ps;
    double da, db, dc, dd, de;
    double df, dg, dh, di, dj, dk, dl, dm;
    FSGraph py;

    memset(&py, 0, sizeof(FSGraph));
    if (vp.gp==NULL) {
        py.state = JBXL_GRAPH_NODATA_ERROR;
        return py;
    }

    xs = vp.xs;
    ys = vp.ys;
    zs = vp.zs;
    ps = xs*ys;

    if (zs<5) { 
        py = make_FSGraph(xs, ys, 1);
        if (py.gp==NULL) return py;

        for (y=2; y<ys-2; y++) {
            cy = y*xs;
            for (x=2; x<xs-2; x++) {
                cx = cy + x;
                da = vp.gp[cx+2*xs-2] - 2*vp.gp[cx-2] + vp.gp[cx-2*xs-2];
                db = vp.gp[cx+2*xs-1] - 2*vp.gp[cx-1] + vp.gp[cx-2*xs-1]; 
                dc = vp.gp[cx+2*xs]   - 2*vp.gp[cx]   + vp.gp[cx-2*xs]; 
                dd = vp.gp[cx+2*xs+1] - 2*vp.gp[cx+1] + vp.gp[cx-2*xs+1]; 
                de = vp.gp[cx+2*xs+2] - 2*vp.gp[cx+2] + vp.gp[cx-2*xs+2];
                py.gp[cx] = (da + 4*db + 6*dc + 4*dd + de)/64.;
            }
        }
    }
    else {
        py = make_FSGraph(xs, ys, zs);
        if (py.gp==NULL) return py;

        for (z=2; z<zs-2; z++) {
            cz = z*ps;
            for (y=2; y<ys-2; y++) {
                cy = cz + y*xs;
                for (x=2; x<xs-2; x++) {
                    cx = cy + x;
                    da = vp.gp[cx+2*xs]      - 2*vp.gp[cx]      + vp.gp[cx-2*xs];    
                    db = vp.gp[cx+1+2*xs]    - 2*vp.gp[cx+1]    + vp.gp[cx+1-2*xs]; 
                    dc = vp.gp[cx-1+2*xs]    - 2*vp.gp[cx-1]    + vp.gp[cx-1-2*xs];
                    dd = vp.gp[cx+ps+2*xs]   - 2*vp.gp[cx+ps]   + vp.gp[cx+ps-2*xs]; 
                    de = vp.gp[cx-ps+2*xs]   - 2*vp.gp[cx-ps]   + vp.gp[cx-ps-2*xs]; 
                    df = vp.gp[cx+1+ps+2*xs] - 2*vp.gp[cx+1+ps] + vp.gp[cx+1+ps-2*xs]; 
                    dg = vp.gp[cx+1-ps+2*xs] - 2*vp.gp[cx+1-ps] + vp.gp[cx+1-ps-2*xs]; 
                    dh = vp.gp[cx-1+ps+2*xs] - 2*vp.gp[cx-1+ps] + vp.gp[cx-1+ps-2*xs]; 
                    di = vp.gp[cx-1-ps+2*xs] - 2*vp.gp[cx-1-ps] + vp.gp[cx-1-ps-2*xs]; 
                    dj = vp.gp[cx+2+2*xs]    - 2*vp.gp[cx+2]    + vp.gp[cx+2-2*xs];
                    dk = vp.gp[cx-2+2*xs]    - 2*vp.gp[cx-2]    + vp.gp[cx-2-2*xs];
                    dl = vp.gp[cx+2*ps+2*xs] - 2*vp.gp[cx+2*ps] + vp.gp[cx+2*ps-2*xs];
                    dm = vp.gp[cx-2*ps+2*xs] - 2*vp.gp[cx-2*ps] + vp.gp[cx-2*ps-2*xs];
                    py.gp[cx] = (8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144.;
                }
            }
        }
    }

    return py;
}



WSGraph  zzSobel(WSGraph vp)
{
    int  x, y, z, xs, ys, zs, cx, cy, cz, ps;
    int  da, db, dc, dd, de;
    int  df, dg, dh, di, dj, dk, dl, dm;
    WSGraph pz;
    
    memset(&pz, 0, sizeof(WSGraph));
    if (vp.gp==NULL) {
        pz.state = JBXL_GRAPH_NODATA_ERROR;
        return pz;
    }

    if (vp.zs<5) {
        //fprintf(stderr,"ZZSOBEL: no 3D data inputed.\n");
        pz.state = JBXL_GRAPH_IVDARG_ERROR;
        return pz;
    }

    xs = vp.xs;
    ys = vp.ys;
    zs = vp.zs;
    ps = xs*ys;
    pz = make_WSGraph(xs, ys, zs);
    if (pz.gp==NULL) return pz;

    for (z=2; z<zs-2; z++) {
        cz = z*ps;
        for (y=2; y<ys-2; y++) {
            cy = cz + y*xs;
            for (x=2; x<xs-2; x++) {
                cx = cy + x;
                da = vp.gp[cx   +2*ps]   - 2*vp.gp[cx]       + vp.gp[cx   -2*ps];    
                db = vp.gp[cx+1 +2*ps]   - 2*vp.gp[cx+1]     + vp.gp[cx+1 -2*ps]; 
                dc = vp.gp[cx-1 +2*ps]   - 2*vp.gp[cx-1]     + vp.gp[cx-1 -2*ps];
                dd = vp.gp[cx+xs+2*ps]   - 2*vp.gp[cx+xs]    + vp.gp[cx+xs-2*ps]; 
                de = vp.gp[cx-xs+2*ps]   - 2*vp.gp[cx-xs]    + vp.gp[cx-xs-2*ps]; 
                df = vp.gp[cx+1+xs+2*ps] - 2*vp.gp[cx+1+xs]  + vp.gp[cx+1+xs-2*ps]; 
                dg = vp.gp[cx+1-xs+2*ps] - 2*vp.gp[cx+1-xs]  + vp.gp[cx+1-xs-2*ps]; 
                dh = vp.gp[cx-1+xs+2*ps] - 2*vp.gp[cx-1+xs]  + vp.gp[cx-1+xs-2*ps]; 
                di = vp.gp[cx-1-xs+2*ps] - 2*vp.gp[cx-1-xs]  + vp.gp[cx-1-xs-2*ps]; 
                dj = vp.gp[cx+2   +2*ps] - 2*vp.gp[cx+2]     + vp.gp[cx+2   -2*ps];
                dk = vp.gp[cx-2   +2*ps] - 2*vp.gp[cx-2]     + vp.gp[cx-2   -2*ps];
                dl = vp.gp[cx+2*xs+2*ps] - 2*vp.gp[cx+2*xs]  + vp.gp[cx+2*xs-2*ps];
                dm = vp.gp[cx-2*xs+2*ps] - 2*vp.gp[cx-2*xs]  + vp.gp[cx-2*xs-2*ps];
                pz.gp[cx] = (sWord)((8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144);
            }
        }
    }

    cz = ps;
    for (y=2; y<ys-2; y++) {
        cy = cz + y*xs;
        for (x=2; x<xs-2; x++) {
            cx = cy + x;
            da = vp.gp[cx   +2*ps]   - 2*vp.gp[cx];  
            db = vp.gp[cx+1 +2*ps]   - 2*vp.gp[cx+1]; 
            dc = vp.gp[cx-1 +2*ps]   - 2*vp.gp[cx-1];
            dd = vp.gp[cx+xs+2*ps]   - 2*vp.gp[cx+xs]; 
            de = vp.gp[cx-xs+2*ps]   - 2*vp.gp[cx-xs]; 
            df = vp.gp[cx+1+xs+2*ps] - 2*vp.gp[cx+1+xs]; 
            dg = vp.gp[cx+1-xs+2*ps] - 2*vp.gp[cx+1-xs]; 
            dh = vp.gp[cx-1+xs+2*ps] - 2*vp.gp[cx-1+xs]; 
            di = vp.gp[cx-1-xs+2*ps] - 2*vp.gp[cx-1-xs]; 
            dj = vp.gp[cx+2   +2*ps] - 2*vp.gp[cx+2];
            dk = vp.gp[cx-2   +2*ps] - 2*vp.gp[cx-2];
            dl = vp.gp[cx+2*xs+2*ps] - 2*vp.gp[cx+2*xs];
            dm = vp.gp[cx-2*xs+2*ps] - 2*vp.gp[cx-2*xs];
            pz.gp[cx] = (sWord)((8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144);
        }
    }

    cz = (zs-2)*ps;
    for (y=2; y<ys-2; y++) {
        cy = cz + y*xs;
        for (x=2; x<xs-2; x++) {
            cx = cy + x;
            da = - 2*vp.gp[cx]      + vp.gp[cx   -2*ps];     
            db = - 2*vp.gp[cx+1]    + vp.gp[cx+1 -2*ps]; 
            dc = - 2*vp.gp[cx-1]    + vp.gp[cx-1 -2*ps];
            dd = - 2*vp.gp[cx+xs]   + vp.gp[cx+xs-2*ps]; 
            de = - 2*vp.gp[cx-xs]   + vp.gp[cx-xs-2*ps]; 
            df = - 2*vp.gp[cx+1+xs] + vp.gp[cx+1+xs-2*ps]; 
            dg = - 2*vp.gp[cx+1-xs] + vp.gp[cx+1-xs-2*ps]; 
            dh = - 2*vp.gp[cx-1+xs] + vp.gp[cx-1+xs-2*ps]; 
            di = - 2*vp.gp[cx-1-xs] + vp.gp[cx-1-xs-2*ps]; 
            dj = - 2*vp.gp[cx+2]    + vp.gp[cx+2   -2*ps];
            dk = - 2*vp.gp[cx-2]    + vp.gp[cx-2   -2*ps];
            dl = - 2*vp.gp[cx+2*xs] + vp.gp[cx+2*xs-2*ps];
            dm = - 2*vp.gp[cx-2*xs] + vp.gp[cx-2*xs-2*ps];
            pz.gp[cx] = (sWord)((8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144);
        }
    }

    return pz;
}



FSGraph  fzzSobel(FSGraph vp)
{
    int  x, y, z, xs, ys, zs, cx, cy, cz, ps;
    double da, db, dc, dd, de;
    double df, dg, dh, di, dj, dk, dl, dm;
    FSGraph pz;
    
    memset(&pz, 0, sizeof(FSGraph));
    if (vp.gp==NULL) {
        pz.state = JBXL_GRAPH_NODATA_ERROR;
        return pz;
    }

    if (vp.zs<5) {
        //fprintf(stderr,"FZZSOBEL: no 3D data inputed.\n");
        pz.state = JBXL_GRAPH_IVDARG_ERROR;
        return pz;
    }

    xs = vp.xs;
    ys = vp.ys;
    zs = vp.zs;
    ps = xs*ys;
    pz = make_FSGraph(xs, ys, zs);
    if (pz.gp==NULL) return pz;

    for (z=2; z<zs-2; z++) {
        cz = z*ps;
        for (y=2; y<ys-2; y++) {
            cy = cz + y*xs;
            for (x=2; x<xs-2; x++) {
                cx = cy + x;
                da = vp.gp[cx     +2*ps] - 2*vp.gp[cx]      + vp.gp[cx     -2*ps];   
                db = vp.gp[cx+1   +2*ps] - 2*vp.gp[cx+1]    + vp.gp[cx+1   -2*ps];
                dc = vp.gp[cx-1   +2*ps] - 2*vp.gp[cx-1]    + vp.gp[cx-1   -2*ps];
                dd = vp.gp[cx  +xs+2*ps] - 2*vp.gp[cx  +xs] + vp.gp[cx  +xs-2*ps]; 
                de = vp.gp[cx  -xs+2*ps] - 2*vp.gp[cx  -xs] + vp.gp[cx  -xs-2*ps]; 
                df = vp.gp[cx+1+xs+2*ps] - 2*vp.gp[cx+1+xs] + vp.gp[cx+1+xs-2*ps]; 
                dg = vp.gp[cx+1-xs+2*ps] - 2*vp.gp[cx+1-xs] + vp.gp[cx+1-xs-2*ps]; 
                dh = vp.gp[cx-1+xs+2*ps] - 2*vp.gp[cx-1+xs] + vp.gp[cx-1+xs-2*ps]; 
                di = vp.gp[cx-1-xs+2*ps] - 2*vp.gp[cx-1-xs] + vp.gp[cx-1-xs-2*ps]; 
                dj = vp.gp[cx+2   +2*ps] - 2*vp.gp[cx+2]    + vp.gp[cx+2   -2*ps];
                dk = vp.gp[cx-2   +2*ps] - 2*vp.gp[cx-2]    + vp.gp[cx-2   -2*ps];
                dl = vp.gp[cx+2*xs+2*ps] - 2*vp.gp[cx+2*xs] + vp.gp[cx+2*xs-2*ps];
                dm = vp.gp[cx-2*xs+2*ps] - 2*vp.gp[cx-2*xs] + vp.gp[cx-2*xs-2*ps];
                pz.gp[cx] = (8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144.;
            }
        }
    }

    cz = ps;
    for (y=2; y<ys-2; y++) {
        cy = cz + y*xs;
        for (x=2; x<xs-2; x++) {
            cx = cy + x;
            da = vp.gp[cx   +2*ps]   - 2*vp.gp[cx];  
            db = vp.gp[cx+1 +2*ps]   - 2*vp.gp[cx+1]; 
            dc = vp.gp[cx-1 +2*ps]   - 2*vp.gp[cx-1];
            dd = vp.gp[cx+xs+2*ps]   - 2*vp.gp[cx+xs]; 
            de = vp.gp[cx-xs+2*ps]   - 2*vp.gp[cx-xs]; 
            df = vp.gp[cx+1+xs+2*ps] - 2*vp.gp[cx+1+xs]; 
            dg = vp.gp[cx+1-xs+2*ps] - 2*vp.gp[cx+1-xs]; 
            dh = vp.gp[cx-1+xs+2*ps] - 2*vp.gp[cx-1+xs]; 
            di = vp.gp[cx-1-xs+2*ps] - 2*vp.gp[cx-1-xs]; 
            dj = vp.gp[cx+2   +2*ps] - 2*vp.gp[cx+2];
            dk = vp.gp[cx-2   +2*ps] - 2*vp.gp[cx-2];
            dl = vp.gp[cx+2*xs+2*ps] - 2*vp.gp[cx+2*xs];
            dm = vp.gp[cx-2*xs+2*ps] - 2*vp.gp[cx-2*xs];
            pz.gp[cx] = (8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144.;
        }
    }

    cz = (zs-2)*ps;
    for (y=2; y<ys-2; y++) {
        cy = cz + y*xs;
        for (x=2; x<xs-2; x++) {
            cx = cy + x;
            da = - 2*vp.gp[cx]      + vp.gp[cx   -2*ps];     
            db = - 2*vp.gp[cx+1]    + vp.gp[cx+1 -2*ps]; 
            dc = - 2*vp.gp[cx-1]    + vp.gp[cx-1 -2*ps];
            dd = - 2*vp.gp[cx+xs]   + vp.gp[cx+xs-2*ps]; 
            de = - 2*vp.gp[cx-xs]   + vp.gp[cx-xs-2*ps]; 
            df = - 2*vp.gp[cx+1+xs] + vp.gp[cx+1+xs-2*ps]; 
            dg = - 2*vp.gp[cx+1-xs] + vp.gp[cx+1-xs-2*ps]; 
            dh = - 2*vp.gp[cx-1+xs] + vp.gp[cx-1+xs-2*ps]; 
            di = - 2*vp.gp[cx-1-xs] + vp.gp[cx-1-xs-2*ps]; 
            dj = - 2*vp.gp[cx+2]    + vp.gp[cx+2   -2*ps];
            dk = - 2*vp.gp[cx-2]    + vp.gp[cx-2   -2*ps];
            dl = - 2*vp.gp[cx+2*xs] + vp.gp[cx+2*xs-2*ps];
            dm = - 2*vp.gp[cx-2*xs] + vp.gp[cx-2*xs-2*ps];
            pz.gp[cx] = (8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144.;
        }
    }

    return pz;
}



VSGraph  vNabra(WSGraph vp)  // Sobel 
{
    int   i, xs, ys, zs;
    double  xx, yy, zz;
    WSGraph px, py, pz;
    VSGraph pn;
     
    memset(&pn, 0, sizeof(VSGraph));
    if (vp.gp==NULL) {
        pn.state = JBXL_GRAPH_NODATA_ERROR;
        return pn;
    }

    xs = vp.xs;
    ys = vp.ys;
    zs = vp.zs;
    pn = make_VSGraph(xs, ys, zs);
    if (pn.gp==NULL) return pn;

    px = xSobel(vp);
    if (px.gp==NULL) {
        free_VSGraph(&pn);
        pn.state = px.state;
        return pn;
    }
    py = ySobel(vp);
    if (py.gp==NULL) {
        free_VSGraph(&pn);
        free(px.gp);
        pn.state = py.state;
        return pn;
    }

    if (vp.zs<3) {
        for (i=0; i<xs*ys; i++) {
            xx = px.gp[i];
            yy = py.gp[i];
            pn.gp[i] = set_vector(xx, yy, 0.0);
            unit_vector(pn.gp[i]);
        }
    }
    else {
        pz = zSobel(vp);
        if (pz.gp==NULL) {
            free_VSGraph(&pn);
            free(px.gp);
            free(py.gp);
            pn.state = pz.state;
            return pn;
        }

        for (i=0; i<xs*ys*zs; i++) {
            xx = px.gp[i];
            yy = py.gp[i];
            zz = pz.gp[i];
            pn.gp[i] = set_vector(xx, yy, zz);
            unit_vector(pn.gp[i]);
        }
        free(pz.gp);
    }

    free(px.gp);
    free(py.gp);

    return pn;
}



VSGraph  vfNabra(FSGraph vp)  // Sobel 
{
    int    i, xs, ys, zs;
    double xx, yy, zz;
    FSGraph px, py, pz;
    VSGraph pn;
     
    memset(&pn, 0, sizeof(VSGraph));
    if (vp.gp==NULL) {
        pn.state = JBXL_GRAPH_NODATA_ERROR;
        return pn;
    }

    xs = vp.xs;
    ys = vp.ys;
    zs = vp.zs;
    pn = make_VSGraph(xs, ys, zs);
    if (pn.gp==NULL) return pn;

    px = fxSobel(vp);
    if (px.gp==NULL) {
        free_VSGraph(&pn);
        pn.state = px.state;
        return pn;
    }
    py = fySobel(vp);
    if (py.gp==NULL) {
        free_VSGraph(&pn);
        free(px.gp);
        pn.state = py.state;
        return pn;
    }

    if (vp.zs<3) {
        for (i=0; i<xs*ys; i++) {
            xx = px.gp[i];
            yy = py.gp[i];
            pn.gp[i] = set_vector(xx, yy, 0.0);
            unit_vector(pn.gp[i]);
        }
    }
    else {
        pz = fzSobel(vp);
        if (pz.gp==NULL) {
            free_VSGraph(&pn);
            free(px.gp);
            free(py.gp);
            pn.state = pz.state;
            return pn;
        }

        for (i=0; i<xs*ys*zs; i++) {
            xx = px.gp[i];
            yy = py.gp[i];
            zz = pz.gp[i];
            pn.gp[i] = set_vector(xx, yy, zz);
            unit_vector(pn.gp[i]);
        }
        free(pz.gp);
    }

    free(px.gp);
    free(py.gp);

    return pn;
}



WSGraph  Nabra(WSGraph vp)  // Sobel 
{
    int   i, xs, ys, zs;
    int   xx, yy, zz;
    WSGraph px, py, pz, pn;
     
    memset(&pn, 0, sizeof(WSGraph));
    if (vp.gp==NULL) {
        pn.state = JBXL_GRAPH_NODATA_ERROR;
        return pn;
    }

    xs = vp.xs;
    ys = vp.ys;
    zs = vp.zs;
    pn = make_WSGraph(xs, ys, zs);
    if (pn.gp==NULL) return pn;

    px = xSobel(vp);
    if (px.gp==NULL) {
        free_WSGraph(&pn);
        pn.state = px.state;
        return pn;
    }
    py = ySobel(vp);
    if (py.gp==NULL) {
        free_WSGraph(&pn);
        free(px.gp);
        pn.state = py.state;
        return pn;
    }

    if (vp.zs<3) {
        for (i=0; i<xs*ys; i++) {
            xx = px.gp[i];
            yy = py.gp[i];
            pn.gp[i] = (sWord)sqrt(xx*xx + yy*yy);
        }
    }
    else {
        pz = zSobel(vp);
        if (pz.gp==NULL) {
            free_WSGraph(&pn);
            free(px.gp);
            free(py.gp);
            pn.state = pz.state;
            return pn;
        }

        for (i=0; i<xs*ys*zs; i++) {
            xx = px.gp[i];
            yy = py.gp[i];
            zz = pz.gp[i];
            pn.gp[i] = (sWord)sqrt(xx*xx + yy*yy + zz*zz);
        }
        free(pz.gp);
    }

    free(px.gp);
    free(py.gp);

    return pn;
}



FSGraph  fNabra(FSGraph vp)  // Sobel 
{
    int    i, xs, ys, zs;
    double xx, yy, zz;
    FSGraph px, py, pz, pn;
     
    memset(&pn, 0, sizeof(FSGraph));
    if (vp.gp==NULL) {
        pn.state = JBXL_GRAPH_NODATA_ERROR;
        return pn;
    }

    xs = vp.xs;
    ys = vp.ys;
    zs = vp.zs;
    pn = make_FSGraph(xs, ys, zs);
    if (pn.gp==NULL) return pn;

    px = fxSobel(vp);
    if (px.gp==NULL) {
        free_FSGraph(&pn);
        pn.state = px.state;
        return pn;
    }
    py = fySobel(vp);
    if (py.gp==NULL) {
        free_FSGraph(&pn);
        free(px.gp);
        pn.state = py.state;
        return pn;
    }

    if (vp.zs<3) {
        for (i=0; i<xs*ys; i++) {
            xx = px.gp[i];
            yy = py.gp[i];
            pn.gp[i] = sqrt(xx*xx + yy*yy);
        }
    }
    else {
        pz = fzSobel(vp);
        if (pz.gp==NULL) {
            free_FSGraph(&pn);
            free(px.gp);
            free(py.gp);
            pn.state = pz.state;
            return pn;
        }

        for (i=0; i<xs*ys*zs; i++) {
            xx = px.gp[i];
            yy = py.gp[i];
            zz = pz.gp[i];
            pn.gp[i] = sqrt(xx*xx + yy*yy + zz*zz);
        }
        free(pz.gp);
    }

    free(px.gp);
    free(py.gp);

    return pn;
}



VSGraph  curvature3D(FSGraph vp)
{
    int  i;
    double alph, beta, gamm, K, H;
    double fx, fy, fz, fxy, fyz, fzx, fxx, fyy, fzz, nb;
    FSGraph px, py, pz, pxy, pyz, pzx, pxx, pyy, pzz, nab;
    VSGraph pp;

    memset(&pp, 0, sizeof(VSGraph));
    if (vp.gp==NULL) {
        pp.state = JBXL_GRAPH_NODATA_ERROR;
        return pp;
    }

    if (vp.zs<5) {
        //fprintf(stderr,"CURVATURE3D: z dimension is < 5.\n");
        pp.state = JBXL_GRAPH_IVDARG_ERROR;
        return pp;
    }

    pp = make_VSGraph(vp.xs, vp.ys, vp.zs);
    if (pp.gp==NULL) return pp;

    nab = fNabra(vp);
    px  = fxSobel(vp);
    py  = fySobel(vp);
    pz  = fzSobel(vp);
    pxy = fySobel(px);
    pyz = fzSobel(py);
    pzx = fxSobel(pz);
    pxx = fxxSobel(vp);
    pyy = fyySobel(vp);
    pzz = fzzSobel(vp);

    if (nab.gp==NULL || px.gp==NULL  || py.gp==NULL  || pz.gp==NULL  ||
                        pxy.gp==NULL || pyz.gp==NULL || pzx.gp==NULL ||
                        pxx.gp==NULL || pyy.gp==NULL || pzz.gp==NULL) {
        freeNull(px.gp);
        freeNull(py.gp);
        freeNull(pz.gp);
        freeNull(pxy.gp);
        freeNull(pyz.gp);
        freeNull(pzx.gp);
        freeNull(pxx.gp);
        freeNull(pyy.gp);
        freeNull(pzz.gp);
        freeNull(nab.gp);
        free_VSGraph(&pp);
        pp.state = JBXL_GRAPH_ERROR;
        return pp;
    }

    for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
        nb  = nab.gp[i];
        fx  = px.gp[i];
        fy  = py.gp[i];
        fz  = pz.gp[i];
        fxy = pxy.gp[i];
        fyz = pyz.gp[i];
        fzx = pzx.gp[i];
        fxx = pxx.gp[i];
        fyy = pyy.gp[i];
        fzz = pzz.gp[i];

        if (nb*(fx*fx+fy*fy) !=0) {
            alph = (2*fx*fy*fxy - fx*fx*fyy - fy*fy*fxx)/(fx*fx+fy*fy);
            beta = (2*fz*(fx*fx+fy*fy)*(fx*fzx+fy*fyz) - 2*fx*fy*fz*fz*fxy
                        - fx*fx*fz*fz*fxx - fy*fy*fz*fz*fyy 
                        - (fx*fx+fy*fy)*(fx*fx+fy*fy)*fzz)/(nb*nb*(fx*fx+fy*fy));
            gamm = ((fx*fx+fy*fy)*(fy*fzx-fx*fyz) + (fx*fx-fy*fy)*fz*fxy
                        - fx*fy*fz*(fxx-fyy))/(nb*(fx*fx+fy*fy));

            K = alph*beta - gamm*gamm;
            H = (alph + beta)/2;
            pp.gp[i] = set_vector(K, H, 0.0);
        }
    }

    free(px.gp);
    free(py.gp);
    free(pz.gp);
    free(pxy.gp);
    free(pyz.gp);
    free(pzx.gp);
    free(pxx.gp);
    free(pyy.gp);
    free(pzz.gp);
    free(nab.gp);

    return pp;
}



VSGraph  curvature(FSGraph vp)
{
    int   i;
    double K, H, d;
    double Ix, Ixx, Iy, Iyy, Ixy;
    FSGraph px, py, pxy, pxx, pyy;
    VSGraph pp;

    memset(&pp, 0, sizeof(VSGraph));
    if (vp.gp==NULL) {
        pp.state = JBXL_GRAPH_NODATA_ERROR;
        return pp;
    }

    if (vp.zs>1) {
        //fprintf(stderr,"CURVATURE: z dimension is > 1.\n");
        pp.state = JBXL_GRAPH_IVDARG_ERROR;
        return pp;
    }

    pp = make_VSGraph(vp.xs, vp.ys, vp.zs);
    if (pp.gp==NULL) return pp;

    px  = fxSobel(vp);
    py  = fySobel(vp);
    pxy = fySobel(px);
    pxx = fxxSobel(vp);
    pyy = fyySobel(vp);

    if (px.gp==NULL||py.gp==NULL||pxy.gp==NULL||pxx.gp==NULL||pyy.gp==NULL) {
        freeNull(px.gp);
        freeNull(py.gp);
        freeNull(pxy.gp);
        freeNull(pxx.gp);
        freeNull(pyy.gp);
        free_VSGraph(&pp);
        pp.state = JBXL_GRAPH_ERROR;
        return pp;
    }

    for (i=0; i<vp.xs*vp.ys; i++) {
        Ix  = px.gp[i];
        Iy  = py.gp[i];
        Ixy = pxy.gp[i];
        Ixx = pxx.gp[i];
        Iyy = pyy.gp[i];
        d   = 1. + Ix*Ix + Iy*Iy;

        K = (Ixx*Iyy-Ixy*Ixy)/(d*d); 
        H = (Ixx+Ixx*Iy*Iy+Iyy+Iyy*Ix*Ix-2*Ix*Ixy*Iy)/(2.*d*sqrt(d));
        pp.gp[i] = set_vector(K, H, 0.0);
    }

    free(px.gp);
    free(py.gp);
    free(pxy.gp);
    free(pxx.gp);
    free(pyy.gp);

    return pp;
}



WSGraph  curv2WSGraph(VSGraph xp)
{
    int  i;
    WSGraph vp;
    double  K, H;
    
    memset(&vp, 0, sizeof(WSGraph));
    if (xp.gp==NULL) {
        vp.state = JBXL_GRAPH_NODATA_ERROR;
        return vp;
    }

    vp = make_WSGraph(xp.xs, xp.ys, xp.zs);
    if (vp.gp==NULL) return vp;

    for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
        K = xp.gp[i].x;
        H = xp.gp[i].y;
        if      (K>0  && H<0)  vp.gp[i] = PEAK;
        else if (K>0  && H>0)  vp.gp[i] = PIT;
        else if (K<0  && H<0)  vp.gp[i] = SADDLE_RIDGE;
        else if (K<0  && H>0)  vp.gp[i] = SADDLE_VALLEY;
        else if (K>0  && H==0) vp.gp[i] = NONE_SHAPE;
        else if (K<0  && H==0) vp.gp[i] = MINIMAL;
        else if (K==0 && H<0)  vp.gp[i] = RIDGE;
        else if (K==0 && H>0)  vp.gp[i] = VALLEY;
        else if (K==0 && H==0) vp.gp[i] = FLAT;
    }

    return vp;
}



WSGraph  WSCurve(WSGraph gx, int mode, int cc)
{
    int  i;
    FSGraph wp;
    VSGraph xp;
    WSGraph vp;

    memset(&vp, 0, sizeof(WSGraph));
    if (gx.gp==NULL) {
        vp.state = JBXL_GRAPH_NODATA_ERROR;
        return vp;
    }

    wp = W2FSGraph(gx);
    if (gx.zs<5) xp = curvature(wp);
    else         xp = curvature3D(wp);

    vp = curv2WSGraph(xp);
    freeNull(wp.gp);
    freeNull(xp.gp);
    if (vp.gp==NULL) return vp;

    if (mode==ALL) {
        for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
            if    (vp.gp[i]==FLAT)            vp.gp[i] = 500;
            else if (vp.gp[i]==PIT)           vp.gp[i] = 1000;
            else if (vp.gp[i]==SADDLE_RIDGE)  vp.gp[i] = 1500;
            else if (vp.gp[i]==SADDLE_VALLEY) vp.gp[i] = 2000;
            else if (vp.gp[i]==NONE_SHAPE)    vp.gp[i] = 2500;
            else if (vp.gp[i]==MINIMAL)       vp.gp[i] = 3000;
            else if (vp.gp[i]==RIDGE)         vp.gp[i] = 3500;
            else if (vp.gp[i]==VALLEY)        vp.gp[i] = 4000;
            else if (vp.gp[i]==PEAK)          vp.gp[i] = 4500;
            else                              vp.gp[i] = 0;
        }
    }
    else {
        for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
            if ((vp.gp[i]&mode)!=0) vp.gp[i] = cc;
            else                    vp.gp[i] = 0;
        }
    }
    return vp;
}



WSGraph  edge_enhance(WSGraph gd, int mode)
{
    int  i;
    WSGraph la, vp;

    memset(&vp, 0, sizeof(WSGraph));
    if (gd.gp==NULL) {
        vp.state = JBXL_GRAPH_NODATA_ERROR;
        return vp;
    }

    la = Laplacian(gd, mode);  
    if (la.gp==NULL) return la;

    vp = make_WSGraph(gd.xs, gd.ys, 1);
    if (vp.gp==NULL) return vp;

    for (i=0; i<vp.xs*vp.ys; i++) vp.gp[i] = gd.gp[i] - la.gp[i];
    return vp;
}




// Filter

FMask  gauss_mask(double sig, int ms, int md) 
{
    int    xx, yy, zz, ns, cp, dx, ps, sw;
    double min, *fm;
    FMask mask;
    
    mask.mode  = 0;
    mask.msize = 0;
    mask.imask = NULL;

    if (md<=2) {    // 2D 
        md = 2;
        ps = ms*ms;
        sw = 0;
    }
    else {          // 3D 
        md = 3;
        ps = ms*ms*ms;
        sw = 1;
    }

    ns  = ms/2;
    min = (double)SINTMAX;
    fm = (double*)malloc(ps*sizeof(double));
    mask.imask = (int*)malloc(ps*sizeof(int));
    if (fm==NULL || mask.imask==NULL) {
        free(fm);
        free(mask.imask);
        memset(&mask, 0, sizeof(FMask));
        return mask;
    }

    for (zz=-ns*sw; zz<=ns*sw; zz++) {
        for (yy=-ns; yy<=ns; yy++) {
            for (xx=-ns; xx<=ns; xx++) {
                cp = (zz+ns)*ms*ms*sw + (yy+ns)*ms + (xx+ns);
                fm[cp] = exp(-(xx*xx+yy*yy+zz*zz)/(sig*sig));
                if (fm[cp]!=0.0) min = Min(min, fm[cp]);
            }
        }
    }

    dx = 0;
    for (xx=0; xx<ps; xx++) {
        mask.imask[xx] = (int)(fm[xx]/min+0.5);
        dx += mask.imask[xx];
    }

    mask.msize = ms;
    mask.nfact = dx;
    mask.mode  = md;

    free(fm);
    return mask; 
}



WSGraph  imask(WSGraph xp, FMask mask)
{
    int    i, x, y, z, cx;
    int    xx, yy, zz, cp, cw, sw;
    int    kc, xc, xs, ps, pm, mz, zc;
    int    ms, nf, min;
    double dd;
    WSGraph vp;

    memset(&vp, 0, sizeof(WSGraph));
    if (xp.gp==NULL) {
        vp.state = JBXL_GRAPH_NODATA_ERROR;
        return vp;
    }

    if (xp.zs<=1 && mask.mode>2) {
        //fprintf(stderr, "IMASK: mismach mask dimension %d %d\n", xp.zs, mask.mode);
        vp.state = JBXL_GRAPH_IVDARG_ERROR;
        return vp;
    }

    nf = mask.nfact;
    ms = mask.msize;
    if (mask.mode==2) {
        sw = 0;
        pm = ms*ms;
    }
    else {
        sw = 1;
        pm = ms*ms*ms;
    }

    mz = Min(ms, xp.zs);
    kc = pm/2;
    zc = mz/2;
    xc = ms/2;
    xs = xp.xs;
    ps = xp.xs*xp.ys;

    min = SINTMAX;
    for (i=0; i<xp.xs*xp.ys; i++) min = Min(min, xp.gp[i]);
    vp = make_WSGraph(xp.xs, xp.ys, 1);
    if (vp.gp==NULL) return vp;

    for (i=0; i<vp.xs*vp.ys; i++) vp.gp[i] = min;

    z = xp.zs/2; 
    for (y=xc; y<xp.ys-xc; y++) 
    for (x=xc; x<xp.xs-xc; x++) {
        cx = z*ps + y*xs + x;
        dd = 0.0;
        for (zz=-zc*sw; zz<=zc*sw; zz++)
        for (yy=-xc; yy<=xc; yy++)
        for (xx=-xc; xx<=xc; xx++) {
            cp = kc + xx + yy*ms + zz*ms*ms;
            cw = cx + xx + yy*xs + zz*ps;
            dd = dd + (double)xp.gp[cw]*mask.imask[cp];
        }
        vp.gp[y*xs + x] = (sWord)(dd/nf);
    }
    return vp;
}



WSGraph  median(WSGraph xp, int ms)  /* 3D */
{
    int   i, j, x, y, z;
    int   xx, yy, zz, cw, ux, mz;
    int   kc, xc, zc, xs, ps, cx;
    WSGraph vp;
    sWord   *me;

    memset(&vp, 0, sizeof(WSGraph));
    if (xp.gp==NULL) {
        vp.state = JBXL_GRAPH_NODATA_ERROR;
        return vp;
    }

    mz = Min(ms, xp.zs);
    me = (sWord*)malloc(ms*ms*mz*sizeof(sWord));
    if (me==NULL) {
        vp.state = JBXL_GRAPH_MEMORY_ERROR;
        return vp;
    }

    kc = ms*ms*mz/2;
    xc = ms/2;
    zc = mz/2;
    xs = xp.xs;
    ps = xp.xs*xp.ys;
    vp = make_WSGraph(vp.xs, vp.ys, vp.zs);
    if (vp.gp==NULL) {
        free(me);
        return vp;
    }
    
    z = xp.zs/2;
    for(y=xc; y<xp.ys-xc; y++) {
        for(x=xc; x<xp.xs-xc; x++) {
            cx = z*ps + y*xs + x;
            i  = 0;
            for (zz=-zc; zz<=zc; zz++) {
                for (yy=-xc; yy<=xc; yy++) {
                    for (xx=-xc; xx<=xc; xx++) {
                        cw = cx + xx + yy*xs + zz*ps;
                        me[i++] = xp.gp[cw];
                    }
                }
            }

            for (i=0; i<ms*ms*mz-1; i++) {
                for (j=i+1; j<ms*ms*mz; j++) {
                    if (me[i]<me[j]) {
                        ux  = me[i];
                        me[i] = me[j];
                        me[j] = ux;
                    }
                }
            }
            vp.gp[cx-z*ps] = me[kc];
        }
    }

    free(me);
    return vp;
}



WSGraph  to2d(WSGraph gd, int mode)
{
    int  i, j, k, psize;
    int  cx, cy, cz, cw, xx, yy;
    WSGraph vp;

    memset(&vp, 0, sizeof(WSGraph));
    if (gd.gp==NULL) {
        vp.state = JBXL_GRAPH_NODATA_ERROR;
        return vp;
    }

    psize = gd.xs*gd.ys;

    if (mode==TOP_VIEW) {
        vp = make_WSGraph(gd.xs, gd.ys, 1);
        if (vp.gp==NULL) return vp;

        for (k=0; k<gd.zs; k++) {
            cz = k*psize;
            for (j=0; j<gd.ys; j++) {
                cy = j*gd.xs;
                for (i=0; i<gd.xs; i++) {
                    cx = cz + cy + i;
                    cw = cy + i;
                    if (gd.gp[cx]!=0) vp.gp[cw] = Max(vp.gp[cw], gd.gp[cx]);
                }
            }
        }
    }

    else if (mode==TOP_VIEW_DEPTH) {
        vp = make_WSGraph(gd.xs, gd.ys, 1);
        if (vp.gp==NULL) return vp;

        for (k=0; k<gd.zs; k++) {
            cz = k*psize;
            for (j=0; j<gd.ys; j++) {
                cy = j*gd.xs;
                for (i=0; i<gd.xs; i++) {
                    cx = cz + cy + i;
                    cw = cy + i;
                    if (gd.gp[cx]!=0) vp.gp[cw] = Max(vp.gp[cw], (gd.zs-k)+100);
                }
            }
        }
    }

    else if (mode==SIDEX_VIEW) {
        vp = make_WSGraph(gd.ys, gd.zs, 1);
        if (vp.gp==NULL) return vp;

        for (k=0; k<gd.zs; k++) {
            cz = k*psize;
            yy = k;
            for (j=0; j<gd.ys; j++) {
                cy = j*gd.xs;
                xx = gd.ys - 1 - j;
                for (i=0; i<gd.xs; i++) {
                    cx = cz + cy + i;
                    cw = yy*vp.xs + xx;
                    if (gd.gp[cx]!=0) vp.gp[cw] = Max(vp.gp[cw], gd.gp[cx]);
                }
            }
        }
    }

    else if (mode==SIDEY_VIEW) {
        vp = make_WSGraph(gd.xs, gd.zs, 1);
        if (vp.gp==NULL) return vp;

        for (k=0; k<gd.zs; k++) {
            cz = k*psize;
            yy = k;
            for (j=0; j<gd.ys; j++) {
                cy = j*gd.xs;
                for (i=0; i<gd.xs; i++) {
                    cx = cz + cy + i;
                    xx = i;
                    cw = yy*vp.xs + xx;
                    if (gd.gp[cx]!=0) vp.gp[cw] = Max(vp.gp[cw], gd.gp[cx]);
                }
            }
        }
    }

    else {
        memset(&vp, 0, sizeof(WSGraph));
        vp.state = JBXL_GRAPH_IVDARG_ERROR;
        //fprintf(stderr,"TO2D: unknown mode = %d.\n",mode);
    }

    return vp;
}



WSGraph  euclid_distance(WSGraph vp, int* rr, int bc)
{
    int  i, j, k, l;
    int  df, db, d, w; 
    int  rmax, rstart, rend;
    WSGraph wp;
    ISGraph pp, buff;

    memset(&wp, 0, sizeof(WSGraph));
    if (vp.gp==NULL) {
        wp.state = JBXL_GRAPH_NODATA_ERROR;
        return wp;
    }

    wp = make_WSGraph(vp.xs, vp.ys, vp.zs);
    if (wp.gp==NULL) return wp;

    pp = make_ISGraph(vp.xs, vp.ys, vp.zs);
    if (pp.gp==NULL) {
        free_WSGraph(&wp);
        wp.state = pp.state;
        return wp;
    }

    for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
         if (vp.gp[i]>=bc) pp.gp[i] = 1;
         else              pp.gp[i] = 0;
    }

    for (k=1; k<=vp.zs; k++) {
        for (j=1; j<=vp.ys; j++) {
            df = vp.xs;
            for (i=1; i<=vp.xs; i++) {
                if (Vxt(pp, i, j, k)!=0) df = df + 1;
                else                     df = 0;
                Vxt(pp, i, j, k) = df*df;
            }
        }
    }
    
    for (k=1; k<=vp.zs; k++) {
        for (j=1; j<=vp.ys; j++) {
            db = vp.xs;
            for (i=vp.xs; i>=1; i--) {
                if (Vxt(pp, i, j, k)!=0) db = db + 1;
                else                     db = 0;
                Vxt(pp, i, j, k) = Min(Vxt(pp, i, j, k), db*db);
            }
        }
    }

    buff = make_ISGraph(vp.ys, 1, 1);
    for (k=1; k<=vp.zs; k++) {
        for (i=1; i<=vp.xs; i++) {
            for (j=1; j<=vp.ys; j++) {
                Lxt(buff, j) = Vxt(pp, i, j, k);
            }
            for (j=1; j<=vp.ys; j++) {
                d = Lxt(buff, j);
                if (d!=0) {
                    rmax   = (int)sqrt((double)d) + 1;
                    rstart = Min(rmax, j-1);
                    rend   = Min(rmax, vp.ys-j);
                    for (l=-rstart; l<=rend; l++) {
                        w = Lxt(buff, j+l) + l*l;
                        if (w<d) d = w;
                    }
                }
                Vxt(pp, i, j, k) = d;
            }
        }
    }
    free(buff.gp);

    *rr = 0;
    buff = make_ISGraph(vp.zs, 1, 1);
    for (j=1; j<=vp.ys; j++) {
        for (i=1; i<=vp.xs; i++) {
            for (k=1; k<=vp.zs; k++) {
                Lxt(buff, k) = Vxt(pp, i, j, k);
            }
            for (k=1; k<=vp.zs; k++) {
                d = Lxt(buff, k);
                if (d!=0) {
                    rmax   = (int)sqrt((double)d) + 1;
                    rstart = Min(rmax, k-1);
                    rend   = Min(rmax, vp.zs-k);
                    for (l=-rstart; l<=rend; l++) {
                        w = Lxt(buff, k+l) + l*l;
                        if (w<d) d = w;
                    }
                    *rr = Max(*rr, d);
                }
                Vxt(pp, i, j, k) = d;
            }
        }
    }
    free(buff.gp);

    for (i=0; i<wp.xs*wp.ys*wp.zs; i++) {
        wp.gp[i] = (sWord)pp.gp[i];
        if (pp.gp[i]>32767) {
            fprintf(stderr,"EUCLID_DISTANCE: WARNING: distance is too long = %d!\n",pp.gp[i]);
        }
    }
    free(pp.gp);

    return wp;
}



int  out_round(WSGraph vp, int x, int y, IRBound* rb, int mode)
{
    int  i, j, sp, cp, w, ll, ss;
    int  xx, yy, vx, vy, ix, eflg=OFF;
    int  r8[8]={-1, 1, -1, -1, 1, -1, 1, 1};
    int  r4[8]={ 0, 1, -1,  0, 0, -1, 1, 0};
    int* cc;

    if (vp.gp==NULL) return JBXL_GRAPH_IVDARG_ERROR;

    i = y*vp.xs + x;
    rb->xmax = rb->xmin = x;
    rb->ymax = rb->ymin = y;
    sp = cp = i;
    ss = 0;
    ll = 0;
    vx = 1;
    vy = 0;

    if (vp.gp[sp]==0 || sp==0) {
        //fprintf(stderr,"OUT_ROUND: irregular start point!! sp = %d\n",sp);
        return JBXL_GRAPH_IVDPARAM_ERROR;
    }

    if (mode==8){
        ix = 8;
        cc = r8;
    }
    else if (mode==4) {
        ix = 4;
        cc = r4;
    }
    else {
        //fprintf(stderr,"OUT_ROUND: invalid mode = %d!!\n",mode);
        return JBXL_GRAPH_IVDMODE_ERROR;
    }

    do {
        w  = abs(vx)+abs(vy);
        xx = (vx*cc[0]+vy*cc[1])/w;
        yy = (vx*cc[2]+vy*cc[3])/w;
        for (j=1; j<=ix; j++) {
            if (vp.gp[cp+yy*vp.xs+xx]!=0) {
                vx = xx;
                vy = yy;
                cp = cp + yy*vp.xs + xx;
                xx = cp%vp.xs;
                yy = (cp-xx)/vp.xs;
                rb->xmax = Max(rb->xmax, xx);
                rb->ymax = Max(rb->ymax, yy);
                rb->xmin = Min(rb->xmin, xx);
                rb->ymin = Min(rb->ymin, yy);
                break;
            }
            else {
                if(sp==cp && xx==-1 && yy==0) {
                    eflg = ON;
                    break;
                }
                w  = abs(xx)+abs(yy);
                vx = (xx*cc[4]+yy*cc[5])/w;
                vy = (xx*cc[6]+yy*cc[7])/w;
                xx = vx;
                yy = vy;
            }
        }
        ll++;
        if (abs(vx)+abs(vy)==2) ss++;
        //
    } while(eflg==OFF);

    if (mode==4) ss = 0;
    (rb->xmax)++;
    (rb->ymax)++;
    (rb->xmin)--;
    (rb->ymin)--;
    rb->misc = ss;

    return ll;
}

---