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

SCM Repository

[diderot] View of /tests/vis15-bugs/isohnd/isofixme.diderot
ViewVC logotype

View of /tests/vis15-bugs/isohnd/isofixme.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4640 - (download) (annotate)
Tue Sep 27 20:54:47 2016 UTC (2 years, 10 months ago) by glk
File size: 3634 byte(s)
initial result of svn export --username anonsvn --password=anonsvn https://svn.smlnj-gforge.cs.uchicago.edu/svn/diderot/branches/vis15/src/tests/
#version 2

// scalar field in which we sample isocontour
field#1(2)[] F = c4hexic ⊛ load_image("handchop.nrrd");
field#0(2)[] D = tent ⊛ load_image("densechop.nrrd");

input vec2 cent = [0,0];  // center of region initially sampled
input real hght = 2; // height or vertical FOV of sampling region
// (size0 and size1 determine the aspect ratio of the initially sampling region)
input int size0; // # horizontal samples
input int size1; // # vertical samples

real wdth = hght*size0/size1;

vec2 spc = [wdth/(size0-1), hght/(size1-1)];
vec2 dir0 = [spc[0], 0.0];
vec2 dir1 = [0.0, spc[1]];
vec2 orig = cent - (dir0*(size0-1) + dir1*(size1-1))/2.0;

input real radius = 2;    // particle interaction radius
input real epsilon = 0.001; // convergence criterion
int iter = 0;               // global iteration count
input int iterLimit = 0;    // upper limit on # iterations

// 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 ID, int ui, int vi) {
  int id = ID;
  output vec2 pos = orig + ui*dir0 + vi*dir1;   // particle position
  vec2 dpos = [0,0];        // change in position
  bool foundIso = false;    // initial isocontour search done
  real tt = radius/2;    // line search step size
  real countme = 1;

  update {
    if (!foundIso) {
      if (!inside(pos, F) || iter == 16) {
        die; // quit if outside field or took too many steps
      }
      if (!inside(pos, D) || D(pos) < 0.01) {
         die;
      }
      // subsequent expressions are undefined if |∇F| is zero
      if (|∇F(pos)| == 0.0) {
         die;
      }
      // Newton-Raphson step
      dpos = -normalize(∇F(pos)) * F(pos)/|∇F(pos)|;
      pos += dpos;
      if (|dpos|/radius < epsilon) {
        foundIso = true;
      }
    } else {
      if (!inside(pos, D) || D(pos) < 7) {
        die;
      }
      real energyLast=0;
      vec2 force=[0,0];
      foreach (point P in sphere(radius)) {
        vec2 r_ij = (pos - P.pos)/radius;
        energyLast += 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);
        vec2 posLast = pos;
        pos += dpos;
        // take Newton-Raphson step back to surface
        if (!inside(pos, F)) die;
        pos -= normalize(∇F(pos)) * F(pos)/|∇F(pos)|;
        if (!inside(pos, F)) die;
        pos -= normalize(∇F(pos)) * F(pos)/|∇F(pos)|;
        if (!inside(pos, F)) die;
        real energy = 0;
        foreach (point P in sphere(radius)) {
real d = |pos - P.pos|;
if (d > radius) { print("dist(", ID, " @ ", pos, ", ", P.id, " @ ", P.pos, ") = ", d, "\n"); }
          energy += phi(|pos - P.pos|/radius);
        }
        if (energy - energyLast > -0.5*(pos - posLast)•force) {
          tt *= 0.5; // backtrack
          pos = posLast;
        } else {
          tt = min(tt*1.1, radius/2); // bigger step on next iter
        }
      }
    }
  }
}

update {
print("***** ", iter, "\n");
  real motion = mean{ |P.dpos|/radius | P in point.all };
  if (motion < epsilon/100)
    stabilize;
  // BUG: try putting next line inside following conditional
  real number = sum{ P.countme | P in point.all };
  if (iter == iterLimit) {
    print("===================, ", number, "\n");
    stabilize;
  }
  iter+=1;
}

create_collection { point(ui + size0*vi, ui, vi) | vi in 0..(size1-1), ui in 0..(size0-1) }

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