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

SCM Repository

[diderot] View of /tests/examples/isopt3d/isopt3d.diderot
ViewVC logotype

View of /tests/examples/isopt3d/isopt3d.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4697 - (download) (annotate)
Tue Oct 4 20:06:26 2016 UTC (2 years, 9 months ago) by glk
File size: 11737 byte(s)
shifting language from isosurface to feature
/* ==========================================
Mutually-repelling particles populating an isosurface

... under construction ...

This is heavily based on the [`sphere.diderot`](../sphere) example; to
understand the functioning of this program it is best to read through that
program and its comments first, as well as the [`circle.diderot`](../circle)
example it is based on.

The main change with this example is that instead of sampling a unit-sphere,
the particles are sampling an isosurface of a field. This means that there are
now two phases for each strand, controlled by bool featFound, initialized to
false. When `!featFound`, strands take Newton-Raphson steps towards the
feature, until they've converged within featEps, at which point
`featFound=true`. Then strands interact and control their population as in the
[`sphere.diderot`](../sphere) example, but with each movement, the feature
must be rediscovered with more Newton-Raphson steps. The initialization of the
system is different as well; instead of starting with a given list of initial
positions, particles are initialized on a grid that the user can set up to
encompass the region of interest.

========================================== */
input image(3)[] img ("data in which to find features") = image("img.nrrd");
input real isoval ("isovalue of isosurface to sample") = 0;

input vec3 cent ("center of region to be sampled") = [0,0,0];
input int sz0 ("# initial samples along X axis") = 10;
input int sz1 ("# initial samples along Y axis") = 10;
input int sz2 ("# initial samples along Z axis") = 10;
input real width ("extent along X axis around center to sample; extents along Y and Z are determined proportionally from sz0, sz1, sz2") = 2;
input real rad ("radius of particle's potential energy support") = 0.1;
input int featStepsMax ("max # steps allowed for initial convergence onto feature, or 0 to make unlimited") = 10;
input real featTravelMax ("maximum distance point is allowed to travel, as multiple of particle radius, during initial iterative search for feature") = 2;
input real featEps ("convergence threshold on initial feature search") = 0.00001;
input real geoEps ("geometric system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.05;
input int pcp ("periodicity of population control (if non-zero)") = 5;
input real hhInit ("initial step size for potential energy gradient descent") = 1;

real newDist = 0.5*rad; // how far away to birth new particles
real stepMax = rad;     // limit on distance to travel per iter
int iter = 0;           // which iteration we're on

// field is defined so that isocontour of interest is the zero levelset
field#1(3)[] F = bspln3 ⊛ clamp(img) - isoval;

// sampling extents along X, Y, Z
vec3 extent = width*[1, real(sz1)/sz0, real(sz2)/sz0];

// Univariate energy functions; see ../circle/circle.diderot for alternatives
function real  phi(real r) = (1 - r)^4;
function real phi'(real r) = -4*(1 - r)^3;
// Energy and force from particle at vec3 x
function real enr(vec3 x) = phi(|x|/rad);
function vec3 frc(vec3 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule

// Returns a non-zero vector perpendicular to given non-zero vector v
function vec3 perp3(vec3 v) {
   int c = 0;
   if (|v[0]| < |v[1]|) {
      c = 1;
   }
   // not v[c] because tensors can only be indexed by constants
   if (|v[1] if 1==c else v[0]| < |v[2]|) {
      c = 2;
   }
   // now c is index of v component with largest absolute value
   vec3 ret = ([v[1] - v[2], -v[0], v[0]] if (c == 0) else
               [-v[1], v[0] - v[2], v[1]] if (c == 1) else
               [-v[2], v[2], v[0] - v[1]]);
   return ret;
}

// From a given vec3, find a 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?
function bool pcIter() = (pcp > 0 && iter > 0 && 0 == iter % pcp);

function vec3 featureNorm(vec3 x) = normalize(∇F(x));
function vec3 featureStep(vec3 x) = -normalize(∇F(x)) * F(x)/|∇F(x)|;
function bool featureLost(vec3 x) = (|∇F(x)| == 0);

strand particle (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 ncount = 0;         // how many neighbors did we have
   bool featFound = false; // whether or not we've found the feature
   real featTravel = 0;    // distance traveled looking for feature
   int featSteps = 0;      // # steps taken looking for feature
   real undone = 1;        // same as in ../sphere
   update {
      /* Whether looking for the feature, or distributing points on it,
         strands not inside the field can't live */
      if (!inside(pos, F)) {
         die;
      }
      posLast = pos;
      if (false) print("hello from ", pos, " (featFound ", featFound, ")--------------\n");
      if (!featFound) {
         // various reasons to call it quits
         if (false) {
            print("   featSteps = ", featSteps, " (vs featStepsMax ", featStepsMax, ")\n");
            print("   featTravel = ", featTravel, " (vs featTravelMax ", featTravelMax, ")\n");
            print("   featureLost = ", featureLost(pos), "\n");
         }
         if ((featStepsMax > 0
                 && featSteps > featStepsMax)   // too many steps
             || (featTravelMax > 0
                 && featTravel > featTravelMax) // too much travel
             || featureLost(pos)) {             // can't compute step
            die;
         }
         if (false) print("   still here; F=", F(pos), "\n");
         // Take one step towards feature of interest
         step = featureStep(pos);
         pos += step;
         if (false) print("   step ", step, " (length ", |step|, ") to ", pos, " where F=", F(pos), "\n");
         // We've converged if step is small enough
         if (|step| < featEps) {
            featFound = true;
            if (false) print("   and we have featFound!\n");
         } else {
            featTravel += |step|/rad;
            featSteps += 1;
         }
      } else {
         // compute energy and forces on us from neighbors
         real energyLast = 0;
         vec3 force = [0,0,0];
         ncount = 0;
         foreach (particle P in sphere(rad)) {
            vec3 x = P.pos - pos;
            if (|x| == 0) { die; }
            energyLast += enr(x);
            force += frc(x);
            ncount += 1;
         }
         if (false) print(" (feat) at ", pos, " where F=", F(pos), ", ncount = ", ncount, "\n");
         vec3 norm = featureNorm(pos);
         if (0 == ncount) {
            if (pcIter()) { // no neighbors, so let's make one
               vec3 npos = pos + newDist*normalize(perp3(norm));
               if (false) print(" (feat) from pos ", pos, " where F=", F(pos), ", delF=", ∇F(pos), ": NEW at ", npos, "\n");
               new particle(npos, hh);
               undone = 1;
            } else {
               undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));
            }
            // set closest to something in case used in global update
            closest = newDist;
            continue;
         }

         /* Else we have interacting neighbors; project force onto
            tangent surface, find step, limit step size */
         if (false) print(" (feat) force = ", force);
         force = (identity[3] - norm⊗norm)•force;
         step = hh*force;
         if (false) print(" ---project--> ", force, "; step = ", step);
         if (|step| > stepMax) {
            hh *= stepMax/|step|;
            step = hh*force;
         }
         if (false) print(" ---limit--> ", step, "\n");

         // Take step, re-apply constraint, find new energy
         pos += step;
         pos += featureStep(pos);
         pos += featureStep(pos);
         pos += featureStep(pos);
         real energy = 0;
         closest = rad;
         ncount = 0;
         foreach (particle P in sphere(rad)) {
            energy += enr(P.pos - pos);
            closest = min(closest, |P.pos - pos|);
            ncount += 1;
         }
         if (false) print(" (feat) moved to pos ", pos, "; ncount = ", ncount, "\n");
         if (false) print(" (feat) energy - energyLast == ", energy, " - ", energyLast, " == ", energy - energyLast, "\n");

         // Armijo-Goldstein sufficient decrease condition
         if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
            hh *= 0.5;  // energy didn't decrease as expected: backtrack
            if (false) print("      OOPS backtrack!  hh == ", hh, "\n");
            pos = posLast;  // try again next iteration
            // no progress, so decrease of undone
         } else {
            hh *= 1.02; // opportunistically increase step size
            // indicate progress; may be over-written below
            undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));
            if (false) print("      progress  hh == ", hh, "; undone = ", undone, "; ncount = ", ncount, "; pcIter = ", pcIter(), "\n");
            // try to have between 5 and 8 neighbors
            if (pcIter()) {
               if (false) print("      considering pop cntl (iter ", iter, ", pcp ", pcp, ") with ncount ", ncount, "\n");
               if (ncount < 5) {
                  //                                                    
                  // new particle(pos + newDist*normalize(force), hh);
                  //                                                    
                  undone = 1;
               } else if (ncount > 8) {
                  if (posrnd(pos) < (ncount - 8.0)/ncount) {
                     die;
                  }
                  // else not done if too many neighbors, w/ popl control
                  undone = 1;
               }
            }
         }
      } // else featFound
      // Record actual step taken, in case used in global update
      step = pos - posLast;
   }
}

