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 4236 - (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 int which ("which function to sample (from 0 to 2)") = 0;
35 :     input vec3 parm0 ("parameters that functions may use") = [0,0,0];
36 :     input vec3 parm1 ("more parameters that functions may use") = [0,0,0];
37 :     input vec3 parm2 ("even more parameters that functions may use") = [0,0,0];
38 :    
39 :     real theta = angle*π/180;
40 :     // unit-length axis of rotation scaled by sin(theta/2)
41 :     vec3 snax = sin(theta/2)*normalize(axis);
42 :     // quaternion representing rotation
43 :     vec4 qq = [cos(theta/2), snax[0], snax[1], snax[2]];
44 :     // create rotation matrix by converstion from quaternion
45 :     // (which aren't currently natively supported in Diderot)
46 :     tensor[3,3] rot = [[qq[0]*qq[0] + qq[1]*qq[1] - qq[2]*qq[2] - qq[3]*qq[3],
47 :     2*(qq[1]*qq[2] - qq[0]*qq[3]),
48 :     2*(qq[1]*qq[3] + qq[0]*qq[2])],
49 :     [2*(qq[1]*qq[2] + qq[0]*qq[3]),
50 :     qq[0]*qq[0] - qq[1]*qq[1] + qq[2]*qq[2] - qq[3]*qq[3],
51 :     2*(qq[2]*qq[3] - qq[0]*qq[1])],
52 :     [2*(qq[1]*qq[3] - qq[0]*qq[2]),
53 :     2*(qq[2]*qq[3] + qq[0]*qq[1]),
54 :     qq[0]*qq[0] - qq[1]*qq[1] - qq[2]*qq[2] + qq[3]*qq[3]]];
55 :     // per axis inter-sample edge vectors
56 :     real space0 = width/(sz0-1);
57 :     real space1 = width/(sz1-1);
58 :     real space2 = width/(sz2-1);
59 :     vec3 edge0 = rot•[space0, 0, 0];
60 :     vec3 edge1 = rot•[0, space1, 0];
61 :     vec3 edge2 = rot•[0, 0, space2];
62 :     // offset in world space
63 :     vec3 offws = [off[0]*space0, off[1]*space1, off[2]*space2];
64 :     // location of first sample point
65 :     vec3 orig = -(edge0*(sz0-1) + edge1*(sz1-1) + edge2*(sz2-1))/2 + offws;
66 :    
67 :     //-- field defined by coefficients --
68 :     // z=1
69 :     input vec2 base1_z0 = [0, 0]; // b*x, c*y
70 :     input vec2 xsq1_z0 = [0, 0]; // d*x^2, g*y*x^2
71 :     input vec2 ysq1_z0 = [0, 0]; // f*y^2, h*x*y^2
72 :     input vec3 diag1_z0 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
73 :     // z=z
74 :     input vec2 base1_z1 = [0, 0]; // b*x, c*y
75 :     input vec2 xsq1_z1 = [0, 0]; // d*x^2, g*y*x^2
76 :     input vec2 ysq1_z1 = [0, 0]; // f*y^2, h*x*y^2
77 :     input vec3 diag1_z1 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
78 :     // z=z*z
79 :     input vec2 base1_z2 = [0, 0]; // b*x, c*y
80 :     input vec2 xsq1_z2 = [0, 0]; // d*x^2, g*y*x^2
81 :     input vec2 ysq1_z2 = [0, 0]; // f*y^2, h*x*y^2
82 :     input vec3 diag1_z2 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
83 :    
84 :     //-- field defined by coefficients --
85 :     // z=1
86 :     input vec2 base2_z0 = [0, 0]; // b*x, c*y
87 :     input vec2 xsq2_z0 = [0, 0]; // d*x^2, g*y*x^2
88 :     input vec2 ysq2_z0 = [0, 0]; // f*y^2, h*x*y^2
89 :     input vec3 diag2_z0 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
90 :     // z=z
91 :     input vec2 base2_z1 = [0, 0]; // b*x, c*y
92 :     input vec2 xsq2_z1 = [0, 0]; // d*x^2, g*y*x^2
93 :     input vec2 ysq2_z1 = [0, 0]; // f*y^2, h*x*y^2
94 :     input vec3 diag2_z1 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
95 :     // z=z*z
96 :     input vec2 base2_z2 = [0, 0]; // b*x, c*y
97 :     input vec2 xsq2_z2 = [0, 0]; // d*x^2, g*y*x^2
98 :     input vec2 ysq2_z2 = [0, 0]; // f*y^2, h*x*y^2
99 :     input vec3 diag2_z2 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
100 :    
101 :     //-- field defined by coefficients --
102 :     // z=1
103 :     input vec2 base3_z0 = [0, 0]; // b*x, c*y
104 :     input vec2 xsq3_z0 = [0, 0]; // d*x^2, g*y*x^2
105 :     input vec2 ysq3_z0 = [0, 0]; // f*y^2, h*x*y^2
106 :     input vec3 diag3_z0 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
107 :     // z=z
108 :     input vec2 base3_z1 = [0, 0]; // b*x, c*y
109 :     input vec2 xsq3_z1 = [0, 0]; // d*x^2, g*y*x^2
110 :     input vec2 ysq3_z1 = [0, 0]; // f*y^2, h*x*y^2
111 :     input vec3 diag3_z1 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
112 :     // z=z*z
113 :     input vec2 base3_z2 = [0, 0]; // b*x, c*y
114 :     input vec2 xsq3_z2 = [0, 0]; // d*x^2, g*y*x^2
115 :     input vec2 ysq3_z2 = [0, 0]; // f*y^2, h*x*y^2
116 :     input vec3 diag3_z2 = [0, 0, 0]; // a, e*x*y, i*y^2*x^2
117 :     input int cat=0;
118 :    
119 :    
120 :    
121 :     function real cvt(real x, real y, real z, vec2 base, vec2 xsq, vec2 ysq, vec3 diag){
122 :     real b = base[0];
123 :     real c = base[1];
124 :     real d = xsq[0];
125 :     real g = xsq[1];
126 :     real f = ysq[0];
127 :     real h = ysq[1];
128 :     real a = diag[0];
129 :     real e = diag[1];
130 :     real i = diag[2];
131 :     real t0 = b*x + c*y;
132 :     real t1 = d*x^2 + g*y*x^2;
133 :     real t2 = f*y^2 + h*x*y^2;
134 :     real t3 = a + e*x*y + i*(y^2)*(x^2);
135 :     return z*(t0+t1+t2+t3);
136 :     }
137 :    
138 :    
139 :     function tensor[3,3] func(vec3 pos) {
140 :     real x = pos[0];
141 :     real y = pos[1];
142 :     real z = pos[2];
143 :     tensor[3,3] ret = [[0,0,0],[0,0,0],[0,0,0]];
144 :     if (0 == which) {
145 :     /*
146 :     real axis1_z0 = cvt(x,y,1,base1_z0,xsq1_z0,ysq1_z0,diag1_z0);
147 :     real axis1_z1 = cvt(x,y,z,base1_z1,xsq1_z1,ysq1_z1,diag1_z1);
148 :     real axis1_z2 = cvt(x,y,z*z,base1_z2,xsq1_z2,ysq1_z2,diag1_z2);
149 :     real axis1_val = axis1_z0+axis1_z1+axis1_z2;
150 :     if(cat==1){
151 :     print(" \naxis1_val", axis1_val);
152 :     }
153 :     real axis2_z0 = cvt(x,y,1,base2_z0,xsq2_z0,ysq2_z0,diag2_z0);
154 :     real axis2_z1 = cvt(x,y,z,base2_z1,xsq2_z1,ysq2_z1,diag2_z1);
155 :     real axis2_z2 = cvt(x,y,z*z,base2_z2,xsq2_z2,ysq2_z2,diag2_z2);
156 :     real axis2_val = axis2_z0+axis2_z1+axis2_z2;
157 :     if(cat==1){
158 :     print(" \naxis2_val", axis2_val);
159 :     }
160 :     real axis3_z0 = cvt(x,y,1,base3_z0,xsq3_z0,ysq3_z0,diag3_z0);
161 :     real axis3_z1 = cvt(x,y,z,base3_z1,xsq3_z1,ysq3_z1,diag3_z1);
162 :     real axis3_z2 = cvt(x,y,z*z,base3_z2,xsq3_z2,ysq3_z2,diag3_z2);
163 :     real axis3_val = axis3_z0+axis3_z1+axis3_z2;
164 :     if(cat==1){
165 :     print(" \naxis3_val", axis3_val);
166 :     }
167 :     // vec3 ret1 = [axis1_val, axis2_val, axis3_val];
168 :     vec3 ret1 = [1,4,2];
169 :     if(cat==1){
170 :     print(" \n ret1", ret1);
171 :     }
172 :     //ret = [ret1, ret1, ret1];
173 :     */
174 :     //fixed return
175 :     ret = [[0,1,2],[4,5,6],[8,9,10]];
176 :     //ret = [[7,1,2],[1,5,2],[8,9,10]];
177 :     } else {
178 :     // update "input int which" annotation above as cases are added
179 :     print("Sorry, no function defined for which = ", which, "\n");
180 :     }
181 :     return ret;
182 :     }
183 :    
184 :     strand sample(int idx0, int idx1, int idx2) {
185 :    
186 :    
187 :    
188 :     output tensor [3,3] out = [[0,0,0],[0,0,0],[0,0,0]];
189 :     update {
190 :     if(cat==1){
191 :     print ("\n\n_------top");
192 :     print ("\n\n\nbase 1");
193 :     print("\nbase1-z0",base1_z0,xsq1_z0,ysq1_z0,diag1_z0);
194 :     print("\nbase1-z1",base1_z1,xsq1_z1,ysq1_z1,diag1_z1);
195 :     print("\nbase1-z2",base1_z2,xsq1_z2,ysq1_z2,diag1_z2);
196 :     print ("\n\n\ndiag2-z0");
197 :     print("\ndiag2_z0",diag2_z0);
198 :     print ("\n\n\n diag3_z0");
199 :     print("diag3_z0",diag3_z0);
200 :    
201 :     }
202 :    
203 :     /* see comment in ../fs2d/fs2d-scl.diderot about the need for these
204 :     conditionals around the print statements */
205 :     if (0 == idx0 && 0 == idx1 && 0 == idx2) {
206 :     print("NRRD0004\n");
207 :     print("# Complete NRRD file format specification at:\n");
208 :     print("# http://teem.sourceforge.net/nrrd/format.html\n");
209 :     /* NOTE: this assumes we haven't been compiled with --double,
210 :     and there isn't currently a way for the program to learn this
211 :     (which in our experience has not been a problem) */
212 :     print("type: float\n");
213 :     print("dimension: 4\n");
214 :     print("sizes: 9 ", sz0, " ", sz1, " ", sz2, "\n");
215 :     print("kinds: 3D-matrix space space space\n");
216 :     // print("centers: none cell cell cell\n");
217 :     // NOTE: this assumes machine endianness
218 :     print("endian: little\n");
219 :     print("encoding: raw\n");
220 :     print("space dimension: 3\n");
221 :     print("space directions: none (", edge0[0], ",", edge0[1], ",", edge0[2],
222 :     ") (", edge1[0], ",", edge1[1], ",", edge1[2],
223 :     ") (", edge2[0], ",", edge2[1], ",", edge2[2], ")\n");
224 :     print("space origin: (", orig[0], ",", orig[1], ",", orig[2], ")\n");
225 :     print("data file: out.nrrd\n");
226 :     print("byte skip: -1\n");
227 :     }
228 :     out = func(orig + idx0*edge0 + idx1*edge1 + idx2*edge2);
229 :     if(cat==1){
230 :     print(" \n out0-", out[1,0],out[1,1],out[1,2]);
231 :     }
232 :    
233 :     stabilize;
234 :     }
235 :     }
236 :     initially [ sample(idx0, idx1, idx2)
237 :     | idx2 in 0..(sz2-1), // slowest axis
238 :     idx1 in 0..(sz1-1), // medium axis
239 :     idx0 in 0..(sz0-1)]; // fastest axis

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