// iso2d-cbrt
//
// Demo of difficulty of finding isocontours via Newton-Raphson method
// with nasty input; in this case it follows the shape of cbrt(r),
// so many of the points dramatically fail to converge.
//
// Process output with:
// unu jhisto -i mip.txt -b 500 500 -min 0 0 -max 499 499 | unu 2op neq - 0 | unu quantize -b 8 -o iso2d-cbrt.png
int gridSize = 150;
real isoval = 0.0;
field#2(2)[] F = image("../data/cbrt-500-unit.nrrd") ⊛ bspln3;
int stepsMax = 40;
real epsilon = 0.001;
strand sample (int ui, int vi) {
output vec2 pos = [lerp(1.0, 499.0, -0.5, real(ui), real(gridSize)-0.5),
lerp(1.0, 499.0, -0.5, real(vi), real(gridSize)-0.5)];
int steps = 0;
update {
// We bail if we're no longer inside or taken too many steps.
// We tag this strand result as garbage by setting it to a
// non-finite position; these will be filtered later.
// This is really what "die" is for, though (see below).
if (!inside(pos, F) || steps >= stepsMax) {
pos = [nan,nan];
// HEY (BUG) we shouldn't even have to use non-finite values;
// we should be able to just say "die;". Currently, using
// "die" leads to the program never finishing...
stabilize;
}
vec2 grad = ∇F(pos);
real gmsq = grad • grad;
// If |∇F|=0 we can't compute update; we have to bail
if (gmsq == 0.0) {
pos = [nan,nan];
stabilize;
}
// delta = Newton-Raphson step
vec2 delta = -(F(pos) - isoval)*grad/gmsq;
// We've converged if the update is small enough
if (|delta| < epsilon) {
stabilize;
}
pos = pos + delta;
steps = steps + 1;
}
}
initially [ sample(ui, vi) | vi in 0..(gridSize-1), ui in 0..(gridSize-1) ];