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

Parent Directory | 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