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

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