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

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

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