// lic-turb.diderot
// 
// demo of line integral convolution (LIC) on 2D turbulent flow
//
// First, need to create padded vector dataset v.nrrd with:
// unu pad -i ../data/turb2d.nrrd -min 0 -5 -5 -max M M+5 M+5 -b pad -v 0.0001 -o v.nrrd
// We pad with value 0.0001; we padded with zero than normalize() would fail
// But if you try take too many steps (increasing stepNum) then eventually
// you will leave the vector field domain, which will crash the program,
// so really the padding should be done with a vector field that 
// pointwards inwards, towards the center of the domain, so that anything
// leaving the domain will get pushed back in.  OR, really, the code
// below should be smart and use inside() or some other test to stop
// the integration when we get too close to the domain boundary.
//
// Then, need to create noise texture with:
// unu slice -i v.nrrd -a 0 -p 0 | unu resample -s 1020 561 | unu 1op nrand -o R.nrrd
// where "1020 561" is copied from imgSizeX and imgSizeY below; we start
// with unu slice in order to inherit the image orientation info, to
// create a noise texture at the same world-space position as the data
//
// process output with:
// unu reshape -i lic-turb2d.txt -s 3 1020 561 | unu quantize -b 8  -o lic-turb2d.ppm


int imgSizeX = 1020;
int imgSizeY = 561;
real h = 0.005;    // step size of integration
int stepNum = 10; // take this many steps both upstream and downstream

field#2(2)[2] boundary = load("../../data/turb2d.nrrd") ⊛ bspln3;
field#2(2)[2] V = load("../../data/v.nrrd") ⊛ bspln3;
field#0(2)[] R = load("../../data/R.nrrd") ⊛ tent;

strand LIC (int xi, int yi) {
	real xx = lerp(0.04, 6.76, 0.0, real(xi), real(imgSizeX)-1.0);
	real yy = lerp(0.04, 3.7, 0.0, real(yi), real(imgSizeY)-1.0);
	vec2 pos0 = [xx,yy];
	vec2 forw = pos0;
	vec2 back = pos0;
	real sum = R(pos0);
	output vec3 RGB = [0.0, 0.0, 0.0];
	real l = 0.0;
	int step = 0;
	int num = 1;
	real crl = 0.0;

	update {
		// Euler integration step
		// forw = forw + h*V(forw);
		// back = back - h*V(back);
		// Midpoint method step
		if(inside(forw, boundary)) {
			forw += h*normalize(V(forw + 0.5*h*normalize(V(forw))));
			sum += R(forw);
			num += 1;
		}
		if(inside(back, boundary)) {
			back -= h*normalize(V(back - 0.5*h*normalize(V(back))));
			sum += R(back);
			num += 1;
		}
		step += 1;
		if((!inside(forw, boundary) && !inside(back, boundary)) || step == stepNum) {
			// the pow() is a way to adjust how much velocity modulates
			// contrast; lower values increase contrast at small velocity
			sum = pow(|V(pos0)|,0.3)*sum/real(num);

			tensor[2,2] M = ∇⊗V(pos0);
			crl = M[0,1] - M[1,0];
			// BUG - curl ("∇×") doesn't work
			//crl = ∇×V(pos0);

			//sum in range -1.56, 2.58
			//curl in rage -131.797974, 131.828476

			if(crl < 0.0) {
				real h = 300.0;
			}
			else {
				real h = 120.0;
			}

			real s = lerp(0.0, 1.0, 0.0, |crl|, 132.0);

			l = lerp(0.0, 1.0, -1.6, sum, 2.6);

			real c = (1.0 - |2.0 * l - 1.0|) * s;

			if(crl < 0.0) {
				int hp = 5;
				real x = c;
				real m = l - c / 2.0;
				RGB = [c + m, 0.0 + m, x + m];
			}
			else {
				int hp = 2;
				real x = 0.0;
				real m = l - c / 2.0;
				RGB = [0.0 + m, c + m, x + m];
			}

			stabilize;
		}
	}
}

initially [ LIC(xi, yi) | yi in 0..(imgSizeY-1), xi in 0..(imgSizeX-1) ];

<script type="text/javascript">//<![CDATA[
		function toggle_ffErrors() {
			var errorsblock = document.getElementById("ffErrorsBlock");
			var errorsgroup = document.getElementById("ffErrors");
			if (errorsblock.style.display == "none") {
				errorsblock.style.display = "block";
				errorsgroup.style.right = "10px";
			} else {
				errorsblock.style.display = "none";
				errorsgroup.style.right = "300px";
			}
		}
//]]></script>
<div id="ffErrors">
<a href="javascript:toggle_ffErrors();">Click to toggle</a>
<div id="ffErrorsBlock">
	<div class="error">does not end with &lt;/html&gt; tag</div>
	<div class="error">does not end with &lt;/body&gt; tag</div>
	<div class="info">The output has ended thus: }

			stabilize;
		}
	}
}

initially [ LIC(xi, yi) | yi in 0..(imgSizeY-1), xi in 0..(imgSizeX-1) ];</div>
</div></div>
</body></html>