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

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