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_fv3d3.diderot
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4456 - (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 : cchiw 4444
35 : cchiw 4236 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 : cchiw 4456 input tensor[4] setA1_z0 =[0,0,0,0];
69 :     input tensor[4] setB1_z0 =[0,0,0,0];
70 :     input tensor[4] setC1_z0 =[0,0,0,0];
71 :     input tensor[4] setD1_z0 =[0,0,0,0];
72 : cchiw 4236
73 : cchiw 4456 input tensor[4] setA1_z1 =[0,0,0,0];
74 :     input tensor[4] setB1_z1 =[0,0,0,0];
75 :     input tensor[4] setC1_z1 =[0,0,0,0];
76 :     input tensor[4] setD1_z1 =[0,0,0,0];
77 : cchiw 4236
78 : cchiw 4456 input tensor[4] setA1_z2 =[0,0,0,0];
79 :     input tensor[4] setB1_z2 =[0,0,0,0];
80 :     input tensor[4] setC1_z2 =[0,0,0,0];
81 :     input tensor[4] setD1_z2 =[0,0,0,0];
82 : cchiw 4236
83 :    
84 : cchiw 4456 input tensor[4] setA2_z0 =[0,0,0,0];
85 :     input tensor[4] setB2_z0 =[0,0,0,0];
86 :     input tensor[4] setC2_z0 =[0,0,0,0];
87 :     input tensor[4] setD2_z0 =[0,0,0,0];
88 :     input tensor[4] setA2_z1 =[0,0,0,0];
89 :     input tensor[4] setB2_z1 =[0,0,0,0];
90 :     input tensor[4] setC2_z1 =[0,0,0,0];
91 :     input tensor[4] setD2_z1 =[0,0,0,0];
92 :     input tensor[4] setA2_z2 =[0,0,0,0];
93 :     input tensor[4] setB2_z2 =[0,0,0,0];
94 :     input tensor[4] setC2_z2 =[0,0,0,0];
95 :     input tensor[4] setD2_z2 =[0,0,0,0];
96 :    
97 :     input tensor[4] setA3_z0 =[0,0,0,0];
98 :     input tensor[4] setB3_z0 =[0,0,0,0];
99 :     input tensor[4] setC3_z0 =[0,0,0,0];
100 :     input tensor[4] setD3_z0 =[0,0,0,0];
101 :     input tensor[4] setA3_z1 =[0,0,0,0];
102 :     input tensor[4] setB3_z1 =[0,0,0,0];
103 :     input tensor[4] setC3_z1 =[0,0,0,0];
104 :     input tensor[4] setD3_z1 =[0,0,0,0];
105 :     input tensor[4] setA3_z2 =[0,0,0,0];
106 :     input tensor[4] setB3_z2 =[0,0,0,0];
107 :     input tensor[4] setC3_z2 =[0,0,0,0];
108 :     input tensor[4] setD3_z2 =[0,0,0,0];
109 :    
110 :     function real cvt(real x, real y, real z, tensor[4] setA, tensor[4] setB, tensor[4] setC, tensor[4] setD) {
111 :     real a = setA[0];
112 :     real b = setA[1];
113 :     real c = setA[2];
114 :     real d = setA[3];
115 :     real e = setB[0];
116 :     real f = setB[1];
117 :     real g = setB[2];
118 :     real h = setB[3];
119 :     real i = setC[0];
120 :     real j = setC[1];
121 :     real k = setC[2];
122 :     real l = setC[3];
123 :     real m = setD[0];
124 :     real n = setD[1];
125 :     real o = setD[2];
126 :     real p = setD[3];
127 :    
128 :     // as intended to be represented
129 :     real tA = a + b*y + c*x*y+ d*x;
130 :     real tB = ((e+f*x+g*(x*x))*y*y) + (h*(x*x)*y);
131 :     real tC = i*(x*x) + (j+k*x+l*(x*x))*y*y*y;
132 :     real tD = (x*x*x)*((m*y*y*y)+(n*y*y)+(o*y)+p);
133 :     return z*(tA+tB+tC+tD);
134 : cchiw 4236 }
135 :    
136 :    
137 : cchiw 4456 function tensor[3] func(vec3 pos) {
138 : cchiw 4444 real x = pos[0];
139 :     real y = pos[1];
140 :     real z = pos[2];
141 : cchiw 4456
142 :     real axis1_z0 = cvt(x,y,1, setA1_z0, setB1_z0, setC1_z0, setD1_z0);
143 :     real axis1_z1 = cvt(x,y,z, setA1_z1, setB1_z1, setC1_z1, setD1_z1);
144 :     real axis1_z2 = cvt(x,y,z*z,setA1_z2, setB1_z2, setC1_z2, setD1_z2);
145 :     real axis1= axis1_z0+axis1_z1+axis1_z2;
146 :    
147 :     real axis2_z0 = cvt(x,y,1, setA2_z0, setB2_z0, setC2_z0, setD2_z0);
148 :     real axis2_z1 = cvt(x,y,z, setA2_z1, setB2_z1, setC2_z1, setD2_z1);
149 :     real axis2_z2 = cvt(x,y,z*z,setA2_z2, setB2_z2, setC2_z2, setD2_z2);
150 :     real axis2= axis2_z0+axis2_z1+axis2_z2;
151 :    
152 :     real axis3_z0 = cvt(x,y,1, setA3_z0, setB3_z0, setC3_z0, setD3_z0);
153 :     real axis3_z1 = cvt(x,y,z, setA3_z1, setB3_z1, setC3_z1, setD3_z1);
154 :     real axis3_z2 = cvt(x,y,z*z,setA3_z2, setB3_z2, setC3_z2, setD3_z2);
155 :     real axis3= axis3_z0+axis3_z1+axis3_z2;
156 :     return [axis1, axis2, axis3];
157 :    
158 : cchiw 4236 }
159 :    
160 :     strand sample(int idx0, int idx1, int idx2) {
161 : cchiw 4456
162 :    
163 :     output tensor [3] out = [0, 0, 0];
164 : cchiw 4236 update {
165 : cchiw 4456
166 :    
167 : cchiw 4236 /* see comment in ../fs2d/fs2d-scl.diderot about the need for these
168 :     conditionals around the print statements */
169 :     if (0 == idx0 && 0 == idx1 && 0 == idx2) {
170 :     print("NRRD0004\n");
171 :     print("# Complete NRRD file format specification at:\n");
172 :     print("# http://teem.sourceforge.net/nrrd/format.html\n");
173 :     /* NOTE: this assumes we haven't been compiled with --double,
174 :     and there isn't currently a way for the program to learn this
175 :     (which in our experience has not been a problem) */
176 :     print("type: float\n");
177 :     print("dimension: 4\n");
178 :     print("sizes: 3 ", sz0, " ", sz1, " ", sz2, "\n");
179 :     print("kinds: 3-vector space space space\n");
180 : cchiw 4456 // print("centers: none cell cell cell\n");
181 : cchiw 4236 // NOTE: this assumes machine endianness
182 :     print("endian: little\n");
183 :     print("encoding: raw\n");
184 :     print("space dimension: 3\n");
185 :     print("space directions: none (", edge0[0], ",", edge0[1], ",", edge0[2],
186 :     ") (", edge1[0], ",", edge1[1], ",", edge1[2],
187 :     ") (", edge2[0], ",", edge2[1], ",", edge2[2], ")\n");
188 :     print("space origin: (", orig[0], ",", orig[1], ",", orig[2], ")\n");
189 :     print("data file: out.nrrd\n");
190 :     print("byte skip: -1\n");
191 :     }
192 :     out = func(orig + idx0*edge0 + idx1*edge1 + idx2*edge2);
193 :     stabilize;
194 :     }
195 :     }
196 :     initially [ sample(idx0, idx1, idx2)
197 :     | idx2 in 0..(sz2-1), // slowest axis
198 :     idx1 in 0..(sz1-1), // medium axis
199 :     idx0 in 0..(sz0-1)]; // fastest axis

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