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

SCM Repository

[diderot] View of /tests/glk-tests/back/isosparse.diderot
ViewVC logotype

View of /tests/glk-tests/back/isosparse.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4640 - (download) (annotate)
Tue Sep 27 20:54:47 2016 UTC (3 years ago) by glk
File size: 4056 byte(s)
initial result of svn export --username anonsvn --password=anonsvn https://svn.smlnj-gforge.cs.uchicago.edu/svn/diderot/branches/vis15/src/tests/
// scalar field in which we sample isocontour
field#1(2)[] F = c4hexic ⊛ image("hand.nrrd");
field#0(2)[] D = tent ⊛ image("dist.nrrd");

input int size0 = 50;
input int size1 = 50;
int pnum = size0*size1;

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

// inter-particle energy, and its derivative
// WOW BUG?: if you change |x| to x, you get different results
function real  phi(real x) =
   (1 - |x|/radius)^4 if (|x| < radius) else 0.0;
function real phi'(real x) =
  -(4/radius)*(1 - |x|/radius)^3 if |x| < radius else 0.0;

strand point (int ID, int ui, int vi) {
  //output vec2 pos = porig + ui*dir0 + vi*dir1;   // particle position
  output vec2 pos = [lerp(15.590633, 106.409363, -0.5, ui, size0-0.5),
                     lerp(14.155009, 67.244995, -0.5, vi, size1-0.5)];
  int myID = ID;
  vec2 delta = [0,0];        // change in position
  bool foundIso = false;    // initial isocontour search done
  real hh = 1;
  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) > 1.4)
         die;
      // subsequent expressions are undefined if |∇F| is zero
      if (|∇F(pos)| == 0.0)
         die;
      // Newton-Raphson step
      delta = -normalize(∇F(pos)) * F(pos)/|∇F(pos)|;
      pos += delta;
      if (|delta|/radius < epsilon) {
        foundIso = true;
      }
    } else {
      if (!inside(pos, D))
        die;
      real energy=0;
      vec2 force=[0,0];
      //int ncount = 0;
      //int idmin = pnum;
      //real radmin = 2.0*radius;
      foreach (P in sphere(SRAD)) {
        vec2 r_ij = pos - P.pos;
        energy += phi(|r_ij|);
        force -= normalize(r_ij)*phi'(|r_ij|);
        //ncount += 1;
        //if (|r_ij| < radmin) {
        //  radmin = |r_ij|;
        //}
        //print(myID, ": P.myID = ", P.myID, ", idmin = ", idmin, " --> ");
        //if (P.myID < idmin) {
        //  idmin = P.myID;
        //}
        //print(idmin, "\n");
      }
      //if (ncount > 1 && radmin/radius < 0.2 && myID < idmin && iter > 10) { //  && iter == 10*(iter/10)) {
      //  die;
      //}
      // project force onto tangent plane
      force -= normalize(∇F(pos))⊗normalize(∇F(pos))•force;
      if (|force| > 0) { // take gradient descent step
        delta = hh*force;
        if (|delta| > limit) {
          //print(myID, ": knocking down hh ", hh, " by ", limit/|delta|, "\n");
          hh *= limit/|delta|;
          delta = hh*force;
        }
        vec2 posLast = pos;
        pos += delta;
        // take Newton-Raphson step back to surface
        if (!inside(pos, D)) die;
        pos -= normalize(∇F(pos)) * F(pos)/|∇F(pos)|;
        if (!inside(pos, D)) die;
        pos -= normalize(∇F(pos)) * F(pos)/|∇F(pos)|;
        if (!inside(pos, D)) die;
        pos -= normalize(∇F(pos)) * F(pos)/|∇F(pos)|;
        if (!inside(pos, D)) die;
        pos -= normalize(∇F(pos)) * F(pos)/|∇F(pos)|;
        if (!inside(pos, D)) die;
        real energyNew = 0;
        foreach (P in sphere(SRAD))
          energyNew += phi(|pos - P.pos|);
        if (energyNew > energy - 0.3*(pos - posLast)•force) {
          hh *= 0.5;  // backtrack
          pos = posLast;
        } else {
          hh *= 1.1; // bigger step on next iter
        }
      }
    }
  }
}
global{
  //real number = sum{ P.countme | P in point.all };
  //print("===================, ", number, "\n");
  real motion = mean{ |P.delta|/radius | P in point.all };
  if ((iter/100)*100 == iter) {
    print(iter, ": motion=", motion, "\n");
  }
  if (motion < epsilon/100)
    stabilize;
  iter+=1;
  if (iter == iterLimit)
    stabilize;
}

initially { 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