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 3026 - (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 :     int gridSize = 10;
8 :     real isoval = 1;
9 :     field#1(2)[] F = bspln3 ⊛ image("data/hex.nrrd");
10 :     input int stepsMax = 50;
11 :     input real stepScale = 1.0;
12 :     real epsilon = 0.0001;
13 :     //////////////// Global Variables for Spacing Particles //////////////////////
14 : glk 2991 input real rr = 0.2; // actual particle radius
15 : glk 2980 real stepMax = rr/2;
16 :     real evariance = ∞;
17 : glk 3000 real RR = rr+0.0; // neighbor query radius (MUST be >= rr; should get same results for any RR >= rr)
18 : glk 2980 real hhInit = 10.0; // initial integration step size; can err too big, will be trimmed down during iterations
19 :     int iterSpacing = 1; // which iteration we're on
20 : glk 2991 input int forceIterMax = 20;
21 : glk 2980 /////////////////////////////////////////////////////////////////////////////
22 :    
23 :     strand Particle (int ID, vec2 pos0) {
24 :     // world is 1x1 centered at (0.5, 0.5)
25 : glk 3026 output vec2 pos = [pos0[0] + pos0[1]/gridSize, pos0[1] + pos0[0]/gridSize];
26 : glk 2980
27 :     bool foundContour = false;
28 :     int contourSteps = 1;
29 :     int pSteps = 0;
30 :     real energy = 0; // HEY: can do convergence test based on variance of energies
31 :     vec2 force = [0,0]; // or can test convergence based on sum of |force|
32 :    
33 :     vec2 posOld1 = pos; // remember last TWO positions
34 :     vec2 posOld2 = pos;
35 :     real hh = hhInit;
36 :    
37 :     update {
38 :     if (!foundContour) {
39 :     // We bail if we're no longer inside or taken too many steps.
40 :     if (!inside(pos, F) || pSteps > stepsMax) {
41 :     die;
42 :     }
43 :     if (|∇F(pos)| == 0.0) { // can't compute step if |∇F|, so have to bail
44 :     die;
45 :     }
46 :     vec2 delta = -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos)); // Newton-Raphson step
47 :     if (|delta| < epsilon) { // we've converged if step is small enough
48 :     foundContour = true;
49 :     } else {
50 :     pos += delta;
51 :     pSteps += 1;
52 :     }
53 :     } else { // if (!foundContour)
54 : glk 2991 if (contourSteps > forceIterMax) {
55 : glk 2980 stabilize;
56 :     }
57 :     energy = 0;
58 :     force = [0,0];
59 :     foreach (Particle p_j in sphere(RR)) {
60 :     if (p_j.foundContour) {
61 :     vec2 r_ij = (pos - p_j.pos)/rr;
62 :     if (|r_ij| < 1) {
63 :     energy += (1 - |r_ij|)^4;
64 : glk 2982 force += - (-4*(1 - |r_ij|)^3) * normalize(r_ij);
65 : glk 2980 }
66 :     }
67 :     }
68 :     force /= rr; // smaller particles make larger forces
69 :     // update position based on force
70 :     posOld2 = posOld1; // shuffle saved positions down
71 :     posOld1 = pos;
72 :     if (energy > 0.0) { // we have neighbors
73 :     tensor[2,2] pten = identity[2] - normalize(∇F(pos))⊗normalize(∇F(pos));
74 :     force = pten•force; // project force onto tangent surface
75 :     vec2 step = hh*force;
76 :     if (|step| > stepMax) {
77 :     // decrease hh by factor by which step was too big
78 :     hh *= stepMax/|step|;
79 :     // and find smaller step
80 :     step = hh*force;
81 :     }
82 :     // take step and re-find implicit surface
83 :     pos += step;
84 :     pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos)); // Newton-Raphson step
85 :     pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos)); // Newton-Raphson step
86 :     real travel = |pos - posOld1| + |posOld1 - posOld2|;
87 :     if (travel > 0) {
88 :     // if we've moved in the past two steps, but we've moved back
89 :     // to where we were two steps ago, we're oscillating ==>
90 :     // okay = 0. Two steps in the same direction ==> okay = 1.
91 :     real okay = |pos - posOld2|/travel;
92 :     // slow down if oscillating, speed up a little if not
93 :     hh *= lerp(0.8, 1.001, 0, okay, 1);
94 :     }
95 :     } // if (energy > 0.0)
96 :     contourSteps+=1;
97 :     } // if (!foundContour)
98 :     }
99 :     }
100 :     global{
101 :     real energyMean = mean{P.energy | P in Particle.all};
102 :     evariance = mean{ (P.energy - energyMean) * (P.energy - energyMean) | P in Particle.all};
103 : glk 3012 // printing here does work
104 :     //print("evariance = ", evariance, "\n");
105 : glk 2980 }
106 :    
107 :     initially { Particle(ui + gridSize*vi,
108 :     [lerp(-1.5, 1.5, -0.5, real(ui), real(gridSize)-0.5),
109 :     lerp(-1.5, 1.5, -0.5, real(vi), real(gridSize)-0.5)])
110 :     | vi in 0..(gridSize-1), ui in 0..(gridSize-1) };

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