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

SCM Repository

[diderot] Annotation of /benchmarks/programs/lic2d/bmark-diderot.diderot
ViewVC logotype

Annotation of /benchmarks/programs/lic2d/bmark-diderot.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1592 - (view) (download)

1 : nseltzer 1592 // lic-turb.diderot
2 :     //
3 :     // demo of line integral convolution (LIC) on 2D turbulent flow
4 :     //
5 :     // First, need to create padded vector dataset v.nrrd with:
6 :     // unu pad -i ../data/turb2d.nrrd -min 0 -5 -5 -max M M+5 M+5 -b pad -v 0.0001 -o v.nrrd
7 :     // We pad with value 0.0001; we padded with zero than normalize() would fail
8 :     // But if you try take too many steps (increasing stepNum) then eventually
9 :     // you will leave the vector field domain, which will crash the program,
10 :     // so really the padding should be done with a vector field that
11 :     // pointwards inwards, towards the center of the domain, so that anything
12 :     // leaving the domain will get pushed back in. OR, really, the code
13 :     // below should be smart and use inside() or some other test to stop
14 :     // the integration when we get too close to the domain boundary.
15 :     //
16 :     // Then, need to create noise texture with:
17 :     // unu slice -i v.nrrd -a 0 -p 0 | unu resample -s 1020 561 | unu 1op nrand -o R.nrrd
18 :     // where "1020 561" is copied from imgSizeX and imgSizeY below; we start
19 :     // with unu slice in order to inherit the image orientation info, to
20 :     // create a noise texture at the same world-space position as the data
21 :     //
22 :     // process output with:
23 :     // unu reshape -i lic-turb2d.txt -s 3 1020 561 | unu quantize -b 8 -o lic-turb2d.ppm
24 :    
25 :    
26 :     int imgSizeX = 1020;
27 :     int imgSizeY = 561;
28 :     real h = 0.005; // step size of integration
29 :     int stepNum = 10; // take this many steps both upstream and downstream
30 :    
31 :     field#2(2)[2] boundary = load("../../data/turb2d.nrrd") ⊛ bspln3;
32 :     field#2(2)[2] V = load("../../data/v.nrrd") ⊛ bspln3;
33 :     field#0(2)[] R = load("../../data/R.nrrd") ⊛ tent;
34 :    
35 :     strand LIC (int xi, int yi) {
36 :     real xx = lerp(0.04, 6.76, 0.0, real(xi), real(imgSizeX)-1.0);
37 :     real yy = lerp(0.04, 3.7, 0.0, real(yi), real(imgSizeY)-1.0);
38 :     vec2 pos0 = [xx,yy];
39 :     vec2 forw = pos0;
40 :     vec2 back = pos0;
41 :     real sum = R(pos0);
42 :     output vec3 RGB = [0.0, 0.0, 0.0];
43 :     real l = 0.0;
44 :     int step = 0;
45 :     int num = 1;
46 :     real crl = 0.0;
47 :    
48 :     update {
49 :     // Euler integration step
50 :     // forw = forw + h*V(forw);
51 :     // back = back - h*V(back);
52 :     // Midpoint method step
53 :     if(inside(forw, boundary)) {
54 :     forw += h*normalize(V(forw + 0.5*h*normalize(V(forw))));
55 :     sum += R(forw);
56 :     num += 1;
57 :     }
58 :     if(inside(back, boundary)) {
59 :     back -= h*normalize(V(back - 0.5*h*normalize(V(back))));
60 :     sum += R(back);
61 :     num += 1;
62 :     }
63 :     step += 1;
64 :     if((!inside(forw, boundary) && !inside(back, boundary)) || step == stepNum) {
65 :     // the pow() is a way to adjust how much velocity modulates
66 :     // contrast; lower values increase contrast at small velocity
67 :     sum = pow(|V(pos0)|,0.3)*sum/real(num);
68 :    
69 :     tensor[2,2] M = ∇⊗V(pos0);
70 :     crl = M[0,1] - M[1,0];
71 :     // BUG - curl ("∇×") doesn't work
72 :     //crl = ∇×V(pos0);
73 :    
74 :     //sum in range -1.56, 2.58
75 :     //curl in rage -131.797974, 131.828476
76 :    
77 :     if(crl < 0.0) {
78 :     real h = 300.0;
79 :     }
80 :     else {
81 :     real h = 120.0;
82 :     }
83 :    
84 :     real s = lerp(0.0, 1.0, 0.0, |crl|, 132.0);
85 :    
86 :     l = lerp(0.0, 1.0, -1.6, sum, 2.6);
87 :    
88 :     real c = (1.0 - |2.0 * l - 1.0|) * s;
89 :    
90 :     if(crl < 0.0) {
91 :     int hp = 5;
92 :     real x = c;
93 :     real m = l - c / 2.0;
94 :     RGB = [c + m, 0.0 + m, x + m];
95 :     }
96 :     else {
97 :     int hp = 2;
98 :     real x = 0.0;
99 :     real m = l - c / 2.0;
100 :     RGB = [0.0 + m, c + m, x + m];
101 :     }
102 :    
103 :     stabilize;
104 :     }
105 :     }
106 :     }
107 :    
108 :     initially [ LIC(xi, yi) | yi in 0..(imgSizeY-1), xi in 0..(imgSizeX-1) ];

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