source: nanovis/trunk/VtkReader.cpp @ 5064

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

Remove some more Rappture deps from nanovis
  • Property svn:eol-style set to native
File size: 14.6 KB
Line 
1/* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */
2/*
3 * Copyright (C) 2004-2013  HUBzero Foundation, LLC
4 *
5 */
6
7#include <iostream>
8#include <float.h>
9
10#include <RpField1D.h>
11#include <RpFieldRect3D.h>
12
13#include <vrmath/Vector3f.h>
14
15#include "ReaderCommon.h"
16
17#include "config.h"
18#include "nanovis.h"
19#include "VtkReader.h"
20#include "Volume.h"
21#include "Trace.h"
22
23using namespace nv;
24
25enum FieldType {
26    FLOAT,
27    DOUBLE,
28    UCHAR,
29    SHORT,
30    USHORT,
31    INT,
32    UINT
33};
34
35static inline int typeSize(FieldType type)
36{
37    switch (type) {
38    case FLOAT:
39        return sizeof(float);
40    case DOUBLE:
41        return sizeof(double);
42    case UCHAR:
43        return sizeof(unsigned char);
44    case SHORT:
45        return sizeof(short);
46    case USHORT:
47        return sizeof(unsigned short);
48    case INT:
49        return sizeof(int);
50    case UINT:
51        return sizeof(unsigned int);
52    default:
53        return 0;
54    }
55}
56
57static inline void allocFieldData(void **data, int npts, FieldType type)
58{
59    switch (type) {
60    case FLOAT:
61        *data = (float *)(new float[npts]);
62        break;
63    case DOUBLE:
64        *data = (double *)(new double[npts]);
65        break;
66    case UCHAR:
67        *data = (unsigned char *)(new unsigned char[npts]);
68        break;
69    case SHORT:
70        *data = (short *)(new short[npts]);
71        break;
72    case USHORT:
73        *data = (unsigned short *)(new unsigned short[npts]);
74        break;
75    case INT:
76        *data = (int *)(new int[npts]);
77        break;
78    case UINT:
79        *data = (unsigned int *)(new unsigned int[npts]);
80        break;
81    default:
82        ;
83    }
84}
85
86static inline void readFieldValue(void *data, int idx, FieldType type, std::iostream& fin)
87{
88    switch (type) {
89    case FLOAT: {
90        float val;
91        fin >> val;
92        ((float *)data)[idx] = val;
93    }
94        break;
95    case DOUBLE: {
96        double val;
97        fin >> val;
98        ((double *)data)[idx] = val;
99    }
100        break;
101    case SHORT: {
102        short val;
103        fin >> val;
104        ((short *)data)[idx] = val;
105    }
106        break;
107    case USHORT: {
108        unsigned short val;
109        fin >> val;
110        ((unsigned short *)data)[idx] = val;
111    }
112        break;
113    case INT: {
114        int val;
115        fin >> val;
116        ((int *)data)[idx] = val;
117    }
118        break;
119    case UINT: {
120        unsigned int val;
121        fin >> val;
122        ((unsigned int *)data)[idx] = val;
123    }
124        break;
125    default:
126        break;
127    }
128
129}
130
131static inline float getFieldData(void *data, int i, FieldType type)
132{
133    switch (type) {
134    case FLOAT:
135        return ((float *)data)[i];
136    case DOUBLE:
137        return (float)((double *)data)[i];
138    case UCHAR:
139        return ((float)((unsigned char *)data)[i]);
140    case SHORT:
141        return (float)((short *)data)[i];
142    case USHORT:
143        return (float)((unsigned short *)data)[i];
144    case INT:
145        return (float)((int *)data)[i];
146    case UINT:
147        return (float)((unsigned int *)data)[i];
148    default:
149        return 0.0f;
150    }
151}
152
153static inline void deleteFieldData(void *data, FieldType type)
154{
155    if (data == NULL)
156        return;
157    switch (type) {
158    case FLOAT:
159        delete [] ((float *)data);
160        break;
161    case DOUBLE:
162        delete [] ((double *)data);
163        break;
164    case UCHAR:
165        delete [] ((unsigned char *)data);
166        break;
167    case SHORT:
168        delete [] ((short *)data);
169        break;
170    case USHORT:
171        delete [] ((unsigned short *)data);
172        break;
173    case INT:
174        delete [] ((int *)data);
175        break;
176    case UINT:
177        delete [] ((unsigned int *)data);
178        break;
179    default:
180        ;
181    }
182}
183
184Volume *
185nv::load_vtk_volume_stream(const char *tag, std::iostream& fin)
186{
187    TRACE("Enter tag:%s", tag);
188
189    int nx, ny, nz, npts;
190    // origin
191    double x0, y0, z0;
192    // deltas (cell dimensions)
193    double dx, dy, dz;
194    // length of volume edges
195    double lx, ly, lz;
196    char line[128], str[128], type[128], *start;
197    int isBinary = -1;
198    int isRect = -1;
199    FieldType ftype = DOUBLE;
200    void *fieldData = NULL;
201    int numRead = 0;
202
203    dx = dy = dz = 0.0;
204    x0 = y0 = z0 = 0.0;
205    nx = ny = nz = npts = 0;
206    while (!fin.eof()) {
207        fin.getline(line, sizeof(line) - 1);
208        if (fin.fail()) {
209            ERROR("Error in data stream");
210            return NULL;
211        }
212        for (start = line; *start == ' ' || *start == '\t'; start++)
213            ;  // skip leading blanks
214
215        if (*start != '#') {  // skip comment lines
216            if (strncmp(start, "BINARY", 6) == 0) {
217                isBinary = 1;
218            } else if (strncmp(start, "ASCII", 5) == 0) {
219                isBinary = 0;
220            } else if (sscanf(start, "DATASET %s", str) == 1) {
221                if (strncmp(str, "STRUCTURED_POINTS", 17) == 0) {
222                    isRect = 1;
223                } else {
224                    ERROR("Unsupported DataSet type '%s'", str);
225                    return NULL;
226                }
227            } else if (sscanf(start, "DIMENSIONS %d %d %d", &nx, &ny, &nz) == 3) {
228                if (nx <= 0 || ny <= 0 || nz <= 0) {
229                    ERROR("Found non-positive dimensions");
230                    return NULL;
231                }
232                // found dimensions
233            } else if (sscanf(start, "ORIGIN %lg %lg %lg", &x0, &y0, &z0) == 3) {
234                // found origin
235            } else if (sscanf(start, "SPACING %lg %lg %lg", &dx, &dy, &dz) == 3) {
236                if ((nx > 1 && dx == 0.0) || (ny > 1 && dy == 0.0) || (nz > 1 && dz == 0.0)) {
237                    ERROR("Found zero spacing for dimension");
238                    return NULL;
239                }
240                // found cell spacing
241            } else if (sscanf(start, "POINT_DATA %d", &npts) == 1) {
242                if (npts < 1 || npts != nx * ny * nz) {
243                    ERROR("Error in data stream: unexpected number of point data: %d", npts);
244                    return NULL;
245                }
246                if (isBinary < 0 || fin.eof()) {
247                    // Should know if ASCII or BINARY by now
248                    ERROR("Error in data stream");
249                    return NULL;
250                }
251                fin.getline(line, sizeof(line) - 1);
252                if (fin.fail()) {
253                    ERROR("Error in data stream");
254                    return NULL;
255                }
256                for (start = line; *start == ' ' || *start == '\t'; start++)
257                    ;  // skip leading blanks
258                if (sscanf(start, "SCALARS %s %s", str, type) != 2) {
259                    ERROR("Error in data stream");
260                    return NULL;
261                }
262                if (strncmp(type, "double", 6) == 0) {
263                    ftype = DOUBLE;
264                } else if (strncmp(type, "float", 5) == 0) {
265                    ftype = FLOAT;
266                } else if (strncmp(type, "unsigned_char", 13) == 0) {
267                    ftype = UCHAR;
268                    if (!isBinary) {
269                        ERROR("unsigned char only supported in binary VTK files");
270                        return NULL;
271                    }
272                } else if (strncmp(type, "short", 5) == 0) {
273                    ftype = SHORT;
274                } else if (strncmp(type, "unsigned_short", 14) == 0) {
275                    ftype = USHORT;
276                } else if (strncmp(type, "int", 3) == 0) {
277                    ftype = INT;
278                } else if (strncmp(type, "unsigned_int", 12) == 0) {
279                    ftype = UINT;
280                } else {
281                    ERROR("Unsupported scalar type: '%s'", type);
282                    return NULL;
283                }
284                if (fin.eof()) {
285                    ERROR("Error in data stream");
286                    return NULL;
287                }
288                fin.getline(line, sizeof(line) - 1);
289                if (fin.fail()) {
290                    ERROR("Error in data stream");
291                    return NULL;
292                }
293                for (start = line; *start == ' ' || *start == '\t'; start++)
294                    ;  // skip leading blanks
295                if (sscanf(start, "LOOKUP_TABLE %s", str) == 1) {
296                    // skip lookup table, but don't read ahead
297                    if (fin.eof()) {
298                        ERROR("Error in data stream");
299                        return NULL;
300                    }
301                } else {
302                    // Lookup table line is required
303                    ERROR("Missing LOOKUP_TABLE");
304                    return NULL;
305                }
306                // found start of field data
307                allocFieldData(&fieldData, npts, ftype);
308                if (isBinary) {
309                    fin.read((char *)fieldData, npts * typeSize(ftype));
310                    if (fin.fail()) {
311                        deleteFieldData(fieldData, ftype);
312                        ERROR("Error in data stream");
313                        return NULL;
314                    }
315                    numRead = npts;
316                } else {
317                    numRead = 0;
318                    while (!fin.eof() && numRead < npts) {
319                        readFieldValue(fieldData, numRead++, ftype, fin);
320                        if (fin.fail()) {
321                            deleteFieldData(fieldData, ftype);
322                            ERROR("Error in data stream");
323                            return NULL;
324                        }
325                    }
326                }
327               
328                break;
329            }
330        }
331    }
332
333    if (isRect < 0 || numRead != npts || npts < 1) {
334        deleteFieldData(fieldData, ftype);
335        ERROR("Error in data stream");
336        return NULL;
337    }
338
339    TRACE("found nx=%d ny=%d nz=%d\ndx=%f dy=%f dz=%f\nx0=%f y0=%f z0=%f",
340          nx, ny, nz, dx, dy, dz, x0, y0, z0);
341
342    lx = (nx - 1) * dx;
343    ly = (ny - 1) * dy;
344    lz = (nz - 1) * dz;
345
346    Volume *volume = NULL;
347    double vmin = DBL_MAX;
348    double nzero_min = DBL_MAX;
349    double vmax = -DBL_MAX;
350    float *data = NULL;
351
352    if (isRect) {
353#ifdef DOWNSAMPLE_DATA
354        Rappture::Mesh1D xgrid(x0, x0 + lx, nx);
355        Rappture::Mesh1D ygrid(y0, y0 + ly, ny);
356        Rappture::Mesh1D zgrid(z0, z0 + lz, nz);
357        Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
358#else // !DOWNSAMPLE_DATA
359        data = new float[npts * 4];
360        memset(data, 0, npts * 4);
361#endif // DOWNSAMPLE_DATA
362        for (int p = 0; p < npts; p++) {
363#ifdef DOWNSAMPLE_DATA
364            field.define(p, getFieldData(fieldData, p, ftype));
365#else // !DOWNSAMPLE_DATA
366            int nindex = p * 4;
367            data[nindex] = getFieldData(fieldData, p, ftype);
368            if (data[nindex] < vmin) {
369                vmin = data[nindex];
370            } else if (data[nindex] > vmax) {
371                vmax = data[nindex];
372            }
373            if (data[nindex] != 0.0 && data[nindex] < nzero_min) {
374                nzero_min = data[nindex];
375            }
376#endif // DOWNSAMPLE_DATA
377        }
378
379        deleteFieldData(fieldData, ftype);
380
381#ifndef DOWNSAMPLE_DATA
382
383        // scale all values [0-1], -1 => out of bounds
384        normalizeScalar(data, npts, 4, vmin, vmax);
385        computeSimpleGradient(data, nx, ny, nz,
386                              dx, dy, dz);
387#else // DOWNSAMPLE_DATA
388        // figure out a good mesh spacing
389        int nsample = 64;
390        double dmin = pow((lx*ly*lz)/((nsample-1)*(nsample-1)*(nsample-1)), 0.333);
391
392        nx = (int)ceil(lx/dmin) + 1;
393        ny = (int)ceil(ly/dmin) + 1;
394        nz = (int)ceil(lz/dmin) + 1;
395
396#ifndef HAVE_NPOT_TEXTURES
397        // must be an even power of 2 for older cards
398        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
399        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
400        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
401#endif
402
403        dx = lx /(double)(nx - 1);
404        dy = ly /(double)(ny - 1);
405        dz = lz /(double)(nz - 1);
406
407        vmin = field.valueMin();
408        vmax = field.valueMax();
409        nzero_min = DBL_MAX;
410        double dv = vmax - vmin;
411        if (dv == 0.0) {
412            dv = 1.0;
413        }
414
415        int ngen = 0;
416#ifdef FILTER_GRADIENTS
417        // Sobel filter expects a single float at
418        // each node
419        const int step = 1;
420        float *cdata = new float[npts * step];
421#else // !FILTER_GRADIENTS
422        // Leave 3 floats of space for gradients
423        // to be filled in by computeSimpleGradient()
424        const int step = 4;
425        data = new float[npts * step];
426#endif // FILTER_GRADIENTS
427
428        for (int iz = 0; iz < nz; iz++) {
429            double zval = z0 + iz*dz;
430            for (int iy = 0; iy < ny; iy++) {
431                double yval = y0 + iy*dy;
432                for (int ix = 0; ix < nx; ix++) {
433                    double xval = x0 + ix*dx;
434                    double v = field.value(xval, yval, zval);
435#ifdef notdef
436                    if (isnan(v)) {
437                        TRACE("NAN at %d,%d,%d: (%g, %g, %g)", ix, iy, iz, xval, yval, zval);
438                    }
439#endif
440                    if (v != 0.0 && v < nzero_min) {
441                        nzero_min = v;
442                    }
443#ifdef FILTER_GRADIENTS
444                    // NOT normalized, may have NaNs (FIXME: how does gradient filter handle NaN?)
445                    cdata[ngen] = v;
446#else // !FILTER_GRADIENTS
447                    // scale all values [0-1], -1 => out of bounds
448                    v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
449                    data[ngen] = v;
450#endif // FILTER_GRADIENTS
451                    ngen += step;
452                }
453            }
454        }
455#ifdef FILTER_GRADIENTS
456        // computeGradient returns a new array with gradients
457        // filled in, so data is now 4 floats per node
458        data = computeGradient(cdata, nx, ny, nz,
459                               dx, dy, dz,
460                               vmin, vmax);
461        delete [] cdata;
462#else // !FILTER_GRADIENTS
463        computeSimpleGradient(data, nx, ny, nz,
464                              dx, dy, dz);
465#endif // FILTER_GRADIENTS
466
467#endif // DOWNSAMPLE_DATA
468    } else {
469        deleteFieldData(fieldData, ftype);
470        ERROR("Unsupported DataSet");
471        return NULL;
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.