// iso2d // // Demo of finding isocontours via Newton-Raphson method. // Initializes positions on a grid, and each update applies one // step of Newton-Raphson. // // Process output with: // unu jhisto -i iso2d.txt -b 512 512 -min 0 0 -max 1 1 | unu 2op neq - 0 | unu quantize -b 8 -o iso2d.png int gridSize = 300; field#1(2)[] F = ctmr ⊛ image("../data/ddro-80.nrrd"); int stepsMax = 10; real epsilon = 0.0001; 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)]; // set the isvalue to 50, 30, or 10, depending on whichever we're closest to real isoval = 50.0 if F(pos) >= 40.0 else 30.0 if F(pos) >= 20.0 else 10.0; int steps = 0; update { // We bail if we're no longer inside or taken too many steps. if (!inside(pos, F) || steps > stepsMax) { die; } // GLKs recent changes (revision 1329) were made to make the code // a more obvious implementation of Newton-Raphson, and also to create // some obvious and not-so-obvious opportunties for optimization vec2 grad = ∇F(pos); if (|grad| == 0.0) { // can't compute step if |∇F|, so have to bail die; } vec2 norm = normalize(grad); vec2 delta = -((F(pos) - isoval)/|grad|)*norm; // Newton-Raphson step if (|delta| < epsilon) { // we've converged if step is small enough stabilize; } pos += delta; 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: steps += 1; } } initially { sample(ui, vi) | vi in 0..(gridSize-1), ui in 0..(gridSize-1) };