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

SCM Repository

[diderot] View of /tests/vis15-bugs/badptclC/pbfsmini.diderot
ViewVC logotype

View of /tests/vis15-bugs/badptclC/pbfsmini.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5432 - (download) (annotate)
Wed Sep 27 17:28:52 2017 UTC (21 months, 2 weeks ago) by glk
File size: 6078 byte(s)
more printing
#version 1.0
input real rad ("inter-particle potential radius");
input real eps ("general convergence threshold");
input real v0 ("which isosurface to sample");
input vec3{} ipos ("initial point positions");
input image(3)[] vol ("data to analyze");
field#2(3)[] F = bspln3 ⊛ clamp(vol); // convolve w/ recon kernel

// Only these feature functions are specific to isosurfaces
function vec3 fStep(vec3 x) = (v0 - F(x))*∇F(x)/(∇F(x)•∇F(x));
function tensor[3,3] fPerp(vec3 x) {
  vec3 norm = normalize(∇F(x));
  return identity[3] - norm⊗norm;
}
function bool fMask(vec3 x) = |∇F(x)| > 0;
int iter = 0;

// Potential energy and force functions
function real phi(real r) = (1 - r)^4;
function real phi'(real r) = -4*(1 - r)^3;
function real enr(vec3 x) = phi(|x|/rad);
function vec3 frc(vec3 x) = phi'(|x|/rad) * (1/rad) * x/|x|;

function vec3 rndrgb(vec3 v) {
   vec3 p = (100000 + iter)*v/rad;
   return 4*[round(lerp(6, 63, |fmod(p[0],1)|)),
             round(lerp(6, 63, |fmod(p[1],1)|)),
             round(lerp(6, 63, |fmod(p[2],1)|))];
}

function vec2 jh(vec3 x) = [
   round(lerp(-0.5, 800-0.5, 0.2, x[0], 0.8)),
   round(lerp(-0.5, 800-0.5, 0.2, x[1], 0.8))];

function bool v3eq(vec3 a, vec3 b) = (a[0]==b[0] && a[1]==b[1] && a[2]==b[2]);

// Strands first find feature, then interact w/ or make neighbors
strand point (int idx, vec3 p0, real hh0, vec3 parentID) {
  output vec3 pos = p0;      // current particle position
  output vec3 ID = rndrgb(p0);
  real hh = hh0;             // energy gradient descent step size
  vec3 step = [0,0,0];       // total movement this iteration
  bool found = false;        // whether feature has been found
  int nstep = 0;             // # steps taken looking for feature
  int nback = 0;             // # times backtracked after found
  bool first = true;
  update {
    vec3 oldpos = pos;           // save previous position
    if (first) {
      print(ID, "(", iter, ") HELLO at pos=", pos, ":", jh(pos), " born of ", parentID, "\n");
      first = false;
    }
    if (!inside(pos, F) || !fMask(pos)) {
      print(ID, "(", iter, ") ", "" if found else "!", "found: inside=", inside(pos, F), "; fMask=", fMask(pos), ": DIE\n");
      die;
    }
    if (!found) {                // looking for feature
      step = fStep(pos);         // one step towards feature
      pos += step;
      print(ID, "(", iter, ")", " !found: step=", step, " to pos=", pos, ":", jh(pos), " where F=", F(pos), "\n");
      if (|step|/rad > eps) {    // if took a big step
        nstep += 1;
        if (nstep > 10) {
          print(ID, "(", iter, ")", " !found: I'm too slow; DIE\n");
          die; // too slow
        }
      } else {
        print(ID, "(", iter, ")", " !found: I've found it\n");
        found = true;
      }
    } else {      // feature found; interact with other points
      real oldE = 0;             // energy at current location
      vec3 force = [0,0,0];      // force on me from neighbors
      int nn = 0;                // number of neighbors
      foreach (point P in sphere(rad)) {
        oldE += enr(P.pos - pos);
        force += frc(P.pos - pos);
        print(ID, "(", iter, ")", ": pos=", pos, ":", jh(pos), " forced from ", P.ID, ":", jh(P.pos), " at distance ", |P.pos - pos|, "\n");
        nn += 1;
      }
      if (0 == nn) {             // no neighbors, so make one
        print(ID, "(", iter, ")", ": pos=", pos, ":", jh(pos), " but zero neighbors, NEW at pos=", pos + [0.5*rad,0,0], ":", jh(pos + [0.5*rad,0,0]), "\n");
        new point(-1, pos + [0.5*rad,0,0], hh, ID);
        print(ID, "(", iter, ")", " " if found else " !", "found: and iter done B\n");
        continue;
      }                          // else interact w/ neighbors
      vec3 away = hh*fPerp(pos)•force; // away from neighbors
      if (|away| > rad) {        // limit motion to radius
        hh *= rad/|away|;
        away = hh*force;
      }
      print(ID, "(", iter, ")", ": pos=", pos, ":", jh(pos), " had nn=", nn, "; oldE=", oldE, " hh=", hh, "\n");
      pos += away;               // take step along force
      print(ID, "(", iter, ")", "      after force step: pos=", pos, ":", jh(pos), "\n");
      pos += fStep(pos);         // move towards feature
      //real newE = sum { enr(P.pos - pos) | P in sphere(rad) };
      real newE = 0;
      foreach (point P in sphere(rad)) {
        newE += enr(P.pos - pos);
        print(ID, "(", iter, ")", ": pos=", pos, ":", jh(pos), " new energy from ", P.ID, ":", jh(P.pos), " at distance ", |P.pos - pos|, "\n");
      }
      print(ID, "(", iter, ")", "      after fStep: pos=", pos, ":", jh(pos), "; newE=", newE, "\n");
      if (newE - oldE < 0.5*(pos - oldpos)•(-force)) {
        hh *= 1.1;               // cautious step size increase
        nback = 0;               // reset backtrack counter
        if (nn < 5 && pos[0] > 0 && pos[1] > 0) {
          print(ID, "(", iter, ")", " pos=", pos, ":", jh(pos), "; nn=", nn, "; NEEW at pos=", pos + 0.5*rad*normalize(away), ":", jh(pos + 0.5*rad*normalize(away)), "\n");
          new point(-10 - nn, pos + 0.5*rad*normalize(away), hh, ID);
        }
      } else {
        pos = oldpos; nback += 1; // backtrack
        print(ID, "(", iter, ")", " newE|oldE=", newE, "|", oldE, " --> backtrack to pos=", pos, ":", jh(pos), "\n");
        if (nback > 10) {
           print(ID, "(", iter, ")", " got stuck DIE\n");
           die;
        }
        hh *= 0.5;                // try again w/ smaller step
      }
      step = pos - oldpos;        // move down force + to feature
      print(ID, "(", iter, ")", " " if found else " !", "found: and iter done C\n");
    } // else found
  } // update
}
global {
  bool allfound = all { P.found | P in point.all};
  real maxstep = max { |P.step| | P in point.all };
  print("================= done with iter ", iter, " w/ ", numActive(),
        "; allfound=", allfound,
        "; max(mvmt)=", maxstep,
        "\n");
  iter += 1;
  if (iter == 40) { stabilize; }
}
initially { point(ii, ipos[ii], 1, [0,0,0]) | ii in 0 .. length(ipos)-1 };

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