Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] Diff of /tests/examples/isopt3d/isopt3d.diderot
 [diderot] / tests / examples / isopt3d / isopt3d.diderot

# Diff of /tests/examples/isopt3d/isopt3d.diderot

revision 4696, Tue Oct 4 19:26:05 2016 UTC revision 4697, Tue Oct 4 20:06:26 2016 UTC
# Line 5  Line 5
5
6  This is heavily based on the [`sphere.diderot`](../sphere) example; to  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  understand the functioning of this program it is best to read through that
8  program and its comments first (as well as the [`circle.diderot`](../circle)  program and its comments first, as well as the [`circle.diderot`](../circle)
9  it based on). The main change with this example is that instead of sampling a  example it is based on.
10  unit-sphere, the particles are sampling an isosurface of a field.  The
11  initialization of the system is different as well; instead of starting with a  The main change with this example is that instead of sampling a unit-sphere,
12  given list of initial positions, particles are initialized on a grid.  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  ========================================== */  ========================================== */
24  input image(3)[] img ("data to isocontour") = image("img.nrrd");  input image(3)[] img ("data in which to find features") = image("img.nrrd");
25  input real isoval ("isovalue of isosurface to sample") = 0;  input real isoval ("isovalue of isosurface to sample") = 0;
26
27  input vec3 cent ("center of region to be sampled") = [0,0,0];  input vec3 cent ("center of region to be sampled") = [0,0,0];
# Line 21  Line 30
30  input int sz2 ("# initial samples along Z axis") = 10;  input int sz2 ("# initial samples along Z axis") = 10;
31  input real width ("extent along X axis around center to sample; extents along Y and Z are determined proportionally from sz0, sz1, sz2") = 2;  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;  input real rad ("radius of particle's potential energy support") = 0.1;
33  input int isoStepsMax ("max # steps allowed for initial convergence, or 0 to make unlimited") = 10;  input int featStepsMax ("max # steps allowed for initial convergence onto feature, or 0 to make unlimited") = 10;
34  input real isoTravelMax ("maximum distance point is allowed to travel, as multiple of particle radius, during initial iterative location of isocontour") = 2;  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 isoEps ("convergence threshold on Newton-based isosurface location") = 0.00001;  input real featEps ("convergence threshold on initial feature search") = 0.00001;
36  input real geoEps ("geometric system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.05;  input real geoEps ("geometric system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.05;
37  input int pcp ("periodicity of population control (if non-zero)") = 5;  input int pcp ("periodicity of population control (if non-zero)") = 5;
38  input real hhInit ("initial step size for potential energy gradient descent") = 1;  input real hhInit ("initial step size for potential energy gradient descent") = 1;
# Line 33  Line 42
42  int iter = 0;           // which iteration we're on  int iter = 0;           // which iteration we're on
43
44  // field is defined so that isocontour of interest is the zero levelset  // field is defined so that isocontour of interest is the zero levelset
45  field#1(3)[] F = c4hexic ⊛ clamp(img) - isoval;  field#1(3)[] F = bspln3 ⊛ clamp(img) - isoval;
46
47  // sampling extents along X, Y, Z  // sampling extents along X, Y, Z
48  vec3 extent = width*[1, real(sz1)/sz0, real(sz2)/sz0];  vec3 extent = width*[1, real(sz1)/sz0, real(sz2)/sz0];
# Line 71  Line 80
80  // Is this an iteration in which to do population control?  // Is this an iteration in which to do population control?
81  function bool pcIter() = (pcp > 0 && iter > 0 && 0 == iter % pcp);  function bool pcIter() = (pcp > 0 && iter > 0 && 0 == iter % pcp);
82
83    function vec3 featureNorm(vec3 x) = normalize(∇F(x));
84  function vec3 featureStep(vec3 x) = -normalize(∇F(x)) * F(x)/|∇F(x)|;  function vec3 featureStep(vec3 x) = -normalize(∇F(x)) * F(x)/|∇F(x)|;
85    function bool featureLost(vec3 x) = (|∇F(x)| == 0);
86
87  strand particle (vec3 pos0, real hh0) {  strand particle (vec3 pos0, real hh0) {
88     output vec3 pos = pos0;     output vec3 pos = pos0;
89     vec3 posLast = pos0;     vec3 posLast = pos0;
90     real hh = hh0;     real hh = hh0;
91     vec3 step = [0,0,0];   // step towards isosurface or away from force     vec3 step = [0,0,0];    // step towards feature or away from force
92     real closest = rad;    // distance to closest neighbor     real closest = rad;    // distance to closest neighbor
93     int ncount = 0;        // how many neighbors did we have     int ncount = 0;        // how many neighbors did we have
94     bool isoFound = false; // whether or not we've found the isosurface     bool featFound = false; // whether or not we've found the feature
95     real isoTravel = 0;    // distance traveled looking for isosurface     real featTravel = 0;    // distance traveled looking for feature
96     int isoSteps = 0;      // # steps taken looking for isosurface     int featSteps = 0;      // # steps taken looking for feature
97     real undone = 1;       // same as in ../sphere     real undone = 1;       // same as in ../sphere
98     update {     update {
99          /* 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        posLast = pos;        posLast = pos;
105        if (false) print("hello from ", pos, " (isoFound ", isoFound, ")--------------\n");        if (false) print("hello from ", pos, " (featFound ", featFound, ")--------------\n");
106        if (!isoFound) {        if (!featFound) {
107           // various reasons to call it quits           // various reasons to call it quits
108           if (false) {           if (false) {
109              print("   inside = ", inside(pos, F), "\n");              print("   featSteps = ", featSteps, " (vs featStepsMax ", featStepsMax, ")\n");
110              print("   isoSteps = ", isoSteps, " (vs isoStepsMax ", isoStepsMax, ")\n");              print("   featTravel = ", featTravel, " (vs featTravelMax ", featTravelMax, ")\n");
111              print("   isoTravel = ", isoTravel, " (vs isoTravelMax ", isoTravelMax, ")\n");              print("   featureLost = ", featureLost(pos), "\n");
112              print("   |del F| = ", |∇F(pos)|, "\n");           }
113           }           if ((featStepsMax > 0
114           if (!inside(pos, F)                  // no longer inside                   && featSteps > featStepsMax)   // too many steps
115               || (isoStepsMax > 0               || (featTravelMax > 0
116                   && isoSteps > isoStepsMax)   // too many steps                   && featTravel > featTravelMax) // too much travel
117               || (isoTravelMax > 0               || featureLost(pos)) {             // can't compute step
&& isoTravel > isoTravelMax) // too much travel
|| |∇F(pos)| == 0) {             // can't compute step
118              die;              die;
119           }           }
120           if (false) print("   still here; F=", F(pos), "\n");           if (false) print("   still here; F=", F(pos), "\n");
# Line 109  Line 123
123           pos += step;           pos += step;
124           if (false) print("   step ", step, " (length ", |step|, ") to ", pos, " where F=", F(pos), "\n");           if (false) print("   step ", step, " (length ", |step|, ") to ", pos, " where F=", F(pos), "\n");
125           // We've converged if step is small enough           // We've converged if step is small enough
126           if (|step| < isoEps) {           if (|step| < featEps) {
127              isoFound = true;              featFound = true;
128              if (false) print("   and we have isoFound!\n");              if (false) print("   and we have featFound!\n");
129           } else {           } else {
131              isoSteps += 1;              featSteps += 1;
132           }           }
133        } else {        } else {
134           // compute energy and forces on us from neighbors           // compute energy and forces on us from neighbors
# Line 128  Line 142
142              force += frc(x);              force += frc(x);
143              ncount += 1;              ncount += 1;
144           }           }
145           if (false) print(" (iso) at ", pos, " where F=", F(pos), ", ncount = ", ncount, "\n");           if (false) print(" (feat) at ", pos, " where F=", F(pos), ", ncount = ", ncount, "\n");
146           vec3 norm = normalize(∇F(pos));           vec3 norm = featureNorm(pos);
147           if (0 == ncount) {           if (0 == ncount) {
148              if (pcIter()) { // no neighbors, so let's make one              if (pcIter()) { // no neighbors, so let's make one
149                 vec3 npos = pos + newDist*normalize(perp3(norm));                 vec3 npos = pos + newDist*normalize(perp3(norm));
150                 if (false) print(" (iso) from pos ", pos, " where F=", F(pos), ", delF=", ∇F(pos), ": NEW at ", npos, "\n");                 if (false) print(" (feat) from pos ", pos, " where F=", F(pos), ", delF=", ∇F(pos), ": NEW at ", npos, "\n");
151                 new particle(npos, hh);                 new particle(npos, hh);
152                 undone = 1;                 undone = 1;
153              } else {              } else {
# Line 146  Line 160
160
161           /* Else we have interacting neighbors; project force onto           /* Else we have interacting neighbors; project force onto
162              tangent surface, find step, limit step size */              tangent surface, find step, limit step size */
163           if (false) print(" (iso) force = ", force);           if (false) print(" (feat) force = ", force);
164           force = (identity[3] - norm⊗norm)•force;           force = (identity[3] - norm⊗norm)•force;
165           step = hh*force;           step = hh*force;
166           if (false) print(" ---project--> ", force, "; step = ", step);           if (false) print(" ---project--> ", force, "; step = ", step);
# Line 169  Line 183
183              closest = min(closest, |P.pos - pos|);              closest = min(closest, |P.pos - pos|);
184              ncount += 1;              ncount += 1;
185           }           }
186           if (false) print(" (iso) moved to pos ", pos, "; ncount = ", ncount, "\n");           if (false) print(" (feat) moved to pos ", pos, "; ncount = ", ncount, "\n");
187           if (false) print(" (iso) energy - energyLast == ", energy, " - ", energyLast, " == ", energy - energyLast, "\n");           if (false) print(" (feat) energy - energyLast == ", energy, " - ", energyLast, " == ", energy - energyLast, "\n");
188
189           // Armijo-Goldstein sufficient decrease condition           // Armijo-Goldstein sufficient decrease condition
190           if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {           if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
# Line 187  Line 201
201              if (pcIter()) {              if (pcIter()) {
202                 if (false) print("      considering pop cntl (iter ", iter, ", pcp ", pcp, ") with ncount ", ncount, "\n");                 if (false) print("      considering pop cntl (iter ", iter, ", pcp ", pcp, ") with ncount ", ncount, "\n");
203                 if (ncount < 5) {                 if (ncount < 5) {
204                    new particle(pos + newDist*normalize(force), hh);                    //
205                      // new particle(pos + newDist*normalize(force), hh);
206                      //
207                    undone = 1;                    undone = 1;
208                 } else if (ncount > 8) {                 } else if (ncount > 8) {
209                    if (posrnd(pos) < (ncount - 8.0)/ncount) {                    if (posrnd(pos) < (ncount - 8.0)/ncount) {
# Line 198  Line 214
214                 }                 }
215              }              }
216           }           }
217        } // else isoFound        } // else featFound
218        // Record actual step taken, in case used in global update        // Record actual step taken, in case used in global update
219        step = pos - posLast;        step = pos - posLast;
220     }     }
# Line 211  Line 227
227     real covcl = sqrt(varicl) / meancl;     real covcl = sqrt(varicl) / meancl;
228     real meanncount = mean { real(P.ncount) | P in particle.all};     real meanncount = mean { real(P.ncount) | P in particle.all};
229     real maxundone = max { P.undone | P in particle.all};     real maxundone = max { P.undone | P in particle.all};
230     bool allfound = (min { 1 if P.isoFound else 0 | P in particle.all} == 1);     bool allfound = (min { 1 if P.featFound else 0 | P in particle.all} == 1);
231     if (false) print("\n\n\n\n");     if (false) print("\n\n\n\n");
232     print("(iter ", iter, " w/ ", numStrands(), " or ", numActive(), ") mean(isofound) = ",     print("(iter ", iter, " w/ ", numStrands(), " or ", numActive(), ") mean(featfound) = ",
233           mean { 1.0 if P.isoFound else 0.0 | P in particle.all },           mean { 1.0 if P.featFound else 0.0 | P in particle.all },
234           "; allfound=", allfound,           "; allfound=", allfound,
235           "; mean(hh)=", mean { P.hh | P in particle.all},           "; mean(hh)=", mean { P.hh | P in particle.all},
236           "; COV(cl)=", covcl,           "; COV(cl)=", covcl,
# Line 222  Line 238
238     if (false) print("\n\n\n\n");     if (false) print("\n\n\n\n");
239
240     if (covcl < geoEps        // seem to be geometrically uniform     if (covcl < geoEps        // seem to be geometrically uniform
241         && allfound           // and all found the isosurface         && allfound           // and all found the feature
242         && maxundone < 0.5) { // and no particle recently set undone=1         && maxundone < 0.5) { // and no particle recently set undone=1
243        print("Stabilizing ", numActive(), " points with COV(closest) ", covcl,        print("Stabilizing ", numActive(), " points with COV(closest) ", covcl,
244              " < ", geoEps, " and maxundone < ", maxundone, " < 0.5 (iter ", iter, ")\n");              " < ", geoEps, " and maxundone < ", maxundone, " < 0.5 (iter ", iter, ")\n");

Legend:
 Removed from v.4696 changed lines Added in v.4697