source: nanovis/branches/1.1/dxReader.cpp @ 4820

Last change on this file since 4820 was 4820, checked in by ldelgass, 8 years ago

remove dead code

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