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

Last change on this file since 2877 was 2877, checked in by ldelgass, 12 years ago

Some minor refactoring, also add some more fine grained config.h defines
(e.g. replace NV40 define with feature defines). Add tests for some required
OpenGL extensions (should always check for extensions or base version before
calling entry points from the extension). Also, clamp diffuse and specular
values on input and warn when they are out of range.
  • Property svn:eol-style set to native
File size: 25.3 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            TRACE("#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 HAVE_NPOT_TEXTURES
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.