Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Annotation of /branches/ein16/synth/d2/symb/symb_fm3x3d3.diderot
ViewVC logotype

Annotation of /branches/ein16/synth/d2/symb/symb_fm3x3d3.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4444 - (view) (download)

1 : cchiw 4236 /* fs3d-vec.diderot: function sampler, 3D vector synthetic fields
2 :    
3 :     This is based on fs3d-scl.diderot in this same directory, after much
4 :     copy-and-pasting. As with that program, the `-which` option will
5 :     determine which function is sampled; look for `(0 == which)` in the
6 :     code to see the start of the function definitions.
7 :    
8 :     Assuming the directions at https://github.com/Diderot-Language/examples
9 :     this program can be compiled with:
10 :    
11 :     ../../vis12/bin/diderotc --exec fs3d-vec.diderot
12 :    
13 :     This program's printed output needs to be captured to generate
14 :     NRRD header with full orientation info. A vector-valued identity
15 :     function is sampled via:
16 :    
17 :     ./fs2d-vec -which 0 | unu save -f nrrd -o ident.nrrd
18 :     rm out.nrrd
19 :    
20 :     Note that like in [`fs2d-scl.diderot`](../fs2d), the NRRD header
21 :     generated by this program assumes:
22 :    
23 :     1. This program was not compiled with `--double`
24 :     2. The program is running on a little-endian machine.
25 :     */
26 :    
27 :     input int sz0 ("# samples on fastest axis") = 52;
28 :     input int sz1 ("# samples on medium axis") = 51;
29 :     input int sz2 ("# samples on slowest axis") = 50;
30 :     input real width ("width of edge of cube region sampled") = 4;
31 :     input vec3 axis ("axis (non-normalized) of rotation of sampling grid") = [1,1,1];
32 :     input real angle ("angle (in degrees) of rotation of sampling grid") = 0;
33 :     input vec3 off ("translation offset, in index space, from origin-centered grid") = [0,0,0];
34 :     input vec3 parm0 ("parameters that functions may use") = [0,0,0];
35 :     input vec3 parm1 ("more parameters that functions may use") = [0,0,0];
36 :     input vec3 parm2 ("even more parameters that functions may use") = [0,0,0];
37 :    
38 :     real theta = angle*π/180;
39 :     // unit-length axis of rotation scaled by sin(theta/2)
40 :     vec3 snax = sin(theta/2)*normalize(axis);
41 :     // quaternion representing rotation
42 :     vec4 qq = [cos(theta/2), snax[0], snax[1], snax[2]];
43 :     // create rotation matrix by converstion from quaternion
44 :     // (which aren't currently natively supported in Diderot)
45 :     tensor[3,3] rot = [[qq[0]*qq[0] + qq[1]*qq[1] - qq[2]*qq[2] - qq[3]*qq[3],
46 :     2*(qq[1]*qq[2] - qq[0]*qq[3]),
47 :     2*(qq[1]*qq[3] + qq[0]*qq[2])],
48 :     [2*(qq[1]*qq[2] + qq[0]*qq[3]),
49 :     qq[0]*qq[0] - qq[1]*qq[1] + qq[2]*qq[2] - qq[3]*qq[3],
50 :     2*(qq[2]*qq[3] - qq[0]*qq[1])],
51 :     [2*(qq[1]*qq[3] - qq[0]*qq[2]),
52 :     2*(qq[2]*qq[3] + qq[0]*qq[1]),
53 :     qq[0]*qq[0] - qq[1]*qq[1] - qq[2]*qq[2] + qq[3]*qq[3]]];
54 :     // per axis inter-sample edge vectors
55 :     real space0 = width/(sz0-1);
56 :     real space1 = width/(sz1-1);
57 :     real space2 = width/(sz2-1);
58 :     vec3 edge0 = rot•[space0, 0, 0];
59 :     vec3 edge1 = rot•[0, space1, 0];
60 :     vec3 edge2 = rot•[0, 0, space2];
61 :     // offset in world space
62 :     vec3 offws = [off[0]*space0, off[1]*space1, off[2]*space2];
63 :     // location of first sample point
64 :     vec3 orig = -(edge0*(sz0-1) + edge1*(sz1-1) + edge2*(sz2-1))/2 + offws;
65 :    
66 :     //-- field defined by coefficients --
67 :     // z=1
68 :     input vec2 base1_z0 = [0, 0]; // b*x, c*y
69 :     input vec2 xsq1_z0 = [0, 0]; // d*x^2, g*y*x^2
70 :     input vec2 ysq1_z0 = [0, 0]; // f*y^2, h*x*y^2
71 :     input vec3 diag1_z0 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
72 :     // z=z
73 :     input vec2 base1_z1 = [0, 0]; // b*x, c*y
74 :     input vec2 xsq1_z1 = [0, 0]; // d*x^2, g*y*x^2
75 :     input vec2 ysq1_z1 = [0, 0]; // f*y^2, h*x*y^2
76 :     input vec3 diag1_z1 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
77 :     // z=z*z
78 :     input vec2 base1_z2 = [0, 0]; // b*x, c*y
79 :     input vec2 xsq1_z2 = [0, 0]; // d*x^2, g*y*x^2
80 :     input vec2 ysq1_z2 = [0, 0]; // f*y^2, h*x*y^2
81 :     input vec3 diag1_z2 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
82 :    
83 :     // z=1
84 :     input vec2 base2_z0 = [0, 0]; // b*x, c*y
85 :     input vec2 xsq2_z0 = [0, 0]; // d*x^2, g*y*x^2
86 :     input vec2 ysq2_z0 = [0, 0]; // f*y^2, h*x*y^2
87 :     input vec3 diag2_z0 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
88 :     // z=z
89 :     input vec2 base2_z1 = [0, 0]; // b*x, c*y
90 :     input vec2 xsq2_z1 = [0, 0]; // d*x^2, g*y*x^2
91 :     input vec2 ysq2_z1 = [0, 0]; // f*y^2, h*x*y^2
92 :     input vec3 diag2_z1 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
93 :     // z=z*z
94 :     input vec2 base2_z2 = [0, 0]; // b*x, c*y
95 :     input vec2 xsq2_z2 = [0, 0]; // d*x^2, g*y*x^2
96 :     input vec2 ysq2_z2 = [0, 0]; // f*y^2, h*x*y^2
97 :     input vec3 diag2_z2 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
98 :    
99 : cchiw 4385
100 : cchiw 4236 // z=1
101 :     input vec2 base3_z0 = [0, 0]; // b*x, c*y
102 :     input vec2 xsq3_z0 = [0, 0]; // d*x^2, g*y*x^2
103 :     input vec2 ysq3_z0 = [0, 0]; // f*y^2, h*x*y^2
104 :     input vec3 diag3_z0 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
105 :     // z=z
106 :     input vec2 base3_z1 = [0, 0]; // b*x, c*y
107 :     input vec2 xsq3_z1 = [0, 0]; // d*x^2, g*y*x^2
108 :     input vec2 ysq3_z1 = [0, 0]; // f*y^2, h*x*y^2
109 :     input vec3 diag3_z1 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
110 :     // z=z*z
111 :     input vec2 base3_z2 = [0, 0]; // b*x, c*y
112 :     input vec2 xsq3_z2 = [0, 0]; // d*x^2, g*y*x^2
113 :     input vec2 ysq3_z2 = [0, 0]; // f*y^2, h*x*y^2
114 :     input vec3 diag3_z2 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
115 : cchiw 4385
116 :     //-- field defined by coefficients --
117 :    
118 :     // z=1
119 :     input vec2 base4_z0 = [0, 0]; // b*x, c*y
120 :     input vec2 xsq4_z0 = [0, 0]; // d*x^2, g*y*x^2
121 :     input vec2 ysq4_z0 = [0, 0]; // f*y^2, h*x*y^2
122 :     input vec3 diag4_z0 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
123 :     // z=z
124 :     input vec2 base4_z1 = [0, 0]; // b*x, c*y
125 :     input vec2 xsq4_z1 = [0, 0]; // d*x^2, g*y*x^2
126 :     input vec2 ysq4_z1 = [0, 0]; // f*y^2, h*x*y^2
127 :     input vec3 diag4_z1 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
128 :     // z=z*z
129 :     input vec2 base4_z2 = [0, 0]; // b*x, c*y
130 :     input vec2 xsq4_z2 = [0, 0]; // d*x^2, g*y*x^2
131 :     input vec2 ysq4_z2 = [0, 0]; // f*y^2, h*x*y^2
132 :     input vec3 diag4_z2 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
133 :    
134 :     input vec2 base5_z0 = [0, 0]; // b*x, c*y
135 :     input vec2 xsq5_z0 = [0, 0]; // d*x^2, g*y*x^2
136 :     input vec2 ysq5_z0 = [0, 0]; // f*y^2, h*x*y^2
137 :     input vec3 diag5_z0 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
138 :     // z=z
139 :     input vec2 base5_z1 = [0, 0]; // b*x, c*y
140 :     input vec2 xsq5_z1 = [0, 0]; // d*x^2, g*y*x^2
141 :     input vec2 ysq5_z1 = [0, 0]; // f*y^2, h*x*y^2
142 :     input vec3 diag5_z1 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
143 :     // z=z*z
144 :     input vec2 base5_z2 = [0, 0]; // b*x, c*y
145 :     input vec2 xsq5_z2 = [0, 0]; // d*x^2, g*y*x^2
146 :     input vec2 ysq5_z2 = [0, 0]; // f*y^2, h*x*y^2
147 :     input vec3 diag5_z2 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
148 :    
149 :     input vec2 base6_z0 = [0, 0]; // b*x, c*y
150 :     input vec2 xsq6_z0 = [0, 0]; // d*x^2, g*y*x^2
151 :     input vec2 ysq6_z0 = [0, 0]; // f*y^2, h*x*y^2
152 :     input vec3 diag6_z0 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
153 :     // z=z
154 :     input vec2 base6_z1 = [0, 0]; // b*x, c*y
155 :     input vec2 xsq6_z1 = [0, 0]; // d*x^2, g*y*x^2
156 :     input vec2 ysq6_z1 = [0, 0]; // f*y^2, h*x*y^2
157 :     input vec3 diag6_z1 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
158 :     // z=z*z
159 :     input vec2 base6_z2 = [0, 0]; // b*x, c*y
160 :     input vec2 xsq6_z2 = [0, 0]; // d*x^2, g*y*x^2
161 :     input vec2 ysq6_z2 = [0, 0]; // f*y^2, h*x*y^2
162 :     input vec3 diag6_z2 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
163 :    
164 :     input vec2 base7_z0 = [0, 0]; // b*x, c*y
165 :     input vec2 xsq7_z0 = [0, 0]; // d*x^2, g*y*x^2
166 :     input vec2 ysq7_z0 = [0, 0]; // f*y^2, h*x*y^2
167 :     input vec3 diag7_z0 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
168 :     // z=z
169 :     input vec2 base7_z1 = [0, 0]; // b*x, c*y
170 :     input vec2 xsq7_z1 = [0, 0]; // d*x^2, g*y*x^2
171 :     input vec2 ysq7_z1 = [0, 0]; // f*y^2, h*x*y^2
172 :     input vec3 diag7_z1 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
173 :     // z=z*z
174 :     input vec2 base7_z2 = [0, 0]; // b*x, c*y
175 :     input vec2 xsq7_z2 = [0, 0]; // d*x^2, g*y*x^2
176 :     input vec2 ysq7_z2 = [0, 0]; // f*y^2, h*x*y^2
177 :     input vec3 diag7_z2 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
178 :    
179 :     input vec2 base8_z0 = [0, 0]; // b*x, c*y
180 :     input vec2 xsq8_z0 = [0, 0]; // d*x^2, g*y*x^2
181 :     input vec2 ysq8_z0 = [0, 0]; // f*y^2, h*x*y^2
182 :     input vec3 diag8_z0 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
183 :     // z=z
184 :     input vec2 base8_z1 = [0, 0]; // b*x, c*y
185 :     input vec2 xsq8_z1 = [0, 0]; // d*x^2, g*y*x^2
186 :     input vec2 ysq8_z1 = [0, 0]; // f*y^2, h*x*y^2
187 :     input vec3 diag8_z1 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
188 :     // z=z*z
189 :     input vec2 base8_z2 = [0, 0]; // b*x, c*y
190 :     input vec2 xsq8_z2 = [0, 0]; // d*x^2, g*y*x^2
191 :     input vec2 ysq8_z2 = [0, 0]; // f*y^2, h*x*y^2
192 :     input vec3 diag8_z2 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
193 :    
194 :     input vec2 base9_z0 = [0, 0]; // b*x, c*y
195 :     input vec2 xsq9_z0 = [0, 0]; // d*x^2, g*y*x^2
196 :     input vec2 ysq9_z0 = [0, 0]; // f*y^2, h*x*y^2
197 :     input vec3 diag9_z0 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
198 :     // z=z
199 :     input vec2 base9_z1 = [0, 0]; // b*x, c*y
200 :     input vec2 xsq9_z1 = [0, 0]; // d*x^2, g*y*x^2
201 :     input vec2 ysq9_z1 = [0, 0]; // f*y^2, h*x*y^2
202 :     input vec3 diag9_z1 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
203 :     // z=z*z
204 :     input vec2 base9_z2 = [0, 0]; // b*x, c*y
205 :     input vec2 xsq9_z2 = [0, 0]; // d*x^2, g*y*x^2
206 :     input vec2 ysq9_z2 = [0, 0]; // f*y^2, h*x*y^2
207 :     input vec3 diag9_z2 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
208 :    
209 : cchiw 4236 input int cat=0;
210 :    
211 :    
212 :    
213 :     function real cvt(real x, real y, real z, vec2 base, vec2 xsq, vec2 ysq, vec3 diag){
214 : cchiw 4444 real b = base[0];
215 :     real c = base[1];
216 :     real d = xsq[0];
217 :     real g = xsq[1];
218 :     real f = ysq[0];
219 :     real h = ysq[1];
220 :     real a = diag[0];
221 :     real e = diag[1];
222 :     real i = diag[2];
223 :     real t0 = b*x + c*y;
224 :     real t1 = d*x^2 + g*y*x^2;
225 :     real t2 = f*y^2 + h*x*y^2;
226 :     real t3 = a + e*x*y + i*(y^2)*(x^2);
227 :     return z*(t0+t1+t2+t3);
228 : cchiw 4236 }
229 :    
230 :    
231 :     strand sample(int idx0, int idx1, int idx2) {
232 :    
233 :    
234 :    
235 :     output tensor [3,3] out = [[0,0,0],[0,0,0],[0,0,0]];
236 :     update {
237 :     if(cat==1){
238 :     print ("\n\n_------top");
239 :     print ("\n\n\nbase 1");
240 :     print("\nbase1-z0",base1_z0,xsq1_z0,ysq1_z0,diag1_z0);
241 :     print("\nbase1-z1",base1_z1,xsq1_z1,ysq1_z1,diag1_z1);
242 :     print("\nbase1-z2",base1_z2,xsq1_z2,ysq1_z2,diag1_z2);
243 :     print ("\n\n\ndiag2-z0");
244 :     print("\ndiag2_z0",diag2_z0);
245 :     print ("\n\n\n diag3_z0");
246 :     print("diag3_z0",diag3_z0);
247 :    
248 :     }
249 :    
250 :     /* see comment in ../fs2d/fs2d-scl.diderot about the need for these
251 :     conditionals around the print statements */
252 :     if (0 == idx0 && 0 == idx1 && 0 == idx2) {
253 :     print("NRRD0004\n");
254 :     print("# Complete NRRD file format specification at:\n");
255 :     print("# http://teem.sourceforge.net/nrrd/format.html\n");
256 :     /* NOTE: this assumes we haven't been compiled with --double,
257 :     and there isn't currently a way for the program to learn this
258 :     (which in our experience has not been a problem) */
259 :     print("type: float\n");
260 :     print("dimension: 4\n");
261 :     print("sizes: 9 ", sz0, " ", sz1, " ", sz2, "\n");
262 :     print("kinds: 3D-matrix space space space\n");
263 :     // print("centers: none cell cell cell\n");
264 :     // NOTE: this assumes machine endianness
265 :     print("endian: little\n");
266 :     print("encoding: raw\n");
267 :     print("space dimension: 3\n");
268 :     print("space directions: none (", edge0[0], ",", edge0[1], ",", edge0[2],
269 :     ") (", edge1[0], ",", edge1[1], ",", edge1[2],
270 :     ") (", edge2[0], ",", edge2[1], ",", edge2[2], ")\n");
271 :     print("space origin: (", orig[0], ",", orig[1], ",", orig[2], ")\n");
272 :     print("data file: out.nrrd\n");
273 :     print("byte skip: -1\n");
274 :     }
275 : cchiw 4385 // out = func(orig + idx0*edge0 + idx1*edge1 + idx2*edge2);
276 :     vec3 pos = orig + idx0*edge0 + idx1*edge1 + idx2*edge2;
277 :     real x = pos[0];
278 :     real y = pos[1];
279 :     real z = pos[2];
280 :     real axis1_z0 = cvt(x,y,1,base1_z0,xsq1_z0,ysq1_z0,diag1_z0);
281 :     real axis1_z1 = cvt(x,y,z,base1_z1,xsq1_z1,ysq1_z1,diag1_z1);
282 :     real axis1_z2 = cvt(x,y,z*z,base1_z2,xsq1_z2,ysq1_z2,diag1_z2);
283 :     real axis1_val = axis1_z0+axis1_z1+axis1_z2;
284 : cchiw 4236
285 : cchiw 4385 real axis2_z0 = cvt(x,y,1,base2_z0,xsq2_z0,ysq2_z0,diag2_z0);
286 :     real axis2_z1 = cvt(x,y,z,base2_z1,xsq2_z1,ysq2_z1,diag2_z1);
287 :     real axis2_z2 = cvt(x,y,z*z,base2_z2,xsq2_z2,ysq2_z2,diag2_z2);
288 :     real axis2_val = axis2_z0+axis2_z1+axis2_z2;
289 :    
290 :     real axis3_z0 = cvt(x,y,1,base3_z0,xsq3_z0,ysq3_z0,diag3_z0);
291 :     real axis3_z1 = cvt(x,y,z,base3_z1,xsq3_z1,ysq3_z1,diag3_z1);
292 :     real axis3_z2 = cvt(x,y,z*z,base3_z2,xsq3_z2,ysq3_z2,diag3_z2);
293 :     real axis3_val = axis3_z0+axis3_z1+axis3_z2;
294 :    
295 :    
296 :     real axis4_z0 = cvt(x,y,1,base4_z0,xsq4_z0,ysq4_z0,diag4_z0);
297 :     real axis4_z1 = cvt(x,y,z,base4_z1,xsq4_z1,ysq4_z1,diag4_z1);
298 :     real axis4_z2 = cvt(x,y,z*z,base4_z2,xsq4_z2,ysq4_z2,diag4_z2);
299 :     real axis4_val = axis4_z0+axis4_z1+axis4_z2;
300 :    
301 :     real axis5_z0 = cvt(x,y,1,base5_z0,xsq5_z0,ysq5_z0,diag5_z0);
302 :     real axis5_z1 = cvt(x,y,z,base5_z1,xsq5_z1,ysq5_z1,diag5_z1);
303 :     real axis5_z2 = cvt(x,y,z*z,base5_z2,xsq5_z2,ysq5_z2,diag5_z2);
304 :     real axis5_val = axis5_z0+axis5_z1+axis5_z2;
305 :    
306 :     real axis6_z0 = cvt(x,y,1,base6_z0,xsq6_z0,ysq6_z0,diag6_z0);
307 :     real axis6_z1 = cvt(x,y,z,base6_z1,xsq6_z1,ysq6_z1,diag6_z1);
308 :     real axis6_z2 = cvt(x,y,z*z,base6_z2,xsq6_z2,ysq6_z2,diag6_z2);
309 :     real axis6_val = axis6_z0+axis6_z1+axis6_z2;
310 :    
311 :     real axis7_z0 = cvt(x,y,1,base7_z0,xsq7_z0,ysq7_z0,diag7_z0);
312 :     real axis7_z1 = cvt(x,y,z,base7_z1,xsq7_z1,ysq7_z1,diag7_z1);
313 :     real axis7_z2 = cvt(x,y,z*z,base7_z2,xsq7_z2,ysq7_z2,diag7_z2);
314 :     real axis7_val = axis7_z0+axis7_z1+axis7_z2;
315 :    
316 :     real axis8_z0 = cvt(x,y,1,base8_z0,xsq8_z0,ysq8_z0,diag8_z0);
317 :     real axis8_z1 = cvt(x,y,z,base8_z1,xsq8_z1,ysq8_z1,diag8_z1);
318 :     real axis8_z2 = cvt(x,y,z*z,base8_z2,xsq8_z2,ysq8_z2,diag8_z2);
319 :     real axis8_val = axis8_z0+axis8_z1+axis8_z2;
320 :    
321 :     real axis9_z0 = cvt(x,y,1,base9_z0,xsq9_z0,ysq9_z0,diag9_z0);
322 :     real axis9_z1 = cvt(x,y,z,base9_z1,xsq9_z1,ysq9_z1,diag9_z1);
323 :     real axis9_z2 = cvt(x,y,z*z,base9_z2,xsq9_z2,ysq9_z2,diag9_z2);
324 :     real axis9_val = axis9_z0+axis9_z1+axis9_z2;
325 :    
326 :    
327 :     out = [[axis1_val, axis2_val, axis3_val],[axis4_val, axis5_val, axis6_val],[axis7_val, axis8_val, axis9_val]];
328 :    
329 : cchiw 4236 stabilize;
330 :     }
331 :     }
332 :     initially [ sample(idx0, idx1, idx2)
333 :     | idx2 in 0..(sz2-1), // slowest axis
334 :     idx1 in 0..(sz1-1), // medium axis
335 :     idx0 in 0..(sz0-1)]; // fastest axis

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0