gvl_calc.c

Go to the documentation of this file.

/*!
   \file lib/ogsf/gvl_calc.c
 
   \brief OGSF library - loading and manipulating volumes (lower level
   functions)
 
   GRASS OpenGL gsurf OGSF Library
 
   (C) 1999-2008 by the GRASS Development Team
 
   This program is free software under the
   GNU General Public License (>=v2).
   Read the file COPYING that comes with GRASS
   for details.
 
   \author Tomas Paudits (February 2004)
   \author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
 */
 
#include <math.h>
 
#include <grass/gis.h>
#include <grass/ogsf.h>
 
#include "rgbpack.h"
#include "mc33_table.h"
 
/*!
   \brief memory buffer for writing
 */
#define BUFFER_SIZE      1000000
 
/* USEFUL MACROS */
 
/* interp. */
#define LINTERP(d, a, b) (a + d * (b - a))
#define TINTERP(d, v)                                                   \
    ((v[0] * (1. - d[0]) * (1. - d[1]) * (1. - d[2])) +                 \
     (v[1] * d[0] * (1. - d[1]) * (1. - d[2])) +                        \
     (v[2] * d[0] * d[1] * (1. - d[2])) +                               \
     (v[3] * (1. - d[0]) * d[1] * (1. - d[2])) +                        \
     (v[4] * (1. - d[0]) * (1. - d[1]) * d[2]) +                        \
     (v[5] * d[0] * (1. - d[1]) * d[2]) + (v[6] * d[0] * d[1] * d[2]) + \
     (v[7] * (1. - d[0]) * d[1] * d[2]))
 
#define FOR_VAR i_for
#define FOR_0_TO_N(n, cmd)                          \
    {                                               \
        int FOR_VAR;                                \
        for (FOR_VAR = 0; FOR_VAR < n; FOR_VAR++) { \
            cmd;                                    \
        }                                           \
    }
 
/*!
   \brief writing and reading isosurface data
 */
#define WRITE(c)         gvl_write_char(dbuff->ndx_new++, &(dbuff->new), c)
#define READ()           gvl_read_char(dbuff->ndx_old++, dbuff->old)
#define SKIP(n)          dbuff->ndx_old = dbuff->ndx_old + n
 
/*!
   \brief check and set data descriptor
 */
#define IS_IN_DATA(att)  ((isosurf->data_desc >> att) & 1)
#define SET_IN_DATA(att) isosurf->data_desc = (isosurf->data_desc | (1 << att))
 
typedef struct {
    unsigned char *old;
    unsigned char *new;
    int ndx_old;
    int ndx_new;
    int num_zero;
} data_buffer;
 
int mc33_process_cube(int c_ndx, float *v);
 
/* global variables */
int Rows, Cols, Depths;
double ResX, ResY, ResZ;
 
/************************************************************************/
/* ISOSURFACES */
 
/*!
   \brief Write cube index
 
   \param ndx
   \param dbuff
 */
void iso_w_cndx(int ndx, data_buffer *dbuff)
{
    /* cube don't contains polys */
    if (ndx == -1) {
        if (dbuff->num_zero == 0) {
            WRITE(0);
            dbuff->num_zero++;
        }
        else if (dbuff->num_zero == 254) {
            WRITE(dbuff->num_zero + 1);
            dbuff->num_zero = 0;
        }
        else {
            dbuff->num_zero++;
        }
    }
    else { /* isosurface cube */
        if (dbuff->num_zero == 0) {
            WRITE((ndx / 256) + 1);
            WRITE(ndx % 256);
        }
        else {
            WRITE(dbuff->num_zero);
            dbuff->num_zero = 0;
            WRITE((ndx / 256) + 1);
            WRITE(ndx % 256);
        }
    }
}
 
/*!
   \brief Read cube index
 
   \param dbuff
 */
int iso_r_cndx(data_buffer *dbuff)
{
    int ndx, ndx2;
 
    if (dbuff->num_zero != 0) {
        dbuff->num_zero--;
        ndx = -1;
    }
    else {
        WRITE(ndx = READ());
        if (ndx == 0) {
            WRITE(dbuff->num_zero = READ());
            dbuff->num_zero--;
            ndx = -1;
        }
        else {
            WRITE(ndx2 = READ());
            ndx = (ndx - 1) * 256 + ndx2;
        }
    }
 
    return ndx;
}
 
