source: trunk/packages/vizservers/nanovis/dxReader.cpp @ 1429

 
/*
 * ----------------------------------------------------------------------
 * Nanovis: Visualization of Nanoelectronics Data
 *
 *  dxReader.cpp 
 *      This module contains openDX readers for 2D and 3D volumes.
 *
 * ======================================================================
 *  AUTHOR:  Wei Qiao <qiaow@purdue.edu>
 *           Michael McLennan <mmclennan@purdue.edu>
 *           Purdue Rendering and Perceptualization Lab (PURPL)
 *
 *  Copyright (c) 2004-2006  Purdue Research Foundation
 *
 *  See the file "license.terms" for information on usage and
 *  redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES.
 * ======================================================================
 */

// common dx functions
#include "dxReaderCommon.h"

#include <stdio.h>
#include <math.h>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <sys/types.h>
#include <unistd.h>

#include "Nv.h"

#include "nanovis.h"
#include "RpField1D.h"
#include "RpFieldRect3D.h"
#include "RpFieldPrism3D.h"
#include <Unirect.h>

//transfer function headers
#include "ZincBlendeVolume.h"
#include "NvZincBlendeReconstructor.h"

#define  _LOCAL_ZINC_TEST_ 0

static INLINE char *
skipspaces(char *string) 
{
    while (isspace(*string)) {
        string++;
    }
    return string;
}

static INLINE char *
getline(char **stringPtr, char *endPtr) 
{
    char *line, *p;

    line = skipspaces(*stringPtr);
    for (p = line; p < endPtr; p++) {
        if (*p == '\n') {
            *p++ = '\0';
            *stringPtr = p;
            return line;
        }
    }
    *stringPtr = p;
    return line;
}

