// iso2d // // Demo of finding isocontours via Newton-Raphson method. // Initializes positions on a grid, and each update applies one // step of Newton-Raphson. 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.2; // 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 forceIterMax = 20; ///////////////////////////////////////////////////////////////////////////// strand Particle (int ID, vec2 pos0) { // world is 1x1 centered at (0.5, 0.5) output vec2 pos = [pos0[0] + pos0[1]/gridSize, pos0[1] + pos0[0]/gridSize]; 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; 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 > forceIterMax) { stabilize; } energy = 0; force = [0,0]; foreach (Particle p_j in sphere(RR)) { if (p_j.foundContour) { vec2 r_ij = (pos - p_j.pos)/rr; if (|r_ij| < 1) { energy += (1 - |r_ij|)^4; force += - (-4*(1 - |r_ij|)^3) * normalize(r_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}; // printing here does work //print("evariance = ", evariance, "\n"); } 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) };