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

Last change on this file since 3596 was 3576, checked in by ldelgass, 12 years ago

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