-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathfs2d-scl.diderot
162 lines (149 loc) · 7.19 KB
/
fs2d-scl.diderot
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
#version 1.0
/* ==========================================
## fs2d-scl.diderot: function sampler, 2D scalar synthetic fields
This programs generates synthetic scalar 2D data on regular grids, and
the grids have specific location and orientation in world space.
The input arguments here make it easy to sample the same underlying
function on grids with differing orientation, resolution, and shear.
By using one strand per function evaluation, this is not an efficient
use of Diderot, but it can be nice to have a controlled place to test
the evaluation of Diderot expressions, especially since it prints out
the NRRD header that contains all the orientation meta-data to locate
the sampling grid in a world-space.
The `-which` option will determine which function is sampled; look
for `("x" == which)` below to see the start of the function definitions.
This program is unusual in that its printed output needs to be captured
in order to have a NRRD header that records the orientation of the
sampling grid, so using the program involves redirection. To
get a self-contained parab.nrrd containing a parabola function
#R
./fs2d-scl -which parab | unu save -f nrrd -o parab.nrrd
rm out.nrrd
As noted in the comments, the NRRD header generated by this program assumes:
1. This program was not compiled with `--double`
2. The program is running on a little-endian machine.
3. The program's `-o` option was not used to change the output filename.
#!
#_ # (here to ensure this program can be compiled, and to be
#_ # available to other tests)
#=diderotc
========================================== */
/* When an input global variable is not initialized, the command-line
interface has no default value for the variable, so it becomes a
required option instead of an optional option. In this case,
the "-which" flag must be used to choose a function to sample */
input string which ("which function to sample, currently one of: x, y, parab, tcubic, sin, circ, scparab, cos, foo, phi");
input int size0 ("# samples on faster axis") = 101;
input int size1 ("# samples on slower axis") = 100;
input real width ("approx width of world-space region sampled") = 1;
input vec2 off ("translation offset from origin-centered grid") = [0,0];
input real shear ("amount of shear in sampling grid") = 0;
input real angle ("orientation (in degrees) of faster axis") = 0;
input real wangle ("orientation (in degrees) of world coords") = 0;
input vec4 parm ("parameters that functions may use") = [0,0,0,0];
real phi = angle*π/180; // π means what you think
// rotation by phi of [1,0] towards [0,1]
tensor[2,2] rot = [[cos(phi),-sin(phi)],[sin(phi),cos(phi)]];
// repeating for world coords
real wphi = wangle*π/180;
tensor[2,2] wrot = [[cos(wphi),-sin(wphi)],[sin(wphi),cos(wphi)]];
// sample spacing on faster and (unsheared) slower axis
vec2 spc = [width/(size0-1), width/(size1-1)];
// inter-sample vector on faster axis
vec2 edge0 = rot•[spc[0], 0];
// inter-sample vector on slower axis
vec2 edge1 = rot•[0, spc[1]] + shear*spc[1]*normalize(edge0);
// location of first sample
vec2 orig = -(edge0*(size0-1) + edge1*(size1-1))/2 + off;
// univariate normal distribution
function real gauss(real stdv, real x) =
exp(-(x*x)/(2*stdv^2))/(stdv*sqrt(2*π));
/* The function to evaluate at each grid point. The parm[] vector is just a
tuple of numbers, rather than a geometric vector, and is made somewhat
confusing by the intent that the function be useful with the default
all-zero parm */
function real func(vec2 pos0) {
vec2 pos = wrot•pos0;
real xx = pos[0];
real yy = pos[1];
real tt = atan2(yy,xx);
real rr = |pos|;
real ret = 0;
/* parenthesized numbers in comments are the integers
previously used to identify the case */
if ("x" == which) { // x ramp (formerly 0)
ret = xx;
} else if ("y" == which) { // y ramp (1)
ret = yy;
} else if ("parab" == which) { // parabola (2)
ret = pos•pos;
} else if ("tcubic" == which) { // Tschirnhausen cubic (3)
// using 4*y^2 instead of y^2 for aesthetic reasons
// expects width of about 8
ret = xx^3 + 3*xx^2 - 4*yy^2;
} else if ("sin" == which) { // sine waves (4)
ret = sin((1 + parm[0])*(xx + parm[2]))
+ sin((1 + parm[1])*(yy + parm[3]));
} else if ("circ" == which) { // circle with soft edge, plus ramps (5)
// expects width of about 4
ret = lerp(1, 0, -1, erf((rr - (1 + parm[0]))*(8*(1+parm[1]))), 1)
+ parm[2]*xx + parm[3]*yy;
} else if ("scparab" == which) { // scaled parabola (6)
ret = xx*xx*(1 + parm[0]) + yy*yy*(1 + parm[1]);
} else if ("cos" == which) { // cosine waves (7)
ret = cos((1 + parm[0])*|pos|)*(1 + parm[1]) + parm[2];
} else if ("foo" == which) { // something w/ convoluted zero-levelset (8)
ret = 0.5 + rr*(sin(4.5*tt)^2 - 7*rr)
- 0.18*gauss(0.04,rr-0.18)*cos(4.5*tt)^2
+ 0.11*gauss(0.08,rr-0.35)*sin(4.5*tt)^2;
} else if ("phi" == which) { // angle around origin, from atan2 (9)
ret = tt;
} else {
// update "input string which" annotation above as cases are added
print("Sorry, no function defined for which = ", which, "\n");
}
return ret;
}
strand sample(int idx0, int idx1) {
output real out = 0.0;
update {
/* This version of the Diderot syntax doesn't allow print statements
from global initialization, so to print something once per
program, you need to test for a condition that will be true
for one strand. Due to the immediate stabilize, this will only
run for one iteration. */
if (0 == idx0 && 0 == idx1) {
print("NRRD0004\n");
print("# Complete NRRD file format specification at:\n");
print("# http://teem.sourceforge.net/nrrd/format.html\n");
// NOTE: this assumes we haven't been compiled with --double
print("type: float\n");
print("dimension: 2\n");
print("sizes: ", size0, " ", size1, "\n");
print("kinds: space space\n");
print("centers: cell cell\n");
// NOTE: this assumes machine endianness
print("endian: little\n");
print("encoding: raw\n");
print("space dimension: 2\n");
// Diderot prints vectors like it parses them, e.g "[0.1,0.3]"
// but this is not how orientation vectors are stored in the NRRD
// header, hence the need to print individual components
print("space directions: (", edge0[0], ",", edge0[1],
") (", edge1[0], ",", edge1[1], ")\n");
print("space origin: (", orig[0], ",", orig[1], ")\n");
// NOTE: this assumes output filename is not explicitly set
print("data file: out.nrrd\n");
print("byte skip: -1\n");
}
out = func(orig + idx0*edge0 + idx1*edge1);
stabilize;
}
}
/* Create one strand per sample point. The "initially []" creates an array of
strands, so the output is saved as a 2D array with axis length size0 and
size1. Had the program been created as a collection of strands (with
"initially {}"), the output would be a 1D array length size0*size1. */
initially [ sample(idx0, idx1)
| idx1 in 0..(size1-1), // SLOWER axis
idx0 in 0..(size0-1) ]; // FASTER axis