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

SCM Repository

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

Annotation of /tests/examples/isopt3d/isopt3d.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4697 - (view) (download)

1 : glk 4691 /* ==========================================
2 :     Mutually-repelling particles populating an isosurface
3 :    
4 :     ... under construction ...
5 :    
6 :     This is heavily based on the [`sphere.diderot`](../sphere) example; to
7 :     understand the functioning of this program it is best to read through that
8 : glk 4697 program and its comments first, as well as the [`circle.diderot`](../circle)
9 :     example it is based on.
10 : glk 4691
11 : glk 4697 The main change with this example is that instead of sampling a unit-sphere,
12 :     the particles are sampling an isosurface of a field. This means that there are
13 :     now two phases for each strand, controlled by bool featFound, initialized to
14 :     false. When `!featFound`, strands take Newton-Raphson steps towards the
15 :     feature, until they've converged within featEps, at which point
16 :     `featFound=true`. Then strands interact and control their population as in the
17 :     [`sphere.diderot`](../sphere) example, but with each movement, the feature
18 :     must be rediscovered with more Newton-Raphson steps. The initialization of the
19 :     system is different as well; instead of starting with a given list of initial
20 :     positions, particles are initialized on a grid that the user can set up to
21 :     encompass the region of interest.
22 :    
23 : glk 4691 ========================================== */
24 : glk 4697 input image(3)[] img ("data in which to find features") = image("img.nrrd");
25 : glk 4695 input real isoval ("isovalue of isosurface to sample") = 0;
26 :    
27 : glk 4691 input vec3 cent ("center of region to be sampled") = [0,0,0];
28 : glk 4695 input int sz0 ("# initial samples along X axis") = 10;
29 :     input int sz1 ("# initial samples along Y axis") = 10;
30 :     input int sz2 ("# initial samples along Z axis") = 10;
31 : glk 4691 input real width ("extent along X axis around center to sample; extents along Y and Z are determined proportionally from sz0, sz1, sz2") = 2;
32 :     input real rad ("radius of particle's potential energy support") = 0.1;
33 : glk 4697 input int featStepsMax ("max # steps allowed for initial convergence onto feature, or 0 to make unlimited") = 10;
34 :     input real featTravelMax ("maximum distance point is allowed to travel, as multiple of particle radius, during initial iterative search for feature") = 2;
35 :     input real featEps ("convergence threshold on initial feature search") = 0.00001;
36 : glk 4695 input real geoEps ("geometric system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.05;
37 : glk 4691 input int pcp ("periodicity of population control (if non-zero)") = 5;
38 :     input real hhInit ("initial step size for potential energy gradient descent") = 1;
39 :    
40 :     real newDist = 0.5*rad; // how far away to birth new particles
41 :     real stepMax = rad; // limit on distance to travel per iter
42 :     int iter = 0; // which iteration we're on
43 :    
44 : glk 4695 // field is defined so that isocontour of interest is the zero levelset
45 : glk 4697 field#1(3)[] F = bspln3 ⊛ clamp(img) - isoval;
46 : glk 4695
47 : glk 4691 // sampling extents along X, Y, Z
48 :     vec3 extent = width*[1, real(sz1)/sz0, real(sz2)/sz0];
49 :    
50 :     // Univariate energy functions; see ../circle/circle.diderot for alternatives
51 :     function real phi(real r) = (1 - r)^4;
52 :     function real phi'(real r) = -4*(1 - r)^3;
53 :     // Energy and force from particle at vec3 x
54 :     function real enr(vec3 x) = phi(|x|/rad);
55 :     function vec3 frc(vec3 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule
56 :    
57 :     // Returns a non-zero vector perpendicular to given non-zero vector v
58 :     function vec3 perp3(vec3 v) {
59 :     int c = 0;
60 :     if (|v[0]| < |v[1]|) {
61 :     c = 1;
62 :     }
63 :     // not v[c] because tensors can only be indexed by constants
64 :     if (|v[1] if 1==c else v[0]| < |v[2]|) {
65 :     c = 2;
66 :     }
67 :     // now c is index of v component with largest absolute value
68 :     vec3 ret = ([v[1] - v[2], -v[0], v[0]] if (c == 0) else
69 :     [-v[1], v[0] - v[2], v[1]] if (c == 1) else
70 :     [-v[2], v[2], v[0] - v[1]]);
71 :     return ret;
72 :     }
73 :    
74 :     // From a given vec3, find a random-ish value uniformly sampling [0,1)
75 :     function real posrnd(vec3 v) {
76 :     vec3 p = 10000*v/rad;
77 :     return fmod(|fmod(p[0],1)| + |fmod(p[1],1)| + |fmod(p[2],1)|, 1);
78 :     }
79 :    
80 :     // Is this an iteration in which to do population control?
81 :     function bool pcIter() = (pcp > 0 && iter > 0 && 0 == iter % pcp);
82 :    
83 : glk 4697 function vec3 featureNorm(vec3 x) = normalize(∇F(x));
84 : glk 4696 function vec3 featureStep(vec3 x) = -normalize(∇F(x)) * F(x)/|∇F(x)|;
85 : glk 4697 function bool featureLost(vec3 x) = (|∇F(x)| == 0);
86 : glk 4696
87 :     strand particle (vec3 pos0, real hh0) {
88 : glk 4691 output vec3 pos = pos0;
89 :     vec3 posLast = pos0;
90 :     real hh = hh0;
91 : glk 4697 vec3 step = [0,0,0]; // step towards feature or away from force
92 :     real closest = rad; // distance to closest neighbor
93 :     int ncount = 0; // how many neighbors did we have
94 :     bool featFound = false; // whether or not we've found the feature
95 :     real featTravel = 0; // distance traveled looking for feature
96 :     int featSteps = 0; // # steps taken looking for feature
97 :     real undone = 1; // same as in ../sphere
98 : glk 4691 update {
99 : glk 4697 /* Whether looking for the feature, or distributing points on it,
100 :     strands not inside the field can't live */
101 :     if (!inside(pos, F)) {
102 :     die;
103 :     }
104 : glk 4691 posLast = pos;
105 : glk 4697 if (false) print("hello from ", pos, " (featFound ", featFound, ")--------------\n");
106 :     if (!featFound) {
107 : glk 4695 // various reasons to call it quits
108 : glk 4696 if (false) {
109 : glk 4697 print(" featSteps = ", featSteps, " (vs featStepsMax ", featStepsMax, ")\n");
110 :     print(" featTravel = ", featTravel, " (vs featTravelMax ", featTravelMax, ")\n");
111 :     print(" featureLost = ", featureLost(pos), "\n");
112 : glk 4695 }
113 : glk 4697 if ((featStepsMax > 0
114 :     && featSteps > featStepsMax) // too many steps
115 :     || (featTravelMax > 0
116 :     && featTravel > featTravelMax) // too much travel
117 :     || featureLost(pos)) { // can't compute step
118 : glk 4695 die;
119 :     }
120 : glk 4696 if (false) print(" still here; F=", F(pos), "\n");
121 :     // Take one step towards feature of interest
122 :     step = featureStep(pos);
123 : glk 4695 pos += step;
124 : glk 4696 if (false) print(" step ", step, " (length ", |step|, ") to ", pos, " where F=", F(pos), "\n");
125 : glk 4695 // We've converged if step is small enough
126 : glk 4697 if (|step| < featEps) {
127 :     featFound = true;
128 :     if (false) print(" and we have featFound!\n");
129 : glk 4695 } else {
130 : glk 4697 featTravel += |step|/rad;
131 :     featSteps += 1;
132 : glk 4695 }
133 : glk 4691 } else {
134 :     // compute energy and forces on us from neighbors
135 :     real energyLast = 0;
136 :     vec3 force = [0,0,0];
137 :     ncount = 0;
138 :     foreach (particle P in sphere(rad)) {
139 :     vec3 x = P.pos - pos;
140 :     if (|x| == 0) { die; }
141 :     energyLast += enr(x);
142 :     force += frc(x);
143 :     ncount += 1;
144 :     }
145 : glk 4697 if (false) print(" (feat) at ", pos, " where F=", F(pos), ", ncount = ", ncount, "\n");
146 :     vec3 norm = featureNorm(pos);
147 : glk 4695 if (0 == ncount) {
148 :     if (pcIter()) { // no neighbors, so let's make one
149 :     vec3 npos = pos + newDist*normalize(perp3(norm));
150 : glk 4697 if (false) print(" (feat) from pos ", pos, " where F=", F(pos), ", delF=", ∇F(pos), ": NEW at ", npos, "\n");
151 : glk 4696 new particle(npos, hh);
152 : glk 4695 undone = 1;
153 :     } else {
154 :     undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));
155 :     }
156 : glk 4691 // set closest to something in case used in global update
157 : glk 4695 closest = newDist;
158 : glk 4691 continue;
159 :     }
160 :    
161 :     /* Else we have interacting neighbors; project force onto
162 :     tangent surface, find step, limit step size */
163 : glk 4697 if (false) print(" (feat) force = ", force);
164 : glk 4691 force = (identity[3] - norm⊗norm)•force;
165 :     step = hh*force;
166 : glk 4696 if (false) print(" ---project--> ", force, "; step = ", step);
167 : glk 4691 if (|step| > stepMax) {
168 :     hh *= stepMax/|step|;
169 :     step = hh*force;
170 :     }
171 : glk 4696 if (false) print(" ---limit--> ", step, "\n");
172 : glk 4691
173 :     // Take step, re-apply constraint, find new energy
174 : glk 4696 pos += step;
175 :     pos += featureStep(pos);
176 :     pos += featureStep(pos);
177 :     pos += featureStep(pos);
178 : glk 4691 real energy = 0;
179 :     closest = rad;
180 :     ncount = 0;
181 :     foreach (particle P in sphere(rad)) {
182 :     energy += enr(P.pos - pos);
183 :     closest = min(closest, |P.pos - pos|);
184 :     ncount += 1;
185 :     }
186 : glk 4697 if (false) print(" (feat) moved to pos ", pos, "; ncount = ", ncount, "\n");
187 :     if (false) print(" (feat) energy - energyLast == ", energy, " - ", energyLast, " == ", energy - energyLast, "\n");
188 : glk 4691
189 :     // Armijo-Goldstein sufficient decrease condition
190 : glk 4696 if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
191 : glk 4695 hh *= 0.5; // energy didn't decrease as expected: backtrack
192 : glk 4696 if (false) print(" OOPS backtrack! hh == ", hh, "\n");
193 : glk 4695 pos = posLast; // try again next iteration
194 :     // no progress, so decrease of undone
195 : glk 4691 } else {
196 : glk 4695 hh *= 1.02; // opportunistically increase step size
197 :     // indicate progress; may be over-written below
198 :     undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));
199 : glk 4696 if (false) print(" progress hh == ", hh, "; undone = ", undone, "; ncount = ", ncount, "; pcIter = ", pcIter(), "\n");
200 : glk 4691 // try to have between 5 and 8 neighbors
201 :     if (pcIter()) {
202 : glk 4696 if (false) print(" considering pop cntl (iter ", iter, ", pcp ", pcp, ") with ncount ", ncount, "\n");
203 : glk 4691 if (ncount < 5) {
204 : glk 4697 //
205 :     // new particle(pos + newDist*normalize(force), hh);
206 :     //
207 : glk 4695 undone = 1;
208 : glk 4691 } else if (ncount > 8) {
209 :     if (posrnd(pos) < (ncount - 8.0)/ncount) {
210 :     die;
211 :     }
212 : glk 4695 // else not done if too many neighbors, w/ popl control
213 :     undone = 1;
214 : glk 4691 }
215 :     }
216 :     }
217 : glk 4697 } // else featFound
218 : glk 4691 // Record actual step taken, in case used in global update
219 :     step = pos - posLast;
220 :     }
221 :     }
222 :    
223 :     global {
224 :     /* Compute coefficient-of-variation of distance to closest neighbor */
225 :     real meancl = mean { P.closest | P in particle.all};
226 :     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
227 :     real covcl = sqrt(varicl) / meancl;
228 : glk 4695 real meanncount = mean { real(P.ncount) | P in particle.all};
229 :     real maxundone = max { P.undone | P in particle.all};
230 : glk 4697 bool allfound = (min { 1 if P.featFound else 0 | P in particle.all} == 1);
231 : glk 4696 if (false) print("\n\n\n\n");
232 : glk 4697 print("(iter ", iter, " w/ ", numStrands(), " or ", numActive(), ") mean(featfound) = ",
233 :     mean { 1.0 if P.featFound else 0.0 | P in particle.all },
234 : glk 4695 "; allfound=", allfound,
235 :     "; mean(hh)=", mean { P.hh | P in particle.all},
236 :     "; COV(cl)=", covcl,
237 : glk 4696 "; mean(ncount)=", meanncount, "; max(undone)=", maxundone, "\n");
238 :     if (false) print("\n\n\n\n");
239 : glk 4691
240 : glk 4695 if (covcl < geoEps // seem to be geometrically uniform
241 : glk 4697 && allfound // and all found the feature
242 : glk 4695 && maxundone < 0.5) { // and no particle recently set undone=1
243 : glk 4691 print("Stabilizing ", numActive(), " points with COV(closest) ", covcl,
244 : glk 4695 " < ", geoEps, " and maxundone < ", maxundone, " < 0.5 (iter ", iter, ")\n");
245 : glk 4691 stabilize;
246 :     }
247 :     iter += 1;
248 :     }
249 :    
250 : glk 4696 initially { particle([cent[0] + lerp(-extent[0], extent[0], 0, ii, sz0-1)/2,
251 : glk 4691 cent[1] + lerp(-extent[1], extent[1], 0, jj, sz1-1)/2,
252 :     cent[2] + lerp(-extent[2], extent[2], 0, kk, sz2-1)/2],
253 :     hhInit)
254 :     | ii in 0 .. sz0-1,
255 :     jj in 0 .. sz1-1,
256 :     kk in 0 .. sz2-1 };

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