1 | #include <stdio.h> |
---|
2 | #include <math.h> |
---|
3 | #include <fstream> |
---|
4 | #include <iostream> |
---|
5 | #include <string> |
---|
6 | #include "RpFieldRect3D.h" |
---|
7 | #include "RpFieldPrism3D.h" |
---|
8 | |
---|
9 | /* Load a 3D volume from a dx-format file |
---|
10 | */ |
---|
11 | void |
---|
12 | load_volume_file(int index, char *fname) { |
---|
13 | Rappture::MeshTri2D xymesh; |
---|
14 | int dummy, nx, ny, nz, nxy, npts; |
---|
15 | double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz; |
---|
16 | char line[128], type[128], *start; |
---|
17 | std::ifstream fin(fname); |
---|
18 | |
---|
19 | int isrect = 1; |
---|
20 | |
---|
21 | do { |
---|
22 | fin.getline(line,sizeof(line)-1); |
---|
23 | for (start=&line[0]; *start == ' ' || *start == '\t'; start++) |
---|
24 | ; // skip leading blanks |
---|
25 | |
---|
26 | if (*start != '#') { // skip comment lines |
---|
27 | if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) { |
---|
28 | // found grid size |
---|
29 | isrect = 1; |
---|
30 | } |
---|
31 | else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) { |
---|
32 | std::ofstream ftmp("tmppts"); |
---|
33 | double xx, yy, zz; |
---|
34 | isrect = 0; |
---|
35 | for (int i=0; i < nxy; i++) { |
---|
36 | fin.getline(line,sizeof(line)-1); |
---|
37 | if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) { |
---|
38 | xymesh.addNode( Rappture::Node2D(xx,yy) ); |
---|
39 | ftmp << xx << " " << yy << std::endl; |
---|
40 | } |
---|
41 | } |
---|
42 | ftmp.close(); |
---|
43 | |
---|
44 | if (system("voronoi -t < tmppts > tmpcells") == 0) { |
---|
45 | int cx, cy, cz; |
---|
46 | std::ifstream ftri("tmpcells"); |
---|
47 | while (!ftri.eof()) { |
---|
48 | ftri.getline(line,sizeof(line)-1); |
---|
49 | if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) { |
---|
50 | xymesh.addCell(cx, cy, cz); |
---|
51 | } |
---|
52 | } |
---|
53 | ftri.close(); |
---|
54 | } else { |
---|
55 | std::cerr << "WARNING: triangularization failed" << std::endl; |
---|
56 | } |
---|
57 | } |
---|
58 | else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) { |
---|
59 | // found z-grid |
---|
60 | } |
---|
61 | else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) { |
---|
62 | // found origin |
---|
63 | } |
---|
64 | else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) { |
---|
65 | // found one of the delta lines |
---|
66 | if (ddx != 0.0) { dx = ddx; } |
---|
67 | else if (ddy != 0.0) { dy = ddy; } |
---|
68 | else if (ddz != 0.0) { dz = ddz; } |
---|
69 | } |
---|
70 | else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) { |
---|
71 | if (isrect && (npts != nx*ny*nz)) { |
---|
72 | std::cerr << "inconsistent data: expected " << nx*ny*nz << " points but found " << npts << " points" << std::endl; |
---|
73 | return; |
---|
74 | } |
---|
75 | else if (!isrect && (npts != nxy*nz)) { |
---|
76 | std::cerr << "inconsistent data: expected " << nxy*nz << " points but found " << npts << " points" << std::endl; |
---|
77 | return; |
---|
78 | } |
---|
79 | break; |
---|
80 | } |
---|
81 | } |
---|
82 | } while (!fin.eof()); |
---|
83 | |
---|
84 | // read data points |
---|
85 | if (!fin.eof()) { |
---|
86 | if (isrect) { |
---|
87 | Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx); |
---|
88 | Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny); |
---|
89 | Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz); |
---|
90 | Rappture::FieldRect3D field(xgrid, ygrid, zgrid); |
---|
91 | |
---|
92 | double dval; |
---|
93 | int nread = 0; |
---|
94 | while (!fin.eof() && nread < npts) { |
---|
95 | if (!(fin >> dval).fail()) { |
---|
96 | field.define(nread++, dval); |
---|
97 | } |
---|
98 | } |
---|
99 | |
---|
100 | // make sure that we read all of the expected points |
---|
101 | if (nread != nx*ny*nz) { |
---|
102 | std::cerr << "inconsistent data: expected " << nx*ny*nz << " points but found " << nread << " points" << std::endl; |
---|
103 | return; |
---|
104 | } |
---|
105 | |
---|
106 | // figure out a good mesh spacing |
---|
107 | int nsample = 100; |
---|
108 | dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis); |
---|
109 | dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis); |
---|
110 | dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis); |
---|
111 | double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333); |
---|
112 | |
---|
113 | nx = (int)ceil(dx/dmin); |
---|
114 | ny = (int)ceil(dy/dmin); |
---|
115 | nz = (int)ceil(dz/dmin); |
---|
116 | float *data = new float[nx*ny*nz]; |
---|
117 | std::cout << "generating " << nx << "x" << ny << "x" << nz << " = " << nx*ny*nz << " points" << std::endl; |
---|
118 | |
---|
119 | // generate the uniformly sampled data that we need for a volume |
---|
120 | int ngen = 0; |
---|
121 | for (int iz=0; iz < nz; iz++) { |
---|
122 | double zval = z0 + iz*dmin; |
---|
123 | for (int iy=0; iy < ny; iy++) { |
---|
124 | double yval = y0 + iy*dmin; |
---|
125 | for (int ix=0; ix < nx; ix++) { |
---|
126 | double xval = x0 + ix*dmin; |
---|
127 | data[ngen++] = field.value(xval,yval,zval); |
---|
128 | } |
---|
129 | } |
---|
130 | } |
---|
131 | } else { |
---|
132 | Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz); |
---|
133 | Rappture::FieldPrism3D field(xymesh, zgrid); |
---|
134 | |
---|
135 | double dval; |
---|
136 | int nread = 0; |
---|
137 | while (!fin.eof() && nread < npts) { |
---|
138 | if (!(fin >> dval).fail()) { |
---|
139 | field.define(nread++, dval); |
---|
140 | } |
---|
141 | } |
---|
142 | |
---|
143 | // make sure that we read all of the expected points |
---|
144 | if (nread != nxy*nz) { |
---|
145 | std::cerr << "inconsistent data: expected " << nxy*nz << " points but found " << nread << " points" << std::endl; |
---|
146 | return; |
---|
147 | } |
---|
148 | |
---|
149 | // figure out a good mesh spacing |
---|
150 | int nsample = 100; |
---|
151 | dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis); |
---|
152 | dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis); |
---|
153 | dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis); |
---|
154 | double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333); |
---|
155 | |
---|
156 | nx = (int)ceil(dx/dmin); |
---|
157 | ny = (int)ceil(dy/dmin); |
---|
158 | nz = (int)ceil(dz/dmin); |
---|
159 | float *data = new float[nx*ny*nz]; |
---|
160 | std::cout << "generating " << nx << "x" << ny << "x" << nz << " = " << nx*ny*nz << " points" << std::endl; |
---|
161 | |
---|
162 | // generate the uniformly sampled data that we need for a volume |
---|
163 | int ngen = 0; |
---|
164 | for (int iz=0; iz < nz; iz++) { |
---|
165 | double zval = z0 + iz*dmin; |
---|
166 | for (int iy=0; iy < ny; iy++) { |
---|
167 | double yval = y0 + iy*dmin; |
---|
168 | for (int ix=0; ix < nx; ix++) { |
---|
169 | double xval = x0 + ix*dmin; |
---|
170 | data[ngen++] = field.value(xval,yval,zval); |
---|
171 | } |
---|
172 | } |
---|
173 | } |
---|
174 | } |
---|
175 | } else { |
---|
176 | std::cerr << "WARNING: data not found in file " << fname << std::endl; |
---|
177 | } |
---|
178 | |
---|
179 | //load_volume(index, nx, ny, nz, 1, data); |
---|
180 | } |
---|
181 | |
---|
182 | int |
---|
183 | main(int argc, char **argv) |
---|
184 | { |
---|
185 | if (argc != 2) { |
---|
186 | std::cerr << "usage: viz file" << std::endl; |
---|
187 | exit(1); |
---|
188 | } |
---|
189 | load_volume_file(0, argv[1]); |
---|
190 | exit(0); |
---|
191 | } |
---|