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

Last change on this file since 3568 was 3558, checked in by gah, 11 years ago

small fix for double type in vtk structure points reader.
  • Property svn:eol-style set to native
File size: 14.9 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                // found dimensions
227            } else if (sscanf(start, "ORIGIN %lg %lg %lg", &x0, &y0, &z0) == 3) {
228                // found origin
229            } else if (sscanf(start, "SPACING %lg %lg %lg", &dx, &dy, &dz) == 3) {
230                // found cell spacing
231            } else if (sscanf(start, "POINT_DATA %d", &npts) == 1) {
232                if (npts != nx * ny * nz) {
233                    result.addError("Error in data stream: unexpected number of point data: %d", npts);
234                    return NULL;
235                }
236                if (isBinary < 0 || fin.eof()) {
237                    // Should know if ASCII or BINARY by now
238                    result.error("Error in data stream");
239                    return NULL;
240                }
241                fin.getline(line, sizeof(line) - 1);
242                if (fin.fail()) {
243                    result.error("Error in data stream");
244                    return NULL;
245                }
246                for (start = line; *start == ' ' || *start == '\t'; start++)
247                    ;  // skip leading blanks
248                if (sscanf(start, "SCALARS %s %s", str, type) != 2) {
249                    result.error("Error in data stream");
250                    return NULL;
251                }
252                if (strncmp(type, "double", 6) == 0) {
253                    ftype = DOUBLE;
254                } else if (strncmp(type, "float", 5) == 0) {
255                    ftype = FLOAT;
256                } else if (strncmp(type, "unsigned_char", 13) == 0) {
257                    ftype = UCHAR;
258                    if (!isBinary) {
259                        result.addError("unsigned char only supported in binary VTK files");
260                        return NULL;
261                    }
262                } else if (strncmp(type, "short", 5) == 0) {
263                    ftype = SHORT;
264                } else if (strncmp(type, "unsigned_short", 14) == 0) {
265                    ftype = USHORT;
266                } else if (strncmp(type, "int", 3) == 0) {
267                    ftype = INT;
268                } else if (strncmp(type, "unsigned_int", 12) == 0) {
269                    ftype = UINT;
270                } else {
271                    result.addError("Unsupported scalar type: '%s'", type);
272                    return NULL;
273                }
274                if (fin.eof()) {
275                    result.error("Error in data stream");
276                    return NULL;
277                }
278                fin.getline(line, sizeof(line) - 1);
279                if (fin.fail()) {
280                    result.error("Error in data stream");
281                    return NULL;
282                }
283                for (start = line; *start == ' ' || *start == '\t'; start++)
284                    ;  // skip leading blanks
285                if (sscanf(start, "LOOKUP_TABLE %s", str) == 1) {
286                    // skip lookup table, but don't read ahead
287                    if (fin.eof()) {
288                        result.error("Error in data stream");
289                        return NULL;
290                    }
291                } else {
292                    // Lookup table line is required
293                    result.error("Missing LOOKUP_TABLE");
294                    return NULL;
295                }
296                // found start of field data
297                allocFieldData(&fieldData, npts, ftype);
298                if (isBinary) {
299                    fin.read((char *)fieldData, npts * typeSize(ftype));
300                    if (fin.fail()) {
301                        deleteFieldData(fieldData, ftype);
302                        result.error("Error in data stream");
303                        return NULL;
304                    }
305                    numRead = npts;
306                } else {
307                    numRead = 0;
308                    while (!fin.eof() && numRead < npts) {
309                        readFieldValue(fieldData, numRead++, ftype, fin);
310                        if (fin.fail()) {
311                            deleteFieldData(fieldData, ftype);
312                            result.error("Error in data stream");
313                            return NULL;
314                        }
315                    }
316                }
317               
318                break;
319            }
320        }
321    }
322
323    if (isRect < 0 || numRead != npts) {
324        deleteFieldData(fieldData, ftype);
325        result.error("Error in data stream");
326        return NULL;
327    }
328
329    TRACE("found nx=%d ny=%d nz=%d\ndx=%f dy=%f dz=%f\nx0=%f y0=%f z0=%f",
330          nx, ny, nz, dx, dy, dz, x0, y0, z0);
331
332    lx = (nx - 1) * dx;
333    ly = (ny - 1) * dy;
334    lz = (nz - 1) * dz;
335
336    Volume *volPtr = NULL;
337    double vmin = DBL_MAX;
338    double nzero_min = DBL_MAX;
339    double vmax = -DBL_MAX;
340    float *data = NULL;
341
342    if (isRect) {
343#ifdef DOWNSAMPLE_DATA
344        Rappture::Mesh1D xgrid(x0, x0 + lx, nx);
345        Rappture::Mesh1D ygrid(y0, y0 + ly, ny);
346        Rappture::Mesh1D zgrid(z0, z0 + lz, nz);
347        Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
348#else // !DOWNSAMPLE_DATA
349        data = new float[npts * 4];
350        memset(data, 0, npts * 4);
351#endif // DOWNSAMPLE_DATA
352        for (int p = 0; p < npts; p++) {
353#ifdef DOWNSAMPLE_DATA
354            field.define(p, getFieldData(fieldData, p, ftype));
355#else // !DOWNSAMPLE_DATA
356            int nindex = p * 4;
357            data[nindex] = getFieldData(fieldData, p, ftype);
358            if (data[nindex] < vmin) {
359                vmin = data[nindex];
360            } else if (data[nindex] > vmax) {
361                vmax = data[nindex];
362            }
363            if (data[nindex] != 0.0 && data[nindex] < nzero_min) {
364                nzero_min = data[nindex];
365            }
366#endif // DOWNSAMPLE_DATA
367        }
368
369        deleteFieldData(fieldData, ftype);
370
371#ifndef DOWNSAMPLE_DATA
372        double dv = vmax - vmin;
373        if (dv == 0.0) {
374            dv = 1.0;
375        }
376
377        int ngen = 0;
378        const int step = 4;
379        for (int i = 0; i < npts; ++i) {
380            double v = data[ngen];
381            // scale all values [0-1], -1 => out of bounds
382            v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
383
384            data[ngen] = v;
385            ngen += step;
386        }
387
388        computeSimpleGradient(data, nx, ny, nz,
389                              dx, dy, dz);
390#else // DOWNSAMPLE_DATA
391        // figure out a good mesh spacing
392        int nsample = 30;
393        double dmin = pow((lx*ly*lz)/((nsample-1)*(nsample-1)*(nsample-1)), 0.333);
394
395        nx = (int)ceil(lx/dmin) + 1;
396        ny = (int)ceil(ly/dmin) + 1;
397        nz = (int)ceil(lz/dmin) + 1;
398
399#ifndef HAVE_NPOT_TEXTURES
400        // must be an even power of 2 for older cards
401        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
402        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
403        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
404#endif
405
406        dx = lx /(double)(nx - 1);
407        dy = ly /(double)(ny - 1);
408        dz = lz /(double)(nz - 1);
409
410        vmin = field.valueMin();
411        vmax = field.valueMax();
412        nzero_min = DBL_MAX;
413        double dv = vmax - vmin;
414        if (dv == 0.0) {
415            dv = 1.0;
416        }
417
418        int ngen = 0;
419#ifdef FILTER_GRADIENTS
420        // Sobel filter expects a single float at
421        // each node
422        const int step = 1;
423        float *cdata = new float[npts * step];
424#else // !FILTER_GRADIENTS
425        // Leave 3 floats of space for gradients
426        // to be filled in by computeSimpleGradient()
427        const int step = 4;
428        data = new float[npts * step];
429#endif // FILTER_GRADIENTS
430
431        for (int iz = 0; iz < nz; iz++) {
432            double zval = z0 + iz*dz;
433            for (int iy = 0; iy < ny; iy++) {
434                double yval = y0 + iy*dy;
435                for (int ix = 0; ix < nx; ix++) {
436                    double xval = x0 + ix*dx;
437                    double v = field.value(xval, yval, zval);
438#ifdef notdef
439                    if (isnan(v)) {
440                        TRACE("NAN at %d,%d,%d: (%g, %g, %g)", ix, iy, iz, xval, yval, zval);
441                    }
442#endif
443                    if (v != 0.0 && v < nzero_min) {
444                        nzero_min = v;
445                    }
446#ifdef FILTER_GRADIENTS
447                    // NOT normalized, -1 => out of bounds
448                    v = (isnan(v)) ? -1.0 : v;
449                    cdata[ngen] = v;
450#else // !FILTER_GRADIENTS
451                    // scale all values [0-1], -1 => out of bounds
452                    v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
453                    data[ngen] = v;
454#endif // FILTER_GRADIENTS
455                    ngen += step;
456                }
457            }
458        }
459#ifdef FILTER_GRADIENTS
460        // computeGradient returns a new array with gradients
461        // filled in, so data is now 4 floats per node
462        data = computeGradient(cdata, nx, ny, nz,
463                               dx, dy, dz,
464                               vmin, vmax);
465        delete [] cdata;
466#else // !FILTER_GRADIENTS
467        computeSimpleGradient(data, nx, ny, nz,
468                              dx, dy, dz);
469#endif // FILTER_GRADIENTS
470
471#endif // DOWNSAMPLE_DATA
472    } else {
473        deleteFieldData(fieldData, ftype);
474        result.error("Unsupported DataSet");
475        return NULL;
476    }
477
478    TRACE("nx = %i ny = %i nz = %i", nx, ny, nz);
479    TRACE("x0 = %lg y0 = %lg z0 = %lg", x0, y0, z0);
480    TRACE("lx = %lg ly = %lg lz = %lg", lx, ly, lz);
481    TRACE("dx = %lg dy = %lg dz = %lg", dx, dy, dz);
482    TRACE("dataMin = %lg dataMax = %lg nzero_min = %lg",
483          vmin, vmax, nzero_min);
484
485    volPtr = NanoVis::loadVolume(tag, nx, ny, nz, 4, data,
486                                 vmin, vmax, nzero_min);
487    volPtr->xAxis.setRange(x0, x0 + lx);
488    volPtr->yAxis.setRange(y0, y0 + ly);
489    volPtr->zAxis.setRange(z0, z0 + lz);
490    volPtr->updatePending = true;
491
492    delete [] data;
493
494    //
495    // Center this new volume on the origin.
496    //
497    float dx0 = -0.5;
498    float dy0 = -0.5*ly/lx;
499    float dz0 = -0.5*lz/lx;
500    if (volPtr) {
501        volPtr->location(vrmath::Vector3f(dx0, dy0, dz0));
502        TRACE("Set volume location to %g %g %g", dx0, dy0, dz0);
503    }
504    return volPtr;
505}
Note: See TracBrowser for help on using the repository browser.