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

SCM Repository

[diderot] Diff of /tests/vis15-bugs/pbfs1.diderot
ViewVC logotype

Diff of /tests/vis15-bugs/pbfs1.diderot

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 5148, Fri Jul 21 17:11:33 2017 UTC revision 5149, Fri Jul 21 17:11:49 2017 UTC
# Line 1  Line 1 
1  input image(3)[] img ("data to analyze") = image("pbfs-vol.nrrd");  #version 1
2  input real isoval ("isovalue of isosurface to sample") = 0;  
3    image(3)[] img = image("pbfs-vol.nrrd");
4    
5  field#2(3)[] F = bspln3 ⊛ clamp(img);  field#2(3)[] F = bspln3 ⊛ clamp(img);
6    
7  /*  vec3 cent = [0,0,0];
8  function vec3 featureStep(vec3 x) =  int sz = 10;
9    (isoval - F(x))*∇F(x)/(∇F(x)•∇F(x));  real width = 2;
10  function tensor[3,3] featureProj(vec3 x) {  
   vec3 norm = normalize(∇F(x));  
   return identity[3] - norm⊗norm;  
 }  
 function bool featureLost(vec3 x) =  
   !inside(x, F) || |∇F(x)| == 0;  
 */  
   
 function vec3 featureStep(vec3 x) {  
   vec3{3} E = evecs(∇⊗∇F(x));  
   real{3} L = evals(∇⊗∇F(x));  
   return -(1/L{2})*E{2}⊗E{2}•∇F(x);  
 }  
 function tensor[3,3] featureProj(vec3 x) {  
   vec3{3} E = evecs(∇⊗∇F(x));  
   return E{0}⊗E{0} + E{1}⊗E{1};  
 }  
