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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1685 - (view) (download)

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

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