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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2979 - (view) (download)

1 : glk 2979 // 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 = 8;
11 :     real isoval = 1;
12 :     image(2)[] diderot = image("data/hex.nrrd");
13 :     field#1(2)[] F = bspln3 ⊛ diderot;
14 :     input int stepsMax = 50;
15 :     input real stepScale = 1.0;
16 :     real epsilon = 0.0001;
17 :     //////////////// Global Variables for Spacing Particles //////////////////////
18 :     input real rr = 0.1; // actual particle radius
19 :     real stdev = ∞;
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 = 0;
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 :     contourSteps+=1;
63 :     if (contourSteps > iterMax) {
64 :     stabilize;
65 :     }
66 :     // compute energy and forces on us
67 :     energy = 0;
68 :     force = [0,0];
69 :     foreach (Particle p_j in sphere(RR)) {
70 :     if (p_j.foundContour) {
71 :     vec2 r_ij = (pos - p_j.pos)/rr;
72 :     vec2 d_ij = normalize(r_ij);
73 :     if (|r_ij| < 1) {
74 :     energy += (1 - |r_ij|)^4;
75 :     force += - (-4*(1 - |r_ij|)^3) * d_ij;
76 :     }
77 :     }
78 :     }
79 :     force /= rr; // smaller particles make larger forces
80 :     real stepMax = rr;
81 :     // update position based on force
82 :     posOld2 = posOld1; // shuffle saved positions down
83 :     posOld1 = pos;
84 :     if (energy > 0.0) { // we have neighbors
85 :     tensor[2,2] pten = identity[2] - normalize(∇F(pos))⊗normalize(∇F(pos));
86 :     force = pten•force; // project force onto tangent surface
87 :     vec2 step = hh*force;
88 :     if (|step| > stepMax) {
89 :     // decrease hh by factor by which step was too big
90 :     hh *= stepMax/|step|;
91 :     // and find smaller step
92 :     step = hh*force;
93 :     }
94 :     // take step and re-find implicit surface
95 :     pos += step;
96 :     pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos)); // Newton-Raphson step
97 :     pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos)); // Newton-Raphson step
98 :     real travel = |pos - posOld1| + |posOld1 - posOld2|;
99 :     if (travel > 0) {
100 :     // if we've moved in the past two steps, but we've moved back
101 :     // to where we were two steps ago, we're oscillating ==>
102 :     // okay = 0. Two steps in the same direction ==> okay = 1.
103 :     real okay = |pos - posOld2|/travel;
104 :     // slow down if oscillating, speed up a little if not
105 :     hh *= lerp(0.8, 1.001, 0, okay, 1);
106 :     }
107 :     } // if (energy > 0.0)
108 :     } // if (!foundContour)
109 :     }
110 :     }
111 :     global{
112 :     real energyMean = mean{P.energy | P in Particle.all};
113 :     real var = mean{ (P.energy - energyMean) * (P.energy - energyMean) | P in Particle.all};
114 :    
115 :     // real var = variance{P.energy | P in Particle.all};
116 :     stdev = sqrt(var);
117 :     }
118 :    
119 :     initially { Particle(ui + gridSize*vi,
120 :     [lerp(-1.5, 1.5, -0.5, real(ui), real(gridSize)-0.5),
121 :     lerp(-1.5, 1.5, -0.5, real(vi), real(gridSize)-0.5)])
122 :     | vi in 0..(gridSize-1), ui in 0..(gridSize-1) };

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