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

Last change on this file since 3502 was 3502, checked in by ldelgass, 7 years ago

Add basic VTK structured points reader to nanovis, update copyright dates.

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