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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4826 - (download) (annotate)
Wed Nov 2 11:39:10 2016 UTC (2 years, 10 months ago) by glk
File size: 7427 byte(s)
another crash
input image(3)[] vol ("volume dataset to render") = image("dvrESbug-vol.nrrd");
input bool doIso ("do isosurface instead of ridge line") = false;
/* see ../mip/mip.diderot for everything about camera set-up */
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 distance along view direction") = 0.1;

input real refStep ("reference (unit) step length, for normalizing opacities") = 0.1;
input real transp0 ("transparency close enough to 0 to terminate ray") = 0.005;
input real isoval ("isovalue at which to render soft isosurface") = 0;
input real thick ("approximate thickness (in world-space) of soft isosurface") = 0.1;
input real maxAlpha ("maximum opacity on rendered surface") = 1;
input real phongKa ("Blinn-Phong ambient component") = 0.1;
input real phongKd ("Blinn-Phong diffuse component") = 0.7;
input real phongKs ("Blinn-Phong specular component") = 0.3;
input real phongSp ("Blinn-Phong specularity exponent") = 60;
input vec3 litdir ("direction (non-normalized) towards light source, in (U,V,N) view-space") = [-1, -2, -1];
input vec3 mcnear ("material color at near clipping plane (for depth cuing)") = [1,1,1];
input vec3 mcfar ("material color at far clipping plane") = [1,1,1];

/* Convolve volume data with one of various possible kernels;
   see ../mip/mip.diderot for more info */
//field#1(3)[] V = c1tent ⊛ vol;
//field#1(3)[] V = ctmr ⊛ vol;
//field#2(3)[] V = bspln3 ⊛ vol;
field#4(3)[] V = c4hexic ⊛ vol;

/* create a field to render from the original volume data field */
field#2(3)[] F = V;

// (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 (right-handed frame)
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
// convert light directions from view-space to world-space
vec3 litwsp = transpose([camU,camV,camN])•normalize(litdir);

/* 2-D opacity function, after Levoy. The 1.4 is a trick to ensure
  that there is a segment of positions receiving maximum opacity. */
function real alpha(real v, real g) = maxAlpha*clamp(0, 1, 1.4 - |v|/(g*thick));

function vec3 nullspace1D(tensor[3,3] m) {
   vec3 c01 = m[0,:]×m[1,:];
   vec3 c02 = m[0,:]×m[2,:];
   vec3 c12 = m[1,:]×m[2,:];
   if (c02•c01 < 0) { c02 = -c02; }
   if (c12•(c01+c02) < 0) { c12 = -c12; }
   return normalize(c01 + c02 + c12);
}

function real mode(real v0_, real v1_, real v2_) {
   real avg = (v0_ + v1_ + v2_)/3;
   real v0 = v0_ - avg;
   real v1 = v1_ - avg;
   real v2 = v2_ - avg;
   real num = (v0 + v1 - 2*v2)*(2*v0 - v1 - v2)*(v0 - 2*v1 + v2);
   real den = v0*v0 + v1*v1 + v2*v2 - v1*v2 - v0*v1 - v0*v2;
   den = sqrt(den);
   return num/(2*den*den*den) if den != 0 else 0;
}

function vec3 ridgeStep(vec3 pp) {
   tensor[3,3] hess = ∇⊗∇F(pp);
   real{3} eval = evals(hess);
   vec3 evec1 = [0,0,0];
   vec3 evec2 = [0,0,0];
   if (false) {
      vec3{3} evec = evecs(hess);
      evec1 = evec{1};
      evec2 = evec{2};
   } else {
      evec1 = nullspace1D(hess - identity[3]*eval{1});
      evec2 = nullspace1D(hess - identity[3]*eval{2});
   }
   return (evec1⊗evec1/eval{1} + evec2⊗evec2/eval{2})•∇F(pp);
}

// how to render ray through (rayU,rayV) on view plane
strand raycast(int ui, int vi /* real rayU, real rayV */) {
   // cell-centered sampling of view plane (intersects look-at)
   real rayU = lerp(-camUmax, camUmax, -0.5, ui, iresU-0.5);
   real rayV = lerp(-camVmax, camVmax, -0.5, vi, iresV-0.5);
   // creation of per-strand ray state based on ../mip/mip.diderot
   real rayN = camNearVsp;
   vec3 UV = rayU*camU + rayV*camV;
   vec3 rayOrig = camEye + (UV if camOrtho else [0,0,0]);
   vec3 rayVec = camN + ([0,0,0] if camOrtho else UV/camDist);
   /* alphaFix is used for opacity correction, so that attenuation
   happens in world-space rather than ray sample index space.
   The actual distance between samples on this ray is |rayVec|*rayStep */
   real alphaFix = |rayVec|*rayStep/refStep;
   vec3 eyeDir = -normalize(rayVec);
   // output for this ray
   output vec4 rgba = [0,0,0,0];
   // state for this ray is current color ...
   vec3 rgb = [0,0,0];
   // ... and current tranparency
   real transp = 1;
   /* example of turning on debuging for one strand (pixel);
      can can then have "if (verb) { print(...); }" */
   bool verb = 200==ui && 220==vi;
   real maxval = -∞;
   update {
      rayN += rayStep;          // increment ray position
      if (rayN > camFarVsp) {   // ray passed the far clipping plane
         stabilize;
      }
      vec3 pos = rayOrig + rayN*rayVec;  // pos is ray sample position
      if (!inside(pos,F)) {              // If not inside field domain,
         continue;                       // then move onto next iteration
      }
      // compute data contribution
      real val = F(pos);
      maxval = max(maxval, val); continue;
      vec3 grad = ∇F(pos);
      real aa = 0;
      if (doIso) {
         aa = alpha(val-isoval, |grad|);
      } else {
         tensor[3,3] hess = ∇⊗∇F(pos);
         real{3} eval = evals(hess);
         vec3 delta = ridgeStep(pos);
         aa = alpha(|delta|, 1) if eval{1} < -2.5 else 0;
         if (aa == 0) {
            continue;
         }
         grad = -delta;
      }
      aa = 1 - (1 - aa)^alphaFix;
      // Note: not standard Phong diffuse (no dark hemisphere)
      real dcomp = lerp(0, 1, -1, normalize(-grad)•litwsp, 1)^3;
      // If phongKs=0, this is conditional expression speeds things slightly
      real scomp = max(0,normalize(-grad)•normalize(eyeDir+litwsp))^phongSp
                   if phongKs != 0 else 0.0;
      // simple depth-cueing
      vec3 dcol = lerp(mcnear, mcfar, camNearVsp, rayN, camFarVsp);
      vec3 mcol = [1,1,1];
      // light color is currently [1,1,1]
      if (verb) {
         print("(", rayN, "|", pos, "): val=", val, "; dcol=", dcol, "; mcol=", mcol, "\n");
      }
      rgb += transp*aa*((phongKa + phongKd*dcomp)*modulate(dcol,mcol)
                        + phongKs*scomp*[1,1,1]);
      transp *= 1 - aa;
      if (transp < transp0) { // early ray termination
         transp = 0;
         stabilize;
      }
   }
   stabilize {
      if (transp < 1) {  // un-pre-multiply opacities
         real aa = 1-transp;
         rgba = [rgb[0]/aa, rgb[1]/aa, rgb[2]/aa, aa];
      } else {
         rgba = [rgb[0], rgb[1], rgb[2], 0];
      }
      rgba = [maxval, 0, 0, 1];
   }
}

initially [ raycast(ui, vi)
            | vi in 0..iresV-1,
              ui in 0..iresU-1];


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