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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3013 - (view) (download)

1 : lamonts 3013 int gridSize = 200;
2 :     real isoval = 30;
3 :     input real rad = 0.3; // actual particle radius
4 :     real hhInit = rad/4;
5 :     real evariance = ∞;
6 : lamonts 2950
7 : lamonts 3013 input int newtonIterMax = 10;
8 :     input int forceIterMax = 20;
9 : lamonts 2950
10 : lamonts 3013 field#1(2)[] F = bspln3 ⊛ image("data/ddro-80.nrrd");
11 : lamonts 2950
12 : lamonts 3013 strand Particle (int ID, vec2 pos0) {
13 :     output vec2 pos = pos0;
14 : lamonts 2950
15 : lamonts 3013 bool foundContour = false;
16 :     int iter = 0;
17 :     int forceIter = 0; // want to get rid of this with global convergence criterion
18 :     real energy = 0; // has to be a member variable so that global reduction can see it
19 :     vec2 force = [0,0];
20 :     real hh = hhInit;
21 : lamonts 2950
22 : lamonts 3013 update {
23 :     if (!foundContour) {
24 :     if (!inside(pos, F) || iter == newtonIterMax) {
25 :     die; // quit if outside field or took too many steps
26 :     }
27 :     // Newton-Raphson step
28 :     vec2 step = -normalize(∇F(pos))*(F(pos) - isoval)/|∇F(pos)|;
29 :     pos += step;
30 :     iter += 1;
31 :     if (|step| < 0.001) {
32 :     foundContour = true;
33 :     }
34 :     } else {
35 :     if (forceIter == forceIterMax) {
36 :     stabilize;
37 :     }
38 :     forceIter +=1;
39 :     energy = 0; force = [0,0];
40 :     foreach (Particle p_j in sphere(rad)) {
41 :     vec2 r_ij = (pos - p_j.pos)/rad;
42 :     energy += (1 - |r_ij|)^4;
43 :     force -= normalize(r_ij)*(-4*(1 - |r_ij|)^3)/rad;
44 :     }
45 :     force -= normalize(∇F(pos))⊗normalize(∇F(pos))•force;
46 :     if (|force| > 0) {
47 :     // project force onto tangent plane
48 :     // take gradient descent step
49 :     vec2 down = normalize(force);
50 :     vec2 posLast = pos;
51 :     pos += hh*down;
52 :     // take Newton-Raphson step back to surface
53 :     pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos));
54 :     pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos));
55 :     real newenergy = 0;
56 :     foreach (Particle p_j2 in sphere(rad)) {
57 :     vec2 r_ij = (pos - p_j2.pos)/rad;
58 :     newenergy += (1 - |r_ij|)^4;
59 : lamonts 2950 }
60 : lamonts 3013 if (newenergy > energy - 0.5*hh*|force|) {
61 :     hh *= 0.5; // backtrack
62 :     pos = posLast;
63 :     } else {
64 :     // re-initialize for new particle neighborhood
65 :     hh = hhInit;
66 :     }
67 :     }
68 :     } // if (!foundContour)
69 :     }
70 : lamonts 2950 }
71 :     global{
72 : lamonts 3013 real energyMean = mean{P.energy | P in Particle.all};
73 :     evariance = mean{ (P.energy - energyMean) * (P.energy - energyMean) | P in Particle.all};
74 : lamonts 2950 }
75 :    
76 : lamonts 3013 initially { Particle(ui + gridSize*vi,
77 :     [lerp(-1.5, 1.5, -0.5, real(ui), real(gridSize)-0.5),
78 :     lerp(-1.5, 1.5, -0.5, real(vi), real(gridSize)-0.5)])
79 :     | vi in 0..(gridSize-1), ui in 0..(gridSize-1) };

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