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

SCM Repository

[diderot] Annotation of /branches/tash/test/iso2d-cbrt.diderot
ViewVC logotype

Annotation of /branches/tash/test/iso2d-cbrt.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1685 - (view) (download)
Original Path: branches/vis12/test/iso2d-cbrt.diderot

1 : jhr 1115 // iso2d-cbrt
2 :     //
3 :     // Demo of difficulty of finding isocontours via Newton-Raphson method
4 :     // with nasty input; in this case it follows the shape of cbrt(r),
5 :     // so many of the points dramatically fail to converge.
6 :     //
7 :     // Process output with:
8 :     // unu jhisto -i mip.txt -b 500 500 -min 0 0 -max 499 499 | unu 2op neq - 0 | unu quantize -b 8 -o iso2d-cbrt.png
9 :    
10 :    
11 :     int gridSize = 150;
12 :     real isoval = 0.0;
13 :     field#2(2)[] F = load("../data/cbrt-500-unit.nrrd") ⊛ bspln3;
14 :     int stepsMax = 40;
15 :     real epsilon = 0.001;
16 :    
17 :     strand sample (int ui, int vi) {
18 :     output vec2 pos = [lerp(1.0, 499.0, -0.5, real(ui), real(gridSize)-0.5),
19 :     lerp(1.0, 499.0, -0.5, real(vi), real(gridSize)-0.5)];
20 :     int steps = 0;
21 :     update {
22 :     // We bail if we're no longer inside or taken too many steps.
23 :     // We tag this strand result as garbage by setting it to a
24 :     // non-finite position; these will be filtered later.
25 :     // This is really what "die" is for, though (see below).
26 :     if (!inside(pos, F) || steps >= stepsMax) {
27 :     pos = [nan,nan];
28 :     // HEY (BUG) we shouldn't even have to use non-finite values;
29 :     // we should be able to just say "die;". Currently, using
30 :     // "die" leads to the program never finishing...
31 :     stabilize;
32 :     }
33 :     vec2 grad = ∇F(pos);
34 :     real gmsq = grad • grad;
35 :     // If |∇F|=0 we can't compute update; we have to bail
36 :     if (gmsq == 0.0) {
37 :     pos = [nan,nan];
38 :     stabilize;
39 :     }
40 :     // delta = Newton-Raphson step
41 :     vec2 delta = -(F(pos) - isoval)*grad/gmsq;
42 :     // We've converged if the update is small enough
43 :     if (|delta| < epsilon) {
44 :     stabilize;
45 :     }
46 :     pos = pos + delta;
47 :     steps = steps + 1;
48 :     }
49 :     }
50 :    
51 :     initially [ sample(ui, vi) | vi in 0..(gridSize-1), ui in 0..(gridSize-1) ];

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