Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] View of /tests/vis15-bugs/pbfs1.diderot
 [diderot] / tests / vis15-bugs / pbfs1.diderot # View of /tests/vis15-bugs/pbfs1.diderot

Fri Jul 21 05:54:36 2017 UTC (2 years ago) by glk
File size: 10492 byte(s)
`adding pbfs1 pbfs2`
```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 - 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 and nnmm 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;
function real enr(vec3 x) = phi(|x|/rad); // energy at x
function vec3 frc(vec3 x) = // force on x (via chain rule)
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) {
return fmod(|fmod(p,1)| + |fmod(p,1)| + |fmod(p,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;
if (false) print("  step ", step, " (length ", |step|, ") to ", pos, " where F=", F(pos), "\n"); // HIDE
// have converged if step is small enough
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;
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 and nnmm neighbors
if (false) print("      considering pop cntl (iter ", iter, ", pcp ", pcp, ") with nn ", nn, "\n"); // HIDE
/* the lower nn is compared to nnmm, the more likely
I'll try to birth a new particle */
if (energy < 0 && posrnd(pos) < (nnmm - nn)/nnmm) {
// 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) {
/* If this particle has nn neighbors, then all those
neighbors probably have a similar # neighbors. To
get down to nnmm neighbors, all die with chance
of nn-nnmm out of nn. */
if (posrnd(pos) < (nn - nnmm)/nn) {
die;
}
}
} // successfully moved to lower energy
} // else found
// Record actual step and movement
step = pos - posLast;
}
}

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 + width*lerp(-1, 1, 0, ii, sz-1)/2,
cent + width*lerp(-1, 1, 0, jj, sz-1)/2,
cent + 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 };
```