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

SCM Repository

[diderot] Annotation of /branches/tash/test/particles-iso.diderot
ViewVC logotype

Annotation of /branches/tash/test/particles-iso.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1941 - (view) (download)
Original Path: branches/vis12/test/particles-iso.diderot

1 : glk 479 // particles-iso.diderot
2 :     //
3 :     // distributes particles on an isosurface
4 :     //
5 :    
6 :     input real radius = 0.2; // radius of particle
7 :     input real isovalue = 0.0; // isocontour value
8 :     input real cnstrEps = 0.01; // epsilon test for constraint satisfaction
9 :     input real travelMax = 5.0; // point can't move further than this
10 :     input int cnstrIterMax = 10; // constraint satisfaction can't iterate more than this
11 :     input real stepLimit = 2; // limit on how far a particle can move in one iteration
12 :     input real stepBack = 0.2; // how to reduce step size if badstep
13 :     input real stepUp = 1.1; // opportunistic scaling of step
14 :     input int pointNum = 100; // number of initial seed points
15 :     input real energyConvg = 0.1; // convergence threshold on system energy
16 :    
17 : jhr 1941 input image(3)[] img;
18 : glk 479 field#2(3)[] F = img ⊛ bspln3;
19 :    
20 :     // should allow functions to take actors, but sensible semantics
21 :     // limit that to a read-only view of actor at start of iteration;
22 :     // compiler can warn if the actor is passed post-killing-definition
23 :    
24 :     // uses Newton-Raphson to find the zero-crossing
25 :     (bool, vec3) constrain (vec3 posOrig)
26 :     {
27 :     vec3 grad, fix, pos = posOrig;
28 :     int iter = 0;
29 :     do {
30 :     // implicitly a (iter,update,pos) tuple is iterated by this
31 :     iter += 1;
32 :     grad = ∇F@pos;
33 :     fix = F@pos*grad / |grad|^2 // should optimize out sqrt()^2
34 :     pos -= fix;
35 :     if (dot(pos - posOrig,pos - posOrig) > travelMax*travelMax) {
36 :     // ?can use "NaN" as a type real literal for IEEE Not-A-Number?
37 :     return (false, [NaN,NaN,NaN])
38 :     }
39 :     if (iter > constrIterMax) {
40 :     return (false, [NaN,NaN,NaN])
41 :     }
42 :     } while (dot(fix,fix) > cnstrEps*cnstrEps)
43 :     // we've converged
44 :     return (true,pos);
45 :     }
46 :    
47 :     // this is the sort of function that would benefit from automatic
48 :     // differentiation; it simply evaluates a piece-wise polynomial
49 :     // along with its derivative
50 :     (real, real) phi (real x)
51 :     {
52 :     // this directly copied from C code (teem/src/pull/energy.c);
53 :     real wx = 0.66, wy = 0.005;
54 :     real enr, denr; // output
55 :    
56 :     real xmo = wx-1;
57 :     real xmoo = xmo*xmo;
58 :     real xmooo = xmoo*xmo;
59 :     if (x < wx) {
60 :     real a = -3*(xmoo + (-1 + 2*wx)*wy)/(xmoo*wx);
61 :     real b = 3*(xmoo + (-1 + wx*(2 + wx))*wy)/(xmoo*wx*wx);
62 :     real c = (-1 + wy - wx*(-2 + wx + 2*(1 + wx)*wy))/(xmoo*wx*wx*wx);
63 :     enr = 1 + x*(a + x*(b + x*c));
64 :     denr = a + x*(2*b + x*3*c);
65 :     } else if (x < 1) {
66 :     real d = ((-1 + 3*wx)*wy)/xmooo;
67 :     real e = -(6*wx*wy)/xmooo;
68 :     real f = (3*(1 + wx)*wy)/xmooo;
69 :     real g = -(2*wy)/xmooo;
70 :     enr = d + x*(e + x*(f + x*g));
71 :     denr = e + x*(2*f + x*3*g);
72 :     } else {
73 :     enr = 0;
74 :     denr = 0;
75 :     }
76 :     return (enr, denr);
77 :     }
78 :    
79 :     // evaluates the inter-particle energy function in 3D
80 :     (real, vec3) pairEnergy(vec3 here, vec3 there) {
81 :     vec3 dir = (there - here)/|there - here|;
82 :     (real enr, real denr) = phi(|there - here|/radius);
83 :     return (enr, -denr*dir);
84 :     }
85 :    
86 :     // computes energy and force from nearby particles
87 :     (real, vec3) totalEnergy(vec3 here) {
88 :     // ?what is our sequence/list/iterable type?
89 :     // ?can use ".pos" to refer to position variable?
90 :     // ?neighbors has to know not return actor at "here" exactly?
91 :     seq(vec3) neighPos = [n.pos for n in neighbors(radius, here)];
92 :     real energy = 0;
93 :     vec3 force = [0,0,0];
94 :     for p in neighPos {
95 :     (energy, force) += pairEnergy(pos, n.pos);
96 :     }
97 :     return (energy, force);
98 :     }
99 :    
100 :     actor IsoPoint (vec3 posInit)
101 :     {
102 :     real step = 1.0; // initial step size for gradient descent
103 :     output real energy;
104 :     output vec3 pos;
105 :    
106 :     initially {
107 :     (bool ok, vec3 pos) = constrain(posInit);
108 :     if (!ok) {
109 :     // we couldn't satisfy contraint from seed location; we're gone
110 :     delete; // ?what's the way for an actor to remove itself?
111 :     }
112 :     }
113 :    
114 :     update
115 :     {
116 :     bool badstep;
117 :     vec3 posNew, force;
118 :     real energNew, energyOrig;
119 :    
120 :     (energyOrig, force) = totalEnergy(pos);
121 : glk 480 // project force onto tangent plane
122 : glk 481 vec3 norm = ∇F@pos/|∇F@pos|;
123 :     force -= outer(norm, norm)*force
124 : glk 479 // ensures that step*force is a reasonable distance
125 :     step = min(step, travelMax/|force|);
126 :     do {
127 :     (bool ok, posNew) = constrain(pos + step*force);
128 :     if (!ok) {
129 :     delete;
130 :     }
131 :     // there should be non-trivial savings from not computing force
132 :     (energyNew,) = pointEnergy(pos, posNew);
133 :     // use of badstep here is rather clunky
134 :     badstep = energyNew > energyOrig;
135 :     if (badstep) {
136 :     step *= stepBack;
137 :     }
138 :     } while badstep;
139 :     step *= stepUp;
140 :     pos = posNew;
141 :     energy = energyNew;
142 :     }
143 :    
144 :     // (stabilization is not determined per-actor)
145 :    
146 :     }
147 :    
148 :     // ?from here on out I'm making things up ...?
149 :    
150 :     // need some way of generating a list of random locations inside image,
151 :     // and then seeding actors there
152 :     initially [ IsoPoint(pos) for pos in randomPositionsInside(F, pointNum) ];
153 :    
154 :     // keep update()ing all actors until total change in system is
155 :     // below some threshold
156 :     real Eold, Enew = reduce(sum, [a.energy for a in actors]);
157 :     do {
158 :     Eold = Enew;
159 :     updateAll;
160 :     Enew = reduce(sum, [a.energy ]);
161 :     } while (Eold - Enew > energyConvg);

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