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

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