source: trunk/packages/vizservers/nanovis/Unirect.cpp @ 3362

Last change on this file since 3362 was 2953, checked in by ldelgass, 12 years ago

Remove unused global origin, make default transfer function a bit more
sensible (used to have full opacity at 0). Fix HeightMap? dtor to use delete[]
instead of free() on array allocated with new[]. Document data response in
protocol.
  • Property svn:eol-style set to native
File size: 25.5 KB
Line 
1/* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */
2#include <float.h>
3#include <tcl.h>
4
5#include <RpField1D.h>
6#include <RpFieldRect3D.h>
7
8#include "Unirect.h"
9#include "Trace.h"
10
11extern int GetFloatFromObj(Tcl_Interp *interp, Tcl_Obj *objPtr,
12                           float *valuePtr);
13
14extern int GetAxisFromObj(Tcl_Interp *interp, Tcl_Obj *objPtr, int *indexPtr);
15
16static inline char *
17skipspaces(char *string)
18{
19    while (isspace(*string)) {
20        string++;
21    }
22    return string;
23}
24
25static inline char *
26getline(char **stringPtr, char *endPtr)
27{
28    char *line, *p;
29
30    line = skipspaces(*stringPtr);
31    for (p = line; p < endPtr; p++) {
32        if (*p == '\n') {
33            *p++ = '\0';
34            *stringPtr = p;
35            return line;
36        }
37    }
38    *stringPtr = p;
39    return line;
40}
41
42int
43Rappture::Unirect2d::parseBuffer(Tcl_Interp *interp, Rappture::Buffer &buf)
44{
45    Tcl_Obj *objPtr;
46    objPtr = Tcl_NewStringObj(buf.bytes(), buf.size());
47    Tcl_Obj **objv;
48    int objc;
49    if (Tcl_ListObjGetElements(interp, objPtr, &objc, &objv) != TCL_OK) {
50        return TCL_ERROR;
51    }
52    int result;
53    result = loadData(interp, objc, objv);
54    Tcl_DecrRefCount(objPtr);
55    if ((result != TCL_OK) || (!isInitialized())) {
56        return TCL_ERROR;
57    }
58    return TCL_OK;
59}
60
61int
62Rappture::Unirect3d::parseBuffer(Tcl_Interp *interp, Rappture::Buffer &buf)
63{
64    Tcl_Obj *objPtr;
65    objPtr = Tcl_NewStringObj(buf.bytes(), buf.size());
66    Tcl_Obj **objv;
67    int objc;
68    if (Tcl_ListObjGetElements(interp, objPtr, &objc, &objv) != TCL_OK) {
69        return TCL_ERROR;
70    }
71    int result;
72    result = loadData(interp, objc, objv);
73    Tcl_DecrRefCount(objPtr);
74    if ((result != TCL_OK) || (!isInitialized())) {
75        return TCL_ERROR;
76    }
77    return TCL_OK;
78}
79
80int
81Rappture::Unirect3d::loadData(Tcl_Interp *interp, int objc,
82                              Tcl_Obj *const *objv)
83{
84    int num[3], nValues;
85    float min[3], max[3];
86    float *values;
87    const char *units[4], *order;
88
89
90    if ((objc-1) & 0x01) {
91        Tcl_AppendResult(interp, Tcl_GetString(objv[0]), ": ",
92                "wrong # of arguments: should be key-value pairs",
93                (char *)NULL);
94        return TCL_ERROR;
95    }
96
97    /* Default order is  z, y, x. */
98    int axis1, axis2, axis3;
99    axis1 = 0;                  /* X-axis */
100    axis2 = 1;                  /* Y-axis */
101    axis3 = 2;                  /* Z-axis */
102
103    values = NULL;
104    num[0] = num[1] = num[2] = nValues = 0;
105    min[0] = min[1] = min[2] = max[0] = max[1] = max[2] = 0.0f;
106    order = units[0] = units[1] = units[2] = units[3] = NULL;
107
108    int i;
109    for (i = 1; i < objc; i += 2) {
110        const char *string;
111        char c;
112
113        string = Tcl_GetString(objv[i]);
114        c = string[0];
115        if ((c == 'x') && (strcmp(string, "xmin") == 0)) {
116            if (GetFloatFromObj(interp, objv[i+1], min+2) != TCL_OK) {
117                return TCL_ERROR;
118            }
119        } else if ((c == 'x') && (strcmp(string, "xmax") == 0)) {
120            if (GetFloatFromObj(interp, objv[i+1], max+2) != TCL_OK) {
121                return TCL_ERROR;
122            }
123        } else if ((c == 'x') && (strcmp(string, "xnum") == 0)) {
124            if (Tcl_GetIntFromObj(interp, objv[i+1], num+2) != TCL_OK) {
125                return TCL_ERROR;
126            }
127            if (num[2] <= 0) {
128                Tcl_AppendResult(interp, "bad xnum value: must be > 0",
129                     (char *)NULL);
130                return TCL_ERROR;
131            }
132        } else if ((c == 'x') && (strcmp(string, "xunits") == 0)) {
133            units[2] = Tcl_GetString(objv[i+1]);
134        } else if ((c == 'y') && (strcmp(string, "ymin") == 0)) {
135            if (GetFloatFromObj(interp, objv[i+1], min+1) != TCL_OK) {
136                return TCL_ERROR;
137            }
138        } else if ((c == 'y') && (strcmp(string, "ymax") == 0)) {
139            if (GetFloatFromObj(interp, objv[i+1], max+1) != TCL_OK) {
140                return TCL_ERROR;
141            }
142        } else if ((c == 'y') && (strcmp(string, "ynum") == 0)) {
143            if (Tcl_GetIntFromObj(interp, objv[i+1], num+1) != TCL_OK) {
144                return TCL_ERROR;
145            }
146            if (num[1] <= 0) {
147                Tcl_AppendResult(interp, "bad ynum value: must be > 0",
148                                 (char *)NULL);
149                return TCL_ERROR;
150            }
151        } else if ((c == 'y') && (strcmp(string, "yunits") == 0)) {
152            units[1] = Tcl_GetString(objv[i+1]);
153        } else if ((c == 'z') && (strcmp(string, "zmin") == 0)) {
154            if (GetFloatFromObj(interp, objv[i+1], min) != TCL_OK) {
155                return TCL_ERROR;
156            }
157        } else if ((c == 'z') && (strcmp(string, "zmax") == 0)) {
158            if (GetFloatFromObj(interp, objv[i+1], max) != TCL_OK) {
159                return TCL_ERROR;
160            }
161        } else if ((c == 'z') && (strcmp(string, "znum") == 0)) {
162            if (Tcl_GetIntFromObj(interp, objv[i+1], num) != TCL_OK) {
163                return TCL_ERROR;
164            }
165            if (num[0] <= 0) {
166                Tcl_AppendResult(interp, "bad znum value: must be > 0",
167                                 (char *)NULL);
168                return TCL_ERROR;
169            }
170        } else if ((c == 'z') && (strcmp(string, "zunits") == 0)) {
171            units[0] = Tcl_GetString(objv[i+1]);
172        } else if ((c == 'v') && (strcmp(string, "values") == 0)) {
173            Tcl_Obj **vobj;
174
175            if (Tcl_ListObjGetElements(interp, objv[i+1], &nValues, &vobj)
176                != TCL_OK) {
177                return TCL_ERROR;
178            }
179            _values = (float *)realloc(_values, sizeof(float) * nValues);
180            int j;
181            for (j = 0; j < nValues; j++) {
182                if (GetFloatFromObj(interp, vobj[j], _values + j)!=TCL_OK) {
183                    return TCL_ERROR;
184                }
185            }
186        } else if ((c == 'u') && (strcmp(string, "units") == 0)) {
187            _vUnits = strdup(Tcl_GetString(objv[i+1]));
188        } else if ((c == 'c') && (strcmp(string, "components") == 0)) {
189            int n;
190
191            if (Tcl_GetIntFromObj(interp, objv[i+1], &n) != TCL_OK) {
192                return TCL_ERROR;
193            }
194            if (n <= 0) {
195                Tcl_AppendResult(interp, "bad extents value: must be > 0",
196                                 (char *)NULL);
197                return TCL_ERROR;
198            }
199            _nComponents = n;
200        } else if ((c == 'a') && (strcmp(string, "axisorder") == 0)) {
201            Tcl_Obj **axes;
202            int n;
203
204            if (Tcl_ListObjGetElements(interp, objv[i+1], &n, &axes)
205                != TCL_OK) {
206                return TCL_ERROR;
207            }
208            if (n != 3) {
209                return TCL_ERROR;
210            }
211            if ((GetAxisFromObj(interp, axes[0], &axis1) != TCL_OK) ||
212                (GetAxisFromObj(interp, axes[1], &axis2) != TCL_OK) ||
213                (GetAxisFromObj(interp, axes[2], &axis3) != TCL_OK)) {
214                return TCL_ERROR;
215            }
216        } else {
217            Tcl_AppendResult(interp, "unknown key \"", string,
218                (char *)NULL);
219            return TCL_ERROR;
220        }
221    }
222    if (_values == NULL) {
223        Tcl_AppendResult(interp, "missing \"values\" key", (char *)NULL);
224        return TCL_ERROR;
225    }
226    if ((size_t)nValues != (num[0] * num[1] * num[2] * _nComponents)) {
227        TRACE("num[0]=%d num[1]=%d num[2]=%d ncomponents=%d nValues=%d\n",
228               num[0], num[1], num[2], _nComponents, nValues);
229        Tcl_AppendResult(interp,
230                "wrong # of values: must be xnum*ynum*znum*extents",
231                         (char *)NULL);
232       return TCL_ERROR;
233    }
234   
235#ifdef notdef
236    if ((axis1 != 0) || (axis2 != 1) || (axis3 != 2)) {
237        // Reorder the data into x, y, z where x varies fastest and so on.
238        int z;
239        float *data, *dp;
240
241        dp = data = new float[nValues];
242        for (z = 0; z < num[0]; z++) {
243            int y;
244
245            for (y = 0; y < num[1]; y++) {
246                int x;
247
248                for (x = 0; x < num[2]; x++) {
249                    int i;
250                   
251                    /* Compute the index from the data's described ordering. */
252                    i = ((z*num[axis2]*num[axis3]) + (y*num[axis3]) + x) * 3;
253                    for(size_t v = 0; v < _nComponents; v++) {
254                        dp[v] = values[i+v];
255                    }
256                    dp += _nComponents;
257                }
258            }
259        }
260        delete [] values;
261        values = data;
262    }
263#endif
264    _nValues = nValues;
265    if (units[3] != NULL) {
266        if (_vUnits != NULL) {
267            free(_vUnits);
268        }
269        _vUnits = strdup(units[3]);
270    }
271    _xMin = min[axis3];
272    _xMax = max[axis3];
273    _xNum = num[axis3];
274    if (units[axis3] != NULL) {
275        if (_xUnits != NULL) {
276            free(_xUnits);
277        }
278        _xUnits = strdup(units[axis3]);
279    }
280    _yMin = min[axis2];
281    _yMax = max[axis2];
282    _yNum = num[axis2];
283    if (units[axis2] != NULL) {
284        if (_yUnits != NULL) {
285            free(_yUnits);
286        }
287        _yUnits = strdup(units[axis2]);
288    }
289    _zMin = min[axis1];
290    _zMax = max[axis1];
291    _zNum = num[axis1];
292    if (units[axis1] != NULL) {
293        if (_zUnits != NULL) {
294            free(_zUnits);
295        }
296        _zUnits = strdup(units[axis1]);
297    }
298    _initialized = true;
299#ifdef notdef
300    {
301        FILE *f;
302        f = fopen("/tmp/unirect3d.txt", "w");
303        fprintf(f, "unirect3d xmin %g xmax %g xnum %d ", _xMin, _xMax, _xNum);
304        fprintf(f, "ymin %g ymax %g ynum %d ", _yMin, _yMax, _yNum);
305        fprintf(f, "zmin %g zmax %g znum %d ", _zMin, _zMax, _zNum);
306        fprintf(f, "components %d values {\n",  _nComponents);
307        for (size_t i = 0; i < _nValues; i+= 3) {
308            fprintf(f, "%g %g %g\n", _values[i], _values[i+1], _values[i+2]);
309        }
310        fprintf(f, "}\n");
311        fclose(f);
312    }
313#endif
314    return TCL_OK;
315}
316
317int
318Rappture::Unirect2d::loadData(Tcl_Interp *interp, int objc,
319                              Tcl_Obj *const *objv)
320{
321    if ((objc & 0x01) == 0) {
322        Tcl_AppendResult(interp, Tcl_GetString(objv[0]), ": ",
323                "wrong number of arguments: should be key-value pairs",
324                (char *)NULL);
325        return TCL_ERROR;
326    }
327
328    int axis[2];
329    axis[0] = 1;                         /* X-axis */
330    axis[1] = 0;                        /* Y-axis */
331
332    _xNum = _yNum = _nValues = 0;
333    _xMin = _yMin = _xMax = _yMax = 0.0f;
334    if (_xUnits != NULL) {
335        free(_xUnits);
336    }
337    if (_yUnits != NULL) {
338        free(_yUnits);
339    }
340    if (_vUnits != NULL) {
341        free(_vUnits);
342    }
343    _xUnits = _yUnits = _vUnits = NULL;
344    _nValues = 0;
345
346    int i;
347    for (i = 1; i < objc; i += 2) {
348        const char *string;
349        char c;
350
351        string = Tcl_GetString(objv[i]);
352        c = string[0];
353        if ((c == 'x') && (strcmp(string, "xmin") == 0)) {
354            if (GetFloatFromObj(interp, objv[i+1], &_xMin) != TCL_OK) {
355                return TCL_ERROR;
356            }
357        } else if ((c == 'x') && (strcmp(string, "xmax") == 0)) {
358            if (GetFloatFromObj(interp, objv[i+1], &_xMax) != TCL_OK) {
359                return TCL_ERROR;
360            }
361        } else if ((c == 'x') && (strcmp(string, "xnum") == 0)) {
362            int n;
363            if (Tcl_GetIntFromObj(interp, objv[i+1], &n) != TCL_OK) {
364                return TCL_ERROR;
365            }
366            if (n <= 0) {
367                Tcl_AppendResult(interp, "bad xnum value: must be > 0",
368                     (char *)NULL);
369                return TCL_ERROR;
370            }
371            _xNum = n;
372        } else if ((c == 'x') && (strcmp(string, "xunits") == 0)) {
373            _xUnits = strdup(Tcl_GetString(objv[i+1]));
374        } else if ((c == 'y') && (strcmp(string, "ymin") == 0)) {
375            if (GetFloatFromObj(interp, objv[i+1], &_yMin) != TCL_OK) {
376                return TCL_ERROR;
377            }
378        } else if ((c == 'y') && (strcmp(string, "ymax") == 0)) {
379            if (GetFloatFromObj(interp, objv[i+1], &_yMax) != TCL_OK) {
380                return TCL_ERROR;
381            }
382        } else if ((c == 'y') && (strcmp(string, "ynum") == 0)) {
383            int n;
384            if (Tcl_GetIntFromObj(interp, objv[i+1], &n) != TCL_OK) {
385                return TCL_ERROR;
386            }
387            if (n <= 0) {
388                Tcl_AppendResult(interp, "bad ynum value: must be > 0",
389                                 (char *)NULL);
390                return TCL_ERROR;
391            }
392            _yNum = n;
393        } else if ((c == 'y') && (strcmp(string, "yunits") == 0)) {
394            _yUnits = strdup(Tcl_GetString(objv[i+1]));
395        } else if ((c == 'v') && (strcmp(string, "values") == 0)) {
396            Tcl_Obj **vobj;
397            int n;
398
399            Tcl_IncrRefCount(objv[i+1]);
400            if (Tcl_ListObjGetElements(interp, objv[i+1], &n, &vobj) != TCL_OK){
401                return TCL_ERROR;
402            }
403            if (n <= 0) {
404                Tcl_AppendResult(interp, "empty values list : must be > 0",
405                                 (char *)NULL);
406                return TCL_ERROR;
407            }
408            _nValues = n;
409            _values = (float *)realloc(_values, sizeof(float) * _nValues);
410            size_t j;
411            for (j = 0; j < _nValues; j++) {
412                if (GetFloatFromObj(interp, vobj[j], _values + j)!=TCL_OK) {
413                    return TCL_ERROR;
414                }
415            }
416            Tcl_DecrRefCount(objv[i+1]);
417        } else if ((c == 'u') && (strcmp(string, "units") == 0)) {
418            _vUnits = strdup(Tcl_GetString(objv[i+1]));
419        } else if ((c == 'c') && (strcmp(string, "components") == 0)) {
420            int n;
421
422            if (Tcl_GetIntFromObj(interp, objv[i+1], &n) != TCL_OK) {
423                return TCL_ERROR;
424            }
425            if (n <= 0) {
426                Tcl_AppendResult(interp, "bad extents value: must be > 0",
427                                 (char *)NULL);
428                return TCL_ERROR;
429            }
430            _nComponents = n;
431        } else if ((c == 'a') && (strcmp(string, "axisorder") == 0)) {
432            Tcl_Obj **order;
433            int n;
434
435            if (Tcl_ListObjGetElements(interp, objv[i+1], &n, &order)
436                != TCL_OK) {
437                return TCL_ERROR;
438            }
439            if (n != 2) {
440                Tcl_AppendResult(interp,
441                        "wrong # of axes defined for ordering the data",
442                        (char *)NULL);
443                return TCL_ERROR;
444            }
445            if ((GetAxisFromObj(interp, order[0], axis) != TCL_OK) ||
446                (GetAxisFromObj(interp, order[1], axis+1) != TCL_OK)) {
447                return TCL_ERROR;
448            }
449        } else {
450            Tcl_AppendResult(interp, "unknown key \"", string,
451                (char *)NULL);
452            return TCL_ERROR;
453        }
454    }
455    if (_nValues == 0) {
456        Tcl_AppendResult(interp, "missing \"values\" key", (char *)NULL);
457        return TCL_ERROR;
458    }
459    if (_nValues != (_xNum * _yNum * _nComponents)) {
460        Tcl_AppendResult(interp,
461                "wrong number of values: must be xnum*ynum*components",
462                         (char *)NULL);
463        return TCL_ERROR;
464    }
465   
466    if ((axis[0] != 1) || (axis[1] != 0)) {
467        TRACE("reordering data\n");
468        /* Reorder the data into x, y where x varies fastest and so on. */
469        size_t y;
470        float *dp;
471
472        dp = _values = (float *)realloc(_values, sizeof(float) * _nValues);
473        for (y = 0; y < _yNum; y++) {
474            size_t x;
475
476            for (x = 0; x < _xNum; x++) {
477                size_t i, v;
478                   
479                /* Compute the index from the data's described ordering. */
480                i = (y + (_yNum * x)) * _nComponents;
481                for(v = 0; v < _nComponents; v++) {
482                    dp[v] = _values[i+v];
483                }
484                dp += _nComponents;
485            }
486        }
487    }
488    _initialized = true;
489    return TCL_OK;
490}
491
492bool
493Rappture::Unirect3d::importDx(Rappture::Outcome &result, size_t nComponents,
494                              size_t length, char *string)
495{
496    int nx, ny, nz, npts;
497    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
498    char *p, *pend;
499
500    dx = dy = dz = 0.0;                 /* Suppress compiler warning. */
501    x0 = y0 = z0 = 0.0;                 /* May not have an origin line. */
502    nx = ny = nz = npts = 0;            /* Suppress compiler warning. */
503    for (p = string, pend = p + length; p < pend; /*empty*/) {
504        char *line;
505
506        line = getline(&p, pend);
507        if (line == pend) {
508            break;                        /* EOF */
509        }
510        if ((line[0] == '#') || (line == '\0')) {
511            continue;                        /* Skip blank or comment lines. */
512        }
513        if (sscanf(line, "object %*d class gridpositions counts %d %d %d",
514                   &nx, &ny, &nz) == 3) {
515            if ((nx < 0) || (ny < 0) || (nz < 0)) {
516                result.addError("invalid grid size: x=%d, y=%d, z=%d",
517                        nx, ny, nz);
518                return false;
519            }
520        } else if (sscanf(line, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
521            /* Found origin. */
522        } else if (sscanf(line, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
523            /* Found one of the delta lines. */
524            if (ddx != 0.0) {
525                dx = ddx;
526            } else if (ddy != 0.0) {
527                dy = ddy;
528            } else if (ddz != 0.0) {
529                dz = ddz;
530            }
531        } else if (sscanf(line, "object %*d class array type %*s shape 3"
532                " rank 1 items %d data follows", &npts) == 1) {
533            if (npts < 0) {
534                result.addError("bad # points %d", npts);
535                return false;
536            }
537            TRACE("#points=%d\n", npts);
538            if (npts != nx*ny*nz) {
539                result.addError("inconsistent data: expected %d points"
540                                " but found %d points", nx*ny*nz, npts);
541                return false;
542            }
543            break;
544        } else if (sscanf(line, "object %*d class array type %*s rank 0"
545                " times %d data follows", &npts) == 1) {
546            if (npts != nx*ny*nz) {
547                result.addError("inconsistent data: expected %d points"
548                                " but found %d points", nx*ny*nz, npts);
549                return false;
550            }
551            break;
552        }
553    }
554    if (npts != nx*ny*nz) {
555        result.addError("inconsistent data: expected %d points"
556                        " but found %d points", nx*ny*nz, npts);
557        return false;
558    }
559
560    _initialized = false;
561    _xValueMin = _yValueMin = _zValueMin = FLT_MAX;
562    _xValueMax = _yValueMax = _zValueMax = -FLT_MAX;
563    _xMin = x0, _yMin = y0, _zMin = z0;
564    _xNum = nx, _yNum = ny, _zNum = nz;
565    _xMax = _xMin + dx * _xNum;
566    _yMax = _yMin + dy * _yNum;
567    _zMax = _zMin + dz * _zNum;
568    _nComponents = nComponents;
569
570    _values = (float *)realloc(_values, sizeof(float) * npts * _nComponents);
571    _nValues = 0;
572    for (size_t ix = 0; ix < _xNum; ix++) {
573        for (size_t iy = 0; iy < _yNum; iy++) {
574            for (size_t iz = 0; iz < _zNum; iz++) {
575                char *line;
576                if ((p == pend) || (_nValues > (size_t)npts)) {
577                    break;
578                }
579                line = getline(&p, pend);
580                if ((line[0] == '#') || (line[0] == '\0')) {
581                    continue;                /* Skip blank or comment lines. */
582                }
583                double vx, vy, vz;
584                if (sscanf(line, "%lg %lg %lg", &vx, &vy, &vz) == 3) {
585                    int nindex = (iz*nx*ny + iy*nx + ix) * 3;
586                    if (vx < _xValueMin) {
587                        _xValueMin = vx;
588                    } else if (vx > _xValueMax) {
589                        _xValueMax = vx;
590                    }
591                    if (vy < _yValueMin) {
592                        _yValueMin = vy;
593                    } else if (vy > _yValueMax) {
594                        _yValueMax = vy;
595                    }
596                    if (vz < _zValueMin) {
597                        _zValueMin = vz;
598                    } else if (vz > _zValueMax) {
599                        _zValueMax = vz;
600                    }
601                    _values[nindex] = vx;
602                    _values[nindex+1] = vy;
603                    _values[nindex+2] = vz;
604                    _nValues++;
605                }
606            }
607        }
608    }
609    /* Make sure that we read all of the expected points. */
610    if (_nValues != (size_t)npts) {
611        result.addError("inconsistent data: expected %d points"
612                        " but found %d points", npts, _nValues);
613        free(_values);
614        _values = NULL;
615        return false;
616    }
617    _nValues *= _nComponents;
618    _initialized = true;
619#ifdef notdef
620    {
621        FILE *f;
622        f = fopen("/tmp/dx.txt", "w");
623        fprintf(f, "unirect3d xmin %g xmax %g xnum %d ", _xMin, _xMax, _xNum);
624        fprintf(f, "ymin %g ymax %g ynum %d ", _yMin, _yMax, _yNum);
625        fprintf(f, "zmin %g zmax %g znum %d ", _zMin, _zMax, _zNum);
626        fprintf(f, "components %d values {\n",  _nComponents);
627        for (size_t i = 0; i < _nValues; i+= 3) {
628            fprintf(f, "%g %g %g\n", _values[i], _values[i+1], _values[i+2]);
629        }
630        fprintf(f, "}\n");
631        fclose(f);
632    }
633#endif
634    return true;
635}
636
637bool
638Rappture::Unirect3d::resample(Rappture::Outcome &result, size_t nSamples)
639{
640    Rappture::Mesh1D xgrid(_xMin, _xMax, _xNum);
641    Rappture::Mesh1D ygrid(_yMin, _yMax, _yNum);
642    Rappture::Mesh1D zgrid(_zMin, _zMax, _zNum);
643    Rappture::FieldRect3D xfield(xgrid, ygrid, zgrid);
644    Rappture::FieldRect3D yfield(xgrid, ygrid, zgrid);
645    Rappture::FieldRect3D zfield(xgrid, ygrid, zgrid);
646
647    size_t i, j;
648    for (i = 0, j = 0; i < _nValues; i += _nComponents, j++) {
649        double vx, vy, vz;
650
651        vx = _values[i];
652        vy = _values[i+1];
653        vz = _values[i+2];
654       
655        xfield.define(j, vx);
656        yfield.define(j, vy);
657        zfield.define(j, vz);
658    }
659    // Figure out a good mesh spacing
660    double dx, dy, dz;
661    double lx, ly, lz;
662    lx = xfield.rangeMax(Rappture::xaxis) - xfield.rangeMin(Rappture::xaxis);
663    ly = xfield.rangeMax(Rappture::yaxis) - xfield.rangeMin(Rappture::yaxis);
664    lz = xfield.rangeMax(Rappture::zaxis) - xfield.rangeMin(Rappture::zaxis);
665
666    double dmin;
667    dmin = pow((lx*ly*lz)/((nSamples-1)*(nSamples-1)*(nSamples-1)), 0.333);
668
669    /* Recompute new number of points for each axis. */
670    _xNum = (size_t)ceil(lx/dmin);
671    _yNum = (size_t)ceil(ly/dmin);
672    _zNum = (size_t)ceil(lz/dmin);
673
674#ifndef HAVE_NPOT_TEXTURES
675    // must be an even power of 2 for older cards
676    _xNum = (int)pow(2.0, ceil(log10((double)_xNum)/log10(2.0)));
677    _yNum = (int)pow(2.0, ceil(log10((double)_yNum)/log10(2.0)));
678    _zNum = (int)pow(2.0, ceil(log10((double)_zNum)/log10(2.0)));
679#endif
680
681    dx = lx/(double)(_xNum-1);
682    dy = ly/(double)(_yNum-1);
683    dz = lz/(double)(_zNum-1);
684
685    TRACE("lx:%lf ly:%lf lz:%lf dmin:%lf dx:%lf dy:%lf dz:%lf\n", lx, ly, lz, dmin, dx, dy, dz);
686
687    size_t n = _nComponents * _xNum * _yNum * _zNum;
688    _values = (float *)realloc(_values, sizeof(float) * n);
689    memset(_values, 0, sizeof(float) * n);
690
691    // Generate the uniformly sampled rectangle that we need for a volume
692    float *destPtr = _values;
693    for (size_t i = 0; i < _zNum; i++) {
694        double z;
695
696        z = _zMin + (i * dx);
697        for (size_t j = 0; j < _yNum; j++) {
698            double y;
699               
700            y = _yMin + (j * dy);
701            for (size_t k = 0; k < _xNum; k++) {
702                double x;
703
704                x = _xMin + (k * dz);
705                destPtr[0] = xfield.value(x, y, z);
706                destPtr[1] = yfield.value(x, y, z);
707                destPtr[2] = zfield.value(x, y, z);
708            }
709        }
710    }
711    _nValues = _xNum * _yNum * _zNum * _nComponents;
712    return true;
713}
714
715void
716Rappture::Unirect3d::getVectorRange()
717{
718    assert(_nComponents == 3);
719    _magMax = -DBL_MAX, _magMin = DBL_MAX;
720    size_t i;
721    for (i = 0; i < _nValues; i += _nComponents) {
722        double vx, vy, vz, vm;
723
724        vx = _values[i];
725        vy = _values[i+1];
726        vz = _values[i+2];
727                   
728        vm = sqrt(vx*vx + vy*vy + vz*vz);
729        if (vm > _magMax) {
730            _magMax = vm;
731        }
732        if (vm < _magMin) {
733            _magMin = vm;
734        }
735    }
736    TRACE("getVectorRange: %g %g\n", _magMin, _magMax);
737}
738
739bool
740Rappture::Unirect3d::convert(Rappture::Unirect2d *dataPtr)
741{
742    _initialized = false;
743
744    _xValueMin = dataPtr->xValueMin();
745    _yValueMin = dataPtr->yValueMin();
746    _zValueMin = 0.0;
747    _xMin = dataPtr->xMin();
748    _yMin = dataPtr->yMin();
749    _zMin = 0.0;
750    _xMax = dataPtr->xMax();
751    _yMax = dataPtr->yMax();
752    _zMax = 0.0;
753    _xNum = dataPtr->xNum();
754    _yNum = dataPtr->yNum();
755    _zNum = 1;
756    switch (dataPtr->nComponents()) {
757    case 1:
758    case 3:
759        _values = (float *)realloc(_values, sizeof(float) * dataPtr->nValues());
760        if (_values == NULL) {
761            return false;
762        }
763        memcpy(_values, dataPtr->values(), dataPtr->nValues());
764        _nValues = dataPtr->nValues();
765        _nComponents = dataPtr->nComponents();
766        break;
767    case 2:
768        float *values;
769        _nValues = 3 * _xNum * _yNum * _zNum;
770        _values = (float *)realloc(_values, sizeof(float) * _nValues);
771        if (_values == NULL) {
772            return false;
773        }
774        values = dataPtr->values();
775        size_t i, j;
776        for (j = i = 0; i < dataPtr->nValues(); i += 2, j+= 3) {
777            _values[j] = values[i];
778            _values[j+1] = values[i+1];
779            _values[j+2] = 0.0f;
780        }
781        _nComponents = 3;
782        break;
783    }
784    return true;
785}
Note: See TracBrowser for help on using the repository browser.