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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3005 - (download) (annotate)
Sat Mar 7 22:25:40 2015 UTC (4 years, 4 months ago) by glk
File size: 3438 byte(s)
latest results
input int isoIterMax = 10;  // limit on isocontour search
input int forceIterMax = 0; // HEY get rid of this with global convergence criterion
input real radius = 0.2;    // particle interaction radius
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) {
  vec3 pos = pos0;       // particle position
  // GLK asks: why can't we just say "output vec2 pos = pos0;"?
  // this leads to "Error: pos needs to be a strand state variable"
  // As is, have to make a separate outPos variable, set via stabilize
  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

  output vec3 outPos = [0,0,0];
  stabilize {
    outPos = pos;
  }

  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;
      if (|dpos|/radius < epsilon) {
        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|);
        force -= normalize(r_ij)*phi'(|r_ij|)/radius;
      }
      // 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.:
        // newenergy = sum{phi(|pos - P.pos|/radius) | Particle P in sphere(radius)}
        foreach (Particle P in sphere(radius))
          newenergy += phi(|pos - P.pos|/radius);
        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) };

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