Ignore:
Timestamp:
May 11, 2009, 7:39:13 PM (11 years ago)
Author:
gah
Message:

Initial commit of new flow visualization command structure

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/packages/vizservers/nanovis/Unirect.cpp

    r1381 r1429  
    11
     2#include <float.h>
    23#include <tcl.h>
    34#include <Unirect.h>
     5#include "RpField1D.h"
     6#include "RpFieldRect3D.h"
    47
    58extern int GetFloatFromObj(Tcl_Interp *interp, Tcl_Obj *objPtr,
    69        float *valuePtr);
    710extern 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}
    837
    938int
     
    3059    axis3 = 0;                  /* X-axis */
    3160
    32     int extents = 1;
    3361    values = NULL;
    3462    num[0] = num[1] = num[2] = nValues = 0;
     
    115143            }
    116144        } else if ((c == 'u') && (strcmp(string, "units") == 0)) {
    117             vUnits_ = strdup(Tcl_GetString(objv[i+1]));
     145            _vUnits = strdup(Tcl_GetString(objv[i+1]));
    118146        } else if ((c == 'o') && (strcmp(string, "order") == 0)) {
    119147            Tcl_Obj **axes;
     
    142170        return TCL_ERROR;
    143171    }
    144     if (nValues != (num[0] * num[1] * num[2] * extents)) {
     172    if ((size_t)nValues != (num[0] * num[1] * num[2] * _nComponents)) {
    145173        Tcl_AppendResult(interp,
    146174                "wrong number of values: must be xnum*ynum*znum*extents",
     
    149177    }
    150178   
    151     fprintf(stderr, "generating %dx%dx%dx%d = %d points\n",
    152             num[0], num[1], num[2], extents, num[0]*num[1]*num[2]*extents);
    153 
    154179    if ((axis1 != 2) || (axis2 != 1) || (axis3 != 0)) {
    155180        // Reorder the data into x, y, z where x varies fastest and so on.
     
    165190
    166191                for (x = 0; x < num[2]; x++) {
    167                     int i, v;
     192                    int i;
    168193                   
    169194                    /* Compute the index from the data's described ordering. */
    170195                    i = ((z*num[axis2]*num[axis3]) + (y*num[axis3]) + x) * 3;
    171                     for(v = 0; v < extents; v++) {
     196                    for(size_t v = 0; v < _nComponents; v++) {
    172197                        dp[v] = values[i+v];
    173198                    }
    174                     dp += extents;
     199                    dp += _nComponents;
    175200                }
    176201            }
     
    180205    }
    181206
    182     values_ = values;
    183     nValues_ = nValues;
    184     extents_ = extents;
     207    _values = values;
     208    _nValues = nValues;
    185209    if (units[3] != NULL) {
    186         vUnits_ = strdup(units[3]);
    187     }
    188     xMin_ = min[axis3];
    189     xMax_ = max[axis3];
    190     xNum_ = num[axis3];
     210        _vUnits = strdup(units[3]);
     211    }
     212    _xMin = min[axis3];
     213    _xMax = max[axis3];
     214    _xNum = num[axis3];
    191215    if (units[axis3] != NULL) {
    192         xUnits_ = strdup(units[axis3]);
    193     }
    194     yMin_ = min[axis2];
    195     yMax_ = max[axis2];
    196     yNum_ = num[axis2];
     216        _xUnits = strdup(units[axis3]);
     217    }
     218    _yMin = min[axis2];
     219    _yMax = max[axis2];
     220    _yNum = num[axis2];
    197221    if (units[axis2] != NULL) {
    198         yUnits_ = strdup(units[axis2]);
    199     }
    200     zMin_ = min[axis1];
    201     zMax_ = max[axis1];
    202     zNum_ = num[axis1];
     222        _yUnits = strdup(units[axis2]);
     223    }
     224    _zMin = min[axis1];
     225    _zMax = max[axis1];
     226    _zNum = num[axis1];
    203227    if (units[axis1] != NULL) {
    204         zUnits_ = strdup(units[axis1]);
    205     }
    206     initialized_ = true;
     228        _zUnits = strdup(units[axis1]);
     229    }
     230    _initialized = true;
    207231    return TCL_OK;
    208232}
     
    225249    axis[1] = 0;                        /* Y-axis */
    226250
    227     extents_ = 1;
    228     xNum_ = yNum_ = nValues_ = 0;
    229     xMin_ = yMin_ = xMax_ = yMax_ = 0.0f;
    230     if (xUnits_ != NULL) {
    231         free(xUnits_);
    232     }
    233     if (yUnits_ != NULL) {
    234         free(yUnits_);
    235     }
    236     if (vUnits_ != NULL) {
    237         free(vUnits_);
    238     }
    239     xUnits_ = yUnits_ = vUnits_ = NULL;
    240     if (values_ != NULL) {
    241         delete [] values_;
    242     }
    243     values_ = NULL;
     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;
    244267
    245268    int i;
     
    251274        c = string[0];
    252275        if ((c == 'x') && (strcmp(string, "xmin") == 0)) {
    253             if (GetFloatFromObj(interp, objv[i+1], &xMin_) != TCL_OK) {
     276            if (GetFloatFromObj(interp, objv[i+1], &_xMin) != TCL_OK) {
    254277                return TCL_ERROR;
    255278            }
    256279        } else if ((c == 'x') && (strcmp(string, "xmax") == 0)) {
    257             if (GetFloatFromObj(interp, objv[i+1], &xMax_) != TCL_OK) {
     280            if (GetFloatFromObj(interp, objv[i+1], &_xMax) != TCL_OK) {
    258281                return TCL_ERROR;
    259282            }
     
    268291                return TCL_ERROR;
    269292            }
    270             xNum_ = n;
     293            _xNum = n;
    271294        } else if ((c == 'x') && (strcmp(string, "xunits") == 0)) {
    272             xUnits_ = strdup(Tcl_GetString(objv[i+1]));
     295            _xUnits = strdup(Tcl_GetString(objv[i+1]));
    273296        } else if ((c == 'y') && (strcmp(string, "ymin") == 0)) {
    274             if (GetFloatFromObj(interp, objv[i+1], &yMin_) != TCL_OK) {
     297            if (GetFloatFromObj(interp, objv[i+1], &_yMin) != TCL_OK) {
    275298                return TCL_ERROR;
    276299            }
    277300        } else if ((c == 'y') && (strcmp(string, "ymax") == 0)) {
    278             if (GetFloatFromObj(interp, objv[i+1], &yMax_) != TCL_OK) {
     301            if (GetFloatFromObj(interp, objv[i+1], &_yMax) != TCL_OK) {
    279302                return TCL_ERROR;
    280303            }
     
    289312                return TCL_ERROR;
    290313            }
    291             yNum_ = n;
     314            _yNum = n;
    292315        } else if ((c == 'y') && (strcmp(string, "yunits") == 0)) {
    293             yUnits_ = strdup(Tcl_GetString(objv[i+1]));
     316            _yUnits = strdup(Tcl_GetString(objv[i+1]));
    294317        } else if ((c == 'v') && (strcmp(string, "values") == 0)) {
    295318            Tcl_Obj **vobj;
     
    304327                return TCL_ERROR;
    305328            }
    306             nValues_ = n;
    307             values_ = new float[nValues_];
     329            _nValues = n;
     330            _values = new float[_nValues];
    308331            size_t j;
    309             for (j = 0; j < nValues_; j++) {
    310                 if (GetFloatFromObj(interp, vobj[j], values_ + j)!=TCL_OK) {
     332            for (j = 0; j < _nValues; j++) {
     333                if (GetFloatFromObj(interp, vobj[j], _values + j)!=TCL_OK) {
    311334                    return TCL_ERROR;
    312335                }
    313336            }
    314337        } else if ((c == 'u') && (strcmp(string, "units") == 0)) {
    315             vUnits_ = strdup(Tcl_GetString(objv[i+1]));
     338            _vUnits = strdup(Tcl_GetString(objv[i+1]));
    316339        } else if ((c == 'e') && (strcmp(string, "extents") == 0)) {
    317340            int n;
     
    325348                return TCL_ERROR;
    326349            }
    327             extents_ = n;
     350            _nComponents = n;
    328351        } else if ((c == 'a') && (strcmp(string, "axisorder") == 0)) {
    329352            Tcl_Obj **order;
     
    350373        }
    351374    }
    352     if (values_ == NULL) {
     375    if (_values == NULL) {
    353376        Tcl_AppendResult(interp, "missing \"values\" key", (char *)NULL);
    354377        return TCL_ERROR;
    355378    }
    356     if (nValues_ != (xNum_ * yNum_ * extents_)) {
     379    if (_nValues != (_xNum * _yNum * _nComponents)) {
    357380        Tcl_AppendResult(interp,
    358381                "wrong number of values: must be xnum*ynum*extents",
     
    361384    }
    362385   
    363     fprintf(stderr, "generating %dx%dx%d = %d points\n",
    364             xNum_, yNum_, extents_, xNum_ * yNum_ * extents_);
    365 
    366386#ifndef notdef
    367387    if ((axis[0] != 1) || (axis[1] != 0)) {
     
    371391        float *data, *dp;
    372392
    373         dp = data = new float[nValues_];
    374         for (y = 0; y < yNum_; y++) {
     393        dp = data = new float[_nValues];
     394        for (y = 0; y < _yNum; y++) {
    375395            size_t x;
    376396
    377             for (x = 0; x < xNum_; x++) {
     397            for (x = 0; x < _xNum; x++) {
    378398                size_t i, v;
    379399                   
    380400                /* Compute the index from the data's described ordering. */
    381                 i = (y + (yNum_ * x)) * extents_;
    382                 for(v = 0; v < extents_; v++) {
    383                     dp[v] = values_[i+v];
     401                i = (y + (_yNum * x)) * _nComponents;
     402                for(v = 0; v < _nComponents; v++) {
     403                    dp[v] = _values[i+v];
    384404                }
    385                 dp += extents_;
    386             }
    387         }
    388         delete [] values_;
    389         values_ = data;
     405                dp += _nComponents;
     406            }
     407        }
     408        delete [] _values;
     409        _values = data;
    390410    }
    391411#endif
    392     initialized_ = true;
     412    _initialized = true;
    393413    return TCL_OK;
    394414}
     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                    if (nindex >= (npts*3)) {
     523                        fprintf(stderr, "nindex=%d, npts=%d, z=%d, y=%d, x=%d\n",
     524                                nindex, npts, iz, iy, ix);
     525                    }
     526                    _values[nindex] = vx;
     527                    _values[nindex+1] = vy;
     528                    _values[nindex+2] = vz;
     529                    _nValues++;
     530                }
     531            }
     532        }
     533    }
     534    // make sure that we read all of the expected points
     535    if (_nValues != npts) {
     536        result.addError("inconsistent data: expected %d points"
     537                        " but found %d points", npts, _nValues);
     538        delete []  _values;
     539        _values = NULL;
     540        return false;
     541    }
     542    _nValues *= _nComponents;
     543    _initialized = true;
     544    return true;
     545}
     546
     547
     548bool
     549Rappture::Unirect3d::Resample(Rappture::Outcome &result, int nSamples)
     550{
     551    Rappture::Mesh1D xgrid(_xMin, _xMax, _xNum);
     552    Rappture::Mesh1D ygrid(_yMin, _yMax, _yNum);
     553    Rappture::Mesh1D zgrid(_zMin, _zMax, _zNum);
     554    Rappture::FieldRect3D xfield(xgrid, ygrid, zgrid);
     555    Rappture::FieldRect3D yfield(xgrid, ygrid, zgrid);
     556    Rappture::FieldRect3D zfield(xgrid, ygrid, zgrid);
     557
     558    size_t i, j;
     559    for (i = 0, j = 0; i < _nValues; i += _nComponents, j++) {
     560        double vx, vy, vz;
     561
     562        vx = _values[i];
     563        vy = _values[i+1];
     564        vz = _values[i+2];
     565       
     566        xfield.define(j, vx);
     567        yfield.define(j, vy);
     568        zfield.define(j, vz);
     569    }
     570    // Figure out a good mesh spacing
     571    double dx, dy, dz;
     572    dx = xfield.rangeMax(Rappture::xaxis) - xfield.rangeMin(Rappture::xaxis);
     573    dy = xfield.rangeMax(Rappture::yaxis) - xfield.rangeMin(Rappture::yaxis);
     574    dz = xfield.rangeMax(Rappture::zaxis) - xfield.rangeMin(Rappture::zaxis);
     575
     576    double dmin;
     577    dmin = pow((dx*dy*dz)/(nSamples*nSamples*nSamples), 0.333);
     578   
     579    printf("dx:%lf dy:%lf dz:%lf dmin:%lf\n", dx, dy, dz, dmin);
     580
     581    /* Recompute new number of points for each axis. */
     582    _xNum = (size_t)ceil(dx/dmin);
     583    _yNum = (size_t)ceil(dy/dmin);
     584    _zNum = (size_t)ceil(dz/dmin);
     585   
     586#ifndef NV40
     587    // must be an even power of 2 for older cards
     588    _xNum = (int)pow(2.0, ceil(log10((double)_xNum)/log10(2.0)));
     589    _yNum = (int)pow(2.0, ceil(log10((double)_yNum)/log10(2.0)));
     590    _zNum = (int)pow(2.0, ceil(log10((double)_zNum)/log10(2.0)));
     591#endif
     592
     593    size_t n = _nComponents * _xNum * _yNum * _zNum;
     594    float *data = new float[n];
     595    memset(data, 0, sizeof(float) * n);
     596   
     597    // Generate the uniformly sampled rectangle that we need for a volume
     598    float *destPtr = data;
     599    for (size_t i = 0; i < _zNum; i++) {
     600        double z;
     601
     602        z = _zMin + (i * dmin);
     603        for (size_t j = 0; j < _yNum; j++) {
     604            double y;
     605               
     606            y = _yMin + (j * dmin);
     607            for (size_t k = 0; k < _xNum; k++) {
     608                double x;
     609
     610                x = _xMin + (k * dmin);
     611                destPtr[0] = xfield.value(x, y, z);
     612                destPtr[1] = yfield.value(x, y, z);
     613                destPtr[2] = zfield.value(x, y, z);
     614            }
     615        }
     616    }
     617    delete [] _values;
     618    _values = data;
     619    _nValues = _xNum * _yNum * _zNum * _nComponents;
     620    return true;
     621}
     622
Note: See TracChangeset for help on using the changeset viewer.