// 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 = 10; real isoval = 1; field#1(2)[] F = bspln3 ⊛ image("data/hex.nrrd"); input int stepsMax = 50; input real stepScale = 1.0; real epsilon = 0.0001; //////////////// Global Variables for Spacing Particles ////////////////////// input real rr = 0.4; // actual particle radius real stepMax = rr/2; real evariance = ∞; real RR = rr+0.0; // neighbor query radius (MUST be >= rr; should get same results for any RR >= rr) real hhInit = 10.0; // initial integration step size; can err too big, will be trimmed down during iterations int iterSpacing = 1; // which iteration we're on input int iterMax = 20; ///////////////////////////////////////////////////////////////////////////// strand Particle (int ID, vec2 pos0) { // world is 1x1 centered at (0.5, 0.5) vec2 pos = [pos0[0] + pos0[1]/gridSize, pos0[1] + pos0[0]/gridSize]; output vec2 outPos = [0,0]; bool foundContour = false; int contourSteps = 1; int pSteps = 0; real energy = 0; // HEY: can do convergence test based on variance of energies vec2 force = [0,0]; // or can test convergence based on sum of |force| vec2 posOld1 = pos; // remember last TWO positions vec2 posOld2 = pos; real hh = hhInit; stabilize { outPos = pos; } update { if (!foundContour) { // We bail if we're no longer inside or taken too many steps. if (!inside(pos, F) || pSteps > stepsMax) { die; } if (|∇F(pos)| == 0.0) { // can't compute step if |∇F|, so have to bail die; } vec2 delta = -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos)); // Newton-Raphson step if (|delta| < epsilon) { // we've converged if step is small enough foundContour = true; } else { pos += delta; pSteps += 1; } } else { // if (!foundContour) if (contourSteps > iterMax) { stabilize; } // GLK asks: can this computation of energy and force // be packaged up in a function or something else that // woule simplify re-computing energy and force at // a second, different position? e.g. after we do the // pos update, we want to test to see if we successfully // lowered energy energy = 0; force = [0,0]; foreach (Particle p_j in sphere(RR)) { if (p_j.foundContour) { vec2 r_ij = (pos - p_j.pos)/rr; vec2 d_ij = normalize(r_ij); if (|r_ij| < 1) { energy += (1 - |r_ij|)^4; force += - (-4*(1 - |r_ij|)^3) * d_ij; } } } force /= rr; // smaller particles make larger forces // update position based on force posOld2 = posOld1; // shuffle saved positions down posOld1 = pos; if (energy > 0.0) { // we have neighbors tensor[2,2] pten = identity[2] - normalize(∇F(pos))⊗normalize(∇F(pos)); force = pten•force; // project force onto tangent surface vec2 step = hh*force; if (|step| > stepMax) { // decrease hh by factor by which step was too big hh *= stepMax/|step|; // and find smaller step step = hh*force; } // take step and re-find implicit surface pos += step; pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos)); // Newton-Raphson step pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos)); // Newton-Raphson step real travel = |pos - posOld1| + |posOld1 - posOld2|; if (travel > 0) { // if we've moved in the past two steps, but we've moved back // to where we were two steps ago, we're oscillating ==> // okay = 0. Two steps in the same direction ==> okay = 1. real okay = |pos - posOld2|/travel; // slow down if oscillating, speed up a little if not hh *= lerp(0.8, 1.001, 0, okay, 1); } } // if (energy > 0.0) contourSteps+=1; } // if (!foundContour) } } global{ real energyMean = mean{P.energy | P in Particle.all}; evariance = mean{ (P.energy - energyMean) * (P.energy - energyMean) | P in Particle.all}; // GLK asks: what can we do with evariance? can we stabilize all strands // once variance goes below some threshold? // GLK asks: can we please make print() work from here? Or how else // can we learn how variance is changing between iterations? It it only // via the C API to the library-compiled-to-program? } initially { Particle(ui + gridSize*vi, [lerp(-1.5, 1.5, -0.5, real(ui), real(gridSize)-0.5), lerp(-1.5, 1.5, -0.5, real(vi), real(gridSize)-0.5)]) | 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: -0.5, real(vi), real(gridSize)-0.5)]) | vi in 0..(gridSize-1), ui in 0..(gridSize-1) };