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

Last change on this file since 1049 was 1000, checked in by vrinside, 16 years ago
File size: 58.0 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            double v;
789            int ngen = 0;
790            if (dv == 0.0)  dv = 1.0;
791
792            for (int iz=0; iz < nz; iz++) {
793                double zval = z0 + iz*dmin;
794                for (int iy=0; iy < ny; iy++) {
795                    double yval = y0 + iy*dmin;
796                    for (int ix=0; ix < nx; ix++) {
797                        double xval = x0 + ix*dmin;
798                        double v = field.value(xval,yval,zval);
799
800                        // scale all values [0-1], -1 => out of bounds
801                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
802
803                        data[ngen] = v;
804                        ngen += 4;
805                    }
806                }
807            }
808
809            computeSimpleGradient(data, nx, ny, nz);
810#endif
811
812            for (int i=0; i<nx*ny*nz; i++) {
813                fprintf(stderr,"enddata[%i] = %lg\n",i,data[i]);
814                fflush(stderr);
815            }
816
817            fprintf(stdout,"End Data Stats index = %i\n",index);
818            fprintf(stdout,"nx = %i ny = %i nz = %i\n",nx,ny,nz);
819            fprintf(stdout,"dx = %lg dy = %lg dz = %lg\n",dx,dy,dz);
820            fprintf(stdout,"dataMin = %lg\tdataMax = %lg\tnzero_min = %lg\n", field.valueMin(),field.valueMax(),nzero_min);
821            fflush(stdout);
822
823            Volume *volPtr;
824            volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
825            field.valueMin(), field.valueMax(), nzero_min);
826            volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
827                   field.rangeMax(Rappture::xaxis));
828            volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
829                   field.rangeMax(Rappture::yaxis));
830            volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
831                   field.rangeMax(Rappture::zaxis));
832            volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
833            volPtr->update_pending = true;
834            // TBD..
835            // POINTSET
836            /*
837              PointSet* pset = new PointSet();
838              pset->initialize(volume[index], (float*) data);
839              pset->setVisible(true);
840              NanoVis::pointSet.push_back(pset);
841              updateColor(pset);
842              NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
843            */
844
845            delete [] data;
846
847        } else {
848            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
849            Rappture::FieldPrism3D field(xymesh, zgrid);
850
851            double dval;
852            int nread = 0;
853            int ixy = 0;
854            int iz = 0;
855            while (!fin.eof() && nread < npts) {
856                fin >> dval;
857                if (fin.fail()) {
858                    char mesg[256];
859                    sprintf(mesg,"after %d of %d points: can't read number",
860                            nread, npts);
861                    return result.error(mesg);
862                } else {
863                    int nid = nxy*iz + ixy;
864                    field.define(nid, dval);
865
866                    nread++;
867                    if (++iz >= nz) {
868                        iz = 0;
869                        ixy++;
870                    }
871                }
872            }
873
874            // make sure that we read all of the expected points
875            if (nread != nxy*nz) {
876                char mesg[256];
877                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, nread);
878                return result.error(mesg);
879            }
880
881            // figure out a good mesh spacing
882            int nsample = 30;
883            x0 = field.rangeMin(Rappture::xaxis);
884            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
885            y0 = field.rangeMin(Rappture::yaxis);
886            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
887            z0 = field.rangeMin(Rappture::zaxis);
888            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
889            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
890
891            nx = (int)ceil(dx/dmin);
892            ny = (int)ceil(dy/dmin);
893            nz = (int)ceil(dz/dmin);
894#ifndef NV40
895            // must be an even power of 2 for older cards
896            nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
897            ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
898            nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
899#endif
900            float *data = new float[4*nx*ny*nz];
901
902            double vmin = field.valueMin();
903            double dv = field.valueMax() - field.valueMin();
904            if (dv == 0.0) { dv = 1.0; }
905
906            // generate the uniformly sampled data that we need for a volume
907            int ngen = 0;
908            double nzero_min = 0.0;
909            for (iz=0; iz < nz; iz++) {
910                double zval = z0 + iz*dmin;
911                for (int iy=0; iy < ny; iy++) {
912                    double yval = y0 + iy*dmin;
913                    for (int ix=0; ix < nx; ix++) {
914                        double xval = x0 + ix*dmin;
915                        double v = field.value(xval,yval,zval);
916
917                        if (v != 0.0f && v < nzero_min)
918                            {
919                                nzero_min = v;
920                            }
921                        // scale all values [0-1], -1 => out of bounds
922                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
923                        data[ngen] = v;
924
925                        ngen += 4;
926                    }
927                }
928            }
929
930            // Compute the gradient of this data.  BE CAREFUL: center
931            // calculation on each node to avoid skew in either direction.
932            ngen = 0;
933            for (int iz=0; iz < nz; iz++) {
934                for (int iy=0; iy < ny; iy++) {
935                    for (int ix=0; ix < nx; ix++) {
936                        // gradient in x-direction
937                        double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
938                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
939                        if (valm1 < 0 || valp1 < 0) {
940                            data[ngen+1] = 0.0;
941                        } else {
942                            data[ngen+1] = valp1-valm1; // assume dx=1
943                            //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1 (ISO)
944                        }
945
946                        // gradient in y-direction
947                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
948                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
949                        if (valm1 < 0 || valp1 < 0) {
950                            data[ngen+2] = 0.0;
951                        } else {
952                            data[ngen+2] = valp1-valm1; // assume dy=1
953                            //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1 (ISO)
954                        }
955
956                        // gradient in z-direction
957                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
958                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
959                        if (valm1 < 0 || valp1 < 0) {
960                            data[ngen+3] = 0.0;
961                        } else {
962                            data[ngen+3] = valp1-valm1; // assume dz=1
963                            //data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1 (ISO)
964                        }
965
966                        ngen += 4;
967                    }
968                }
969            }
970
971            Volume *volPtr;
972            volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
973            field.valueMin(), field.valueMax(), nzero_min);
974            volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
975                   field.rangeMax(Rappture::xaxis));
976            volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
977                   field.rangeMax(Rappture::yaxis));
978            volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
979                   field.rangeMax(Rappture::zaxis));
980            volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
981            volPtr->update_pending = true;
982            // TBD..
983            // POINTSET
984            /*
985              PointSet* pset = new PointSet();
986              pset->initialize(volume[index], (float*) data);
987              pset->setVisible(true);
988              NanoVis::pointSet.push_back(pset);
989              updateColor(pset);
990              NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
991            */
992
993
994            delete [] data;
995        }
996    } else {
997        return result.error("data not found in stream");
998    }
999
1000    //
1001    // Center this new volume on the origin.
1002    //
1003    float dx0 = -0.5;
1004    float dy0 = -0.5*dy/dx;
1005    float dz0 = -0.5*dz/dx;
1006    NanoVis::volume[index]->move(Vector3(dx0, dy0, dz0));
1007    return result;
1008}
1009
1010Rappture::Outcome
1011load_volume_stream_insoo(int index, std::iostream& fin)
1012{
1013    printf("load_volume_stream\n");
1014    Rappture::Outcome result;
1015
1016    Rappture::MeshTri2D xymesh;
1017    int dummy, nx, ny, nz, nxy, npts;
1018    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
1019    char line[128], type[128], *start;
1020
1021    int isrect = 1;
1022
1023    dx = dy = dz = 0.0;         // Suppress compiler warning.
1024    x0 = y0 = z0 = 0.0;         // May not have an origin line.
1025    while (!fin.eof()) {
1026        fin.getline(line, sizeof(line) - 1);
1027        if (fin.fail()) {
1028            return result.error("error in data stream");
1029        }
1030        for (start=line; *start == ' ' || *start == '\t'; start++)
1031            ;  // skip leading blanks
1032
1033        if (*start != '#') {  // skip comment lines
1034            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
1035                // found grid size
1036                isrect = 1;
1037            } else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
1038                isrect = 0;
1039
1040                double xx, yy, zz;
1041                for (int i=0; i < nxy; i++) {
1042                    fin.getline(line,sizeof(line)-1);
1043                    if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
1044                        xymesh.addNode( Rappture::Node2D(xx,yy) );
1045                    }
1046                }
1047
1048                char fpts[128];
1049                sprintf(fpts, "/tmp/tmppts%d", getpid());
1050                char fcells[128];
1051                sprintf(fcells, "/tmp/tmpcells%d", getpid());
1052
1053                std::ofstream ftmp(fpts);
1054                // save corners of bounding box first, to work around meshing
1055                // problems in voronoi utility
1056                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
1057                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
1058                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
1059                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
1060                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
1061                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
1062                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
1063                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
1064                for (int i=0; i < nxy; i++) {
1065                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
1066
1067                }
1068                ftmp.close();
1069
1070                char cmdstr[512];
1071                sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
1072                if (system(cmdstr) == 0) {
1073                    int cx, cy, cz;
1074                    std::ifstream ftri(fcells);
1075                    while (!ftri.eof()) {
1076                        ftri.getline(line,sizeof(line)-1);
1077                        if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
1078                            if (cx >= 4 && cy >= 4 && cz >= 4) {
1079                                // skip first 4 boundary points
1080                                xymesh.addCell(cx-4, cy-4, cz-4);
1081                            }
1082                        }
1083                    }
1084                    ftri.close();
1085                } else {
1086                    return result.error("triangularization failed");
1087                }
1088
1089                sprintf(cmdstr, "rm -f %s %s", fpts, fcells);
1090                system(cmdstr);
1091            } else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
1092                // found z-grid
1093            } else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
1094                // found origin
1095            } else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
1096                // found one of the delta lines
1097                if (ddx != 0.0) { dx = ddx; }
1098                else if (ddy != 0.0) { dy = ddy; }
1099                else if (ddz != 0.0) { dz = ddz; }
1100            } else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
1101                if (isrect && (npts != nx*ny*nz)) {
1102                    char mesg[256];
1103                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
1104                    return result.error(mesg);
1105                } else if (!isrect && (npts != nxy*nz)) {
1106                    char mesg[256];
1107                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, npts);
1108                    return result.error(mesg);
1109                }
1110                break;
1111            } else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
1112                if (npts != nx*ny*nz) {
1113                    char mesg[256];
1114                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
1115                    return result.error(mesg);
1116                }
1117                break;
1118            }
1119        }
1120    }
1121
1122    // read data points
1123    if (!fin.eof()) {
1124        if (isrect) {
1125            Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
1126            Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
1127            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
1128            Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
1129
1130            double dval[6];
1131            int nread = 0;
1132            int ix = 0;
1133            int iy = 0;
1134            int iz = 0;
1135            while (!fin.eof() && nread < npts) {
1136                fin.getline(line,sizeof(line)-1);
1137                if (fin.fail()) {
1138                    return result.error("error reading data points");
1139                }
1140                int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
1141
1142                for (int p=0; p < n; p++) {
1143                    int nindex = iz*nx*ny + iy*nx + ix;
1144                    field.define(nindex, dval[p]);
1145                    fprintf(stderr,"nindex = %i\tdval[%i] = %lg\n", nindex, p,
1146                        dval[p]);
1147                    fflush(stderr);
1148                    nread++;
1149                    if (++iz >= nz) {
1150                        iz = 0;
1151                        if (++iy >= ny) {
1152                            iy = 0;
1153                            ++ix;
1154                        }
1155                    }
1156                }
1157            }
1158
1159            // make sure that we read all of the expected points
1160            if (nread != nx*ny*nz) {
1161                char mesg[256];
1162                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, nread);
1163                result.error(mesg);
1164                return result;
1165            }
1166
1167            // figure out a good mesh spacing
1168            int nsample = 30;
1169            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
1170            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
1171            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
1172            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
1173
1174            nx = (int)ceil(dx/dmin);
1175            ny = (int)ceil(dy/dmin);
1176            nz = (int)ceil(dz/dmin);
1177
1178#ifndef NV40
1179            // must be an even power of 2 for older cards
1180            nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
1181            ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
1182            nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
1183#endif
1184
1185            float *data = new float[4*nx*ny*nz];
1186
1187            double vmin = field.valueMin();
1188            double dv = field.valueMax() - field.valueMin();
1189            if (dv == 0.0) { dv = 1.0; }
1190
1191            int ngen = 0;
1192            double nzero_min = 0.0;
1193            for (iz=0; iz < nz; iz++) {
1194                double zval = z0 + iz*dmin;
1195                for (int iy=0; iy < ny; iy++) {
1196                    double yval = y0 + iy*dmin;
1197                    for (int ix=0; ix < nx; ix++) {
1198                        double xval = x0 + ix*dmin;
1199                        double v = field.value(xval,yval,zval);
1200
1201                        if (v != 0.0f && v < nzero_min)
1202                            {
1203                                nzero_min = v;
1204                            }
1205                        // scale all values [0-1], -1 => out of bounds
1206                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
1207                        data[ngen] = v;
1208
1209                        ngen += 4;
1210                    }
1211                }
1212            }
1213
1214            // Compute the gradient of this data.  BE CAREFUL: center
1215            // calculation on each node to avoid skew in either direction.
1216            ngen = 0;
1217            for (int iz=0; iz < nz; iz++) {
1218                for (int iy=0; iy < ny; iy++) {
1219                    for (int ix=0; ix < nx; ix++) {
1220                        // gradient in x-direction
1221                        double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
1222                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
1223                        if (valm1 < 0 || valp1 < 0) {
1224                            data[ngen+1] = 0.0;
1225                        } else {
1226                            data[ngen+1] = valp1-valm1; // assume dx=1
1227                            //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1 (ISO)
1228                        }
1229
1230                        // gradient in y-direction
1231                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
1232                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
1233                        if (valm1 < 0 || valp1 < 0) {
1234                            data[ngen+2] = 0.0;
1235                        } else {
1236                            data[ngen+2] = valp1-valm1; // assume dy=1
1237                            //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1 (ISO)
1238                        }
1239
1240                        // gradient in z-direction
1241                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
1242                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
1243                        if (valm1 < 0 || valp1 < 0) {
1244                            data[ngen+3] = 0.0;
1245                        } else {
1246                            data[ngen+3] = valp1-valm1; // assume dz=1
1247                            //data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1 (ISO)
1248                        }
1249
1250                        ngen += 4;
1251                    }
1252                }
1253            }
1254
1255/*
1256            float *cdata = new float[nx*ny*nz];
1257            int ngen = 0;
1258            double nzero_min = 0.0;
1259            for (int iz=0; iz < nz; iz++) {
1260                double zval = z0 + iz*dmin;
1261                for (int iy=0; iy < ny; iy++) {
1262                    double yval = y0 + iy*dmin;
1263                    for (int ix=0; ix < nx; ix++) {
1264                        double xval = x0 + ix*dmin;
1265                        double v = field.value(xval,yval,zval);
1266
1267                        if (v != 0.0f && v < nzero_min) {
1268                            nzero_min = v;
1269                        }
1270
1271                        // scale all values [0-1], -1 => out of bounds
1272                        v = (isnan(v)) ? -1.0 : v;
1273
1274                        cdata[ngen] = v;
1275                        ++ngen;
1276                    }
1277                }
1278            }
1279
1280            float* data = computeGradient(cdata, nx, ny, nz, field.valueMin(),
1281                                          field.valueMax());
1282
1283            for (int i=0; i<nx*ny*nz; i++) {
1284                fprintf(stderr,"enddata[%i] = %lg\n",i,data[i]);
1285                fflush(stderr);
1286            }
1287
1288            fprintf(stdout,"End Data Stats index = %i\n",index);
1289            fprintf(stdout,"nx = %i ny = %i nz = %i\n",nx,ny,nz);
1290            fprintf(stdout,"dx = %lg dy = %lg dz = %lg\n",dx,dy,dz);
1291            fprintf(stdout,"dataMin = %lg\tdataMax = %lg\tnzero_min = %lg\n", field.valueMin(),field.valueMax(),nzero_min);
1292            fflush(stdout);
1293            */
1294
1295            Volume *volPtr;
1296            volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
1297            field.valueMin(), field.valueMax(), nzero_min);
1298            volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
1299                   field.rangeMax(Rappture::xaxis));
1300            volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
1301                   field.rangeMax(Rappture::yaxis));
1302            volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
1303                   field.rangeMax(Rappture::zaxis));
1304            volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
1305            volPtr->update_pending = true;
1306            // TBD..
1307            // POINTSET
1308            /*
1309              PointSet* pset = new PointSet();
1310              pset->initialize(volume[index], (float*) data);
1311              pset->setVisible(true);
1312              NanoVis::pointSet.push_back(pset);
1313              updateColor(pset);
1314              NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
1315            */
1316
1317            delete [] data;
1318
1319        } else {
1320            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
1321            Rappture::FieldPrism3D field(xymesh, zgrid);
1322
1323            double dval;
1324            int nread = 0;
1325            int ixy = 0;
1326            int iz = 0;
1327            while (!fin.eof() && nread < npts) {
1328                fin >> dval;
1329                if (fin.fail()) {
1330                    char mesg[256];
1331                    sprintf(mesg,"after %d of %d points: can't read number",
1332                            nread, npts);
1333                    return result.error(mesg);
1334                } else {
1335                    int nid = nxy*iz + ixy;
1336                    field.define(nid, dval);
1337
1338                    nread++;
1339                    if (++iz >= nz) {
1340                        iz = 0;
1341                        ixy++;
1342                    }
1343                }
1344            }
1345
1346            // make sure that we read all of the expected points
1347            if (nread != nxy*nz) {
1348                char mesg[256];
1349                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, nread);
1350                return result.error(mesg);
1351            }
1352
1353            // figure out a good mesh spacing
1354            int nsample = 30;
1355            x0 = field.rangeMin(Rappture::xaxis);
1356            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
1357            y0 = field.rangeMin(Rappture::yaxis);
1358            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
1359            z0 = field.rangeMin(Rappture::zaxis);
1360            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
1361            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
1362
1363            nx = (int)ceil(dx/dmin);
1364            ny = (int)ceil(dy/dmin);
1365            nz = (int)ceil(dz/dmin);
1366#ifndef NV40
1367            // must be an even power of 2 for older cards
1368            nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
1369            ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
1370            nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
1371#endif
1372            float *data = new float[4*nx*ny*nz];
1373
1374            double vmin = field.valueMin();
1375            double dv = field.valueMax() - field.valueMin();
1376            if (dv == 0.0) { dv = 1.0; }
1377
1378            // generate the uniformly sampled data that we need for a volume
1379            int ngen = 0;
1380            double nzero_min = 0.0;
1381            for (iz=0; iz < nz; iz++) {
1382                double zval = z0 + iz*dmin;
1383                for (int iy=0; iy < ny; iy++) {
1384                    double yval = y0 + iy*dmin;
1385                    for (int ix=0; ix < nx; ix++) {
1386                        double xval = x0 + ix*dmin;
1387                        double v = field.value(xval,yval,zval);
1388
1389                        if (v != 0.0f && v < nzero_min)
1390                            {
1391                                nzero_min = v;
1392                            }
1393                        // scale all values [0-1], -1 => out of bounds
1394                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
1395                        data[ngen] = v;
1396
1397                        ngen += 4;
1398                    }
1399                }
1400            }
1401
1402            // Compute the gradient of this data.  BE CAREFUL: center
1403            // calculation on each node to avoid skew in either direction.
1404            ngen = 0;
1405            for (int iz=0; iz < nz; iz++) {
1406                for (int iy=0; iy < ny; iy++) {
1407                    for (int ix=0; ix < nx; ix++) {
1408                        // gradient in x-direction
1409                        double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
1410                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
1411                        if (valm1 < 0 || valp1 < 0) {
1412                            data[ngen+1] = 0.0;
1413                        } else {
1414                            data[ngen+1] = valp1-valm1; // assume dx=1
1415                            //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1 (ISO)
1416                        }
1417
1418                        // gradient in y-direction
1419                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
1420                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
1421                        if (valm1 < 0 || valp1 < 0) {
1422                            data[ngen+2] = 0.0;
1423                        } else {
1424                            data[ngen+2] = valp1-valm1; // assume dy=1
1425                            //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1 (ISO)
1426                        }
1427
1428                        // gradient in z-direction
1429                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
1430                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
1431                        if (valm1 < 0 || valp1 < 0) {
1432                            data[ngen+3] = 0.0;
1433                        } else {
1434                            data[ngen+3] = valp1-valm1; // assume dz=1
1435                            //data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1 (ISO)
1436                        }
1437
1438                        ngen += 4;
1439                    }
1440                }
1441            }
1442
1443            Volume *volPtr;
1444            volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
1445            field.valueMin(), field.valueMax(), nzero_min);
1446            volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
1447                   field.rangeMax(Rappture::xaxis));
1448            volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
1449                   field.rangeMax(Rappture::yaxis));
1450            volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
1451                   field.rangeMax(Rappture::zaxis));
1452            volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
1453            volPtr->update_pending = true;
1454            // TBD..
1455            // POINTSET
1456            /*
1457              PointSet* pset = new PointSet();
1458              pset->initialize(volume[index], (float*) data);
1459              pset->setVisible(true);
1460              NanoVis::pointSet.push_back(pset);
1461              updateColor(pset);
1462              NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
1463            */
1464
1465
1466            delete [] data;
1467        }
1468    } else {
1469        return result.error("data not found in stream");
1470    }
1471
1472    //
1473    // Center this new volume on the origin.
1474    //
1475    float dx0 = -0.5;
1476    float dy0 = -0.5*dy/dx;
1477    float dz0 = -0.5*dz/dx;
1478    NanoVis::volume[index]->move(Vector3(dx0, dy0, dz0));
1479    return result;
1480}
Note: See TracBrowser for help on using the repository browser.