Rappture::Unirect3d *
ReadDxVectorFieldData(Rappture::Outcome &result, size_t length, char *string) 
{
    int nx, ny, nz, npts;
    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
    char *p, *endPtr;


    dx = dy = dz = 0.0;         // Suppress compiler warning.
    x0 = y0 = z0 = 0.0;         // May not have an origin line.
    for (p = string, endPtr = p + length; p < endPtr; /*empty*/) {
        char *line;

        line = getline(&p, endPtr);
        if (line == endPtr) {
            break;
        }
        if ((line[0] == '#') || (line == '\0')) {
            continue;           // Skip blank or comment lines.
        }
        if (sscanf(line, "object %*d class gridpositions counts %d %d %d", 
                   &nx, &ny, &nz) == 3) {
            printf("w:%d h:%d d:%d\n", nx, ny, nz);
            // found grid size
        } else if (sscanf(line, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
            // found origin
        } else if (sscanf(line, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
            // found one of the delta lines
            if (ddx != 0.0) {
                dx = ddx;
            } else if (ddy != 0.0) {
                dy = ddy;
            } else if (ddz != 0.0) {
                dz = ddz;
            }
        } else if (sscanf(line, "object %*d class array type %*s shape 3"
                " rank 1 items %d data follows", &npts) == 1) {
            printf("point %d\n", npts);
            if (npts != nx*ny*nz) {
                result.addError("inconsistent data: expected %d points"
                                " but found %d points", nx*ny*nz, npts);
                return NULL;
            }
            break;
        } else if (sscanf(line, "object %*d class array type %*s rank 0"
                " times %d data follows", &npts) == 1) {
            if (npts != nx*ny*nz) {
                result.addError("inconsistent data: expected %d points"
                                " but found %d points", nx*ny*nz, npts);
                return NULL;
            }
            break;
        }
    }
    if (npts != nx*ny*nz) {
        result.addError("inconsistent data: expected %d points"
                        " but found %d points", nx*ny*nz, npts);
        return NULL;
    }
    float *values = new float[npts];
    int nValues = 0;
    for (int ix=0; ix < nx; ix++) {
        for (int iy=0; iy < ny; iy++) {
            for (int iz=0; iz < nz; iz++) {
                char *line;
                if ((p == endPtr) || (nValues > npts)) {
                    break;
                }
                line = getline(&p, endPtr);
                if ((line[0] == '#') || (line[0] == '\0')) {
                    continue;   // Skip blank or comment lines.
                }
                double vx, vy, vz;
                if (sscanf(line, "%lg %lg %lg", &vx, &vy, &vz) == 3) {
                    int nindex = (iz*nx*ny + iy*nx + ix) * 3;
                    values[nindex] = vx;
                    values[nindex+1] = vy;
                    values[nindex+2] = vz;
                    nValues++;
                }
            }
        }
    }
    // make sure that we read all of the expected points
    if (nValues != npts) {
        result.addError("inconsistent data: expected %d points"
                        " but found %d points", npts, nValues);
        delete values;
        return NULL;
    }
    return new Rappture::Unirect3d(x0, x0 + dx * nx, nx, y0, y0 + dy * ny, ny,
                                   z0, z0 + dz * nz, nz, nValues, values, 3);
}


#ifdef notdef
/*
 * Load a 3D vector field from a dx-format file
 */
bool
load_vector_stream2(Rappture::Outcome &result, int ivol, size_t length,
                    char *string)
{
    Rappture::Unirect3d *dataPtr;
    dataPtr = ReadDxVectorFieldData(result, length, string);
    if (dataPtr == NULL) {
        return false;
    }
    Vector3 physicalMin;
    Vector3 physicalMax;

    physicalMin.set(dataPtr->xMin(), dataPtr->yMin(), dataPtr->zMin());
    physicalMax.set(dataPtr->xMax(), dataPtr->yMax(), dataPtr->zMax());

    double max_mag = -DBL_MAX, min_mag = DBL_MAX;
    float *values = dataPtr->values();
    for (size_t ix=0; ix < dataPtr->xNum(); ix++) {
        for (size_t iy=0; iy < dataPtr->yNum(); iy++) {
            for (size_t iz=0; iz < dataPtr->zNum(); iz++) {
                double vx, vy, vz, vm;
                vx = values[0];
                vy = values[1];
                vz = values[2];
                vm = sqrt(vx*vx + vy*vy + vz*vz);
                if (vm > max_mag) {
                    max_mag = vm;
                }
                if (vm < min_mag) {
                    min_mag = vm;
                }
                values += 3;
            }
        }
    }
    
    float *data = new float[4*dataPtr->nValues()];
    memset(data, 0, sizeof(float) * 4 * dataPtr->nValues());
    fprintf(stderr, "generating %dx%dx%d = %d points\n", 
            dataPtr->xNum(), dataPtr->yNum(), dataPtr->zNum(), 
            dataPtr->nValues());

    float *destPtr = data;
    float *srcPtr = dataPtr->values();
    for (size_t iz=0; iz < dataPtr->zNum(); iz++) {
        for (size_t iy=0; iy < dataPtr->yNum(); iy++) {
            for (size_t ix=0; ix < dataPtr->xNum(); ix++) {
                double vx, vy, vz, vm;
                vx = values[0];
                vy = values[1];
                vz = values[2];
                vm = sqrt(vx*vx + vy*vy + vz*vz);
                destPtr[0] = vm / max_mag; 
                destPtr[1] = vx /(2.0*max_mag) + 0.5; 
                destPtr[2] = vy /(2.0*max_mag) + 0.5; 
                destPtr[3] = vz /(2.0*max_mag) + 0.5; 
                srcPtr += 3;
                destPtr += 4;
            }
        }
    }
    
    Volume *volPtr;
    volPtr = NanoVis::load_volume(ivol, dataPtr->xNum(), dataPtr->yNum(), 
                dataPtr->zNum(), 4, data, min_mag, max_mag, 0);
    volPtr->xAxis.SetRange(dataPtr->xMin(), dataPtr->xMax());
    volPtr->yAxis.SetRange(dataPtr->yMin(), dataPtr->yMax());
    volPtr->zAxis.SetRange(dataPtr->zMin(), dataPtr->zMax());
    volPtr->wAxis.SetRange(min_mag, max_mag);
    volPtr->update_pending = true;
    volPtr->setPhysicalBBox(physicalMin, physicalMax);
    delete dataPtr;
    delete data;
    return true;
}

bool
load_vector_stream(Rappture::Outcome result, int index, size_t length,
                    char *string)
{
    Rappture::Unirect3d *dataPtr;
    dataPtr = ReadDxVectorFieldData(result, length, string);
    if (dataPtr == NULL) {
        return false;
    }
    Vector3 physicalMin;
    Vector3 physicalMax;

    physicalMin.set(dataPtr->xMin(), dataPtr->yMin(), dataPtr->zMin());
    physicalMax.set(dataPtr->xMax(), dataPtr->yMax(), dataPtr->zMax());

    Rappture::Mesh1D xgrid(dataPtr->xMin(), dataPtr->yMin(), dataPtr->xNum());
    Rappture::Mesh1D ygrid(dataPtr->yMin(), dataPtr->yMax(), dataPtr->yNum());
    Rappture::Mesh1D zgrid(dataPtr->zMin(), dataPtr->zMax(), dataPtr->zNum());
    Rappture::FieldRect3D xfield(xgrid, ygrid, zgrid);
    Rappture::FieldRect3D yfield(xgrid, ygrid, zgrid);
    Rappture::FieldRect3D zfield(xgrid, ygrid, zgrid);
    
    float *values = dataPtr->values();
    size_t npts = 0;
    for (size_t ix=0; ix < dataPtr->xNum(); ix++) {
        for (size_t iy=0; iy < dataPtr->yNum(); iy++) {
            for (size_t iz=0; iz < dataPtr->zNum(); iz++) {
                xfield.define(npts, values[0]);
                yfield.define(npts, values[1]);
                zfield.define(npts, values[2]);
                npts++;
                values += 3;
            }
        }
    }

    double dx, dy, dz;
    // figure out a good mesh spacing
    int nsample = 30;
    dx = xfield.rangeMax(Rappture::xaxis) - xfield.rangeMin(Rappture::xaxis);
    dy = xfield.rangeMax(Rappture::yaxis) - xfield.rangeMin(Rappture::yaxis);
    dz = xfield.rangeMax(Rappture::zaxis) - xfield.rangeMin(Rappture::zaxis);
    double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);

    printf("dx:%lf dy:%lf dz:%lf dmin:%lf\n", dx, dy, dz, dmin);

    size_t nx, ny, nz;
    nx = (int)ceil(dx/dmin);
    ny = (int)ceil(dy/dmin);
    nz = (int)ceil(dz/dmin);
    
#ifndef NV40
    // must be an even power of 2 for older cards
    nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
    ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
    nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
#endif
    
    float *data = new float[4*nx*ny*nz];
    memset(data, 0, sizeof(float) * 4 * nx * ny * nz);
    
    std::cout << "generating " << nx << "x" << ny << "x" << nz << " = " << nx*ny*nz << " points" << std::endl;
    
    // generate the uniformly sampled data that we need for a volume
    double vmin = 1e21;
    double vmax = -1e21;
    double nzero_min = 0.0;
    size_t ngen = 0;
    for (size_t iz=0; iz < nz; iz++) {
        double zval = dataPtr->zMin() + iz*dmin;
        for (size_t iy=0; iy < ny; iy++) {
            double yval = dataPtr->yMin() + iy*dmin;
            for (size_t ix=0; ix < nx; ix++) {
                double xval = dataPtr->xMin() + ix*dmin;
                double vx, vy, vz;

                vx = xfield.value(xval,yval,zval);
                vy = yfield.value(xval,yval,zval);
                vz = zfield.value(xval,yval,zval);
                
                double vm;
                vm = sqrt(vx*vx + vy*vy + vz*vz);
                if (vm < vmin) {
                    vmin = vm;
                } else if (vm > vmax) {
                    vmax = vm;
                }
                if ((vm != 0.0f) && (vm < nzero_min)) {
                    nzero_min = vm;
                }
                data[ngen++] = vm;
                data[ngen++] = vx;
                data[ngen++] = vy;
                data[ngen++] = vz;
            }
        }
    }
    
    ngen = 0;
    // scale should be accounted.
    for (ngen=0; ngen < npts; /*empty*/) {
        data[ngen] = data[ngen] / vmax; ++ngen;
        data[ngen] = (data[ngen]/(2.0*vmax) + 0.5); ++ngen;
        data[ngen] = (data[ngen]/(2.0*vmax) + 0.5); ++ngen;
        data[ngen] = (data[ngen]/(2.0*vmax) + 0.5); ++ngen;
    }
    Volume *volPtr;
    volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data, vmin, vmax,
                                  nzero_min);
    
    volPtr->xAxis.SetRange(dataPtr->xMin(), dataPtr->xMin() + (nx * dx));
    volPtr->yAxis.SetRange(dataPtr->yMin(), dataPtr->yMin() + (ny * dy));
    volPtr->zAxis.SetRange(dataPtr->zMin(), dataPtr->zMin() + (nz * dz));
    volPtr->wAxis.SetRange(vmin, vmax);
    volPtr->update_pending = true;
    delete [] data;
    return true;
}
#endif

