// lic.diderot
// 
// simple demo of the line integral convolution (LIC) method
// of vector field viz
//
// process output with:
// unu reshape -i mip.txt -s 300 200 | unu quantize -b 8  -o lic.png

int imgSizeX = 300;
int imgSizeY = 200;
real h = 0.01;    // step size of integration
int stepNum = 25; // take this many steps both upstream and downstream

//image(2)[2] Vimg = load("../data/vorttest.nrrd");
//field#1(2)[2] V = Vimg ⊛ ctmr;
//image(2)[] Rimg = load("../data/vorttest-rand.nrrd");
//field#0(2)[] R = Rimg ⊛ tent;

// HEY (BUG): the program compiles and run with these two field definitions
// but the results are NOT correct!!
// Correct result obtained with four-line field definition above
field#1(2)[2] V = load("../data/vorttest.nrrd") ⊛ ctmr;
field#0(2)[] R = load("../data/vorttest-rand.nrrd") ⊛ tent;


strand LIC (int xi, int yi) {
    real xx = lerp(0.0, 3.0, -0.5, real(xi), real(imgSizeX)-0.5);
    real yy = lerp(0.0, 2.0, -0.5, real(yi), real(imgSizeY)-0.5);
    vec2 pos0 = [xx,yy];
    vec2 forw = pos0;
    vec2 back = pos0;
    output real sum = R(pos0);
    int step = 0;

    update {
        // Euler integration step
        // forw = forw + h*V(forw);
        // back = back - h*V(back);
        // Midpoint method step
        forw = forw + h*V(forw + 0.5*h*V(forw));
        back = back - h*V(back - 0.5*h*V(back));
        sum = sum + R(forw) + R(back);
        step = step + 1;
        if (step == stepNum) {
          // modulate output by velocity at initial position
          sum = |V(pos0)|*sum/real(1 + 2*stepNum);
          stabilize;
        }
    }
}

initially [ LIC(xi, yi) | yi in 0..(imgSizeY-1), xi in 0..(imgSizeX-1) ];