SCM Repository
[diderot] / tests / examples / mip / mip.diderot |
View of /tests/examples/mip/mip.diderot
Parent Directory | Revision Log
Revision 4640 -
(download)
(annotate)
Tue Sep 27 20:54:47 2016 UTC (2 years, 10 months ago) by glk
File size: 11849 byte(s)
Tue Sep 27 20:54:47 2016 UTC (2 years, 10 months ago) by glk
File size: 11849 byte(s)
initial result of svn export --username anonsvn --password=anonsvn https://svn.smlnj-gforge.cs.uchicago.edu/svn/diderot/branches/vis15/src/tests/
/* ========================================== ## 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,z)=(+-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), and then we can compile `mip.diderot`: 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` step. 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 done This makes orthographic MIPs of the four low-res parab?.nrrd datasets, producing images that reveal the underlying sampling grid: [parab0-tent.png](parab0-tent.png), [parab1-tent.png](parab1-tent.png), [parab2-tent.png](parab2-tent.png), [parab3-tent.png](parab3-tent.png). 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 itself. 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 done This generates four new images: [parab0-ctmr.png](parab0-ctmr.png), [parab1-ctmr.png](parab1-ctmr.png), [parab2-ctmr.png](parab2-ctmr.png), [parab3-ctmr.png](parab3-ctmr.png). 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 (in degrees) 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; // Currently, Unicode is not allowed in Diderot strings, // hence using "-inf" instead of "-∞" in this usage info. 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 stabilize; } 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 ];
root@smlnj-gforge.cs.uchicago.edu | ViewVC Help |
Powered by ViewVC 1.0.0 |