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

SCM Repository

[diderot] Annotation of /tests/examples/mip/mip.diderot
ViewVC logotype

Annotation of /tests/examples/mip/mip.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4640 - (view) (download)

1 : glk 4640 /* ==========================================
2 :     ## mip.diderot: basic maximum intensity projection (MIP) volume rendering
3 :    
4 :     Maximum intensity projection is the minimalist volume visualization
5 :     tool. This implementation is about as short as possible; much of the code is
6 :     spent on specifying and setting up the ray geometry. It can also be adapted
7 :     (by changing a line or two) to do other kinds of projections. The code for
8 :     camera geometry and ray sampling will be re-used verbatim in other volume
9 :     rendering programs. This may later become a target for evolving Diderot to
10 :     support libraries that contain common functionality.
11 :    
12 :     Note that the explicitly pedagogical nature of this example means that some
13 :     unusual usages are required (like recompiling the program multiple times for
14 :     different datasets). In addition, the interaction loop of changing the
15 :     program parameters via command-line option, generating (on the command-line)
16 :     output images, and viewing the images, is very clumsy. Diderot supports
17 :     compiling programs to libraries, around which one can generate nicer GUIs,
18 :     but these examples are intended to be as stand-alone and minimal as possible.
19 :    
20 :     Just like [`iso2d`](../iso2d) depends on first creating a dataset with
21 :     [`fs2d-scl`](../fs2d), we need to create a volume dataset with [`fs3d-scl`](../fs3d)
22 :     in order to compile this program `mip.diderot`:
23 :    
24 :     ../fs3d/fs3d-scl -which 13 -width 3 -sz0 73 -sz1 73 -sz2 73 | unu save -f nrrd -o cube.nrrd
25 :    
26 :     In this case the volume size is chosen to ensure that the local maxima of
27 :     this -which 13 synthetic function, a cube frame with maxima at
28 :     (x,y,z)=(+-1,+1,+1), are actually hit by grid sample points, which helps
29 :     reason about subsequent debugging. We can inspect the volume by tiled slicing
30 :     in `unu` to make [cube-tile.png](cube-tile.png):
31 :    
32 :     unu pad -i cube.nrrd -min 0 0 0 -max M M 79 |
33 :     unu tile -a 2 0 1 -s 10 8 |
34 :     unu quantize -b 8 -o cube-tile.png
35 :    
36 :     These examples are mainly to demonstrate Diderot, but some learning some
37 :     `unu` is a useful side-effect. Type `unu` to get a list of all `unu` commands,
38 :     and type, for example, `unu pad` to get the usage info for `unu pad`.
39 :    
40 :     We copy `cube.nrrd` to `vol.nrrd`, so that `mip.diderot` can see it in its
41 :     default location (in the current directory), and then we can compile
42 :     `mip.diderot`:
43 :    
44 :     cp cube.nrrd vol.nrrd
45 :     ../../vis12/bin/diderotc --exec mip.diderot
46 :    
47 :     And then make some MIP renderings:
48 :    
49 :     ./mip -out0 0 -camFOV 20 -rayStep 0.03 -iresU 300 -iresV 300
50 :     unu quantize -b 8 -i out.nrrd -o cube-persp.png
51 :     ./mip -out0 0 -camFOV 20 -rayStep 0.03 -iresU 300 -iresV 300 -camOrtho true
52 :     unu quantize -b 8 -i out.nrrd -o cube-ortho.png
53 :    
54 :     Make sure you can run these steps to get the same
55 :     [cube-persp.png](cube-persp.png) and [cube-ortho.png](cube-ortho.png). The
56 :     only difference in the second image is the use of orthographic instead of
57 :     perspective projection, which is apparent in the result. The `-out0 0`
58 :     initialization of the MIP accumulator is safe because `unu minmax cube.nrrd`
59 :     reports that the minimum value is 0; this avoids an extra `unu 2op exists`
60 :     step.
61 :    
62 :     The usual parameters in volume rendering can be set by command-line options.
63 :     This sampling of the various options is not exhaustive; you should try more
64 :     yourself. For the purposes of this demo we set a (shell) variable to store
65 :     repeatedly used options, which can be over-ridden later in the command-line.
66 :    
67 :     OPTS="-out0 0 -camFOV 20 -rayStep 0.03 -iresU 300 -iresV 300"
68 :     ./mip $OPTS -camAt 1 1 1; unu quantize -b 8 -i out.nrrd -o cube-111.png
69 :     ./mip $OPTS -camFOV 30; unu quantize -b 8 -i out.nrrd -o cube-wide.png
70 :     ./mip $OPTS -rayStep 0.3; unu quantize -b 8 -i out.nrrd -o cube-sparse.png
71 :     ./mip $OPTS -camFar 0; unu quantize -b 8 -i out.nrrd -o cube-farclip.png
72 :     ./mip $OPTS -camEye 7 7 7; unu quantize -b 8 -i out.nrrd -o cube-eyefar.png
73 :    
74 :     Make sure that these commands generate similar [cube-111.png](cube-111.png),
75 :     [cube-wide.png](cube-wide.png), [cube-sparse.png](cube-sparse.png),
76 :     [cube-farclip.png](cube-farclip.png), and [cube-eyefar.png](cube-eyefar.png)
77 :     for you, and also make sure that you understand why the results look the way they do.
78 :    
79 :     Aside from camera and ray parameters, we now have a basis for testing that
80 :     field reconstruction can be invariant with respect to the details of the
81 :     sampling grid. We first sample a quadratic function in a few different ways
82 :     on a very low resolution grid (which is translating and rotating across these
83 :     four datasets):
84 :    
85 :     ../fs3d/fs3d-scl -which 4 -sz0 10 -sz1 10 -sz2 10 | unu save -f nrrd -o parab0.nrrd
86 :     ../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
87 :     ../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
88 :     ../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
89 :    
90 :     This demo (of volume rendering with sampling grid invariance) involves a very
91 :     unusual invocation of the Diderot compiler and execution of the resulting
92 :     program, so the following steps are unfortunately rather clumsy: each
93 :     different grid orientation requires a recompilation. First, we use the
94 :     program as is, with the
95 :    
96 :     field#0(3)[] F = vol ⊛ tent;
97 :    
98 :     line **uncommented out**, and then run some sh/bash commands:
99 :    
100 :     OPTS="-out0 0 -inSphere true -camOrtho true -camEye 8 0 0 -camFOV 15 -rayStep 0.01 -iresU 300 -iresV 300"
101 :     for I in 0 1 2 3; do
102 :     echo === $I ====
103 :     cp parab${I}.nrrd vol.nrrd
104 :     ../../vis12/bin/diderotc --exec mip.diderot
105 :     ./mip $OPTS
106 :     unu quantize -b 8 -min 0 -max 1 -i out.nrrd -o parab${I}-tent.png
107 :     done
108 :    
109 :     This makes orthographic MIPs of the four low-res parab?.nrrd datasets, producing
110 :     images that reveal the underlying sampling grid:
111 :     [parab0-tent.png](parab0-tent.png),
112 :     [parab1-tent.png](parab1-tent.png),
113 :     [parab2-tent.png](parab2-tent.png),
114 :     [parab3-tent.png](parab3-tent.png).
115 :    
116 :     The fact that these images differ means that the quadratic function sampled
117 :     by `fs3d-scl -which 4` is not exactly reconstructed by the `tent`
118 :     reconstruction kernel. If the reconstruction had been exact, then the
119 :     rendering would have been entirely determined by the underlying quadratic
120 :     function, rather than the particulars of the sampling grid. Note that the
121 :     `inSphere` parameter is important here: it ensures that the region rendered
122 :     (the sphere) has the same rotational symmetry as the quadratic function
123 :     itself.
124 :    
125 :     Now **we change one line of code below**, so that only the
126 :    
127 :     field#1(3)[] F = vol ⊛ ctmr;
128 :    
129 :     field definition is uncommented (commenting out `field#0(3)[] F = vol ⊛ tent;`).
130 :     This changes the reconstruction kernel the Catmull-Rom cubic interpolating spline,
131 :     which can accurately reconstruct quadratic functions.
132 :     Then we re-run the commands above, but saving the results to differently-named images.
133 :    
134 :     OPTS="-out0 0 -inSphere true -camOrtho true -camEye 8 0 0 -camFOV 15 -rayStep 0.01 -iresU 300 -iresV 300"
135 :     for I in 0 1 2 3; do
136 :     echo === $I ====
137 :     cp parab${I}.nrrd vol.nrrd
138 :     ../../vis12/bin/diderotc --exec mip.diderot
139 :     ./mip $OPTS
140 :     unu quantize -b 8 -min 0 -max 1 -i out.nrrd -o parab${I}-ctmr.png
141 :     done
142 :    
143 :     This generates four new images:
144 :     [parab0-ctmr.png](parab0-ctmr.png),
145 :     [parab1-ctmr.png](parab1-ctmr.png),
146 :     [parab2-ctmr.png](parab2-ctmr.png),
147 :     [parab3-ctmr.png](parab3-ctmr.png).
148 :    
149 :     As should be visually clear, these images are all the same (or very nearly so), demonstrating the
150 :     accuracy of the Catmull-Rom spline for quadratic functions.
151 :    
152 :     Things to try by further modifications of this program:
153 :     * Change the rendering from being a maximum intensity, to a summation intensity
154 :     (e.g. `out += F(pos)`). Make sure `out0` is zero.
155 :     * Change rendering to be mean intensity, which requires adding a counter to
156 :     count the number of samples of current ray inside the volume.
157 :     * Change the rendering to summating of gradient magnitude
158 :     (e.g. `out += |∇F(pos)|`).
159 :     * Change the rendering to generate a `vec3` output summing gradient vectors
160 :     (e.g. `out += ∇F(pos)`). Type of `out` must be `vec3`, initialized to `[0,0,0]`.
161 :    
162 :     ========================================== */
163 :    
164 :     input image(3)[] vol ("volume dataset to MIP") = image("vol.nrrd");
165 :     /* Look-from and look-at are points, but up is a vector;
166 :     all are stored in vec3s. Diderot has not adopted the
167 :     mathematical vector-vs-point distinction into its type system */
168 :     input vec3 camEye ("camera look-from point") = [6, 9, 2];
169 :     input vec3 camAt ("camera look-at point") = [0, 0, 0];
170 :     input vec3 camUp ("camera pseudo-up vector") = [0, 0, 1];
171 :     input real camNear ("relative to look-at point, distance to near clipping plane (where rays start from)") = -3;
172 :     input real camFar ("relative to look-at point, distance to far clipping plane") = 3;
173 :     input real camFOV ("field-of-view angle (in degrees) subtended vertically by view plane") = 15;
174 :     input bool camOrtho ("whether to use orthographic, instead of perspective, projection") = false;
175 :     input int iresU ("# samples on horizontal axis of image") = 640;
176 :     input int iresV ("# samples on vertical axis of image") = 480;
177 :     input real rayStep ("inter-sample step along view direction") = 0.1;
178 :     // Currently, Unicode is not allowed in Diderot strings,
179 :     // hence using "-inf" instead of "-∞" in this usage info.
180 :     input real out0 ("value at which to initilize output max accumulator; using -inf ensures that output will stay -inf if ray misses volume entirely") = -∞;
181 :    
182 :     // this is important for this example, but not other volume renderers
183 :     input bool inSphere ("only render samples inside a unit sphere") = false;
184 :    
185 :     /* Convolve volume data with one of various possible kernels. Typically
186 :     Diderot programs use one particular kernel suited for the task; knowing
187 :     the kernel at compile time permits the evaluation of the kernel and its
188 :     derivatives to be expressed in terms of specific constants rather via
189 :     further indirection. There is as yet not a good way to specify a kernel
190 :     on the command-line, or assign one field to another at runtime, hence
191 :     the need to uncomment one of these at a time */
192 :     //field#0(3)[] F = vol ⊛ tent;
193 :     field#1(3)[] F = vol ⊛ ctmr;
194 :     //field#2(3)[] F = vol ⊛ bspln3;
195 :     //field#4(3)[] F = vol ⊛ c4hexic;
196 :    
197 :     // (boilerplate) computing ray parameters and view-space basis
198 :     vec3 camN = normalize(camAt - camEye); // N: away from eye
199 :     vec3 camU = normalize(camN × camUp); // U: right
200 :     vec3 camV = camN × camU; // V: down
201 :     real camDist = |camAt - camEye|;
202 :     real camVmax = tan(camFOV*π/360)*camDist;
203 :     real camUmax = camVmax*iresU/iresV;
204 :     real camNearVsp = camNear + camDist; // near clip, view space
205 :     real camFarVsp = camFar + camDist; // far clip, view space
206 :    
207 :     // how to compute MIP of ray through (rayU,rayV) on view plane
208 :     strand raycast(real rayU, real rayV) {
209 :     real rayN = camNearVsp;
210 :     // offset from view-plane center to where ray hits it
211 :     vec3 UV = rayU*camU + rayV*camV;
212 :     // where ray starts (ray position at N=0)
213 :     vec3 rayOrig = camEye + (UV if camOrtho else [0,0,0]);
214 :     // the vector to parameterize position along ray for this pixel
215 :     vec3 rayVec = camN + ([0,0,0] if camOrtho else UV/camDist);
216 :     // initialize output value
217 :     output real out = out0;
218 :    
219 :     update { // how to compute one sample of a MIP
220 :     vec3 pos = rayOrig + rayN*rayVec; // pos = ray sample position
221 :     if (inside(pos,F) // can only query field inside its domain
222 :     && (!inSphere || |pos| < 1)) { // and maybe only in unit sphere
223 :     out = max(out, F(pos)); // update output based on last sample
224 :     }
225 :     if (rayN > camFarVsp) { // ray hit the far clipping plane
226 :     stabilize;
227 :     }
228 :     rayN += rayStep; // increment ray position
229 :     }
230 :     }
231 :    
232 :     /* this creates a cell-centered sampling of the view plane */
233 :     initially [ raycast(lerp(-camUmax, camUmax, -0.5, ui, iresU-0.5),
234 :     lerp(-camVmax, camVmax, -0.5, vi, iresV-0.5))
235 :     | vi in 0..iresV-1, ui in 0..iresU-1 ];

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