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

SCM Repository

[diderot] View of /examples/iso2d-spatial/iso2d-glk4.diderot
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3031 - (download) (annotate)
Tue Mar 10 09:45:24 2015 UTC (4 years, 2 months ago) by glk
File size: 3287 byte(s)
now with basic population control
input real radius = 0.2;    // particle interaction radius
input real epsilon = 0.001; // convergence criterion
input int res = 5;         // initialization sampling resolution
input real isoval = 1;      // isovalue
real motion = ∞;            // mean particle motion in last iter
int giter = 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 point (int theID, 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
  int liter = 0;
  int ID = theID;

  update {
    liter += 1;
    if (!foundIso) {
      if (!inside(pos, F) || liter == 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;
      if (|dpos| < epsilon*radius) {
        foundIso = true;
        liter = 0;
      }
    } else {
      if (motion < epsilon)
        stabilize;
      if (liter == iterLimit)
        stabilize;
      real energy=0;
      vec2 force=[0,0];
      int ncount = 0;
      foreach (point P in sphere(radius)) {
        ncount += 1;
        if (true || P.foundIso) {
          vec2 r_ij = (pos - P.pos)/radius;
          energy += phi(|r_ij|);
          force -= normalize(r_ij)*phi'(|r_ij|)/radius;
        }
      }
      vec2 norm = -normalize(∇F(pos));
      if (ncount < 4 && (liter/10)*10 == liter) {
        vec2 npos = pos + 0.2*radius*[norm[1],-norm[0]];
        print(ID, ": birthing new at ", npos, "\n");
        new point(-1, npos);
      }
      // project force onto tangent plane
      force -= norm⊗norm•force;
      if (|force| > 0) { // take gradient descent step
        vec2 posLast = pos;
        pos += tt*normalize(force);
        // take Newton-Raphson steps back to surface
        pos -= normalize(∇F(pos)) * (F(pos) - isoval)/|∇F(pos)|;
        pos -= normalize(∇F(pos)) * (F(pos) - isoval)/|∇F(pos)|;
        if ((F(pos) - isoval)/|∇F(pos)| > epsilon*radius) {
          die; // this does occasionally happen
        }
        dpos = pos - posLast;
        real energyNew = 0;
        foreach (point P in sphere(radius))
          energyNew += phi(|pos - P.pos|/radius) if P.foundIso else 0.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 point.all };
  real pcount = sum{ 1.0 if P.foundIso else 0.0 | P in point.active };
  print(giter, ": count=", pcount, " motion=", motion, "\n");
  giter+=1;
}

initially { point(ui + res*vi,
                  [lerp(-0.3, 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) };

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