Material - SimpleOcean.osl

// SimpleOcean by Philippe Groarke
// Copyright 2024 Philippe Groarke, All rights reserved. This file is licensed under Apache 2.0 license
//    https://github.com/ADN-DevTech/3dsMax-OSL-Shaders/blob/master/LICENSE.txt

/*
TODO :
- debug cell noise
- test vray displace
- figure direction utility map
- debugging


Help

Arnold Properties
- General > Opaque : off
- Displacement > Enable
- Displacement > Bounds pads : 1
- Displacement > Use Map
- Connect Displacement output to that map.

Optional, Auto-Bump can replace subdiv and is much faster.
- Subdivision > Catclark
- Subdivision > Iterations : ~2-3

Arnold Standard Surface Setup
- Diffuse : 0
- IOR : 1.33
- Specular : ~1
- Transmission : 1
- Transmission Color : ~blue
- Transmission Depth : 0.1+
- Scatter : White (or other color)
- Scatter Anisotropy : 1
*/

#define MAX_OCTAVES 12
#define MAX_WAVES 32
#define ARR_SIZE 384

struct wave_info {
	float wavelength;
	float amp;
	vector dir;
	vector offset;

	// computed
	float w;
	float q;
	float phi;
	float cos_term;
	float sin_term;
};

float wave_w(float wavelength) {
	return 2.0 / wavelength;
}

float wave_q(float steepness, float w, float amp, float num_waves) {
	return steepness / (amp * w * num_waves);
}

float wave_phi(float speed, float wavelength) {
	return speed * sqrt(9.8 * (M_2PI / wavelength));
}

float max_height(int num_octaves, int num_waves, wave_info waves[ARR_SIZE]) {
	float ret = 0.0;
	for (int j = 0; j < num_octaves; ++j) {
		int count = j * num_waves + num_waves;
		for (int i = j * num_waves; i < count; ++i) {
			ret += waves[i].amp;
		}
	}
	return ret / 2.0;
}

float fract(float x) {
	return x - floor(x);
}

// https://www.shadertoy.com/view/Xt23Ry
float rand(float co) {
	return fract(sin(co * 91.3458) * 47453.5453);
}

float remap(float val, float low1, float high1, float low2, float high2) {
	return low2 + (val - low1) * (high2 - low2) / (high1 - low1);
}

void compute_terms(point pos, float t, int num_octaves, int num_waves, output wave_info waves[ARR_SIZE]) {
	for (int j = 0; j < num_octaves; ++j) {
		int count = j * num_waves + num_waves;
		for (int i = j * num_waves; i < count; ++i) {
			float term = waves[i].w * dot(waves[i].dir, pos + waves[i].offset) + waves[i].phi * t;

			waves[i].cos_term = cos(term);
			waves[i].sin_term = sin(term);
		}
	}
}

vector compute_displacement(int num_octaves, int num_waves, wave_info waves[ARR_SIZE]) {
	vector ret = 0;
	for (int j = 0; j < num_octaves; ++j) {
		int count = j * num_waves + num_waves;
		for (int i = j * num_waves; i < count; ++i) {
			float xy_term = waves[i].q * waves[i].amp * waves[i].cos_term;

			ret[0] += xy_term * waves[i].dir[0];
			ret[1] += xy_term * waves[i].dir[1];
			ret[2] += waves[i].amp * waves[i].sin_term;
		}
	}

	return ret;
}

vector compute_normal(int num_octaves, int num_waves, wave_info waves[ARR_SIZE]) {
	vector ret = 0;
	for (int j = 0; j < num_octaves; ++j) {
		int count = j * num_waves + num_waves;
		for (int i = j * num_waves; i < count; ++i) {
			float xy_term = waves[i].w * waves[i].amp * waves[i].cos_term;

			ret[0] += waves[i].dir[0] * xy_term;
			ret[1] += waves[i].dir[1] * xy_term;
			ret[2] += waves[i].q * waves[i].w * waves[i].amp * waves[i].sin_term;
		}
	}

	ret[0] = -ret[0];
	ret[1] = -ret[1];
	ret[2] = 1.0 - ret[2];
	return ret;
}

