input image(3)[] img ("data to analyze") = image("pbfs-vol.nrrd"); input real isoval ("isovalue of isosurface to sample") = 0; field#2(3)[] F = bspln3 ⊛ clamp(img); /* function vec3 featureStep(vec3 x) = (isoval - F(x))*∇F(x)/(∇F(x)•∇F(x)); 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}; } function bool featureLost(vec3 x) { if (!inside(x, F)) { return true; } real{3} L = evals(∇⊗∇F(x)); return -L{2}/(|∇F(x)| + 2) < 10; } 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 update { // can't proceed if outside field or have zero gradient if (featureLost(pos)) { die; } if (false) print("hello from ", pos, " (found ", found, ")--------------\n"); // HIDE 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"); stabilize; } iter += 1; } initially { point([cent[0] + width*lerp(-1, 1, 0, ii, sz-1)/2, cent[1] + width*lerp(-1, 1, 0, jj, sz-1)/2, cent[2] + width*lerp(-1, 1, 0, kk, sz-1)/2], hhInit) | ii in 0 .. sz-1, jj in 0 .. sz-1, kk in 0 .. sz-1 };
Click to toggle
does not end with </html> tag
does not end with </body> tag
The output has ended thus: Init) | ii in 0 .. sz-1, jj in 0 .. sz-1, kk in 0 .. sz-1 };