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

SCM Repository

[diderot] View of /branches/vis15/src/tests/examples/mip/mip.diderot
ViewVC logotype

View of /branches/vis15/src/tests/examples/mip/mip.diderot

Parent Directory Parent Directory | Revision Log Revision Log

Revision 3588 - (download) (annotate)
Thu Jan 14 15:25:22 2016 UTC (4 years, 2 months ago) by jhr
File size: 11717 byte(s)
adding tests to merge branch
/* ==========================================
## mip.diderot: basic maximum intensity projection (MIP) volume rendering

Maximum intensity projection is the minimalist volume visualization
tool. This implementation is about as short as possible; much of the code is
spent on specifying and setting up the ray geometry. It can also be adapted
(by changing a line or two) to do other kinds of projections.  The code for
camera geometry and ray sampling will be re-used verbatim in other volume
rendering programs. This may later become a target for evolving Diderot to
support libraries that contain common functionality.

Note that the explicitly pedagogical nature of this example means that some
unusual usages are required (like recompiling the program multiple times for
different datasets).  In addition, the interaction loop of changing the
program parameters via command-line option, generating (on the command-line)
output images, and viewing the images, is very clumsy.  Diderot supports
compiling programs to libraries, around which one can generate nicer GUIs,
but these examples are intended to be as stand-alone and minimal as possible.

Just like [`iso2d`](../iso2d) depends on first creating a dataset with
[`fs2d-scl`](../fs2d), we need to create a volume dataset with [`fs3d-scl`](../fs3d)
in order to compile this program `mip.diderot`:

	../fs3d/fs3d-scl -which 13 -width 3 -sz0 73 -sz1 73 -sz2 73 | unu save -f nrrd -o cube.nrrd

In this case the volume size is chosen to ensure that the local maxima of
this -which 13 synthetic function, a cube frame with maxima at
(x,y,)=(+-1,+1,+1), are actually hit by grid sample points, which helps
reason about subsequent debugging.  We can inspect the volume by tiled slicing
in `unu` to make [cube-tile.png](cube-tile.png):

	unu pad -i cube.nrrd -min 0 0 0 -max M M 79 |
	  unu tile -a 2 0 1 -s 10 8 |
	  unu quantize -b 8 -o cube-tile.png

These examples are mainly to demonstrate Diderot, but some learning some
`unu` is a useful side-effect.  Type `unu` to get a list of all `unu` commands,
and type, for example, `unu pad` to get the usage info for `unu pad`.

We copy `cube.nrrd` to `vol.nrrd`, so that `mip.diderot` can see it in its
default location (in the current directory), we then we can compile

	cp cube.nrrd vol.nrrd
	../../vis12/bin/diderotc --exec mip.diderot

And then make some MIP renderings:

	./mip -out0 0 -camFOV 20 -rayStep 0.03 -iresU 300 -iresV 300
	unu quantize -b 8 -i out.nrrd -o cube-persp.png
	./mip -out0 0 -camFOV 20 -rayStep 0.03 -iresU 300 -iresV 300 -camOrtho true
	unu quantize -b 8 -i out.nrrd -o cube-ortho.png

Make sure you can run these steps to get the same
[cube-persp.png](cube-persp.png) and [cube-ortho.png](cube-ortho.png).  The
only difference in the second image is the use of orthographic instead of
perspective projection, which is apparent in the result.  The `-out0 0`
initialization of the MIP accumulator is safe because `unu minmax cube.nrrd`
reports that the minimum value is 0; this avoids an extra `unu 2op exists`

The usual parameters in volume rendering can be set by command-line options.
This sampling of the various options is not exhaustive; you should try more
yourself.  For the purposes of this demo we set a (shell) variable to store
repeatedly used options, which can be over-ridden later in the command-line.

	OPTS="-out0 0 -camFOV 20 -rayStep 0.03 -iresU 300 -iresV 300"
	./mip $OPTS -camAt 1 1 1;  unu quantize -b 8 -i out.nrrd -o cube-111.png
	./mip $OPTS -camFOV 30;    unu quantize -b 8 -i out.nrrd -o cube-wide.png
	./mip $OPTS -rayStep 0.3;  unu quantize -b 8 -i out.nrrd -o cube-sparse.png
	./mip $OPTS -camFar 0;     unu quantize -b 8 -i out.nrrd -o cube-farclip.png
	./mip $OPTS -camEye 7 7 7; unu quantize -b 8 -i out.nrrd -o cube-eyefar.png

Make sure that these commands generate similar [cube-111.png](cube-111.png),
[cube-wide.png](cube-wide.png), [cube-sparse.png](cube-sparse.png),
[cube-farclip.png](cube-farclip.png), and [cube-eyefar.png](cube-eyefar.png)
for you, and also make sure that you understand why the results look the way they do.

Aside from camera and ray parameters, we now have a basis for testing that
field reconstruction can be invariant with respect to the details of the
sampling grid. We first sample a quadratic function in a few different ways
on a very low resolution grid (which is translating and rotating across these
four datasets):

	../fs3d/fs3d-scl -which 4 -sz0 10 -sz1 10 -sz2 10 | unu save -f nrrd -o parab0.nrrd
	../fs3d/fs3d-scl -which 4 -sz0 10 -sz1 10 -sz2 10 -angle 30 -off 0.25 0.25 0.25 | unu save -f nrrd -o parab1.nrrd
	../fs3d/fs3d-scl -which 4 -sz0 10 -sz1 10 -sz2 10 -angle 60 -off 0.05 0.50 0.50 | unu save -f nrrd -o parab2.nrrd
	../fs3d/fs3d-scl -which 4 -sz0 10 -sz1 10 -sz2 10 -angle 90 -off 0.75 0.75 0.75 | unu save -f nrrd -o parab3.nrrd

This demo (of volume rendering with sampling grid invariance) involves a very
unusual invocation of the Diderot compiler and execution of the resulting
program, so the following steps are unfortunately rather clumsy: each
different grid orientation requires a recompilation.  First, we use the
program as is, with the

	field#0(3)[] F = vol ⊛ tent;

line **uncommented out**, and then run some sh/bash commands:

	OPTS="-out0 0 -inSphere true -camOrtho true -camEye 8 0 0 -camFOV 15 -rayStep 0.01 -iresU 300 -iresV 300"
	for I in 0 1 2 3; do
	  echo === $I ====
	  cp parab${I}.nrrd vol.nrrd
	  ../../vis12/bin/diderotc --exec mip.diderot
	  ./mip $OPTS
	  unu quantize -b 8 -min 0 -max 1 -i out.nrrd -o parab${I}-tent.png

This makes orthographic MIPs of the four low-res parab?.nrrd datasets, producing
images that reveal the underlying sampling grid:

The fact that these images differ means that the quadratic function sampled
by `fs3d-scl -which 4` is not exactly reconstructed by the `tent`
reconstruction kernel.  If the reconstruction had been exact, then the
rendering would have been entirely determined by the underlying quadratic
function, rather than the particulars of the sampling grid.  Note that the
`inSphere` parameter is important here: it ensures that the region rendered
(the sphere) has the same rotational symmetry as the quadratic function

Now **we change one line of code below**, so that only the

	field#1(3)[] F = vol ⊛ ctmr;

field definition is uncommented (commenting out `field#0(3)[] F = vol ⊛ tent;`).
This changes the reconstruction kernel the Catmull-Rom cubic interpolating spline,
which can accurately reconstruct quadratic functions.
Then we re-run the commands above, but saving the results to differently-named images.

	OPTS="-out0 0 -inSphere true -camOrtho true -camEye 8 0 0 -camFOV 15 -rayStep 0.01 -iresU 300 -iresV 300"
	for I in 0 1 2 3; do
	  echo === $I ====
	  cp parab${I}.nrrd vol.nrrd
	  ../../vis12/bin/diderotc --exec mip.diderot
	  ./mip $OPTS
	  unu quantize -b 8 -min 0 -max 1 -i out.nrrd -o parab${I}-ctmr.png

This generates four new images:

As should be visually clear, these images are all the same (or very nearly so), demonstrating the
accuracy of the Catmull-Rom spline for quadratic functions.

Things to try by further modifications of this program:
* Change the rendering from being a maximum intensity, to a summation intensity
(e.g. `out += F(pos)`). Make sure `out0` is zero.
* Change rendering to be mean intensity, which requires adding a counter to
count the number of samples of current ray inside the volume.
* Change the rendering to summating of gradient magnitude
(e.g. `out += |∇F(pos)|`).
* Change the rendering to generate a `vec3` output summing gradient vectors
(e.g. `out += ∇F(pos)`).  Type of `out` must be `vec3`, initialized to `[0,0,0]`.

========================================== */

