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

SCM Repository

[diderot] Annotation of /branches/vis12/test/iso2d.diderot
ViewVC logotype

Annotation of /branches/vis12/test/iso2d.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2174 - (view) (download)

1 : jhr 1115 // 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 : jhr 2173 // unu jhisto -i pos.nrrd -b 512 512 -min 0 0 -max 1 1 | unu 2op neq - 0 | unu quantize -b 8 -o iso2d.png
9 : jhr 1115
10 :     int gridSize = 300;
11 : jhr 1927 field#1(2)[] F = ctmr ⊛ image("../data/ddro-80.nrrd");
12 : jhr 1115 int stepsMax = 10;
13 :     real epsilon = 0.0001;
14 :    
15 :     strand sample (int ui, int vi) {
16 :     output vec2 pos = [lerp(0.0, 1.0, -0.5, real(ui), real(gridSize)-0.5),
17 :     lerp(0.0, 1.0, -0.5, real(vi), real(gridSize)-0.5)];
18 :     // set the isvalue to 50, 30, or 10, depending on whichever we're closest to
19 : glk 1134 real isoval = 50.0 if F(pos) >= 40.0
20 :     else 30.0 if F(pos) >= 20.0
21 :     else 10.0;
22 : jhr 2174 // the following won't work, because the compiler doesn't know how to deal with state vars that
23 :     // are fields (really just a local binding)!
24 :     // field#1(2)[] G = F - isoval;
25 : jhr 1115 int steps = 0;
26 :     update {
27 : jhr 2173 field#1(2)[] G = F - isoval;
28 : jhr 1115 // We bail if we're no longer inside or taken too many steps.
29 : jhr 2173 if (!inside(pos, G) || steps > stepsMax) {
30 : jhr 1131 die;
31 : jhr 1115 }
32 : glk 1329 // GLKs recent changes (revision 1329) were made to make the code
33 :     // a more obvious implementation of Newton-Raphson, and also to create
34 :     // some obvious and not-so-obvious opportunties for optimization
35 : jhr 2173 vec2 grad = ∇G(pos);
36 : glk 1329 if (|grad| == 0.0) { // can't compute step if |∇F|, so have to bail
37 : jhr 1131 die;
38 : jhr 1115 }
39 : glk 1329 vec2 norm = normalize(grad);
40 : jhr 2173 vec2 delta = -(G(pos)/|grad|)*norm; // Newton-Raphson step
41 : glk 1329 if (|delta| < epsilon) { // we've converged if step is small enough
42 : jhr 1115 stabilize;
43 :     }
44 : glk 1329 pos += delta;
45 :     steps += 1;
46 : jhr 1115 }
47 :     }
48 :    
49 :     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