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

SCM Repository

[diderot] Annotation of /branches/vis12/test/lzc2d-brute.diderot
ViewVC logotype

Annotation of /branches/vis12/test/lzc2d-brute.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1955 - (view) (download)

1 : glk 1955 // lzc2d-brute
2 :     //
3 :     // Demo of laplacian zero-crossing edges via brute-force search
4 :     // and then Illinois false position
5 :     //
6 :     // Process output with:
7 :     // unu jhisto -i pos.nrrd -b 512 512 -min 0 0 -max 1 1 | unu 2op neq - 0 | unu quantize -b 8 -o lzc2d.png
8 :    
9 :     int gridSize = 300;
10 :     field#3(2)[] F = bspln5 ⊛ image("../data/ddro-80.nrrd");
11 :     int stepsMax = 30;
12 :     real epsilon = 0.00001;
13 :    
14 :     strand sample (int ui, int vi) {
15 :     output vec2 pos = [lerp(0.0, 1.0, -0.5, real(ui), real(gridSize)-0.5),
16 :     lerp(0.0, 1.0, -0.5, real(vi), real(gridSize)-0.5)];
17 :     vec2 posA=pos;
18 :     vec2 posB=pos;
19 :     int stepsNum = 0;
20 :     int side = 0;
21 :     real step = 0.004;
22 :     real fA=0.0;
23 :     real fB=0.0;
24 :     real A=0.0;
25 :     real B=1.0;
26 :     bool searching = true;
27 :     update {
28 :     // We bail if we're no longer inside or taken too many steps.
29 :     if (!inside(pos, F) || stepsNum > stepsMax) {
30 :     die;
31 :     }
32 :     if (searching) {
33 :     posA = pos;
34 :     real ll = trace(∇⊗∇F(pos));
35 :     vec2 dir = normalize(∇F(pos));
36 :     real sgn = 1.0 if ll > 0.0 else -1.0;
37 :     posB = pos + sgn*step*dir;
38 :     searching = trace(∇⊗∇F(posB))*ll > 0.0;
39 :     if (!searching) {
40 :     /* set up next phase */
41 :     fA = trace(∇⊗∇F(posA));
42 :     fB = trace(∇⊗∇F(posB));
43 :     }
44 :     } else {
45 :     pos = lerp(posA, posB, fA, 0.0, fB);
46 :     real fs = trace(∇⊗∇F(pos));
47 :     if (0.0 == fs) {
48 :     stabilize;
49 :     }
50 :     if (fs*fB > 0.0) { /* not between s and b */
51 :     posB = pos;
52 :     fB = fs;
53 :     if (1 == side) {
54 :     fA /= 2.0;
55 :     }
56 :     side = 1;
57 :     } else { /* not between a and s */
58 :     posA = pos;
59 :     fA = fs;
60 :     if (-1 == side) {
61 :     fB /= 2.0;
62 :     }
63 :     side = -1;
64 :     }
65 :     if (|posA - posB| < epsilon ) {
66 :     if (|∇F(pos)| < 300.0) {
67 :     die;
68 :     }
69 :     stabilize;
70 :     }
71 :     }
72 :     stepsNum += 1;
73 :     }
74 :     }
75 :    
76 :     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