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

SCM Repository

[diderot] Annotation of /branches/pure-cfg/test/iso2d.diderot
ViewVC logotype

Annotation of /branches/pure-cfg/test/iso2d.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 995 - (view) (download)

1 : glk 953 // iso2d
2 :     //
3 :     // Demo of finding isocontours via Newton-Raphson method.
4 :     // Initializes positions on a grid, and each update applies one
5 :     // step of Newton-Raphson.
6 :     //
7 :     // Process output with:
8 : glk 995 // unu jhisto -i mip.txt -b 512 512 -min 0 0 -max 1 1 | unu 2op neq - 0 | unu quantize -b 8 -o iso2d.png
9 : glk 953
10 : glk 995 int gridSize = 300;
11 :     image(2)[] img = load("../data/ddro-80.nrrd");
12 : glk 953 field#1(2)[] F = img ⊛ ctmr;
13 : glk 954 int stepsMax = 10;
14 : glk 995 real epsilon = 0.0001;
15 : glk 953
16 :     strand sample (int ui, int vi) {
17 : glk 995 output vec2 pos = [lerp(0.0, 1.0, -0.5, real(ui), real(gridSize)-0.5),
18 :     lerp(0.0, 1.0, -0.5, real(vi), real(gridSize)-0.5)];
19 :     // set the isvalue to 50, 30, or 10, depending on whichever we're closest to
20 :     real isoval = 50.0 if F(pos) >= 40.0 else 30.0 if F(pos) >= 20.0 else 10.0;
21 : glk 953 int steps = 0;
22 :     update {
23 :     // We bail if we're no longer inside or taken too many steps.
24 :     // We tag this strand result as garbage by setting it to a
25 :     // non-finite position; these will be filtered later.
26 :     // This is really what "die" is for, though (see below).
27 : glk 954 if (!inside(pos, F) || steps > stepsMax) {
28 : jhr 966 pos = [nan,nan];
29 :     // HEY (BUG) we shouldn't even have to use non-finite values;
30 :     // we should be able to just say "die;". Currently, using
31 :     // "die" leads to the program never finishing...
32 :     stabilize;
33 : glk 953 }
34 :     vec2 grad = ∇F(pos);
35 :     real gmsq = grad • grad;
36 :     // If |∇F|=0 we can't compute update; we have to bail
37 : glk 976 if (gmsq == 0.0) {
38 : jhr 966 pos = [nan,nan];
39 :     stabilize;
40 : glk 953 }
41 :     // delta = Newton-Raphson step
42 :     vec2 delta = -(F(pos) - isoval)*grad/gmsq;
43 :     // We've converged if the update is small enough
44 :     if (|delta| < epsilon) {
45 : jhr 966 stabilize;
46 : glk 953 }
47 :     pos = pos + delta;
48 :     steps = steps + 1;
49 :     }
50 :     }
51 :    
52 :     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