Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] View of /tests/examples/vimg/vimg.diderot
ViewVC logotype

View of /tests/examples/vimg/vimg.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4640 - (download) (annotate)
Tue Sep 27 20:54:47 2016 UTC (2 years, 10 months ago) by glk
File size: 9638 byte(s)
initial result of svn export --username anonsvn --password=anonsvn https://svn.smlnj-gforge.cs.uchicago.edu/svn/diderot/branches/vis15/src/tests/
/* ==========================================
## 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) ];

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0