Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

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

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

Mon Mar 9 22:00:35 2015 UTC (4 years, 2 months ago) by glk
File size: 3173 byte(s)
`no longer need outPos`
```input int isoIterMax = 10;  // limit on isocontour search
input int forceIterMax = 0; // HEY get rid of this with global convergence criterion
input real epsilon = 0.001; // convergence criterion
input int res = 14;         // initialization sampling resolution
input real isoval = 1;      // isovalue
// scalar field in which we sample isocontour at isoval
field#1(3)[] F = bspln3 ⊛ image("data/hex3.nrrd");

real motion = ∞;            // mean normalized particle motion

function real phi(real x) = (1 - |x|)^4;
function real phi'(real x) = -4*(1 - |x|)^3;

strand Particle (int ID, vec3 pos0) {
output vec3 pos = pos0;       // particle position
vec3 dpos = [0,0,0];   // change in position
bool foundIso = false; // initial isocontour search done
int isoIter = 0;       // iteration of isocontour search
int forceIter = 0;  // HEY get rid of this with global convergence criterion
real tt = radius/2;    // line search step size

update {
if (!foundIso) {
if (!inside(pos, F) || isoIter == isoIterMax) {
die; // quit if outside field or took too many steps
}
// Newton-Raphson step
dpos = -normalize(∇F(pos)) * (F(pos) - isoval)/|∇F(pos)|;
pos += dpos;
isoIter += 1;
foundIso = true;
}
} else {
// if (0 == ID) print(ID, " tt=", tt, ", motion=", motion, "\n");
if (motion < epsilon) {
stabilize;
}
forceIter += 1;
if (forceIter == forceIterMax) {
stabilize;
}
real energy=0;
vec3 force=[0,0,0];
foreach (Particle P in sphere(radius)) {
vec3 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);
vec3 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 newenergy = 0;
// GLK asks: is there a way to write this summation with the
// same syntax as used in the global reductions below? e.g.:
if (newenergy > 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 };
}

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