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

Last change on this file since 1321 was 1312, checked in by dkearney, 16 years ago

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