/*!
   \brief Get value from data input
 
   \param isosurf
   \param desc
   \param x,y,z
   \param[out] v value
 
   \return 0
   \return ?
 */
int iso_get_cube_value(geovol_isosurf *isosurf, int desc, int x, int y, int z,
                       float *v)
{
    double d;
    geovol_file *vf;
    int type, ret = 1;
 
    /* get volume file from attribute handle */
    vf = gvl_file_get_volfile(isosurf->att[desc].hfile);
    type = gvl_file_get_data_type(vf);
 
    /* get value from volume file */
    if (type == VOL_DTYPE_FLOAT) {
        gvl_file_get_value(vf, (int)(x * ResX), (int)(y * ResY),
                           (int)(z * ResZ), v);
    }
    else if (type == VOL_DTYPE_DOUBLE) {
        gvl_file_get_value(vf, (int)(x * ResX), (int)(y * ResY),
                           (int)(z * ResZ), &d);
        *v = (float)d;
    }
    else {
        return 0;
    }
 
    /* null check */
    if (gvl_file_is_null_value(vf, v))
        ret = 0;
 
    /* adjust data */
    switch (desc) {
    case (ATT_TOPO):
        *v = (*v) - isosurf->att[desc].constant;
        break;
    case (ATT_MASK):
        if (isosurf->att[desc].constant)
            ret = !ret;
        break;
    }
 
    return ret;
}
 
/*!
   \brief Get volume file values range
 
   \param isosurf
   \param desc
   \param[out] min
   \param[out] max
 */
void iso_get_range(geovol_isosurf *isosurf, int desc, double *min, double *max)
{
    gvl_file_get_min_max(gvl_file_get_volfile(isosurf->att[desc].hfile), min,
                         max);
}
 
/*!
   \brief Read values for cube
 
   \param isosurf
   \param desc
   \param x,y,z
   \param[out] v
 
   \return
 */
int iso_get_cube_values(geovol_isosurf *isosurf, int desc, int x, int y, int z,
                        float *v)
{
    int p, ret = 1;
 
    for (p = 0; p < 8; ++p) {
        if (iso_get_cube_value(isosurf, desc, x + ((p ^ (p >> 1)) & 1),
                               y + ((p >> 1) & 1), z + ((p >> 2) & 1),
                               &v[p]) == 0) {
            ret = 0;
        }
    }
 
    return ret;
}
 
/*!
   \brief Calculate cube grads
 
   \param isosurf
   \param x,y,z
   \param grad
 */
void iso_get_cube_grads(geovol_isosurf *isosurf, int x, int y, int z,
                        float (*grad)[3])
{
    float v[3];
    int i, j, k, p;
 
    for (p = 0; p < 8; ++p) {
        i = x + ((p ^ (p >> 1)) & 1);
        j = y + ((p >> 1) & 1);
        k = z + ((p >> 2) & 1);
 
        /* x */
        if (i == 0) {
            iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
            iso_get_cube_value(isosurf, ATT_TOPO, i + 1, j, k, &v[2]);
            grad[p][0] = v[2] - v[1];
        }
        else {
            if (i == (Cols - 1)) {
                iso_get_cube_value(isosurf, ATT_TOPO, i - 1, j, k, &v[0]);
                iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
                grad[p][0] = v[1] - v[0];
            }
            else {
                iso_get_cube_value(isosurf, ATT_TOPO, i - 1, j, k, &v[0]);
                iso_get_cube_value(isosurf, ATT_TOPO, i + 1, j, k, &v[2]);
                grad[p][0] = (v[2] - v[0]) / 2;
            }
        }
 
        /* y */
        if (j == 0) {
            iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
            iso_get_cube_value(isosurf, ATT_TOPO, i, j + 1, k, &v[2]);
            grad[p][1] = v[2] - v[1];
        }
        else {
            if (j == (Rows - 1)) {
                iso_get_cube_value(isosurf, ATT_TOPO, i, j - 1, k, &v[0]);
                iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
                grad[p][1] = v[1] - v[0];
            }
            else {
                iso_get_cube_value(isosurf, ATT_TOPO, i, j - 1, k, &v[0]);
                iso_get_cube_value(isosurf, ATT_TOPO, i, j + 1, k, &v[2]);
                grad[p][1] = (v[2] - v[0]) / 2;
            }
        }
 
        /* z */
        if (k == 0) {
            iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
            iso_get_cube_value(isosurf, ATT_TOPO, i, j, k + 1, &v[2]);
            grad[p][2] = v[2] - v[1];
        }
        else {
            if (k == (Depths - 1)) {
                iso_get_cube_value(isosurf, ATT_TOPO, i, j, k - 1, &v[0]);
                iso_get_cube_value(isosurf, ATT_TOPO, i, j, k, &v[1]);
                grad[p][2] = v[1] - v[0];
            }
            else {
                iso_get_cube_value(isosurf, ATT_TOPO, i, j, k - 1, &v[0]);
                iso_get_cube_value(isosurf, ATT_TOPO, i, j, k + 1, &v[2]);
                grad[p][2] = (v[2] - v[0]) / 2;
            }
        }
    }
}
 
