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

Last change on this file since 4830 was 4612, checked in by ldelgass, 10 years ago

merge r3597 from trunk

  • 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
18static inline char *
19skipspaces(char *string)
20{
21    while (isspace(*string)) {
22        string++;
23    }
24    return string;
25}
26
27static inline char *
28getline(char **stringPtr, char *endPtr)
29{
30    char *line, *p;
31
32    line = skipspaces(*stringPtr);
33    for (p = line; p < endPtr; p++) {
34        if (*p == '\n') {
35            *p++ = '\0';
36            *stringPtr = p;
37            return line;
38        }
39    }
40    *stringPtr = p;
41    return line;
42}
43
44int
45Rappture::Unirect2d::parseBuffer(Tcl_Interp *interp, Rappture::Buffer &buf)
46{
47    Tcl_Obj *objPtr;
48    objPtr = Tcl_NewStringObj(buf.bytes(), buf.size());
49    Tcl_Obj **objv;
50    int objc;
51    if (Tcl_ListObjGetElements(interp, objPtr, &objc, &objv) != TCL_OK) {
52        return TCL_ERROR;
53    }
54    int result;
55    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.