11  function bool featureLost(vec3 x) {  function bool featureLost(vec3 x) {
12    if (!inside(x, F)) { return true; }    if (!inside(x, F)) { return true; }
13    real{3} L = evals(∇⊗∇F(x));    return F(x) < 10;
   return -L{2}/(|∇F(x)| + 2) < 10;  
14  }  }
15    
16    strand point (vec3 pos0)
17    {
18      output real out = -1.0;
19      vec3 pos = pos0;
20    
 input vec3 cent ("center of region to initially sampled") = [0,0,0];  
 input int sz ("# initial samples along each axis") = 10;  
 input real width ("extent along each axis around center to sample") = 2;  
   
 input real tipd ("target inter-particle distance") = 0.1;  
 input real travMax ("max travel allowed to seek feature") = 2;  
 input int seekMax ("max # steps allowed to seek feature") = 10;  
 input real seekEps ("conv. threshold on feature seeking") = 0.00001;  
 input real mvmtEps ("conv. threshold on point movement") = 0.05;  
 input real geoEps ("conv. threshold on system geometry") = 0.05;  
 input int pcp ("periodicity of population control") = 5;  
 input real hhInit ("initial step size for energy descent") = 1;  
   
 /* Each particle wants between nnmm[0] and nnmm[1] neighbors */  
 input vec2 nnmm ("desired # neighbors min max") = [6.0, 8.0];  
   
 /* Potential function phi(r) has phi(0)=1, phi(r)=0 for r >= 1,  
    and a minima at phi'(0.666)=0 and phi(0.666)=-0.005, with C^3  
    continuity across the well and with 0 for r >= 1. The  
    potential well enables using inter-particle energy to create  
    the desired packing density */  
 function real phi(real r) =  
 (1 + r*(-4.248582222661734 + r*(5.991287388535527 + (-2.864766048887071 + 0.06790595504541913*r)*r)))  
 if r < 0.666 else  
 (-0.005 + ((0.4482053856358993 + (-2.6838645846460842 + (6.026642031390917 - 4.811690244623482*(-0.666 + r))*(-0.666 + r))*(-0.666 + r))*(-0.666 + r))*(-0.666 + r))  
 if r < 1 else 0;  
 // phi'(r) is d phi(r) / dr  
 function real phi'(real r) =  
 (-4.248582222661734 + r*(11.982574777071054 + (-8.594298146661213 + 0.2716238201816765*r)*r))  
 if r < 0.666 else  
 (-15.42591894092919 + 0.4482053856358993*(-1.332 + 2*r) + r*(-2.6838645846460842*(-3.996 + 3*r) + 6.026642031390917*(5.322672 + r*(-7.992 + 4*r)) - 4.811690244623482*(-5.90816592 + r*(13.30668 + r*(-13.32 + 5*r)))))  
 if r < 1 else 0;  
 real phiWellRad = 0.666;    // radius of potential well  
 real rad = tipd/phiWellRad; // actual radius of potential support  
 function real enr(vec3 x) = phi(|x|/rad); // energy at x  
 function vec3 frc(vec3 x) = // force on x (via chain rule)  
    phi'(|x|/rad) * (1/rad) * x/|x|;  
 real newDist = tipd;        // how far away to birth  
 real stepMax = rad;         // limit on travel per iter  
 int iter = 0;               // which iteration I'm on  
   
 // generate random-ish value uniformly sampling [0,1)  
 function real posrnd(vec3 v) {  
   vec3 p = 10000*v/rad;  
   return fmod(|fmod(p[0],1)| + |fmod(p[1],1)| + |fmod(p[2],1)|, 1);  
 }  
   
 /* Is this an iteration in which to do population control?  If  
    not, pcIter() returns 0. Otherwise, returns 1 when should  
    birth new particles, and -1 when should kill then off.  This  
    alternation is not due to any language limitations; just a  
    strategy for aiding the population control heuristics. */  
 function int pcIter() = ((iter/pcp)%2)*2 - 1  
                         if (pcp > 0 && iter > 0 && 0 == iter % pcp)  
                         else 0;  
   
 /* each strand seeks the feature (if !found), else interacts  
    with other strands to uniformly sample the feature */  
 strand point (vec3 pos0, real hh0) {  
   output vec3 pos = pos0;  
   vec3 posLast = pos0;  
   real hh = hh0;  
   vec3 step = [0,0,0]; // step towards feature or away from force  
   real closest = rad;  // distance to closest neighbor  
   int nn = 0;          // how many neighbors do I have  
   bool found = false;  // whether or not I've found the feature  
   real trav = 0;       // distance traveled looking for feature  
   int steps = 0;       // # steps taken looking for feature  
   real mvmt = 1;       // recent movement  
21    update {    update {
22      // can't proceed if outside field or have zero gradient      // can't proceed if outside field or have zero gradient
23      if (featureLost(pos)) { die; }      if (featureLost(pos)) { die; }
24      if (false) print("hello from ", pos, " (found ", found, ")--------------\n"); // HIDE      out = F (pos);
     posLast = pos;  
     if (!found) {  
       // various reasons to call it quits  
       if (false) { // HIDE  
         print("   steps = ", steps, " (vs seekMax ", seekMax, ")\n"); // HIDE  
         print("   trav = ", trav, " (vs travMax ", travMax, ")\n"); // HIDE  
         print("   featureLost = ", featureLost(pos), "\n"); // HIDE  
       } // HIDE  
       if ((seekMax > 0 && steps > seekMax)     // too many steps  
          || (travMax > 0 && trav > travMax)) { // too much travel  
         die;  
       }  
       if (false) print("   still here; F=", F(pos), "\n"); // HIDE  
       // Take one step towards feature of interest  
       step = featureStep(pos);  
       pos += step;  
       mvmt = |step|/rad;  
       if (false) print("  step ", step, " (length ", |step|, ") to ", pos, " where F=", F(pos), "\n"); // HIDE  
       // have converged if step is small enough  
       if (|step| < seekEps*rad) {  
         found = true;  
         if (false) print("   and I have found!\n"); // HIDE  
       } else {  
         trav += |step|;  
         steps += 1;  
       }  
     } else { // feature has been found  
       // find energy and forces on me from neighbors  
       real energyLast = 0;  
       vec3 force = [0,0,0];  
       nn = 0;  
       foreach (point P in sphere(rad)) {  
         vec3 d = P.pos - pos;  
         if (|d| == 0) { die; }  
         energyLast += enr(d);  
         force += frc(d);  
         nn += 1;  
       }  
       if (false) print(" (feat) at ", pos, " where F=", F(pos), ", nn = ", nn, "\n"); // HIDE  
       if (0 == nn) {  
         // wait until I have neighbors  
         continue;  
       }  
       /* Else I have interacting neighbors; project force onto  
         tangent surface, find step, limit step size */  
       if (false) print(" (feat) force = ", force);  
       force = featureProj(pos)•force;  
       step = hh*force;  
       if (false) print(" ---project--> ", force, "; step = ", step); // HIDE  
       if (|step| > stepMax) { // rescale hh to limit step taken  
         hh *= stepMax/|step|;  
         step = hh*force;  
       }  
       if (false) print(" ---limit--> ", step, "\n"); // HIDE  
   
       // Take step, re-apply constraint, find new energy  
       pos += step;  
       pos += featureStep(pos);  
       pos += featureStep(pos);  
       real energy = 0;  
       closest = rad;  
       nn = 0;  
       /* compute a mean neighbor offset in order to know (opposite)  
          direction in which to add new particles with pop. cntl. */  
       vec3 mno = [0,0,0];  
       foreach (point P in sphere(rad)) {  
         vec3 d = P.pos - pos;  
         energy += enr(d);  
         closest = min(closest, |d|);  
         mno += d;  
         nn += 1;  
       }  
       mno /= nn;  
       if (false) print(" (feat) moved to pos ", pos, "; nn = ", nn, "\n"); // HIDE  
       if (false) print(" (feat) energy - energyLast == ", energy, " - ", energyLast, " == ", energy - energyLast, "\n"); // HIDE  
   
       /* Armijo-Goldstein sufficient decrease condition: change  
          in energy should roughly match change predicted by force.  
          If not, backtrack and try again next iter with smaller hh */  
       if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {  
         hh *= 0.5;  
         if (false) print("      OOPS backtrack!  hh == ", hh, "\n"); //HIDE  
         pos = posLast;  // try again next iteration //HIDE  
       } else {  
         hh *= 1.1; // opportunistic step increase  
         if (false) print("      progress  hh == ", hh, "; nn = ", nn, "; pcIter = ", pcIter(), "\n"); // HIDE  
         if (pcIter() > 0) {  
           // try to have between nnmm[0] and nnmm[1] neighbors  
           if (false) print("      considering pop cntl (iter ", iter, ", pcp ", pcp, ") with nn ", nn, "\n"); // HIDE  
           /* the lower nn is compared to nnmm[0], the more likely  
              I'll try to birth a new particle */  
           if (energy < 0 && posrnd(pos) < (nnmm[0] - nn)/nnmm[0]) {  
             // stagger npos to avoid coincident particles */  
             vec3 noff = newDist*(posrnd(pos)/30 + 1)*normalize(-mno);  
             if (!featureLost(pos + noff)) { new point(pos + noff, hh); }  
           }  
         } else if (pcIter() < 0 && energy > 0 && nn > nnmm[1]) {  
           /* If this particle has nn neighbors, then all those  
              neighbors probably have a similar # neighbors. To  
              get down to nnmm[1] neighbors, all die with chance  
              of nn-nnmm[1] out of nn. */  
           if (posrnd(pos) < (nn - nnmm[1])/nn) {  
             die;  
           }  
         }  
       } // successfully moved to lower energy  
     } // else found  
     // Record actual step and movement  
     step = pos - posLast;  
     mvmt = |step|/rad;  
   }  
 }  
   
 global {  
   /* use coefficient-of-variation of distance to closest  
      neighbor as way of measuring geometric uniformity */  
   real meancl = mean { P.closest | P in point.all };  
   real varicl = mean { (P.closest - meancl)^2 | P in point.all };  
   real covcl = sqrt(varicl) / meancl;  
   bool allfound = (min { 1 if P.found else 0 | P in point.all} == 1);  
   real maxmv = max { P.mvmt | P in point.all };  
   print("(iter ", iter, " w/ ", numActive(), ")",  
         "; allfound=", allfound,  
         "; mean(hh)=", mean { P.hh | P in point.all}, //HIDE  
         "; mean(cl)=", meancl,  
         "; COV(cl)=", covcl,  
         "; max(mvmt)=", maxmv, "\n");  
   
   if (allfound              // all particles have found the feature  
       && covcl < geoEps     // is geometrically uniform  
       && maxmv < mvmtEps) { // nothing's moving much  
     print("Stabilizing ", numActive(), "(iter ", iter, "); COV(cl)=", covcl,  
         " < ", geoEps, "; max(mvmt)=", maxmv, " < ", mvmtEps, "\n");  
25        stabilize;        stabilize;
26    }    }
27    iter += 1;  
28  }  }
29    
30  initially { point([cent[0] + width*lerp(-1, 1, 0, ii, sz-1)/2,  initially { point([cent[0] + width*lerp(-1, 1, 0, ii, sz-1)/2,
31                     cent[1] + width*lerp(-1, 1, 0, jj, sz-1)/2,                     cent[1] + width*lerp(-1, 1, 0, jj, sz-1)/2,
32                     cent[2] + width*lerp(-1, 1, 0, kk, sz-1)/2],                     cent[2] + width*lerp(-1, 1, 0, kk, sz-1)/2])
                    hhInit)  
33              | ii in 0 .. sz-1,              | ii in 0 .. sz-1,
34                jj in 0 .. sz-1,                jj in 0 .. sz-1,
35                kk in 0 .. sz-1 };                kk in 0 .. sz-1 };

Legend:
Removed from v.5148  
changed lines
  Added in v.5149

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