/*!
   \brief Process cube
 
   \param isosurf
   \param x,y,z
   \param dbuff
 */
void iso_calc_cube(geovol_isosurf *isosurf, int x, int y, int z,
                   data_buffer *dbuff)
{
    int i, c_ndx;
    int crnt, v1, v2, c;
    float val[MAX_ATTS][8], grad[8][3];
    float d, d3[3], d_sum[3], n[3], n_sum[3], tv;
    double min, max;
 
    if (isosurf->att[ATT_TOPO].changed) {
        /* read topo values, if there are NULL values then return */
        if (!iso_get_cube_values(isosurf, ATT_TOPO, x, y, z, val[ATT_TOPO])) {
            iso_w_cndx(-1, dbuff);
            return;
        }
 
        /* mask */
        if (isosurf->att[ATT_MASK].att_src == MAP_ATT) {
            if (!iso_get_cube_values(isosurf, ATT_MASK, x, y, z,
                                     val[ATT_MASK])) {
                iso_w_cndx(-1, dbuff);
                return;
            }
        }
 
        /* index to precalculated table */
        c_ndx = 0;
        for (i = 0; i < 8; i++) {
            if (val[ATT_TOPO][i] > 0)
                c_ndx |= 1 << i;
        }
        c_ndx = mc33_process_cube(c_ndx, val[ATT_TOPO]);
 
        iso_w_cndx(c_ndx, dbuff);
 
        if (c_ndx == -1)
            return;
 
        /* calc cube grads */
        iso_get_cube_grads(isosurf, x, y, z, grad);
    }
    else {
        /* read cube index */
        if ((c_ndx = iso_r_cndx(dbuff)) == -1)
            return;
    }
 
    /* get color values */
    if (isosurf->att[ATT_COLOR].changed &&
        isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
        iso_get_cube_values(isosurf, ATT_COLOR, x, y, z, val[ATT_COLOR]);
    }
 
    /* get transparency values */
    if (isosurf->att[ATT_TRANSP].changed &&
        isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
        iso_get_cube_values(isosurf, ATT_TRANSP, x, y, z, val[ATT_TRANSP]);
    }
 
    /* get shine values */
    if (isosurf->att[ATT_SHINE].changed &&
        isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
        iso_get_cube_values(isosurf, ATT_SHINE, x, y, z, val[ATT_SHINE]);
    }
 
    /* get emit values */
    if (isosurf->att[ATT_EMIT].changed &&
        isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
        iso_get_cube_values(isosurf, ATT_EMIT, x, y, z, val[ATT_EMIT]);
    }
 
    FOR_0_TO_N(3, d_sum[FOR_VAR] = 0.; n_sum[FOR_VAR] = 0.);
 
    /* loop in edges */
    for (i = 0; i < cell_table[c_ndx].nedges; i++) {
        /* get edge number */
        crnt = cell_table[c_ndx].edges[i];
 
        /* set topo */
        if (isosurf->att[ATT_TOPO].changed) {
            /* interior vertex */
            if (crnt == 12) {
                FOR_0_TO_N(3, WRITE((d3[FOR_VAR] =
                                         d_sum[FOR_VAR] /
                                         ((float)(cell_table[c_ndx].nedges))) *
                                    255));
                GS_v3norm(n_sum);
                FOR_0_TO_N(3, WRITE((n_sum[FOR_VAR] /
                                         ((float)(cell_table[c_ndx].nedges)) +
                                     1.) *
                                    127));
                /* edge vertex */
            }
            else {
                /* set edges verts */
                v1 = edge_vert[crnt][0];
                v2 = edge_vert[crnt][1];
 
                /* calc intersection point - edge and isosurf */
                d = val[ATT_TOPO][v1] / (val[ATT_TOPO][v1] - val[ATT_TOPO][v2]);
 
                d_sum[edge_vert_pos[crnt][0]] += d;
                d_sum[edge_vert_pos[crnt][1]] += edge_vert_pos[crnt][2];
                d_sum[edge_vert_pos[crnt][3]] += edge_vert_pos[crnt][4];
 
                WRITE(d * 255);
 
                /* set normal for intersect. point */
                FOR_0_TO_N(3, n[FOR_VAR] = LINTERP(d, grad[v1][FOR_VAR],
                                                   grad[v2][FOR_VAR]));
                GS_v3norm(n);
                FOR_0_TO_N(3, n_sum[FOR_VAR] += n[FOR_VAR]);
                FOR_0_TO_N(3, WRITE((n[FOR_VAR] + 1.) * 127));
            }
        }
        else {
            /* read x,y,z of intersection point in cube coords */
            if (crnt == 12) {
                WRITE(c = READ());
                d3[0] = ((float)c) / 255.0;
                WRITE(c = READ());
                d3[1] = ((float)c) / 255.0;
                WRITE(c = READ());
                d3[2] = ((float)c) / 255.0;
            }
            else {
                /* set edges verts */
                v1 = edge_vert[crnt][0];
                v2 = edge_vert[crnt][1];
 
                WRITE(c = READ());
                d = ((float)c) / 255.0;
            }
 
            /* set normals */
            FOR_0_TO_N(3, WRITE(READ()));
        }
 
        /* set color */
        if (isosurf->att[ATT_COLOR].changed &&
            isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
            if (crnt == 12) {
                tv = TINTERP(d3, val[ATT_COLOR]);
            }
            else {
                tv = LINTERP(d, val[ATT_COLOR][v1], val[ATT_COLOR][v2]);
            }
 
            c = Gvl_get_color_for_value(isosurf->att[ATT_COLOR].att_data, &tv);
 
            WRITE(c & RED_MASK);
            WRITE((c & GRN_MASK) >> 8);
            WRITE((c & BLU_MASK) >> 16);
 
            if (IS_IN_DATA(ATT_COLOR))
                SKIP(3);
        }
        else {
            if (isosurf->att[ATT_COLOR].att_src == MAP_ATT) {
                FOR_0_TO_N(3, WRITE(READ()));
            }
            else {
                if (IS_IN_DATA(ATT_COLOR))
                    SKIP(3);
            }
        }
 
        /* set transparency */
        if (isosurf->att[ATT_TRANSP].changed &&
            isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
            if (crnt == 12) {
                tv = TINTERP(d3, val[ATT_TRANSP]);
            }
            else {
                tv = LINTERP(d, val[ATT_TRANSP][v1], val[ATT_TRANSP][v2]);
            }
 
            iso_get_range(isosurf, ATT_TRANSP, &min, &max);
            c = (min != max) ? 255 - (tv - min) / (max - min) * 255 : 0;
 
            WRITE(c);
            if (IS_IN_DATA(ATT_TRANSP))
                SKIP(1);
        }
        else {
            if (isosurf->att[ATT_TRANSP].att_src == MAP_ATT) {
                WRITE(READ());
            }
            else {
                if (IS_IN_DATA(ATT_TRANSP))
                    SKIP(1);
            }
        }
 
        /* set shin */
        if (isosurf->att[ATT_SHINE].changed &&
            isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
            if (crnt == 12) {
                tv = TINTERP(d3, val[ATT_SHINE]);
            }
            else {
                tv = LINTERP(d, val[ATT_SHINE][v1], val[ATT_SHINE][v2]);
            }
 
            iso_get_range(isosurf, ATT_SHINE, &min, &max);
            c = (min != max) ? (tv - min) / (max - min) * 255 : 0;
 
            WRITE(c);
            if (IS_IN_DATA(ATT_SHINE))
                SKIP(1);
        }
        else {
            if (isosurf->att[ATT_SHINE].att_src == MAP_ATT) {
                WRITE(READ());
            }
            else {
                if (IS_IN_DATA(ATT_SHINE))
                    SKIP(1);
            }
        }
 
        /* set emit */
        if (isosurf->att[ATT_EMIT].changed &&
            isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
            if (crnt == 12) {
                tv = TINTERP(d3, val[ATT_EMIT]);
            }
            else {
                tv = LINTERP(d, val[ATT_EMIT][v1], val[ATT_EMIT][v2]);
            }
 
            iso_get_range(isosurf, ATT_EMIT, &min, &max);
            c = (min != max) ? (tv - min) / (max - min) * 255 : 0;
 
            WRITE(c);
            if (IS_IN_DATA(ATT_SHINE))
                SKIP(1);
        }
        else {
            if (isosurf->att[ATT_EMIT].att_src == MAP_ATT) {
                WRITE(READ());
            }
            else {
                if (IS_IN_DATA(ATT_EMIT))
                    SKIP(1);
            }
        }
    }
}
 
