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

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