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
 [diderot] / tests / vis15-bugs / badptclC / pbfsmini.diderot

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

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