# SCM Repository

# View of /branches/vis15/src/tests/examples/vimg/vimg.diderot

Parent Directory | Revision Log

Revision

File size: 9617 byte(s)

**4178**- (**download**) (**annotate**)*Fri Jul 8 18:47:11 2016 UTC*(3 years, 1 month ago) by*jhr*File size: 9617 byte(s)

updating examples in merge to latest put from github

/* ========================================== ## 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("img.nrrd"); input image(1)[3] cmimg ("colormap to use") = image("cmap.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) ];

root@smlnj-gforge.cs.uchicago.edu | ViewVC Help |

Powered by ViewVC 1.0.0 |