/* Load a 3D volume from a dx-format file
 */
bool
load_volume_stream2(Rappture::Outcome &result, int index, std::iostream& fin)
{
    printf("load_volume_stream2\n");
    Rappture::MeshTri2D xymesh;
    int dummy, nx, ny, nz, nxy, npts;
    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
    char line[128], type[128], *start;

    int isrect = 1;
    dx = dy = dz = 0.0;         // Suppress compiler warning.
    x0 = y0 = z0 = 0.0;         // May not have an origin line.
    do {
        fin.getline(line,sizeof(line)-1);
        for (start=&line[0]; *start == ' ' || *start == '\t'; start++)
            ;  // skip leading blanks

        if (*start != '#') {  // skip comment lines
            if (sscanf(start, "object %d class gridpositions counts %d %d %d", 
                       &dummy, &nx, &ny, &nz) == 4) {
                // found grid size
                isrect = 1;
            } else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
                isrect = 0;
                double xx, yy, zz;
                for (int i=0; i < nxy; i++) {
                    fin.getline(line,sizeof(line)-1);
                    if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
                        xymesh.addNode( Rappture::Node2D(xx,yy) );
                    }
                }

                char fpts[128];
                sprintf(fpts, "/tmp/tmppts%d", getpid());
                char fcells[128];
                sprintf(fcells, "/tmp/tmpcells%d", getpid());

                std::ofstream ftmp(fpts);
                // save corners of bounding box first, to work around meshing
                // problems in voronoi utility
                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
                for (int i=0; i < nxy; i++) {
                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
                }
                ftmp.close();

                char cmdstr[512];
                sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
                if (system(cmdstr) == 0) {
                    int cx, cy, cz;
                    std::ifstream ftri(fcells);
                    while (!ftri.eof()) {
                        ftri.getline(line,sizeof(line)-1);
                        if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
                            if (cx >= 4 && cy >= 4 && cz >= 4) {
                                // skip first 4 boundary points
                                xymesh.addCell(cx-4, cy-4, cz-4);
                            }
                        }
                    }
                    ftri.close();
                } else {
                    result.error("triangularization failed");
                    return false;
                }
                unlink(fpts), unlink(fcells);
            } else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
                // found z-grid
            } else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
                // found origin
            } else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
        int count = 0;
                // found one of the delta lines
                if (ddx != 0.0) {
            dx = ddx;
            count++;
        }
        if (ddy != 0.0) {
            dy = ddy;
            count++;
        }
        if (ddz != 0.0) {
            dz = ddz;
            count++;
        }
        if (count > 1) {
            result.addError("don't know how to handle multiple non-zero"
                            " delta values");
            return false;
        }
            } else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
                if (isrect && (npts != nx*ny*nz)) {
                    result.addError("inconsistent data: expected %d points "
                                    " but found %d points", nx*ny*nz, npts);
                    return false;
                } else if (!isrect && (npts != nxy*nz)) {
                    result.addError("inconsistent data: expected %d points "
                                    " but found %d points", nxy*nz, npts);
                    return false;
                }
                break;
            } else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
                if (npts != nx*ny*nz) {
                    result.addError("inconsistent data: expected %d points "
                                    " but found %d points", nx*ny*nz, npts);
                    return false;
                }
                break;
            }
        }
    } while (!fin.eof());

    fprintf(stderr, "found nx=%d ny=%d, nz=%d, x0=%f, y0=%f, z0=%f\n",
            nx, ny, nz, x0, y0, z0);
    // read data points
    if (fin.eof()) {
        result.addError("EOF found: expecting %d points", npts);
        return false;
    }
    if (isrect) {
        double dval[6];
        int nread = 0;
        int ix = 0;
        int iy = 0;
        int iz = 0;
        float* data = new float[nx *  ny *  nz * 4];
        memset(data, 0, nx*ny*nz*4);
        double vmin = 1e21;
        double nzero_min = 1e21;
        double vmax = -1e21;
        
        
        while (!fin.eof() && nread < npts) {
            fin.getline(line,sizeof(line)-1);
            int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
            
            for (int p=0; p < n; p++) {
                int nindex = (iz*nx*ny + iy*nx + ix) * 4;
                data[nindex] = dval[p];
                
                if (dval[p] < vmin) {
                    vmin = dval[p];
                } else if (dval[p] > vmax) {
                    vmax = dval[p];
                }
                if (dval[p] != 0.0f && dval[p] < nzero_min) {
                    nzero_min = dval[p];
                }
                
                nread++;
                if (++iz >= nz) {
                    iz = 0;
                    if (++iy >= ny) {
                        iy = 0;
                        ++ix;
                    }
                }
            }
        }
        
        // make sure that we read all of the expected points
        if (nread != nx*ny*nz) {
            result.addError("inconsistent data: expected %d points "
                            " but found %d points", nx*ny*nz, nread);
            return false;
        }
        
        double dv = vmax - vmin;
        int count = nx*ny*nz;
        int ngen = 0;
        double v;
        if (dv == 0.0) {
            dv = 1.0;
        }
        
        for (int i = 0; i < count; ++i) {
            v = data[ngen];
            // scale all values [0-1], -1 => out of bounds
            //
            // INSOO
            v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
            data[ngen] = v;
            ngen += 4;
        }
        
        computeSimpleGradient(data, nx, ny, nz);
        
        dx = nx;
        dy = ny;
        dz = nz;
        
        Volume *volPtr;
        volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
                                      vmin, vmax, nzero_min);
        volPtr->xAxis.SetRange(x0, x0 + (nx * dx));
        volPtr->yAxis.SetRange(y0, y0 + (ny * dy));
        volPtr->zAxis.SetRange(z0, z0 + (nz * dz)); 
        volPtr->wAxis.SetRange(vmin, vmax);
        volPtr->update_pending = true;
        delete [] data;
        
    } else {
        Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
        Rappture::FieldPrism3D field(xymesh, zgrid);
        
        double dval;
        int nread = 0;
        int ixy = 0;
        int iz = 0;
        while (!fin.eof() && nread < npts) {
            fin >> dval;
            if (fin.fail()) {
                result.addError("after %d of %d points: can't read number", 
                                nread, npts);
                return false;
            } else {
                int nid = nxy*iz + ixy;
                field.define(nid, dval);
                
                nread++;
                if (++iz >= nz) {
                    iz = 0;
                    ixy++;
                }
            }
        }
        
        // make sure that we read all of the expected points
        if (nread != nxy*nz) {
            result.addError("inconsistent data: expected %d points "
                            "but found %d points", nxy*nz, nread);
            return false;
        }
        
        // figure out a good mesh spacing
        int nsample = 30;
        x0 = field.rangeMin(Rappture::xaxis);
        dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
        y0 = field.rangeMin(Rappture::yaxis);
        dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
        z0 = field.rangeMin(Rappture::zaxis);
        dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
        double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
        
        nx = (int)ceil(dx/dmin);
        ny = (int)ceil(dy/dmin);
        nz = (int)ceil(dz/dmin);
#ifndef NV40
        // must be an even power of 2 for older cards
        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
#endif
        float *data = new float[4*nx*ny*nz];
        
        double vmin = field.valueMin();
        double dv = field.valueMax() - field.valueMin();
        if (dv == 0.0) {
            dv = 1.0; 
        }
        // generate the uniformly sampled data that we need for a volume
        int ngen = 0;
        double nzero_min = 0.0;
        for (iz=0; iz < nz; iz++) {
            double zval = z0 + iz*dmin;
            for (int iy=0; iy < ny; iy++) {
                double yval = y0 + iy*dmin;
                for (int ix=0; ix < nx; ix++) {
                    double xval = x0 + ix*dmin;
                    double v = field.value(xval,yval,zval);
                    
                    if (v != 0.0f && v < nzero_min) {
                        nzero_min = v;
                    }
                    // scale all values [0-1], -1 => out of bounds
                    v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
                    data[ngen] = v;
                    
                    ngen += 4;
                }
            }
        }
        
        // FIXME: This next section of code should be replaced by a
        // call to the computeSimpleGradient() function. There is a slight
        // difference in the code below and the aforementioned function
        // in that the commented out lines in the else statements are
        // different.
        //
        // Compute the gradient of this data.  BE CAREFUL: center
        // calculation on each node to avoid skew in either direction.
        ngen = 0;
        for (int iz=0; iz < nz; iz++) {
            for (int iy=0; iy < ny; iy++) {
                for (int ix=0; ix < nx; ix++) {
                    // gradient in x-direction
                    //double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
                    //double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
                    double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
                    double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
                    if (valm1 < 0 || valp1 < 0) {
                        data[ngen+1] = 0.0;
                    } else {
                        data[ngen+1] = valp1-valm1; // assume dx=1
                        //data[ngen+1] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
                    }
                    
                    // gradient in y-direction
                    valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
                    valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
                    if (valm1 < 0 || valp1 < 0) {
                        data[ngen+2] = 0.0;
                    } else {
                        data[ngen+2] = valp1-valm1; // assume dy=1
                        //data[ngen+2] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
                    }
                    
                    // gradient in z-direction
                    valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
                    valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
                    if (valm1 < 0 || valp1 < 0) {
                        data[ngen+3] = 0.0;
                    } else {
                        data[ngen+3] = valp1-valm1; // assume dz=1
                        //data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
                    }
                    
                    ngen += 4;
                }
            }
        }
        
        Volume *volPtr;
        volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
                                      field.valueMin(), field.valueMax(), nzero_min);
        volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
                               field.rangeMax(Rappture::xaxis));
        volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
                               field.rangeMax(Rappture::yaxis));
        volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
                               field.rangeMax(Rappture::zaxis));
        volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
        volPtr->update_pending = true;
        delete [] data;
    }
    //
    // Center this new volume on the origin.
    //
    float dx0 = -0.5;
    float dy0 = -0.5*dy/dx;
    float dz0 = -0.5*dz/dx;
    NanoVis::volume[index]->move(Vector3(dx0, dy0, dz0));
    return true;
}