/*!
   \brief Fill data structure with computed isosurfaces polygons
 
   \param gvol pointer to geovol struct
 
   \return 1
 */
int gvl_isosurf_calc(geovol *gvol)
{
    int x, y, z;
    int i, a, read;
    geovol_file *vf;
    geovol_isosurf *isosurf;
 
    data_buffer *dbuff;
    int *need_update, need_update_global;
 
    dbuff = G_malloc(gvol->n_isosurfs * sizeof(data_buffer));
    need_update = G_malloc(gvol->n_isosurfs * sizeof(int));
 
    /* flag - changed any isosurface */
    need_update_global = 0;
 
    /* initialize */
    for (i = 0; i < gvol->n_isosurfs; i++) {
        isosurf = gvol->isosurf[i];
 
        /* initialize read/write buffers */
        dbuff[i].old = NULL;
        dbuff[i].new = NULL;
        dbuff[i].ndx_old = 0;
        dbuff[i].ndx_new = 0;
        dbuff[i].num_zero = 0;
 
        need_update[i] = 0;
        for (a = 1; a < MAX_ATTS; a++) {
            if (isosurf->att[a].changed) {
                read = 0;
                /* changed to map attribute */
                if (isosurf->att[a].att_src == MAP_ATT) {
                    vf = gvl_file_get_volfile(isosurf->att[a].hfile);
                    read = 1;
                }
                /* changed threshold value */
                if (a == ATT_TOPO) {
                    isosurf->att[a].hfile = gvol->hfile;
                    vf = gvl_file_get_volfile(gvol->hfile);
                    read = 1;
                }
                /* initialize reading in selected mode */
                if (read) {
                    gvl_file_set_mode(vf, 3);
                    gvl_file_start_read(vf);
                }
 
                /* set update flag - isosurface will be calc */
                if (read || IS_IN_DATA(a)) {
                    need_update[i] = 1;
                    need_update_global = 1;
                }
            }
        }
 
        if (need_update[i]) {
            /* set data buffer */
            dbuff[i].old = isosurf->data;
        }
    }
 
    /* calculate if only some isosurface changed */
    if (need_update_global) {
 
        ResX = gvol->isosurf_x_mod;
        ResY = gvol->isosurf_y_mod;
        ResZ = gvol->isosurf_z_mod;
 
        Cols = gvol->cols / ResX;
        Rows = gvol->rows / ResY;
        Depths = gvol->depths / ResZ;
 
        /* calc isosurface - marching cubes - start */
 
        for (z = 0; z < Depths - 1; z++) {
            for (y = 0; y < Rows - 1; y++) {
                for (x = 0; x < Cols - 1; x++) {
                    for (i = 0; i < gvol->n_isosurfs; i++) {
                        /* recalculate only changed isosurfaces */
                        if (need_update[i]) {
                            iso_calc_cube(gvol->isosurf[i], x, y, z, &dbuff[i]);
                        }
                    }
                }
            }
        }
    }
    /* end */
 
    /* deinitialize */
    for (i = 0; i < gvol->n_isosurfs; i++) {
        isosurf = gvol->isosurf[i];
 
        /* set new isosurface data */
        if (need_update[i]) {
            if (dbuff[i].num_zero != 0)
                gvl_write_char(dbuff[i].ndx_new++, &(dbuff[i].new),
                               dbuff[i].num_zero);
 
            if (dbuff[i].old == isosurf->data)
                dbuff[i].old = NULL;
            G_free(isosurf->data);
            gvl_align_data(dbuff[i].ndx_new, &(dbuff[i].new));
            isosurf->data = dbuff[i].new;
            isosurf->data_desc = 0;
        }
 
        for (a = 1; a < MAX_ATTS; a++) {
            if (isosurf->att[a].changed) {
                read = 0;
                /* changed map attribute */
                if (isosurf->att[a].att_src == MAP_ATT) {
                    vf = gvl_file_get_volfile(isosurf->att[a].hfile);
                    read = 1;
                }
                /* changed threshold value */
                if (a == ATT_TOPO) {
                    isosurf->att[a].hfile = gvol->hfile;
                    vf = gvl_file_get_volfile(gvol->hfile);
                    read = 1;
                }
                /* deinitialize reading */
                if (read) {
                    gvl_file_end_read(vf);
 
                    /* set data description */
                    SET_IN_DATA(a);
                }
                isosurf->att[a].changed = 0;
            }
            else if (isosurf->att[a].att_src == MAP_ATT) {
                /* set data description */
                SET_IN_DATA(a);
            }
        }
    }
 
    G_free(dbuff);
    G_free(need_update);
 
    return (1);
}
 
