int gridSize = 200; real isoval = 30; input real rad = 0.3; // actual particle radius real hhInit = rad/4; real evariance = ∞; input int newtonIterMax = 10; input int forceIterMax = 20; field#1(2)[] F = bspln3 ⊛ image("data/ddro-80.nrrd"); strand Particle (int ID, vec2 pos0) { output vec2 pos = pos0; bool foundContour = false; int iter = 0; int forceIter = 0; // want to get rid of this with global convergence criterion real energy = 0; // has to be a member variable so that global reduction can see it vec2 force = [0,0]; real hh = hhInit; update { if (!foundContour) { if (!inside(pos, F) || iter == newtonIterMax) { die; // quit if outside field or took too many steps } // Newton-Raphson step vec2 step = -normalize(∇F(pos))*(F(pos) - isoval)/|∇F(pos)|; pos += step; iter += 1; if (|step| < 0.001) { foundContour = true; } } else { if (forceIter == forceIterMax) { stabilize; } forceIter +=1; energy = 0; force = [0,0]; foreach (Particle p_j in sphere(rad)) { vec2 r_ij = (pos - p_j.pos)/rad; energy += (1 - |r_ij|)^4; force -= normalize(r_ij)*(-4*(1 - |r_ij|)^3)/rad; } force -= normalize(∇F(pos))⊗normalize(∇F(pos))•force; if (|force| > 0) { // project force onto tangent plane // take gradient descent step vec2 down = normalize(force); vec2 posLast = pos; pos += hh*down; // take Newton-Raphson step back to surface pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos)); pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos)); real newenergy = 0; foreach (Particle p_j2 in sphere(rad)) { vec2 r_ij = (pos - p_j2.pos)/rad; newenergy += (1 - |r_ij|)^4; } if (newenergy > energy - 0.5*hh*|force|) { hh *= 0.5; // backtrack pos = posLast; } else { // re-initialize for new particle neighborhood hh = hhInit; } } } // if (!foundContour) } } global{ real energyMean = mean{P.energy | P in Particle.all}; evariance = mean{ (P.energy - energyMean) * (P.energy - energyMean) | P in Particle.all}; } 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) };