// lzc2d-brute
//
// Demo of laplacian zero-crossing edges via brute-force search
// and then Illinois false position
//
// Process output with:
// unu jhisto -i pos.nrrd -b 512 512 -min 0 0 -max 1 1 | unu 2op neq - 0 | unu quantize -b 8  -o lzc2d.png

int gridSize = 300;
field#3(2)[] F = bspln5 ⊛ image("../data/ddro-80.nrrd");
int stepsMax = 30;
real epsilon = 0.00001;

strand sample (int ui, int vi) {
    output 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)];
    vec2 posA=pos;
    vec2 posB=pos;
    int stepsNum = 0;
    int side = 0;
    real step = 0.004;
    real fA=0.0;
    real fB=0.0;
    real A=0.0;
    real B=1.0;
    bool searching = true;
    update {
        // We bail if we're no longer inside or taken too many steps.
        if (!inside(pos, F) || stepsNum > stepsMax) {
	    die;
        }
        if (searching) {
          posA = pos;
          real ll = trace(∇⊗∇F(pos));
          // vec2 dir = normalize(∇F(pos));
          vec2 dir = normalize(identity[2]:∇⊗∇⊗∇F(pos));
          /* BUG: this third-derivative-based setting of dir
            (gradient-based setting works fine).  Compiling gives:
            high-il DFA: cpu = 0.000 seconds, gc = 0.000 seconds, 11 nodes, 10 visits, 1 iterations
            high-il DFA: cpu = 0.001 seconds, gc = 0.000 seconds, 48 nodes, 47 visits, 1 iterations
            high-il DFA: cpu = 0.004 seconds, gc = 0.000 seconds, 157 nodes, 156 visits, 1 iterations
            ***** Internal error after translation to HighIL: see log file for details
          */
          real sgn = 1.0 if ll > 0.0 else -1.0;
          posB = pos + sgn*step*dir;
          searching = trace(∇⊗∇F(posB))*ll > 0.0;
          if (!searching) {
             /* set up next phase */
             fA = trace(∇⊗∇F(posA));
             fB = trace(∇⊗∇F(posB));
          }
        } else {
          pos = lerp(posA, posB, fA, 0.0, fB);
          real fs = trace(∇⊗∇F(pos));
          if (0.0 == fs) {
            stabilize;
          }
          if (fs*fB > 0.0) { /* not between s and b */
            posB = pos;
            fB = fs;
            if (1 == side) {
              fA /= 2.0;
            }
            side = 1;
          } else { /* not between a and s */
            posA = pos;
            fA = fs;
            if (-1 == side) {
              fB /= 2.0;
            }
            side = -1;
          }
          if (|posA - posB| < epsilon ) {
            if (|∇F(pos)| < 180.0) {
              pos = [-1.0,-1.0];
            }
            stabilize;
          }
        }
        stepsNum += 1;
    }
}

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