// zslice-hess2
//
// Demo of field of tensor[2,2] from image of tensor[2,2]
//
// Process output with:
// unu reshape -s 2 200 200 -i zslice-hess2.txt | unu pad -min 0 0 0 -max 2 M M -b wrap | unu quantize -b 8  -o zslice-hess2.png

int gridSize = 200;
field#1(2)[2] G = ctmr ⊛ image("../data/ddro-80-grad.nrrd");


/* HEY: Not able to load nrrd as tensor[2,2] image; error is:

[zslice-hess2.diderot] Error: image file "../data/ddro-80-hess.nrrd" has range tensor[4], expected tensor[2,2]
uncaught exception Fail [Fail: Error in compiling zslice-hess2.diderot]
   raised at common/phase-timer.sml:76.50-76.52
   raised at common/phase-timer.sml:76.50-76.52
   raised at driver/main.sml:27.39-27.76

Can comment out following two lines, and below replace "hess = H(pos)"
with "hess = ∇⊗G(pos)", to get program to run and generate similar output
*/
image(2)[2,2] Himg = image("../data/ddro-80-hess.nrrd");
field#1(2)[2,2] H = ctmr ⊛ Himg;

strand sample (int ui, int vi) {
    vec2 pos = [lerp(0.0, 1.0, -0.5, real(ui), real(gridSize)-0.5),
                lerp(0.0, 1.0, -0.5, real(vi), real(gridSize)-0.5)];
    output vec2 val = [0.0,0.0];
    vec2 grad = [0.0,0.0];
    tensor[2,2] hess = [[0.0,0.0],[0.0,0.0]];
    update {
        grad = G(pos);
        //hess = ∇⊗G(pos);  // use this as long as tensor[2,2] fields not working
        hess = H(pos);    // use this when tensor[2,2] fields working
        val = hess • grad;
        stabilize;
    }
}

initially { sample(ui, vi) | vi in 0..(gridSize-1), ui in 0..(gridSize-1) };