global {
   /* Compute coefficient-of-variation of distance to closest neighbor */
   real meancl = mean { P.closest | P in particle.all};
   real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
   real covcl = sqrt(varicl) / meancl;
   real meanncount = mean { real(P.ncount) | P in particle.all};
   real maxundone = max { P.undone | P in particle.all};
   bool allfound = (min { 1 if P.featFound else 0 | P in particle.all} == 1);
   if (false) print("\n\n\n\n");
   print("(iter ", iter, " w/ ", numStrands(), " or ", numActive(), ") mean(featfound) = ",
         mean { 1.0 if P.featFound else 0.0 | P in particle.all },
         "; allfound=", allfound,
         "; mean(hh)=", mean { P.hh | P in particle.all},
         "; COV(cl)=", covcl,
         "; mean(ncount)=", meanncount, "; max(undone)=", maxundone, "\n");
   if (false) print("\n\n\n\n");

   if (covcl < geoEps        // seem to be geometrically uniform
       && allfound           // and all found the feature
       && maxundone < 0.5) { // and no particle recently set undone=1
      print("Stabilizing ", numActive(), " points with COV(closest) ", covcl,
            " < ", geoEps, " and maxundone < ", maxundone, " < 0.5 (iter ", iter, ")\n");
      stabilize;
   }
   iter += 1;
}

initially { particle([cent[0] + lerp(-extent[0], extent[0], 0, ii, sz0-1)/2,
                      cent[1] + lerp(-extent[1], extent[1], 0, jj, sz1-1)/2,
                      cent[2] + lerp(-extent[2], extent[2], 0, kk, sz2-1)/2],
                     hhInit)
            | ii in 0 .. sz0-1,
              jj in 0 .. sz1-1,
              kk in 0 .. sz2-1 };

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