Annotation of /trunk/test/particles-iso.diderot

Parent Directory | Revision Log

---

Revision 479 - (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 :			// ensures that step*force is a reasonable distance
123 :			step = min(step, travelMax/|force|);
124 :			if (step*|force| > travelMax) {
125 :			step = travelMax/|force|;
126 :			}
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