input image(2)[] img ("2D image dataset to view") = image("test.nrrd");
input image(1)[3] cmimg ("colormap to use") = image("cmap.nrrd");
input vec2 cent ("center of viewing window") = [290,414];
input real fov ("height of viewing window") = 20000;
input int ss ("# of samples") = 500;
input real angle ("orientation (in counter-clockwise degrees) of viewing window") = 0;
input vec3 bkgd ("RGB color to show outside image domain") = [0.2,0.1,0.1];
input vec3 fcol ("RGB color to show inside features") = [0.0,1.0,0.0];
input real iso ("value at which to show approximate isocontour") = 0;
input real th ("apparent thickness of feature (meaning varies depending on which)") = 0.1;
input real cmin ("value to use at min end of colormap") = -1;
input real cmax ("value to use at max end of colormap") = 1;
/*
The field to sample and view; needs to be more or less differentiable
depending on what we want to do with it.
C1 kernels include: c1tent, ctmr, bspln3, c4hexic
C2 kernels include: c2ctmr, bspln3, c4hexic
*/
field#2(2)[] F = c4hexic ⊛ mirror(img);
/*
The colormap can always use tent; it is not differentiated. Note the
clamp() here: it means that before evaluating cmap(x), x will be clamped to
inside the domain of cmap, as defined by the location of the samples in the
cmimg, and the support of the reconstruction filter tent. From a
functional sense, it might make more sense for the notation to use
composition: "tent⊛cming◦clamp": clamp first, then lookup into cmimg and
reconstruct with tent, but this looks a bit weird.
*/
field#0(1)[3] cmap = tent ⊛ clamp(cmimg);
// compute the image-to-world transform based on specification
// of the viewing window
real wdth = fov;
real phi = angle*π/180;
tensor[2,2] rot = [[cos(phi),sin(phi)],[-sin(phi),cos(phi)]];
vec2 spc = [wdth/(ss-1), fov/(ss-1)];
vec2 dir0 = rot•[spc[0], 0];
vec2 dir1 = rot•[0, spc[1]];
vec2 orig = cent - (dir0*(ss-1) + dir1*(ss-1))/2;
// A function that is 1 at 0, and goes down linearly to 0 according to
// thickness parameter th. Note that if a function is a simple expression like
// this it can be defined with a single `=`, with no need for a `return`
// inside a sequence of statements. You might want to call this function
// "tent" but that name is already taken by a kernel.
function real bump(real z) = max(0, 1 - |z|/th);
strand sample(int ui, int vi) {
// The output type for a Diderot program is fixed at compile-time,
// regardless of input variable values or execution outcomes. So
// this program sometimes saves what is really just grayscale
// information in an RGB color image.
output vec3 rgb = bkgd;
real ll = 0;
update {
vec2 x = orig + ui*dir0 + vi*dir1;
ll = bump(F(x) - iso);
rgb = cmap(lerp(0,1,cmin,F(x),cmax));
rgb = lerp(rgb, fcol, ll);
stabilize;
}
}
initially [ sample(ui, vi) | vi in 0..(ss-1),
ui in 0..(ss-1) ];