source: nanovis/branches/1.1/Unirect.cpp @ 4891

Last change on this file since 4891 was 4877, checked in by ldelgass, 9 years ago

some cleanups
  • 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 = Tcl_NewStringObj(buf.bytes(), buf.size());
67    Tcl_Obj **objv;
68    int objc;
69    if (Tcl_ListObjGetElements(interp, objPtr, &objc, &objv) != TCL_OK) {
70        return TCL_ERROR;
71    }
72    int 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    const char *units[4];
87
88    if ((objc-1) & 0x01) {
89        Tcl_AppendResult(interp, Tcl_GetString(objv[0]), ": ",
90                "wrong # of arguments: should be key-value pairs",
91                (char *)NULL);
92        return TCL_ERROR;
93    }
94
95    /* Default order is  z, y, x. */
96    int axis1, axis2, axis3;
97    axis1 = 0;                  /* X-axis */
98    axis2 = 1;                  /* Y-axis */
99    axis3 = 2;                  /* Z-axis */
100
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    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 = (float *)realloc(_values, sizeof(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",
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    _nValues = nValues;
262    if (units[3] != NULL) {
263        if (_vUnits != NULL) {
264            free(_vUnits);
265        }
266        _vUnits = strdup(units[3]);
267    }
268    _xMin = min[axis3];
269    _xMax = max[axis3];
270    _xNum = num[axis3];
271    if (units[axis3] != NULL) {
272        if (_xUnits != NULL) {
273            free(_xUnits);
274        }
275        _xUnits = strdup(units[axis3]);
276    }
277    _yMin = min[axis2];
278    _yMax = max[axis2];
279    _yNum = num[axis2];
280    if (units[axis2] != NULL) {
281        if (_yUnits != NULL) {
282            free(_yUnits);
283        }
284        _yUnits = strdup(units[axis2]);
285    }
286    _zMin = min[axis1];
287    _zMax = max[axis1];
288    _zNum = num[axis1];
289    if (units[axis1] != NULL) {
290        if (_zUnits != NULL) {
291            free(_zUnits);
292        }
293        _zUnits = strdup(units[axis1]);
294    }
295    _initialized = true;
296#ifdef notdef
297    {
298        FILE *f;
299        f = fopen("/tmp/unirect3d.txt", "w");
300        fprintf(f, "unirect3d xmin %g xmax %g xnum %d ", _xMin, _xMax, _xNum);
301        fprintf(f, "ymin %g ymax %g ynum %d ", _yMin, _yMax, _yNum);
302        fprintf(f, "zmin %g zmax %g znum %d ", _zMin, _zMax, _zNum);
303        fprintf(f, "components %d values {\n",  _nComponents);
304        for (size_t i = 0; i < _nValues; i+= 3) {
305            fprintf(f, "%g %g %g\n", _values[i], _values[i+1], _values[i+2]);
306        }
307        fprintf(f, "}\n");
308        fclose(f);
309    }
310#endif
311    return TCL_OK;
312}
313
314int
315Rappture::Unirect2d::loadData(Tcl_Interp *interp, int objc,
316                              Tcl_Obj *const *objv)
317{
318    if ((objc & 0x01) == 0) {
319        Tcl_AppendResult(interp, Tcl_GetString(objv[0]), ": ",
320                "wrong number of arguments: should be key-value pairs",
321                (char *)NULL);
322        return TCL_ERROR;
323    }
324
325    int axis[2];
326    axis[0] = 1;                         /* X-axis */
327    axis[1] = 0;                        /* Y-axis */
328
329    _xNum = _yNum = _nValues = 0;
330    _xMin = _yMin = _xMax = _yMax = 0.0f;
331    if (_xUnits != NULL) {
332        free(_xUnits);
333    }
334    if (_yUnits != NULL) {
335        free(_yUnits);
336    }
337    if (_vUnits != NULL) {
338        free(_vUnits);
339    }
340    _xUnits = _yUnits = _vUnits = NULL;
341    _nValues = 0;
342
343    int i;
344    for (i = 1; i < objc; i += 2) {
345        const char *string;
346        char c;
347
348        string = Tcl_GetString(objv[i]);
349        c = string[0];
350        if ((c == 'x') && (strcmp(string, "xmin") == 0)) {
351            if (GetFloatFromObj(interp, objv[i+1], &_xMin) != TCL_OK) {
352                return TCL_ERROR;
353            }
354        } else if ((c == 'x') && (strcmp(string, "xmax") == 0)) {
355            if (GetFloatFromObj(interp, objv[i+1], &_xMax) != TCL_OK) {
356                return TCL_ERROR;
357            }
358        } else if ((c == 'x') && (strcmp(string, "xnum") == 0)) {
359            int n;
360            if (Tcl_GetIntFromObj(interp, objv[i+1], &n) != TCL_OK) {
361                return TCL_ERROR;
362            }
363            if (n <= 0) {
364                Tcl_AppendResult(interp, "bad xnum value: must be > 0",
365                     (char *)NULL);
366                return TCL_ERROR;
367            }
368            _xNum = n;
369        } else if ((c == 'x') && (strcmp(string, "xunits") == 0)) {
370            _xUnits = strdup(Tcl_GetString(objv[i+1]));
371        } else if ((c == 'y') && (strcmp(string, "ymin") == 0)) {
372            if (GetFloatFromObj(interp, objv[i+1], &_yMin) != TCL_OK) {
373                return TCL_ERROR;
374            }
375        } else if ((c == 'y') && (strcmp(string, "ymax") == 0)) {
376            if (GetFloatFromObj(interp, objv[i+1], &_yMax) != TCL_OK) {
377                return TCL_ERROR;
378            }
379        } else if ((c == 'y') && (strcmp(string, "ynum") == 0)) {
380            int n;
381            if (Tcl_GetIntFromObj(interp, objv[i+1], &n) != TCL_OK) {
382                return TCL_ERROR;
383            }
384            if (n <= 0) {
385                Tcl_AppendResult(interp, "bad ynum value: must be > 0",
386                                 (char *)NULL);
387                return TCL_ERROR;
388            }
389            _yNum = n;
390        } else if ((c == 'y') && (strcmp(string, "yunits") == 0)) {
391            _yUnits = strdup(Tcl_GetString(objv[i+1]));
392        } else if ((c == 'v') && (strcmp(string, "values") == 0)) {
393            Tcl_Obj **vobj;
394            int n;
395
396            Tcl_IncrRefCount(objv[i+1]);
397            if (Tcl_ListObjGetElements(interp, objv[i+1], &n, &vobj) != TCL_OK){
398                return TCL_ERROR;
399            }
400            if (n <= 0) {
401                Tcl_AppendResult(interp, "empty values list : must be > 0",
402                                 (char *)NULL);
403                return TCL_ERROR;
404            }
405            _nValues = n;
406            _values = (float *)realloc(_values, sizeof(float) * _nValues);
407            size_t j;
408            for (j = 0; j < _nValues; j++) {
409                if (GetFloatFromObj(interp, vobj[j], _values + j)!=TCL_OK) {
410                    return TCL_ERROR;
411                }
412            }
413            Tcl_DecrRefCount(objv[i+1]);
414        } else if ((c == 'u') && (strcmp(string, "units") == 0)) {
415            _vUnits = strdup(Tcl_GetString(objv[i+1]));
416        } else if ((c == 'c') && (strcmp(string, "components") == 0)) {
417            int n;
418
419            if (Tcl_GetIntFromObj(interp, objv[i+1], &n) != TCL_OK) {
420                return TCL_ERROR;
421            }
422            if (n <= 0) {
423                Tcl_AppendResult(interp, "bad extents value: must be > 0",
424                                 (char *)NULL);
425                return TCL_ERROR;
426            }
427            _nComponents = n;
428        } else if ((c == 'a') && (strcmp(string, "axisorder") == 0)) {
429            Tcl_Obj **order;
430            int n;
431
432            if (Tcl_ListObjGetElements(interp, objv[i+1], &n, &order)
433                != TCL_OK) {
434                return TCL_ERROR;
435            }
436            if (n != 2) {
437                Tcl_AppendResult(interp,
438                        "wrong # of axes defined for ordering the data",
439                        (char *)NULL);
440                return TCL_ERROR;
441            }
442            if ((GetAxisFromObj(interp, order[0], axis) != TCL_OK) ||
443                (GetAxisFromObj(interp, order[1], axis+1) != TCL_OK)) {
444                return TCL_ERROR;
445            }
446        } else {
447            Tcl_AppendResult(interp, "unknown key \"", string,
448                (char *)NULL);
449            return TCL_ERROR;
450        }
451    }
452    if (_nValues == 0) {
453        Tcl_AppendResult(interp, "missing \"values\" key", (char *)NULL);
454        return TCL_ERROR;
455    }
456    if (_nValues != (_xNum * _yNum * _nComponents)) {
457        Tcl_AppendResult(interp,
458                "wrong number of values: must be xnum*ynum*components",
459                         (char *)NULL);
460        return TCL_ERROR;
461    }
462   
463    if ((axis[0] != 1) || (axis[1] != 0)) {
464        TRACE("reordering data");
465        /* Reorder the data into x, y where x varies fastest and so on. */
466        size_t y;
467        float *dp;
468
469        dp = _values = (float *)realloc(_values, sizeof(float) * _nValues);
470        for (y = 0; y < _yNum; y++) {
471            size_t x;
472
473            for (x = 0; x < _xNum; x++) {
474                size_t i, v;
475                   
476                /* Compute the index from the data's described ordering. */
477                i = (y + (_yNum * x)) * _nComponents;
478                for(v = 0; v < _nComponents; v++) {
479                    dp[v] = _values[i+v];
480                }
481                dp += _nComponents;
482            }
483        }
484    }
485    _initialized = true;
486    return TCL_OK;
487}
488
489bool
490Rappture::Unirect3d::importDx(Rappture::Outcome &result, size_t nComponents,
491                              size_t length, char *string)
492{
493    int nx, ny, nz, npts;
494    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
495    char *p, *pend;
496
497    dx = dy = dz = 0.0;                 /* Suppress compiler warning. */
498    x0 = y0 = z0 = 0.0;                 /* May not have an origin line. */
499    nx = ny = nz = npts = 0;            /* Suppress compiler warning. */
500    for (p = string, pend = p + length; p < pend; /*empty*/) {
501        char *line;
502
503        line = getline(&p, pend);
504        if (line == pend) {
505            break;                        /* EOF */
506        }
507        if ((line[0] == '#') || (line == '\0')) {
508            continue;                        /* Skip blank or comment lines. */
509        }
510        if (sscanf(line, "object %*d class gridpositions counts %d %d %d",
511                   &nx, &ny, &nz) == 3) {
512            if ((nx < 0) || (ny < 0) || (nz < 0)) {
513                result.addError("invalid grid size: x=%d, y=%d, z=%d",
514                        nx, ny, nz);
515                return false;
516            }
517        } else if (sscanf(line, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
518            /* Found origin. */
519        } else if (sscanf(line, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
520            /* Found one of the delta lines. */
521            if (ddx != 0.0) {
522                dx = ddx;
523            } else if (ddy != 0.0) {
524                dy = ddy;
525            } else if (ddz != 0.0) {
526                dz = ddz;
527            }
528        } else if (sscanf(line, "object %*d class array type %*s shape 3"
529                " rank 1 items %d data follows", &npts) == 1) {
530            if (npts < 0) {
531                result.addError("bad # points %d", npts);
532                return false;
533            }
534            TRACE("#points=%d", npts);
535            if (npts != nx*ny*nz) {
536                result.addError("inconsistent data: expected %d points"
537                                " but found %d points", nx*ny*nz, npts);
538                return false;
539            }
540            break;
541        } else if (sscanf(line, "object %*d class array type %*s rank 0"
542                " times %d data follows", &npts) == 1) {
543            if (npts != nx*ny*nz) {
544                result.addError("inconsistent data: expected %d points"
545                                " but found %d points", nx*ny*nz, npts);
546                return false;
547            }
548            break;
549        }
550    }
551    if (npts != nx*ny*nz) {
552        result.addError("inconsistent data: expected %d points"
553                        " but found %d points", nx*ny*nz, npts);
554        return false;
555    }
556
557    _initialized = false;
558    _xValueMin = _yValueMin = _zValueMin = FLT_MAX;
559    _xValueMax = _yValueMax = _zValueMax = -FLT_MAX;
560    _xMin = x0, _yMin = y0, _zMin = z0;
561    _xNum = nx, _yNum = ny, _zNum = nz;
562    _xMax = _xMin + dx * _xNum;
563    _yMax = _yMin + dy * _yNum;
564    _zMax = _zMin + dz * _zNum;
565    _nComponents = nComponents;
566
567    _values = (float *)realloc(_values, sizeof(float) * npts * _nComponents);
568    _nValues = 0;
569    for (size_t ix = 0; ix < _xNum; ix++) {
570        for (size_t iy = 0; iy < _yNum; iy++) {
571            for (size_t iz = 0; iz < _zNum; iz++) {
572                char *line;
573                if ((p == pend) || (_nValues > (size_t)npts)) {
574                    break;
575                }
576                line = getline(&p, pend);
577                if ((line[0] == '#') || (line[0] == '\0')) {
578                    continue;                /* Skip blank or comment lines. */
579                }
580                double vx, vy, vz;
581                if (sscanf(line, "%lg %lg %lg", &vx, &vy, &vz) == 3) {
582                    int nindex = (iz*nx*ny + iy*nx + ix) * 3;
583                    if (vx < _xValueMin) {
584                        _xValueMin = vx;
585                    } else if (vx > _xValueMax) {
586                        _xValueMax = vx;
587                    }
588                    if (vy < _yValueMin) {
589                        _yValueMin = vy;
590                    } else if (vy > _yValueMax) {
591                        _yValueMax = vy;
592                    }
593                    if (vz < _zValueMin) {
594                        _zValueMin = vz;
595                    } else if (vz > _zValueMax) {
596                        _zValueMax = vz;
597                    }
598                    _values[nindex] = vx;
599                    _values[nindex+1] = vy;
600                    _values[nindex+2] = vz;
601                    _nValues++;
602                }
603            }
604        }
605    }
606    /* Make sure that we read all of the expected points. */
607    if (_nValues != (size_t)npts) {
608        result.addError("inconsistent data: expected %d points"
609                        " but found %d points", npts, _nValues);
610        free(_values);
611        _values = NULL;
612        return false;
613    }
614    _nValues *= _nComponents;
615    _initialized = true;
616#ifdef notdef
617    {
618        FILE *f;
619        f = fopen("/tmp/dx.txt", "w");
620        fprintf(f, "unirect3d xmin %g xmax %g xnum %d ", _xMin, _xMax, _xNum);
621        fprintf(f, "ymin %g ymax %g ynum %d ", _yMin, _yMax, _yNum);
622        fprintf(f, "zmin %g zmax %g znum %d ", _zMin, _zMax, _zNum);
623        fprintf(f, "components %d values {\n",  _nComponents);
624        for (size_t i = 0; i < _nValues; i+= 3) {
625            fprintf(f, "%g %g %g\n", _values[i], _values[i+1], _values[i+2]);
626        }
627        fprintf(f, "}\n");
628        fclose(f);
629    }
630#endif
631    return true;
632}
633
634bool
635Rappture::Unirect3d::resample(Rappture::Outcome &result, size_t nSamples)
636{
637    Rappture::Mesh1D xgrid(_xMin, _xMax, _xNum);
638    Rappture::Mesh1D ygrid(_yMin, _yMax, _yNum);
639    Rappture::Mesh1D zgrid(_zMin, _zMax, _zNum);
640    Rappture::FieldRect3D xfield(xgrid, ygrid, zgrid);
641    Rappture::FieldRect3D yfield(xgrid, ygrid, zgrid);
642    Rappture::FieldRect3D zfield(xgrid, ygrid, zgrid);
643
644    size_t i, j;
645    for (i = 0, j = 0; i < _nValues; i += _nComponents, j++) {
646        double vx, vy, vz;
647
648        vx = _values[i];
649        vy = _values[i+1];
650        vz = _values[i+2];
651       
652        xfield.define(j, vx);
653        yfield.define(j, vy);
654        zfield.define(j, vz);
655    }
656    // Figure out a good mesh spacing
657    double dx, dy, dz;
658    double lx, ly, lz;
659    lx = xfield.rangeMax(Rappture::xaxis) - xfield.rangeMin(Rappture::xaxis);
660    ly = xfield.rangeMax(Rappture::yaxis) - xfield.rangeMin(Rappture::yaxis);
661    lz = xfield.rangeMax(Rappture::zaxis) - xfield.rangeMin(Rappture::zaxis);
662
663    double dmin;
664    dmin = pow((lx*ly*lz)/((nSamples-1)*(nSamples-1)*(nSamples-1)), 0.333);
665
666    /* Recompute new number of points for each axis. */
667    _xNum = (size_t)ceil(lx/dmin);
668    _yNum = (size_t)ceil(ly/dmin);
669    _zNum = (size_t)ceil(lz/dmin);
670
671#ifndef HAVE_NPOT_TEXTURES
672    // must be an even power of 2 for older cards
673    _xNum = (int)pow(2.0, ceil(log10((double)_xNum)/log10(2.0)));
674    _yNum = (int)pow(2.0, ceil(log10((double)_yNum)/log10(2.0)));
675    _zNum = (int)pow(2.0, ceil(log10((double)_zNum)/log10(2.0)));
676#endif
677
678    dx = lx/(double)(_xNum-1);
679    dy = ly/(double)(_yNum-1);
680    dz = lz/(double)(_zNum-1);
681
682    TRACE("lx:%lf ly:%lf lz:%lf dmin:%lf dx:%lf dy:%lf dz:%lf", lx, ly, lz, dmin, dx, dy, dz);
683
684    size_t n = _nComponents * _xNum * _yNum * _zNum;
685    _values = (float *)realloc(_values, sizeof(float) * n);
686    memset(_values, 0, sizeof(float) * n);
687
688    // Generate the uniformly sampled rectangle that we need for a volume
689    float *destPtr = _values;
690    for (size_t i = 0; i < _zNum; i++) {
691        double z;
692
693        z = _zMin + (i * dx);
694        for (size_t j = 0; j < _yNum; j++) {
695            double y;
696               
697            y = _yMin + (j * dy);
698            for (size_t k = 0; k < _xNum; k++) {
699                double x;
700
701                x = _xMin + (k * dz);
702                destPtr[0] = xfield.value(x, y, z);
703                destPtr[1] = yfield.value(x, y, z);
704                destPtr[2] = zfield.value(x, y, z);
705            }
706        }
707    }
708    _nValues = _xNum * _yNum * _zNum * _nComponents;
709    return true;
710}
711
712void
713Rappture::Unirect3d::getVectorRange()
714{
715    assert(_nComponents == 3);
716    _magMax = -DBL_MAX, _magMin = DBL_MAX;
717    size_t i;
718    for (i = 0; i < _nValues; i += _nComponents) {
719        double vx, vy, vz, vm;
720
721        vx = _values[i];
722        vy = _values[i+1];
723        vz = _values[i+2];
724                   
725        vm = sqrt(vx*vx + vy*vy + vz*vz);
726        if (vm > _magMax) {
727            _magMax = vm;
728        }
729        if (vm < _magMin) {
730            _magMin = vm;
731        }
732    }
733    TRACE("getVectorRange: %g %g", _magMin, _magMax);
734}
735
736bool
737Rappture::Unirect3d::convert(Rappture::Unirect2d *dataPtr)
738{
739    _initialized = false;
740
741    _xValueMin = dataPtr->xValueMin();
742    _yValueMin = dataPtr->yValueMin();
743    _zValueMin = 0.0;
744    _xMin = dataPtr->xMin();
745    _yMin = dataPtr->yMin();
746    _zMin = 0.0;
747    _xMax = dataPtr->xMax();
748    _yMax = dataPtr->yMax();
749    _zMax = 0.0;
750    _xNum = dataPtr->xNum();
751    _yNum = dataPtr->yNum();
752    _zNum = 1;
753    switch (dataPtr->nComponents()) {
754    case 1:
755    case 3:
756        _values = (float *)realloc(_values, sizeof(float) * dataPtr->nValues());
757        if (_values == NULL) {
758            return false;
759        }
760        memcpy(_values, dataPtr->values(), dataPtr->nValues());
761        _nValues = dataPtr->nValues();
762        _nComponents = dataPtr->nComponents();
763        break;
764    case 2:
765        float *values;
766        _nValues = 3 * _xNum * _yNum * _zNum;
767        _values = (float *)realloc(_values, sizeof(float) * _nValues);
768        if (_values == NULL) {
769            return false;
770        }
771        values = dataPtr->values();
772        size_t i, j;
773        for (j = i = 0; i < dataPtr->nValues(); i += 2, j+= 3) {
774            _values[j] = values[i];
775            _values[j+1] = values[i+1];
776            _values[j+2] = 0.0f;
777        }
778        _nComponents = 3;
779        break;
780    }
781    return true;
782}
Note: See TracBrowser for help on using the repository browser.