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

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