source: branches/1.3/packages/vizservers/nanovis/VtkReader.cpp @ 3877

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

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