Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] View of /examples/iso2d-spatial/iso2d-glk.diderot
 [diderot] / examples / iso2d-spatial / iso2d-glk.diderot # View of /examples/iso2d-spatial/iso2d-glk.diderot

Sat Mar 7 22:31:02 2015 UTC (4 years, 4 months ago) by glk
File size: 5603 byte(s)
`some questions answered`
```// 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.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)
// GLK asks: this says "should get same results for any RR >= rr",
// which is in fact true, but if you try, for example "RR = rr+0.3",
// then the results are NOT the same: the particles are not moving the
// same, and they don't converge correctly.  Why not?
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)
vec2 pos = [pos0 + pos0/gridSize, pos0 + pos0/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 > forceIterMax) {
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;
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 - 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: can we please make print() work from here?  Or how else
// can we learn how variance is changing between iterations?  Is it only
// via the C API to the library-compiled-to-program that we can learn
// this value at run-time?
}

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) };
```