/*!
   \brief ADD
 
   \param pos
   \param data
   \param c
 */
void gvl_write_char(int pos, unsigned char **data, unsigned char c)
{
    /* check to need allocation memory */
    if ((pos % BUFFER_SIZE) == 0) {
        *data = G_realloc(*data, sizeof(char) * ((pos / BUFFER_SIZE) + 1) *
                                     BUFFER_SIZE);
        if (!(*data)) {
            return;
        }
 
        G_debug(3,
                "gvl_write_char(): reallocate memory for pos : %d to : %lu B",
                pos, sizeof(char) * ((pos / BUFFER_SIZE) + 1) * BUFFER_SIZE);
    }
 
    (*data)[pos] = c;
}
 
/*!
   \brief Read char
 
   \param pos position index
   \param data data buffer
 
   \return char on success
   \return NULL on failure
 */
unsigned char gvl_read_char(int pos, const unsigned char *data)
{
    if (!data)
        return '\0';
 
    return data[pos];
}
 
/*!
   \brief Append data to buffer
 
   \param pos position index
   \param data data buffer
 */
void gvl_align_data(int pos, unsigned char **data)
{
    if (pos <= 0) {
        if (*data) {
            G_free(*data);
            *data = NULL;
        }
        return;
    }
    /* realloc memory to fit in data length */
    unsigned char *p;
    p = (unsigned char *)G_realloc(*data, sizeof(unsigned char) *
                                              pos); /* G_fatal_error */
    if (!p) {
        return;
    }
 
    G_debug(3, "gvl_align_data(): reallocate memory finally to : %d B", pos);
 
    *data = p;
 
    return;
}
 
