source: nanovis/trunk/dxReader.cpp @ 5399

Last change on this file since 5399 was 5326, checked in by ldelgass, 9 years ago

nanovis can't handle negative deltas in DX files, so return an error if one is
found. I'm not even sure this is valid in the DX format, but there is a tool
generating negative deltas.
  • Property svn:eol-style set to native
File size: 17.3 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 = %lg y0 = %lg z0 = %lg", x0, y0, z0);
488    TRACE("lx = %lg ly = %lg lz = %lg", lx, ly, lz);
489    TRACE("dx = %lg dy = %lg dz = %lg", dx, dy, dz);
490    TRACE("dataMin = %lg dataMax = %lg nzero_min = %lg",
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    return volume;
503}
Note: See TracBrowser for help on using the repository browser.