// Dx
vector compute_bitangent(int num_octaves, int num_waves, wave_info waves[ARR_SIZE]) {
	vector ret = 0;
	for (int j = 0; j < num_octaves; ++j) {
		int count = j * num_waves + num_waves;
		for (int i = j * num_waves; i < count; ++i) {
			float xy_term = waves[i].q * waves[i].w * waves[i].amp * waves[i].sin_term;

			ret[0] += pow(waves[i].dir[0], 2.0) * xy_term;
			ret[1] += waves[i].dir[0] * waves[i].dir[1] * xy_term;
			ret[2] += waves[i].dir[0] * waves[i].w * waves[i].amp * waves[i].cos_term;
		}
	}

	ret[0] = 1.0 - ret[0];
	ret[1] = -ret[1];
	return ret;
}

// Dy
vector compute_tangent(int num_octaves, int num_waves, wave_info waves[ARR_SIZE]) {
	vector ret = 0;
	for (int j = 0; j < num_octaves; ++j) {
		int count = j * num_waves + num_waves;
		for (int i = j * num_waves; i < count; ++i) {
			float xy_term = waves[i].q * waves[i].w * waves[i].amp * waves[i].sin_term;

			ret[0] += waves[i].dir[0] * waves[i].dir[1] * xy_term;
			ret[1] += pow(waves[i].dir[1], 2.0) * xy_term;
			ret[2] += waves[i].dir[1] * waves[i].w * waves[i].amp * waves[i].cos_term;
		}
	}

	ret[0] = -ret[0];
	ret[1] = 1.0 - ret[1];
	return ret;
}

// Play with k,
// https://www.desmos.com/calculator/og836nvwmx
float interp(float k, float percent) {
	float ret = 0.0;
	float epsilon = 0.0001;
	if (fabs(k) < epsilon) {
		// Actual k = 0 == 0 always. Just do linear interp.
		ret = percent;
	} else {
		ret = (exp(k * percent) - 1.0) / (exp(k) - 1.0);
	}
	return ret;
}

