/*

../bin/diderotc --exec isosurf-lapl.diderot \
 && ./isosurf-lapl \
 && overrgb -i rayRGBA.nrrd -b 0.1 0.15 0.2 -o isosurf-lapl.png 

*/

field#4(3)[] FF = bspln5 ⊛ image("../data/mobius.nrrd");

// set camera, image, and rendering parameters
vec3 camEye = [-104.957, -249.445, 297.334];
vec3 camAt = [39.5087, 37.2972, 40.6096];
vec3 camUp = [0.16498, 0.610459, 0.774675];
real camNear = -46.0;
real camFar = 46.0;
real camFOV = 6.56;
int imgResU = 480;
int imgResV = 345;
real rayStep = 1;
int stepsMax = 20;
real epsilon = 0.003;
vec3 lightVspDir = [0.9, -2.5, -2.2];
real phongKa = 0.05;
real phongKd = 0.6;
real phongKs = 0.4;
real phongSp = 60.0;

// boilerplate computation of camera and light info
real camDist = |camAt - camEye|;
real camVspNear = camNear + camDist;
real camVspFar = camFar + camDist;
vec3 camN = normalize(camAt - camEye);
vec3 camU = normalize(camN × camUp);
vec3 camV = camN × camU;
real camVmax = tan(camFOV*π/360.0)*camDist;
real camUmax = camVmax*real(imgResU)/real(imgResV);
vec3 lightDir = normalize(lightVspDir[0]*camU + 
                          lightVspDir[1]*camV + 
                          lightVspDir[2]*camN);

// RayCast(ui,vi) computes rendered color for pixel (ui,vi)
strand RayCast (int ui, int vi) {
    real rayU = lerp(-camUmax, camUmax, -0.5, real(ui), real(imgResU)-0.5);
    real rayV = lerp(-camVmax, camVmax, -0.5, real(vi), real(imgResV)-0.5);
    vec3 rayVec = (camDist*camN + rayU*camU + rayV*camV)/camDist;
    vec3 rayDir = normalize(rayVec);
    bool skipping = true;
    real rayN = camVspNear;
    output vec4 rayRGBA = [0.0, 0.0, 0.0, 0.0];
    int steps = 0;

    vec3 pos0 = [0.0,0.0,0.0]; // HEY unwanted initialization 
    vec3 pos1 = [0.0,0.0,0.0]; // HEY unwanted initialization 
    vec3 pos = [0.0,0.0,0.0]; // HEY unwanted initialization 

    update {
       if (skipping) {
          // skipping along ray, looking for sign change
          if (rayN > camVspFar) {
             stabilize;
          }
          pos0 = camEye + rayN*rayVec;
          pos1 = camEye + (rayN + rayStep)*rayVec;
          pos = (pos0 + pos1)/2.0;
          if (inside(pos0, FF) && inside(pos1, FF)) {
             if (trace(∇⊗∇FF(pos0))*trace(∇⊗∇FF(pos1)) < 0.0
                 && |∇FF(pos)| > 10.0) {
                // segment between pos0 and pos1 intersects isosurface;
                // we're done skipping, initialize root finding
                skipping = false;
             }
          }
          if (skipping) {
             rayN = rayN + rayStep;
          }
       } else {
          // root finding starting at pos
          if (steps > stepsMax) {
             // failure, code green
             rayRGBA = [0.0, 1.0, 0.0, 1.0];
             stabilize;
          }
          if (!inside(pos, FF)) {
             // failure, code red
             rayRGBA = [1.0, 0.0, 0.0, 1.0];
             stabilize;
          }
          vec3 grad = identity[3]:∇⊗∇⊗∇FF(pos);
          if (|grad| == 0.0) {
             // failure, code blue
             rayRGBA = [0.0, 0.0, 1.0, 1.0];
             stabilize;
          }
          real delta = -trace(∇⊗∇FF(pos))/(grad•rayDir);
          if (|delta| < epsilon) {    // we've converged; color the point
             vec3 norm = normalize(grad);
             real ld = max(0.0, norm • lightDir);
             real hd = max(0.0, norm • normalize(lightDir - rayDir));
             real gray = phongKa + phongKd*ld + phongKs*hd^phongSp;
             rayRGBA = [gray, gray, gray, 1.0];
             stabilize;
          }
          pos += delta*rayDir;
          steps += 1;
       }
    }
}

initially [ RayCast(ui, vi) | vi in 0..(imgResV-1), ui in 0..(imgResU-1) ];