// Copyright (c) 2023 Lachlan McDonald
// This work is licensed under the MIT License (MIT)
// https://github.com/lachlanmcdonald/magicavoxel-shaders
//
// This script utilises or modifies code from other projects or publications.
// Please see the attributions below for more information:
//
// 1. Copyright (c) 2011 Stefan Gustavson <https://github.com/ashima/webgl-noise>
//    MIT License (MIT)
//    https://github.com/ashima/webgl-noise/blob/master/LICENSE
//
// xs noise/cellular3D [Mode] [Scale] [Jitter] [Noise] [Power] [Cavity] [Seed] [Tile X] [Tile Y] [Tile Z]
//
// xs_begin
// author : '@lmcdx.bsky.social'
// arg : { name = 'Mode'  var = 'm_mode'  range = '0 5'  value = '0'  step = '1'  precision = '0' }
// arg : { name = 'Scale'  var = 'm_scale'  range = '1 100'  value = '24'  step = '1'  precision = '0' }
// arg : { name = 'Jitter'  var = 'm_jitter'  range = '-100 100'  value = '100'  step = '5'  precision = '0' }
// arg : { name = 'Noise'  var = 'm_noise'  range = '0 100'  value = '0'  step = '1'  precision = '0' }
// arg : { name = 'Power'  var = 'm_power'  range = '0 10'  value = '1'  step = '0.05'  precision = '2' }
// arg : { name = 'Cavity'  var = 'm_cavity'  range = '0 100'  value = '50'  step = '1'  precision = '0' }
// arg : { name = 'Seed'  var = 'global_seed'  range = '1 100'  value = '1'  step = '1'  precision = '0' }
// arg : { name = 'Tile X'  var = 'tile_x'  range = '0 40'  value = '0'  step = '1'  precision = '0' }
// arg : { name = 'Tile Y'  var = 'tile_y'  range = '0 40'  value = '0'  step = '1'  precision = '0' }
// arg : { name = 'Tile Z'  var = 'tile_z'  range = '0 40'  value = '0'  step = '1'  precision = '0' }
// xs_end

int mode = int(m_mode);
float noise = m_noise / 100.0;
float jitter = m_jitter / 100.0;
float cavity = m_cavity / 100.0;
vec3 tile = vec3(tile_x, tile_y, tile_z) * i_volume_size;

float pal(float p) {
	float f = floor(mix(0.0, float(i_num_color_sels), p));
	return color_sel(f);
}

vec3 mod289(vec3 x) {
	return x - floor(x * (1.0 / 289.0)) * 289.0;
}

vec3 mod7(vec3 x) {
	return x - floor(x * (1.0 / 7.0)) * 7.0;
}

vec3 permute(vec3 x) {
	return mod289((34.0 * x + 10.0) * x);
}

