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 | |
---|
14 | static INLINE char * |
---|
15 | skipspaces(char *string) |
---|
16 | { |
---|
17 | while (isspace(*string)) { |
---|
18 | string++; |
---|
19 | } |
---|
20 | return string; |
---|
21 | } |
---|
22 | |
---|
23 | static INLINE char * |
---|
24 | getline(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 | |
---|
40 | bool |
---|
41 | SetVectorFieldData(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 | |
---|
114 | bool |
---|
115 | SetResampledVectorFieldData(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 | |
---|
235 | bool |
---|
236 | SetVectorFieldDataFromUnirect3d(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 | */ |
---|
252 | bool |
---|
253 | load_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 |
---|