// FIXED [JHR; 23-03-2012]
// BUG:
// this program does not terminate. The issue is that the body of the
// "if (dir > 0.0)" is empty in the generated C code.
// The problem is that the SSA form constructed for the successor of the block
// is missing the PHI node, which makes the assignments in the block unused.
// It appears related to the fact that one arm of the if is a stabilize, which
// so there is a single-predecessor join node.
//
int imgSize = 200;
real h = 0.25; // step size of integration
int stepNum = 5; // take this many steps both upstream and downstream
field#1(2)[2] V = load("tor2slice-vec0.nrrd") ⊛ ctmr;
field#0(2)[] R = load("tor2slice-rand.nrrd") ⊛ tent;
strand LIC (int xi, int yi) {
real xx = lerp(0.0, 79.0, -0.5, real(xi), real(imgSize)-0.5);
real yy = lerp(0.0, 79.0, -0.5, real(yi), real(imgSize)-0.5);
vec2 pos0 = [xx,yy];
vec2 pos0backup = [xx,yy];
vec2 pos = pos0;
vec2 step = [0.0,0.0];
output vec3 out = [0.0,0.0,0.0];
real sum = 0.0;
int num = 0;
real dir = 1.0;
update {
step = h*dir*V(pos);
if (inside(pos + step, V) && inside(pos + step, R)) {
pos += step;
}
if (num == stepNum) {
if (dir > 0.0) {
num = 0;
pos = pos0backup;
dir = -1.0;
}
else
stabilize;
}
sum += R(pos);
num += 1;
}
stabilize {
out = [pos[0], pos[1], sum];
}
}
initially [ LIC(xi, yi) | yi in 0..(imgSize-1), xi in 0..(imgSize-1) ];