Home My Page Projects Code Snippets Project Openings diderot
 Summary Activity Tracker Tasks SCM

# SCM Repository

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

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

revision 3012, Mon Mar 9 11:24:46 2015 UTC revision 3013, Mon Mar 9 13:06:51 2015 UTC
# Line 1  Line 1
1  // iso2d  int gridSize = 200;
2  //  real isoval = 30;
3  // Demo of finding isocontours via Newton-Raphson method.  input real rad = 0.3;      // actual particle radius
4  // Initializes positions on a grid, and each update applies one  real hhInit = rad/4;
5  // step of Newton-Raphson.  real evariance = ∞;
//
// 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 = 100;
image(2)[] diderot = image("data/ddro-80.nrrd");
field#1(2)[] F = ctmr ⊛ diderot;
input int stepsMax = 50;
input real stepScale = 1.0;
real epsilon = 0.0001;
//////////////// Global Variables for Spacing Particles //////////////////////
input real rr = 2.0;      // actual particle radius
real stdev = ∞;
real RR = rr+0.0;         // neighbor query radius (MUST be >= rr; should get same results for any RR >= rr)
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
int iterMax = 200;
real iosVal50 = 0.0;
real isoVal30 = 0.0;
real isoVal10 = 0.0;
/////////////////////////////////////////////////////////////////////////////