shader SimpleOcean
[[
	string help  =
		"<h3>Simple Ocean</h3>"
		"A vector displacement shader to generate basic ocean waves.<br>"
		"<h4>Instructions</h4>"
		"The shader operates in Tangent Space. To get any useful results, use a 'UVW Map' modifier and apply 'XYZ to UVW' to your target object.<br>"
		"The output Displacement will only work with 'Arnold Properties' modifier vector displacement input.",
	string label = "Material - Simple Ocean"
]]
(
	int NumOctaves = 7
	[[
		string help = "Wave batches are generated per octave, going down from the wave length specified.",
		int min = 1,
		int max = MAX_OCTAVES,
		int connectable = 0
	]],

	int NumWaves = 8
	[[
		string help = "Number of waves to generate, per octave.",
		int min = 1,
		int max = MAX_WAVES,
		int connectable = 0
	]],

	int UVSet = 1
	[[
		string label="UV Map Channel",
		string help ="The 3ds max Map channel to look up. Must have been primed with 'XYZ to UVW'.",
		int min = 0,
		int max = 99
	]],

	vector UVOffset = 0.0
	[[
		string help = "The UV lookup offset. Can be used to move waves.",
		int connectable = 0
	]],

	string Randomness = "Hash 1"
	[[
		string widget= "popup",
		string help  = "Various noises used to gather entropy.",
		string options="Hash 1|Hash 2|Hash 3|Randy Mc. Randface|Cell 1|Rand 1|Rand 2|Rand 3|Perlin|Even Distribution|Simplex",
		int connectable = 0
	]],

	int Seed = 0
	[[
		string help  = "Randomness seed. Doesn't affect all random engines.",
		int connectable = 0
	]],

	float WaveLength = 40.0
	[[
		string help = "The largest wave length. Waves will be generated starting with this value.",
		int connectable = 0
	]],

	float Amplitude = 1.0
	[[
		string help = "The height of the waves.",
		int connectable = 0
	]],

	float AmplitudeDecay = 5.0
	[[
		string help = "Amplitude is scaled according to the octave frequencies. You can increase the scaling rate (making smaller waves less tall).",
		float min = 0.0,
		int connectable = 0
	]],

	int RandomizeAmplitude = 0
	[[
		string help  = "Randomize amplitude with values around your provided amplitude.",
		string widget= "checkBox",
		int connectable = 0
	]],

	float Steepness = 7.0
	[[
		string help = "How spiky the wave is. Too much steepness will cause crests to curl." \
			" Depending on your material settings, this can result in a fake foam look.",
		int connectable = 0
	]],

	float SteepnessDecay = 5.0
	[[
		string help = "Steepness is scaled according to the octave frequencies. You can increase or decrease the steepness scaling if you wish.",
		float min = 0.0,
		// float max = 1.0,
		int connectable = 0
	]],

	int RandomizeSteepness = 0
	[[
		string help  = "Randomize Steepness with values around your provided steepness.",
		string widget= "checkBox",
		int connectable = 0
	]],

	float Direction = 45.0
	[[
		string help = "Direction of the waves in degrees. 0 degrees is +x.",
		string units = "degrees",
		float min = 0.0,
		float max = 360.0,
		int connectable = 0
	]],

	float Deviation = 120.0
	[[
		string help = "The max angle to deviate from direction, in degrees.",
		string units = "degrees",
		float min = 0.0,
		float max = 360.0,
		int connectable = 0
	]],

	float DeviationGrowth = 5.0
	[[
		string help = "This parameter affects Deviation increase of smaller waves. "
		"It helps reduce banding and tiling, while keeping the general wave direction intact.",
		int connectable = 0
	]],

	float Speed = 1.0
	[[
		string help = "Wave speed multiplier.",
		int connectable = 0
	]],

	float Time = 0.0
	[[
		string label = "Time (s)",
		float timeValue = 0.0
	]],

	// Sea Level Settings
	float SeaLevel = 25.0
	[[
		string help = "Increase or decrease the Sea Level height."
	]],

	int SeaLevelRamp = 1
	[[
		string help  = "Use interpolation to smooth the sea level map.",
		string widget= "checkBox",
		int connectable = 0
	]],

	float SeaLevelRampSize = 25.0
	[[
		string help = "How wide is the ramp from under sea level to above sea level. Only valid when using SeaLevelRamp.",
		float min = 0.1,
		int connectable = 0
	]],

	float SeaLevelRampSteepness = 5.0
	[[
		string help = "How steep is the sea level ramp. Positive is exponential, negative is logarithmic and '0' is linear interpolation.",
		int connectable = 0
	]],


	// Foam Settings
	float FoamThreshold = 0.2
	[[
		string help = "Controls how many peaks are considered foam. Higher values give result in more foam.",
		float min = 0.0
	]],

	int FoamRamp = 1
	[[
		string help  = "Use interpolation to smooth the foam map.",
		string widget = "checkBox",
		int connectable = 0
		// string label = "Sea Level Ramp"
	]],

	float FoamRampSize = 0.4
	[[
		string help = "How wide is the ramp from no foam to all foam. Only valid when using FoamRamp.",
		float min = 0.1,
		int connectable = 0
	]],

	float FoamRampSteepness = 5.0
	[[
		string help = "How steep is the foam ramp. Positive is exponential, negative is logarithmic and '0' is linear interpolation.",
		int connectable = 0
	]],

	// ModulatedFoam Settings
	float FoamTwoThreshold = 0.5
	[[
		string help = "Controls how many peaks are considered foam, for Modulated Foam Map. Higher values give result in more foam.",
		float min = 0.0
	]],

	float FoamTwoBrightness = 2.0
	[[
		string help = "Affects the FoamMapTwo output brightness.",
		float min = 0.0
	]],

	// int Debug = 0
	// [[
	// 	string help = "",
	// 	string widget = "checkBox",
	// 	string label = "Debug",
	// 	int connectable = 0
	// ]],


	// Outputs
	output vector Displacement = 0
	[[
		string help = "The vector displacement. This must be used with 'Arnold Properties' vector displacement input."
	]],

	output color SeaLevelMap = 0
	[[
		string help = "Outputs a utility map which colors white above provided sea level."
	]],

	output color BreakMap = 0
	[[
		string help = "Outputs a utility map which identifies where peaks 'break' or 'curl'. You can use this to hide curls when using extreme steepness."
	]],

// Disable for now. Needs more work and investigation.
#if 0
	output vector BreakDirection = 0
	[[
		string help = "EXPERIMENTAL : The local wave break direction."
	]],
#endif

	output color FoamMap = 0
	[[
		string help = "Outputs a utility map which identifies wave peaks, according to the provided parameters."
	]],

	output color FoamMapTwo = 0
	[[
		string help = "An alternative Foam map with a different smooth algorithm."
	]]

	// output color DebugMap = 0
)
{

	// Screw it.
	float randy_mcrandface[ARR_SIZE] = {
		0.7613283154, 0.6417486594, 0.7077160716, 0.2767145632, 0.3459236354, 0.4595805590,
		0.4221504986, 0.5430591227, 0.7397216781, 0.6009534408, 0.7567760852, 0.4356619047,
		0.7305973256, 0.8078904319, 0.9670611668, 0.8112846347, 0.8832108729, 0.2371557390,
		0.9638971629, 0.1849049127, 0.5075650226, 0.6080068042, 0.7992122890, 0.9957855124,
		0.7726620573, 0.5271973617, 0.6167741389, 0.1633856410, 0.5368961641, 0.9393507032,
		0.8341412603, 0.1431013143, 0.1222931448, 0.9149053562, 0.6667115294, 0.7362864532,
		0.6042048595, 0.1710390230, 0.6100750959, 0.9878706709, 0.9662866655, 0.4458209901,
		0.8420691186, 0.0717683244, 0.9372985664, 0.2519136408, 0.3386565954, 0.5898405778,
		0.7636390582, 0.0590747309, 0.3200221082, 0.8101507442, 0.0922381526, 0.1751451629,
		0.8260694506, 0.4569933524, 0.9231153214, 0.2806759626, 0.8814646580, 0.2871386945,
		0.1294133616, 0.1840314361, 0.7272801375, 0.5209047709, 0.4468515935, 0.0739586178,
		0.7808635264, 0.1385085642, 0.5818542706, 0.1419560811, 0.4914149198, 0.3984367663,
		0.5534051739, 0.9021820860, 0.0457077229, 0.9430585064, 0.5832536408, 0.9134246462,
		0.5584909940, 0.0115561512, 0.5102350494, 0.0843100883, 0.6761218329, 0.8182395543,
		0.0778287485, 0.7653847158, 0.5508998159, 0.2821854233, 0.2219167320, 0.9900815780,
		0.0790549319, 0.9048603125, 0.0935980279, 0.9857482777, 0.8467555816, 0.1699877931,
		0.0481672073, 0.2302571620, 0.4959770668, 0.5960198219, 0.4312765748, 0.1591156525,
		0.0687268269, 0.4639244974, 0.0731531983, 0.2757107353, 0.6317098928, 0.2759345973,
		0.0107261248, 0.8062241995, 0.3991566052, 0.4470180022, 0.4447787828, 0.7814391798,
		0.3507617008, 0.2148204368, 0.0951097680, 0.7792412681, 0.6008316251, 0.4973092279,
		0.1175454673, 0.8624851465, 0.3561905606, 0.7163313779, 0.1186382648, 0.9134720187,
		0.6145221594, 0.1339072762, 0.3994272383, 0.9292100772, 0.8215033421, 0.6260928175,
		0.1578142001, 0.6831931992, 0.4945404445, 0.3332440578, 0.2528690529, 0.7808411693,
		0.7977090604, 0.6156838816, 0.9647099684, 0.2205307551, 0.1242177648, 0.5828725769,
		0.6933682418, 0.7172738086, 0.7519134569, 0.2997492017, 0.3270603535, 0.0220112573,
		0.0187798525, 0.3919452145, 0.5571065984, 0.9738236155, 0.6471761395, 0.7259059581,
		0.0263854616, 0.8480789135, 0.4626754320, 0.0731667099, 0.8714149142, 0.0847650250,
		0.0097162705, 0.0560917960, 0.8960117455, 0.6216735959, 0.8982684282, 0.8750739734,
		0.1650972796, 0.3906805203, 0.3209709014, 0.4877145134, 0.0306574656, 0.5971331151,
		0.6189258743, 0.4184951607, 0.8335943660, 0.1242530933, 0.8798338352, 0.5833763266,
		0.8094074573, 0.5254715167, 0.3668876261, 0.1422191387, 0.1721910074, 0.2026817602,
		0.2119211370, 0.6253000589, 0.4904723281, 0.5218543811, 0.9301133550, 0.9683309118,
		0.3661083587, 0.2307221251, 0.3787893454, 0.4751772670, 0.0579906878, 0.1768354991,
		0.8891576681, 0.1431432921, 0.4797410106, 0.2663625361, 0.8325811059, 0.5140320869,
		0.7597127482, 0.0470517628, 0.5803547234, 0.3548459410, 0.2100508372, 0.5428815615,
		0.1411484061, 0.3602227092, 0.7215295438, 0.5434041909, 0.1822010627, 0.8677546826,
		0.7585891515, 0.9008217104, 0.6470395283, 0.8681916740, 0.4534938072, 0.3995711447,
		0.8564553250, 0.0570914895, 0.6910691710, 0.5583367680, 0.5398101111, 0.7521228693,
		0.7568702395, 0.8925501029, 0.6770392693, 0.9093505238, 0.9102231200, 0.2661670667,
		0.1226216037, 0.8111049737, 0.2353793834, 0.7364151392, 0.4747462044, 0.9299663536,
		0.8123320081, 0.6555540549, 0.7453898179, 0.3053406749, 0.3159677258, 0.7146935034,
		0.4320777924, 0.8925243566, 0.1110892016, 0.4212926955, 0.8471286552, 0.4607157468,
		0.5123215652, 0.1875337625, 0.7745056160, 0.3593490972, 0.5003064330, 0.9421407414,
		0.6629999877, 0.1430016877, 0.3685406819, 0.5161234146, 0.8825954808, 0.8835394491,
		0.8222370408, 0.2594687195, 0.4559979994, 0.5882772874, 0.0754071303, 0.1083060122,
		0.7392239567, 0.6197596061, 0.5819020307, 0.7951468252, 0.7117926842, 0.2702588942,
		0.1472675702, 0.0637894921, 0.9838483921, 0.7359211919, 0.0214477716, 0.9105100154,
		0.3542747739, 0.6448500090, 0.9145919663, 0.4144151693, 0.0563616620, 0.8628780555,
		0.4093985005, 0.9681249187, 0.5956522112, 0.4876563354, 0.2623010908, 0.2294213032,
		0.1300411096, 0.8833217846, 0.1880329714, 0.1318175852, 0.4766515362, 0.7480518772,
		0.6049902499, 0.5670014017, 0.8094345290, 0.0488128479, 0.3806065262, 0.5129413053,
		0.8827513878, 0.3742173559, 0.5897817464, 0.1265311938, 0.1666639114, 0.3513735293,
		0.0268498170, 0.5559057612, 0.2938191986, 0.4754748694, 0.1886570996, 0.1506464087,
		0.8847426608, 0.4884451589, 0.4350855332, 0.9345445426, 0.4099686316, 0.7942126964,
		0.9765494371, 0.0297738076, 0.0198415258, 0.3841140226, 0.6416234156, 0.9758983813,
		0.6127427981, 0.7106013024, 0.1153512188, 0.6440086521, 0.9619477967, 0.4960318780,
		0.1374903689, 0.0664052156, 0.3181687044, 0.5267333783, 0.3895587626, 0.3595304375,
		0.0330008479, 0.6608198678, 0.7672469086, 0.3995902480, 0.7524112037, 0.7605942560,
		0.1179058104, 0.0814233157, 0.7398095277, 0.4125950086, 0.8702069902, 0.4078816518,
		0.1060588817, 0.5912999079, 0.3719537622, 0.4171108244, 0.3215406658, 0.1411436386,
		0.5832609558, 0.4554747267, 0.0141743945, 0.4657320091, 0.6624735639, 0.5576015787,
		0.7421591684, 0.3185219243, 0.6550715958, 0.3824203749, 0.2708995102, 0.7972102891,
		0.3792130014, 0.2413789324, 0.3362235952, 0.4771132808, 0.7399085797, 0.0936726877,
		0.1454969690, 0.7006940519, 0.5992551756, 0.2845851767, 0.4452032488, 0.5779639096
	};

	// Array of waves settings.
	wave_info waves[ARR_SIZE];

	vector right = { 1.0, 0.0, 0.0 };
	int total_waves = NumOctaves * NumWaves;

	float current_wavelength = WaveLength;
	for (int j = 0; j < NumOctaves; ++j) {
		float min_wavelength = current_wavelength / 1.5;
		float max_wavelength = current_wavelength * 1.5;

		int count = j * NumWaves + NumWaves;
		for (int i = j * NumWaves; i < count; ++i) {
			float abs_idx = float(i + Seed);
			float local_idx = float((i + Seed) % NumWaves);
			float norm_idx = local_idx / float(NumWaves);

			float abs_seed = float(Seed);
			float local_seed = float(Seed % NumWaves);
			float norm_seed = local_seed / float(NumWaves);

			float rng = 0.0;
			if (Randomness == "Hash 1") {
				rng = noise("hash", local_idx, local_seed);
			} else if (Randomness == "Hash 2") {
				rng = noise("hash", norm_idx, norm_seed + norm_idx);
			} else if (Randomness == "Hash 3") {
				rng = noise("hash", abs_idx, abs_seed);
			} else if (Randomness == "Randy Mc. Randface") {
				rng = randy_mcrandface[(i + Seed) % ARR_SIZE];
			} else if (Randomness == "Cell 1") {
				rng = noise("cell", abs_idx, abs_seed);
			} else if (Randomness == "Rand 1") {
				rng = rand(abs_idx);
			} else if (Randomness == "Rand 2") {
				rng = rand(norm_idx);
			} else if (Randomness == "Rand 3") {
				rng = rand(local_idx);
			} else if (Randomness == "Perlin") {
				rng = noise("uperlin", norm_idx, norm_seed);
			} else if (Randomness == "Even Distribution") {
				rng = norm_idx;
			} else if (Randomness == "Simplex") {
				rng = noise("usimplex", norm_idx, norm_seed);
			}
			// too heavy :/ provide noise map input?
			// float rng = fabs(noise("gabor", norm_idx));

			// rng = mod(rng + randy_mcrandface[i], 1.0);
			rng = mod(rng, 1.0);

			float wlength = mix(min_wavelength, max_wavelength, rng);
			float decay_per = float(j) / float(NumOctaves);

			// As waves get smaller, make their amplitude smaller.
			float amp_decay = 1.0 / (1.0 + decay_per * AmplitudeDecay);

			// Add randomization to amplitude.
			float in_amp = Amplitude;
			if (RandomizeAmplitude) {
				in_amp *= (rng * 2.0);
			}

			// Amp 1 ~= Wavelength
			float amp = (wlength * in_amp * amp_decay) / NumWaves;

			float w = wave_w(wlength);

			// As waves get smaller, so does decay.
			float steep_decay = 1.0 / (1.0 + decay_per * SteepnessDecay);

			// Add randomization to steepness.
			float in_steepness = Steepness;
			if (RandomizeSteepness) {
				in_steepness *= (rng * 2.0);
			}

			float steep = in_steepness * steep_decay;
			float q = wave_q(steep, w, amp, float(total_waves));

			float phi = wave_phi(Speed, wlength);

			// The angle is randomly rotated between the allowed range.
			// The smaller the waves, the higher the angle range (to reduce banding).
			float angle = Direction;
			float dev_growth = 1.0 + decay_per * DeviationGrowth;
			float dev = Deviation * dev_growth;

			angle += mix(0.0, dev, rng) - (dev / 2.0);
			vector my_dir = rotate(right, radians(angle), vector(0.0, 0.0, 1.0));

			vector my_offset = vector(0,0,0);
			float off = i * wlength + rng * wlength;
			my_offset = vector(off, off, 0.0);

			waves[i].wavelength = wlength;
			waves[i].amp = amp;
			waves[i].dir = my_dir;
			waves[i].offset = my_offset;
			waves[i].w = w;
			waves[i].q = q;
			waves[i].phi = phi;
		}

		current_wavelength /= 2.0;
	}



	// Position to compute comes from uv coords, aka tangent space xy.
	point pos = point(u, v, 0.0);

	// Get UVW channel
	if (UVSet != 1) {
		if (!getattribute(format("UV%d", UVSet), pos))
		{
			// Default UV set handled in
			// a special way in MAXtoA.
			// But we don't care since w is ignored.
			getattribute(format("uv_%d", UVSet), pos);
		}
	}

	// Just in case, reset w.
	pos[0] += UVOffset[0];
	pos[1] += UVOffset[1];
	pos[2] = 0.0;

	compute_terms(pos, Time, NumOctaves, NumWaves, waves);

	Displacement = compute_displacement(NumOctaves, NumWaves, waves);

	// TODO : bump()?
	// vector Normal = compute_normal(NumOctaves, NumWaves, waves);

	// Trivial over/under z.
	{
		float height = Displacement[2];
		if (height >= SeaLevel) {
			SeaLevelMap = color(1);
		}

		if (SeaLevelRamp) {
			float ramp_start = SeaLevel - SeaLevelRampSize;

			if (height >= ramp_start && height < SeaLevel) {
				float per = (height - ramp_start) / SeaLevelRampSize;
				float r = interp(SeaLevelRampSteepness, per);
				SeaLevelMap = mix(color(0), color(1), r);
			}
		}
	}


	// Jacobian crest detect
	{
		vector Bitangent = compute_bitangent(NumOctaves, NumWaves, waves);
		vector Tangent = compute_tangent(NumOctaves, NumWaves, waves);
		float Jxx = Bitangent[0];
		float Jxy = Bitangent[1];
		float Jyy = Tangent[1];

		float Jmin = 0.5 * (Jxx + Jyy) - 0.5 * sqrt(pow(Jxx - Jyy, 2.0) + 4.0 * pow(Jxy, 2.0));
		float Jplus = 0.5 * (Jxx + Jyy) + 0.5 * sqrt(pow(Jxx - Jyy, 2.0) + 4.0 * pow(Jxy, 2.0));
		float J = Jmin * Jplus;

// The paper states you can only check with Jmin, but its unclear to me how
// that interacts with the ramp stuff.
#if 0
		if (Jmin < 0.0) {
			BreakMap = color(1,1,1);

			float qmin = (Jmin - Jxx) / Jxy;
			vector emin = vector(1, qmin, 0.0) / sqrt(1 + pow(qmin, 2.0));
			BreakDirection = emin;
		}

		if (Jmin < FoamThreshold) {
			FoamMap = color(1,1,1);
		}

		if (FoamRamp) {
			float ramp_start = FoamThreshold + FoamRampSize;

			if (FoamThreshold < Jmin && Jmin <= ramp_start) {
				float per = 1.0 - ((Jmin - FoamThreshold) / FoamRampSize);
				float r = interp(FoamRampSteepness, per);
				FoamMap = mix(color(0), color(1), r);
			}
		}
#else

// What does this actually give us? What is the use case?
#if 0
		// if (J < 0.0) {
			// BreakMap = color(1,1,1);

			float qmin = (Jmin - Jxx) / Jxy;
			vector emin = vector(1, qmin, 0.0) / sqrt(1 + pow(qmin, 2.0));
			BreakDirection = emin;
		// }
#endif

		if (J < 0.0) {
			BreakMap = color(1,1,1);
		}

		if (J < FoamThreshold) {
			FoamMap = color(1,1,1);
		}

		if (FoamRamp) {
			float ramp_start = FoamThreshold + FoamRampSize;

			if (FoamThreshold < J && J <= ramp_start) {
				float per = 1.0 - ((J - FoamThreshold) / FoamRampSize);
				float r = interp(FoamRampSteepness, per);
				FoamMap = mix(color(0), color(1), r);
			}
		}

		// This is another way than ramp to have nicer foam.
		{
			float m_col = FoamTwoBrightness * (-J + FoamTwoThreshold);
			FoamMapTwo = clamp(color(m_col, m_col, m_col), color(0,0,0), color(1,1,1));
		}


#endif

		// TODO : Experiment with bending curls.
	}
}