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

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