input image(3)[] vol ("volume dataset to MIP") = image("vol.nrrd");
/* Look-from and look-at are points, but up is a vector;
   all are stored in vec3s.  Diderot has not adopted the
   mathematical vector-vs-point distinction into its type system */
input vec3 camEye ("camera look-from point") = [6, 9, 2];
input vec3 camAt ("camera look-at point") = [0, 0, 0];
input vec3 camUp ("camera pseudo-up vector") = [0, 0, 1];
input real camNear ("relative to look-at point, distance to near clipping plane (where rays start from)") = -3;
input real camFar ("relative to look-at point, distance to far clipping plane") = 3;
input real camFOV ("field-of-view angle subtended vertically by view plane") = 15;
input bool camOrtho ("whether to use orthographic, instead of perspective, projection") = false;
input int iresU ("# samples on horizontal axis of image") = 640;
input int iresV ("# samples on vertical axis of image") = 480;
input real rayStep ("inter-sample step along view direction") = 0.1;
input real out0 ("value at which to initilize output max accumulator; using -inf ensures that output will stay -inf if ray misses volume entirely") = -∞;

// this is important for this example, but not other volume renderers
input bool inSphere ("only render samples inside a unit sphere") = false;

/* Convolve volume data with one of various possible kernels. Typically
   Diderot programs use one particular kernel suited for the task; knowing
   the kernel at compile time permits the evaluation of the kernel and its
   derivatives to be expressed in terms of specific constants rather via
   further indirection.  There is as yet not a good way to specify a kernel
   on the command-line, or assign one field to another at runtime, hence
   the need to uncomment one of these at a time */
