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

SCM Repository

[diderot] Annotation of /tests/vis15-bugs/dvrMIPbug.diderot
ViewVC logotype

Annotation of /tests/vis15-bugs/dvrMIPbug.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4828 - (view) (download)

1 : glk 4826 input image(3)[] vol ("volume dataset to render") = image("dvrESbug-vol.nrrd");
2 :     /* see ../mip/mip.diderot for everything about camera set-up */
3 :     input vec3 camEye ("camera look-from point") = [6, 9, 2];
4 :     input vec3 camAt ("camera look-at point") = [0, 0, 0];
5 :     input vec3 camUp ("camera pseudo-up vector") = [0, 0, 1];
6 :     input real camNear ("relative to look-at point, distance to near clipping plane (where rays start from)") = -3;
7 :     input real camFar ("relative to look-at point, distance to far clipping plane") = 3;
8 :     input real camFOV ("field-of-view angle (in degrees) subtended vertically by view plane") = 15;
9 :     input bool camOrtho ("whether to use orthographic, instead of perspective, projection") = false;
10 :     input int iresU ("# samples on horizontal axis of image") = 640;
11 :     input int iresV ("# samples on vertical axis of image") = 480;
12 :     input real rayStep ("inter-sample distance along view direction") = 0.1;
13 :    
14 :     input real refStep ("reference (unit) step length, for normalizing opacities") = 0.1;
15 :     input real transp0 ("transparency close enough to 0 to terminate ray") = 0.005;
16 :     input real thick ("approximate thickness (in world-space) of soft isosurface") = 0.1;
17 :     input real maxAlpha ("maximum opacity on rendered surface") = 1;
18 :     input real phongKa ("Blinn-Phong ambient component") = 0.1;
19 :     input real phongKd ("Blinn-Phong diffuse component") = 0.7;
20 :     input real phongKs ("Blinn-Phong specular component") = 0.3;
21 :     input real phongSp ("Blinn-Phong specularity exponent") = 60;
22 :     input vec3 litdir ("direction (non-normalized) towards light source, in (U,V,N) view-space") = [-1, -2, -1];
23 :     input vec3 mcnear ("material color at near clipping plane (for depth cuing)") = [1,1,1];
24 :     input vec3 mcfar ("material color at far clipping plane") = [1,1,1];
25 :    
26 :     /* Convolve volume data with one of various possible kernels;
27 :     see ../mip/mip.diderot for more info */
28 :     //field#1(3)[] V = c1tent ⊛ vol;
29 :     //field#1(3)[] V = ctmr ⊛ vol;
30 :     //field#2(3)[] V = bspln3 ⊛ vol;
31 :     field#4(3)[] V = c4hexic ⊛ vol;
32 :    
33 :     /* create a field to render from the original volume data field */
34 :     field#2(3)[] F = V;
35 :    
36 :     // (boilerplate) computing ray parameters and view-space basis
37 :     vec3 camN = normalize(camAt - camEye); // N: away from eye
38 :     vec3 camU = normalize(camN × camUp); // U: right
39 :     vec3 camV = camN × camU; // V: down (right-handed frame)
40 :     real camDist = |camAt - camEye|;
41 :     real camVmax = tan(camFOV*π/360)*camDist;
42 :     real camUmax = camVmax*iresU/iresV;
43 :     real camNearVsp = camNear + camDist; // near clip, view space
44 :     real camFarVsp = camFar + camDist; // far clip, view space
45 :     // convert light directions from view-space to world-space
46 :     vec3 litwsp = transpose([camU,camV,camN])•normalize(litdir);
47 :    
48 :     /* 2-D opacity function, after Levoy. The 1.4 is a trick to ensure
49 :     that there is a segment of positions receiving maximum opacity. */
50 :     function real alpha(real v, real g) = maxAlpha*clamp(0, 1, 1.4 - |v|/(g*thick));
51 :    
52 :     function vec3 nullspace1D(tensor[3,3] m) {
53 :     vec3 c01 = m[0,:]×m[1,:];
54 :     vec3 c02 = m[0,:]×m[2,:];
55 :     vec3 c12 = m[1,:]×m[2,:];
56 :     if (c02•c01 < 0) { c02 = -c02; }
57 :     if (c12•(c01+c02) < 0) { c12 = -c12; }
58 :     return normalize(c01 + c02 + c12);
59 :     }
60 :    
61 :     function real mode(real v0_, real v1_, real v2_) {
62 :     real avg = (v0_ + v1_ + v2_)/3;
63 :     real v0 = v0_ - avg;
64 :     real v1 = v1_ - avg;
65 :     real v2 = v2_ - avg;
66 :     real num = (v0 + v1 - 2*v2)*(2*v0 - v1 - v2)*(v0 - 2*v1 + v2);
67 :     real den = v0*v0 + v1*v1 + v2*v2 - v1*v2 - v0*v1 - v0*v2;
68 :     den = sqrt(den);
69 :     return num/(2*den*den*den) if den != 0 else 0;
70 :     }
71 :    
72 :     function vec3 ridgeStep(vec3 pp) {
73 :     tensor[3,3] hess = ∇⊗∇F(pp);
74 :     real{3} eval = evals(hess);
75 :     vec3 evec1 = [0,0,0];
76 :     vec3 evec2 = [0,0,0];
77 :     if (false) {
78 :     vec3{3} evec = evecs(hess);
79 :     evec1 = evec{1};
80 :     evec2 = evec{2};
81 :     } else {
82 :     evec1 = nullspace1D(hess - identity[3]*eval{1});
83 :     evec2 = nullspace1D(hess - identity[3]*eval{2});
84 :     }
85 :     return (evec1⊗evec1/eval{1} + evec2⊗evec2/eval{2})•∇F(pp);
86 :     }
87 :    
88 :     // how to render ray through (rayU,rayV) on view plane
89 :     strand raycast(int ui, int vi /* real rayU, real rayV */) {
90 :     // cell-centered sampling of view plane (intersects look-at)
91 :     real rayU = lerp(-camUmax, camUmax, -0.5, ui, iresU-0.5);
92 :     real rayV = lerp(-camVmax, camVmax, -0.5, vi, iresV-0.5);
93 :     // creation of per-strand ray state based on ../mip/mip.diderot
94 :     real rayN = camNearVsp;
95 :     vec3 UV = rayU*camU + rayV*camV;
96 :     vec3 rayOrig = camEye + (UV if camOrtho else [0,0,0]);
97 :     vec3 rayVec = camN + ([0,0,0] if camOrtho else UV/camDist);
98 :     /* alphaFix is used for opacity correction, so that attenuation
99 :     happens in world-space rather than ray sample index space.
100 :     The actual distance between samples on this ray is |rayVec|*rayStep */
101 :     real alphaFix = |rayVec|*rayStep/refStep;
102 :     vec3 eyeDir = -normalize(rayVec);
103 :     // output for this ray
104 :     output vec4 rgba = [0,0,0,0];
105 :     // state for this ray is current color ...
106 :     vec3 rgb = [0,0,0];
107 :     // ... and current tranparency
108 :     real transp = 1;
109 :     real maxval = -∞;
110 :     update {
111 :     rayN += rayStep; // increment ray position
112 :     if (rayN > camFarVsp) { // ray passed the far clipping plane
113 :     stabilize;
114 :     }
115 :     vec3 pos = rayOrig + rayN*rayVec; // pos is ray sample position
116 :     if (!inside(pos,F)) { // If not inside field domain,
117 :     continue; // then move onto next iteration
118 :     }
119 :     // compute data contribution
120 :     real val = F(pos);
121 : jhr 4828 maxval = max(maxval, val);
122 :     continue; /* This continue causes a compiler crash */
123 :     tensor[3,3] hess = ∇⊗∇F(pos);
124 :     real{3} eval = evals(hess);
125 :     vec3 delta = ridgeStep(pos);
126 :     real aa = alpha(|delta|, 1) if eval{1} < -2.5 else 0;
127 :     if (aa == 0) {
128 :     continue;
129 : glk 4826 }
130 : jhr 4828 vec3 grad = -delta;
131 : glk 4826 aa = 1 - (1 - aa)^alphaFix;
132 :     // Note: not standard Phong diffuse (no dark hemisphere)
133 :     real dcomp = lerp(0, 1, -1, normalize(-grad)•litwsp, 1)^3;
134 :     // If phongKs=0, this is conditional expression speeds things slightly
135 :     real scomp = max(0,normalize(-grad)•normalize(eyeDir+litwsp))^phongSp
136 :     if phongKs != 0 else 0.0;
137 :     // simple depth-cueing
138 :     vec3 dcol = lerp(mcnear, mcfar, camNearVsp, rayN, camFarVsp);
139 :     vec3 mcol = [1,1,1];
140 :     // light color is currently [1,1,1]
141 :     rgb += transp*aa*((phongKa + phongKd*dcomp)*modulate(dcol,mcol)
142 :     + phongKs*scomp*[1,1,1]);
143 :     transp *= 1 - aa;
144 :     if (transp < transp0) { // early ray termination
145 :     transp = 0;
146 :     stabilize;
147 :     }
148 :     }
149 :     }
150 :    
151 :     initially [ raycast(ui, vi)
152 :     | vi in 0..iresV-1,
153 :     ui in 0..iresU-1];
154 :    

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