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

SCM Repository

[diderot] Annotation of /tests/vis15-bugs/halftone-bug2.diderot
ViewVC logotype

Annotation of /tests/vis15-bugs/halftone-bug2.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4929 - (view) (download)

1 : glk 4710 /* ==========================================
2 :     Half-toning via interacting repulsive particles
3 :    
4 :     This example is based on the [`sphere.diderot`](../sphere) example,
5 :     which in turn is based on the [`circle.diderot`](../circle) example;
6 :     both of these should be read first. While in those cases a surface
7 :     (circle or sphere) is sampled uniformly, in this case the flat surface
8 :     of an image is sampled non-uniformly, according to the image
9 :     intensity.
10 :    
11 :     ... TODO: example of running program and looking at output
12 :    
13 :     ========================================== */
14 :    
15 :     input image(2)[] img ("data in which to find features") = image("halftone-bug2-img.nrrd");
16 :     field#0(2)[] F = ctmr ⊛ clamp(img);
17 :    
18 :     input vec2{} ipos ("initial positions for all particles") = load("halftone-bug2-vec2.nrrd");
19 :     input vec2 radmm ("particle minimum > 0 and maximum radius") = [0.01, 0.1];
20 :     input real eps ("system convergence threshold, computed as the coefficient-of-variation of normalized distances to nearest neighbors") = 0.05;
21 :     input int pcp ("periodicity of population control (if greater than zero)") = 2;
22 :     input real hhInit ("initial step size for potential energy gradient descent") = 1;
23 :    
24 :     real newDist = ((0.666)); // how far away to birth, as multiple of radius
25 :     real stepMax = ((0.666)); // limit on travel per iter, as multiple of radius
26 :     int iter = 0; // which iteration we're on
27 :    
28 :     // potential function phi(r) found via Mathematica: this is the lowest-degree
29 :     // piecewise polynomial with phi(0)=1 and a potential well at phi(0.666)=-0.005.
30 :     // phi(r) is C^3 continuous across the well and with 0 at r=1.
31 :     function real phi(real r) =
32 :     1 + r*(-4.248582222661734 + r*(5.991287388535527 + (-2.864766048887071 + 0.06790595504541913*r)*r))
33 :     if r < 0.666 else
34 :     -0.005 + ((0.4482053856358993 + (-2.6838645846460842 + (6.026642031390917 - 4.811690244623482*(-0.666 + r))*(-0.666 + r))*(-0.666 + r))*(-0.666 + r))*(-0.666 + r)
35 :     if r < 1 else 0;
36 :     function real phi'(real r) =
37 :     -4.248582222661734 + r*(11.982574777071054 + (-8.594298146661213 + 0.2716238201816765*r)*r)
38 :     if r < 0.666 else
39 :     -15.42591894092919 + 0.4482053856358993*(-1.332 + 2*r) + r*(-2.6838645846460842*(-3.996 + 3*r) + 6.026642031390917*(5.322672 + r*(-7.992 + 4*r)) - 4.811690244623482*(-5.90816592 + r*(13.30668 + r*(-13.32 + 5*r))))
40 :     if r < 1 else 0;
41 :    
42 :     // Energy and force from particle (with radius rad) at vec2 x
43 :     function real enr(vec2 x, real rad) = rad*rad*phi(|x|/rad);
44 :     function vec2 frc(vec2 x, real rad) = rad*rad*phi'(|x|/rad) * (1/rad) * x/|x|;
45 :    
46 :     // From a given vec2, find a random-ish value uniformly sampling [0,1)
47 :     function real posrnd(vec2 v, real rad) {
48 :     vec2 p = 10000*v/rad;
49 :     return fmod(|fmod(p[0],1)| + |fmod(p[1],1)|, 1);
50 :     }
51 :    
52 :     // Is this an iteration in which to do population control?
53 :     function bool pcIter() = (pcp > 0 && iter > 0 && 0 == iter % pcp);
54 :     function bool birth() = ( (iter/pcp) % 2 == 0);
55 :    
56 :     function real radius(real ff) =
57 :     min(radmm[1], radmm[0]*sqrt(1/clamp(0.0000001, 1, ff)));
58 :    
59 :    
60 :     /* The particle is initialized at position pos0, with initial stepsize hh0.
61 :     The first set of particles gets hhInit for the initial stepsize, but new
62 :     particles created by population control benefit from getting the stepsize
63 :     that was adaptively learned by the parent. */
64 :     strand particle (int gen_, vec2 pos0, real hh0) {
65 :     int gen = gen_;
66 :     vec2 pos = pos0;
67 :     output vec3 posn = [pos[0], pos[1], 3];
68 :     real hh = hh0;
69 :     vec2 step = [0,0]; // step along force
70 :     real rad = 0; // zero isn't a valid value
71 :     real closest = 1; // distance to closest neighbor (as mult of radius)
72 :     int ncount = 0; // how many neighbors did we have
73 :     real mvmt = 1;
74 :     real energy=0;
75 :     update {
76 :     if (!inside(pos, F)) { die; }
77 :     rad = radius(F(pos));
78 :     // compute energy and forces on us from neighbors
79 :     real energyLast = 0;
80 :     vec2 force = [0,0];
81 :     ncount = 0;
82 :     foreach (particle P in sphere(2*rad)) {
83 :     vec2 x = P.pos - pos;
84 :     real rr = radius(F(lerp(pos, P.pos, 0.5)));
85 :     energyLast += enr(x, rr);
86 :     force += frc(x, rr);
87 :     ncount += 1 if |x| < rr else 0;
88 :     }
89 :    
90 :     /* Else have interacting neighbors; find step, limit step size */
91 :     step = hh*force;
92 :     if (|step| > stepMax*rad) {
93 :     hh *= stepMax*rad/|step|;
94 :     step = hh*force;
95 :     }
96 :     mvmt = lerp(mvmt, |step|/(rad*((0.666))), 0.5);
97 :    
98 :     vec2 posLast = pos;
99 :     pos += step;
100 :     energy = 0;
101 :     closest = 1;
102 :     ncount = 0;
103 :     vec2 mno = [0,0]; // mean neighbor offset
104 :     if (|pos| < 0.025 && (iter==5 || iter==6)) {
105 : jhr 4929 print("\n(iter ", iter, ") HELLO (gen ", gen, ") at pos ", pos, " (|pos|=",
106 :     |pos|, "), looking for neighbors with base radius ", rad, "\n");
107 : glk 4710 }
108 :     foreach (particle P in sphere(2*rad)) {
109 :     real rr = radius(F(lerp(pos, P.pos, 0.5)));
110 :     vec2 x = P.pos - pos;
111 : jhr 4929 // if (|pos| < 0.025 && (iter==5 || iter==6)) {
112 :     // print(" neighbor at ", P.pos, " ==> halfway ", lerp(pos, P.pos, 0.5),
113 :     // " F=", F(lerp(pos, P.pos, 0.5)), " ==> rr=", rr, "; |x|=", |x|,
114 :     // "<" if |x| < rr else ">=", rr, "\n");
115 :     // }
116 : glk 4710 if (|x| < rr) {
117 :     energy += enr(x, rr);
118 :     closest = min(closest, |x|/rr);
119 :     mno += x;
120 :     ncount += 1;
121 :     if (|pos| < 0.025 && (iter==5 || iter==6)) {
122 : jhr 4929 print("\n REAL neighbor P at ", P.pos, " (|P.pos|==", |P.pos|,
123 :     "==" if |P.pos|==0 else "!=", "0; P.gen ", P.gen, "); now, energy=",
124 :     energy, ", closest=", closest, "; ncount=", ncount, "\n");
125 : glk 4710 }
126 :     }
127 :     }
128 :     mno /= ncount;
129 :    
130 :     // Armijo-Goldstein sufficient decrease condition
131 :     if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
132 :     hh *= 0.5; // energy didn't decrease as expected: backtrack
133 :     pos = posLast; // try again next iteration
134 :     } else {
135 :     hh *= 1.05; // opportunistically increase step size
136 :     // try to have between 6 and 8 neighbors
137 :     if (pcIter()) {
138 :     if (birth()) {
139 :     if (energy < 0 && posrnd(pos, rad) < (6.0 - ncount)/6.0) {
140 :     real a = atan2(-mno[1],-mno[0]);
141 :     a += lerp(-1, 1, gen % 2)*π/6.0;
142 :     vec2 npos = pos + newDist*rad*[cos(a),sin(a)];
143 :     print("\n(iter ", iter, ") NEW particle(", gen+1, ",", npos, ",", hh, ")\n");
144 :     new particle(gen+1, npos, hh);
145 :     }
146 :     } else if (energy > 0 && posrnd(pos, rad) < (ncount - 8.0)/ncount) {
147 :     if (|pos| < 0.1 && (iter==5 || iter==6)) {
148 : jhr 4929 print("(iter ", iter, ") DIE at pos ", pos, " (|pos|=", |pos|, "): energy ",
149 :     energy, "; ncount ", ncount, "; rad ", rad, "\n");
150 : glk 4710 print(" ", posrnd(pos, rad), " < ", (ncount - 8.0)/ncount, "\n");
151 :     }
152 :     die;
153 :     }
154 :     }
155 :     }
156 :    
157 :     posn = [pos[0], pos[1], ncount];
158 :     if (|pos| == 0) {
159 :     print("(iter ", iter, ") WHOA pos exactly at 0!\n");
160 :     }
161 :    
162 :     // Record actual step taken, in case used in global update
163 :     step = pos - posLast;
164 :     }
165 :     }
166 :    
167 :     global {
168 :     /* Compute coefficient-of-variation of distance to closest neighbor */
169 :     real totenr = sum { P.energy | P in particle.all};
170 :     real meancl = mean { P.closest | P in particle.all};
171 :     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
172 :     real covcl = sqrt(varicl) / meancl;
173 :     print("\n=-=-=\n=-=-= iter ", iter,
174 :     "; tot#=", numActive(),
175 :     "; totenr=", totenr,
176 :     "; COV(cl)=", covcl,
177 :     "; mean(cl)=", meancl,
178 :     "; mean(hh)=", mean { P.hh | P in particle.all},
179 :     "; mean(ncount)=", mean { real(P.ncount) | P in particle.all},
180 :     "; mean(mvmt)=", mean { P.mvmt | P in particle.all },
181 :     "\n");
182 :     iter += 1;
183 :     }
184 :    
185 :     initially { particle(2, ipos[ii], hhInit) | ii in 0 .. length(ipos)-1 };

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