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

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