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

Last change on this file since 1382 was 1382, checked in by gah, 15 years ago

fixed encoding problems

File size: 50.7 KB
Line 
1 
2/*
3 * ----------------------------------------------------------------------
4 * Nanovis: Visualization of Nanoelectronics Data
5 *
6 *  dxReader.cpp
7 *      This module contains openDX readers for 2D and 3D volumes.
8 *
9 * ======================================================================
10 *  AUTHOR:  Wei Qiao <qiaow@purdue.edu>
11 *           Michael McLennan <mmclennan@purdue.edu>
12 *           Purdue Rendering and Perceptualization Lab (PURPL)
13 *
14 *  Copyright (c) 2004-2006  Purdue Research Foundation
15 *
16 *  See the file "license.terms" for information on usage and
17 *  redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES.
18 * ======================================================================
19 */
20
21// common dx functions
22#include "dxReaderCommon.h"
23
24#include <stdio.h>
25#include <math.h>
26#include <fstream>
27#include <iostream>
28#include <sstream>
29#include <string>
30#include <sys/types.h>
31#include <unistd.h>
32
33#include "Nv.h"
34
35#include "nanovis.h"
36#include "RpField1D.h"
37#include "RpFieldRect3D.h"
38#include "RpFieldPrism3D.h"
39#include <Unirect.h>
40
41//transfer function headers
42#include "ZincBlendeVolume.h"
43#include "NvZincBlendeReconstructor.h"
44
45#define  _LOCAL_ZINC_TEST_ 0
46
47static INLINE char *
48skipspaces(char *string)
49{
50    while (isspace(*string)) {
51        string++;
52    }
53    return string;
54}
55
56static INLINE char *
57getline(char **stringPtr, char *endPtr)
58{
59    char *line, *p;
60
61    line = skipspaces(*stringPtr);
62    for (p = line; p < endPtr; p++) {
63        if (*p == '\n') {
64            *p++ = '\0';
65            *stringPtr = p;
66            return line;
67        }
68    }
69    *stringPtr = p;
70    return line;
71}
72
73Rappture::Unirect3d *
74ReadDxVectorFieldData(Rappture::Outcome &result, size_t length, char *string)
75{
76    int nx, ny, nz, npts;
77    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
78    char *p, *endPtr;
79
80
81    dx = dy = dz = 0.0;         // Suppress compiler warning.
82    x0 = y0 = z0 = 0.0;         // May not have an origin line.
83    for (p = string, endPtr = p + length; p < endPtr; /*empty*/) {
84        char *line;
85
86        line = getline(&p, endPtr);
87        if (line == endPtr) {
88            break;
89        }
90        if ((line[0] == '#') || (line == '\0')) {
91            continue;           // Skip blank or comment lines.
92        }
93        if (sscanf(line, "object %*d class gridpositions counts %d %d %d",
94                   &nx, &ny, &nz) == 4) {
95            printf("w:%d h:%d d:%d\n", nx, ny, nz);
96            // found grid size
97        } else if (sscanf(line, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
98            // found origin
99        } else if (sscanf(line, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
100            // found one of the delta lines
101            if (ddx != 0.0) {
102                dx = ddx;
103            } else if (ddy != 0.0) {
104                dy = ddy;
105            } else if (ddz != 0.0) {
106                dz = ddz;
107            }
108        } else if (sscanf(line, "object %*d class array type %*s shape 3"
109                " rank 1 items %d data follows", &npts) == 3) {
110            printf("point %d\n", npts);
111            if (npts != nx*ny*nz) {
112                result.addError("inconsistent data: expected %d points"
113                                " but found %d points", nx*ny*nz, npts);
114                return NULL;
115            }
116            break;
117        } else if (sscanf(line, "object %*d class array type %*s rank 0"
118                " times %d data follows", &npts) == 3) {
119            if (npts != nx*ny*nz) {
120                result.addError("inconsistent data: expected %d points"
121                                " but found %d points", nx*ny*nz, npts);
122                return NULL;
123            }
124            break;
125        }
126    }
127    if (npts != nx*ny*nz) {
128        result.addError("inconsistent data: expected %d points"
129                        " but found %d points", nx*ny*nz, npts);
130        return NULL;
131    }
132    float *values = new float[npts];
133    int nValues = 0;
134    for (int ix=0; ix < nx; ix++) {
135        for (int iy=0; iy < ny; iy++) {
136            for (int iz=0; iz < nz; iz++) {
137                char *line;
138                if ((p == endPtr) || (nValues > npts)) {
139                    break;
140                }
141                line = getline(&p, endPtr);
142                if ((line[0] == '#') || (line[0] == '\0')) {
143                    continue;   // Skip blank or comment lines.
144                }
145                double vx, vy, vz;
146                if (sscanf(line, "%lg %lg %lg", &vx, &vy, &vz) == 3) {
147                    int nindex = (iz*nx*ny + iy*nx + ix) * 3;
148                    values[nindex] = vx;
149                    values[nindex+1] = vy;
150                    values[nindex+2] = vz;
151                    nValues++;
152                }
153            }
154        }
155    }
156    // make sure that we read all of the expected points
157    if (nValues != npts) {
158        result.addError("inconsistent data: expected %d points"
159                        " but found %d points", npts, nValues);
160        delete values;
161        return NULL;
162    }
163    return new Rappture::Unirect3d(x0, x0 + dx * nx, nx, y0, y0 + dy * ny, ny,
164                                   z0, z0 + dz * nz, nz, nValues, values);
165}
166
167
168/*
169 * Load a 3D vector field from a dx-format file
170 */
171bool
172load_vector_stream2(Rappture::Outcome &result, int ivol, size_t length,
173                    char *string)
174{
175    Rappture::Unirect3d *dataPtr;
176    dataPtr = ReadDxVectorFieldData(result, length, string);
177    if (dataPtr == NULL) {
178        return false;
179    }
180    Vector3 physicalMin;
181    Vector3 physicalMax;
182
183    physicalMin.set(dataPtr->xMin(), dataPtr->yMin(), dataPtr->zMin());
184    physicalMax.set(dataPtr->xMax(), dataPtr->yMax(), dataPtr->zMax());
185
186    double max_mag = -DBL_MAX, min_mag = DBL_MAX;
187    float *values = dataPtr->values();
188    for (size_t ix=0; ix < dataPtr->xNum(); ix++) {
189        for (size_t iy=0; iy < dataPtr->yNum(); iy++) {
190            for (size_t iz=0; iz < dataPtr->zNum(); iz++) {
191                double vx, vy, vz, vm;
192                vx = values[0];
193                vy = values[1];
194                vz = values[2];
195                vm = sqrt(vx*vx + vy*vy + vz*vz);
196                if (vm > max_mag) {
197                    max_mag = vm;
198                }
199                if (vm < min_mag) {
200                    min_mag = vm;
201                }
202                values += 3;
203            }
204        }
205    }
206   
207    float *data = new float[4*dataPtr->nValues()];
208    memset(data, 0, sizeof(float) * 4 * dataPtr->nValues());
209    fprintf(stderr, "generating %dx%dx%d = %d points\n",
210            dataPtr->xNum(), dataPtr->yNum(), dataPtr->zNum(),
211            dataPtr->nValues());
212
213    float *destPtr = data;
214    float *srcPtr = dataPtr->values();
215    for (size_t iz=0; iz < dataPtr->zNum(); iz++) {
216        for (size_t iy=0; iy < dataPtr->yNum(); iy++) {
217            for (size_t ix=0; ix < dataPtr->xNum(); ix++) {
218                double vx, vy, vz, vm;
219                vx = values[0];
220                vy = values[1];
221                vz = values[2];
222                vm = sqrt(vx*vx + vy*vy + vz*vz);
223                destPtr[0] = vm / max_mag;
224                destPtr[1] = vx /(2.0*max_mag) + 0.5;
225                destPtr[2] = vy /(2.0*max_mag) + 0.5;
226                destPtr[3] = vz /(2.0*max_mag) + 0.5;
227                srcPtr += 3;
228                destPtr += 4;
229            }
230        }
231    }
232   
233    Volume *volPtr;
234    volPtr = NanoVis::load_volume(ivol, dataPtr->xNum(), dataPtr->yNum(),
235                dataPtr->zNum(), 4, data, min_mag, max_mag, 0);
236    volPtr->xAxis.SetRange(dataPtr->xMin(), dataPtr->xMax());
237    volPtr->yAxis.SetRange(dataPtr->yMin(), dataPtr->yMax());
238    volPtr->zAxis.SetRange(dataPtr->zMin(), dataPtr->zMax());
239    volPtr->wAxis.SetRange(min_mag, max_mag);
240    volPtr->update_pending = true;
241    volPtr->setPhysicalBBox(physicalMin, physicalMax);
242    delete dataPtr;
243    delete data;
244    return true;
245}
246
247bool
248load_vector_stream(Rappture::Outcome result, int index, size_t length,
249                    char *string)
250{
251    Rappture::Unirect3d *dataPtr;
252    dataPtr = ReadDxVectorFieldData(result, length, string);
253    if (dataPtr == NULL) {
254        return false;
255    }
256    Vector3 physicalMin;
257    Vector3 physicalMax;
258
259    physicalMin.set(dataPtr->xMin(), dataPtr->yMin(), dataPtr->zMin());
260    physicalMax.set(dataPtr->xMax(), dataPtr->yMax(), dataPtr->zMax());
261
262    Rappture::Mesh1D xgrid(dataPtr->xMin(), dataPtr->yMin(), dataPtr->xNum());
263    Rappture::Mesh1D ygrid(dataPtr->yMin(), dataPtr->yMax(), dataPtr->yNum());
264    Rappture::Mesh1D zgrid(dataPtr->zMin(), dataPtr->zMax(), dataPtr->zNum());
265    Rappture::FieldRect3D xfield(xgrid, ygrid, zgrid);
266    Rappture::FieldRect3D yfield(xgrid, ygrid, zgrid);
267    Rappture::FieldRect3D zfield(xgrid, ygrid, zgrid);
268   
269    float *values = dataPtr->values();
270    size_t npts = 0;
271    for (size_t ix=0; ix < dataPtr->xNum(); ix++) {
272        for (size_t iy=0; iy < dataPtr->yNum(); iy++) {
273            for (size_t iz=0; iz < dataPtr->zNum(); iz++) {
274                xfield.define(npts, values[0]);
275                yfield.define(npts, values[1]);
276                zfield.define(npts, values[2]);
277                npts++;
278                values += 3;
279            }
280        }
281    }
282
283    double dx, dy, dz;
284    // figure out a good mesh spacing
285    int nsample = 30;
286    dx = xfield.rangeMax(Rappture::xaxis) - xfield.rangeMin(Rappture::xaxis);
287    dy = xfield.rangeMax(Rappture::yaxis) - xfield.rangeMin(Rappture::yaxis);
288    dz = xfield.rangeMax(Rappture::zaxis) - xfield.rangeMin(Rappture::zaxis);
289    double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
290
291    printf("dx:%lf dy:%lf dz:%lf dmin:%lf\n", dx, dy, dz, dmin);
292
293    size_t nx, ny, nz;
294    nx = (int)ceil(dx/dmin);
295    ny = (int)ceil(dy/dmin);
296    nz = (int)ceil(dz/dmin);
297   
298#ifndef NV40
299    // must be an even power of 2 for older cards
300    nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
301    ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
302    nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
303#endif
304   
305    float *data = new float[4*nx*ny*nz];
306    memset(data, 0, sizeof(float) * 4 * nx * ny * nz);
307   
308    std::cout << "generating " << nx << "x" << ny << "x" << nz << " = " << nx*ny*nz << " points" << std::endl;
309   
310    // generate the uniformly sampled data that we need for a volume
311    double vmin = 1e21;
312    double vmax = -1e21;
313    double nzero_min = 0.0;
314    size_t ngen = 0;
315    for (size_t iz=0; iz < nz; iz++) {
316        double zval = dataPtr->zMin() + iz*dmin;
317        for (size_t iy=0; iy < ny; iy++) {
318            double yval = dataPtr->yMin() + iy*dmin;
319            for (size_t ix=0; ix < nx; ix++) {
320                double xval = dataPtr->xMin() + ix*dmin;
321                double vx, vy, vz;
322
323                vx = xfield.value(xval,yval,zval);
324                vy = yfield.value(xval,yval,zval);
325                vz = zfield.value(xval,yval,zval);
326               
327                double vm;
328                vm = sqrt(vx*vx + vy*vy + vz*vz);
329                if (vm < vmin) {
330                    vmin = vm;
331                } else if (vm > vmax) {
332                    vmax = vm;
333                }
334                if ((vm != 0.0f) && (vm < nzero_min)) {
335                    nzero_min = vm;
336                }
337                data[ngen++] = vm;
338                data[ngen++] = vx;
339                data[ngen++] = vy;
340                data[ngen++] = vz;
341            }
342        }
343    }
344   
345    ngen = 0;
346    // scale should be accounted.
347    for (ngen=0; ngen < npts; /*empty*/) {
348        data[ngen] = data[ngen] / vmax; ++ngen;
349        data[ngen] = (data[ngen]/(2.0*vmax) + 0.5); ++ngen;
350        data[ngen] = (data[ngen]/(2.0*vmax) + 0.5); ++ngen;
351        data[ngen] = (data[ngen]/(2.0*vmax) + 0.5); ++ngen;
352    }
353    Volume *volPtr;
354    volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data, vmin, vmax,
355                                  nzero_min);
356   
357    volPtr->xAxis.SetRange(dataPtr->xMin(), dataPtr->xMin() + (nx * dx));
358    volPtr->yAxis.SetRange(dataPtr->yMin(), dataPtr->yMin() + (ny * dy));
359    volPtr->zAxis.SetRange(dataPtr->zMin(), dataPtr->zMin() + (nz * dz));
360    volPtr->wAxis.SetRange(vmin, vmax);
361    volPtr->update_pending = true;
362    delete [] data;
363    return true;
364}
365
366/* Load a 3D volume from a dx-format file
367 */
368bool
369load_volume_stream2(Rappture::Outcome &result, int index, std::iostream& fin)
370{
371    printf("load_volume_stream2\n");
372    Rappture::MeshTri2D xymesh;
373    int dummy, nx, ny, nz, nxy, npts;
374    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
375    char line[128], type[128], *start;
376
377    int isrect = 1;
378    dx = dy = dz = 0.0;         // Suppress compiler warning.
379    x0 = y0 = z0 = 0.0;         // May not have an origin line.
380    do {
381        fin.getline(line,sizeof(line)-1);
382        for (start=&line[0]; *start == ' ' || *start == '\t'; start++)
383            ;  // skip leading blanks
384
385        if (*start != '#') {  // skip comment lines
386            if (sscanf(start, "object %d class gridpositions counts %d %d %d",
387                       &dummy, &nx, &ny, &nz) == 4) {
388                // found grid size
389                isrect = 1;
390            } else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
391                isrect = 0;
392                double xx, yy, zz;
393                for (int i=0; i < nxy; i++) {
394                    fin.getline(line,sizeof(line)-1);
395                    if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
396                        xymesh.addNode( Rappture::Node2D(xx,yy) );
397                    }
398                }
399
400                char fpts[128];
401                sprintf(fpts, "/tmp/tmppts%d", getpid());
402                char fcells[128];
403                sprintf(fcells, "/tmp/tmpcells%d", getpid());
404
405                std::ofstream ftmp(fpts);
406                // save corners of bounding box first, to work around meshing
407                // problems in voronoi utility
408                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
409                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
410                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
411                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
412                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
413                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
414                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
415                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
416                for (int i=0; i < nxy; i++) {
417                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
418                }
419                ftmp.close();
420
421                char cmdstr[512];
422                sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
423                if (system(cmdstr) == 0) {
424                    int cx, cy, cz;
425                    std::ifstream ftri(fcells);
426                    while (!ftri.eof()) {
427                        ftri.getline(line,sizeof(line)-1);
428                        if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
429                            if (cx >= 4 && cy >= 4 && cz >= 4) {
430                                // skip first 4 boundary points
431                                xymesh.addCell(cx-4, cy-4, cz-4);
432                            }
433                        }
434                    }
435                    ftri.close();
436                } else {
437                    result.error("triangularization failed");
438                    return false;
439                }
440                unlink(fpts), unlink(fcells);
441            } else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
442                // found z-grid
443            } else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
444                // found origin
445            } else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
446        int count = 0;
447                // found one of the delta lines
448                if (ddx != 0.0) {
449            dx = ddx;
450            count++;
451        }
452        if (ddy != 0.0) {
453            dy = ddy;
454            count++;
455        }
456        if (ddz != 0.0) {
457            dz = ddz;
458            count++;
459        }
460        if (count > 1) {
461            result.addError("don't know how to handle multiple non-zero"
462                            " delta values");
463            return false;
464        }
465            } else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
466                if (isrect && (npts != nx*ny*nz)) {
467                    result.addError("inconsistent data: expected %d points "
468                                    " but found %d points", nx*ny*nz, npts);
469                    return false;
470                } else if (!isrect && (npts != nxy*nz)) {
471                    result.addError("inconsistent data: expected %d points "
472                                    " but found %d points", nxy*nz, npts);
473                    return false;
474                }
475                break;
476            } else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
477                if (npts != nx*ny*nz) {
478                    result.addError("inconsistent data: expected %d points "
479                                    " but found %d points", nx*ny*nz, npts);
480                    return false;
481                }
482                break;
483            }
484        }
485    } while (!fin.eof());
486
487    fprintf(stderr, "found nx=%d ny=%d, nz=%d, x0=%f, y0=%f, z0=%f\n",
488            nx, ny, nz, x0, y0, z0);
489    // read data points
490    if (fin.eof()) {
491        result.addError("EOF found: expecting %d points", npts);
492        return false;
493    }
494    if (isrect) {
495        double dval[6];
496        int nread = 0;
497        int ix = 0;
498        int iy = 0;
499        int iz = 0;
500        float* data = new float[nx *  ny *  nz * 4];
501        memset(data, 0, nx*ny*nz*4);
502        double vmin = 1e21;
503        double nzero_min = 1e21;
504        double vmax = -1e21;
505       
506       
507        while (!fin.eof() && nread < npts) {
508            fin.getline(line,sizeof(line)-1);
509            int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
510           
511            for (int p=0; p < n; p++) {
512                int nindex = (iz*nx*ny + iy*nx + ix) * 4;
513                data[nindex] = dval[p];
514               
515                if (dval[p] < vmin) {
516                    vmin = dval[p];
517                } else if (dval[p] > vmax) {
518                    vmax = dval[p];
519                }
520                if (dval[p] != 0.0f && dval[p] < nzero_min) {
521                    nzero_min = dval[p];
522                }
523               
524                nread++;
525                if (++iz >= nz) {
526                    iz = 0;
527                    if (++iy >= ny) {
528                        iy = 0;
529                        ++ix;
530                    }
531                }
532            }
533        }
534       
535        // make sure that we read all of the expected points
536        if (nread != nx*ny*nz) {
537            result.addError("inconsistent data: expected %d points "
538                            " but found %d points", nx*ny*nz, nread);
539            return false;
540        }
541       
542        double dv = vmax - vmin;
543        int count = nx*ny*nz;
544        int ngen = 0;
545        double v;
546        if (dv == 0.0) {
547            dv = 1.0;
548        }
549       
550        for (int i = 0; i < count; ++i) {
551            v = data[ngen];
552            // scale all values [0-1], -1 => out of bounds
553            //
554            // INSOO
555            v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
556            data[ngen] = v;
557            ngen += 4;
558        }
559       
560        computeSimpleGradient(data, nx, ny, nz);
561       
562        dx = nx;
563        dy = ny;
564        dz = nz;
565       
566        Volume *volPtr;
567        volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
568                                      vmin, vmax, nzero_min);
569        volPtr->xAxis.SetRange(x0, x0 + (nx * dx));
570        volPtr->yAxis.SetRange(y0, y0 + (ny * dy));
571        volPtr->zAxis.SetRange(z0, z0 + (nz * dz));
572        volPtr->wAxis.SetRange(vmin, vmax);
573        volPtr->update_pending = true;
574        delete [] data;
575       
576    } else {
577        Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
578        Rappture::FieldPrism3D field(xymesh, zgrid);
579       
580        double dval;
581        int nread = 0;
582        int ixy = 0;
583        int iz = 0;
584        while (!fin.eof() && nread < npts) {
585            fin >> dval;
586            if (fin.fail()) {
587                result.addError("after %d of %d points: can't read number",
588                                nread, npts);
589                return false;
590            } else {
591                int nid = nxy*iz + ixy;
592                field.define(nid, dval);
593               
594                nread++;
595                if (++iz >= nz) {
596                    iz = 0;
597                    ixy++;
598                }
599            }
600        }
601       
602        // make sure that we read all of the expected points
603        if (nread != nxy*nz) {
604            result.addError("inconsistent data: expected %d points "
605                            "but found %d points", nxy*nz, nread);
606            return false;
607        }
608       
609        // figure out a good mesh spacing
610        int nsample = 30;
611        x0 = field.rangeMin(Rappture::xaxis);
612        dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
613        y0 = field.rangeMin(Rappture::yaxis);
614        dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
615        z0 = field.rangeMin(Rappture::zaxis);
616        dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
617        double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
618       
619        nx = (int)ceil(dx/dmin);
620        ny = (int)ceil(dy/dmin);
621        nz = (int)ceil(dz/dmin);
622#ifndef NV40
623        // must be an even power of 2 for older cards
624        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
625        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
626        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
627#endif
628        float *data = new float[4*nx*ny*nz];
629       
630        double vmin = field.valueMin();
631        double dv = field.valueMax() - field.valueMin();
632        if (dv == 0.0) {
633            dv = 1.0;
634        }
635        // generate the uniformly sampled data that we need for a volume
636        int ngen = 0;
637        double nzero_min = 0.0;
638        for (iz=0; iz < nz; iz++) {
639            double zval = z0 + iz*dmin;
640            for (int iy=0; iy < ny; iy++) {
641                double yval = y0 + iy*dmin;
642                for (int ix=0; ix < nx; ix++) {
643                    double xval = x0 + ix*dmin;
644                    double v = field.value(xval,yval,zval);
645                   
646                    if (v != 0.0f && v < nzero_min) {
647                        nzero_min = v;
648                    }
649                    // scale all values [0-1], -1 => out of bounds
650                    v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
651                    data[ngen] = v;
652                   
653                    ngen += 4;
654                }
655            }
656        }
657       
658        // FIXME: This next section of code should be replaced by a
659        // call to the computeSimpleGradient() function. There is a slight
660        // difference in the code below and the aforementioned function
661        // in that the commented out lines in the else statements are
662        // different.
663        //
664        // Compute the gradient of this data.  BE CAREFUL: center
665        // calculation on each node to avoid skew in either direction.
666        ngen = 0;
667        for (int iz=0; iz < nz; iz++) {
668            for (int iy=0; iy < ny; iy++) {
669                for (int ix=0; ix < nx; ix++) {
670                    // gradient in x-direction
671                    //double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
672                    //double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
673                    double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
674                    double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
675                    if (valm1 < 0 || valp1 < 0) {
676                        data[ngen+1] = 0.0;
677                    } else {
678                        data[ngen+1] = valp1-valm1; // assume dx=1
679                        //data[ngen+1] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
680                    }
681                   
682                    // gradient in y-direction
683                    valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
684                    valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
685                    if (valm1 < 0 || valp1 < 0) {
686                        data[ngen+2] = 0.0;
687                    } else {
688                        data[ngen+2] = valp1-valm1; // assume dy=1
689                        //data[ngen+2] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
690                    }
691                   
692                    // gradient in z-direction
693                    valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
694                    valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
695                    if (valm1 < 0 || valp1 < 0) {
696                        data[ngen+3] = 0.0;
697                    } else {
698                        data[ngen+3] = valp1-valm1; // assume dz=1
699                        //data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
700                    }
701                   
702                    ngen += 4;
703                }
704            }
705        }
706       
707        Volume *volPtr;
708        volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
709                                      field.valueMin(), field.valueMax(), nzero_min);
710        volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
711                               field.rangeMax(Rappture::xaxis));
712        volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
713                               field.rangeMax(Rappture::yaxis));
714        volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
715                               field.rangeMax(Rappture::zaxis));
716        volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
717        volPtr->update_pending = true;
718        delete [] data;
719    }
720    //
721    // Center this new volume on the origin.
722    //
723    float dx0 = -0.5;
724    float dy0 = -0.5*dy/dx;
725    float dz0 = -0.5*dz/dx;
726    NanoVis::volume[index]->move(Vector3(dx0, dy0, dz0));
727    return true;
728}
729
730bool
731load_volume_stream(Rappture::Outcome &result, int index, std::iostream& fin)
732{
733    printf("load_volume_stream\n");
734
735    Rappture::MeshTri2D xymesh;
736    int dummy, nx, ny, nz, nxy, npts;
737    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
738    char line[128], type[128], *start;
739
740    int isrect = 1;
741
742    dx = dy = dz = 0.0;         // Suppress compiler warning.
743    x0 = y0 = z0 = 0.0;         // May not have an origin line.
744    while (!fin.eof()) {
745        fin.getline(line, sizeof(line) - 1);
746        fprintf(stderr, "line is %s\n", line);
747        if (fin.fail()) {
748            result.error("error in data stream");
749            return false;
750        }
751        for (start=line; *start == ' ' || *start == '\t'; start++)
752            ;  // skip leading blanks
753
754        if (*start != '#') {  // skip comment lines
755            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
756                // found grid size
757                isrect = 1;
758            } else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
759                isrect = 0;
760
761                double xx, yy, zz;
762                for (int i=0; i < nxy; i++) {
763                    fin.getline(line,sizeof(line)-1);
764                    if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
765                        xymesh.addNode( Rappture::Node2D(xx,yy) );
766                    }
767                }
768
769                char fpts[128];
770                sprintf(fpts, "/tmp/tmppts%d", getpid());
771                char fcells[128];
772                sprintf(fcells, "/tmp/tmpcells%d", getpid());
773
774                std::ofstream ftmp(fpts);
775                // save corners of bounding box first, to work around meshing
776                // problems in voronoi utility
777                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
778                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
779                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
780                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
781                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
782                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
783                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
784                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
785                for (int i=0; i < nxy; i++) {
786                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
787
788                }
789                ftmp.close();
790
791                char cmdstr[512];
792                sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
793                if (system(cmdstr) == 0) {
794                    int cx, cy, cz;
795                    std::ifstream ftri(fcells);
796                    while (!ftri.eof()) {
797                        ftri.getline(line,sizeof(line)-1);
798                        if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
799                            if (cx >= 4 && cy >= 4 && cz >= 4) {
800                                // skip first 4 boundary points
801                                xymesh.addCell(cx-4, cy-4, cz-4);
802                            }
803                        }
804                    }
805                    ftri.close();
806                } else {
807                    result.error("triangularization failed");
808                    return false;
809                }
810                unlink(fpts);
811                unlink(fcells);
812            } else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
813                // found z-grid
814            } else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
815                // found origin
816            } else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
817                // found one of the delta lines
818                if (ddx != 0.0) { dx = ddx; }
819                else if (ddy != 0.0) { dy = ddy; }
820                else if (ddz != 0.0) { dz = ddz; }
821            } else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
822                if (isrect && (npts != nx*ny*nz)) {
823                    result.addError("inconsistent data: expected %d points"
824                                    " but found %d points", nx*ny*nz, npts);
825                    return false;
826                } else if (!isrect && (npts != nxy*nz)) {
827                    result.addError("inconsistent data: expected %d points"
828                                    " but found %d points", nx*ny*nz, npts);
829                    return false;
830                }
831                break;
832            } else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
833                if (npts != nx*ny*nz) {
834                    result.addError("inconsistent data: expected %d points"
835                                    " but found %d points", nx*ny*nz, npts);
836                    return false;
837                }
838                break;
839            }
840        }
841    }
842    fprintf(stderr, "found nx=%d ny=%d, nz=%d, x0=%f, y0=%f, z0=%f\n",
843            nx, ny, nz, x0, y0, z0);
844    // read data points
845    if (fin.eof()) {
846        result.error("data not found in stream");
847        return false;
848    }
849    if (isrect) {
850        Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
851        Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
852        Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
853        Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
854
855        double dval[6];
856        int nread = 0;
857        int ix = 0;
858        int iy = 0;
859        int iz = 0;
860        while (!fin.eof() && nread < npts) {
861            fin.getline(line,sizeof(line)-1);
862            if (fin.fail()) {
863                result.addError("error reading data points");
864                return false;
865            }
866            int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
867           
868            for (int p=0; p < n; p++) {
869                int nindex = iz*nx*ny + iy*nx + ix;
870                field.define(nindex, dval[p]);
871                fprintf(stderr,"nindex = %i\tdval[%i] = %lg\n", nindex, p,
872                        dval[p]);
873                fflush(stderr);
874                nread++;
875                if (++iz >= nz) {
876                    iz = 0;
877                    if (++iy >= ny) {
878                        iy = 0;
879                        ++ix;
880                    }
881                }
882            }
883        }
884       
885        // make sure that we read all of the expected points
886        if (nread != nx*ny*nz) {
887            result.addError("inconsistent data: expected %d points"
888                            " but found %d points", nx*ny*nz, npts);
889            return false;
890        }
891       
892        // figure out a good mesh spacing
893        int nsample = 30;
894        dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
895        dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
896        dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
897        double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
898       
899        nx = (int)ceil(dx/dmin);
900        ny = (int)ceil(dy/dmin);
901        nz = (int)ceil(dz/dmin);
902       
903#ifndef NV40
904        // must be an even power of 2 for older cards
905        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
906        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
907        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
908#endif
909       
910        //#define _SOBEL
911#ifdef _SOBEL_
912        const int step = 1;
913        float *cdata = new float[nx*ny*nz * step];
914        int ngen = 0;
915        double nzero_min = 0.0;
916        for (int iz=0; iz < nz; iz++) {
917            double zval = z0 + iz*dmin;
918            for (int iy=0; iy < ny; iy++) {
919                double yval = y0 + iy*dmin;
920                for (int ix=0; ix < nx; ix++) {
921                    double xval = x0 + ix*dmin;
922                    double v = field.value(xval,yval,zval);
923                   
924                    if (v != 0.0f && v < nzero_min) {
925                        nzero_min = v;
926                    }
927                   
928                    // scale all values [0-1], -1 => out of bounds
929                    v = (isnan(v)) ? -1.0 : v;
930                   
931                    cdata[ngen] = v;
932                    ngen += step;
933                }
934            }
935        }
936       
937        float* data = computeGradient(cdata, nx, ny, nz, field.valueMin(),
938                                      field.valueMax());
939#else
940        double vmin = field.valueMin();
941        double vmax = field.valueMax();
942        double nzero_min = 0;
943        float *data = new float[nx*ny*nz * 4];
944        double dv = vmax - vmin;
945        int ngen = 0;
946        if (dv == 0.0)  dv = 1.0;
947       
948        for (int iz=0; iz < nz; iz++) {
949            double zval = z0 + iz*dmin;
950            for (int iy=0; iy < ny; iy++) {
951                double yval = y0 + iy*dmin;
952                for (int ix=0; ix < nx; ix++) {
953                    double xval = x0 + ix*dmin;
954                    double v = field.value(xval,yval,zval);
955                   
956                    // scale all values [0-1], -1 => out of bounds
957                    v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
958                   
959                    data[ngen] = v;
960                    ngen += 4;
961                }
962            }
963        }
964       
965        computeSimpleGradient(data, nx, ny, nz);
966#endif
967       
968        for (int i=0; i<nx*ny*nz; i++) {
969            fprintf(stderr,"enddata[%i] = %lg\n",i,data[i]);
970            fflush(stderr);
971        }
972       
973        fprintf(stdout,"End Data Stats index = %i\n",index);
974        fprintf(stdout,"nx = %i ny = %i nz = %i\n",nx,ny,nz);
975        fprintf(stdout,"dx = %lg dy = %lg dz = %lg\n",dx,dy,dz);
976        fprintf(stdout,"dataMin = %lg\tdataMax = %lg\tnzero_min = %lg\n", field.valueMin(),field.valueMax(),nzero_min);
977        fflush(stdout);
978       
979        Volume *volPtr;
980        volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
981                                      field.valueMin(), field.valueMax(), nzero_min);
982        volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
983                               field.rangeMax(Rappture::xaxis));
984        volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
985                               field.rangeMax(Rappture::yaxis));
986        volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
987                               field.rangeMax(Rappture::zaxis));
988        volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
989        volPtr->update_pending = true;
990        // TBD..
991        // POINTSET
992        /*
993          PointSet* pset = new PointSet();
994          pset->initialize(volume[index], (float*) data);
995          pset->setVisible(true);
996          NanoVis::pointSet.push_back(pset);
997          updateColor(pset);
998          NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
999        */
1000       
1001        delete [] data;
1002       
1003    } else {
1004        Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
1005        Rappture::FieldPrism3D field(xymesh, zgrid);
1006       
1007        double dval;
1008        int nread = 0;
1009        int ixy = 0;
1010        int iz = 0;
1011        while (!fin.eof() && nread < npts) {
1012            fin >> dval;
1013            if (fin.fail()) {
1014                result.addError("after %d of %d points: can't read number",
1015                                nread, npts);
1016                return false;
1017            } else {
1018                int nid = nxy*iz + ixy;
1019                field.define(nid, dval);
1020               
1021                nread++;
1022                if (++iz >= nz) {
1023                    iz = 0;
1024                    ixy++;
1025                }
1026            }
1027        }
1028       
1029        // make sure that we read all of the expected points
1030        if (nread != nxy*nz) {
1031            result.addError("inconsistent data: expected %d points"
1032                            " but found %d points", nx*ny*nz, npts);
1033            return false;
1034        }
1035       
1036        // figure out a good mesh spacing
1037        int nsample = 30;
1038        x0 = field.rangeMin(Rappture::xaxis);
1039        dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
1040        y0 = field.rangeMin(Rappture::yaxis);
1041        dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
1042        z0 = field.rangeMin(Rappture::zaxis);
1043        dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
1044        double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
1045       
1046        nx = (int)ceil(dx/dmin);
1047        ny = (int)ceil(dy/dmin);
1048        nz = (int)ceil(dz/dmin);
1049#ifndef NV40
1050        // must be an even power of 2 for older cards
1051        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
1052        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
1053        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
1054#endif
1055        float *data = new float[4*nx*ny*nz];
1056       
1057        double vmin = field.valueMin();
1058        double dv = field.valueMax() - field.valueMin();
1059        if (dv == 0.0) { dv = 1.0; }
1060       
1061        // generate the uniformly sampled data that we need for a volume
1062        int ngen = 0;
1063        double nzero_min = 0.0;
1064        for (iz=0; iz < nz; iz++) {
1065            double zval = z0 + iz*dmin;
1066            for (int iy=0; iy < ny; iy++) {
1067                double yval = y0 + iy*dmin;
1068                for (int ix=0; ix < nx; ix++) {
1069                    double xval = x0 + ix*dmin;
1070                    double v = field.value(xval,yval,zval);
1071                   
1072                    if (v != 0.0f && v < nzero_min) {
1073                        nzero_min = v;
1074                    }
1075                    // scale all values [0-1], -1 => out of bounds
1076                    v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
1077                    data[ngen] = v;
1078                   
1079                    ngen += 4;
1080                }
1081            }
1082        }
1083       
1084        // Compute the gradient of this data.  BE CAREFUL: center
1085        // calculation on each node to avoid skew in either direction.
1086        ngen = 0;
1087        for (int iz=0; iz < nz; iz++) {
1088            for (int iy=0; iy < ny; iy++) {
1089                for (int ix=0; ix < nx; ix++) {
1090                    // gradient in x-direction
1091                    double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
1092                    double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
1093                    if (valm1 < 0 || valp1 < 0) {
1094                        data[ngen+1] = 0.0;
1095                    } else {
1096                        data[ngen+1] = valp1-valm1; // assume dx=1
1097                        //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1 (ISO)
1098                    }
1099                   
1100                    // gradient in y-direction
1101                    valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
1102                    valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
1103                    if (valm1 < 0 || valp1 < 0) {
1104                        data[ngen+2] = 0.0;
1105                    } else {
1106                        data[ngen+2] = valp1-valm1; // assume dy=1
1107                        //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1 (ISO)
1108                    }
1109                   
1110                    // gradient in z-direction
1111                    valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
1112                    valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
1113                    if (valm1 < 0 || valp1 < 0) {
1114                        data[ngen+3] = 0.0;
1115                    } else {
1116                        data[ngen+3] = valp1-valm1; // assume dz=1
1117                        //data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1 (ISO)
1118                    }
1119                   
1120                    ngen += 4;
1121                }
1122            }
1123        }
1124       
1125        Volume *volPtr;
1126        volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
1127                                      field.valueMin(), field.valueMax(), nzero_min);
1128        volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
1129                               field.rangeMax(Rappture::xaxis));
1130        volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
1131                               field.rangeMax(Rappture::yaxis));
1132        volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
1133                               field.rangeMax(Rappture::zaxis));
1134        volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
1135        volPtr->update_pending = true;
1136        // TBD..
1137        // POINTSET
1138        /*
1139          PointSet* pset = new PointSet();
1140          pset->initialize(volume[index], (float*) data);
1141          pset->setVisible(true);
1142          NanoVis::pointSet.push_back(pset);
1143          updateColor(pset);
1144          NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
1145        */
1146       
1147       
1148        delete [] data;
1149    }
1150
1151    //
1152    // Center this new volume on the origin.
1153    //
1154    float dx0 = -0.5;
1155    float dy0 = -0.5*dy/dx;
1156    float dz0 = -0.5*dz/dx;
1157    NanoVis::volume[index]->move(Vector3(dx0, dy0, dz0));
1158    return true;
1159}
1160
1161
1162bool
1163load_volume_stream_insoo(Rappture::Outcome &result, int index,
1164                         std::iostream& fin)
1165{
1166    printf("load_volume_stream\n");
1167 
1168    Rappture::MeshTri2D xymesh;
1169    int dummy, nx, ny, nz, nxy, npts;
1170    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
1171    char line[128], type[128], *start;
1172
1173    int isrect = 1;
1174
1175    dx = dy = dz = 0.0;         // Suppress compiler warning.
1176    x0 = y0 = z0 = 0.0;         // May not have an origin line.
1177    while (!fin.eof()) {
1178        fin.getline(line, sizeof(line) - 1);
1179        fprintf(stderr, "line is %s\n", line);
1180        if (fin.fail()) {
1181            result.addError("line \"%s\"error in data stream");
1182            return false;
1183        }
1184        for (start=line; *start == ' ' || *start == '\t'; start++)
1185            ;  // skip leading blanks
1186
1187        if (*start != '#') {  // skip comment lines
1188            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
1189                // found grid size
1190                isrect = 1;
1191            } else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
1192                isrect = 0;
1193
1194                double xx, yy, zz;
1195                for (int i=0; i < nxy; i++) {
1196                    fin.getline(line,sizeof(line)-1);
1197                    if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
1198                        xymesh.addNode( Rappture::Node2D(xx,yy) );
1199                    }
1200                }
1201
1202                char fpts[128];
1203                sprintf(fpts, "/tmp/tmppts%d", getpid());
1204                char fcells[128];
1205                sprintf(fcells, "/tmp/tmpcells%d", getpid());
1206
1207                std::ofstream ftmp(fpts);
1208                // save corners of bounding box first, to work around meshing
1209                // problems in voronoi utility
1210                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
1211                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
1212                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
1213                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
1214                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
1215                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
1216                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
1217                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
1218                for (int i=0; i < nxy; i++) {
1219                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
1220
1221                }
1222                ftmp.close();
1223
1224                char cmdstr[512];
1225                sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
1226                if (system(cmdstr) == 0) {
1227                    int cx, cy, cz;
1228                    std::ifstream ftri(fcells);
1229                    while (!ftri.eof()) {
1230                        ftri.getline(line,sizeof(line)-1);
1231                        if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
1232                            if (cx >= 4 && cy >= 4 && cz >= 4) {
1233                                // skip first 4 boundary points
1234                                xymesh.addCell(cx-4, cy-4, cz-4);
1235                            }
1236                        }
1237                    }
1238                    ftri.close();
1239                } else {
1240                    result.error("triangularization failed");
1241                    return false;
1242                }
1243                unlink(fpts), unlink(fcells);
1244            } else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
1245                // found z-grid
1246            } else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
1247                // found origin
1248            } else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
1249                // found one of the delta lines
1250                if (ddx != 0.0) { dx = ddx; }
1251                else if (ddy != 0.0) { dy = ddy; }
1252                else if (ddz != 0.0) { dz = ddz; }
1253            } else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
1254                if (isrect && (npts != nx*ny*nz)) {
1255                    result.addError("inconsistent data: expected %d points"
1256                                    " but found %d points", nx*ny*nz, npts);
1257                    return false;
1258                } else if (!isrect && (npts != nxy*nz)) {
1259                    result.addError("inconsistent data: expected %d points"
1260                                    " but found %d points", nx*ny*nz, npts);
1261                    return false;
1262                }
1263                break;
1264            } else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
1265                if (npts != nx*ny*nz) {
1266                    result.addError("inconsistent data: expected %d points"
1267                                    " but found %d points", nx*ny*nz, npts);
1268                    return false;
1269                }
1270                break;
1271            }
1272        }
1273    }
1274
1275    fprintf(stderr, "found nx=%d ny=%d, nz=%d, x0=%f, y0=%f, z0=%f\n",
1276            nx, ny, nz, x0, y0, z0);
1277    // read data points
1278    if (fin.eof()) {
1279        result.error("data not found in stream");
1280        return false;
1281    }
1282
1283    if (isrect) {
1284        Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
1285        Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
1286        Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
1287        Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
1288       
1289        double dval[6];
1290        int nread = 0;
1291        int ix = 0;
1292        int iy = 0;
1293        int iz = 0;
1294        while (!fin.eof() && nread < npts) {
1295            fin.getline(line,sizeof(line)-1);
1296            if (fin.fail()) {
1297                result.error("error reading data points");
1298                return false;
1299            }
1300            int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
1301           
1302            for (int p=0; p < n; p++) {
1303                int nindex = iz*nx*ny + iy*nx + ix;
1304                field.define(nindex, dval[p]);
1305                fprintf(stderr,"nindex = %i\tdval[%i] = %lg\n", nindex, p,
1306                        dval[p]);
1307                fflush(stderr);
1308                nread++;
1309                if (++iz >= nz) {
1310                    iz = 0;
1311                    if (++iy >= ny) {
1312                        iy = 0;
1313                        ++ix;
1314                    }
1315                }
1316            }
1317        }
1318       
1319        // make sure that we read all of the expected points
1320        if (nread != nx*ny*nz) {
1321            result.addError("inconsistent data: expected %d points"
1322                            " but found %d points", nx*ny*nz, npts);
1323            return false;
1324        }
1325       
1326        // figure out a good mesh spacing
1327        int nsample = 30;
1328        dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
1329        dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
1330        dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
1331        double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
1332       
1333        nx = (int)ceil(dx/dmin);
1334        ny = (int)ceil(dy/dmin);
1335        nz = (int)ceil(dz/dmin);
1336       
1337#ifndef NV40
1338        // must be an even power of 2 for older cards
1339        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
1340        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
1341        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
1342#endif
1343       
1344        float *data = new float[4*nx*ny*nz];
1345       
1346        double vmin = field.valueMin();
1347        double dv = field.valueMax() - field.valueMin();
1348        if (dv == 0.0) { dv = 1.0; }
1349       
1350        int ngen = 0;
1351        double nzero_min = 0.0;
1352        for (iz=0; iz < nz; iz++) {
1353            double zval = z0 + iz*dmin;
1354            for (int iy=0; iy < ny; iy++) {
1355                double yval = y0 + iy*dmin;
1356                for (int ix=0; ix < nx; ix++) {
1357                    double xval = x0 + ix*dmin;
1358                    double v = field.value(xval,yval,zval);
1359                   
1360                    if (v != 0.0f && v < nzero_min) {
1361                        nzero_min = v;
1362                    }
1363                    // scale all values [0-1], -1 => out of bounds
1364                    v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
1365                    data[ngen] = v;
1366                   
1367                    ngen += 4;
1368                }
1369            }
1370        }
1371       
1372        // Compute the gradient of this data.  BE CAREFUL: center
1373        // calculation on each node to avoid skew in either direction.
1374        ngen = 0;
1375        for (int iz=0; iz < nz; iz++) {
1376            for (int iy=0; iy < ny; iy++) {
1377                for (int ix=0; ix < nx; ix++) {
1378                    // gradient in x-direction
1379                    double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
1380                    double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
1381                    if (valm1 < 0 || valp1 < 0) {
1382                        data[ngen+1] = 0.0;
1383                    } else {
1384                        data[ngen+1] = valp1-valm1; // assume dx=1
1385                        //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1 (ISO)
1386                    }
1387                   
1388                    // gradient in y-direction
1389                    valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
1390                    valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
1391                    if (valm1 < 0 || valp1 < 0) {
1392                        data[ngen+2] = 0.0;
1393                    } else {
1394                        data[ngen+2] = valp1-valm1; // assume dy=1
1395                        //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1 (ISO)
1396                    }
1397                   
1398                    // gradient in z-direction
1399                    valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
1400                    valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
1401                    if (valm1 < 0 || valp1 < 0) {
1402                        data[ngen+3] = 0.0;
1403                    } else {
1404                        data[ngen+3] = valp1-valm1; // assume dz=1
1405                        //data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1 (ISO)
1406                    }
1407                   
1408                    ngen += 4;
1409                }
1410            }
1411        }
1412       
1413        /*
1414          float *cdata = new float[nx*ny*nz];
1415          int ngen = 0;
1416          double nzero_min = 0.0;
1417          for (int iz=0; iz < nz; iz++) {
1418          double zval = z0 + iz*dmin;
1419          for (int iy=0; iy < ny; iy++) {
1420          double yval = y0 + iy*dmin;
1421          for (int ix=0; ix < nx; ix++) {
1422          double xval = x0 + ix*dmin;
1423          double v = field.value(xval,yval,zval);
1424         
1425          if (v != 0.0f && v < nzero_min) {
1426          nzero_min = v;
1427          }
1428         
1429          // scale all values [0-1], -1 => out of bounds
1430          v = (isnan(v)) ? -1.0 : v;
1431         
1432          cdata[ngen] = v;
1433          ++ngen;
1434          }
1435                }
1436            }
1437           
1438            float* data = computeGradient(cdata, nx, ny, nz, field.valueMin(),
1439                                          field.valueMax());
1440                                         
1441            for (int i=0; i<nx*ny*nz; i++) {
1442                fprintf(stderr,"enddata[%i] = %lg\n",i,data[i]);
1443                fflush(stderr);
1444            }
1445
1446            fprintf(stdout,"End Data Stats index = %i\n",index);
1447            fprintf(stdout,"nx = %i ny = %i nz = %i\n",nx,ny,nz);
1448            fprintf(stdout,"dx = %lg dy = %lg dz = %lg\n",dx,dy,dz);
1449            fprintf(stdout,"dataMin = %lg\tdataMax = %lg\tnzero_min = %lg\n", field.valueMin(),field.valueMax(),nzero_min);
1450            fflush(stdout);
1451            */
1452
1453        Volume *volPtr;
1454        volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
1455                                      field.valueMin(), field.valueMax(), nzero_min);
1456        volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
1457                               field.rangeMax(Rappture::xaxis));
1458        volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
1459                               field.rangeMax(Rappture::yaxis));
1460        volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
1461                               field.rangeMax(Rappture::zaxis));
1462        volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
1463        volPtr->update_pending = true;
1464        // TBD..
1465        // POINTSET
1466        /*
1467          PointSet* pset = new PointSet();
1468          pset->initialize(volume[index], (float*) data);
1469          pset->setVisible(true);
1470          NanoVis::pointSet.push_back(pset);
1471          updateColor(pset);
1472          NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
1473        */
1474       
1475        delete [] data;
1476       
1477    } else {
1478        Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
1479        Rappture::FieldPrism3D field(xymesh, zgrid);
1480       
1481        double dval;
1482        int nread = 0;
1483        int ixy = 0;
1484        int iz = 0;
1485        while (!fin.eof() && nread < npts) {
1486            fin >> dval;
1487            if (fin.fail()) {
1488                char mesg[256];
1489                sprintf(mesg,"after %d of %d points: can't read number",
1490                        nread, npts);
1491                result.error(mesg);
1492                return false;
1493            } else {
1494                int nid = nxy*iz + ixy;
1495                field.define(nid, dval);
1496               
1497                nread++;
1498                if (++iz >= nz) {
1499                    iz = 0;
1500                    ixy++;
1501                }
1502            }
1503        }
1504       
1505        // make sure that we read all of the expected points
1506        if (nread != nxy*nz) {
1507            result.addError("inconsistent data: expected %d points"
1508                            " but found %d points", nx*ny*nz, npts);
1509            return false;
1510        }
1511       
1512        // figure out a good mesh spacing
1513        int nsample = 30;
1514        x0 = field.rangeMin(Rappture::xaxis);
1515        dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
1516        y0 = field.rangeMin(Rappture::yaxis);
1517        dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
1518        z0 = field.rangeMin(Rappture::zaxis);
1519        dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
1520        double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
1521       
1522        nx = (int)ceil(dx/dmin);
1523        ny = (int)ceil(dy/dmin);
1524        nz = (int)ceil(dz/dmin);
1525#ifndef NV40
1526        // must be an even power of 2 for older cards
1527        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
1528        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
1529        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
1530#endif
1531        float *data = new float[4*nx*ny*nz];
1532       
1533        double vmin = field.valueMin();
1534        double dv = field.valueMax() - field.valueMin();
1535        if (dv == 0.0) { dv = 1.0; }
1536       
1537        // generate the uniformly sampled data that we need for a volume
1538        int ngen = 0;
1539        double nzero_min = 0.0;
1540        for (iz=0; iz < nz; iz++) {
1541            double zval = z0 + iz*dmin;
1542            for (int iy=0; iy < ny; iy++) {
1543                double yval = y0 + iy*dmin;
1544                for (int ix=0; ix < nx; ix++) {
1545                    double xval = x0 + ix*dmin;
1546                    double v = field.value(xval,yval,zval);
1547                   
1548                    if (v != 0.0f && v < nzero_min) {
1549                        nzero_min = v;
1550                    }
1551                    // scale all values [0-1], -1 => out of bounds
1552                    v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
1553                    data[ngen] = v;
1554                   
1555                    ngen += 4;
1556                }
1557            }
1558        }
1559       
1560        // Compute the gradient of this data.  BE CAREFUL: center
1561        // calculation on each node to avoid skew in either direction.
1562        ngen = 0;
1563        for (int iz=0; iz < nz; iz++) {
1564            for (int iy=0; iy < ny; iy++) {
1565                for (int ix=0; ix < nx; ix++) {
1566                    // gradient in x-direction
1567                    double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
1568                    double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
1569                    if (valm1 < 0 || valp1 < 0) {
1570                        data[ngen+1] = 0.0;
1571                    } else {
1572                        data[ngen+1] = valp1-valm1; // assume dx=1
1573                        //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1 (ISO)
1574                    }
1575                   
1576                    // gradient in y-direction
1577                    valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
1578                    valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
1579                    if (valm1 < 0 || valp1 < 0) {
1580                        data[ngen+2] = 0.0;
1581                    } else {
1582                        data[ngen+2] = valp1-valm1; // assume dy=1
1583                        //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1 (ISO)
1584                    }
1585                   
1586                    // gradient in z-direction
1587                    valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
1588                    valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
1589                    if (valm1 < 0 || valp1 < 0) {
1590                        data[ngen+3] = 0.0;
1591                    } else {
1592                        data[ngen+3] = valp1-valm1; // assume dz=1
1593                        //data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1 (ISO)
1594                    }
1595                   
1596                    ngen += 4;
1597                }
1598            }
1599        }
1600       
1601        Volume *volPtr;
1602        volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
1603                                      field.valueMin(), field.valueMax(), nzero_min);
1604        volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
1605                               field.rangeMax(Rappture::xaxis));
1606        volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
1607                               field.rangeMax(Rappture::yaxis));
1608        volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
1609                               field.rangeMax(Rappture::zaxis));
1610        volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
1611        volPtr->update_pending = true;
1612        // TBD..
1613        // POINTSET
1614        /*
1615          PointSet* pset = new PointSet();
1616          pset->initialize(volume[index], (float*) data);
1617          pset->setVisible(true);
1618          NanoVis::pointSet.push_back(pset);
1619          updateColor(pset);
1620          NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
1621        */
1622       
1623       
1624        delete [] data;
1625    }
1626
1627    //
1628    // Center this new volume on the origin.
1629    //
1630    float dx0 = -0.5;
1631    float dy0 = -0.5*dy/dx;
1632    float dz0 = -0.5*dz/dx;
1633    NanoVis::volume[index]->move(Vector3(dx0, dy0, dz0));
1634    return true;
1635}
Note: See TracBrowser for help on using the repository browser.