source: trunk/packages/vizservers/nanovis/VtkReader.cpp @ 3628

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

Use nv namespace for classes in nanovis rather than prefixing class names with
Nv (still need to convert shader classes).
  • Property svn:eol-style set to native
File size: 15.0 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 <RpOutcome.h>
11#include <RpField1D.h>
12#include <RpFieldRect3D.h>
13
14#include <vrmath/Vector3f.h>
15
16#include "ReaderCommon.h"
17
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(Rappture::Outcome& result, 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            result.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                    result.addError("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                    result.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                    result.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                    result.addError("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                    result.error("Error in data stream");
249                    return NULL;
250                }
251                fin.getline(line, sizeof(line) - 1);
252                if (fin.fail()) {
253                    result.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                    result.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                        result.addError("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                    result.addError("Unsupported scalar type: '%s'", type);
282                    return NULL;
283                }
284                if (fin.eof()) {
285                    result.error("Error in data stream");
286                    return NULL;
287                }
288                fin.getline(line, sizeof(line) - 1);
289                if (fin.fail()) {
290                    result.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                        result.error("Error in data stream");
299                        return NULL;
300                    }
301                } else {
302                    // Lookup table line is required
303                    result.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                        result.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                            result.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        result.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        result.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    //
491    // Center this new volume on the origin.
492    //
493    float dx0 = -0.5;
494    float dy0 = -0.5*ly/lx;
495    float dz0 = -0.5*lz/lx;
496    if (volume) {
497        volume->location(vrmath::Vector3f(dx0, dy0, dz0));
498        TRACE("Set volume location to %g %g %g", dx0, dy0, dz0);
499    }
500    return volume;
501}
Note: See TracBrowser for help on using the repository browser.