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