Changeset 956 for trunk


Ignore:
Timestamp:
Mar 21, 2008, 11:12:26 AM (17 years ago)
Author:
dkearney
Message:

updated the process of querying data from dx files but using native dx library calls instead of calculating grid positions and data points myself.
the Rappture::DX::interpolate() function does not quite work as intended, but if you do not change the axis lengths you can get the original data values back.
created a new function !computeSimpleGradient, located in dxReaderCommon.cpp, as part of the effort to simplify !dxReader.cpp.
removed old code from !dxReader2.cpp

Location:
trunk/vizservers/nanovis
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/vizservers/nanovis/Command.cpp

    r934 r956  
    9393
    9494extern Rappture::Outcome load_volume_stream(int index, std::iostream& fin);
    95 extern Rappture::Outcome load_volume_stream_odx(int index, const char *buf, 
    96                                                 int nBytes);
     95extern Rappture::Outcome load_volume_stream_odx(int index, const char *buf,
     96                        int nBytes);
    9797extern Rappture::Outcome load_volume_stream2(int index, std::iostream& fin);
    98 extern void load_volume(int index, int width, int height, int depth, 
    99                         int n_component, float* data, double vmin, double vmax,
    100                         double nzero_min);
     98extern void load_volume(int index, int width, int height, int depth,
     99            int n_component, float* data, double vmin, double vmax,
     100            double nzero_min);
    101101extern void load_vector_stream(int index, std::istream& fin);
    102102
    … …  
    12641264            NanoVis::n_volumes++;
    12651265        }
    1266        
     1266
    12671267        if (NanoVis::volume[n] != NULL) {
    12681268            delete NanoVis::volume[n];
    12691269            NanoVis::volume[n] = NULL;
    12701270        }
    1271        
     1271
    12721272        float dx0 = -0.5;
    12731273        float dy0 = -0.5*vol->height/vol->width;
    12741274        float dz0 = -0.5*vol->depth/vol->width;
    12751275        vol->move(Vector3(dx0, dy0, dz0));
    1276        
     1276
    12771277        NanoVis::volume[n] = vol;
    12781278#if __TEST_CODE__
    … …  
    12901290    } else if (strcmp(header, "<ODX>") == 0) {
    12911291        Rappture::Outcome err;
    1292        
     1292
    12931293        printf("Loading DX using OpenDX library...\n");
    12941294        fflush(stdout);
    12951295        err = load_volume_stream_odx(n, buf.bytes()+5, buf.size()-5);
    1296         //err = load_volume_stream2(n, fdata);
    12971296        if (err) {
    12981297            Tcl_AppendResult(interp, err.remark().c_str(), (char*)NULL);
    … …  
    13011300    } else {
    13021301        Rappture::Outcome err;
    1303        
     1302
    13041303        printf("OpenDX loading...\n");
    13051304        fflush(stdout);
    13061305        std::stringstream fdata;
    13071306        fdata.write(buf.bytes(),buf.size());
    1308        
     1307
    13091308#if ISO_TEST
    13101309        err = load_volume_stream2(n, fdata);
    … …  
    13171316        }
    13181317    }
    1319    
     1318
    13201319    //
    13211320    // BE CAREFUL: Set the number of slices to something slightly different

  • trunk/vizservers/nanovis/RpDX.cpp

    r932 r956  
    2121using namespace Rappture;
    2222
    23 DX::DX(const char* filename) : 
    24     _dataMin(0),
    25     _dataMax(0),
    26     _nzero_min(0),
     23DX::DX(const char* filename) :
     24    _dataMin(1e21),
     25    _dataMax(1e-21),
     26    _nzero_min(1e21),
    2727    _numAxis(0),
    2828    _axisLen(NULL),
    … …  
    3636    _origin(NULL)
    3737{
    38     Array dxarr;
     38    Array dxpos;
     39    Array dxdata;
     40    Array dxgrid;
    3941    // category and type are probably not needed
    4042    // we keep them around in case they hold useful info for the future
    … …  
    4244    Category category;
    4345    Type type;
    44    
     46
    4547    if (filename == NULL) {
    4648        // error
    4749    }
    48    
     50
    4951    // open the file with libdx
    5052    fprintf(stdout, "Calling DXImportDX(%s)\n", filename);
    … …  
    5658
    5759    // parse out the positions array
    58     dxarr = (Array) DXGetComponentValue((Field) _dxobj, (char *)"positions");
    59     DXGetArrayInfo(dxarr, &_n, &type, &category, &_rank, &_shape);
    60 
    61     _positions = (float*) DXGetArrayData(dxarr);
     60    // FIXME: nanowire will need a different way to parse out the positions array
     61    // since it uses a productarray to store its positions.
     62    // possibly use DXGetProductArray()
     63    dxpos = (Array) DXGetComponentValue((Field) _dxobj, (char *)"positions");
     64    DXGetArrayInfo(dxpos, &_n, &type, &category, &_rank, &_shape);
     65
     66    fprintf(stdout, "_n = %d\n",_n);
     67    fprintf(stdout, "_rank = %d\n",_rank);
     68    fprintf(stdout, "_shape = %d\n",_shape);
     69
     70    float* pos = NULL;
     71    pos = (float*) DXGetArrayData(dxpos);
     72    if (pos == NULL) {
     73        fprintf(stdout, "DXGetArrayData failed to return positions array\n");
     74        fflush(stdout);
     75    }
     76
     77    // first call to get the number of axis needed
     78    dxgrid = (Array) DXGetComponentValue((Field) _dxobj, (char *)"connections");
     79    if (DXQueryGridConnections(dxgrid, &_numAxis, NULL) == NULL) {
     80        // raise error, data is not a regular grid and we cannot handle it
     81        fprintf(stdout,"DX says our grid is not regular, we cannot handle this data\n");
     82        fflush(stdout);
     83    }
     84
     85    _positions = new float[_n*_numAxis];
    6286    if (_positions == NULL) {
    63         fprintf(stdout, "DXGetArrayData failed to return positions array\n");
    64         fflush(stdout);
    65     }
    66 
    67     _numAxis = _rank*_shape;
     87        // malloc failed, raise error
     88        fprintf(stdout, "malloc of _positions array failed");
     89        fflush(stdout);
     90    }
     91    memcpy(_positions,pos,sizeof(float)*_n*_numAxis);
     92
    6893    _axisLen = new int[_numAxis];
    6994    if (_axisLen == NULL) {
    … …  
    7297        fflush(stdout);
    7398    }
    74 
    75     _delta = new float[_numAxis];
     99    memset(_axisLen, 0, _numAxis);
     100
     101    _delta = new float[_numAxis*_numAxis];
    76102    if (_delta == NULL) {
    77103        // malloc failed, raise error
    … …  
    79105        fflush(stdout);
    80106    }
     107    memset(_delta, 0, _numAxis*_numAxis);
     108
     109    _origin = new float[_numAxis];
     110    if (_origin == NULL) {
     111        // malloc failed, raise error
     112        fprintf(stdout, "malloc of _origin array failed");
     113        fflush(stdout);
     114    }
     115    memset(_origin, 0, _numAxis);
    81116
    82117    _max = new float[_numAxis];
    … …  
    86121        fflush(stdout);
    87122    }
    88 
    89     __findPosDeltaMax();
    90 
    91     for (int lcv = 0; lcv < _numAxis; lcv++) {
    92         if (_delta[lcv] == 0) {
    93             // div by 0, raise error
    94             fprintf(stdout, "delta is 0, can't divide by zero\n");
    95             fflush(stdout);
    96         }
    97         // FIXME: find a way to grab the number of points per axis in dx
    98         // this is a horrible way to find the number of points
    99         // per axis and only works when each axis has equal number of pts
    100         _axisLen[lcv] = (int)ceil(pow((double)_n,(double)(1.0/_numAxis)));
    101 
    102     }
    103 
    104     fprintf(stdout, "_max = [%lg,%lg,%lg]\n",_max[0],_max[1],_max[2]);
    105     fprintf(stdout, "_delta = [%lg,%lg,%lg]\n",_delta[0],_delta[1],_delta[2]);
     123    memset(_max, 0, _numAxis);
     124
     125    __findPosMax();
     126
     127    // parse out the gridconnections (length of each axis) array
     128    // dxgrid = (Array) DXQueryGridConnections(dxpos, &_numAxis, _axisLen);
     129    DXQueryGridPositions(dxpos, NULL, _axisLen, _origin, _delta);
     130
     131    fprintf(stdout, "_max = [%g,%g,%g]\n",_max[0],_max[1],_max[2]);
     132    fprintf(stdout, "_delta = [%g,%g,%g]\n",_delta[0],_delta[1],_delta[2]);
     133    fprintf(stdout, "         [%g,%g,%g]\n",_delta[3],_delta[4],_delta[5]);
     134    fprintf(stdout, "         [%g,%g,%g]\n",_delta[6],_delta[7],_delta[8]);
     135    fprintf(stdout, "_origin = [%g,%g,%g]\n",_origin[0],_origin[1],_origin[2]);
    106136    fprintf(stdout, "_axisLen = [%i,%i,%i]\n",_axisLen[0],_axisLen[1],_axisLen[2]);
    107137    fflush(stdout);
    108138
     139    // grab the data array from the dx object and store it in _data
     140    float *data = NULL;
     141    dxdata = (Array) DXGetComponentValue((Field) _dxobj, (char *)"data");
     142    data = (float*) DXGetArrayData(dxdata);
    109143    _data = new float[_n];
    110144    if (_data == NULL) {
    … …  
    113147        fflush(stdout);
    114148    }
    115 
    116     __getInterpData2();
     149    memcpy(_data,data,sizeof(float)*_n);
     150
     151    // print debug info
     152    for (int lcv = 0, pt = 0; lcv < _n; lcv+=3,pt+=9) {
     153        fprintf(stdout,
     154            "(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n",
     155            _positions[pt],_positions[pt+1],_positions[pt+2], _data[lcv],
     156            _positions[pt+3],_positions[pt+4],_positions[pt+5],_data[lcv+1],
     157            _positions[pt+6],_positions[pt+7],_positions[pt+8],_data[lcv+2]);
     158        fflush(stdout);
     159    }
     160
     161    __collectDataStats();
    117162
    118163}
    … …  
    122167    delete[] _axisLen;
    123168    delete[] _delta;
     169    delete[] _origin;
    124170    delete[] _max;
    125171    delete[] _data;
    … …  
    128174
    129175void
    130 DX::__findPosDeltaMax()
     176DX::__findPosMax()
    131177{
    132178    int lcv = 0;
    … …  
    134180    // initialize the max array and delta array
    135181    // max holds the maximum value found for each index
    136     // delta holds the difference between each entry's value
    137182    for (lcv = 0; lcv < _numAxis; lcv++) {
    138183        _max[lcv] = _positions[lcv];
    139         _delta[lcv] = _positions[lcv];
    140184    }
    141185
    142186    for (lcv=lcv; lcv < _n*_numAxis; lcv++) {
    143         if (_positions[lcv] > _max[lcv%_numAxis]) {
    144             _max[lcv%_numAxis] = _positions[lcv];
    145         }
    146         if (_delta[lcv%_numAxis] == _positions[lcv%_numAxis]) {
    147             if (_positions[lcv] != _positions[lcv-_numAxis]) {
    148                 _delta[lcv%_numAxis] = _positions[lcv] - _positions[lcv-_numAxis];
    149             }
     187        int maxIdx = lcv%_numAxis;
     188        if (_positions[lcv] > _max[maxIdx]) {
     189            _max[maxIdx] = _positions[lcv];
    150190        }
    151191    }
    … …  
    153193
    154194void
    155 DX::__getInterpData2()
     195DX::__collectDataStats()
     196{
     197    _dataMin = 1e21;
     198    _dataMax = 1e-21;
     199    _nzero_min = 1e21;
     200
     201    // populate _dataMin, _dataMax, _nzero_min
     202    for (int lcv = 0; lcv < _n; lcv=lcv+3) {
     203        if (_data[lcv] < _dataMin) {
     204            _dataMin = _data[lcv];
     205        }
     206        if (_data[lcv] > _dataMax) {
     207            _dataMax = _data[lcv];
     208        }
     209        if ((_data[lcv] != 0) && (_data[lcv] < _nzero_min)) {
     210            _nzero_min = _data[lcv];
     211        }
     212    }
     213}
     214
     215/*
     216 * getInterpPos()
     217 *
     218 * creates a new grid of positions which can be used for interpolation.
     219 * we create the positions array using the function DXMakeGridPositionsV
     220 * which creates an n-dimensional grid of regularly spaced positions.
     221 * This function overwrites the original positions array
     222 */
     223void
     224DX::__getInterpPos()
     225{
     226    Array dxpos;
     227    float* pos = NULL;
     228
     229    // gather the positions we want to interpolate over
     230    dxpos = DXMakeGridPositionsV(_numAxis, _axisLen, _origin, _delta);
     231    DXGetArrayInfo(dxpos, &_n, NULL, NULL, &_rank, &_shape);
     232    pos = (float*) DXGetArrayData(dxpos);
     233    if (pos == NULL) {
     234        fprintf(stdout, "DXGetArrayData failed to return positions array\n");
     235        fflush(stdout);
     236    }
     237
     238    if (_positions != NULL) {
     239        delete[] _positions;
     240    }
     241    _positions = new float[_n*_numAxis];
     242    if (_positions == NULL) {
     243        // malloc failed, raise error
     244        fprintf(stdout, "malloc of _axisLen array failed");
     245        fflush(stdout);
     246    }
     247    memcpy(_positions,pos,sizeof(float)*_n*_numAxis);
     248
     249    pos = NULL;
     250    DXDelete((object*)dxpos);
     251}
     252
     253/*
     254 * getInterpData()
     255 *
     256 * this function interpolates over a positions array to produce data for each
     257 * point in the positions array. we use the position data stored in _positions
     258 * array.
     259 */
     260void
     261DX::__getInterpData()
    156262{
    157263    int pts = _n;
    158264    int interppts = pts;
    159265    Interpolator interpolator;
     266
     267    _data = new float[_n];
     268    if (_data == NULL) {
     269        // malloc failed, raise error
     270        fprintf(stdout, "malloc of _data array failed");
     271        fflush(stdout);
     272    }
     273    memset(_data,0,_n);
    160274
    161275    // build the interpolator and interpolate
    … …  
    183297    }
    184298
    185     for (int lcv = 0; lcv < pts; lcv=lcv+3) {
    186         if (_data[lcv] < _dataMin) {
    187             _dataMin = _data[lcv];
    188             if (_dataMin != 0) {
    189                 _nzero_min = _dataMin;
    190             }
    191         }
    192         if (_data[lcv] > _dataMax) {
    193             _dataMax = _data[lcv];
    194         }
    195     }
    196 }
    197 
     299    __collectDataStats();
     300}
     301
     302/*
     303 * interpolate()
     304 *
     305 * generate a positions array with optional new axis length and
     306 * interpolate to get the new data values at each point in
     307 * the positions array. this function currently only works if you
     308 * do not change the axis length (i.e. newAxisLen == NULL).
     309 */
     310DX&
     311DX::interpolate(int* newAxisLen)
     312{
     313    fprintf(stdout, "----begin interpolation----\n");
     314    fflush(stdout);
     315    if (newAxisLen != NULL) {
     316        for (int i = 0; i < _numAxis; i++) {
     317            _axisLen[i] = newAxisLen[i];
     318        }
     319    }
     320    __getInterpPos();
     321    __getInterpData();
     322    fprintf(stdout, "----end interpolation----\n");
     323    fflush(stdout);
     324    return *this;
     325}
    198326
    199327int
    … …  
    268396    return _nzero_min;
    269397}
    270 
    271 
    272 

  • trunk/vizservers/nanovis/RpDX.h

    r931 r956  
    4646    virtual float nzero_min() const;
    4747
     48    virtual DX& interpolate(int* newAxisLen);
    4849    virtual int n() const;
    4950    virtual int rank() const;
    … …  
    7273    Object _dxobj;
    7374
    74     void __findPosDeltaMax();
    75     void __getInterpData2();
     75    void __findPosMax();
     76    void __collectDataStats();
     77    void __getInterpPos();
     78    void __getInterpData();
    7679
    7780

  • trunk/vizservers/nanovis/dxReader.cpp

    r932 r956  
    170170                    double vm;
    171171                    vm = sqrt(vx*vx + vy*vy + vz*vz);
    172                     if (vm < vmin) { 
    173                         vmin = vm; 
    174                     } else if (vm > vmax) { 
    175                         vmax = vm; 
     172                    if (vm < vmin) {
     173                        vmin = vm;
     174                    } else if (vm > vmax) {
     175                        vmax = vm;
    176176                    }
    177177                    if ((vm != 0.0f) && (vm < nzero_min)) {
    … …  
    195195                    nzero_min);
    196196
    197         volPtr->xAxis.SetRange(x0, x0 + (nx * ddx));
     197        volPtr->xAxis.SetRange(x0, x0 + (nx * ddx));
    198198        volPtr->yAxis.SetRange(y0, y0 + (ny * ddy));
    199199        volPtr->zAxis.SetRange(z0, z0 + (nz * ddz));
    200200        volPtr->wAxis.SetRange(y0, y0 + (ny * ddy));
    201         volPtr->update_pending = true;
     201        volPtr->update_pending = true;
    202202        delete [] data;
    203203    } else {
    … …  
    357357
    358358                    if (dval[p] < vmin) {
    359             vmin = dval[p];
    360             } else if (dval[p] > vmax) {
    361             vmax = dval[p];
    362             }
     359                        vmin = dval[p];
     360                    } else if (dval[p] > vmax) {
     361                        vmax = dval[p];
     362                    }
    363363                    if (dval[p] != 0.0f && dval[p] < nzero_min) {
    364364                        nzero_min = dval[p];
    … …  
    391391            fflush(stdout);
    392392            if (dv == 0.0) {
    393         dv = 1.0;
    394         }
     393                dv = 1.0;
     394            }
     395
    395396            for (int i = 0; i < count; ++i) {
    396         v = data[ngen];
    397         // scale all values [0-1], -1 => out of bounds
    398         v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
    399         data[ngen] = v;
    400         ngen += 4;
    401         }
    402             // Compute the gradient of this data.  BE CAREFUL: center
    403             // calculation on each node to avoid skew in either direction.
    404             ngen = 0;
    405             for (int iz=0; iz < nz; iz++) {
    406                 for (int iy=0; iy < ny; iy++) {
    407                     for (int ix=0; ix < nx; ix++) {
    408                         // gradient in x-direction
    409                         double valm1 = (ix == 0) ? 0.0 : data[ngen - 4];
    410                         double valp1 = (ix == nx-1) ? 0.0 : data[ngen + 4];
    411                         if (valm1 < 0 || valp1 < 0) {
    412                             data[ngen+1] = 0.0;
    413                         } else {
    414                             data[ngen+1] = valp1-valm1; // assume dx=1
    415                             //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
    416                         }
    417 
    418                         // gradient in y-direction
    419                         valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
    420                         valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
    421                         if (valm1 < 0 || valp1 < 0) {
    422                             data[ngen+2] = 0.0;
    423                         } else {
    424                             data[ngen+2] = valp1-valm1; // assume dx=1
    425                             //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1
    426                         }
    427 
    428                         // gradient in z-direction
    429                         valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
    430                         valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
    431                         if (valm1 < 0 || valp1 < 0) {
    432                             data[ngen+3] = 0.0;
    433                         } else {
    434                             data[ngen+3] = valp1-valm1; // assume dx=1
    435                             //data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
    436                         }
    437 
    438                         ngen += 4;
    439                     }
    440                 }
    441             }
     397                v = data[ngen];
     398                // scale all values [0-1], -1 => out of bounds
     399                v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
     400                data[ngen] = v;
     401                ngen += 4;
     402            }
     403
     404            computeSimpleGradient(data, nx, ny, nz);
    442405
    443406            dx = nx;
    444407            dy = ny;
    445408            dz = nz;
     409
    446410            Volume *volPtr;
    447411            volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
    448                                           vmin, vmax, nzero_min);
    449             fprintf(stderr, "x nx=%d ddx=%g min=%g max=%g\n", nx, ddx,
    450                     x0, x0 + (nx * ddx));
    451             fflush(stderr);
     412                      vmin, vmax, nzero_min);
     413            fprintf(stderr, "x nx=%d ddx=%g min=%g max=%g\n", nx, ddx,
     414                    x0, x0 + (nx * ddx));
     415            fflush(stderr);
    452416#ifdef notdef
    453             volPtr->xAxis.SetRange(x0, x0 + (nx * ddx));
    454             volPtr->yAxis.SetRange(y0, y0 + (ny * ddy));
    455             volPtr->zAxis.SetRange(z0, z0 + (nz * ddz));
     417            volPtr->xAxis.SetRange(x0, x0 + (nx * ddx));
     418            volPtr->yAxis.SetRange(y0, y0 + (ny * ddy));
     419            volPtr->zAxis.SetRange(z0, z0 + (nz * ddz));
    456420#endif
    457             volPtr->xAxis.SetRange(-20.0, 108.0);
    458             volPtr->yAxis.SetRange(0.001, 102.2);
    459             volPtr->zAxis.SetRange(-21, 19);
    460             volPtr->update_pending = true;
     421            volPtr->xAxis.SetRange(-20.0, 108.0);
     422            volPtr->yAxis.SetRange(0.001, 102.2);
     423            volPtr->zAxis.SetRange(-21, 19);
     424            volPtr->update_pending = true;
    461425            delete [] data;
    462426
    … …  
    545509            }
    546510
     511            // FIXME: This next section of code should be replaced by a
     512            // call to the computeSimpleGradient() function. There is a slight
     513            // difference in the code below and the aforementioned function
     514            // in that the commented out lines in the else statements are
     515            // different.
     516            //
    547517            // Compute the gradient of this data.  BE CAREFUL: center
    548518            // calculation on each node to avoid skew in either direction.
    … …  
    588558            Volume *volPtr;
    589559            volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
    590                 field.valueMin(), field.valueMax(), nzero_min);
    591             fprintf(stderr, "x min=%g max=%g\n",
    592                     field.rangeMin(Rappture::xaxis),
    593                     field.rangeMax(Rappture::xaxis));
    594             fflush(stderr);
    595             volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
    596                                field.rangeMax(Rappture::xaxis));
    597             volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
    598                                field.rangeMax(Rappture::yaxis));
    599             volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
    600                                field.rangeMax(Rappture::zaxis));
    601             volPtr->update_pending = true;
     560                        field.valueMin(), field.valueMax(), nzero_min);
     561            fprintf(stderr, "x min=%g max=%g\n",
     562                    field.rangeMin(Rappture::xaxis),
     563                    field.rangeMax(Rappture::xaxis));
     564            fflush(stderr);
     565            volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
     566                   field.rangeMax(Rappture::xaxis));
     567            volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
     568                   field.rangeMax(Rappture::yaxis));
     569            volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
     570                   field.rangeMax(Rappture::zaxis));
     571            volPtr->update_pending = true;
    602572            delete [] data;
    603573        }
    … …  
    750720                    int nindex = iz*nx*ny + iy*nx + ix;
    751721                    field.define(nindex, dval[p]);
    752                     fprintf(stderr,"nindex = %i\tdval[%i] = %lg\n", nindex, p, 
    753                             dval[p]);
     722                    fprintf(stderr,"nindex = %i\tdval[%i] = %lg\n", nindex, p,
     723                        dval[p]);
    754724                    fflush(stderr);
    755725                    nread++;
    … …  
    821791                fflush(stderr);
    822792            }
    823             // Compute the gradient of this data.  BE CAREFUL: center
    824             /*
    825               float *data = new float[4*nx*ny*nz];
    826 
    827               double vmin = field.valueMin();
    828               double dv = field.valueMax() - field.valueMin();
    829               if (dv == 0.0) { dv = 1.0; }
    830 
    831               // generate the uniformly sampled data that we need for a volume
    832               int ngen = 0;
    833               double nzero_min = 0.0;
    834               for (int iz=0; iz < nz; iz++) {
    835               double zval = z0 + iz*dmin;
    836               for (int iy=0; iy < ny; iy++) {
    837               double yval = y0 + iy*dmin;
    838               for (int ix=0; ix < nx; ix++) {
    839               double xval = x0 + ix*dmin;
    840               double v = field.value(xval,yval,zval);
    841 
    842               if (v != 0.0f && v < nzero_min)
    843               {
    844               nzero_min = v;
    845               }
    846 
    847               // scale all values [0-1], -1 => out of bounds
    848               v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
    849 
    850               data[ngen] = v;
    851               ngen += 4;
    852               }
    853               }
    854               }
    855               // Compute the gradient of this data.  BE CAREFUL: center
    856               // calculation on each node to avoid skew in either direction.
    857               ngen = 0;
    858               for (int iz=0; iz < nz; iz++) {
    859               for (int iy=0; iy < ny; iy++) {
    860               for (int ix=0; ix < nx; ix++) {
    861               // gradient in x-direction
    862               double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
    863               double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
    864               if (valm1 < 0 || valp1 < 0) {
    865               data[ngen+1] = 0.0;
    866               } else {
    867               data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
    868               }
    869 
    870               // gradient in y-direction
    871               valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
    872               valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
    873               if (valm1 < 0 || valp1 < 0) {
    874               data[ngen+2] = 0.0;
    875               } else {
    876               //data[ngen+2] = valp1-valm1; // assume dy=1
    877               data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
    878               }
    879 
    880               // gradient in z-direction
    881               valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
    882               valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
    883               if (valm1 < 0 || valp1 < 0) {
    884               data[ngen+3] = 0.0;
    885               } else {
    886               //data[ngen+3] = valp1-valm1; // assume dz=1
    887               data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1
    888               }
    889 
    890               ngen += 4;
    891               }
    892               }
    893               }
    894             */
     793
     794            fprintf(stdout,"End Data Stats index = %i\n",index);
     795            fprintf(stdout,"nx = %i ny = %i nz = %i\n",nx,ny,nz);
     796            fprintf(stdout,"dx = %lg dy = %lg dz = %lg\n",dx,dy,dz);
     797            fprintf(stdout,"dataMin = %lg\tdataMax = %lg\tnzero_min = %lg\n", field.valueMin(),field.valueMax(),nzero_min);
     798            fflush(stdout);
    895799
    896800            Volume *volPtr;
    … …  
    898802            field.valueMin(), field.valueMax(), nzero_min);
    899803            volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
    900                                field.rangeMax(Rappture::xaxis));
     804                   field.rangeMax(Rappture::xaxis));
    901805            volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
    902                                field.rangeMax(Rappture::yaxis));
     806                   field.rangeMax(Rappture::yaxis));
    903807            volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
    904                                field.rangeMax(Rappture::zaxis));
    905             volPtr->update_pending = true;
     808                   field.rangeMax(Rappture::zaxis));
     809            volPtr->update_pending = true;
    906810            // TBD..
    907811            // POINTSET
    … …  
    1045949            field.valueMin(), field.valueMax(), nzero_min);
    1046950            volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
    1047                                field.rangeMax(Rappture::xaxis));
     951                   field.rangeMax(Rappture::xaxis));
    1048952            volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
    1049                                field.rangeMax(Rappture::yaxis));
     953                   field.rangeMax(Rappture::yaxis));
    1050954            volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
    1051                                field.rangeMax(Rappture::zaxis));
    1052             volPtr->update_pending = true;
     955                   field.rangeMax(Rappture::zaxis));
     956            volPtr->update_pending = true;
    1053957            // TBD..
    1054958            // POINTSET

  • trunk/vizservers/nanovis/dxReader2.cpp

    r934 r956  
    1313// rappture headers
    1414#include "RpEncode.h"
    15 #include "RpField1D.h"
    16 #include "RpFieldRect3D.h"
    17 #include "RpFieldPrism3D.h"
    18 
    19 
    20 // FIXME: move newCoolInterplation to the Rappture::DX object
    21 // this definition is kept only so we can test newCoolInterpolation
    22 typedef struct {
    23     int nx;
    24     int ny;
    25     int nz;
    26     double dataMin;
    27     double dataMax;
    28     double nzero_min;
    29     float* data;
    30 } dataInfo;
    31 
    32 // This is legacy code that can be removed during code cleanup
    33 int getInterpPos(  float* numPts, float *origin,
    34                     float* max, int rank, float* arr)
    35 {
    36     float dn[rank];
    37     int i = 0;
    38     float x = 0;
    39     float y = 0;
    40     float z = 0;
    41 
    42     // dn holds the delta you want for interpolation
    43     for (i = 0; i < rank; i++) {
    44         // dn[i] = (max[i] - origin[i]) / (numPts[i] - 1);
    45         dn[i] = (max[i] - origin[i]) / (numPts[i]);
    46     }
    47 
    48     // points are calculated based on the new delta
    49     i = 0;
    50     for (x = origin[0]; x < numPts[0]; x++) {
    51         for (y = origin[1]; y < numPts[1]; y++) {
    52             for (z = origin[2]; z < numPts[2]; z++) {
    53                 arr[i++] = x*dn[0];
    54                 arr[i++] = y*dn[1];
    55                 arr[i++] = z*dn[2];
    56                 fprintf(stderr, "(%f,%f,%f)\n",arr[i-3],arr[i-2],arr[i-1]);
    57             }
    58         }
    59     }
    60 
    61     return 0;
    62 }
    63 
    64 // This is legacy code that can be removed during code cleanup
    65 int
    66 getInterpData(int rank, float* numPts, float* max, float* dxpos, Object* dxobj,
    67               double* dataMin, double* dataMax, float** result)
    68 {
    69     int pts = int(numPts[0]*numPts[1]*numPts[2]);
    70     int interppts = pts;
    71     Interpolator interpolator;
    72 
    73 #ifdef test
    74     int arrSize = interppts*rank;
    75     float interppos[arrSize];
    76 //     commented this out to see what happens when we use the dxpos array
    77 //     which holds the positions dx generated for interpolation.
    78 //
    79 //    generate the positions we want to interpolate on
    80 //    getInterpPos(numPts,dxpos,max,rank,interppos);
    81 //    fprintf(stdout, "after getInterpPos\n");
    82 //    fflush(stdout);
    83 #endif
    84     // build the interpolator and interpolate
    85     interpolator = DXNewInterpolator(*dxobj,INTERP_INIT_IMMEDIATE,-1.0);
    86     fprintf(stdout, "after DXNewInterpolator=%x\n", (unsigned int)interpolator);
    87     fflush(stdout);
    88 //    DXInterpolate(interpolator,&interppts,interppos,*result);
    89     DXInterpolate(interpolator,&interppts,dxpos,*result);
    90     fprintf(stdout,"interppts = %i\n",interppts);
    91     fflush(stdout);
    92 
    93     // print debug info
    94     for (int lcv = 0, pt = 0; lcv < pts; lcv+=3,pt+=9) {
    95         fprintf(stdout,
    96             "(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n",
    97 #ifdef test
    98 //            interppos[pt],interppos[pt+1],interppos[pt+2], (*result)[lcv],
    99 //            interppos[pt+3],interppos[pt+4],interppos[pt+5],(*result)[lcv+1],
    100 //            interppos[pt+6],interppos[pt+7],interppos[pt+8],(*result)[lcv+2]);
    101 #endif
    102             dxpos[pt],dxpos[pt+1],dxpos[pt+2], (*result)[lcv],
    103             dxpos[pt+3],dxpos[pt+4],dxpos[pt+5],(*result)[lcv+1],
    104             dxpos[pt+6],dxpos[pt+7],dxpos[pt+8],(*result)[lcv+2]);
    105         fflush(stdout);
    106     }
    107 
    108     for (int lcv = 0; lcv < pts; lcv=lcv+3) {
    109         if ((*result)[lcv] < *dataMin) {
    110             *dataMin = (*result)[lcv];
    111         }
    112         if ((*result)[lcv] > *dataMax) {
    113             *dataMax = (*result)[lcv];
    114         }
    115     }
    116 
    117     return 0;
    118 }
    119 
    120 /*
    121  * Generalized linear to cartesian transformation function
    122  * This function should work with n dimensional matricies
    123  * with varying axis lengths.
    124  *
    125  * Inputs:
    126  * axisLen - array containing the length of each axis
    127  * idx - the index in a linear array from which x,y,z values are derived
    128  * Outputs:
    129  * xyz - array containing the values of the x, y, and z axis
    130  */
    131 int idx2xyz(const int *axisLen, int idx, int *xyz) {
    132     int mult = 1;
    133     int axis = 0;
    134 
    135     for (axis = 0; axis < 3; axis++) {
    136         xyz[axis] = (idx/mult) % axisLen[axis];
    137         mult = mult*axisLen[axis];
    138     }
    139 
    140     return 0;
    141 }
    142 
    143 /*
    144  * Generalized cartesian to linear transformation function
    145  * This function should work with n dimensional matricies
    146  * with varying axis lengths.
    147  *
    148  * Inputs:
    149  * axisLen - array containing the length of each axis
    150  * xyz - array containing the values of the x, y, and z axis
    151  * Outputs:
    152  * idx - the index that represents xyz in a linear array
    153  */
    154 int xyz2idx(const int *axisLen, int *xyz, int *idx) {
    155     int mult = 1;
    156     int axis = 0;
    157 
    158     *idx = 0;
    159 
    160     for (axis = 0; axis < 3; axis++) {
    161         *idx = *idx + xyz[axis]*mult;
    162         mult = mult*axisLen[axis];
    163     }
    164 
    165     return 0;
    166 }
    167 
    168 Rappture::FieldRect3D* fillRectMesh(const int* numPts, const float* delta, const float* dxpos, const float* data)
    169 {
    170     Rappture::Mesh1D xgrid(dxpos[0], dxpos[0]+numPts[0]*delta[0], (int)(ceil(numPts[0])));
    171     Rappture::Mesh1D ygrid(dxpos[1], dxpos[1]+numPts[1]*delta[1], (int)(ceil(numPts[1])));
    172     Rappture::Mesh1D zgrid(dxpos[2], dxpos[2]+numPts[2]*delta[2], (int)(ceil(numPts[2])));
    173     Rappture::FieldRect3D* field = new Rappture::FieldRect3D(xgrid, ygrid, zgrid);
    174 
    175     int i = 0;
    176     int idx = 0;
    177     int xyz[] = {0,0,0};
    178 
    179     //FIXME: This can be switched to 3 for loops to avoid the if checks
    180     for (i=0; i<numPts[0]*numPts[1]*numPts[2]; i++,xyz[2]++) {
    181         if (xyz[2] >= numPts[2]) {
    182             xyz[2] = 0;
    183             xyz[1]++;
    184             if (xyz[1] >= numPts[1]) {
    185                 xyz[1] = 0;
    186                 xyz[0]++;
    187                 if (xyz[0] >= numPts[0]) {
    188                     xyz[0] = 0;
    189                 }
    190             }
    191         }
    192         xyz2idx(numPts,xyz,&idx);
    193         field->define(idx, data[i]);
    194         fprintf(stdout,"xyz = {%i,%i,%i}\tidx = %i\tdata[%i] = % 8e\n",xyz[0],xyz[1],xyz[2],idx,i,data[i]);
    195         fflush(stdout);
    196     }
    197 
    198     return field;
    199 }
    200 
    201 float* oldCrustyInterpolation(int* nx, int* ny, int* nz, double* dx, double* dy, double* dz, const float* dxpos, Rappture::FieldRect3D* field, double* nzero_min)
    202 {
    203     // figure out a good mesh spacing
    204     int nsample = 30;
    205     *dx = field->rangeMax(Rappture::xaxis) - field->rangeMin(Rappture::xaxis);
    206     *dy = field->rangeMax(Rappture::yaxis) - field->rangeMin(Rappture::yaxis);
    207     *dz = field->rangeMax(Rappture::zaxis) - field->rangeMin(Rappture::zaxis);
    208     double dmin = pow((*dx**dy**dz)/(nsample*nsample*nsample), 0.333);
    209 
    210     *nx = (int)ceil(*dx/dmin);
    211     *ny = (int)ceil(*dy/dmin);
    212     *nz = (int)ceil(*dz/dmin);
    213 
    214 #ifndef NV40
    215     // must be an even power of 2 for older cards
    216     *nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
    217     *ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
    218     *nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
    219 #endif
    220 
    221     float *cdata = new float[*nx**ny**nz];
    222     int ngen = 0;
    223     *nzero_min = 0.0;
    224     for (int iz=0; iz < *nz; iz++) {
    225         double zval = dxpos[2] + iz*dmin;
    226         for (int iy=0; iy < *ny; iy++) {
    227             double yval = dxpos[1] + iy*dmin;
    228             for (int ix=0; ix < *nx; ix++)
    229             {
    230                 double xval = dxpos[0] + ix*dmin;
    231                 double v = field->value(xval,yval,zval);
    232 
    233                 if (v != 0.0f && v < *nzero_min)
    234                 {
    235                     *nzero_min = v;
    236                 }
    237 
    238                 // scale all values [0-1], -1 => out of bounds
    239                 v = (isnan(v)) ? -1.0 : v;
    240 
    241                 cdata[ngen] = v;
    242                 ++ngen;
    243             }
    244         }
    245     }
    246 
    247     float* data = computeGradient(cdata, *nx, *ny, *nz, field->valueMin(), field->valueMax());
    248 
    249     delete[] cdata;
    250     return data;
    251 }
    252 
    253 // FIXME: move newCoolInterplation to the Rappture::DX object
    254 int
    255 newCoolInterpolation(int pts, dataInfo* inData, dataInfo* outData, float* max, float* dxpos, float* delta, float* numPts, double* dx, double* dy, double* dz)
    256 {
    257     // wouldnt dataMin have the non-zero min of the data?
    258     //double nzero_min = 0.0;
    259     outData->nzero_min = inData->dataMin;
    260 
    261     // need to rewrite this starting here!!! dsk
    262     // figure out a good mesh spacing
    263     // dxpos[0,1,2] are the origins for x,y,z axis
    264     *dx = max[0] - dxpos[0];
    265     *dy = max[1] - dxpos[1];
    266     *dz = max[2] - dxpos[2];
    267     double dmin = pow((*dx**dy**dz)/(pts), (double)(1.0/3.0));
    268 
    269     fprintf(stdout, "dx = %f    dy = %f    dz = %f\n",*dx,*dy,*dz);
    270     fprintf(stdout, "dmin = %f\n",dmin);
    271     fprintf(stdout, "max = [%f,%f,%f]\n",max[0],max[1],max[2]);
    272     fprintf(stdout, "delta = [%f,%f,%f]\n",delta[0],delta[1],delta[2]);
    273     fprintf(stdout, "numPts = [%f,%f,%f]\n",numPts[0],numPts[1],numPts[2]);
    274     fflush(stdout);
    275 
    276 
    277     outData->nx = (int)ceil(numPts[0]);
    278     outData->ny = (int)ceil(numPts[1]);
    279     outData->nz = (int)ceil(numPts[2]);
    280 
    281     fprintf(stdout, "nx = %i    ny = %i    nz = %i\n",outData->nx,outData->ny,outData->nz);
    282 
    283     outData->data = new float[outData->nx*outData->ny*outData->nz];
    284     for (int i=0; i<(outData->nx*outData->ny*outData->nz); i+=1) {
    285         // scale all values [0,1], -1 => out of bounds
    286         outData->data[i] = (isnan(inData->data[i])) ? -1.0 :
    287             (inData->data[i] - inData->dataMin)/(inData->dataMax-inData->dataMin);
    288         fprintf(stdout, "outData[%i] = % 8e\tinData[%i] = % 8e\n",i,outData->data[i],i,inData->data[i]);
    289     }
    290 
    291     return 0;
    292 }
    29315
    29416/* Load a 3D volume from a dx-format file the new way
    29517 */
    29618Rappture::Outcome
    297 load_volume_stream_odx(int index, const char *buf, int nBytes) {
     19load_volume_stream_odx(int index, const char *buf, int nBytes)
     20{
    29821    Rappture::Outcome outcome;
    29922
    … …  
    32144    }
    32245
    323     // store our data in a rappture structure to be interpolated
    324     Rappture::FieldRect3D* field = fillRectMesh(dxObj.axisLen(),
    325                                                 dxObj.delta(),
    326                                                 dxObj.positions(),
    327                                                 dxObj.data());
     46    int nx = dxObj.axisLen()[0];
     47    int ny = dxObj.axisLen()[1];
     48    int nz = dxObj.axisLen()[2];
     49    float dx = nx;
     50    float dy = ny;
     51    float dz = nz;
    32852
    329     int nx = 0;
    330     int ny = 0;
    331     int nz = 0;
    332     double dx = 0.0;
    333     double dy = 0.0;
    334     double dz = 0.0;
    335     double nzero_min = 0.0;
    336     float* data = oldCrustyInterpolation(&nx,&ny,&nz,&dx,&dy,&dz,
    337                                         dxObj.positions(),field,&nzero_min);
     53    const float* data1 = dxObj.data();
     54    float *data = new float[nx*ny*nz*4];
     55    memset(data, 0, nx*ny*nz*4);
     56    int iz=0, ix=0, iy=0;
     57    float dv = dxObj.dataMax() - dxObj.dataMin();
     58    float vmin = dxObj.dataMin();
    33859
    339     double dataMin = field->valueMin();
    340     double dataMax = field->valueMax();
     60    for (int i=0; i < nx*ny*nz; i++) {
     61        int nindex = (iz*nx*ny + iy*nx + ix) * 4;
     62        float v = data1[i];
    34163
     64        // scale all values [0-1], -1 => out of bounds
     65        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
    34266
    343     for (int i=0; i<nx*ny*nz; i++) {
    344         fprintf(stdout,"enddata[%i] = %lg\n",i,data[i]);
    345         fflush(stdout);
     67        // place the value in the correct index in cdata
     68        data[nindex] = v;
     69
     70        // calculate next x,y,z coordinate
     71        if (++iz >= nz) {
     72            iz = 0;
     73            if (++iy >= ny) {
     74                iy = 0;
     75                ++ix;
     76            }
     77        }
    34678    }
     79
     80    computeSimpleGradient(data, nx, ny, nz);
    34781
    34882    fprintf(stdout,"End Data Stats index = %i\n",index);
    34983    fprintf(stdout,"nx = %i ny = %i nz = %i\n",nx,ny,nz);
    35084    fprintf(stdout,"dx = %lg dy = %lg dz = %lg\n",dx,dy,dz);
    351     fprintf(stdout,"dataMin = %lg\tdataMax = %lg\tnzero_min = %lg\n", dataMin,dataMax,nzero_min);
     85    fprintf(stdout,"dataMin = %lg\tdataMax = %lg\tnzero_min = %lg\n", dxObj.dataMin(),dxObj.dataMax(),dxObj.nzero_min());
    35286    fflush(stdout);
    353     NanoVis::load_volume(index, nx, ny, nz, 4, data, dataMin, dataMax,
    354                              nzero_min);
    35587
    356     // free the result array
    357     // delete[] data;
    358     // delete[] result;
    359     // delete field;
     88    Volume *volPtr;
     89    volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
     90                                  dxObj.dataMin(), dxObj.dataMax(), dxObj.nzero_min());
     91    /*
     92    volPtr->SetRange(AxisRange::X, x0, x0 + (nx * ddx));
     93    volPtr->SetRange(AxisRange::Y, y0, y0 + (ny * ddy));
     94    volPtr->SetRange(AxisRange::Z, z0, z0 + (nz * ddz));
     95    */
     96
     97    volPtr->update_pending = true;
     98    delete [] data;
    36099
    361100    //

  • trunk/vizservers/nanovis/dxReaderCommon.cpp

    r931 r956  
    44#include "Vector3.h"
    55#include "stdlib.h"
    6 
    7 /*
    8 float* merge(float* scalar, float* gradient, int size)
    9 {
    10     float* data = (float*) malloc(sizeof(float) * 4 * size);
    11 
    12     Vector3* g = (Vector3*) gradient;
    13 
    14     int ngen = 0, sindex = 0;
    15     for (sindex = 0; sindex <size; ++sindex)
    16     {
    17         data[ngen++] = scalar[sindex];
    18         data[ngen++] = g[sindex].x;
    19         data[ngen++] = g[sindex].y;
    20         data[ngen++] = g[sindex].z;
    21     }
    22     return data;
    23 }
    24 
    25 void normalizeScalar(float* fdata, int count, float min, float max)
    26 {
    27     float v = max - min;
    28     if (v != 0.0f)
    29     {
    30         for (int i = 0; i < count; ++i)
    31             fdata[i] = fdata[i] / v;
    32     }
    33 }
    34 
    35 float* computeGradient(float* fdata, int width, int height, int depth, float min, float max)
    36 {
    37     float* gradients = (float *)malloc(width * height * depth * 3 * sizeof(float));
    38     float* tempGradients = (float *)malloc(width * height * depth * 3 * sizeof(float));
    39     int sizes[3] = { width, height, depth };
    40     computeGradients(tempGradients, fdata, sizes, DATRAW_FLOAT);
    41     filterGradients(tempGradients, sizes);
    42     quantizeGradients(tempGradients, gradients, sizes, DATRAW_FLOAT);
    43     normalizeScalar(fdata, width * height * depth, min, max);
    44     float* data = merge(fdata, gradients, width * height * depth);
    45 
    46     return data;
    47 }
    48 */
    496
    507float *
    … …  
    9350}
    9451
     52void
     53computeSimpleGradient(float* data, int nx, int ny, int nz)
     54{
     55    // Compute the gradient of this data.  BE CAREFUL: center
     56    // calculation on each node to avoid skew in either direction.
     57    int ngen = 0;
     58    for (int iz=0; iz < nz; iz++) {
     59        for (int iy=0; iy < ny; iy++) {
     60            for (int ix=0; ix < nx; ix++) {
     61                // gradient in x-direction
     62                double valm1 = (ix == 0) ? 0.0 : data[ngen - 4];
     63                double valp1 = (ix == nx-1) ? 0.0 : data[ngen + 4];
     64                if (valm1 < 0 || valp1 < 0) {
     65                    data[ngen+1] = 0.0;
     66                } else {
     67                    data[ngen+1] = valp1-valm1; // assume dx=1
     68                    //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
     69                }
     70
     71                // gradient in y-direction
     72                valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
     73                valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
     74                if (valm1 < 0 || valp1 < 0) {
     75                    data[ngen+2] = 0.0;
     76                } else {
     77                    data[ngen+2] = valp1-valm1; // assume dx=1
     78                    //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1
     79                }
     80
     81                // gradient in z-direction
     82                valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
     83                valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
     84                if (valm1 < 0 || valp1 < 0) {
     85                    data[ngen+3] = 0.0;
     86                } else {
     87                    data[ngen+3] = valp1-valm1; // assume dx=1
     88                    //data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
     89                }
     90
     91                ngen += 4;
     92            }
     93        }
     94    }
     95}

  • trunk/vizservers/nanovis/dxReaderCommon.h

    r932 r956  
    22#define _DX_READER_COMMON_H
    33
    4 
    54float* merge(float* scalar, float* gradient, int size);
    65void normalizeScalar(float* fdata, int count, float min, float max);
    7 float* computeGradient(float* fdata, int width, int height, int depth,
    8         float min, float max);
     6float* computeGradient(float* fdata, int width, int height, int depth,
     7        float min, float max);
     8void computeSimpleGradient(float* data, int nx, int ny, int nz);
    99
    1010#endif /*_DX_READER_COMMON_H*/
Note: See TracChangeset for help on using the changeset viewer.