source: nanovis/tags/1.2.3/dxReader.cpp @ 6369

Last change on this file since 6369 was 5588, checked in by ldelgass, 9 years ago

Merge support for VTK vector fields in nanovis to release branch
  • Property svn:eol-style set to native
File size: 17.6 KB
Line 
1 /* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */
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-2013  HUBzero Foundation, LLC
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#include <stdio.h>
21#include <math.h>
22#include <sys/types.h>
23#include <unistd.h>
24#include <float.h>
25
26#include <fstream>
27#include <iostream>
28#include <sstream>
29#include <string>
30
31#include <RpMesh1D.h>
32#include <RpFieldRect3D.h>
33#include <RpFieldPrism3D.h>
34
35// common file/data reader functions
36#include "ReaderCommon.h"
37
38#include "config.h"
39#include "nanovis.h"
40#include "dxReader.h"
41#include "Unirect.h"
42#include "Volume.h"
43
44using namespace nv;
45
46/**
47 * \brief Load a 3D volume from a dx-format file
48 *
49 * In DX format:
50 *  rank 0 means scalars,
51 *  rank 1 means vectors,
52 *  rank 2 means matrices,
53 *  rank 3 means tensors
54 *
55 *  For rank 1, shape is a single number equal to the number of dimensions.
56 *  e.g. rank 1 shape 3 means a 3-component vector field
57 */
58Volume *
59nv::load_dx_volume_stream(const char *tag, std::iostream& fin)
60{
61    TRACE("Enter tag: %s", tag);
62
63    Rappture::MeshTri2D xymesh;
64    int dummy, nx, ny, nz, nxy, npts, naxes;
65    // origin
66    double x0, y0, z0;
67    // deltas (cell dimensions)
68    double dx, dy, dz;
69    // length of volume edges
70    double lx, ly, lz;
71    // temps for delta values
72    double ddx, ddy, ddz;
73    char line[128], type[128], *start;
74
75    int isrect = 1;
76
77    dx = dy = dz = 0.0;
78    x0 = y0 = z0 = 0.0; // May not have an origin line.
79    nx = ny = nz = npts = nxy = naxes = 0;
80    while (!fin.eof()) {
81        fin.getline(line, sizeof(line) - 1);
82        if (fin.fail()) {
83            ERROR("error in data stream");
84            return NULL;
85        }
86        for (start = line; *start == ' ' || *start == '\t'; start++)
87            ;  // skip leading blanks
88
89        if (*start != '#') {  // skip comment lines
90            if (sscanf(start, "object %d class gridpositions counts %d %d %d",
91                       &dummy, &nx, &ny, &nz) == 4) {
92                // found grid size
93                isrect = 1;
94            } else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows",
95                              &dummy, &nxy) == 2) {
96                isrect = 0;
97
98                double xx, yy, zz;
99                for (int i = 0; i < nxy; i++) {
100                    fin.getline(line,sizeof(line)-1);
101                    if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
102                        xymesh.addNode(Rappture::Node2D(xx, yy));
103                    }
104                }
105
106                char fpts[128];
107                sprintf(fpts, "/tmp/tmppts%d", getpid());
108                char fcells[128];
109                sprintf(fcells, "/tmp/tmpcells%d", getpid());
110
111                std::ofstream ftmp(fpts);
112                // save corners of bounding box first, to work around meshing
113                // problems in voronoi utility
114                int numBoundaryPoints = 4;
115
116                // Bump out the bounds by an epsilon to avoid problem when
117                // corner points are already nodes
118                double XEPS = (xymesh.rangeMax(Rappture::xaxis) -
119                               xymesh.rangeMin(Rappture::xaxis)) / 10.0f;
120
121                double YEPS = (xymesh.rangeMax(Rappture::yaxis) -
122                               xymesh.rangeMin(Rappture::yaxis)) / 10.0f;
123
124                ftmp << xymesh.rangeMin(Rappture::xaxis) - XEPS << " "
125                     << xymesh.rangeMin(Rappture::yaxis) - YEPS << std::endl;
126                ftmp << xymesh.rangeMax(Rappture::xaxis) + XEPS << " "
127                     << xymesh.rangeMin(Rappture::yaxis) - YEPS << std::endl;
128                ftmp << xymesh.rangeMax(Rappture::xaxis) + XEPS << " "
129                     << xymesh.rangeMax(Rappture::yaxis) + YEPS << std::endl;
130                ftmp << xymesh.rangeMin(Rappture::xaxis) - XEPS << " "
131                     << xymesh.rangeMax(Rappture::yaxis) + YEPS << std::endl;
132
133                for (int i = 0; i < nxy; i++) {
134                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
135                }
136                ftmp.close();
137
138                char cmdstr[512];
139                sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
140                if (system(cmdstr) == 0) {
141                    int cx, cy, cz;
142                    std::ifstream ftri(fcells);
143                    while (!ftri.eof()) {
144                        ftri.getline(line,sizeof(line)-1);
145                        if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
146                            if (cx >= numBoundaryPoints &&
147                                cy >= numBoundaryPoints &&
148                                cz >= numBoundaryPoints) {
149                                // skip boundary points we added
150                                xymesh.addCell(cx - numBoundaryPoints,
151                                               cy - numBoundaryPoints,
152                                               cz - numBoundaryPoints);
153                            }
154                        }
155                    }
156                    ftri.close();
157                } else {
158                    ERROR("triangularization failed");
159                    return NULL;
160                }
161                unlink(fpts);
162                unlink(fcells);
163            } else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
164                // found z-grid
165            } else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
166                // found origin
167            } else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
168                // found one of the delta lines
169                int count = 0;
170                if (ddx != 0.0) {
171                    dx = ddx;
172                    count++;
173                }
174                if (ddy != 0.0) {
175                    dy = ddy;
176                    count++;
177                }
178                if (ddz != 0.0) {
179                    dz = ddz;
180                    count++;
181                }
182                if (count > 1) {
183                    ERROR("don't know how to handle multiple non-zero"
184                          " delta values");
185                    return NULL;
186                }
187                naxes++;
188            } else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows",
189                              &dummy, type, &npts) == 3) {
190                if (isrect && (npts != nx*ny*nz)) {
191                    ERROR("inconsistent data: expected %d points"
192                          " but found %d points", nx*ny*nz, npts);
193                    return NULL;
194                } else if (!isrect && (npts != nxy*nz)) {
195                    ERROR("inconsistent data: expected %d points"
196                          " but found %d points", nxy*nz, npts);
197                    return NULL;
198                }
199                break;
200            } else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows",
201                              &dummy, type, &npts) == 3) {
202                if (npts != nx*ny*nz) {
203                    ERROR("inconsistent data: expected %d points"
204                          " but found %d points", nx*ny*nz, npts);
205                    return NULL;
206                }
207                break;
208            }
209        }
210    }
211
212    TRACE("found nx=%d ny=%d nxy=%d nz=%d\nnaxes=%d dx=%g dy=%g dz=%g\nx0=%g y0=%g z0=%g",
213          nx, ny, nxy, nz, naxes, dx, dy, dz, x0, y0, z0);
214
215    if (npts > 1 && naxes == 0) {
216        ERROR("Invalid DX file: no deltas found");
217        return NULL;
218    }
219    if (npts > 1 && ((dx == dy) && (dx == dz) && (dx == 0.0))) {
220        ERROR("Invalid deltas in DX file: %g %g %g", dx, dy, dz);
221        return NULL;
222    }
223    if (dx < 0.0 || dy < 0.0 || dz < 0.0) {
224        ERROR("Negative deltas not supported in DX file: %g %g %g", dx, dy, dz);
225        return NULL;
226    }
227
228    lx = (nx - 1) * dx;
229    ly = (ny - 1) * dy;
230    lz = (nz - 1) * dz;
231
232    // read data points
233    if (fin.eof()) {
234        ERROR("data not found in stream");
235        return NULL;
236    }
237    Volume *volume = NULL;
238    float *data = NULL;
239    double vmin = DBL_MAX;
240    double nzero_min = DBL_MAX;
241    double vmax = -DBL_MAX;
242    if (isrect) {
243#ifdef DOWNSAMPLE_DATA
244        Rappture::Mesh1D xgrid(x0, x0 + lx, nx);
245        Rappture::Mesh1D ygrid(y0, y0 + ly, ny);
246        Rappture::Mesh1D zgrid(z0, z0 + lz, nz);
247        Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
248#else // !DOWNSAMPLE_DATA
249        data = new float[nx *  ny *  nz * 4];
250        memset(data, 0, nx*ny*nz*4);
251#endif // DOWNSAMPLE_DATA
252        double dval[6];
253        int nread = 0;
254        int ix = 0;
255        int iy = 0;
256        int iz = 0;
257
258        while (!fin.eof() && nread < npts) {
259            fin.getline(line,sizeof(line)-1);
260            if (fin.fail()) {
261                ERROR("error reading data points");
262                return NULL;
263            }
264            int n = sscanf(line, "%lg %lg %lg %lg %lg %lg",
265                           &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
266            for (int p = 0; p < n; p++) {
267#ifdef notdef
268                if (isnan(dval[p])) {
269                    TRACE("Found NAN in input at %d,%d,%d", ix, iy, iz);
270                }
271#endif
272#ifdef DOWNSAMPLE_DATA
273                int nindex = iz*nx*ny + iy*nx + ix;
274                field.define(nindex, dval[p]);
275#else // !DOWNSAMPLE_DATA
276                int nindex = (iz*nx*ny + iy*nx + ix) * 4;
277                data[nindex] = dval[p];
278
279                if (dval[p] < vmin) {
280                    vmin = dval[p];
281                } else if (dval[p] > vmax) {
282                    vmax = dval[p];
283                }
284                if (dval[p] != 0.0 && dval[p] < nzero_min) {
285                    nzero_min = dval[p];
286                }
287#endif // DOWNSAMPLE_DATA
288                nread++;
289                if (++iz >= nz) {
290                    iz = 0;
291                    if (++iy >= ny) {
292                        iy = 0;
293                        ++ix;
294                    }
295                }
296            }
297        }
298
299        // make sure that we read all of the expected points
300        if (nread != npts) {
301            ERROR("inconsistent data: expected %d points"
302                  " but found %d points", npts, nread);
303            return NULL;
304        }
305
306#ifndef DOWNSAMPLE_DATA
307
308        // scale all values [0-1], -1 => out of bounds
309        normalizeScalar(data, npts, 4, vmin, vmax);
310        computeSimpleGradient(data, nx, ny, nz,
311                              dx, dy, dz);
312#else // DOWNSAMPLE_DATA
313        // figure out a good mesh spacing
314        int nsample = 64;
315        double dmin = pow((lx*ly*lz)/((nsample-1)*(nsample-1)*(nsample-1)), 0.333);
316
317        nx = (int)ceil(lx/dmin) + 1;
318        ny = (int)ceil(ly/dmin) + 1;
319        nz = (int)ceil(lz/dmin) + 1;
320
321#ifndef HAVE_NPOT_TEXTURES
322        // must be an even power of 2 for older cards
323        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
324        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
325        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
326#endif
327
328        dx = lx /(double)(nx - 1);
329        dy = ly /(double)(ny - 1);
330        dz = lz /(double)(nz - 1);
331
332        vmin = field.valueMin();
333        vmax = field.valueMax();
334        nzero_min = DBL_MAX;
335        double dv = vmax - vmin;
336        if (dv == 0.0) {
337            dv = 1.0;
338        }
339
340        int ngen = 0;
341#ifdef FILTER_GRADIENTS
342        // Sobel filter expects a single float at
343        // each node
344        const int step = 1;
345        float *cdata = new float[nx*ny*nz * step];
346#else // !FILTER_GRADIENTS
347        // Leave 3 floats of space for gradients
348        // to be filled in by computeSimpleGradient()
349        const int step = 4;
350        data = new float[nx*ny*nz * step];
351#endif // FILTER_GRADIENTS
352
353        for (int iz = 0; iz < nz; iz++) {
354            double zval = z0 + iz*dz;
355            for (int iy = 0; iy < ny; iy++) {
356                double yval = y0 + iy*dy;
357                for (int ix = 0; ix < nx; ix++) {
358                    double xval = x0 + ix*dx;
359                    double v = field.value(xval, yval, zval);
360#ifdef notdef
361                    if (isnan(v)) {
362                        TRACE("NAN at %d,%d,%d: (%g, %g, %g)", ix, iy, iz, xval, yval, zval);
363                    }
364#endif
365                    if (v != 0.0 && v < nzero_min) {
366                        nzero_min = v;
367                    }
368#ifdef FILTER_GRADIENTS
369                    // NOT normalized, may have NaNs (FIXME: how does gradient filter handle NaN?)
370                    cdata[ngen] = v;
371#else // !FILTER_GRADIENTS
372                    // scale all values [0-1], -1 => out of bounds
373                    v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
374                    data[ngen] = v;
375#endif // FILTER_GRADIENTS
376                    ngen += step;
377                }
378            }
379        }
380#ifdef FILTER_GRADIENTS
381        // computeGradient returns a new array with gradients
382        // filled in, so data is now 4 floats per node
383        data = computeGradient(cdata, nx, ny, nz,
384                               dx, dy, dz,
385                               vmin, vmax);
386        delete [] cdata;
387#else // !FILTER_GRADIENTS
388        computeSimpleGradient(data, nx, ny, nz,
389                              dx, dy, dz);
390#endif // FILTER_GRADIENTS
391
392#endif // DOWNSAMPLE_DATA
393    } else {
394        Rappture::Mesh1D zgrid(z0, z0 + (nz-1)*dz, nz);
395        Rappture::FieldPrism3D field(xymesh, zgrid);
396
397        double dval;
398        int nread = 0;
399        int ixy = 0;
400        int iz = 0;
401        while (!fin.eof() && nread < npts) {
402            fin >> dval;
403            if (fin.fail()) {
404                ERROR("after %d of %d points: can't read number",
405                      nread, npts);
406                return NULL;
407            } else {
408                int nid = nxy*iz + ixy;
409                field.define(nid, dval);
410
411                nread++;
412                if (++iz >= nz) {
413                    iz = 0;
414                    ixy++;
415                }
416            }
417        }
418
419        // make sure that we read all of the expected points
420        if (nread != npts) {
421            ERROR("inconsistent data: expected %d points"
422                  " but found %d points", npts, nread);
423            return NULL;
424        }
425
426        x0 = field.rangeMin(Rappture::xaxis);
427        lx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
428        y0 = field.rangeMin(Rappture::yaxis);
429        ly = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
430        z0 = field.rangeMin(Rappture::zaxis);
431        lz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
432
433        // figure out a good mesh spacing
434        int nsample = 64;
435        double dmin = pow((lx*ly*lz)/((nsample-1)*(nsample-1)*(nsample-1)), 0.333);
436
437        nx = (int)ceil(lx/dmin) + 1;
438        ny = (int)ceil(ly/dmin) + 1;
439        nz = (int)ceil(lz/dmin) + 1;
440#ifndef HAVE_NPOT_TEXTURES
441        // must be an even power of 2 for older cards
442        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
443        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
444        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
445#endif
446
447        dx = lx /(double)(nx - 1);
448        dy = ly /(double)(ny - 1);
449        dz = lz /(double)(nz - 1);
450
451        vmin = field.valueMin();
452        vmax = field.valueMax();
453        nzero_min = DBL_MAX;
454        double dv = vmax - vmin;
455        if (dv == 0.0) {
456            dv = 1.0;
457        }
458
459        data = new float[4*nx*ny*nz];
460        // generate the uniformly sampled data that we need for a volume
461        int ngen = 0;
462        for (int iz = 0; iz < nz; iz++) {
463            double zval = z0 + iz*dz;
464            for (int iy = 0; iy < ny; iy++) {
465                double yval = y0 + iy*dy;
466                for (int ix = 0; ix < nx; ix++) {
467                    double xval = x0 + ix*dx;
468                    double v = field.value(xval, yval, zval);
469
470                    if (v != 0.0 && v < nzero_min) {
471                        nzero_min = v;
472                    }
473                    // scale all values [0-1], -1 => out of bounds
474                    v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
475                    data[ngen] = v;
476
477                    ngen += 4;
478                }
479            }
480        }
481
482        computeSimpleGradient(data, nx, ny, nz,
483                              dx, dy, dz);
484    }
485
486    TRACE("nx = %i ny = %i nz = %i", nx, ny, nz);
487    TRACE("x0 = %g y0 = %g z0 = %g", x0, y0, z0);
488    TRACE("lx = %g ly = %g lz = %g", lx, ly, lz);
489    TRACE("dx = %g dy = %g dz = %g", dx, dy, dz);
490    TRACE("dataMin = %g dataMax = %g nzero_min = %g",
491          vmin, vmax, nzero_min);
492
493    volume = NanoVis::loadVolume(tag, nx, ny, nz, 4, data,
494                                 vmin, vmax, nzero_min);
495    volume->xAxis.setRange(x0, x0 + lx);
496    volume->yAxis.setRange(y0, y0 + ly);
497    volume->zAxis.setRange(z0, z0 + lz);
498    volume->updatePending = true;
499
500    delete [] data;
501
502    //
503    // Center this new volume on the origin.
504    //
505    float dx0 = -0.5;
506    float dy0 = -0.5*ly/lx;
507    float dz0 = -0.5*lz/lx;
508    if (volume) {
509        volume->setPosition(vrmath::Vector3f(dx0, dy0, dz0));
510        TRACE("Set volume position to %g %g %g", dx0, dy0, dz0);
511    }
512    return volume;
513}
Note: See TracBrowser for help on using the repository browser.