source: nanovis/trunk/dxReader.cpp @ 6369

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

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