bool
load_volume_stream(Rappture::Outcome &result, int index, std::iostream& fin)
{
    printf("load_volume_stream\n");

    Rappture::MeshTri2D xymesh;
    int dummy, nx, ny, nz, nxy, npts;
    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
    char line[128], type[128], *start;

    int isrect = 1;

    dx = dy = dz = 0.0;         // Suppress compiler warning.
    x0 = y0 = z0 = 0.0;         // May not have an origin line.
    while (!fin.eof()) {
        fin.getline(line, sizeof(line) - 1);
        if (fin.fail()) {
            result.error("error in data stream");
            return false;
        }
        for (start=line; *start == ' ' || *start == '\t'; start++)
            ;  // skip leading blanks

        if (*start != '#') {  // skip comment lines
            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
                // found grid size
                isrect = 1;
            } else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
                isrect = 0;

                double xx, yy, zz;
                for (int i=0; i < nxy; i++) {
                    fin.getline(line,sizeof(line)-1);
                    if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
                        xymesh.addNode( Rappture::Node2D(xx,yy) );
                    }
                }

                char fpts[128];
                sprintf(fpts, "/tmp/tmppts%d", getpid());
                char fcells[128];
                sprintf(fcells, "/tmp/tmpcells%d", getpid());

                std::ofstream ftmp(fpts);
                // save corners of bounding box first, to work around meshing
                // problems in voronoi utility
                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
                for (int i=0; i < nxy; i++) {
                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;

                }
                ftmp.close();

                char cmdstr[512];
                sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
                if (system(cmdstr) == 0) {
                    int cx, cy, cz;
                    std::ifstream ftri(fcells);
                    while (!ftri.eof()) {
                        ftri.getline(line,sizeof(line)-1);
                        if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
                            if (cx >= 4 && cy >= 4 && cz >= 4) {
                                // skip first 4 boundary points
                                xymesh.addCell(cx-4, cy-4, cz-4);
                            }
                        }
                    }
                    ftri.close();
                } else {
                    result.error("triangularization failed");
                    return false;
                }
                unlink(fpts);
                unlink(fcells);
            } else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
                // found z-grid
            } else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
                // found origin
            } else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
                // found one of the delta lines
                if (ddx != 0.0) { dx = ddx; }
                else if (ddy != 0.0) { dy = ddy; }
                else if (ddz != 0.0) { dz = ddz; }
            } else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
                if (isrect && (npts != nx*ny*nz)) {
                    result.addError("inconsistent data: expected %d points"
                                    " but found %d points", nx*ny*nz, npts);
                    return false;
                } else if (!isrect && (npts != nxy*nz)) {
                    result.addError("inconsistent data: expected %d points"
                                    " but found %d points", nx*ny*nz, npts);
                    return false;
                }
                break;
            } else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
                if (npts != nx*ny*nz) {
                    result.addError("inconsistent data: expected %d points"
                                    " but found %d points", nx*ny*nz, npts);
                    return false;
                }
                break;
            }
        }
    }
    // read data points
    if (fin.eof()) {
        result.error("data not found in stream");
        return false;
    }
    if (isrect) {
        Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
        Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
        Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
        Rappture::FieldRect3D field(xgrid, ygrid, zgrid);

        double dval[6];
        int nread = 0;
        int ix = 0;
        int iy = 0;
        int iz = 0;
        while (!fin.eof() && nread < npts) {
            fin.getline(line,sizeof(line)-1);
            if (fin.fail()) {
                result.addError("error reading data points");
                return false;
            }
            int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
            
            for (int p=0; p < n; p++) {
                int nindex = iz*nx*ny + iy*nx + ix;
                field.define(nindex, dval[p]);
                fflush(stderr);
                nread++;
                if (++iz >= nz) {
                    iz = 0;
                    if (++iy >= ny) {
                        iy = 0;
                        ++ix;
                    }
                }
            }
        }
        
        // make sure that we read all of the expected points
        if (nread != nx*ny*nz) {
            result.addError("inconsistent data: expected %d points"
                            " but found %d points", nx*ny*nz, npts);
            return false;
        }
        
        // figure out a good mesh spacing
        int nsample = 30;
        dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
        dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
        dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
        double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
        
        nx = (int)ceil(dx/dmin);
        ny = (int)ceil(dy/dmin);
        nz = (int)ceil(dz/dmin);
        
#ifndef NV40
        // must be an even power of 2 for older cards
        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
#endif
        
        //#define _SOBEL
#ifdef _SOBEL_
        const int step = 1;
        float *cdata = new float[nx*ny*nz * step];
        int ngen = 0;
        double nzero_min = 0.0;
        for (int iz=0; iz < nz; iz++) {
            double zval = z0 + iz*dmin;
            for (int iy=0; iy < ny; iy++) {
                double yval = y0 + iy*dmin;
                for (int ix=0; ix < nx; ix++) {
                    double xval = x0 + ix*dmin;
                    double v = field.value(xval,yval,zval);
                    
                    if (v != 0.0f && v < nzero_min) {
                        nzero_min = v;
                    }
                    
                    // scale all values [0-1], -1 => out of bounds
                    v = (isnan(v)) ? -1.0 : v;
                    
                    cdata[ngen] = v;
                    ngen += step;
                }
            }
        }
        
        float* data = computeGradient(cdata, nx, ny, nz, field.valueMin(),
                                      field.valueMax());
#else
        double vmin = field.valueMin();
        double vmax = field.valueMax();
        double nzero_min = 0;
        float *data = new float[nx*ny*nz * 4];
        double dv = vmax - vmin;
        int ngen = 0;
        if (dv == 0.0)  dv = 1.0; 
        
        for (int iz=0; iz < nz; iz++) {
            double zval = z0 + iz*dmin;
            for (int iy=0; iy < ny; iy++) {
                double yval = y0 + iy*dmin;
                for (int ix=0; ix < nx; ix++) {
                    double xval = x0 + ix*dmin;
                    double v = field.value(xval,yval,zval);
                    
                    // scale all values [0-1], -1 => out of bounds
                    v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
                    
                    data[ngen] = v;
                    ngen += 4;
                }
            }
        }
        
        computeSimpleGradient(data, nx, ny, nz);
#endif
        
        for (int i=0; i<nx*ny*nz; i++) {
            fprintf(stderr,"enddata[%i] = %lg\n",i,data[i]);
            fflush(stderr);
        }
        
        fprintf(stdout,"End Data Stats index = %i\n",index);
        fprintf(stdout,"nx = %i ny = %i nz = %i\n",nx,ny,nz);
        fprintf(stdout,"dx = %lg dy = %lg dz = %lg\n",dx,dy,dz);
        fprintf(stdout,"dataMin = %lg\tdataMax = %lg\tnzero_min = %lg\n", field.valueMin(),field.valueMax(),nzero_min);
        fflush(stdout);
        
        Volume *volPtr;
        volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
                                      field.valueMin(), field.valueMax(), nzero_min);
        volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
                               field.rangeMax(Rappture::xaxis));
        volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
                               field.rangeMax(Rappture::yaxis));
        volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
                               field.rangeMax(Rappture::zaxis));
        volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
        volPtr->update_pending = true;
        // TBD..
        // POINTSET
        /*
          PointSet* pset = new PointSet();
          pset->initialize(volume[index], (float*) data);
          pset->setVisible(true);
          NanoVis::pointSet.push_back(pset);
          updateColor(pset);
          NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
        */
        
        delete [] data;
        
    } else {
        Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
        Rappture::FieldPrism3D field(xymesh, zgrid);
        
        double dval;
        int nread = 0;
        int ixy = 0;
        int iz = 0;
        while (!fin.eof() && nread < npts) {
            fin >> dval;
            if (fin.fail()) {
                result.addError("after %d of %d points: can't read number", 
                                nread, npts);
                return false;
            } else {
                int nid = nxy*iz + ixy;
                field.define(nid, dval);
                
                nread++;
                if (++iz >= nz) {
                    iz = 0;
                    ixy++;
                }
            }
        }
        
        // make sure that we read all of the expected points
        if (nread != nxy*nz) {
            result.addError("inconsistent data: expected %d points"
                            " but found %d points", nx*ny*nz, npts);
            return false;
        }
        
        // figure out a good mesh spacing
        int nsample = 30;
        x0 = field.rangeMin(Rappture::xaxis);
        dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
        y0 = field.rangeMin(Rappture::yaxis);
        dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
        z0 = field.rangeMin(Rappture::zaxis);
        dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
        double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
        
        nx = (int)ceil(dx/dmin);
        ny = (int)ceil(dy/dmin);
        nz = (int)ceil(dz/dmin);
#ifndef NV40
        // must be an even power of 2 for older cards
        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
#endif
        float *data = new float[4*nx*ny*nz];
        
        double vmin = field.valueMin();
        double dv = field.valueMax() - field.valueMin();
        if (dv == 0.0) { dv = 1.0; }
        
        // generate the uniformly sampled data that we need for a volume
        int ngen = 0;
        double nzero_min = 0.0;
        for (iz=0; iz < nz; iz++) {
            double zval = z0 + iz*dmin;
            for (int iy=0; iy < ny; iy++) {
                double yval = y0 + iy*dmin;
                for (int ix=0; ix < nx; ix++) {
                    double xval = x0 + ix*dmin;
                    double v = field.value(xval,yval,zval);
                    
                    if (v != 0.0f && v < nzero_min) {
                        nzero_min = v;
                    }
                    // scale all values [0-1], -1 => out of bounds
                    v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
                    data[ngen] = v;
                    
                    ngen += 4;
                }
            }
        }
        
        // Compute the gradient of this data.  BE CAREFUL: center
        // calculation on each node to avoid skew in either direction.
        ngen = 0;
        for (int iz=0; iz < nz; iz++) {
            for (int iy=0; iy < ny; iy++) {
                for (int ix=0; ix < nx; ix++) {
                    // gradient in x-direction
                    double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
                    double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
                    if (valm1 < 0 || valp1 < 0) {
                        data[ngen+1] = 0.0;
                    } else {
                        data[ngen+1] = valp1-valm1; // assume dx=1
                        //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1 (ISO)
                    }
                    
                    // gradient in y-direction
                    valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
                    valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
                    if (valm1 < 0 || valp1 < 0) {
                        data[ngen+2] = 0.0;
                    } else {
                        data[ngen+2] = valp1-valm1; // assume dy=1
                        //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1 (ISO)
                    }
                    
                    // gradient in z-direction
                    valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
                    valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
                    if (valm1 < 0 || valp1 < 0) {
                        data[ngen+3] = 0.0;
                    } else {
                        data[ngen+3] = valp1-valm1; // assume dz=1
                        //data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1 (ISO)
                    }
                    
                    ngen += 4;
                }
            }
        }
        
        Volume *volPtr;
        volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
                                      field.valueMin(), field.valueMax(), nzero_min);
        volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
                               field.rangeMax(Rappture::xaxis));
        volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
                               field.rangeMax(Rappture::yaxis));
        volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
                               field.rangeMax(Rappture::zaxis));
        volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
        volPtr->update_pending = true;
        // TBD..
        // POINTSET
        /*
          PointSet* pset = new PointSet();
          pset->initialize(volume[index], (float*) data);
          pset->setVisible(true);
          NanoVis::pointSet.push_back(pset);
          updateColor(pset);
          NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
        */
        
        
        delete [] data;
    }

    //
    // Center this new volume on the origin.
    //
    float dx0 = -0.5;
    float dy0 = -0.5*dy/dx;
    float dz0 = -0.5*dz/dx;
    NanoVis::volume[index]->move(Vector3(dx0, dy0, dz0));
    return true;
}


