/* ==========================================
## 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.
Assuming the directions at https://github.com/Diderot-Language/examples
this program can be compiled with:
../../vis12/bin/diderotc --exec fs2d-scl.diderot
The `-which` option will determine which function is sampled; look
for `(0 == 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
./fs2d-scl -which 2 | unu save -f nrrd -o parab.nrrd
rm out.nrrd
Note that 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.
========================================== */
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 int which ("which function to sample, currently from 0 to 5") = 0;
input vec4 parm ("parameters that functions may use") = [0,0,0,0];
real theta = angle*π/180;
// rotation by theta of [1,0] towards [0,1]
tensor[2,2] rot = [[cos(theta),-sin(theta)],[sin(theta),cos(theta)]];
// 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*edge0;
// location of first sample
vec2 orig = -(edge0*(size0-1) + edge1*(size1-1))/2 + off;
/*
The function to evaluate at each grid point The usage of the parm[]
vector elements is just as a tuple of numbers, rather than a geometric
vector, and is made somewhat confusing by the intent that the the
function be useful with the default all-zero parm
*/
function real func(vec2 pos) {
real x = pos[0];
real y = pos[1];
real ret = 0;
if (0 == which) { // 0: x ramp
ret = x;
} else if (1 == which) { // 1: y ramp
ret = y;
} else if (2 == which) { // 2: parabola
ret = pos•pos;
} else if (3 == which) { // 3: Tschirnhausen cubic
// using 4*y^2 instead of y^2 for aesthetic reasons
// expects width of about 8
ret = x^3 + 3*x^2 - 4*y^2;
} else if (4 == which) { // 4: sine waves
ret = sin((1 + parm[0])*(x + parm[2]))
+ sin((1 + parm[1])*(y + parm[3]));
} else if (5 == which) { // 5: circle with soft edge, plus ramps
real r = |pos|;
// expects width of about 4
ret = lerp(1, 0, -1, erf((r - (1 + parm[0]))*(8*(1+parm[1]))), 1)
+ parm[2]*x + parm[3]*y;
} else {
// update "input int 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 {
/* Diderot doesn't (currently) 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. By 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 a grid
** of strands (as opposed to a collection).
*/
initially [ sample(idx0, idx1)
| idx1 in 0..(size1-1), // SLOWER axis
idx0 in 0..(size0-1) ]; // FASTER axis