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 4826 - (view) (download)

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

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