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

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

Cleanups, no functional changes

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