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

SCM Repository

[diderot] Annotation of /examples/iso2d-spatial/iso2d-glk.diderot
ViewVC logotype

Annotation of /examples/iso2d-spatial/iso2d-glk.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2980 - (view) (download)

1 : glk 2980 // iso2d
2 :     //
3 :     // Demo of finding isocontours via Newton-Raphson method.
4 :     // Initializes positions on a grid, and each update applies one
5 :     // step of Newton-Raphson.
6 :     //
7 :     // Process output with:
8 :     // unu jhisto -i iso2d.txt -b 512 512 -min 0 0 -max 1 1 | unu 2op neq - 0 | unu quantize -b 8 -o iso2d.png
9 :    
10 :     int gridSize = 10;
11 :     real isoval = 1;
12 :     field#1(2)[] F = bspln3 ⊛ image("data/hex.nrrd");
13 :     input int stepsMax = 50;
14 :     input real stepScale = 1.0;
15 :     real epsilon = 0.0001;
16 :     //////////////// Global Variables for Spacing Particles //////////////////////
17 :     input real rr = 0.4; // actual particle radius
18 :     real stepMax = rr/2;
19 :     real evariance = ∞;
20 :     real RR = rr+0.0; // neighbor query radius (MUST be >= rr; should get same results for any RR >= rr)
21 :     real hhInit = 10.0; // initial integration step size; can err too big, will be trimmed down during iterations
22 :     int iterSpacing = 1; // which iteration we're on
23 :     input int iterMax = 20;
24 :     /////////////////////////////////////////////////////////////////////////////
25 :    
26 :     strand Particle (int ID, vec2 pos0) {
27 :     // world is 1x1 centered at (0.5, 0.5)
28 :     vec2 pos = [pos0[0] + pos0[1]/gridSize, pos0[1] + pos0[0]/gridSize];
29 :     output vec2 outPos = [0,0];
30 :    
31 :     bool foundContour = false;
32 :     int contourSteps = 1;
33 :     int pSteps = 0;
34 :     real energy = 0; // HEY: can do convergence test based on variance of energies
35 :     vec2 force = [0,0]; // or can test convergence based on sum of |force|
36 :    
37 :     vec2 posOld1 = pos; // remember last TWO positions
38 :     vec2 posOld2 = pos;
39 :     real hh = hhInit;
40 :    
41 :     stabilize {
42 :     outPos = pos;
43 :     }
44 :    
45 :     update {
46 :     if (!foundContour) {
47 :     // We bail if we're no longer inside or taken too many steps.
48 :     if (!inside(pos, F) || pSteps > stepsMax) {
49 :     die;
50 :     }
51 :     if (|∇F(pos)| == 0.0) { // can't compute step if |∇F|, so have to bail
52 :     die;
53 :     }
54 :     vec2 delta = -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos)); // Newton-Raphson step
55 :     if (|delta| < epsilon) { // we've converged if step is small enough
56 :     foundContour = true;
57 :     } else {
58 :     pos += delta;
59 :     pSteps += 1;
60 :     }
61 :     } else { // if (!foundContour)
62 :     if (contourSteps > iterMax) {
63 :     stabilize;
64 :     }
65 :     // GLK asks: can this computation of energy and force
66 :     // be packaged up in a function or something else that
67 :     // woule simplify re-computing energy and force at
68 :     // a second, different position? e.g. after we do the
69 :     // pos update, we want to test to see if we successfully
70 :     // lowered energy
71 :     energy = 0;
72 :     force = [0,0];
73 :     foreach (Particle p_j in sphere(RR)) {
74 :     if (p_j.foundContour) {
75 :     vec2 r_ij = (pos - p_j.pos)/rr;
76 :     vec2 d_ij = normalize(r_ij);
77 :     if (|r_ij| < 1) {
78 :     energy += (1 - |r_ij|)^4;
79 :     force += - (-4*(1 - |r_ij|)^3) * d_ij;
80 :     }
81 :     }
82 :     }
83 :     force /= rr; // smaller particles make larger forces
84 :     // update position based on force
85 :     posOld2 = posOld1; // shuffle saved positions down
86 :     posOld1 = pos;
87 :     if (energy > 0.0) { // we have neighbors
88 :     tensor[2,2] pten = identity[2] - normalize(∇F(pos))⊗normalize(∇F(pos));
89 :     force = pten•force; // project force onto tangent surface
90 :     vec2 step = hh*force;
91 :     if (|step| > stepMax) {
92 :     // decrease hh by factor by which step was too big
93 :     hh *= stepMax/|step|;
94 :     // and find smaller step
95 :     step = hh*force;
96 :     }
97 :     // take step and re-find implicit surface
98 :     pos += step;
99 :     pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos)); // Newton-Raphson step
100 :     pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos)); // Newton-Raphson step
101 :     real travel = |pos - posOld1| + |posOld1 - posOld2|;
102 :     if (travel > 0) {
103 :     // if we've moved in the past two steps, but we've moved back
104 :     // to where we were two steps ago, we're oscillating ==>
105 :     // okay = 0. Two steps in the same direction ==> okay = 1.
106 :     real okay = |pos - posOld2|/travel;
107 :     // slow down if oscillating, speed up a little if not
108 :     hh *= lerp(0.8, 1.001, 0, okay, 1);
109 :     }
110 :     } // if (energy > 0.0)
111 :     contourSteps+=1;
112 :     } // if (!foundContour)
113 :     }
114 :     }
115 :     global{
116 :     real energyMean = mean{P.energy | P in Particle.all};
117 :     evariance = mean{ (P.energy - energyMean) * (P.energy - energyMean) | P in Particle.all};
118 :    
119 :     // GLK asks: what can we do with evariance? can we stabilize all strands
120 :     // once variance goes below some threshold?
121 :    
122 :     // GLK asks: can we please make print() work from here? Or how else
123 :     // can we learn how variance is changing between iterations? It it only
124 :     // via the C API to the library-compiled-to-program?
125 :     }
126 :    
127 :     initially { Particle(ui + gridSize*vi,
128 :     [lerp(-1.5, 1.5, -0.5, real(ui), real(gridSize)-0.5),
129 :     lerp(-1.5, 1.5, -0.5, real(vi), real(gridSize)-0.5)])
130 :     | vi in 0..(gridSize-1), ui in 0..(gridSize-1) };

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