Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] View of /examples/iso2d-spatial/iso2d-glk3.diderot
 [diderot] / examples / iso2d-spatial / iso2d-glk3.diderot

# View of /examples/iso2d-spatial/iso2d-glk3.diderot

Mon Mar 9 21:54:11 2015 UTC (4 years, 2 months ago) by glk
File size: 2648 byte(s)
`good version, just before trying population control`
```input real radius = 0.1;    // particle interaction radius
input real epsilon = 0.001; // convergence criterion
input int res = 13;         // initialization sampling resolution
input real isoval = 1;      // isovalue
real motion = ∞;            // mean particle motion in last iter
int iter = 0;               // global iteration count
input int iterLimit = 0;    // upper limit on # iterations

// scalar field in which we sample isocontour at isoval
field#1(2)[] F = bspln3 ⊛ image("data/hex.nrrd");

// inter-particle energy, and its derivative
function real  phi(real x)  =   (1 - |x|)^4 if |x| < 1 else 0.0;
function real phi'(real x) = -4*(1 - |x|)^3 if |x| < 1 else 0.0;

strand Particle (vec2 pos0) {
output vec2 pos = pos0;   // particle position
vec2 dpos = [0,0];        // change in position
bool foundIso = false;    // initial isocontour search done
real tt = radius/2;    // line search step size

update {
if (!foundIso) {
if (!inside(pos, F) || iter == 10) {
die; // quit if outside field or took too many steps
}
// Newton-Raphson step
dpos = -normalize(∇F(pos)) * (F(pos) - isoval)/|∇F(pos)|;
pos += dpos;
foundIso = true;
}
} else {
if (motion < epsilon)
stabilize;
if (iter == iterLimit)
stabilize;
real energy=0;
vec2 force=[0,0];
foreach (Particle P in sphere(radius)) {
vec2 r_ij = (pos - P.pos)/radius;
energy += phi(|r_ij|);
}
// project force onto tangent plane
force -= normalize(∇F(pos))⊗normalize(∇F(pos))•force;
if (|force| > 0) { // take gradient descent step
dpos = tt*normalize(force);
vec2 posLast = pos;
pos += dpos;
// take Newton-Raphson step back to surface
pos -= normalize(∇F(pos)) * (F(pos) - isoval)/|∇F(pos)|;
pos -= normalize(∇F(pos)) * (F(pos) - isoval)/|∇F(pos)|;
real energyNew = 0;
if (energyNew > energy - 0.5*tt*|force|) {
tt *= 0.5; // backtrack
pos = posLast;
} else {
tt *= 2;   // bigger step for next iteration
}
}
}
}
}
global{
motion = mean{ |P.dpos|/radius | P in Particle.all };
print(iter, ": motion=", motion, "\n");
iter+=1;
}

initially { Particle([lerp(-1.5, 1.5, 0, ui, res-1),
lerp(-1.5, 1.5, 0, vi, res-1)])
| vi in 0..(res-1), ui in 0..(res-1) };
```