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

Last change on this file since 3492 was 3492, checked in by ldelgass, 11 years ago

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