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

Last change on this file since 3874 was 3874, checked in by ldelgass, 6 years ago

fix trace output

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