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

Last change on this file since 1431 was 1431, checked in by gah, 11 years ago

fixup new flow visualization command structure

File size: 17.9 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 = 2;                  /* Z-axis */
58    axis2 = 1;                  /* Y-axis */
59    axis3 = 0;                  /* X-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 == 'o') && (strcmp(string, "order") == 0)) {
147            Tcl_Obj **axes;
148            int n;
149
150            if (Tcl_ListObjGetElements(interp, objv[i+1], &n, &axes) 
151                != TCL_OK) {
152                return TCL_ERROR;
153            }
154            if (n != 3) {
155                return TCL_ERROR;
156            }
157            if ((GetAxisFromObj(interp, axes[0], &axis1) != TCL_OK) ||
158                (GetAxisFromObj(interp, axes[1], &axis2) != TCL_OK) ||
159                (GetAxisFromObj(interp, axes[2], &axis3) != TCL_OK)) {
160                return TCL_ERROR;
161            }
162        } else {
163            Tcl_AppendResult(interp, "unknown key \"", string,
164                (char *)NULL);
165            return TCL_ERROR;
166        }
167    }
168    if (values == NULL) {
169        Tcl_AppendResult(interp, "missing \"values\" key", (char *)NULL);
170        return TCL_ERROR;
171    }
172    if ((size_t)nValues != (num[0] * num[1] * num[2] * _nComponents)) {
173        Tcl_AppendResult(interp, 
174                "wrong number of values: must be xnum*ynum*znum*extents", 
175                         (char *)NULL);
176        return TCL_ERROR;
177    }
178   
179    if ((axis1 != 2) || (axis2 != 1) || (axis3 != 0)) {
180        // Reorder the data into x, y, z where x varies fastest and so on.
181        int z;
182        float *data, *dp;
183
184        dp = data = new float[nValues];
185        for (z = 0; z < num[0]; z++) {
186            int y;
187
188            for (y = 0; y < num[1]; y++) {
189                int x;
190
191                for (x = 0; x < num[2]; x++) {
192                    int i;
193                   
194                    /* Compute the index from the data's described ordering. */
195                    i = ((z*num[axis2]*num[axis3]) + (y*num[axis3]) + x) * 3;
196                    for(size_t v = 0; v < _nComponents; v++) {
197                        dp[v] = values[i+v];
198                    }
199                    dp += _nComponents;
200                }
201            }
202        }
203        delete [] values;
204        values = data;
205    }
206
207    _values = values;
208    _nValues = nValues;
209    if (units[3] != NULL) {
210        _vUnits = strdup(units[3]);
211    }
212    _xMin = min[axis3];
213    _xMax = max[axis3];
214    _xNum = num[axis3];
215    if (units[axis3] != NULL) {
216        _xUnits = strdup(units[axis3]);
217    }
218    _yMin = min[axis2];
219    _yMax = max[axis2];
220    _yNum = num[axis2];
221    if (units[axis2] != NULL) {
222        _yUnits = strdup(units[axis2]);
223    }
224    _zMin = min[axis1];
225    _zMax = max[axis1];
226    _zNum = num[axis1];
227    if (units[axis1] != NULL) {
228        _zUnits = strdup(units[axis1]);
229    }
230    _initialized = true;
231    return TCL_OK;
232}
233
234
235int 
236Rappture::Unirect2d::LoadData(Tcl_Interp *interp, int objc, 
237                              Tcl_Obj *const *objv)
238{
239
240    if ((objc & 0x01) == 0) {
241        Tcl_AppendResult(interp, Tcl_GetString(objv[0]), ": ",
242                "wrong number of arguments: should be key-value pairs",
243                (char *)NULL);
244        return TCL_ERROR;
245    }
246
247    int axis[2];
248    axis[0] = 1;                        /* X-axis */
249    axis[1] = 0;                        /* Y-axis */
250
251    _xNum = _yNum = _nValues = 0;
252    _xMin = _yMin = _xMax = _yMax = 0.0f;
253    if (_xUnits != NULL) {
254        free(_xUnits);
255    }
256    if (_yUnits != NULL) {
257        free(_yUnits);
258    }
259    if (_vUnits != NULL) {
260        free(_vUnits);
261    }
262    _xUnits = _yUnits = _vUnits = NULL;
263    if (_values != NULL) {
264        delete [] _values;
265    }
266    _values = NULL;
267
268    int i;
269    for (i = 1; i < objc; i += 2) {
270        const char *string;
271        char c;
272
273        string = Tcl_GetString(objv[i]);
274        c = string[0];
275        if ((c == 'x') && (strcmp(string, "xmin") == 0)) {
276            if (GetFloatFromObj(interp, objv[i+1], &_xMin) != TCL_OK) {
277                return TCL_ERROR;
278            }
279        } else if ((c == 'x') && (strcmp(string, "xmax") == 0)) {
280            if (GetFloatFromObj(interp, objv[i+1], &_xMax) != TCL_OK) {
281                return TCL_ERROR;
282            }
283        } else if ((c == 'x') && (strcmp(string, "xnum") == 0)) {
284            int n;
285            if (Tcl_GetIntFromObj(interp, objv[i+1], &n) != TCL_OK) {
286                return TCL_ERROR;
287            }
288            if (n <= 0) {
289                Tcl_AppendResult(interp, "bad xnum value: must be > 0",
290                     (char *)NULL);
291                return TCL_ERROR;
292            }
293            _xNum = n;
294        } else if ((c == 'x') && (strcmp(string, "xunits") == 0)) {
295            _xUnits = strdup(Tcl_GetString(objv[i+1]));
296        } else if ((c == 'y') && (strcmp(string, "ymin") == 0)) {
297            if (GetFloatFromObj(interp, objv[i+1], &_yMin) != TCL_OK) {
298                return TCL_ERROR;
299            }
300        } else if ((c == 'y') && (strcmp(string, "ymax") == 0)) {
301            if (GetFloatFromObj(interp, objv[i+1], &_yMax) != TCL_OK) {
302                return TCL_ERROR;
303            }
304        } else if ((c == 'y') && (strcmp(string, "ynum") == 0)) {
305            int n;
306            if (Tcl_GetIntFromObj(interp, objv[i+1], &n) != TCL_OK) {
307                return TCL_ERROR;
308            }
309            if (n <= 0) {
310                Tcl_AppendResult(interp, "bad ynum value: must be > 0",
311                                 (char *)NULL);
312                return TCL_ERROR;
313            }
314            _yNum = n;
315        } else if ((c == 'y') && (strcmp(string, "yunits") == 0)) {
316            _yUnits = strdup(Tcl_GetString(objv[i+1]));
317        } else if ((c == 'v') && (strcmp(string, "values") == 0)) {
318            Tcl_Obj **vobj;
319            int n;
320
321            if (Tcl_ListObjGetElements(interp, objv[i+1], &n, &vobj) != TCL_OK){
322                return TCL_ERROR;
323            }
324            if (n <= 0) {
325                Tcl_AppendResult(interp, "empty values list : must be > 0",
326                                 (char *)NULL);
327                return TCL_ERROR;
328            }
329            _nValues = n;
330            _values = new float[_nValues];
331            size_t j;
332            for (j = 0; j < _nValues; j++) {
333                if (GetFloatFromObj(interp, vobj[j], _values + j)!=TCL_OK) {
334                    return TCL_ERROR;
335                }
336            }
337        } else if ((c == 'u') && (strcmp(string, "units") == 0)) {
338            _vUnits = strdup(Tcl_GetString(objv[i+1]));
339        } else if ((c == 'e') && (strcmp(string, "extents") == 0)) {
340            int n;
341
342            if (Tcl_GetIntFromObj(interp, objv[i+1], &n) != TCL_OK) {
343                return TCL_ERROR;
344            }
345            if (n <= 0) {
346                Tcl_AppendResult(interp, "bad extents value: must be > 0",
347                                 (char *)NULL);
348                return TCL_ERROR;
349            }
350            _nComponents = n;
351        } else if ((c == 'a') && (strcmp(string, "axisorder") == 0)) {
352            Tcl_Obj **order;
353            int n;
354
355            if (Tcl_ListObjGetElements(interp, objv[i+1], &n, &order) 
356                != TCL_OK) {
357                return TCL_ERROR;
358            }
359            if (n != 2) {
360                Tcl_AppendResult(interp, 
361                        "wrong # of axes defined for ordering the data",
362                        (char *)NULL);
363                return TCL_ERROR;
364            }
365            if ((GetAxisFromObj(interp, order[0], axis) != TCL_OK) ||
366                (GetAxisFromObj(interp, order[1], axis+1) != TCL_OK)) {
367                return TCL_ERROR;
368            }
369        } else {
370            Tcl_AppendResult(interp, "unknown key \"", string,
371                (char *)NULL);
372            return TCL_ERROR;
373        }
374    }
375    if (_values == NULL) {
376        Tcl_AppendResult(interp, "missing \"values\" key", (char *)NULL);
377        return TCL_ERROR;
378    }
379    if (_nValues != (_xNum * _yNum * _nComponents)) {
380        Tcl_AppendResult(interp, 
381                "wrong number of values: must be xnum*ynum*extents", 
382                         (char *)NULL);
383        return TCL_ERROR;
384    }
385   
386#ifndef notdef
387    if ((axis[0] != 1) || (axis[1] != 0)) {
388        fprintf(stderr, "reordering data\n");
389        // Reorder the data into x, y where x varies fastest and so on.
390        size_t y;
391        float *data, *dp;
392
393        dp = data = new float[_nValues];
394        for (y = 0; y < _yNum; y++) {
395            size_t x;
396
397            for (x = 0; x < _xNum; x++) {
398                size_t i, v;
399                   
400                /* Compute the index from the data's described ordering. */
401                i = (y + (_yNum * x)) * _nComponents;
402                for(v = 0; v < _nComponents; v++) {
403                    dp[v] = _values[i+v];
404                }
405                dp += _nComponents;
406            }
407        }
408        delete [] _values;
409        _values = data;
410    }
411#endif
412    _initialized = true;
413    return TCL_OK;
414}
415
416
417bool
418Rappture::Unirect3d::ReadVectorDataFromDx(Rappture::Outcome &result, 
419        size_t length, char *string) 
420{
421    size_t nx, ny, nz, npts;
422    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
423    char *p, *endPtr;
424
425    dx = dy = dz = 0.0;         // Suppress compiler warning.
426    x0 = y0 = z0 = 0.0;         // May not have an origin line.
427    for (p = string, endPtr = p + length; p < endPtr; /*empty*/) {
428        char *line;
429
430        line = getline(&p, endPtr);
431        if (line == endPtr) {
432            break;
433        }
434        if ((line[0] == '#') || (line == '\0')) {
435            continue;           // Skip blank or comment lines.
436        }
437        if (sscanf(line, "object %*d class gridpositions counts %d %d %d", 
438                   &nx, &ny, &nz) == 3) {
439            if ((nx < 0) || (ny < 0) || (nz < 0)) {
440                result.addError("invalid grid size: x=%d, y=%d, z=%d", 
441                        nx, ny, nz);
442                return false;
443            }
444        } else if (sscanf(line, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
445            // found origin
446        } else if (sscanf(line, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
447            // found one of the delta lines
448            if (ddx != 0.0) {
449                dx = ddx;
450            } else if (ddy != 0.0) {
451                dy = ddy;
452            } else if (ddz != 0.0) {
453                dz = ddz;
454            }
455        } else if (sscanf(line, "object %*d class array type %*s shape 3"
456                " rank 1 items %d data follows", &npts) == 1) {
457            printf("#points=%d\n", npts);
458            if (npts != nx*ny*nz) {
459                result.addError("inconsistent data: expected %d points"
460                                " but found %d points", nx*ny*nz, npts);
461                return false;
462            }
463            break;
464        } else if (sscanf(line, "object %*d class array type %*s rank 0"
465                " times %d data follows", &npts) == 1) {
466            if (npts != nx*ny*nz) {
467                result.addError("inconsistent data: expected %d points"
468                                " but found %d points", nx*ny*nz, npts);
469                return false;
470            }
471            break;
472        }
473    }
474    if (npts != nx*ny*nz) {
475        result.addError("inconsistent data: expected %d points"
476                        " but found %d points", nx*ny*nz, npts);
477        return false;
478    }
479
480    _initialized = false;
481    _xValueMin = _yValueMin = _zValueMin = FLT_MAX;
482    _xValueMax = _yValueMax = _zValueMax = -FLT_MAX;
483    _xMin = x0, _yMin = y0, _zMin = z0;
484    _xNum = nx, _yNum = ny, _zNum = nz;
485    _xMax = _xMin + dx * _xNum;
486    _yMax = _yMin + dy * _yNum;
487    _zMax = _zMin + dz * _zNum;
488    if (_values != NULL) {
489        delete [] _values;
490    }
491    _values = new float[npts * _nComponents];
492    _nValues = 0;
493    for (size_t ix = 0; ix < _xNum; ix++) {
494        for (size_t iy = 0; iy < _yNum; iy++) {
495            for (size_t iz = 0; iz < _zNum; iz++) {
496                char *line;
497                if ((p == endPtr) || (_nValues > (size_t)npts)) {
498                    break;
499                }
500                line = getline(&p, endPtr);
501                if ((line[0] == '#') || (line[0] == '\0')) {
502                    continue;   // Skip blank or comment lines.
503                }
504                double vx, vy, vz;
505                if (sscanf(line, "%lg %lg %lg", &vx, &vy, &vz) == 3) {
506                    int nindex = (iz*nx*ny + iy*nx + ix) * 3;
507                    if (vx < _xValueMin) {
508                        _xValueMin = vx;
509                    } else if (vx > _xValueMax) {
510                        _xValueMax = vx;
511                    }
512                    if (vy < _yValueMin) {
513                        _yValueMin = vy;
514                    } else if (vy > _yValueMax) {
515                        _yValueMax = vy;
516                    }
517                    if (vz < _zValueMin) {
518                        _zValueMin = vz;
519                    } else if (vz > _zValueMax) {
520                        _zValueMax = vz;
521                    }
522                    _values[nindex] = vx;
523                    _values[nindex+1] = vy;
524                    _values[nindex+2] = vz;
525                    _nValues++;
526                }
527            }
528        }
529    }
530    // make sure that we read all of the expected points
531    if (_nValues != npts) {
532        result.addError("inconsistent data: expected %d points"
533                        " but found %d points", npts, _nValues);
534        delete []  _values;
535        _values = NULL;
536        return false;
537    }
538    _nValues *= _nComponents;
539    _initialized = true;
540    return true;
541}
542
543
544bool
545Rappture::Unirect3d::Resample(Rappture::Outcome &result, int nSamples)
546{
547    Rappture::Mesh1D xgrid(_xMin, _xMax, _xNum);
548    Rappture::Mesh1D ygrid(_yMin, _yMax, _yNum);
549    Rappture::Mesh1D zgrid(_zMin, _zMax, _zNum);
550    Rappture::FieldRect3D xfield(xgrid, ygrid, zgrid);
551    Rappture::FieldRect3D yfield(xgrid, ygrid, zgrid);
552    Rappture::FieldRect3D zfield(xgrid, ygrid, zgrid);
553
554    size_t i, j;
555    for (i = 0, j = 0; i < _nValues; i += _nComponents, j++) {
556        double vx, vy, vz;
557
558        vx = _values[i];
559        vy = _values[i+1];
560        vz = _values[i+2];
561       
562        xfield.define(j, vx);
563        yfield.define(j, vy);
564        zfield.define(j, vz);
565    }
566    // Figure out a good mesh spacing
567    double dx, dy, dz;
568    dx = xfield.rangeMax(Rappture::xaxis) - xfield.rangeMin(Rappture::xaxis);
569    dy = xfield.rangeMax(Rappture::yaxis) - xfield.rangeMin(Rappture::yaxis);
570    dz = xfield.rangeMax(Rappture::zaxis) - xfield.rangeMin(Rappture::zaxis);
571
572    double dmin;
573    dmin = pow((dx*dy*dz)/(nSamples*nSamples*nSamples), 0.333);
574   
575    printf("dx:%lf dy:%lf dz:%lf dmin:%lf\n", dx, dy, dz, dmin);
576
577    /* Recompute new number of points for each axis. */
578    _xNum = (size_t)ceil(dx/dmin);
579    _yNum = (size_t)ceil(dy/dmin);
580    _zNum = (size_t)ceil(dz/dmin);
581   
582#ifndef NV40
583    // must be an even power of 2 for older cards
584    _xNum = (int)pow(2.0, ceil(log10((double)_xNum)/log10(2.0)));
585    _yNum = (int)pow(2.0, ceil(log10((double)_yNum)/log10(2.0)));
586    _zNum = (int)pow(2.0, ceil(log10((double)_zNum)/log10(2.0)));
587#endif
588
589    size_t n = _nComponents * _xNum * _yNum * _zNum;
590    float *data = new float[n];
591    memset(data, 0, sizeof(float) * n);
592   
593    // Generate the uniformly sampled rectangle that we need for a volume
594    float *destPtr = data;
595    for (size_t i = 0; i < _zNum; i++) {
596        double z;
597
598        z = _zMin + (i * dmin);
599        for (size_t j = 0; j < _yNum; j++) {
600            double y;
601               
602            y = _yMin + (j * dmin);
603            for (size_t k = 0; k < _xNum; k++) {
604                double x;
605
606                x = _xMin + (k * dmin);
607                destPtr[0] = xfield.value(x, y, z);
608                destPtr[1] = yfield.value(x, y, z);
609                destPtr[2] = zfield.value(x, y, z);
610            }
611        }
612    }
613    delete [] _values;
614    _values = data;
615    _nValues = _xNum * _yNum * _zNum * _nComponents;
616    return true;
617}
618
Note: See TracBrowser for help on using the repository browser.