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

Last change on this file since 1053 was 1053, checked in by gah, 16 years ago

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