//field#0(3)[] F = vol ⊛ tent;
field#1(3)[] F = vol ⊛ ctmr;
//field#2(3)[] F = vol ⊛ bspln3;
//field#4(3)[] F = vol ⊛ c4hexic;

// (boilerplate) computing ray parameters and view-space basis
vec3 camN = normalize(camAt - camEye);  // N: away from eye
vec3 camU = normalize(camN × camUp);    // U: right
vec3 camV = camN × camU;                // V: down
real camDist = |camAt - camEye|;
real camVmax = tan(camFOV*π/360)*camDist;
real camUmax = camVmax*iresU/iresV;
real camNearVsp = camNear + camDist; // near clip, view space
real camFarVsp = camFar + camDist;   // far clip, view space

// how to compute MIP of ray through (rayU,rayV) on view plane
strand raycast(real rayU, real rayV) {
   real rayN = camNearVsp;
   // offset from view-plane center to where ray hits it
   vec3 UV = rayU*camU + rayV*camV;
   // where ray starts (ray position at N=0)
   vec3 rayOrig = camEye + (UV if camOrtho else [0,0,0]);
   // the vector to parameterize position along ray for this pixel
   vec3 rayVec = camN + ([0,0,0] if camOrtho else UV/camDist);
   // initialize output value
   output real out = out0;

   update { // how to compute one sample of a MIP
      vec3 pos = rayOrig + rayN*rayVec;  // pos = ray sample position
      if (inside(pos,F)     // can only query field inside its domain
          && (!inSphere || |pos| < 1)) { // and maybe only in unit sphere
         out = max(out, F(pos));  // update output based on last sample
      if (rayN > camFarVsp) {  // ray hit the far clipping plane
      rayN += rayStep;  // increment ray position

/* this creates a cell-centered sampling of the view plane */
initially [ raycast(lerp(-camUmax, camUmax, -0.5, ui, iresU-0.5),
                    lerp(-camVmax, camVmax, -0.5, vi, iresV-0.5))
            | vi in 0..iresV-1, ui in 0..iresU-1 ];

ViewVC Help
Powered by ViewVC 1.0.0