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

Last change on this file since 1028 was 1028, checked in by gah, 13 years ago

various cleanups

File size: 10.9 KB
Line 
1/*
2 * ----------------------------------------------------------------------
3 *  Rappture::DX
4 *
5 *  Rappture DX object for file reading and interacting
6 *  with libDX and friends.
7 *
8 * ======================================================================
9 *  AUTHOR:  Derrick S. Kearney, Purdue University
10 *  Copyright (c) 2005-2008  Purdue Research Foundation
11 *
12 *  See the file "license.terms" for information on usage and
13 *  redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES.
14 * ======================================================================
15 */
16#include "RpDX.h"
17#include <math.h>
18#include <stdio.h>
19#include <stdlib.h>
20#include <float.h>
21
22using namespace Rappture;
23
24DX::DX(const char* filename, Outcome *resultPtr) :
25    _dataMin(FLT_MAX),
26    _dataMax(-FLT_MAX),
27    _nzero_min(FLT_MAX),
28    _numAxis(0),
29    _axisLen(NULL),
30    _data(NULL),
31    _n(0),
32    _rank(0),
33    _shape(0),
34    _positions(NULL),
35    _delta(NULL),
36    _max(NULL),
37    _origin(NULL)
38{
39    Array dxpos;
40    Array dxdata;
41    Array dxgrid;
42    // category and type are probably not needed
43    // we keep them around in case they hold useful info for the future
44    // we could replace them with NULL inthe DXGetArrayInfo fxn call
45    Category category;
46    Type type;
47
48    if (filename == NULL) {
49        resultPtr->AddError("filename is NULL");
50        return;
51    }
52    // open the file with libdx
53    fprintf(stdout, "Calling DXImportDX(%s)\n", filename);
54    fflush(stdout);
55    DXenable_locks(0);
56    _dxobj = DXImportDX((char*)filename,NULL,NULL,NULL,NULL);
57    if (_dxobj == NULL) {
58        resultPtr->AddError("can't import DX file \"%s\"", filename);
59        return;
60    }
61
62    fprintf(stdout, "dxobj=%x\n", (unsigned int)_dxobj);
63    fflush(stdout);
64
65    // parse out the positions array
66
67    // FIXME: nanowire will need a different way to parse out the positions
68    //        array since it uses a productarray to store its positions. 
69    //        Possibly use DXGetProductArray().
70    dxpos = (Array) DXGetComponentValue((Field) _dxobj, (char *)"positions");
71    if (dxpos == NULL) {
72        resultPtr->AddError("can't get component value of \"positions\"");
73        return;
74    }
75    DXGetArrayInfo(dxpos, &_n, &type, &category, &_rank, &_shape);
76
77    fprintf(stdout, "_n = %d\n",_n);
78    if (type != TYPE_FLOAT) {
79        resultPtr->AddError("\"positions\" is not type float (type=%d)\n", type);
80        return;
81    }
82    fprintf(stdout, "_rank = %d\n",_rank);
83    fprintf(stdout, "_shape = %d\n",_shape);
84
85    float* pos = NULL;
86    pos = (float*) DXGetArrayData(dxpos);
87    if (pos == NULL) {
88        resultPtr->AddError("DXGetArrayData failed to return positions array");
89        return;
90    }
91
92    // first call to get the number of axis needed
93    dxgrid = (Array) DXGetComponentValue((Field) _dxobj, (char *)"connections");
94    if (DXQueryGridConnections(dxgrid, &_numAxis, NULL) == NULL) {
95        // raise error, data is not a regular grid and we cannot handle it
96        resultPtr->AddError("DX says our grid is not regular, we cannot handle this data");
97        return;
98    }
99
100    _positions = new float[_n*_numAxis];
101    if (_positions == NULL) {
102        // malloc failed, raise error
103        resultPtr->AddError("malloc of _positions array failed");
104        return;
105    }
106    memcpy(_positions,pos,sizeof(float)*_n*_numAxis);
107
108    _axisLen = new int[_numAxis];
109    if (_axisLen == NULL) {
110        // malloc failed, raise error
111        resultPtr->AddError("malloc of _axisLen array failed");
112        return;
113    }
114    memset(_axisLen, 0, _numAxis);
115
116    _delta = new float[_numAxis*_numAxis];
117    if (_delta == NULL) {
118        resultPtr->AddError("malloc of _delta array failed");
119        return;
120    }
121    memset(_delta, 0, _numAxis*_numAxis);
122
123    _origin = new float[_numAxis];
124    if (_origin == NULL) {
125        resultPtr->AddError("malloc of _origin array failed");
126        return;
127    }
128    memset(_origin, 0, _numAxis);
129
130    _max = new float[_numAxis];
131    if (_max == NULL) {
132        resultPtr->AddError("malloc of _max array failed");
133        return;
134    }
135    memset(_max, 0, _numAxis);
136
137    __findPosMax();
138
139    // parse out the gridconnections (length of each axis) array
140    // dxgrid = (Array) DXQueryGridConnections(dxpos, &_numAxis, _axisLen);
141    DXQueryGridPositions(dxpos, NULL, _axisLen, _origin, _delta);
142
143    fprintf(stdout, "_max = [%g,%g,%g]\n",_max[0],_max[1],_max[2]);
144    fprintf(stdout, "_delta = [%g,%g,%g]\n",_delta[0],_delta[1],_delta[2]);
145    fprintf(stdout, "         [%g,%g,%g]\n",_delta[3],_delta[4],_delta[5]);
146    fprintf(stdout, "         [%g,%g,%g]\n",_delta[6],_delta[7],_delta[8]);
147    fprintf(stdout, "_origin = [%g,%g,%g]\n",_origin[0],_origin[1],_origin[2]);
148    fprintf(stdout, "_axisLen = [%i,%i,%i]\n",_axisLen[0],_axisLen[1],_axisLen[2]);
149    fflush(stdout);
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        resultPtr->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    case TYPE_DOUBLE:
168        double *double_data;
169        double_data = (double*) DXGetArrayData(dxdata);
170        for (int i = 0; i < _n; i++) {
171            _data[i] = double_data[i];
172        }
173    default:
174        resultPtr->AddError("don't know how to handle data of type %d\n", type);
175        return;
176    }
177
178    // print debug info
179    for (int lcv = 0, pt = 0; lcv < _n; lcv +=3, pt+=9) {
180        fprintf(stdout,
181            "(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n",
182            _positions[pt],_positions[pt+1],_positions[pt+2], _data[lcv],
183            _positions[pt+3],_positions[pt+4],_positions[pt+5],_data[lcv+1],
184            _positions[pt+6],_positions[pt+7],_positions[pt+8],_data[lcv+2]);
185        fflush(stdout);
186    }
187    __collectDataStats();
188}
189
190DX::~DX()
191{
192    delete[] _axisLen;
193    delete[] _delta;
194    delete[] _origin;
195    delete[] _max;
196    delete[] _data;
197    delete[] _positions;
198}
199
200void
201DX::__findPosMax()
202{
203    int lcv = 0;
204
205    // initialize the max array and delta array
206    // max holds the maximum value found for each index
207    for (lcv = 0; lcv < _numAxis; lcv++) {
208        _max[lcv] = _positions[lcv];
209    }
210
211    for (lcv=lcv; lcv < _n*_numAxis; lcv++) {
212        int maxIdx = lcv%_numAxis;
213        if (_positions[lcv] > _max[maxIdx]) {
214            _max[maxIdx] = _positions[lcv];
215        }
216    }
217}
218
219void
220DX::__collectDataStats()
221{
222    _dataMin = FLT_MAX;
223    _dataMax = -FLT_MAX;
224    _nzero_min = FLT_MAX;
225
226    // populate _dataMin, _dataMax, _nzero_min
227    for (int lcv = 0; lcv < _n; lcv++) {
228        if (_data[lcv] < _dataMin) {
229            _dataMin = _data[lcv];
230        }
231        if (_data[lcv] > _dataMax) {
232            _dataMax = _data[lcv];
233        }
234        if ((_data[lcv] > 0.0f) && (_data[lcv] < _nzero_min)) {
235            _nzero_min = _data[lcv];
236        }
237    }
238    if (_nzero_min == FLT_MAX) {
239        fprintf(stderr, "could not find a positive minimum value\n");
240        fflush(stderr);
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        fprintf(stdout, "DXGetArrayData failed to return positions array\n");
264        fflush(stdout);
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        fprintf(stdout, "malloc of _axisLen array failed");
274        fflush(stdout);
275    }
276    memcpy(_positions,pos,sizeof(float)*_n*_numAxis);
277
278    pos = NULL;
279    DXDelete((object*)dxpos);
280}
281
282/*
283 * getInterpData()
284 *
285 * this function interpolates over a positions array to produce data for each
286 * point in the positions array. we use the position data stored in _positions
287 * array.
288 */
289void
290DX::__getInterpData()
291{
292    int pts = _n;
293    int interppts = pts;
294    Interpolator interpolator;
295
296    _data = new float[_n];
297    if (_data == NULL) {
298        // malloc failed, raise error
299        fprintf(stdout, "malloc of _data array failed");
300        fflush(stdout);
301    }
302    memset(_data,0,_n);
303
304    // build the interpolator and interpolate
305    fprintf(stdout, "creating DXNewInterpolator...\n");
306    fflush(stdout);
307    interpolator = DXNewInterpolator(_dxobj,INTERP_INIT_IMMEDIATE,-1.0);
308    fprintf(stdout, "DXNewInterpolator=%x\n", (unsigned int)interpolator);
309    fprintf(stdout,"_rank = %i\n",_rank);
310    fprintf(stdout,"_shape = %i\n",_shape);
311    fprintf(stdout,"_n = %i\n",_n);
312    fprintf(stdout,"start interppts = %i\n",interppts);
313    fflush(stdout);
314    DXInterpolate(interpolator,&interppts,_positions,_data);
315    fprintf(stdout,"interppts = %i\n",interppts);
316    fflush(stdout);
317
318    // print debug info
319    for (int lcv = 0, pt = 0; lcv < pts; lcv+=3,pt+=9) {
320        fprintf(stdout,
321            "(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n",
322            _positions[pt],_positions[pt+1],_positions[pt+2], _data[lcv],
323            _positions[pt+3],_positions[pt+4],_positions[pt+5],_data[lcv+1],
324            _positions[pt+6],_positions[pt+7],_positions[pt+8],_data[lcv+2]);
325        fflush(stdout);
326    }
327
328    __collectDataStats();
329}
330
331/*
332 * interpolate()
333 *
334 * generate a positions array with optional new axis length and
335 * interpolate to get the new data values at each point in
336 * the positions array. this function currently only works if you
337 * do not change the axis length (i.e. newAxisLen == NULL).
338 */
339DX&
340DX::interpolate(int* newAxisLen)
341{
342    fprintf(stdout, "----begin interpolation----\n");
343    fflush(stdout);
344    if (newAxisLen != NULL) {
345        for (int i = 0; i < _numAxis; i++) {
346            _axisLen[i] = newAxisLen[i];
347        }
348    }
349    __getInterpPos();
350    __getInterpData();
351    fprintf(stdout, "----end interpolation----\n");
352    fflush(stdout);
353    return *this;
354}
355
356int
357DX::n() const
358{
359    return _n;
360}
361
362int
363DX::rank() const
364{
365    return _rank;
366}
367
368int
369DX::shape() const
370{
371    return _shape;
372}
373
374const float*
375DX::delta() const
376{
377    return _delta;
378}
379
380const float*
381DX::max() const
382{
383    return _max;
384}
385
386const float*
387DX::origin() const
388{
389    return _origin;
390}
391
392const float*
393DX::positions() const
394{
395    return _positions;
396}
397
398const int*
399DX::axisLen() const
400{
401    return _axisLen;
402}
403
404const float*
405DX::data() const
406{
407    return _data;
408}
409
410float
411DX::dataMin() const
412{
413    return _dataMin;
414}
415
416float
417DX::dataMax() const
418{
419    return _dataMax;
420}
421
422float
423DX::nzero_min() const
424{
425    return _nzero_min;
426}
Note: See TracBrowser for help on using the repository browser.