strand Particle (int ui, int vi) {
// world is 1x1 centered at (0.5, 0.5)
vec2 pos = [lerp(0.0, 1.0, -0.5, real(ui), real(gridSize)-0.5),
lerp(0.0, 1.0, -0.5, real(vi), real(gridSize)-0.5)];
// set the isvalue to 50, 30, or 10, depending on whichever we're closest to
real isoval = 50.0 if F(pos) >= 40.0
else 30.0 if F(pos) >= 20.0
else 10.0;

real pIsoVal50 = 0.0;
real pIsoVal30 = 0.0;
real pIsoVal10 = 0.0;
6
7      bool foundContour = false;  input int newtonIterMax = 10;
8      int contourSteps = 0;  input int forceIterMax = 20;
9      int pSteps = 0;
10      real energy = 0;     // HEY: can do convergence test based on variance of energies  field#1(2)[] F = bspln3 ⊛ image("data/ddro-80.nrrd");
11      vec2 force = [0,0];   // or can test convergence based on sum of |force|
12      output vec2 outPos = pos;  strand Particle (int ID, vec2 pos0) {
13      output vec2 pos = pos0;
14
15      vec2 posOld1 = pos;   // remember last TWO positions    bool foundContour = false;
16      vec2 posOld2 = pos;    int iter = 0;
17      int forceIter = 0;  // want to get rid of this with global convergence criterion
18      real energy = 0;    // has to be a member variable so that global reduction can see it
19      vec2 force = [0,0];
20      real hh = hhInit;      real hh = hhInit;
21
22      update {      update {
23        if (!foundContour) {
24          if(foundContour)        if (!inside(pos, F) || iter == newtonIterMax) {
25          {          die; // quit if outside field or took too many steps
26            if( contourSteps > iterMax)  {        }
27          // Newton-Raphson step
28          vec2 step = -normalize(∇F(pos))*(F(pos) - isoval)/|∇F(pos)|;
29          pos += step;
30          iter += 1;
31          if (|step| < 0.001) {
32            foundContour = true;
33          }
34        } else {
35          if (forceIter == forceIterMax) {
36              stabilize;              stabilize;
37            }            }
38          forceIter +=1;
39            // compute energy and forces on us        energy = 0; force = [0,0];
40            energy = 0;        foreach (Particle p_j in sphere(rad)) {
41            force = [0,0];          vec2 r_ij = (pos - p_j.pos)/rad;

foreach (Particle p_j in sphere(RR)) {

if(p_j.foundContour)
{
vec2 r_ij = (pos - p_j.pos)/rr;
vec2 d_ij = normalize(r_ij);
if (|r_ij| < 1) {
42                  energy += (1 - |r_ij|)^4;                  energy += (1 - |r_ij|)^4;
43                  force += - (-4*(1 - |r_ij|)^3) * d_ij;          force -= normalize(r_ij)*(-4*(1 - |r_ij|)^3)/rad;
44                }                }
45              }        force -= normalize(∇F(pos))⊗normalize(∇F(pos))•force;
46            }        if (|force| > 0) {
47            force /= rr;     // smaller particles make larger forces          // project force onto tangent plane
48            real stepMax = rr;  // iosVal50 if isoval == 5.0          // take gradient descent step
49                        // else isoVal10 if isoval == 30.0          vec2 down = normalize(force);
50                               //    else isoVal10;  // limit on distance to travel per iter          vec2 posLast = pos;
51            pos += hh*down;
52            // update position based on force          // take Newton-Raphson step back to surface
53            posOld2 = posOld1;   // shuffle saved positions down          pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos));
54            posOld1 = pos;          pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos));
55            if (energy > 0.0) {  // we have neighbors          real newenergy = 0;
56                tensor[2,2] pten = identity[2] - normalize(pos)⊗normalize(pos);          foreach (Particle p_j2 in sphere(rad)) {
57                force = pten•force;  // project force onto tangent surface            vec2 r_ij = (pos - p_j2.pos)/rad;
58                vec2 step = hh*force;            newenergy += (1 - |r_ij|)^4;
59                if (|step| > stepMax) {          }
60                    // decrease hh by factor by which step was too big          if (newenergy > energy - 0.5*hh*|force|) {
61                    hh *= stepMax/|step|;            hh *= 0.5; // backtrack
62                    // and find smaller step            pos = posLast;
step = hh*force;
}
// take step and re-find implicit surface
pos = normalize(pos + 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);
}
}
contourSteps+=1;
63          }else {          }else {
64              // We bail if we're no longer inside or taken too many steps.            // re-initialize for new particle neighborhood
65              if (!inside(pos, F) || pSteps > stepsMax) {            hh = hhInit;
die;
}
vec2 grad = ∇F(pos);
if (|grad| == 0.0) {    // can't compute step if |∇F|, so have to bail
die;
}
vec2 norm = normalize(grad);
vec2 delta = -((F(pos) - isoval)/|grad|)*norm;  // Newton-Raphson step
if (|delta| < epsilon) {    // we've converged if step is small enough
foundContour = true;
print("converged\n");
if(isoval == 50) {
pIsoVal50 = 1.0;
}else if (isoval == 30) {
pIsoVal30 = 1.0;
}else {
pIsoVal10 = 1.0;
66                   }                   }
67              }              }
68              pos += stepScale * delta;      } // if (!foundContour)
pSteps += 1;
}
outPos = pos;
69      }      }
70  }  }
71  global{  global{
72     real energyMean = mean{P.energy | P in Particle.all};     real energyMean = mean{P.energy | P in Particle.all};
73     real var = mean{ (P.energy - energyMean) * (P.energy - energyMean) | P in Particle.all};    evariance = mean{ (P.energy - energyMean) * (P.energy - energyMean) | P in Particle.all};
iosVal50 = sum{P.pIsoVal50 | P in Particle.all};
isoVal30 = sum{P.pIsoVal30 | P in Particle.all};
isoVal10 = sum{P.pIsoVal10 | P in Particle.all};

// real var = variance{P.energy | P in Particle.all};
stdev = sqrt(var);
74  }  }
75
initially { Particle(ui, vi) | vi in 0..(gridSize-1), ui in 0..(gridSize-1) };
76    initially { Particle(ui + gridSize*vi,
77                         [lerp(-1.5, 1.5, -0.5, real(ui), real(gridSize)-0.5),
78                          lerp(-1.5, 1.5, -0.5, real(vi), real(gridSize)-0.5)])
79                 | vi in 0..(gridSize-1), ui in 0..(gridSize-1) };

Legend:
 Removed from v.3012 changed lines Added in v.3013

 root@smlnj-gforge.cs.uchicago.edu ViewVC Help Powered by ViewVC 1.0.0