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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4456 - (view) (download)

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

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