// 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[2]; // 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[2]; 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) ];
Click to toggle
does not end with </html> tag
does not end with </body> tag
The output has ended thus: s = steps + 1; } } initially [ sample(ui, vi) | vi in 0..(gridSize-1), ui in 0..(gridSize-1) ];