source: trunk/packages/vizservers/nanovis/LoadVector.cpp @ 1374

Last change on this file since 1374 was 1374, checked in by gah, 13 years ago
File size: 9.4 KB
Line 
1
2#include <stdio.h>
3#include <string.h>
4#include <rappture.h>
5#include <tcl.h>
6#include "Unirect.h"
7#include "RpField1D.h"
8#include "RpFieldRect3D.h"
9#include "RpFieldPrism3D.h"
10#include "Nv.h"
11#include "nanovis.h"
12
13static INLINE char *
14skipspaces(char *string)
15{
16    while (isspace(*string)) {
17        string++;
18    }
19    return string;
20}
21
22static INLINE char *
23getline(char **stringPtr, char *endPtr)
24{
25    char *line, *p;
26
27    line = skipspaces(*stringPtr);
28    for (p = line; p < endPtr; p++) {
29        if (*p == '\n') {
30            *p++ = '\0';
31            *stringPtr = p;
32            return line;
33        }
34    }
35    *stringPtr = p;
36    return line;
37}
38bool
39MakeVectorField(Rappture::Outcome &result,
40                float xMin, float xMax, size_t xNum,
41                float yMin, float yMax, size_t yNum,
42                float zMin, float zMax, size_t zNum,
43                size_t nValues, float *values)
44{
45#ifdef notdef
46    double max_x = -1e21, min_x = 1e21;
47    double max_y = -1e21, min_y = 1e21;
48    double max_z = -1e21, min_z = 1e21;
49#endif
50    double max_mag = -1e21, min_mag = 1e21;
51    for (size_t i = 0; i < nValues; i += 3) {
52        double vx, vy, vz, vm;
53
54        vx = values[i];
55        vy = values[i+1];
56        vz = values[i+2];
57                   
58        vm = sqrt(vx*vx + vy*vy + vz*vz);
59        if (vm > max_mag) {
60            max_mag = vm;
61        }
62        if (vm < min_mag) {
63            min_mag = vm;
64        }
65    }
66   
67    size_t n = 4* xNum * yNum * zNum;
68    float *data = new float[n];
69    memset(data, 0, sizeof(float) * n);
70   
71    fprintf(stderr, "generating %dx%dx%d = %d points\n",
72            xNum, yNum, zNum, xNum * yNum * zNum);
73   
74    // Generate the uniformly sampled rectangle that we need for a volume
75    float *destPtr = data;
76    float *srcPtr = values;
77    for (size_t i = 0; i < zNum; i++) {
78        for (size_t j = 0; j < yNum; j++) {
79            for (size_t k = 0; k < xNum; k++) {
80                double vx, vy, vz, vm;
81
82                vx = srcPtr[0];
83                vy = srcPtr[1];
84                vz = srcPtr[2];
85               
86                vm = sqrt(vx*vx + vy*vy + vz*vz);
87               
88                destPtr[0] = vm / max_mag;
89                destPtr[1] = vx /(2.0*max_mag) + 0.5;
90                destPtr[2] = vy /(2.0*max_mag) + 0.5;
91                destPtr[3] = vz /(2.0*max_mag) + 0.5;
92                srcPtr += 3;
93                destPtr += 4;
94            }
95        }
96    }
97   
98    Volume *volPtr;
99    int ivol = NanoVis::n_volumes;
100    volPtr = NanoVis::load_volume(ivol, xNum, yNum, zNum, 4, data,
101        min_mag, max_mag, 0);
102    volPtr->xAxis.SetRange(xMin, xMax);
103    volPtr->yAxis.SetRange(yMin, yMax);
104    volPtr->zAxis.SetRange(zMin, zMax);
105    volPtr->wAxis.SetRange(min_mag, max_mag);
106    volPtr->update_pending = true;
107    delete [] data;
108    return true;
109}
110
111bool
112MakeResampledVectorField(Rappture::Outcome &result,
113                         float xMin, float xMax, size_t xNum,
114                         float yMin, float yMax, size_t yNum,
115                         float zMin, float zMax, size_t zNum,
116                         size_t nValues, float *values)
117{
118    Rappture::Mesh1D xgrid(xMin, xMax, xNum);
119    Rappture::Mesh1D ygrid(yMin, yMax, yNum);
120    Rappture::Mesh1D zgrid(zMin, zMax, zNum);
121    Rappture::FieldRect3D xfield(xgrid, ygrid, zgrid);
122    Rappture::FieldRect3D yfield(xgrid, ygrid, zgrid);
123    Rappture::FieldRect3D zfield(xgrid, ygrid, zgrid);
124
125#ifdef notdef
126    double max_x = -1e21, min_x = 1e21;
127    double max_y = -1e21, min_y = 1e21;
128    double max_z = -1e21, min_z = 1e21;
129#endif
130
131    double max_mag, min_mag, nzero_min;
132    max_mag = -1e21, nzero_min = min_mag = 1e21;
133    size_t i, j;
134    for (i = j = 0; i < nValues; i += 3, j++) {
135        double vx, vy, vz, vm;
136
137        vx = values[i];
138        vy = values[i+1];
139        vz = values[i+2];
140       
141        xfield.define(j, vx);
142        yfield.define(j, vy);
143        zfield.define(j, vz);
144
145        vm = sqrt(vx*vx + vy*vy + vz*vz);
146        if (vm > max_mag) {
147            max_mag = vm;
148        }
149        if (vm < min_mag) {
150            min_mag = vm;
151        }
152        if ((vm != 0.0f) && (vm < nzero_min)) {
153            nzero_min = vm;
154        }
155    }
156   
157    // figure out a good mesh spacing
158    int nSamples = 30;
159    double dx, dy, dz;
160    dx = xfield.rangeMax(Rappture::xaxis) - xfield.rangeMin(Rappture::xaxis);
161    dy = xfield.rangeMax(Rappture::yaxis) - xfield.rangeMin(Rappture::yaxis);
162    dz = xfield.rangeMax(Rappture::zaxis) - xfield.rangeMin(Rappture::zaxis);
163
164    double dmin;
165    dmin = pow((dx*dy*dz)/(nSamples*nSamples*nSamples), 0.333);
166   
167    printf("dx:%lf dy:%lf dz:%lf dmin:%lf\n", dx, dy, dz, dmin);
168   
169    /* Recompute new number of points for each axis. */
170    xNum = (int)ceil(dx/dmin);
171    yNum = (int)ceil(dy/dmin);
172    zNum = (int)ceil(dz/dmin);
173   
174#ifndef NV40
175    // must be an even power of 2 for older cards
176    xNum = (int)pow(2.0, ceil(log10((double)xNum)/log10(2.0)));
177    yNum = (int)pow(2.0, ceil(log10((double)yNum)/log10(2.0)));
178    zNum = (int)pow(2.0, ceil(log10((double)zNum)/log10(2.0)));
179#endif
180
181    size_t n = 4 * xNum * yNum * zNum;
182    float *data = new float[n];
183    memset(data, 0, sizeof(float) * n);
184   
185    fprintf(stderr, "generating %dx%dx%d = %d points\n",
186            xNum, yNum, zNum, xNum * yNum * zNum);
187   
188    // Generate the uniformly sampled rectangle that we need for a volume
189    float *destPtr = data;
190    for (size_t i = 0; i < zNum; i++) {
191        double z;
192
193        z = zMin + (i * dmin);
194        for (size_t j = 0; j < yNum; j++) {
195            double y;
196               
197            y = yMin + (j * dmin);
198            for (size_t k = 0; k < xNum; k++) {
199                double x;
200                double vx, vy, vz, vm;
201
202                x = xMin + (k * dmin);
203                vx = xfield.value(x, y, z);
204                vy = yfield.value(x, y, z);
205                vz = zfield.value(x, y, z);
206                vm = sqrt(vx*vx + vy*vy + vz*vz);
207               
208                destPtr[0] = vm / max_mag;
209                destPtr[1] = vx /(2.0*max_mag) + 0.5;
210                destPtr[2] = vy /(2.0*max_mag) + 0.5;
211                destPtr[3] = vz /(2.0*max_mag) + 0.5;
212                destPtr += 4;
213            }
214        }
215    }
216   
217    Volume *volPtr;
218    int ivol = NanoVis::n_volumes;
219    volPtr = NanoVis::load_volume(ivol, xNum, yNum, zNum, 4, data,
220        min_mag, max_mag, nzero_min);
221    volPtr->xAxis.SetRange(xMin, xMax);
222    volPtr->yAxis.SetRange(yMin, yMax);
223    volPtr->zAxis.SetRange(zMin, zMax);
224    volPtr->wAxis.SetRange(min_mag, max_mag);
225    volPtr->update_pending = true;
226    delete [] data;
227    return true;
228}
229
230bool
231MakeVectorFieldFromUnirect3d(Rappture::Outcome &result,
232                             Rappture::Unirect3d &grid)
233{
234    return MakeVectorField(result,
235        grid.xMin(), grid.xMax(), grid.xNum(),
236        grid.yMin(), grid.yMax(), grid.yNum(),
237        grid.zMin(), grid.zMax(), grid.zNum(),
238        grid.nValues(), grid.values());
239}
240
241
242#ifdef notdef
243/*
244 * Load a 3D vector field from a dx-format file
245 */
246bool
247load_vector_stream2(Rappture::Outcome &result, int ivol, size_t length,
248                    char *string)
249{
250    int nx, ny, nz, npts;
251    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
252    char *p, *endPtr;
253
254    dx = dy = dz = 0.0;         // Suppress compiler warning.
255    x0 = y0 = z0 = 0.0;         // May not have an origin line.
256    for (p = string, endPtr = p + length; p < endPtr; /*empty*/) {
257        char *line;
258
259        line = getline(&p, endPtr);
260        if (line == endPtr) {
261            break;
262        }
263        if (*line == '#') {  // skip comment lines
264            continue;
265        }
266        if (sscanf(line, "object %*d class gridpositions counts %d %d %d",
267                   &nx, &ny, &nz) == 4) {
268            printf("w:%d h:%d d:%d\n", nx, ny, nz);
269            // found grid size
270        } else if (sscanf(line, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
271            // found origin
272        } else if (sscanf(line, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
273            // found one of the delta lines
274            if (ddx != 0.0) {
275                dx = ddx;
276            } else if (ddy != 0.0) {
277                dy = ddy;
278            } else if (ddz != 0.0) {
279                dz = ddz;
280            }
281        } else if (sscanf(line, "object %*d class array type %*s shape 3"
282                " rank 1 items %d data follows", &npts) == 3) {
283            printf("point %d\n", npts);
284            if (npts != nx*ny*nz) {
285                result.addError("inconsistent data: expected %d points"
286                                " but found %d points", nx*ny*nz, npts);
287                return false;
288            }
289            break;
290        } else if (sscanf(line, "object %*d class array type %*s rank 0"
291                " times %d data follows", &npts) == 3) {
292            if (npts != nx*ny*nz) {
293                result.addError("inconsistent data: expected %d points"
294                                " but found %d points", nx*ny*nz, npts);
295                return false;
296            }
297            break;
298        }
299    }
300
301    // read data points
302    float* srcData = new float[nx * ny * nz * 3];
303    if (p >= endPtr) {
304        std::cerr << "WARNING: data not found in stream" << std::endl;
305        return true;
306    }
307#ifdef notdef
308    double max_x = -1e21, min_x = 1e21;
309    double max_y = -1e21, min_y = 1e21;
310    double max_z = -1e21, min_z = 1e21;
311#endif
312    double max_mag = -1e21, min_mag = 1e21;
313    int nRead = 0;
314    for (int ix=0; ix < nx; ix++) {
315        for (int iy=0; iy < ny; iy++) {
316            for (int iz=0; iz < nz; iz++) {
317                char *line;
318                double vx, vy, vz, vm;
319
320                if ((p == endPtr) || nRead > npts) {
321                    break;
322                }
323                line = getline(&p, endPtr);
324                if (sscanf(line, "%lg %lg %lg", &vx, &vy, &vz) == 3) {
325                    int nindex = (iz*nx*ny + iy*nx + ix) * 3;
326                    srcData[nindex] = vx;
327                    //if (srcData[nindex] > max_x) max_x = srcData[nindex];
328                    //if (srcData[nindex] < min_x) min_x = srcData[nindex];
329                    ++nindex;
330                   
331                    srcData[nindex] = vy;
332                    //if (srcData[nindex] > max_y) max_y = srcData[nindex];
333                    //if (srcData[nindex] < min_y) min_y = srcData[nindex];
334                    ++nindex;
335                   
336                    srcData[nindex] = vz;
337                    //if (srcData[nindex] > max_z) max_z = srcData[nindex];
338                    //if (srcData[nindex] < min_z) min_z = srcData[nindex];
339                   
340                    vm = sqrt(vx*vx + vy*vy + vz*vz);
341                    if (vm > max_mag) {
342                        max_mag = vm;
343                    }
344                    if (vm < min_mag) {
345                        min_mag = vm;
346                    }
347                    ++nRead;
348                }
349            }
350        }
351    }
352   
353    // make sure that we read all of the expected points
354    if (nRead != nx*ny*nz) {
355        result.addError("inconsistent data: expected %d points"
356                        " but found %d points", nx*ny*nz, npts);
357        return false;
358    }
359    bool status;
360    status = MakeVectorField(result, x0, x0 + nx, nx, y0, y0 + ny, ny,
361        z0, z0 + ny, ny, nRead, srcData);
362    delete [] srcData;
363    return status;
364}
365
366#endif
Note: See TracBrowser for help on using the repository browser.