source: trunk/vizservers/nanovis/RpDX.cpp @ 973

Last change on this file since 973 was 973, checked in by gah, 13 years ago
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) :
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        // error
50    }
51    // open the file with libdx
52    fprintf(stdout, "Calling DXImportDX(%s)\n", filename);
53    fflush(stdout);
54    DXenable_locks(0);
55    _dxobj = DXImportDX((char*)filename,NULL,NULL,NULL,NULL);
56
57    fprintf(stdout, "dxobj=%x\n", (unsigned int)_dxobj);
58    fflush(stdout);
59
60    // parse out the positions array
61    // FIXME: nanowire will need a different way to parse out the positions array
62    // since it uses a productarray to store its positions.
63    // possibly use DXGetProductArray()
64    dxpos = (Array) DXGetComponentValue((Field) _dxobj, (char *)"positions");
65    DXGetArrayInfo(dxpos, &_n, &type, &category, &_rank, &_shape);
66
67    fprintf(stdout, "_n = %d\n",_n);
68    if (type != TYPE_FLOAT) {
69        fprintf(stderr, "WARNING: positions is not type float (type=%d)\n",
70            type);
71        fflush(stderr);
72    }
73    fprintf(stdout, "_rank = %d\n",_rank);
74    fprintf(stdout, "_shape = %d\n",_shape);
75
76    float* pos = NULL;
77    pos = (float*) DXGetArrayData(dxpos);
78    if (pos == NULL) {
79        fprintf(stdout, "DXGetArrayData failed to return positions array\n");
80        fflush(stdout);
81    }
82
83    // first call to get the number of axis needed
84    dxgrid = (Array) DXGetComponentValue((Field) _dxobj, (char *)"connections");
85    if (DXQueryGridConnections(dxgrid, &_numAxis, NULL) == NULL) {
86        // raise error, data is not a regular grid and we cannot handle it
87        fprintf(stdout,"DX says our grid is not regular, we cannot handle this data\n");
88        fflush(stdout);
89    }
90
91    _positions = new float[_n*_numAxis];
92    if (_positions == NULL) {
93        // malloc failed, raise error
94        fprintf(stdout, "malloc of _positions array failed");
95        fflush(stdout);
96    }
97    memcpy(_positions,pos,sizeof(float)*_n*_numAxis);
98
99    _axisLen = new int[_numAxis];
100    if (_axisLen == NULL) {
101        // malloc failed, raise error
102        fprintf(stdout, "malloc of _axisLen array failed");
103        fflush(stdout);
104    }
105    memset(_axisLen, 0, _numAxis);
106
107    _delta = new float[_numAxis*_numAxis];
108    if (_delta == NULL) {
109        // malloc failed, raise error
110        fprintf(stdout, "malloc of _delta array failed");
111        fflush(stdout);
112    }
113    memset(_delta, 0, _numAxis*_numAxis);
114
115    _origin = new float[_numAxis];
116    if (_origin == NULL) {
117        // malloc failed, raise error
118        fprintf(stdout, "malloc of _origin array failed");
119        fflush(stdout);
120    }
121    memset(_origin, 0, _numAxis);
122
123    _max = new float[_numAxis];
124    if (_max == NULL) {
125        // malloc failed, raise error
126        fprintf(stdout, "malloc of _max array failed");
127        fflush(stdout);
128    }
129    memset(_max, 0, _numAxis);
130
131    __findPosMax();
132
133    // parse out the gridconnections (length of each axis) array
134    // dxgrid = (Array) DXQueryGridConnections(dxpos, &_numAxis, _axisLen);
135    DXQueryGridPositions(dxpos, NULL, _axisLen, _origin, _delta);
136
137    fprintf(stdout, "_max = [%g,%g,%g]\n",_max[0],_max[1],_max[2]);
138    fprintf(stdout, "_delta = [%g,%g,%g]\n",_delta[0],_delta[1],_delta[2]);
139    fprintf(stdout, "         [%g,%g,%g]\n",_delta[3],_delta[4],_delta[5]);
140    fprintf(stdout, "         [%g,%g,%g]\n",_delta[6],_delta[7],_delta[8]);
141    fprintf(stdout, "_origin = [%g,%g,%g]\n",_origin[0],_origin[1],_origin[2]);
142    fprintf(stdout, "_axisLen = [%i,%i,%i]\n",_axisLen[0],_axisLen[1],_axisLen[2]);
143    fflush(stdout);
144
145    // grab the data array from the dx object and store it in _data
146    dxdata = (Array) DXGetComponentValue((Field) _dxobj, (char *)"data");
147    DXGetArrayInfo(dxdata, NULL, &type, NULL, NULL, NULL);
148    _data = new float[_n];
149    if (_data == NULL) {
150        // malloc failed, raise error
151        fprintf(stdout, "malloc of _data array failed");
152        fflush(stdout);
153    }
154
155    switch (type) {
156    case TYPE_FLOAT:
157        float *float_data;
158
159        float_data = (float*) DXGetArrayData(dxdata);
160        memcpy(_data, float_data, sizeof(float)*_n);
161        break;
162    case TYPE_DOUBLE:
163        double *double_data;
164        double_data = (double*) DXGetArrayData(dxdata);
165        for (int i = 0; i < _n; i++) {
166            _data[i] = double_data[i];
167        }
168    default:
169        fprintf(stdout, "don't know how to handle data of type %d\n", type);
170        fflush(stdout);
171    }
172
173    // print debug info
174    for (int lcv = 0, pt = 0; lcv < _n; lcv +=3, pt+=9) {
175        fprintf(stdout,
176            "(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n",
177            _positions[pt],_positions[pt+1],_positions[pt+2], _data[lcv],
178            _positions[pt+3],_positions[pt+4],_positions[pt+5],_data[lcv+1],
179            _positions[pt+6],_positions[pt+7],_positions[pt+8],_data[lcv+2]);
180        fflush(stdout);
181    }
182    __collectDataStats();
183}
184
185DX::~DX()
186{
187    delete[] _axisLen;
188    delete[] _delta;
189    delete[] _origin;
190    delete[] _max;
191    delete[] _data;
192    delete[] _positions;
193}
194
195void
196DX::__findPosMax()
197{
198    int lcv = 0;
199
200    // initialize the max array and delta array
201    // max holds the maximum value found for each index
202    for (lcv = 0; lcv < _numAxis; lcv++) {
203        _max[lcv] = _positions[lcv];
204    }
205
206    for (lcv=lcv; lcv < _n*_numAxis; lcv++) {
207        int maxIdx = lcv%_numAxis;
208        if (_positions[lcv] > _max[maxIdx]) {
209            _max[maxIdx] = _positions[lcv];
210        }
211    }
212}
213
214void
215DX::__collectDataStats()
216{
217    _dataMin = FLT_MAX;
218    _dataMax = -FLT_MAX;
219    _nzero_min = FLT_MAX;
220
221    // populate _dataMin, _dataMax, _nzero_min
222    for (int lcv = 0; lcv < _n; lcv++) {
223        if (_data[lcv] < _dataMin) {
224            _dataMin = _data[lcv];
225        }
226        if (_data[lcv] > _dataMax) {
227            _dataMax = _data[lcv];
228        }
229        if ((_data[lcv] > 0.0f) && (_data[lcv] < _nzero_min)) {
230            _nzero_min = _data[lcv];
231        }
232    }
233    if (_nzero_min == FLT_MAX) {
234        fprintf(stderr, "could not find a positive minimum value\n");
235        fflush(stderr);
236    }
237}
238
239/*
240 * getInterpPos()
241 *
242 * creates a new grid of positions which can be used for interpolation.
243 * we create the positions array using the function DXMakeGridPositionsV
244 * which creates an n-dimensional grid of regularly spaced positions.
245 * This function overwrites the original positions array
246 */
247void
248DX::__getInterpPos()
249{
250    Array dxpos;
251    float* pos = NULL;
252
253    // gather the positions we want to interpolate over
254    dxpos = DXMakeGridPositionsV(_numAxis, _axisLen, _origin, _delta);
255    DXGetArrayInfo(dxpos, &_n, NULL, NULL, &_rank, &_shape);
256    pos = (float*) DXGetArrayData(dxpos);
257    if (pos == NULL) {
258        fprintf(stdout, "DXGetArrayData failed to return positions array\n");
259        fflush(stdout);
260    }
261
262    if (_positions != NULL) {
263        delete[] _positions;
264    }
265    _positions = new float[_n*_numAxis];
266    if (_positions == NULL) {
267        // malloc failed, raise error
268        fprintf(stdout, "malloc of _axisLen array failed");
269        fflush(stdout);
270    }
271    memcpy(_positions,pos,sizeof(float)*_n*_numAxis);
272
273    pos = NULL;
274    DXDelete((object*)dxpos);
275}
276
277/*
278 * getInterpData()
279 *
280 * this function interpolates over a positions array to produce data for each
281 * point in the positions array. we use the position data stored in _positions
282 * array.
283 */
284void
285DX::__getInterpData()
286{
287    int pts = _n;
288    int interppts = pts;
289    Interpolator interpolator;
290
291    _data = new float[_n];
292    if (_data == NULL) {
293        // malloc failed, raise error
294        fprintf(stdout, "malloc of _data array failed");
295        fflush(stdout);
296    }
297    memset(_data,0,_n);
298
299    // build the interpolator and interpolate
300    fprintf(stdout, "creating DXNewInterpolator...\n");
301    fflush(stdout);
302    interpolator = DXNewInterpolator(_dxobj,INTERP_INIT_IMMEDIATE,-1.0);
303    fprintf(stdout, "DXNewInterpolator=%x\n", (unsigned int)interpolator);
304    fprintf(stdout,"_rank = %i\n",_rank);
305    fprintf(stdout,"_shape = %i\n",_shape);
306    fprintf(stdout,"_n = %i\n",_n);
307    fprintf(stdout,"start interppts = %i\n",interppts);
308    fflush(stdout);
309    DXInterpolate(interpolator,&interppts,_positions,_data);
310    fprintf(stdout,"interppts = %i\n",interppts);
311    fflush(stdout);
312
313    // print debug info
314    for (int lcv = 0, pt = 0; lcv < pts; lcv+=3,pt+=9) {
315        fprintf(stdout,
316            "(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n",
317            _positions[pt],_positions[pt+1],_positions[pt+2], _data[lcv],
318            _positions[pt+3],_positions[pt+4],_positions[pt+5],_data[lcv+1],
319            _positions[pt+6],_positions[pt+7],_positions[pt+8],_data[lcv+2]);
320        fflush(stdout);
321    }
322
323    __collectDataStats();
324}
325
326/*
327 * interpolate()
328 *
329 * generate a positions array with optional new axis length and
330 * interpolate to get the new data values at each point in
331 * the positions array. this function currently only works if you
332 * do not change the axis length (i.e. newAxisLen == NULL).
333 */
334DX&
335DX::interpolate(int* newAxisLen)
336{
337    fprintf(stdout, "----begin interpolation----\n");
338    fflush(stdout);
339    if (newAxisLen != NULL) {
340        for (int i = 0; i < _numAxis; i++) {
341            _axisLen[i] = newAxisLen[i];
342        }
343    }
344    __getInterpPos();
345    __getInterpData();
346    fprintf(stdout, "----end interpolation----\n");
347    fflush(stdout);
348    return *this;
349}
350
351int
352DX::n() const
353{
354    return _n;
355}
356
357int
358DX::rank() const
359{
360    return _rank;
361}
362
363int
364DX::shape() const
365{
366    return _shape;
367}
368
369const float*
370DX::delta() const
371{
372    return _delta;
373}
374
375const float*
376DX::max() const
377{
378    return _max;
379}
380
381const float*
382DX::origin() const
383{
384    return _origin;
385}
386
387const float*
388DX::positions() const
389{
390    return _positions;
391}
392
393const int*
394DX::axisLen() const
395{
396    return _axisLen;
397}
398
399const float*
400DX::data() const
401{
402    return _data;
403}
404
405float
406DX::dataMin() const
407{
408    return _dataMin;
409}
410
411float
412DX::dataMax() const
413{
414    return _dataMax;
415}
416
417float
418DX::nzero_min() const
419{
420    return _nzero_min;
421}
Note: See TracBrowser for help on using the repository browser.