source: trunk/packages/vizservers/nanovis/RpDX.cpp @ 2798

Last change on this file since 2798 was 2798, checked in by ldelgass, 8 years ago

Add emacs mode magic line in preparation for indentation cleanup

  • Property svn:eol-style set to native
File size: 9.8 KB
Line 
1/* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */
2/*
3 * ----------------------------------------------------------------------
4 *  Rappture::DX
5 *
6 *  Rappture DX object for file reading and interacting
7 *  with libDX and friends.
8 *
9 * ======================================================================
10 *  AUTHOR:  Derrick S. Kearney, Purdue University
11 *  Copyright (c) 2005-2009  Purdue Research Foundation
12 *
13 *  See the file "license.terms" for information on usage and
14 *  redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES.
15 * ======================================================================
16 */
17#include "nvconf.h"
18#ifdef HAVE_DX_DX_H
19#include "RpDX.h"
20#undef ERROR
21#include "Trace.h"
22#include <math.h>
23#include <stdio.h>
24#include <stdlib.h>
25#include <float.h>
26
27using namespace Rappture;
28
29DX::DX(Outcome &result, const char* filename) :
30    _dataMin(FLT_MAX),
31    _dataMax(-FLT_MAX),
32    _nzero_min(FLT_MAX),
33    _numAxis(0),
34    _axisLen(NULL),
35    _data(NULL),
36    _n(0),
37    _rank(0),
38    _shape(0),
39    _positions(NULL),
40    _delta(NULL),
41    _max(NULL),
42    _origin(NULL)
43{
44    Array dxpos;
45    Array dxdata;
46    Array dxgrid;
47    // category and type are probably not needed
48    // we keep them around in case they hold useful info for the future
49    // we could replace them with NULL inthe DXGetArrayInfo fxn call
50    Category category;
51    Type type;
52
53    if (filename == NULL) {
54        result.addError("filename is NULL");
55        return;
56    }
57    // open the file with libdx
58    TRACE("Calling DXImportDX(%s)\n", filename);
59    DXenable_locks(0);
60    _dxobj = DXImportDX((char*)filename,NULL,NULL,NULL,NULL);
61    if (_dxobj == NULL) {
62        result.addError("can't import DX file \"%s\"", filename);
63        return;
64    }
65
66    // parse out the positions array
67
68    // FIXME: nanowire will need a different way to parse out the positions
69    //        array since it uses a productarray to store its positions.
70    //        Possibly use DXGetProductArray().
71    dxpos = (Array) DXGetComponentValue((Field) _dxobj, (char *)"positions");
72    if (dxpos == NULL) {
73        result.addError("can't get component value of \"positions\"");
74        return;
75    }
76    DXGetArrayInfo(dxpos, &_n, &type, &category, &_rank, &_shape);
77
78    TRACE("_n = %d\n",_n);
79    if (type != TYPE_FLOAT) {
80        result.addError("\"positions\" is not type float (type=%d)\n", type);
81        return;
82    }
83    TRACE("_rank = %d\n",_rank);
84    TRACE("_shape = %d\n",_shape);
85
86    float* pos = NULL;
87    pos = (float*) DXGetArrayData(dxpos);
88    if (pos == NULL) {
89        result.addError("DXGetArrayData failed to return positions array");
90        return;
91    }
92
93    // first call to get the number of axis needed
94    dxgrid = (Array) DXGetComponentValue((Field) _dxobj, (char *)"connections");
95    if (DXQueryGridConnections(dxgrid, &_numAxis, NULL) == NULL) {
96        // raise error, data is not a regular grid and we cannot handle it
97        result.addError("DX says our grid is not regular, we cannot handle this data");
98        return;
99    }
100
101    _positions = new float[_n*_numAxis];
102    if (_positions == NULL) {
103        // malloc failed, raise error
104        result.addError("malloc of _positions array failed");
105        return;
106    }
107    memcpy(_positions,pos,sizeof(float)*_n*_numAxis);
108
109    _axisLen = new int[_numAxis];
110    if (_axisLen == NULL) {
111        // malloc failed, raise error
112        result.addError("malloc of _axisLen array failed");
113        return;
114    }
115    memset(_axisLen, 0, _numAxis);
116
117    _delta = new float[_numAxis*_numAxis];
118    if (_delta == NULL) {
119        result.addError("malloc of _delta array failed");
120        return;
121    }
122    memset(_delta, 0, _numAxis*_numAxis);
123
124    _origin = new float[_numAxis];
125    if (_origin == NULL) {
126        result.addError("malloc of _origin array failed");
127        return;
128    }
129    memset(_origin, 0, _numAxis);
130
131    _max = new float[_numAxis];
132    if (_max == NULL) {
133        result.addError("malloc of _max array failed");
134        return;
135    }
136    memset(_max, 0, _numAxis);
137
138    __findPosMax();
139
140    // parse out the gridconnections (length of each axis) array
141    // dxgrid = (Array) DXQueryGridConnections(dxpos, &_numAxis, _axisLen);
142    DXQueryGridPositions(dxpos, NULL, _axisLen, _origin, _delta);
143
144    TRACE("_max = [%g,%g,%g]\n",_max[0],_max[1],_max[2]);
145    TRACE("_delta = [%g,%g,%g]\n",_delta[0],_delta[1],_delta[2]);
146    TRACE("         [%g,%g,%g]\n",_delta[3],_delta[4],_delta[5]);
147    TRACE("         [%g,%g,%g]\n",_delta[6],_delta[7],_delta[8]);
148    TRACE("_origin = [%g,%g,%g]\n",_origin[0],_origin[1],_origin[2]);
149    TRACE("_axisLen = [%i,%i,%i]\n",_axisLen[0],_axisLen[1],_axisLen[2]);
150
151    // grab the data array from the dx object and store it in _data
152    dxdata = (Array) DXGetComponentValue((Field) _dxobj, (char *)"data");
153    DXGetArrayInfo(dxdata, NULL, &type, NULL, NULL, NULL);
154    _data = new float[_n];
155    if (_data == NULL) {
156        result.addError("malloc of _data array failed");
157        return;
158    }
159
160    switch (type) {
161    case TYPE_FLOAT:
162        float *float_data;
163
164        float_data = (float *)DXGetArrayData(dxdata);
165        memcpy(_data, float_data, sizeof(float)*_n);
166        break;
167
168    case TYPE_DOUBLE:
169        double *double_data;
170        double_data = (double *)DXGetArrayData(dxdata);
171        for (int i = 0; i < _n; i++) {
172            _data[i] = double_data[i];
173        }
174        break;
175
176    default:
177        result.addError("don't know how to handle data of type %d\n", type);
178        return;
179    }
180
181    // print debug info
182    for (int lcv = 0, pt = 0; lcv < _n; lcv +=3, pt+=9) {
183        TRACE("(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n",
184            _positions[pt],_positions[pt+1],_positions[pt+2], _data[lcv],
185            _positions[pt+3],_positions[pt+4],_positions[pt+5],_data[lcv+1],
186            _positions[pt+6],_positions[pt+7],_positions[pt+8],_data[lcv+2]);
187    }
188    __collectDataStats();
189}
190
191DX::~DX()
192{
193    delete[] _axisLen;
194    delete[] _delta;
195    delete[] _origin;
196    delete[] _max;
197    delete[] _data;
198    delete[] _positions;
199}
200
201void
202DX::__findPosMax()
203{
204    int lcv = 0;
205
206    // initialize the max array and delta array
207    // max holds the maximum value found for each index
208    for (lcv = 0; lcv < _numAxis; lcv++) {
209        _max[lcv] = _positions[lcv];
210    }
211
212    for (lcv=lcv; lcv < _n*_numAxis; lcv++) {
213        int maxIdx = lcv%_numAxis;
214        if (_positions[lcv] > _max[maxIdx]) {
215            _max[maxIdx] = _positions[lcv];
216        }
217    }
218}
219
220void
221DX::__collectDataStats()
222{
223    _dataMin = FLT_MAX;
224    _dataMax = -FLT_MAX;
225    _nzero_min = FLT_MAX;
226
227    // populate _dataMin, _dataMax, _nzero_min
228    for (int lcv = 0; lcv < _n; lcv++) {
229        if (_data[lcv] < _dataMin) {
230            _dataMin = _data[lcv];
231        }
232        if (_data[lcv] > _dataMax) {
233            _dataMax = _data[lcv];
234        }
235        if ((_data[lcv] > 0.0f) && (_data[lcv] < _nzero_min)) {
236            _nzero_min = _data[lcv];
237        }
238    }
239    if (_nzero_min == FLT_MAX) {
240        ERROR("could not find a positive minimum value\n");
241    }
242}
243
244/*
245 * getInterpPos()
246 *
247 * creates a new grid of positions which can be used for interpolation.
248 * we create the positions array using the function DXMakeGridPositionsV
249 * which creates an n-dimensional grid of regularly spaced positions.
250 * This function overwrites the original positions array
251 */
252void
253DX::__getInterpPos()
254{
255    Array dxpos;
256    float* pos = NULL;
257
258    // gather the positions we want to interpolate over
259    dxpos = DXMakeGridPositionsV(_numAxis, _axisLen, _origin, _delta);
260    DXGetArrayInfo(dxpos, &_n, NULL, NULL, &_rank, &_shape);
261    pos = (float*) DXGetArrayData(dxpos);
262    if (pos == NULL) {
263        TRACE("DXGetArrayData failed to return positions array\n");
264    }
265
266    if (_positions != NULL) {
267        delete[] _positions;
268    }
269    _positions = new float[_n*_numAxis];
270    if (_positions == NULL) {
271        // malloc failed, raise error
272        TRACE("malloc of _axisLen array failed");
273    }
274    memcpy(_positions,pos,sizeof(float)*_n*_numAxis);
275
276    pos = NULL;
277    DXDelete((object*)dxpos);
278}
279
280/*
281 * getInterpData()
282 *
283 * this function interpolates over a positions array to produce data for each
284 * point in the positions array. we use the position data stored in _positions
285 * array.
286 */
287void
288DX::__getInterpData()
289{
290    int pts = _n;
291    int interppts = pts;
292    Interpolator interpolator;
293
294    _data = new float[_n];
295    if (_data == NULL) {
296        // malloc failed, raise error
297        TRACE("malloc of _data array failed");
298    }
299    memset(_data,0,_n);
300
301    // build the interpolator and interpolate
302    TRACE("creating DXNewInterpolator...\n");
303    interpolator = DXNewInterpolator(_dxobj,INTERP_INIT_IMMEDIATE,-1.0);
304    TRACE("_rank = %i\n",_rank);
305    TRACE("_shape = %i\n",_shape);
306    TRACE("_n = %i\n",_n);
307    TRACE("start interppts = %i\n",interppts);
308    DXInterpolate(interpolator,&interppts,_positions,_data);
309    TRACE("interppts = %i\n",interppts);
310
311    // print debug info
312    for (int lcv = 0, pt = 0; lcv < pts; lcv+=3,pt+=9) {
313        TRACE("(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n",
314            _positions[pt],_positions[pt+1],_positions[pt+2], _data[lcv],
315            _positions[pt+3],_positions[pt+4],_positions[pt+5],_data[lcv+1],
316            _positions[pt+6],_positions[pt+7],_positions[pt+8],_data[lcv+2]);
317    }
318
319    __collectDataStats();
320}
321
322/*
323 * interpolate()
324 *
325 * generate a positions array with optional new axis length and
326 * interpolate to get the new data values at each point in
327 * the positions array. this function currently only works if you
328 * do not change the axis length (i.e. newAxisLen == NULL).
329 */
330DX&
331DX::interpolate(int* newAxisLen)
332{
333    TRACE("----begin interpolation----\n");
334    if (newAxisLen != NULL) {
335        for (int i = 0; i < _numAxis; i++) {
336            _axisLen[i] = newAxisLen[i];
337        }
338    }
339    __getInterpPos();
340    __getInterpData();
341    TRACE("----end interpolation----\n");
342    return *this;
343}
344
345#endif /*HAVE_DX_DX_H*/
Note: See TracBrowser for help on using the repository browser.