vec2 cellular(vec3 co, float jitter) {
	float K = 0.142857142857;
	float Ko = 0.428571428571;
	float K2 = 0.020408163265306;
	float Kz = 0.166666666667;
	float Kzo = 0.416666666667;

	vec3 Pi = mod289(floor(co));
 	vec3 Pf = fract(co) - 0.5;

	vec3 Pfx = Pf.x + vec3(1.0, 0.0, -1.0);
	vec3 Pfy = Pf.y + vec3(1.0, 0.0, -1.0);
	vec3 Pfz = Pf.z + vec3(1.0, 0.0, -1.0);

	vec3 p = permute(Pi.x + vec3(-1.0, 0.0, 1.0));
	vec3 p1 = permute(p + Pi.y - 1.0);
	vec3 p2 = permute(p + Pi.y);
	vec3 p3 = permute(p + Pi.y + 1.0);

	vec3 p11 = permute(p1 + Pi.z - 1.0);
	vec3 p12 = permute(p1 + Pi.z);
	vec3 p13 = permute(p1 + Pi.z + 1.0);

	vec3 p21 = permute(p2 + Pi.z - 1.0);
	vec3 p22 = permute(p2 + Pi.z);
	vec3 p23 = permute(p2 + Pi.z + 1.0);

	vec3 p31 = permute(p3 + Pi.z - 1.0);
	vec3 p32 = permute(p3 + Pi.z);
	vec3 p33 = permute(p3 + Pi.z + 1.0);

	vec3 ox11 = fract(p11 * K) - Ko;
	vec3 oy11 = mod7(floor(p11 * K)) * K - Ko;
	vec3 oz11 = floor(p11 * K2) * Kz - Kzo;

	vec3 ox12 = fract(p12 * K) - Ko;
	vec3 oy12 = mod7(floor(p12 * K)) * K - Ko;
	vec3 oz12 = floor(p12 * K2) * Kz - Kzo;

	vec3 ox13 = fract(p13 * K) - Ko;
	vec3 oy13 = mod7(floor(p13 * K)) * K - Ko;
	vec3 oz13 = floor(p13 * K2) * Kz - Kzo;

	vec3 ox21 = fract(p21 * K) - Ko;
	vec3 oy21 = mod7(floor(p21 * K)) * K - Ko;
	vec3 oz21 = floor(p21 * K2) * Kz - Kzo;

	vec3 ox22 = fract(p22 * K) - Ko;
	vec3 oy22 = mod7(floor(p22 * K)) * K - Ko;
	vec3 oz22 = floor(p22 * K2) * Kz - Kzo;

	vec3 ox23 = fract(p23 * K) - Ko;
	vec3 oy23 = mod7(floor(p23 * K)) * K - Ko;
	vec3 oz23 = floor(p23 * K2) * Kz - Kzo;

	vec3 ox31 = fract(p31 * K) - Ko;
	vec3 oy31 = mod7(floor(p31 * K)) * K - Ko;
	vec3 oz31 = floor(p31 * K2) * Kz - Kzo;

	vec3 ox32 = fract(p32 * K) - Ko;
	vec3 oy32 = mod7(floor(p32 * K)) * K - Ko;
	vec3 oz32 = floor(p32 * K2) * Kz - Kzo;

	vec3 ox33 = fract(p33 * K) - Ko;
	vec3 oy33 = mod7(floor(p33 * K)) * K - Ko;
	vec3 oz33 = floor(p33 * K2) * Kz - Kzo;

	vec3 dx11 = Pfx + jitter * ox11;
	vec3 dy11 = Pfy.x + jitter * oy11;
	vec3 dz11 = Pfz.x + jitter * oz11;

	vec3 dx12 = Pfx + jitter * ox12;
	vec3 dy12 = Pfy.x + jitter * oy12;
	vec3 dz12 = Pfz.y + jitter * oz12;

	vec3 dx13 = Pfx + jitter * ox13;
	vec3 dy13 = Pfy.x + jitter * oy13;
	vec3 dz13 = Pfz.z + jitter * oz13;

	vec3 dx21 = Pfx + jitter * ox21;
	vec3 dy21 = Pfy.y + jitter * oy21;
	vec3 dz21 = Pfz.x + jitter * oz21;

	vec3 dx22 = Pfx + jitter * ox22;
	vec3 dy22 = Pfy.y + jitter * oy22;
	vec3 dz22 = Pfz.y + jitter * oz22;

	vec3 dx23 = Pfx + jitter * ox23;
	vec3 dy23 = Pfy.y + jitter * oy23;
	vec3 dz23 = Pfz.z + jitter * oz23;

	vec3 dx31 = Pfx + jitter * ox31;
	vec3 dy31 = Pfy.z + jitter * oy31;
	vec3 dz31 = Pfz.x + jitter * oz31;

	vec3 dx32 = Pfx + jitter * ox32;
	vec3 dy32 = Pfy.z + jitter * oy32;
	vec3 dz32 = Pfz.y + jitter * oz32;

	vec3 dx33 = Pfx + jitter * ox33;
	vec3 dy33 = Pfy.z + jitter * oy33;
	vec3 dz33 = Pfz.z + jitter * oz33;

	vec3 d11 = dx11 * dx11 + dy11 * dy11 + dz11 * dz11;
	vec3 d12 = dx12 * dx12 + dy12 * dy12 + dz12 * dz12;
	vec3 d13 = dx13 * dx13 + dy13 * dy13 + dz13 * dz13;
	vec3 d21 = dx21 * dx21 + dy21 * dy21 + dz21 * dz21;
	vec3 d22 = dx22 * dx22 + dy22 * dy22 + dz22 * dz22;
	vec3 d23 = dx23 * dx23 + dy23 * dy23 + dz23 * dz23;
	vec3 d31 = dx31 * dx31 + dy31 * dy31 + dz31 * dz31;
	vec3 d32 = dx32 * dx32 + dy32 * dy32 + dz32 * dz32;
	vec3 d33 = dx33 * dx33 + dy33 * dy33 + dz33 * dz33;

	vec3 d1a = min(d11, d12);
	d12 = max(d11, d12);
	d11 = min(d1a, d13);
	d13 = max(d1a, d13);
	d12 = min(d12, d13);
	vec3 d2a = min(d21, d22);
	d22 = max(d21, d22);
	d21 = min(d2a, d23);
	d23 = max(d2a, d23);
	d22 = min(d22, d23);
	vec3 d3a = min(d31, d32);
	d32 = max(d31, d32);
	d31 = min(d3a, d33);
	d33 = max(d3a, d33);
	d32 = min(d32, d33);
	vec3 da = min(d11, d21);
	d21 = max(d11, d21);
	d11 = min(da, d31);
	d31 = max(da, d31);
	d11.xy = (d11.x < d11.y) ? d11.xy : d11.yx;
	d11.xz = (d11.x < d11.z) ? d11.xz : d11.zx;
	d12 = min(d12, d21);
	d12 = min(d12, d22);
	d12 = min(d12, d31);
	d12 = min(d12, d32);
	d11.yz = min(d11.yz,d12.xy);
	d11.y = min(d11.y,d12.z);
	d11.y = min(d11.y,d11.z);

	return sqrt(d11.xy);
}

float map(vec3 v) {
	vec3 offset = global_seed * i_volume_size + tile;
	vec3 k = floor(v) / i_volume_size;
	vec2 p = cellular((v + offset) / m_scale, jitter);
	vec2 q = cellular((v + offset), 1.0);
	float r;
	float z;

	if (mode <= 2) {
		r = clamp(pow(p.x + q.x * noise, m_power), 0.0, 1.0);
	} else {
		r = clamp(pow(p.y + q.y * noise, m_power), 0.0, 1.0);
	}

	if (mode == 1 || mode == 4) {
		z = r;
	}  else if (mode == 2 || mode == 5) {
		z = 1.0 - r;
	} else {
		z = 1.0;
	}

	if (mode == 0 || mode == 3) {
		return (r > cavity) ? pal((r - cavity) / (1.0 - cavity)) : 0.0;
	} else {
		return (k.z <= z) ? pal(r) : 0.0;
	}
}