/************************************************************************/
/* SLICES */
 
/************************************************************************/
 
#define DISTANCE_2(x1, y1, x2, y2) \
    sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2))
 
#define SLICE_MODE_INTERP_NO  0
#define SLICE_MODE_INTERP_YES 1
 
/*!
   \brief Get volume value
 
   \param gvl pointer to geovol struct
   \param x,y,z
 
   \return value
 */
float slice_get_value(geovol *gvl, int x, int y, int z)
{
    static double d;
    static geovol_file *vf;
    static int type;
    static float value;
 
    if (x < 0 || y < 0 || z < 0 || (x > gvl->cols - 1) || (y > gvl->rows - 1) ||
        (z > gvl->depths - 1))
        return 0.;
 
    /* get volume file from attribute handle */
    vf = gvl_file_get_volfile(gvl->hfile);
    type = gvl_file_get_data_type(vf);
 
    /* get value from volume file */
    if (type == VOL_DTYPE_FLOAT) {
        gvl_file_get_value(vf, x, y, z, &value);
    }
    else if (type == VOL_DTYPE_DOUBLE) {
        gvl_file_get_value(vf, x, y, z, &d);
        value = (float)d;
    }
    else {
        return 0.;
    }
 
    return value;
}
 
/*!
   \brief Calculate slices
 
   \param gvl pointer to geovol struct
   \param ndx_slc
   \param colors
 
   \return 1
 */
