Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

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

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

Mon Mar 9 13:06:51 2015 UTC (4 years, 4 months ago) by lamonts
File size: 2555 byte(s)
`Made slight changes to iso2d.diderot, main.c`
```int gridSize = 200;
real isoval = 30;
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(∇F(pos))⊗normalize(∇F(pos))•force;
if (|force| > 0) {
// project force onto tangent plane
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) };
```