// quad-k1k2: demo of curvature calculation on quad patch dataset
// to process output:
// unu reshape -i mip.txt -s 3 200 200 | unu quantize -b 8 -o zslice-k1k2.png
//

int imgRes = 200;
int dataRes = 110;
vec3 orig = [0.0, 0.0, 10.0];
vec3 uVec = [1.0, 0.0, 0.0];
vec3 vVec = [0.0, 1.0, 0.0];
image(3)[] vol = load("../data/quad-patches-pad.nrrd");
field#2(3)[] F = vol ⊛ bspln3;
// HEY (want) some kind of error upon redeclaration
// of variables; was getting really confused by error messages
// associated with having used G both for green image (here) 
// and geometry tensor (below)
image(2)[] imgR = load("../data/hue-0.nrrd"); field#0(2)[] RR = imgR ⊛ tent;
image(2)[] imgG = load("../data/hue-1.nrrd"); field#0(2)[] GG = imgG ⊛ tent;
image(2)[] imgB = load("../data/hue-2.nrrd"); field#0(2)[] BB = imgB ⊛ tent;
//image(2)[3] imgRGB = load("../data/hue.nhdr"); field#0(2)[3] RGB = imgRGB ⊛ tent;

strand RayCast (int ui, int vi)
{
    real uu = lerp(-0.5, real(dataRes)-0.5, -0.5, real(ui), real(imgRes)-0.5);
    real vv = lerp(-0.5, real(dataRes)-0.5, -0.5, real(vi), real(imgRes)-0.5);
    output vec3 rgb = [0.0,0.0,0.0];

    update
    {
        vec3 pos = orig + uu*uVec + vv*vVec;
        if (inside (pos,F)) {
           vec3 g = -∇F@pos;
           // the "if (g[2] > 0.10)" test is specific to this dataset;
           // it assesses if we're within one of the quadratic patches
           if (g[2] > 0.10) {
              vec3 n = normalize(g);
              // HEY (BUG) using "-∇(∇F)@pos" instead of "-1.0*∇(∇F)@pos" causes:
              // uncaught exception Fail [Fail: Error in compiling zslice-k1k2.diderot]
              //   raised at driver/main.sml:31.39-31.76
              tensor[3,3] H = -1.0*∇(∇F)@pos; // replace with -∇(∇F)@pos
              tensor[3,3] P = identity[3] - n⊗n;
              tensor[3,3] G = (P•H•P)/|g|;
              real disc = max(0.0, sqrt(2.0*|G|^2 - trace(G)^2));
              real k1 = (trace(G) + disc)/2.0;
              real k2 = (trace(G) - disc)/2.0;
              // correct curvature computation ends here
              k1 = max(-0.99, min(0.99, 4.2*k1));
              k2 = max(-0.99, min(0.99, 4.2*k2));
              vec2 kk = [k1,k2];
              rgb = [RR@kk,GG@kk,BB@kk] if inside(kk,RR) 
                    else [0.0,0.0,0.0];
              //rgb = RGB@kk if inside(kk,RGB) else [0.0,0.0,0.0];
           }
        }
        stabilize;
    }
}

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