bool
load_volume_stream_insoo(Rappture::Outcome &result, int index, 
                         std::iostream& fin)
{
    printf("load_volume_stream\n");
 
    Rappture::MeshTri2D xymesh;
    int dummy, nx, ny, nz, nxy, npts;
    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
    char line[128], type[128], *start;

    int isrect = 1;

    dx = dy = dz = 0.0;         // Suppress compiler warning.
    x0 = y0 = z0 = 0.0;         // May not have an origin line.
    while (!fin.eof()) {
        fin.getline(line, sizeof(line) - 1);
        if (fin.fail()) {
            result.addError("line \"%s\"error in data stream");
            return false;
        }
        for (start=line; *start == ' ' || *start == '\t'; start++)
            ;  // skip leading blanks

        if (*start != '#') {  // skip comment lines
            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
                // found grid size
                isrect = 1;
            } else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
                isrect = 0;

                double xx, yy, zz;
                for (int i=0; i < nxy; i++) {
                    fin.getline(line,sizeof(line)-1);
                    if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
                        xymesh.addNode( Rappture::Node2D(xx,yy) );
                    }
                }

                char fpts[128];
                sprintf(fpts, "/tmp/tmppts%d", getpid());
                char fcells[128];
                sprintf(fcells, "/tmp/tmpcells%d", getpid());

                std::ofstream ftmp(fpts);
                // save corners of bounding box first, to work around meshing
                // problems in voronoi utility
                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
                for (int i=0; i < nxy; i++) {
                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;

                }
                ftmp.close();

                char cmdstr[512];
                sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
                if (system(cmdstr) == 0) {
                    int cx, cy, cz;
                    std::ifstream ftri(fcells);
                    while (!ftri.eof()) {
                        ftri.getline(line,sizeof(line)-1);
                        if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
                            if (cx >= 4 && cy >= 4 && cz >= 4) {
                                // skip first 4 boundary points
                                xymesh.addCell(cx-4, cy-4, cz-4);
                            }
                        }
                    }
                    ftri.close();
                } else {
                    result.error("triangularization failed");
                    return false;
                }
                unlink(fpts), unlink(fcells);
            } else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
                // found z-grid
            } else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
                // found origin
            } else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
                // found one of the delta lines
                if (ddx != 0.0) { dx = ddx; }
                else if (ddy != 0.0) { dy = ddy; }
                else if (ddz != 0.0) { dz = ddz; }
            } else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
                if (isrect && (npts != nx*ny*nz)) {
                    result.addError("inconsistent data: expected %d points"
                                    " but found %d points", nx*ny*nz, npts);
                    return false;
                } else if (!isrect && (npts != nxy*nz)) {
                    result.addError("inconsistent data: expected %d points"
                                    " but found %d points", nx*ny*nz, npts);
                    return false;
                }
                break;
            } else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
                if (npts != nx*ny*nz) {
                    result.addError("inconsistent data: expected %d points"
                                    " but found %d points", nx*ny*nz, npts);
                    return false;
                }
                break;
            }
        }
    }

    // read data points
    if (fin.eof()) {
        result.error("data not found in stream");
        return false;
    }

    if (isrect) {
        Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
        Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
        Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
        Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
        
        double dval[6];
        int nread = 0;
        int ix = 0;
        int iy = 0;
        int iz = 0;
        while (!fin.eof() && nread < npts) {
            fin.getline(line,sizeof(line)-1);
            if (fin.fail()) {
                result.error("error reading data points");
                return false;
            }
            int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
            
            for (int p=0; p < n; p++) {
                int nindex = iz*nx*ny + iy*nx + ix;
                field.define(nindex, dval[p]);
                fprintf(stderr,"nindex = %i\tdval[%i] = %lg\n", nindex, p,
                        dval[p]);
                fflush(stderr);
                nread++;
                if (++iz >= nz) {
                    iz = 0;
                    if (++iy >= ny) {
                        iy = 0;
                        ++ix;
                    }
                }
            }
        }
        
        // make sure that we read all of the expected points
        if (nread != nx*ny*nz) {
            result.addError("inconsistent data: expected %d points"
                            " but found %d points", nx*ny*nz, npts);
            return false;
        }
        
        // figure out a good mesh spacing
        int nsample = 30;
        dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
        dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
        dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
        double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
        
        nx = (int)ceil(dx/dmin);
        ny = (int)ceil(dy/dmin);
        nz = (int)ceil(dz/dmin);
        
#ifndef NV40
        // must be an even power of 2 for older cards
        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
#endif
        
        float *data = new float[4*nx*ny*nz];
        
        double vmin = field.valueMin();
        double dv = field.valueMax() - field.valueMin();
        if (dv == 0.0) { dv = 1.0; }
        
        int ngen = 0;
        double nzero_min = 0.0;
        for (iz=0; iz < nz; iz++) {
            double zval = z0 + iz*dmin;
            for (int iy=0; iy < ny; iy++) {
                double yval = y0 + iy*dmin;
                for (int ix=0; ix < nx; ix++) {
                    double xval = x0 + ix*dmin;
                    double v = field.value(xval,yval,zval);
                    
                    if (v != 0.0f && v < nzero_min) {
                        nzero_min = v;
                    }
                    // scale all values [0-1], -1 => out of bounds
                    v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
                    data[ngen] = v;
                    
                    ngen += 4;
                }
            }
        }
        
        // Compute the gradient of this data.  BE CAREFUL: center
        // calculation on each node to avoid skew in either direction.
        ngen = 0;
        for (int iz=0; iz < nz; iz++) {
            for (int iy=0; iy < ny; iy++) {
                for (int ix=0; ix < nx; ix++) {
                    // gradient in x-direction
                    double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
                    double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
                    if (valm1 < 0 || valp1 < 0) {
                        data[ngen+1] = 0.0;
                    } else {
                        data[ngen+1] = valp1-valm1; // assume dx=1
                        //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1 (ISO)
                    }
                    
                    // gradient in y-direction
                    valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
                    valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
                    if (valm1 < 0 || valp1 < 0) {
                        data[ngen+2] = 0.0;
                    } else {
                        data[ngen+2] = valp1-valm1; // assume dy=1
                        //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1 (ISO)
                    }
                    
                    // gradient in z-direction
                    valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
                    valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
                    if (valm1 < 0 || valp1 < 0) {
                        data[ngen+3] = 0.0;
                    } else {
                        data[ngen+3] = valp1-valm1; // assume dz=1
                        //data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1 (ISO)
                    }
                    
                    ngen += 4;
                }
            }
        }
        
        /*
          float *cdata = new float[nx*ny*nz];
          int ngen = 0;
          double nzero_min = 0.0;
          for (int iz=0; iz < nz; iz++) {
          double zval = z0 + iz*dmin;
          for (int iy=0; iy < ny; iy++) {
          double yval = y0 + iy*dmin;
          for (int ix=0; ix < nx; ix++) {
          double xval = x0 + ix*dmin;
          double v = field.value(xval,yval,zval);
          
          if (v != 0.0f && v < nzero_min) {
          nzero_min = v;
          }
          
          // scale all values [0-1], -1 => out of bounds
          v = (isnan(v)) ? -1.0 : v;
          
          cdata[ngen] = v;
          ++ngen;
          }
                }
            }
            
            float* data = computeGradient(cdata, nx, ny, nz, field.valueMin(),
                                          field.valueMax());
                                          
            for (int i=0; i<nx*ny*nz; i++) {
                fprintf(stderr,"enddata[%i] = %lg\n",i,data[i]);
                fflush(stderr);
            }

            fprintf(stdout,"End Data Stats index = %i\n",index);
            fprintf(stdout,"nx = %i ny = %i nz = %i\n",nx,ny,nz);
            fprintf(stdout,"dx = %lg dy = %lg dz = %lg\n",dx,dy,dz);
            fprintf(stdout,"dataMin = %lg\tdataMax = %lg\tnzero_min = %lg\n", field.valueMin(),field.valueMax(),nzero_min);
            fflush(stdout);
            */

        Volume *volPtr;
        volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
                                      field.valueMin(), field.valueMax(), nzero_min);
        volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
                               field.rangeMax(Rappture::xaxis));
        volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
                               field.rangeMax(Rappture::yaxis));
        volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
                               field.rangeMax(Rappture::zaxis));
        volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
        volPtr->update_pending = true;
        // TBD..
        // POINTSET
        /*
          PointSet* pset = new PointSet();
          pset->initialize(volume[index], (float*) data);
          pset->setVisible(true);
          NanoVis::pointSet.push_back(pset);
          updateColor(pset);
          NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
        */
        
        delete [] data;
        
    } else {
        Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
        Rappture::FieldPrism3D field(xymesh, zgrid);
        
        double dval;
        int nread = 0;
        int ixy = 0;
        int iz = 0;
        while (!fin.eof() && nread < npts) {
            fin >> dval;
            if (fin.fail()) {
                char mesg[256];
                sprintf(mesg,"after %d of %d points: can't read number", 
                        nread, npts);
                result.error(mesg);
                return false;
            } else {
                int nid = nxy*iz + ixy;
                field.define(nid, dval);
                
                nread++;
                if (++iz >= nz) {
                    iz = 0;
                    ixy++;
                }
            }
        }
        
        // make sure that we read all of the expected points
        if (nread != nxy*nz) {
            result.addError("inconsistent data: expected %d points"
                            " but found %d points", nx*ny*nz, npts);
            return false;
        }
        
        // figure out a good mesh spacing
        int nsample = 30;
        x0 = field.rangeMin(Rappture::xaxis);
        dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
        y0 = field.rangeMin(Rappture::yaxis);
        dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
        z0 = field.rangeMin(Rappture::zaxis);
        dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
        double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
        
        nx = (int)ceil(dx/dmin);
        ny = (int)ceil(dy/dmin);
        nz = (int)ceil(dz/dmin);
#ifndef NV40
        // must be an even power of 2 for older cards
        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
#endif
        float *data = new float[4*nx*ny*nz];
        
        double vmin = field.valueMin();
        double dv = field.valueMax() - field.valueMin();
        if (dv == 0.0) { dv = 1.0; }
        
        // generate the uniformly sampled data that we need for a volume
        int ngen = 0;
        double nzero_min = 0.0;
        for (iz=0; iz < nz; iz++) {
            double zval = z0 + iz*dmin;
            for (int iy=0; iy < ny; iy++) {
                double yval = y0 + iy*dmin;
                for (int ix=0; ix < nx; ix++) {
                    double xval = x0 + ix*dmin;
                    double v = field.value(xval,yval,zval);
                    
                    if (v != 0.0f && v < nzero_min) {
                        nzero_min = v;
                    }
                    // scale all values [0-1], -1 => out of bounds
                    v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
                    data[ngen] = v;
                    
                    ngen += 4;
                }
            }
        }
        
        // Compute the gradient of this data.  BE CAREFUL: center
        // calculation on each node to avoid skew in either direction.
        ngen = 0;
        for (int iz=0; iz < nz; iz++) {
            for (int iy=0; iy < ny; iy++) {
                for (int ix=0; ix < nx; ix++) {
                    // gradient in x-direction
                    double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
                    double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
                    if (valm1 < 0 || valp1 < 0) {
                        data[ngen+1] = 0.0;
                    } else {
                        data[ngen+1] = valp1-valm1; // assume dx=1
                        //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1 (ISO)
                    }
                    
                    // gradient in y-direction
                    valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
                    valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
                    if (valm1 < 0 || valp1 < 0) {
                        data[ngen+2] = 0.0;
                    } else {
                        data[ngen+2] = valp1-valm1; // assume dy=1
                        //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1 (ISO)
                    }
                    
                    // gradient in z-direction
                    valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
                    valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
                    if (valm1 < 0 || valp1 < 0) {
                        data[ngen+3] = 0.0;
                    } else {
                        data[ngen+3] = valp1-valm1; // assume dz=1
                        //data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1 (ISO)
                    }
                    
                    ngen += 4;
                }
            }
        }
        
        Volume *volPtr;
        volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
                                      field.valueMin(), field.valueMax(), nzero_min);
        volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
                               field.rangeMax(Rappture::xaxis));
        volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
                               field.rangeMax(Rappture::yaxis));
        volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
                               field.rangeMax(Rappture::zaxis));
        volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
        volPtr->update_pending = true;
        // TBD..
        // POINTSET
        /*
          PointSet* pset = new PointSet();
          pset->initialize(volume[index], (float*) data);
          pset->setVisible(true);
          NanoVis::pointSet.push_back(pset);
          updateColor(pset);
          NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
        */
        
        
        delete [] data;
    }

    //
    // Center this new volume on the origin.
    //
    float dx0 = -0.5;
    float dy0 = -0.5*dy/dx;
    float dz0 = -0.5*dz/dx;
    NanoVis::volume[index]->move(Vector3(dx0, dy0, dz0));
    return true;
}