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

SCM Repository

[diderot] Annotation of /tests/vis15-bugs/pbfs2.diderot
ViewVC logotype

Annotation of /tests/vis15-bugs/pbfs2.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5149 - (view) (download)

1 : jhr 5149 #version 1
2 :    
3 : glk 5146 input image(3)[] img ("data to analyze") = image("pbfs-vol.nrrd");
4 :     input real isoval ("isovalue of isosurface to sample") = 0;
5 :     input vec3{} ipos ("initial point positions") = load("pbfs-pos.txt");
6 :    
7 :     field#2(3)[] F = bspln3 ⊛ clamp(img);
8 :     //field#2(3)[] F = c4hexic ⊛ clamp(img);
9 :    
10 :     /*
11 :     function vec3 featureStep(vec3 x) =
12 :     (isoval - F(x))*∇F(x)/(∇F(x)•∇F(x));
13 :     function tensor[3,3] featureProj(vec3 x) {
14 :     vec3 norm = normalize(∇F(x));
15 :     return identity[3] - norm⊗norm;
16 :     }
17 :     function bool featureOut(vec3 x) =
18 :     |∇F(x)| == 0;
19 :     */
20 :     /*
21 :     function vec3 featureStep(vec3 x) {
22 :     vec3{3} E = evecs(∇⊗∇F(x));
23 :     real{3} L = evals(∇⊗∇F(x));
24 :     return -(1/L{2})*E{2}⊗E{2}•∇F(x);
25 :     }
26 :     function tensor[3,3] featureProj(vec3 x) {
27 :     vec3{3} E = evecs(∇⊗∇F(x));
28 :     return E{0}⊗E{0} + E{1}⊗E{1};
29 :     }
30 :     function bool featureOut(vec3 x) {
31 :     real{3} L = evals(∇⊗∇F(x));
32 :     return -L{2}/(|∇F(x)| + 2) < 10;
33 :     }
34 :     */
35 :     function vec3 featureStep(vec3 x) {
36 :     vec3{3} E = evecs(∇⊗∇F(x));
37 :     real{3} L = evals(∇⊗∇F(x));
38 :     return -(E{1}⊗E{1}/L{1} + E{2}⊗E{2}/L{2})•∇F(x);
39 :     }
40 :     function tensor[3,3] featureProj(vec3 x) {
41 :     vec3{3} E = evecs(∇⊗∇F(x));
42 :     return E{0}⊗E{0};
43 :     }
44 :     function bool featureOut(vec3 x) {
45 :     real{3} L = evals(∇⊗∇F(x));
46 :     return -L{1}/(∇F(x)•∇F(x) + 240) < 0.001 /* 0.1 */;
47 :     }
48 :    
49 :     input real tipd ("target inter-particle distance") = 0.1;
50 :     input real travMax ("max travel allowed to seek feature") = 2;
51 :     input int seekMax ("max # steps allowed to seek feature") = 10;
52 :     input real seekEps ("conv. threshold on feature seeking") = 0.001;
53 :     input real mvmtEps ("conv. threshold on point movement") = 0.01;
54 :     input real geoEps ("conv. threshold on system geometry") = 0.05;
55 :     input int pcp ("periodicity of population control (PC)") = 5;
56 :     input real hhInit ("initial step size for energy descent") = 1;
57 :     input bool verb ("verbose") = false; // HIDE
58 :    
59 :     /* Each particle wants between nnmm[0] and nnmm[1] neighbors */
60 :     input vec2 nnmm ("desired # neighbors min max") = [6.0, 8.0];
61 :    
62 :     /* Potential function phi(r) has phi(0)=1, phi(r)=0 for r >= 1,
63 :     and a minima at phi'(0.666)=0 and phi(0.666)=-0.005, with C^3
64 :     continuity across the well and with 0 for r >= 1. The
65 :     potential well enables using inter-particle energy to create
66 :     the desired packing density */
67 :     function real phi(real r) =
68 :     (1 + r*(-4.248582222661734 + r*(5.991287388535527 + (-2.864766048887071 + 0.06790595504541913*r)*r)))
69 :     if r < 0.666 else
70 :     (-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))
71 :     if r < 1 else 0;
72 :     // phi'(r) is d phi(r) / dr
73 :     function real phi'(real r) =
74 :     (-4.248582222661734 + r*(11.982574777071054 + (-8.594298146661213 + 0.2716238201816765*r)*r))
75 :     if r < 0.666 else
76 :     (-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)))))
77 :     if r < 1 else 0;
78 :     real phiWellRad = 0.666; // radius of potential well
79 :     real rad = tipd/phiWellRad; // actual radius of potential support
80 :     function real enr(vec3 x) = phi(|x|/rad); // energy at x
81 :     function vec3 frc(vec3 x) = // force on x (via chain rule)
82 :     phi'(|x|/rad) * (1/rad) * x/|x|;
83 :     real newDist = tipd; // how far away to birth
84 :     real stepMax = rad; // limit on travel per iter
85 :     int iter = 0; // current iteration of system
86 :     /* how much history matters in averaging movement; higher
87 :     means more memory and more stringent test of convergence */
88 :     real mvlerp = 0.5;
89 :    
90 :     /* a number in [0,1) roughly based on low 32 bits of significand
91 :     of given real x. Useful only when compiling with --double */
92 :     function real urnd(real x) {
93 :     if (x==0) return 0;
94 :     real l2 = log2(|x|);
95 :     real frxp = 2^(l2 - floor(l2) - 1); // like frexp(x)
96 :     real tmp = fmod(frxp*1048576, 1);
97 :     return fmod(2*tmp, 1);
98 :     // return fmod((frxp-0.5)*1048576,1); // 1048576==2^20; 2097152 == 2^21
99 :     }
100 :    
101 :     // generate random-ish value uniformly sampling [0,1)
102 :     function real posrnd(vec3 v) {
103 :     //real ret = fmod(urnd(v[0]) + urnd(v[1]) + urnd(v[2]), 1);
104 :     real ret = urnd(v[0]);
105 :     print("66666666 ", v[0], " ", v[1], " ", v[2], " ", ret, "\n");
106 :     return ret;
107 :     }
108 :    
109 :     /*
110 :     function int posid(vec3 v) {
111 :     vec3 p = 10000000*(v/rad);
112 :     return fmod(|fmod(p[0],1)| + |fmod(p[1],1)| + |fmod(p[2],1)|, 1);
113 :     }
114 :     */
115 :    
116 :     /* Is this an iteration in which to do population control (PC)?
117 :     If not, pcIter() returns 0. Otherwise, returns 1 when should
118 :     birth new particles, and -1 when should kill then off. This
119 :     alternation is not due to any language limitations; just a
120 :     strategy for aiding the PC heuristics. */
121 :     function int pcIter() = ((iter/pcp)%2)*2 - 1
122 :     if (pcp > 0 && iter > 0 && 0 == iter % pcp)
123 :     else 0;
124 :    
125 :     /* each strand seeks the feature (if !found), else interacts
126 :     with other strands to uniformly sample the feature */
127 :     strand point (vec3 pos0, real hh0, int pnn) {
128 :     output vec4 posn = [pos0[0], pos0[1], pos0[2], 0];
129 :     vec3 pos = pos0;
130 :     vec3 posLast = pos0; // previous position
131 :     real hh = hh0; // integration step size
132 :     vec3 step = [0,0,0]; // step towards feature or lower energy
133 :     real closest = rad; // distance to closest neighbor
134 :     int nn = 0; // how many neighbors do I have
135 :     bool found = false; // whether or not I've found the feature
136 :     real trav = 0; // distance traveled looking for feature
137 :     int steps = 0; // # steps taken looking for feature
138 :     real mvmt = 1; // recent movement avg, to test convergence
139 :     int fresh = 1; // 1: just created, else 0
140 :     update {
141 :     bool vv = false; // pos•[cos(1.658),sin(1.658),0] > 0.82; // HIDE
142 :     real tmp = posrnd(pos0);
143 :     stabilize;
144 :     if (verb && vv) print("hello from ", pos, " (found ", found, ")--------------------------------------\n"); // HIDE
145 :     // can't proceed if outside field or feature isn't there
146 :     if (!inside(pos, F) || featureOut(pos)) {
147 :     if (verb && vv) print("!inside || featureOut DIE\n"); // HIDE
148 :     die;
149 :     }
150 :     posLast = pos;
151 :     if (!found) {
152 :     // various reasons to call it quits
153 :     if ((seekMax > 0 && steps > seekMax) // too many steps
154 :     || (travMax > 0 && trav > travMax)) { // too much travel
155 :     if (1==pnn /* verb && vv */) { // HIDE
156 :     print(" steps = ", steps, " (vs seekMax ", seekMax, ")\n"); // HIDE
157 :     print(" trav = ", trav, " (vs travMax ", travMax, ")\n"); // HIDE
158 :     print(" featureOut = ", featureOut(pos), "\n"); // HIDE
159 :     } // HIDE
160 :     die;
161 :     }
162 :     if (verb && vv) print(" still here; F=", F(pos), "\n"); // HIDE
163 :     // Take one step towards feature of interest
164 :     step = featureStep(pos);
165 :     pos += step;
166 :     posn = [pos[0], pos[1], pos[2], 0];
167 :     if (verb && vv) print(" step ", step, " (length ", |step|, ") to ", pos, " where F=", F(pos), "\n"); // HIDE
168 :     // have converged if step is small enough
169 :     if (|step|/rad < seekEps) {
170 :     found = true;
171 :     if (verb && vv) print(" and I have found!\n"); // HIDE
172 :     } else {
173 :     trav += |step|;
174 :     steps += 1;
175 :     }
176 :     // leave mvmt at 1; only change with interactions below
177 :     } else { // feature has been found
178 :     // every path through following should update mvmt
179 :     real energyLast = 0; // energy at current location
180 :     vec3 force = [0,0,0]; // force on me from neighbors
181 :     nn = 0;
182 :     foreach (point P in sphere(rad)) {
183 :     /* interact with point P even when !P.found to
184 :     avoid disruption from change in P.found */
185 :     vec3 off = P.pos - pos;
186 :     //if (|off| == 0) { die; } // something has gone wrong
187 :     energyLast += enr(off);
188 :     force += frc(off);
189 :     nn += 1;
190 :     }
191 :     if (verb && vv) print(" (feat) at ", pos, " where F=", F(pos), ", nn = ", nn, "\n"); // HIDE
192 :     if (0 == nn) {
193 :     /* No neighbors. Can just wait with "continue", or, can
194 :     make a neighbor, but it takes some work to make sure
195 :     that the new position is well-defined for all
196 :     feature codimensions and orientations. */
197 :     vec3 noff0 = featureProj(pos)•[newDist,0,0];
198 :     vec3 noff1 = featureProj(pos)•[0,newDist,0];
199 :     vec3 noff2 = featureProj(pos)•[0,0,newDist];
200 :     vec3 noff = noff0;
201 :     noff = noff if |noff| > |noff1| else noff1;
202 :     noff = noff if |noff| > |noff2| else noff2;
203 :     if (verb && vv) print(" (feat) at ", pos, ": creating new at ", newDist*normalize(noff) + pos, "\n"); // HIDE
204 :     new point(newDist*normalize(noff) + pos, hh, nn);
205 :     mvmt = 1; // reset convergence
206 :     continue;
207 :     }
208 :     /* Else I have interacting neighbors */
209 :     if (verb && vv) print(" (feat) force = ", force); // HIDE
210 :     force = featureProj(pos)•force;
211 :     step = hh*force;
212 :     if (verb && vv) print(" ---project--> ", force, "; step = ", step); // HIDE
213 :     if (|step| > stepMax) { // rescale hh to limit step taken
214 :     hh *= stepMax/|step|;
215 :     step = hh*force;
216 :     }
217 :     if (verb && vv) print(" ---limit--> ", step, "\n"); // HIDE
218 :    
219 :     // Take step, re-apply constraint, find new energy
220 :     pos += step;
221 :     pos += featureStep(pos);
222 :     pos += featureStep(pos);
223 :     real energy = 0;
224 :     closest = rad;
225 :     nn = 0;
226 :     /* find mean neighbor offset to know (opposite) direction
227 :     in which to add new particles with PC */
228 :     vec3 mno = [0,0,0];
229 :     foreach (point P in sphere(rad)) {
230 :     vec3 off = P.pos - pos;
231 :     energy += enr(off);
232 :     closest = min(closest, |off|);
233 :     mno += off;
234 :     nn += 1;
235 :     }
236 :     mno /= nn;
237 :     posn = [pos[0], pos[1], pos[2], nn];
238 :     if (verb && vv) print(" (feat) moved to pos ", pos, "; nn = ", nn, "\n"); // HIDE
239 :     if (verb && vv) print(" (feat) energy - energyLast == ", energy, " - ", energyLast, " == ", energy - energyLast, "\n"); // HIDE
240 :    
241 :     /* Armijo-Goldstein sufficient decrease condition: change
242 :     in energy should roughly match change predicted by force.
243 :     If not, backtrack and try again next iter with smaller hh */
244 :     step = pos - posLast;
245 :     if (energy - energyLast > 0.5*step•(-force)) {
246 :     hh *= 0.5;
247 :     if (verb && vv) print(" OOPS backtrack! hh == ", hh, "\n"); //HIDE
248 :     pos = posLast; // try again next iteration
249 :     posn = [pos[0], pos[1], pos[2], nn];
250 :     } else {
251 :     hh *= 1.1; // opportunistic increase in step size
252 :     if (verb && vv) print("(", pos, ") progress hh == ", hh, "; nn = ", nn, "; pcIter = ", pcIter(), "\n"); // HIDE
253 :     }
254 :     if (|step|/rad < seekEps) {
255 :     /* even after backtracking, want to permit PC as long as
256 :     there is small (possibly tentative) motion */
257 :     if (verb && vv) print("(", pos, ") |step|/rad=", |step|/rad, " < ", seekEps, "=seekEps: maybe PC iter=", pcIter(), ", nn=", nn, ", energy=", energy, "\n"); // HIDE
258 :     if (pcIter() > 0 && energy < 0 && nn < nnmm[0]) {
259 :     /* the lower nn is compared to nnmm[0], the more likely
260 :     I'll try to birth a new particle */
261 :     vec3 npos = pos + newDist*normalize(-mno);
262 :     /* don't want new point being killed right away */
263 :     if (inside(npos, F) && !featureOut(npos)) {
264 :     if (verb && vv) print("(", pos, ") npos=", npos, "; in=", inside(npos,F), "; !fout=", !featureOut(npos), "; posrnd=", posrnd(pos), "\n"); // HIDE
265 :     if (posrnd(pos) < (nnmm[0] - nn)/nnmm[0]) {
266 :     new point(npos, hh, nn);
267 :     }
268 :     }
269 :     } else if (pcIter() < 0 && energy > 0 && nn > nnmm[1]) {
270 :     /* If this particle has nn neighbors, then all those
271 :     neighbors probably have a similar # neighbors. To
272 :     get down to nnmm[1] neighbors, all die with chance
273 :     of nn-nnmm[1] out of nn. */
274 :     if (posrnd(pos) < (nn - nnmm[1])/nn) {
275 :     if (1==pnn) { print("BBB pnn=1 dying\n"); }
276 :     die;
277 :     }
278 :     }
279 :     } // successfully moved to lower energy
280 :     // update movement history
281 :     mvmt = lerp(|pos - posLast|/rad, mvmt, mvlerp);
282 :     } // else found
283 :     } // update
284 :     }
285 :    
286 :     global {
287 :     /* use coefficient-of-variation of distance to closest
288 :     neighbor as way of measuring geometric uniformity */
289 :     real meancl = mean { P.closest | P in point.all };
290 :     real varicl = mean { (P.closest - meancl)^2 | P in point.all };
291 :     real covcl = sqrt(varicl) / meancl;
292 :     bool allfound = (min { 1 if P.found else 0 | P in point.all} == 1);
293 :     real maxmv = max { P.mvmt | P in point.all };
294 :     if (verb) { print("\n\n"); } // HIDE
295 :     print("(iter ", iter, " w/ ", numActive(), ")",
296 :     "; allfound=", allfound,
297 :     "; mean(hh)=", mean { P.hh | P in point.all}, //HIDE
298 :     "; mean(cl)=", meancl,
299 :     "; COV(cl)=", covcl,
300 :     "; range(nn)=[", min { P.nn | P in point.all }, ",",
301 :     max { P.nn | P in point.all }, "]",
302 :     "; max(mvmt)=", maxmv, "\n");
303 :     if (verb) { print("==============================================\n\n"); } // HIDE
304 :     if (allfound // all particles have found the feature
305 :     && covcl < geoEps // is geometrically uniform
306 :     && maxmv < mvmtEps) { // nothing's moving much
307 :     print("Stabilizing ", numActive(), " (iter ", iter, ")",
308 :     "; COV(cl)=", covcl, " < ", geoEps,
309 :     "; max(mvmt)=", maxmv, " < ", mvmtEps, "\n");
310 :     stabilize;
311 :     }
312 :     iter += 1;
313 :     }
314 :    
315 :     initially { point(ipos[ii], hhInit, 0) | ii in 0 .. length(ipos)-1 };

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