int slice_calc(geovol *gvl, int ndx_slc, void *colors)
{
    int cols, rows, c, r;
    int i, j, k, pos, color;
    int *p_x, *p_y, *p_z;
    float *p_ex, *p_ey, *p_ez;
    float value, v[8];
    float x, y, z, ei, ej, ek, stepx, stepy, stepz;
    float f_cols, f_rows, distxy, distz, modxy, modx, mody, modz;
 
    geovol_slice *slice;
    geovol_file *vf;
 
    slice = gvl->slice[ndx_slc];
 
    /* set mods, pointer to x, y, z step value */
    if (slice->dir == X) {
        modx = ResY;
        mody = ResZ;
        modz = ResX;
        p_x = &k;
        p_y = &i;
        p_z = &j;
        p_ex = &ek;
        p_ey = &ei;
        p_ez = &ej;
    }
    else if (slice->dir == Y) {
        modx = ResX;
        mody = ResZ;
        modz = ResY;
        p_x = &i;
        p_y = &k;
        p_z = &j;
        p_ex = &ei;
        p_ey = &ek;
        p_ez = &ej;
    }
    else {
        modx = ResX;
        mody = ResY;
        modz = ResZ;
        p_x = &i;
        p_y = &j;
        p_z = &k;
        p_ex = &ei;
        p_ey = &ej;
        p_ez = &ek;
    }
 
    /* distance between slice def. points */
    distxy = DISTANCE_2(slice->x2, slice->y2, slice->x1, slice->y1);
    distz = fabsf(slice->z2 - slice->z1);
 
    /* distance between slice def points is zero - nothing to do */
    if (distxy == 0. || distz == 0.) {
        return (1);
    }
 
    /* start reading volume file */
    vf = gvl_file_get_volfile(gvl->hfile);
    gvl_file_set_mode(vf, 3);
    gvl_file_start_read(vf);
 
    /* set xy resolution */
    modxy = DISTANCE_2((slice->x2 - slice->x1) / distxy * modx,
                       (slice->y2 - slice->y1) / distxy * mody, 0., 0.);
 
    /* cols/rows of slice */
    f_cols = distxy / modxy;
    cols = f_cols > (int)f_cols ? (int)f_cols + 1 : (int)f_cols;
 
    f_rows = distz / modz;
    rows = f_rows > (int)f_rows ? (int)f_rows + 1 : (int)f_rows;
 
    /* set x,y initially to first slice point */
    x = slice->x1;
    y = slice->y1;
 
    /* set x,y step */
    stepx = (slice->x2 - slice->x1) / f_cols;
    stepy = (slice->y2 - slice->y1) / f_cols;
    stepz = (slice->z2 - slice->z1) / f_rows;
 
    /* set position in slice data */
    pos = 0;
 
    /* loop in slice cols */
    for (c = 0; c < cols + 1; c++) {
 
        /* convert x, y to integer - index in grid */
        i = (int)x;
        j = (int)y;
 
        /* distance between index and real position */
        ei = x - (float)i;
        ej = y - (float)j;
 
        /* set z to slice z1 point */
        z = slice->z1;
 
        /* loop in slice rows */
        for (r = 0; r < rows + 1; r++) {
 
            /* distance between index and real position */
            k = (int)z;
            ek = z - (float)k;
 
            /* get interpolated value */
            if (slice->mode == SLICE_MODE_INTERP_YES) {
                /* get grid values */
                v[0] = slice_get_value(gvl, *p_x, *p_y, *p_z);
                v[1] = slice_get_value(gvl, *p_x + 1, *p_y, *p_z);
                v[2] = slice_get_value(gvl, *p_x, *p_y + 1, *p_z);
                v[3] = slice_get_value(gvl, *p_x + 1, *p_y + 1, *p_z);
 
                v[4] = slice_get_value(gvl, *p_x, *p_y, *p_z + 1);
                v[5] = slice_get_value(gvl, *p_x + 1, *p_y, *p_z + 1);
                v[6] = slice_get_value(gvl, *p_x, *p_y + 1, *p_z + 1);
                v[7] = slice_get_value(gvl, *p_x + 1, *p_y + 1, *p_z + 1);
 
                /* get interpolated value */
                value = v[0] * (1. - *p_ex) * (1. - *p_ey) * (1. - *p_ez) +
                        v[1] * (*p_ex) * (1. - *p_ey) * (1. - *p_ez) +
                        v[2] * (1. - *p_ex) * (*p_ey) * (1. - *p_ez) +
                        v[3] * (*p_ex) * (*p_ey) * (1. - *p_ez) +
                        v[4] * (1. - *p_ex) * (1. - *p_ey) * (*p_ez) +
                        v[5] * (*p_ex) * (1. - *p_ey) * (*p_ez) +
                        v[6] * (1. - *p_ex) * (*p_ey) * (*p_ez) +
                        v[7] * (*p_ex) * (*p_ey) * (*p_ez);
 
                /* no interp value */
            }
            else {
                value = slice_get_value(gvl, *p_x, *p_y, *p_z);
            }
 
            /* translate value to color */
            color = Gvl_get_color_for_value(colors, &value);
 
            /* write color to slice data */
            gvl_write_char(pos++, &(slice->data), color & RED_MASK);
            gvl_write_char(pos++, &(slice->data), (color & GRN_MASK) >> 8);
            gvl_write_char(pos++, &(slice->data), (color & BLU_MASK) >> 16);
 
            /* step in z */
            if (r + 1 > f_rows) {
                z += stepz * (f_rows - (float)r);
            }
            else {
                z += stepz;
            }
        }
 
        /* step in x,y */
        if (c + 1 > f_cols) {
            x += stepx * (f_cols - (float)c);
            y += stepy * (f_cols - (float)c);
        }
        else {
            x += stepx;
            y += stepy;
        }
    }
 
    /* end reading volume file */
    gvl_file_end_read(vf);
    gvl_align_data(pos, &(slice->data));
 
    return (1);
}
 
/*!
   \brief Calculate slices for given volume set
 
   \param gvol pointer to geovol struct
 
   \return 1
 */
int gvl_slices_calc(geovol *gvol)
{
    int i;
    void *colors;
 
    G_debug(5, "gvl_slices_calc(): id=%d", gvol->gvol_id);
 
    /* set current resolution */
    ResX = gvol->slice_x_mod;
    ResY = gvol->slice_y_mod;
    ResZ = gvol->slice_z_mod;
 
    /* set current num of cols, rows, depths */
    Cols = gvol->cols / ResX;
    Rows = gvol->rows / ResY;
    Depths = gvol->depths / ResZ;
 
    /* load colors for geovol file */
    Gvl_load_colors_data(&colors, gvl_file_get_name(gvol->hfile));
 
    /* calc changed slices */
    for (i = 0; i < gvol->n_slices; i++) {
        if (gvol->slice[i]->changed) {
            slice_calc(gvol, i, colors);
 
            /* set changed flag */
            gvol->slice[i]->changed = 0;
        }
    }
 
    /* free color */
    Gvl_unload_colors_data(colors);
 
    return (1);
}