/* volume.c */
/*
* Vis5D system for visualizing five dimensional gridded data sets.
* Copyright (C) 1990 - 2000 Bill Hibbard, Johan Kellum, Brian Paul,
* Dave Santek, and Andre Battaiola.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* As a special exception to the terms of the GNU General Public
* License, you are permitted to link Vis5D with (and distribute the
* resulting source and executables) the LUI library (copyright by
* Stellar Computer Inc. and licensed for distribution with Vis5D),
* the McIDAS library, and/or the NetCDF library, where those
* libraries are governed by the terms of their own licenses.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
#include "../config.h"
/*
* Volume rendering. Requires alpha blending capability. This is
* supported in on SGI hardware with VGX or higher graphics. Otherwise,
* it's done in software on some other systems.
*
* Volume rendering is also supported on everthing running OpenGL!
*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "globals.h"
#include "graphics.h"
#include "grid.h"
#include "memory.h"
#include "proj.h"
#include "volume.h"
#if HAVE_OPENGL
# include <GL/gl.h>
#elif HAVE_SGI_GL
# include <gl/gl.h>
#endif
#define ABS(A) ( (A) < 0 ? -(A) : (A) )
/* Direction of slices: */
#define BOTTOM_TO_TOP 0
#define TOP_TO_BOTTOM 1
#define EAST_TO_WEST 2
#define WEST_TO_EAST 3
#define NORTH_TO_SOUTH 4
#define SOUTH_TO_NORTH 5
#define LUT_SIZE 255
/*
* Allocate a new volume structure and return a pointer to it. If we
* can't do volume rendering return NULL.
* Input: nr, nc, nl - dimensions of largest volume we'll encounter.
* Return: pointer to volume struct or NULL
*/
struct volume *alloc_volume( Context ctx, int nr, int nc, int nl )
{
int volren = 0;
struct volume *v = NULL;
if (ctx->dpy_ctx->Projection == PROJ_CYLINDRICAL ||
ctx->dpy_ctx->Projection == PROJ_SPHERICAL){
ctx->dpy_ctx->VolRender = 0;
return 0;
}
if (nl <= 1){
ctx->dpy_ctx->VolRender = 0;
return 0; /* no volume variables */
}
#if defined (HAVE_SGI_GL) || defined (DENALI)
alphabits = getgdesc( GD_BLEND );
if (alphabits==0) {
ctx->dpy_ctx->VolRender = 0;
/* no alpha support means no volume rendering */
return NULL;
}
volren = 1;
#endif
#if defined(HAVE_OPENGL)
volren = 1;
#endif
/* MJK 12.15.98 */
#ifdef HAVE_PEX
volren = 1;
#endif
if (volren) {
v = (struct volume *) malloc( sizeof(struct volume) );
/* initialize the volume struct */
v->valid = 0;
/* No matter which way we slice it, we need the same size arrays for */
/* storing the vertices and colors. */
v->vertex = (float *) allocate( ctx, nl*nr*nc*3*sizeof(float) );
v->index = (uint_1 *) allocate( ctx, nl*nr*nc*sizeof(uint_1) );
if (!v->vertex || !v->index) {
printf("WARNING: insufficient memory for volume rendering (%d bytes needed)\n",
nl * nr * nc * ((int) (3*sizeof(float)+sizeof(uint_1))) );
ctx->dpy_ctx->VolRender = 0;
return NULL;
}
v->oldnr = nr;
v->oldnc = nc;
v->oldnl = nl;
/* MJK 12.15.98 */
#ifdef HAVE_PEX
v->dir = -1;
#endif
}
if (v){
ctx->dpy_ctx->VolRender = 1;
}
else{
ctx->dpy_ctx->VolRender = 0;
}
return v;
}
void free_volume( Context ctx)
{
deallocate( ctx, ctx->Volume->vertex,
ctx->Volume->oldnl*ctx->Volume->oldnr*ctx->Volume->oldnc*3*sizeof(float));
deallocate( ctx, ctx->Volume->index,
ctx->Volume->oldnl*ctx->Volume->oldnr*ctx->Volume->oldnc*sizeof(uint_1));
free(ctx->Volume);
ctx->Volume = NULL;
}
/*
* Compute a volume rendering of the given grid.
* Input: data - 3-D data grid.
* time, var - time step and variable
* nr, nc, nl - size of 3-D grid.
* min, max - min and max data value in data grid.
* dir - direction to do slicing
* v - pointer to a volume struct with the vertex and index
* fields pointing to sufficiently large buffers.
* Output: v - volume struct describing the computed volume.
*/
static int compute_volume( Context ctx, float data[],
int time, int var,
int nr, int nc, int nl, int lowlev,
float min, float max,
int dir,
struct volume *v )
{
Display_Context dtx;
float zs[MAXLEVELS];
register int ir, ic, il, i, j;
register float x, y, dx, dy;
register float dscale, val;
register int ilnrnc, icnr; /* to help speed up array indexing */
dtx = ctx->dpy_ctx;
/* compute graphics Z for each grid level */
for (il=0; il<nl; il++) {
zs[il] = gridlevel_to_z(ctx, time, var, (float) (il + lowlev));
}
/* compute some useful values */
dx = (dtx->Xmax-dtx->Xmin) / (nc-1);
dy = (dtx->Ymax-dtx->Ymin) / (nr-1);
dscale = (float) (LUT_SIZE-1) / (max-min);
v->dir = dir;
switch (dir) {
case BOTTOM_TO_TOP:
/* compute slices from bottom to top */
v->slices = nl;
v->rows = nr;
v->cols = nc;
i = 0; /* index into vertex array */
j = 0; /* index into index array */
for (il=0; il<nl; il++) {
y = dtx->Ymax;
ilnrnc = il * nr * nc;
for (ir=0;ir<nr;ir++) {
x = dtx->Xmin;
for (ic=0;ic<nc;ic++) {
/* compute vertex */
v->vertex[i++] = x;
v->vertex[i++] = y;
v->vertex[i++] = zs[il];
/* compute color table index */
val = data[ ilnrnc + ic * nr + ir ];
if (IS_MISSING(val) ||
val < min || val > max)
v->index[j++] = 255;
else
v->index[j++] = (uint_1) (int) ((val-min) * dscale);
x += dx;
}
y -= dy;
}
}
break;
case TOP_TO_BOTTOM:
/* compute slices from top to bottom */
v->slices = nl;
v->rows = nr;
v->cols = nc;
i = 0; /* index into vertex array */
j = 0; /* index into index array */
for (il=nl-1; il>=0; il--) {
y = dtx->Ymax;
ilnrnc = il * nr * nc;
for (ir=0;ir<nr;ir++) {
x = dtx->Xmin;
for (ic=0;ic<nc;ic++) {
/* compute vertex */
v->vertex[i++] = x;
v->vertex[i++] = y;
v->vertex[i++] = zs[il];
/* compute color table index */
val = data[ ilnrnc + ic * nr + ir ];
if (IS_MISSING(val) ||
val < min || val > max)
v->index[j++] = 255;
else
v->index[j++] = (uint_1) (int) ((val-min) * dscale);
x += dx;
}
y -= dy;
}
}
break;
case WEST_TO_EAST:
/* compute slices from west to east */
v->slices = nc;
v->rows = nl;
v->cols = nr;
i = 0; /* index into vertex array */
j = 0; /* index into index array */
x = dtx->Xmin;
for (ic=0; ic<nc; ic++) {
icnr = ic * nr;
for (il=nl-1;il>=0;il--) {
y = dtx->Ymin;
ilnrnc = il * nr * nc + icnr;
for (ir=nr-1;ir>=0;ir--) {
/* compute vertex */
v->vertex[i++] = x;
v->vertex[i++] = y;
v->vertex[i++] = zs[il];
/* compute color table index */
val = data[ ilnrnc + ir ];
if (IS_MISSING(val) ||
val < min || val > max)
v->index[j++] = 255;
else
v->index[j++] = (uint_1) (int) ((val-min) * dscale);
y += dy;
}
}
x += dx;
}
break;
case EAST_TO_WEST:
/* compute slices from east to west */
v->slices = nc;
v->rows = nl;
v->cols = nr;
i = 0; /* index into vertex array */
j = 0; /* index into index array */
x = dtx->Xmax;
for (ic=nc-1; ic>=0; ic--) {
icnr = ic*nr;
for (il=nl-1;il>=0;il--) {
y = dtx->Ymin;
ilnrnc = il * nr * nc + icnr;
for (ir=nr-1;ir>=0;ir--) {
/* compute vertex */
v->vertex[i++] = x;
v->vertex[i++] = y;
v->vertex[i++] = zs[il];
/* compute color table index */
val = data[ ilnrnc + ir ];
if (IS_MISSING(val) ||
val < min || val > max)
v->index[j++] = 255;
else
v->index[j++] = (uint_1) (int) ((val-min) * dscale);
y += dy;
}
}
x -= dx;
}
break;
case NORTH_TO_SOUTH:
/* compute slices from north to south */
v->slices = nr;
v->rows = nl;
v->cols = nc;
i = 0; /* index into vertex array */
j = 0; /* index into index array */
y = dtx->Ymax;
for (ir=0; ir<nr; ir++) {
for (il=nl-1;il>=0;il--) {
x = dtx->Xmin;
ilnrnc = il * nr * nc;
for (ic=0;ic<nc;ic++) {
/* compute vertex */
v->vertex[i++] = x;
v->vertex[i++] = y;
v->vertex[i++] = zs[il];
/* compute color table index */
val = data[ ilnrnc + ic*nr + ir ];
if (IS_MISSING(val) ||
val < min || val > max)
v->index[j++] = 255;
else
v->index[j++] = (uint_1) (int) ((val-min) * dscale);
x += dx;
}
}
y -= dy;
}
break;
case SOUTH_TO_NORTH:
/* compute slices from south to north */
v->slices = nr;
v->rows = nl;
v->cols = nc;
i = 0; /* index into vertex array */
j = 0; /* index into index array */
y = dtx->Ymin;
for (ir=nr-1; ir>=0; ir--) {
for (il=nl-1;il>=0;il--) {
x = dtx->Xmin;
ilnrnc = il * nr * nc;
for (ic=0;ic<nc;ic++) {
/* compute vertex */
v->vertex[i++] = x;
v->vertex[i++] = y;
v->vertex[i++] = zs[il];
/* compute color table index */
val = data[ ilnrnc + ic*nr + ir ];
if (IS_MISSING(val) ||
val < min || val > max)
v->index[j++] = 255;
else
v->index[j++] = (uint_1) (int) ((val-min) * dscale);
x += dx;
}
}
y += dy;
}
break;
default:
printf("Error in compute_volume!\n");
} /* switch */
v->valid = 1;
return 1;
}
static int compute_volumePRIME( Context ctx, float data[],
int time, int var,
int nr, int nc, int nl, int lowlev,
float min, float max,
int dir,
struct volume *v )
{
Display_Context dtx;
float zs[MAXLEVELS];
register int ir, ic, il, i, j;
register float x, y, dx, dy;
register float dscale, val;
register int ilnrnc, icnr; /* to help speed up array indexing */
float s1,s2,s3,s4,s5,s6,s7,s8;
float grow, gcol, glev;
float row, col, lev;
int gr0,gr1,gc0,gc1,gl0,gl1;
float ger, gec, gel;
dtx = ctx->dpy_ctx;
/* compute graphics Z for each grid level */
for (il=0; il<nl; il++) {
zs[il] = gridlevelPRIME_to_zPRIME(dtx, time, var, (float) (il + lowlev));
}
/* compute some useful values */
dx = (dtx->Xmax-dtx->Xmin) / (nc-1);
dy = (dtx->Ymax-dtx->Ymin) / (nr-1);
dscale = (float) (LUT_SIZE-1) / (max-min);
v->dir = dir;
switch (dir) {
case BOTTOM_TO_TOP:
/* compute slices from bottom to top */
v->slices = nl;
v->rows = nr;
v->cols = nc;
i = 0; /* index into vertex array */
j = 0; /* index into index array */
for (il=0; il<nl; il++) {
y = dtx->Ymax;
ilnrnc = il * nr * nc;
for (ir=0;ir<nr;ir++) {
x = dtx->Xmin;
for (ic=0;ic<nc;ic++) {
/* compute vertex */
v->vertex[i++] = x;
v->vertex[i++] = y;
v->vertex[i++] = zs[il];
/* compute color table index */
row = ir;
col = ic;
lev = il;
gridPRIME_to_grid( ctx, time, var, 1, &row, &col, &lev,
&grow, &gcol, &glev);
if ( grow < 0 || gcol < 0 || glev < 0 ||
grow >= ctx->Nr || gcol >= ctx->Nc || glev >= ctx->Nl[var]){
v->index[j++] = 255;
}
else{
gr0= (int) grow;
gr1= gr0+1;
if (gr0 == ctx->Nr-1){
gr1 = gr0;
}
gc0 = (int) gcol;
gc1 = gc0 + 1;
if (gc0 == ctx->Nc-1){
gc1 = gc0;
}
gl0 = (int) glev;
gl1 = gl0+1;
if (gl0 == ctx->Nl[var]-1){
gl1 = gl0;
}
ger = grow - (float) gr0; /* in [0,1) */
gec = gcol - (float) gc0; /* in [0,1) */
gel = glev - (float) gl0; /* in [0,1) */
if (ger==0.0){
gr1 = gr0;
}
if (gec==0.0){
gc1 = gc0;
}
if (gel==0.0){
gl1 = gl0;
}
s1 = data[(gl0*ctx->Nc+gc0)*ctx->Nr+gr0];
s2 = data[(gl0*ctx->Nc+gc0)*ctx->Nr+gr1];
s3 = data[(gl0*ctx->Nc+gc1)*ctx->Nr+gr0];
s4 = data[(gl0*ctx->Nc+gc1)*ctx->Nr+gr1];
s5 = data[(gl1*ctx->Nc+gc0)*ctx->Nr+gr0];
s6 = data[(gl1*ctx->Nc+gc0)*ctx->Nr+gr1];
s7 = data[(gl1*ctx->Nc+gc1)*ctx->Nr+gr0];
s8 = data[(gl1*ctx->Nc+gc1)*ctx->Nr+gr1];
if (IS_MISSING(s1) || IS_MISSING(s2) ||
IS_MISSING(s3) || IS_MISSING(s4) ||
IS_MISSING(s5) || IS_MISSING(s6) ||
IS_MISSING(s7) || IS_MISSING(s8)){
v->index[j++] = 255;
}
else{
val = ( s1 * (1.0-ger) * (1.0-gec)
+ s2 * ger * (1.0-gec)
+ s3 * (1.0-ger) * gec
+ s4 * ger * gec ) * (1.0-gel)
+ ( s5 * (1.0-ger) * (1.0-gec)
+ s6 * ger * (1.0-gec)
+ s7 * (1.0-ger) * gec
+ s8 * ger * gec ) * gel;
if (val < min || val > max){
v->index[j++] = 255;
}
else{
v->index[j++] = (uint_1) (int) ((val-min) * dscale);
}
}
}
x += dx;
}
y -= dy;
}
}
break;
case TOP_TO_BOTTOM:
/* compute slices from top to bottom */
v->slices = nl;
v->rows = nr;
v->cols = nc;
i = 0; /* index into vertex array */
j = 0; /* index into index array */
for (il=nl-1; il>=0; il--) {
y = dtx->Ymax;
ilnrnc = il * nr * nc;
for (ir=0;ir<nr;ir++) {
x = dtx->Xmin;
for (ic=0;ic<nc;ic++) {
/* compute vertex */
v->vertex[i++] = x;
v->vertex[i++] = y;
v->vertex[i++] = zs[il];
/* compute color table index */
row = ir;
col = ic;
lev = il;
gridPRIME_to_grid( ctx, time, var, 1, &row, &col, &lev,
&grow, &gcol, &glev);
if ( grow < 0 || gcol < 0 || glev < 0 ||
grow >= ctx->Nr || gcol >= ctx->Nc || glev >= ctx->Nl[var]){
v->index[j++] = 255;
}
else{
gr0= (int) grow;
gr1= gr0+1;
if (gr0 == ctx->Nr-1){
gr1 = gr0;
}
gc0 = (int) gcol;
gc1 = gc0 + 1;
if (gc0 == ctx->Nc-1){
gc1 = gc0;
}
gl0 = (int) glev;
gl1 = gl0+1;
if (gl0 == ctx->Nl[var]-1){
gl1 = gl0;
}
ger = grow - (float) gr0; /* in [0,1) */
gec = gcol - (float) gc0; /* in [0,1) */
gel = glev - (float) gl0; /* in [0,1) */
if (ger==0.0){
gr1 = gr0;
}
if (gec==0.0){
gc1 = gc0;
}
if (gel==0.0){
gl1 = gl0;
}
s1 = data[(gl0*ctx->Nc+gc0)*ctx->Nr+gr0];
s2 = data[(gl0*ctx->Nc+gc0)*ctx->Nr+gr1];
s3 = data[(gl0*ctx->Nc+gc1)*ctx->Nr+gr0];
s4 = data[(gl0*ctx->Nc+gc1)*ctx->Nr+gr1];
s5 = data[(gl1*ctx->Nc+gc0)*ctx->Nr+gr0];
s6 = data[(gl1*ctx->Nc+gc0)*ctx->Nr+gr1];
s7 = data[(gl1*ctx->Nc+gc1)*ctx->Nr+gr0];
s8 = data[(gl1*ctx->Nc+gc1)*ctx->Nr+gr1];
if (IS_MISSING(s1) || IS_MISSING(s2) ||
IS_MISSING(s3) || IS_MISSING(s4) ||
IS_MISSING(s5) || IS_MISSING(s6) ||
IS_MISSING(s7) || IS_MISSING(s8)){
v->index[j++] = 255;
}
else{
val = ( s1 * (1.0-ger) * (1.0-gec)
+ s2 * ger * (1.0-gec)
+ s3 * (1.0-ger) * gec
+ s4 * ger * gec ) * (1.0-gel)
+ ( s5 * (1.0-ger) * (1.0-gec)
+ s6 * ger * (1.0-gec)
+ s7 * (1.0-ger) * gec
+ s8 * ger * gec ) * gel;
if (val < min || val > max){
v->index[j++] = 255;
}
else{
v->index[j++] = (uint_1) (int) ((val-min) * dscale);
}
}
}
x += dx;
}
y -= dy;
}
}
break;
case WEST_TO_EAST:
/* compute slices from west to east */
v->slices = nc;
v->rows = nl;
v->cols = nr;
i = 0; /* index into vertex array */
j = 0; /* index into index array */
x = dtx->Xmin;
for (ic=0; ic<nc; ic++) {
icnr = ic * nr;
for (il=nl-1;il>=0;il--) {
y = dtx->Ymin;
ilnrnc = il * nr * nc + icnr;
for (ir=nr-1;ir>=0;ir--) {
/* compute vertex */
v->vertex[i++] = x;
v->vertex[i++] = y;
v->vertex[i++] = zs[il];
/* compute color table index */
row = ir;
col = ic;
lev = il;
gridPRIME_to_grid( ctx, time, var, 1, &row, &col, &lev,
&grow, &gcol, &glev);
if ( grow < 0 || gcol < 0 || glev < 0 ||
grow >= ctx->Nr || gcol >= ctx->Nc || glev >= ctx->Nl[var]){
v->index[j++] = 255;
}
else{
gr0= (int) grow;
gr1= gr0+1;
if (gr0 == ctx->Nr-1){
gr1 = gr0;
}
gc0 = (int) gcol;
gc1 = gc0 + 1;
if (gc0 == ctx->Nc-1){
gc1 = gc0;
}
gl0 = (int) glev;
gl1 = gl0+1;
if (gl0 == ctx->Nl[var]-1){
gl1 = gl0;
}
ger = grow - (float) gr0; /* in [0,1) */
gec = gcol - (float) gc0; /* in [0,1) */
gel = glev - (float) gl0; /* in [0,1) */
if (ger==0.0){
gr1 = gr0;
}
if (gec==0.0){
gc1 = gc0;
}
if (gel==0.0){
gl1 = gl0;
}
s1 = data[(gl0*ctx->Nc+gc0)*ctx->Nr+gr0];
s2 = data[(gl0*ctx->Nc+gc0)*ctx->Nr+gr1];
s3 = data[(gl0*ctx->Nc+gc1)*ctx->Nr+gr0];
s4 = data[(gl0*ctx->Nc+gc1)*ctx->Nr+gr1];
s5 = data[(gl1*ctx->Nc+gc0)*ctx->Nr+gr0];
s6 = data[(gl1*ctx->Nc+gc0)*ctx->Nr+gr1];
s7 = data[(gl1*ctx->Nc+gc1)*ctx->Nr+gr0];
s8 = data[(gl1*ctx->Nc+gc1)*ctx->Nr+gr1];
if (IS_MISSING(s1) || IS_MISSING(s2) ||
IS_MISSING(s3) || IS_MISSING(s4) ||
IS_MISSING(s5) || IS_MISSING(s6) ||
IS_MISSING(s7) || IS_MISSING(s8)){
v->index[j++] = 255;
}
else{
val = ( s1 * (1.0-ger) * (1.0-gec)
+ s2 * ger * (1.0-gec)
+ s3 * (1.0-ger) * gec
+ s4 * ger * gec ) * (1.0-gel)
+ ( s5 * (1.0-ger) * (1.0-gec)
+ s6 * ger * (1.0-gec)
+ s7 * (1.0-ger) * gec
+ s8 * ger * gec ) * gel;
if (val < min || val > max){
v->index[j++] = 255;
}
else{
v->index[j++] = (uint_1) (int) ((val-min) * dscale);
}
}
}
y += dy;
}
}
x += dx;
}
break;
case EAST_TO_WEST:
/* compute slices from east to west */
v->slices = nc;
v->rows = nl;
v->cols = nr;
i = 0; /* index into vertex array */
j = 0; /* index into index array */
x = dtx->Xmax;
for (ic=nc-1; ic>=0; ic--) {
icnr = ic*nr;
for (il=nl-1;il>=0;il--) {
y = dtx->Ymin;
ilnrnc = il * nr * nc + icnr;
for (ir=nr-1;ir>=0;ir--) {
/* compute vertex */
v->vertex[i++] = x;
v->vertex[i++] = y;
v->vertex[i++] = zs[il];
/* compute color table index */
row = ir;
col = ic;
lev = il;
gridPRIME_to_grid( ctx, time, var, 1, &row, &col, &lev,
&grow, &gcol, &glev);
if ( grow < 0 || gcol < 0 || glev < 0 ||
grow >= ctx->Nr || gcol >= ctx->Nc || glev >= ctx->Nl[var]){
v->index[j++] = 255;
}
else{
gr0= (int) grow;
gr1= gr0+1;
if (gr0 == ctx->Nr-1){
gr1 = gr0;
}
gc0 = (int) gcol;
gc1 = gc0 + 1;
if (gc0 == ctx->Nc-1){
gc1 = gc0;
}
gl0 = (int) glev;
gl1 = gl0+1;
if (gl0 == ctx->Nl[var]-1){
gl1 = gl0;
}
ger = grow - (float) gr0; /* in [0,1) */
gec = gcol - (float) gc0; /* in [0,1) */
gel = glev - (float) gl0; /* in [0,1) */
if (ger==0.0){
gr1 = gr0;
}
if (gec==0.0){
gc1 = gc0;
}
if (gel==0.0){
gl1 = gl0;
}
s1 = data[(gl0*ctx->Nc+gc0)*ctx->Nr+gr0];
s2 = data[(gl0*ctx->Nc+gc0)*ctx->Nr+gr1];
s3 = data[(gl0*ctx->Nc+gc1)*ctx->Nr+gr0];
s4 = data[(gl0*ctx->Nc+gc1)*ctx->Nr+gr1];
s5 = data[(gl1*ctx->Nc+gc0)*ctx->Nr+gr0];
s6 = data[(gl1*ctx->Nc+gc0)*ctx->Nr+gr1];
s7 = data[(gl1*ctx->Nc+gc1)*ctx->Nr+gr0];
s8 = data[(gl1*ctx->Nc+gc1)*ctx->Nr+gr1];
if (IS_MISSING(s1) || IS_MISSING(s2) ||
IS_MISSING(s3) || IS_MISSING(s4) ||
IS_MISSING(s5) || IS_MISSING(s6) ||
IS_MISSING(s7) || IS_MISSING(s8)){
v->index[j++] = 255;
}
else{
val = ( s1 * (1.0-ger) * (1.0-gec)
+ s2 * ger * (1.0-gec)
+ s3 * (1.0-ger) * gec
+ s4 * ger * gec ) * (1.0-gel)
+ ( s5 * (1.0-ger) * (1.0-gec)
+ s6 * ger * (1.0-gec)
+ s7 * (1.0-ger) * gec
+ s8 * ger * gec ) * gel;
if (val < min || val > max){
v->index[j++] = 255;
}
else{
v->index[j++] = (uint_1) (int) ((val-min) * dscale);
}
}
}
y += dy;
}
}
x -= dx;
}
break;
case NORTH_TO_SOUTH:
/* compute slices from north to south */
v->slices = nr;
v->rows = nl;
v->cols = nc;
i = 0; /* index into vertex array */
j = 0; /* index into index array */
y = dtx->Ymax;
for (ir=0; ir<nr; ir++) {
for (il=nl-1;il>=0;il--) {
x = dtx->Xmin;
ilnrnc = il * nr * nc;
for (ic=0;ic<nc;ic++) {
/* compute vertex */
v->vertex[i++] = x;
v->vertex[i++] = y;
v->vertex[i++] = zs[il];
/* compute color table index */
row = ir;
col = ic;
lev = il;
gridPRIME_to_grid( ctx, time, var, 1, &row, &col, &lev,
&grow, &gcol, &glev);
if ( grow < 0 || gcol < 0 || glev < 0 ||
grow >= ctx->Nr || gcol >= ctx->Nc || glev >= ctx->Nl[var]){
v->index[j++] = 255;
}
else{
gr0= (int) grow;
gr1= gr0+1;
if (gr0 == ctx->Nr-1){
gr1 = gr0;
}
gc0 = (int) gcol;
gc1 = gc0 + 1;
if (gc0 == ctx->Nc-1){
gc1 = gc0;
}
gl0 = (int) glev;
gl1 = gl0+1;
if (gl0 == ctx->Nl[var]-1){
gl1 = gl0;
}
ger = grow - (float) gr0; /* in [0,1) */
gec = gcol - (float) gc0; /* in [0,1) */
gel = glev - (float) gl0; /* in [0,1) */
if (ger==0.0){
gr1 = gr0;
}
if (gec==0.0){
gc1 = gc0;
}
if (gel==0.0){
gl1 = gl0;
}
s1 = data[(gl0*ctx->Nc+gc0)*ctx->Nr+gr0];
s2 = data[(gl0*ctx->Nc+gc0)*ctx->Nr+gr1];
s3 = data[(gl0*ctx->Nc+gc1)*ctx->Nr+gr0];
s4 = data[(gl0*ctx->Nc+gc1)*ctx->Nr+gr1];
s5 = data[(gl1*ctx->Nc+gc0)*ctx->Nr+gr0];
s6 = data[(gl1*ctx->Nc+gc0)*ctx->Nr+gr1];
s7 = data[(gl1*ctx->Nc+gc1)*ctx->Nr+gr0];
s8 = data[(gl1*ctx->Nc+gc1)*ctx->Nr+gr1];
if (IS_MISSING(s1) || IS_MISSING(s2) ||
IS_MISSING(s3) || IS_MISSING(s4) ||
IS_MISSING(s5) || IS_MISSING(s6) ||
IS_MISSING(s7) || IS_MISSING(s8)){
v->index[j++] = 255;
}
else{
val = ( s1 * (1.0-ger) * (1.0-gec)
+ s2 * ger * (1.0-gec)
+ s3 * (1.0-ger) * gec
+ s4 * ger * gec ) * (1.0-gel)
+ ( s5 * (1.0-ger) * (1.0-gec)
+ s6 * ger * (1.0-gec)
+ s7 * (1.0-ger) * gec
+ s8 * ger * gec ) * gel;
if (val < min || val > max){
v->index[j++] = 255;
}
else{
v->index[j++] = (uint_1) (int) ((val-min) * dscale);
}
}
}
x += dx;
}
}
y -= dy;
}
break;
case SOUTH_TO_NORTH:
/* compute slices from south to north */
v->slices = nr;
v->rows = nl;
v->cols = nc;
i = 0; /* index into vertex array */
j = 0; /* index into index array */
y = dtx->Ymin;
for (ir=nr-1; ir>=0; ir--) {
for (il=nl-1;il>=0;il--) {
x = dtx->Xmin;
ilnrnc = il * nr * nc;
for (ic=0;ic<nc;ic++) {
/* compute vertex */
v->vertex[i++] = x;
v->vertex[i++] = y;
v->vertex[i++] = zs[il];
/* compute color table index */
row = ir;
col = ic;
lev = il;
gridPRIME_to_grid( ctx, time, var, 1, &row, &col, &lev,
&grow, &gcol, &glev);
if ( grow < 0 || gcol < 0 || glev < 0 ||
grow >= ctx->Nr || gcol >= ctx->Nc || glev >= ctx->Nl[var]){
v->index[j++] = 255;
}
else{
gr0= (int) grow;
gr1= gr0+1;
if (gr0 == ctx->Nr-1){
gr1 = gr0;
}
gc0 = (int) gcol;
gc1 = gc0 + 1;
if (gc0 == ctx->Nc-1){
gc1 = gc0;
}
gl0 = (int) glev;
gl1 = gl0+1;
if (gl0 == ctx->Nl[var]-1){
gl1 = gl0;
}
ger = grow - (float) gr0; /* in [0,1) */
gec = gcol - (float) gc0; /* in [0,1) */
gel = glev - (float) gl0; /* in [0,1) */
if (ger==0.0){
gr1 = gr0;
}
if (gec==0.0){
gc1 = gc0;
}
if (gel==0.0){
gl1 = gl0;
}
s1 = data[(gl0*ctx->Nc+gc0)*ctx->Nr+gr0];
s2 = data[(gl0*ctx->Nc+gc0)*ctx->Nr+gr1];
s3 = data[(gl0*ctx->Nc+gc1)*ctx->Nr+gr0];
s4 = data[(gl0*ctx->Nc+gc1)*ctx->Nr+gr1];
s5 = data[(gl1*ctx->Nc+gc0)*ctx->Nr+gr0];
s6 = data[(gl1*ctx->Nc+gc0)*ctx->Nr+gr1];
s7 = data[(gl1*ctx->Nc+gc1)*ctx->Nr+gr0];
s8 = data[(gl1*ctx->Nc+gc1)*ctx->Nr+gr1];
if (IS_MISSING(s1) || IS_MISSING(s2) ||
IS_MISSING(s3) || IS_MISSING(s4) ||
IS_MISSING(s5) || IS_MISSING(s6) ||
IS_MISSING(s7) || IS_MISSING(s8)){
v->index[j++] = 255;
}
else{
val = ( s1 * (1.0-ger) * (1.0-gec)
+ s2 * ger * (1.0-gec)
+ s3 * (1.0-ger) * gec
+ s4 * ger * gec ) * (1.0-gel)
+ ( s5 * (1.0-ger) * (1.0-gec)
+ s6 * ger * (1.0-gec)
+ s7 * (1.0-ger) * gec
+ s8 * ger * gec ) * gel;
if (val < min || val > max){
v->index[j++] = 255;
}
else{
v->index[j++] = (uint_1) (int) ((val-min) * dscale);
}
}
}
x += dx;
}
}
y += dy;
}
break;
default:
printf("Error in compute_volumePRIME!\n");
} /* switch */
v->valid = 1;
return 1;
}
/*
* Render the volume described by the given volume struct.
* Return: 1 = ok
* 0 = bad volume struct.
*/
static int render_volume( Context ctx,
struct volume *v, unsigned int ctable[] )
{
register int rows, cols, slices, i, j, s;
register int rows1, cols1;
register uint_1 *cp0, *cp1;
register float *vp0, *vp1;
int fastdraw;
int stride = 1;
if (!v || !v->slices)
return 0;
#if defined (HAVE_SGI_GL) || defined (DENALI)
lmcolor( LMC_COLOR ); /* no shading */
blendfunction( BF_SA, BF_MSA ); /* enable alpha blending */
#endif
#ifdef HAVE_OPENGL
glBlendFunc( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA );
glEnable( GL_BLEND );
check_gl_error( "render_volume (glBlendFunc)" );
#endif
/* put rows, cols, slices values into registers */
rows = v->rows-1; /* number of quad strips = number of data rows - 1 */
cols = v->cols;
slices = v->slices;
/* setup color and vertex pointers */
cp0 = v->index;
cp1 = cp0 + cols;
vp0 = v->vertex;
vp1 = vp0 + cols * 3; /* 3 floats per vertex */
/* MJK 12.15.98 */
#ifdef HAVE_PEX
rows++;
j = rows * cols;
for (s = 0; s < slices; s++)
{
draw_volume_quadmesh (rows, cols, vp0, cp0, ctable);
vp0 += j * 3;
cp0 += j;
}
return 1;
#endif
vis5d_check_fastdraw(ctx->dpy_ctx->dpy_context_index, &fastdraw);
if (fastdraw) {
stride = ctx->dpy_ctx->VStride;
}
/* sanity check */
if(stride<=0)
stride = 1;
/*
** adjust rows and cols based on stride. N.B. appears to be one more
** row than we actually use
*/
rows1 = (rows + 1 - 1) / stride;
cols1 = ((cols - 1) / stride) + 1;
/* loop over slices */
for (s=0;s<slices;s+=stride) {
cp0 = v->index + (s * rows * cols) +
(s * cols); /* skip a row after each slice */
vp0 = v->vertex + (s * rows * cols * 3) +
(s * cols * 3); /* skip a row after each slice */
cp1 = cp0 + (cols * stride);
vp1 = vp0 + (cols * stride * 3); /* 3 floats per vertex */
/* draw 'rows' quadrilateral strips */
for (i=0;i<rows1;i++) {
#if defined(SGI_GL) || defined(DENALI)
bgnqstrip();
for (j=0;j<cols1;j++) {
cpack( ctable[cp0[i*stride*cols+j*stride]] );
v3f( &vp0[i*stride*cols+j*stride] );
cpack( ctable[cp1[i*stride*cols+j*stride]] );
v3f( &vp1[i*stride*cols+j*stride] );
}
endqstrip();
#endif
#ifdef HAVE_OPENGL
glBegin( GL_QUAD_STRIP );
for (j=0;j<cols1;j++) {
glColor4ubv( (GLubyte *) &ctable[cp0[i*stride*cols+j*stride]] );
glVertex3fv( &vp0[(i*stride*cols+j*stride)*3] );
glColor4ubv( (GLubyte *) &ctable[cp1[i*stride*cols+j*stride]] );
glVertex3fv( &vp1[(i*stride*cols+j*stride)*3] );
}
glEnd();
#endif
}
}
#if defined(HAVE_SGI_GL) || defined(DENALI)
blendfunction( BF_ONE, BF_ZERO ); /* disable alpha blending */
#endif
#ifdef HAVE_OPENGL
glDisable( GL_BLEND );
check_gl_error( "render_volume (glDisable)" );
#endif
return 1;
}
/* MJK 12.15.98 */
#ifdef HAVE_PEX
void draw_volume( Context ctx, int it, int ip, unsigned int *ctable )
{
float *data;
static int prev_it[VIS5D_MAX_CONTEXTS];
static int prev_ip[VIS5D_MAX_CONTEXTS];
static int do_once = 1;
int dir;
float x, y, z, ax, ay, az;
float xyz[3], xy[3][2], xy0[2];
Display_Context dtx;
dtx = ctx->dpy_ctx;
if (do_once){
int yo;
for (yo=0; yo<VIS5D_MAX_CONTEXTS; yo++){
prev_it[yo] = -1;
prev_ip[yo] = -1;
}
do_once = 0;
}
xyz[0] = xyz[1] = xyz[2] = 0.0;
project (xyz, &xy0[0], &xy0[1]);
xyz[0] = 1.0, xyz[1] = xyz[2] = 0.0;
project (xyz, &xy[0][0], &xy[0][1]);
xy[0][0] -= xy0[0], xy[0][1] -= xy0[1];
xyz[1] = 1.0, xyz[0] = xyz[2] = 0.0;
project (xyz, &xy[1][0], &xy[1][1]);
xy[1][0] -= xy0[0], xy[1][1] -= xy0[1];
xyz[2] = 1.0, xyz[0] = xyz[1] = 0.0;
project (xyz, &xy[2][0], &xy[2][1]);
xy[2][0] -= xy0[0], xy[2][1] -= xy0[1];
ax = (xy[0][0] * xy[0][0]) + (xy[0][1] * xy[0][1]);
ay = (xy[1][0] * xy[1][0]) + (xy[1][1] * xy[1][1]);
az = (xy[2][0] * xy[2][0]) + (xy[2][1] * xy[2][1]);
if ((ax <= ay) && (ax <= az))
{
x = (xy[1][1] * xy[2][0]) - (xy[1][0] * xy[2][1]);
dir = (x > 0.0) ? WEST_TO_EAST : EAST_TO_WEST;
}
else if ((ay <= ax) && (ay <= az))
{
x = (xy[2][1] * xy[0][0]) - (xy[2][0] * xy[0][1]);
dir = (x > 0.0) ? SOUTH_TO_NORTH : NORTH_TO_SOUTH;
}
else
{
x = (xy[0][1] * xy[1][0]) - (xy[0][0] * xy[1][1]);
dir = (x > 0.0) ? BOTTOM_TO_TOP : TOP_TO_BOTTOM;
}
/* If this is a new time step or variable then invalidate old volumes */
if (it!=prev_it[ctx->context_index] || ip!=prev_ip[ctx->context_index]) {
ctx->Volume->valid = 0;
prev_it[ctx->context_index] = it;
prev_ip[ctx->context_index] = ip;
}
/* Determine if we have to compute a set of slices for the direction. */
if (ctx->Volume->dir!=dir || ctx->Volume->valid==0) {
data = get_grid( ctx, it, ip );
if (data) {
if (ctx->GridSameAsGridPRIME){
compute_volume( ctx, data, it, ip, ctx->Nr, ctx->Nc, ctx->Nl[ip],
ctx->Variable[ip]->LowLev, ctx->Variable[ip]->MinVal, ctx->Variable[ip]->MaxVal,
dir, ctx->Volume );
}
else{
compute_volumePRIME( ctx, data, it, ip, dtx->Nr, dtx->Nc, dtx->Nl,
dtx->LowLev, ctx->Variable[ip]->MinVal, ctx->Variable[ip]->MaxVal,
dir, ctx->Volume );
}
release_grid( ctx, it, ip, data );
}
}
render_volume( ctx, ctx->Volume, ctable );
}
#else
/*
* Draw a volumetric rendering of the grid for timestep it and variable ip.
* Input: it - timestep
* ip - variable
*/
void draw_volume( Context ctx, int it, int ip, unsigned int *ctable )
{
float *data;
static int prev_it[VIS5D_MAX_CONTEXTS];
static int prev_ip[VIS5D_MAX_CONTEXTS];
static int do_once = 1;
int dir;
float x, y, z, ax, ay, az;
MATRIX ctm, proj;
Display_Context dtx;
dtx = ctx->dpy_ctx;
if (do_once){
int yo;
for (yo=0; yo<VIS5D_MAX_CONTEXTS; yo++){
prev_it[yo] = -1;
prev_ip[yo] = -1;
}
do_once = 0;
}
/* Get 3rd column values from transformation matrix */
#if defined (HAVE_SGI_GL) || defined (DENALI)
/* Compute orientation of 3-D box with respect to current matrices */
/* with no assumptions about the location of the camera. This was */
/* done for the CAVE. */
mmode( MPROJECTION );
getmatrix( proj );
mmode( MVIEWING );
getmatrix( ctm );
#endif
#if defined(HAVE_OPENGL)
glGetFloatv( GL_PROJECTION_MATRIX, (GLfloat *) proj );
glGetFloatv( GL_MODELVIEW_MATRIX, (GLfloat *) ctm );
check_gl_error( "draw_volume" );
#endif
/* compute third column values in the product of ctm*proj */
x = ctm[0][0]*proj[0][2] + ctm[0][1]*proj[1][2]
+ ctm[0][2]*proj[2][2] + ctm[0][3]*proj[3][2];
y = ctm[1][0]*proj[0][2] + ctm[1][1]*proj[1][2]
+ ctm[1][2]*proj[2][2] + ctm[1][3]*proj[3][2];
z = ctm[2][0]*proj[0][2] + ctm[2][1]*proj[1][2]
+ ctm[2][2]*proj[2][2] + ctm[2][3]*proj[3][2];
/* examine values to determine how to draw slices */
ax = ABS(x);
ay = ABS(y);
az = ABS(z);
if (ax>=ay && ax>=az) {
/* draw x-axis slices */
dir = (x<0.0) ? WEST_TO_EAST : EAST_TO_WEST;
}
else if (ay>=ax && ay>=az) {
/* draw y-axis slices */
dir = (y<0.0) ? SOUTH_TO_NORTH : NORTH_TO_SOUTH;
}
else {
/* draw z-axis slices */
dir = (z<0.0) ? BOTTOM_TO_TOP : TOP_TO_BOTTOM;
}
/* If this is a new time step or variable then invalidate old volumes */
if (it!=prev_it[ctx->context_index] || ip!=prev_ip[ctx->context_index]) {
ctx->Volume->valid = 0;
prev_it[ctx->context_index] = it;
prev_ip[ctx->context_index] = ip;
}
/* Determine if we have to compute a set of slices for the direction. */
if (ctx->Volume->dir!=dir || ctx->Volume->valid==0) {
data = get_grid( ctx, it, ip );
if (data) {
if (ctx->GridSameAsGridPRIME){
compute_volume( ctx, data, it, ip, ctx->Nr, ctx->Nc, ctx->Nl[ip],
ctx->Variable[ip]->LowLev, ctx->Variable[ip]->MinVal, ctx->Variable[ip]->MaxVal,
dir, ctx->Volume );
}
else{
compute_volumePRIME( ctx, data, it, ip, dtx->Nr, dtx->Nc, dtx->Nl,
dtx->LowLev, ctx->Variable[ip]->MinVal, ctx->Variable[ip]->MaxVal,
dir, ctx->Volume );
}
release_grid( ctx, it, ip, data );
}
}
render_volume( ctx, ctx->Volume, ctable );
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1