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

Last change on this file since 2831 was 2831, checked in by ldelgass, 13 years ago

Refactor texture classes, misc. cleanups, cut down on header pollution -- still
need to fix header deps in Makefile.in
  • Property svn:eol-style set to native
File size: 22.4 KB
Line 
1/* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */
2#include <float.h>
3#include <tcl.h>
4
5#include <RpField1D.h>
6#include <RpFieldRect3D.h>
7
8#include "Unirect.h"
9#include "Trace.h"
10
11extern int GetFloatFromObj(Tcl_Interp *interp, Tcl_Obj *objPtr,
12                           float *valuePtr);
13
14extern int GetAxisFromObj(Tcl_Interp *interp, Tcl_Obj *objPtr, int *indexPtr);
15
16static INLINE char *
17skipspaces(char *string)
18{
19    while (isspace(*string)) {
20        string++;
21    }
22    return string;
23}
24
25static INLINE char *
26getline(char **stringPtr, char *endPtr)
27{
28    char *line, *p;
29
30    line = skipspaces(*stringPtr);
31    for (p = line; p < endPtr; p++) {
32        if (*p == '\n') {
33            *p++ = '\0';
34            *stringPtr = p;
35            return line;
36        }
37    }
38    *stringPtr = p;
39    return line;
40}
41
42int
43Rappture::Unirect2d::ParseBuffer(Tcl_Interp *interp, Rappture::Buffer &buf)
44{
45    Tcl_Obj *objPtr;
46    objPtr = Tcl_NewStringObj(buf.bytes(), buf.size());
47    Tcl_Obj **objv;
48    int objc;
49    if (Tcl_ListObjGetElements(interp, objPtr, &objc, &objv) != TCL_OK) {
50        return TCL_ERROR;
51    }
52    int result;
53    result = LoadData(interp, objc, objv);
54    Tcl_DecrRefCount(objPtr);
55    if ((result != TCL_OK) || (!isInitialized())) {
56        return TCL_ERROR;
57    }
58    return TCL_OK;
59}
60
61int
62Rappture::Unirect3d::ParseBuffer(Tcl_Interp *interp, Rappture::Buffer &buf)
63{
64    Tcl_Obj *objPtr;
65    objPtr = Tcl_NewStringObj(buf.bytes(), buf.size());
66    Tcl_Obj **objv;
67    int objc;
68    if (Tcl_ListObjGetElements(interp, objPtr, &objc, &objv) != TCL_OK) {
69        return TCL_ERROR;
70    }
71    int result;
72    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    float *values;
87    const char *units[4], *order;
88
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    values = NULL;
104    num[0] = num[1] = num[2] = nValues = 0;
105    min[0] = min[1] = min[2] = max[0] = max[1] = max[2] = 0.0f;
106    order = units[0] = units[1] = units[2] = units[3] = NULL;
107
108    int i;
109    for (i = 1; i < objc; i += 2) {
110        const char *string;
111        char c;
112
113        string = Tcl_GetString(objv[i]);
114        c = string[0];
115        if ((c == 'x') && (strcmp(string, "xmin") == 0)) {
116            if (GetFloatFromObj(interp, objv[i+1], min+2) != TCL_OK) {
117                return TCL_ERROR;
118            }
119        } else if ((c == 'x') && (strcmp(string, "xmax") == 0)) {
120            if (GetFloatFromObj(interp, objv[i+1], max+2) != TCL_OK) {
121                return TCL_ERROR;
122            }
123        } else if ((c == 'x') && (strcmp(string, "xnum") == 0)) {
124            if (Tcl_GetIntFromObj(interp, objv[i+1], num+2) != TCL_OK) {
125                return TCL_ERROR;
126            }
127            if (num[2] <= 0) {
128                Tcl_AppendResult(interp, "bad xnum value: must be > 0",
129                     (char *)NULL);
130                return TCL_ERROR;
131            }
132        } else if ((c == 'x') && (strcmp(string, "xunits") == 0)) {
133            units[2] = Tcl_GetString(objv[i+1]);
134        } else if ((c == 'y') && (strcmp(string, "ymin") == 0)) {
135            if (GetFloatFromObj(interp, objv[i+1], min+1) != TCL_OK) {
136                return TCL_ERROR;
137            }
138        } else if ((c == 'y') && (strcmp(string, "ymax") == 0)) {
139            if (GetFloatFromObj(interp, objv[i+1], max+1) != TCL_OK) {
140                return TCL_ERROR;
141            }
142        } else if ((c == 'y') && (strcmp(string, "ynum") == 0)) {
143            if (Tcl_GetIntFromObj(interp, objv[i+1], num+1) != TCL_OK) {
144                return TCL_ERROR;
145            }
146            if (num[1] <= 0) {
147                Tcl_AppendResult(interp, "bad ynum value: must be > 0",
148                                 (char *)NULL);
149                return TCL_ERROR;
150            }
151        } else if ((c == 'y') && (strcmp(string, "yunits") == 0)) {
152            units[1] = Tcl_GetString(objv[i+1]);
153        } else if ((c == 'z') && (strcmp(string, "zmin") == 0)) {
154            if (GetFloatFromObj(interp, objv[i+1], min) != TCL_OK) {
155                return TCL_ERROR;
156            }
157        } else if ((c == 'z') && (strcmp(string, "zmax") == 0)) {
158            if (GetFloatFromObj(interp, objv[i+1], max) != TCL_OK) {
159                return TCL_ERROR;
160            }
161        } else if ((c == 'z') && (strcmp(string, "znum") == 0)) {
162            if (Tcl_GetIntFromObj(interp, objv[i+1], num) != TCL_OK) {
163                return TCL_ERROR;
164            }
165            if (num[0] <= 0) {
166                Tcl_AppendResult(interp, "bad znum value: must be > 0",
167                                 (char *)NULL);
168                return TCL_ERROR;
169            }
170        } else if ((c == 'z') && (strcmp(string, "zunits") == 0)) {
171            units[0] = Tcl_GetString(objv[i+1]);
172        } else if ((c == 'v') && (strcmp(string, "values") == 0)) {
173            Tcl_Obj **vobj;
174
175            if (Tcl_ListObjGetElements(interp, objv[i+1], &nValues, &vobj)
176                != TCL_OK) {
177                return TCL_ERROR;
178            }
179            _values = (float *)realloc(_values, sizeof(float) * nValues);
180            int j;
181            for (j = 0; j < nValues; j++) {
182                if (GetFloatFromObj(interp, vobj[j], _values + j)!=TCL_OK) {
183                    return TCL_ERROR;
184                }
185            }
186        } else if ((c == 'u') && (strcmp(string, "units") == 0)) {
187            _vUnits = strdup(Tcl_GetString(objv[i+1]));
188        } else if ((c == 'c') && (strcmp(string, "components") == 0)) {
189            int n;
190
191            if (Tcl_GetIntFromObj(interp, objv[i+1], &n) != TCL_OK) {
192                return TCL_ERROR;
193            }
194            if (n <= 0) {
195                Tcl_AppendResult(interp, "bad extents value: must be > 0",
196                                 (char *)NULL);
197                return TCL_ERROR;
198            }
199            _nComponents = n;
200        } else if ((c == 'a') && (strcmp(string, "axisorder") == 0)) {
201            Tcl_Obj **axes;
202            int n;
203
204            if (Tcl_ListObjGetElements(interp, objv[i+1], &n, &axes)
205                != TCL_OK) {
206                return TCL_ERROR;
207            }
208            if (n != 3) {
209                return TCL_ERROR;
210            }
211            if ((GetAxisFromObj(interp, axes[0], &axis1) != TCL_OK) ||
212                (GetAxisFromObj(interp, axes[1], &axis2) != TCL_OK) ||
213                (GetAxisFromObj(interp, axes[2], &axis3) != TCL_OK)) {
214                return TCL_ERROR;
215            }
216        } else {
217            Tcl_AppendResult(interp, "unknown key \"", string,
218                (char *)NULL);
219            return TCL_ERROR;
220        }
221    }
222    if (_values == NULL) {
223        Tcl_AppendResult(interp, "missing \"values\" key", (char *)NULL);
224        return TCL_ERROR;
225    }
226    if ((size_t)nValues != (num[0] * num[1] * num[2] * _nComponents)) {
227        TRACE("num[0]=%d num[1]=%d num[2]=%d ncomponents=%d nValues=%d\n",
228               num[0], num[1], num[2], _nComponents, nValues);
229        Tcl_AppendResult(interp,
230                "wrong # of values: must be xnum*ynum*znum*extents",
231                         (char *)NULL);
232       return TCL_ERROR;
233    }
234   
235#ifdef notdef
236    if ((axis1 != 0) || (axis2 != 1) || (axis3 != 2)) {
237        // Reorder the data into x, y, z where x varies fastest and so on.
238        int z;
239        float *data, *dp;
240
241        dp = data = new float[nValues];
242        for (z = 0; z < num[0]; z++) {
243            int y;
244
245            for (y = 0; y < num[1]; y++) {
246                int x;
247
248                for (x = 0; x < num[2]; x++) {
249                    int i;
250                   
251                    /* Compute the index from the data's described ordering. */
252                    i = ((z*num[axis2]*num[axis3]) + (y*num[axis3]) + x) * 3;
253                    for(size_t v = 0; v < _nComponents; v++) {
254                        dp[v] = values[i+v];
255                    }
256                    dp += _nComponents;
257                }
258            }
259        }
260        delete [] values;
261        values = data;
262    }
263#endif
264    _nValues = nValues;
265    if (units[3] != NULL) {
266        if (_vUnits != NULL) {
267            free(_vUnits);
268        }
269        _vUnits = strdup(units[3]);
270    }
271    _xMin = min[axis3];
272    _xMax = max[axis3];
273    _xNum = num[axis3];
274    if (units[axis3] != NULL) {
275        if (_xUnits != NULL) {
276            free(_xUnits);
277        }
278        _xUnits = strdup(units[axis3]);
279    }
280    _yMin = min[axis2];
281    _yMax = max[axis2];
282    _yNum = num[axis2];
283    if (units[axis2] != NULL) {
284        if (_yUnits != NULL) {
285            free(_yUnits);
286        }
287        _yUnits = strdup(units[axis2]);
288    }
289    _zMin = min[axis1];
290    _zMax = max[axis1];
291    _zNum = num[axis1];
292    if (units[axis1] != NULL) {
293        if (_zUnits != NULL) {
294            free(_zUnits);
295        }
296        _zUnits = strdup(units[axis1]);
297    }
298    _initialized = true;
299#ifdef notdef
300    {
301        FILE *f;
302        f = fopen("/tmp/unirect3d.txt", "w");
303        fprintf(f, "unirect3d xmin %g xmax %g xnum %d ", _xMin, _xMax, _xNum);
304        fprintf(f, "ymin %g ymax %g ynum %d ", _yMin, _yMax, _yNum);
305        fprintf(f, "zmin %g zmax %g znum %d ", _zMin, _zMax, _zNum);
306        fprintf(f, "components %d values {\n",  _nComponents);
307        for (size_t i = 0; i < _nValues; i+= 3) {
308            fprintf(f, "%g %g %g\n", _values[i], _values[i+1], _values[i+2]);
309        }
310        fprintf(f, "}\n");
311        fclose(f);
312    }
313#endif
314    return TCL_OK;
315}
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\n");
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
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            printf("#points=%d\n", 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
639
640bool
641Rappture::Unirect3d::Resample(Rappture::Outcome &result, size_t nSamples)
642{
643    Rappture::Mesh1D xgrid(_xMin, _xMax, _xNum);
644    Rappture::Mesh1D ygrid(_yMin, _yMax, _yNum);
645    Rappture::Mesh1D zgrid(_zMin, _zMax, _zNum);
646    Rappture::FieldRect3D xfield(xgrid, ygrid, zgrid);
647    Rappture::FieldRect3D yfield(xgrid, ygrid, zgrid);
648    Rappture::FieldRect3D zfield(xgrid, ygrid, zgrid);
649
650    size_t i, j;
651    for (i = 0, j = 0; i < _nValues; i += _nComponents, j++) {
652        double vx, vy, vz;
653
654        vx = _values[i];
655        vy = _values[i+1];
656        vz = _values[i+2];
657       
658        xfield.define(j, vx);
659        yfield.define(j, vy);
660        zfield.define(j, vz);
661    }
662    // Figure out a good mesh spacing
663    double dx, dy, dz;
664    dx = xfield.rangeMax(Rappture::xaxis) - xfield.rangeMin(Rappture::xaxis);
665    dy = xfield.rangeMax(Rappture::yaxis) - xfield.rangeMin(Rappture::yaxis);
666    dz = xfield.rangeMax(Rappture::zaxis) - xfield.rangeMin(Rappture::zaxis);
667
668    double dmin;
669    dmin = pow((dx*dy*dz)/(nSamples*nSamples*nSamples), 0.333);
670   
671    TRACE("dx:%lf dy:%lf dz:%lf dmin:%lf\n", dx, dy, dz, dmin);
672
673    /* Recompute new number of points for each axis. */
674    _xNum = (size_t)ceil(dx/dmin);
675    _yNum = (size_t)ceil(dy/dmin);
676    _zNum = (size_t)ceil(dz/dmin);
677   
678#ifndef NV40
679    // must be an even power of 2 for older cards
680    _xNum = (int)pow(2.0, ceil(log10((double)_xNum)/log10(2.0)));
681    _yNum = (int)pow(2.0, ceil(log10((double)_yNum)/log10(2.0)));
682    _zNum = (int)pow(2.0, ceil(log10((double)_zNum)/log10(2.0)));
683#endif
684
685    size_t n = _nComponents * _xNum * _yNum * _zNum;
686    _values = (float *)realloc(_values, sizeof(float) * n);
687    memset(_values, 0, sizeof(float) * n);
688   
689    // Generate the uniformly sampled rectangle that we need for a volume
690    float *destPtr = _values;
691    for (size_t i = 0; i < _zNum; i++) {
692        double z;
693
694        z = _zMin + (i * dmin);
695        for (size_t j = 0; j < _yNum; j++) {
696            double y;
697               
698            y = _yMin + (j * dmin);
699            for (size_t k = 0; k < _xNum; k++) {
700                double x;
701
702                x = _xMin + (k * dmin);
703                destPtr[0] = xfield.value(x, y, z);
704                destPtr[1] = yfield.value(x, y, z);
705                destPtr[2] = zfield.value(x, y, z);
706            }
707        }
708    }
709    _nValues = _xNum * _yNum * _zNum * _nComponents;
710    return true;
711}
712
713
714void
715Rappture::Unirect3d::GetVectorRange(void)
716{
717    assert(_nComponents == 3);
718    TRACE("GetVectorRange\n");
719    _magMax = -DBL_MAX, _magMin = DBL_MAX;
720    size_t i;
721    for (i = 0; i < _nValues; i += _nComponents) {
722        double vx, vy, vz, vm;
723
724        vx = _values[i];
725        vy = _values[i+1];
726        vz = _values[i+2];
727                   
728        vm = sqrt(vx*vx + vy*vy + vz*vz);
729        if (vm > _magMax) {
730            _magMax = vm;
731        }
732        if (vm < _magMin) {
733            _magMin = vm;
734        }
735    }
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.