/* ==========================================
## vimg.diderot: 2D image sampler/viewer
This program needs a dataset to render, and a colormap, for example:
ln -s ../data/sscand.nrrd img.nrrd
ln -s ../cmap/spiral.nrrd cmap.nrrd
Then compile this program; assuming the directions at
https://github.com/Diderot-Language/examples you can:
../../vis12/bin/diderotc --exec vimg.diderot
If the needed `img.nrrd` file is missing, the error message looks something like:
uncaught exception Fail [Fail: Nrrd file "img.nrrd" does not exist]
raised at common/phase-timer.sml:78.50-78.52
raised at common/phase-timer.sml:78.50-78.52
raised at nrrd/nrrd-info.sml:146.15-146.74
in which case you should run the `ln -s` command above, or link `img.nrrd`
to some other 2D scalar nrrd file to view. The same applies to the need
for `cmap.nrrd` to link to a colormap.
The `-which` option will determine which function is sampled; look
for `(0 == which)` below to see the start of the function definitions.
Assuming the symbolic links given above (of `img.nrrd`
and `cmap.nrrd` to `../data/sscand.nrrd` and `../cmap/spiral.nrrd`
respectively), some examples usages are:
* Grayscale image of reconstructed field:
`./vimg -cent 290 414 -fov 45 -which 0`
* Grayscale image of gradient magnitude:
`./vimg -cent 290 414 -fov 45 -which 1`
* Color image of gradient vector (blue is always zero):
`./vimg -cent 290 414 -fov 45 -which 2`
* Colormapped field, with naive pseudo-isocontour:
`./vimg -cent 290 414 -fov 45 -which 3 -cmin -500 -cmax 1900 -iso 1210 -th 20`
* Colormapped field, with smarter pseudo-isocontour:
`./vimg -cent 290 414 -fov 45 -which 4 -cmin -500 -cmax 1900 -iso 1210 -th 0.2`
* Colormapped field, with ridge lines:
`./vimg -cent 290 414 -fov 45 -which 5 -cmin -500 -cmax 1900 -th 0.2 -sthr 2`
* Colormapped field, with valley lines:
`./vimg -cent 290 414 -fov 45 -which 6 -cmin -500 -cmax 1900 -th 0.2 -sthr 2`
* Colormapped field, with blue maxima:
`./vimg -cent 290 414 -fov 45 -which 7 -cmin -500 -cmax 1900 -th 0.25 -sthr 25 -fcol 0 0 0.8`
* Colormapped field, with green saddle points:
`./vimg -cent 290 414 -fov 45 -which 8 -cmin -500 -cmax 1900 -th 0.25 -fcol 0 1 0`
* Colormapped field, with cyan minima:
`./vimg -cent 290 414 -fov 45 -which 9 -cmin -500 -cmax 1900 -th 0.25 -sthr 25 -fcol 0 1 1`
Each command can be followed by `unu quantize -b 8 -i rgb.nrrd -o rgb.png` to create
an 8-bit image version of the output. The `-which 3` and `-which 4` commands
show an important comparison, demonstrating how knowing the gradient permits
drawing of equal-thickness isocontours (according to the first-order Taylor
expansion).
Viewing `../data/sscand.nrrd` with the parameters above gives a roughly
100km view of the area around Geilo, Norway, site of the
[Winter School](http://www.sintef.no/projectweb/geilowinterschool/2016-scientific-visualization/)
for which this program was originally created.
========================================== */
input image(2)[] img ("2D image dataset to view") = image("../data/sscand.nrrd");
input image(1)[3] cmimg ("colormap to use") = image("../cmap/spiral.nrrd");
input vec2 cent ("center of viewing window") = [0,0];
input real fov ("height of viewing window") = 2;
input int sx ("# of horizontal samples in viewing window") = 640;
input int sy ("# of vertical samples in viewing window") = 480;
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,0.0,0.0];
input int which ("what to show about the image, currently from 0 to 9") = 0;
input real iso ("value at which to show approximate isocontour") = 0;
input real th ("apparent thickness of feature (meaning varies depending on which)") = 1;
input real cmin ("value to use at min end of colormap") = 0;
input real cmax ("value to use at max end of colormap") = 1;
input real seps ("epsilon >= 0 value to use in divisor when computing strength") = 0.001;
input real sthr ("strength threshold for showing ridges and valleys") = 0;
/*
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 ⊛ 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*sx/sy;
real phi = angle*π/180;
tensor[2,2] rot = [[cos(phi),sin(phi)],[-sin(phi),cos(phi)]];
vec2 spc = [wdth/(sx-1), fov/(sy-1)];
vec2 dir0 = rot•[spc[0], 0];
vec2 dir1 = rot•[0, spc[1]];
vec2 orig = cent - (dir0*(sx-1) + dir1*(sy-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);
// 2x2 matrix inversion should be part of the language; until
// then it is pretty easy to create a function for it.
function tensor[2,2] inv(tensor[2,2] m) {
real d = m[0,0]*m[1,1] - m[0,1]*m[1,0];
return [[m[1,1],-m[0,1]],[-m[1,0],m[0,0]]]/d;
}
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;
if (inside(x, F)) {
if (0 == which) {
// the image value itself
real y = F(x);
rgb = [y,y,y];
} else if (1 == which) {
// the gradient magnitude
real y = |∇F(x)|;
rgb = [y,y,y];
} else if (2 == which) {
// the gradient vector, saved in
// first two components of output
vec2 gg = ∇F(x);
rgb = [gg[0],gg[1],0];
} else if (3 == which) {
/* Colormapped, with simple value-based isocontour. Note that
there isn't (currently) a way to learn within the language what
the domain of colormap (or of any field) is, so rescaling the
function range to [0,1] here reflects either prior knowledge
or an assumption about the orientation meta-data in whatever
cmap.nrrd links to. There is no need to clamp before evaluating
cmap(); that is handled by the clamp in cmap's declaration */
rgb = cmap(lerp(0,1,cmin,F(x),cmax));
ll = bump(F(x) - iso);
} else if (4 == which) {
// colormapped, with Newton-based isocontour
rgb = cmap(lerp(0,1,cmin,F(x),cmax));
ll = bump((F(x) - iso)/|∇F(x)|);
} else if (5 == which) {
// colormapped, with Newton-based ridge line
rgb = cmap(lerp(0,1,cmin,F(x),cmax));
vec2{2} evec = evecs(∇⊗∇F(x));
real{2} eval = evals(∇⊗∇F(x));
real strn = -eval{1}/(|∇F(x)| + seps);
if (eval{1} < 0 && strn > sthr) {
vec2 delta = (evec{1}⊗evec{1}/eval{1})•∇F(x);
ll = bump(|delta|);
}
} else if (6 == which) {
// colormapped, with Newton-based valley line
rgb = cmap(lerp(0,1,cmin,F(x),cmax));
vec2{2} evec = evecs(∇⊗∇F(x));
real{2} eval = evals(∇⊗∇F(x));
real strn = eval{0}/(|∇F(x)| + seps);
if (eval{0} > 0 && strn > sthr) {
vec2 delta = (evec{0}⊗evec{0}/eval{0})•∇F(x);
ll = bump(|delta|);
}
} else if (7 == which || 8 == which || 9 == which) {
// colormapped, with critical points
rgb = cmap(lerp(0,1,cmin,F(x),cmax));
vec2 delta = inv(∇⊗∇F(x))•∇F(x);
real{2} eval = evals(∇⊗∇F(x));
if (7 == which) { // maxima
if (eval{0} < 0 && eval{1} < 0 && eval{0} < -sthr) {
ll = bump(|delta|);
}
} else if (8 == which) { // saddle
if (eval{0} > 0 && eval{1} < 0) {
ll = bump(|delta|);
}
} else { // 9 == which: minima
if (eval{0} > 0 && eval{1} > 0 && eval{1} > sthr) {
ll = bump(|delta|);
}
}
} else {
// update "input int which" annotation above as cases are added
print("Sorry, no function defined for which = ", which, "\n");
}
rgb = lerp(rgb, fcol, ll);
}
stabilize;
}
}
initially [ sample(ui, vi) | vi in 0..(sy-1),
ui in 0..(sx-1) ];