using System;
using System.IO;
using System.Text;
using System.Runtime.InteropServices;

namespace CharacterEditor;
/// <summary>Fits loaded human clothes against shared body surfaces, without garment authoring files.</summary>
public static class AutomaticWardrobeFit
{
	const string ReferencePath = "models/characterkit/wardrobe-reference.bin";
	public const string Morph = "ce_body_female";
	// Cache by resource identity and authored cut, never by the current slider value.
	static readonly Dictionary<(Model, bool, bool), Model> fittedModels = new();
	static readonly Dictionary<(Model, bool, bool), string> fittingFailures = new();
	static Surface originalMale, unifiedMale, unifiedFemale;
	static Surface headSurface;
	static Model replacementNeck;
	public static int CachedModels => fittedModels.Count;

	public static void ClearCache()
	{
		fittedModels.Clear();
		fittingFailures.Clear();
		replacementNeck = null;
		headSurface = null;
	}

	/// <summary>Keep the lower neck when clothing replaces the head bodygroup.</summary>
	public static Model ReplacementNeck()
	{
		if (replacementNeck is not null) return replacementNeck;
		References();
		return replacementNeck = Build(Model.Load(UnifiedHuman.ModelPath), false, neckOnly: true);
	}

	public static bool IsFitted(Model model) => model?.Name.StartsWith("automatic-human-fit:", StringComparison.Ordinal) == true;
	public static Model Fit(Model source, bool female, bool fitFace = true)
	{
		if (source is null || source.IsError || IsFitted(source) || source.BoneCount == 0 || source.Bounds.Mins.z > 60)
			return source;
		if (fittedModels.TryGetValue((source, female, fitFace), out var ready))
			return ready;
		if (fittingFailures.TryGetValue((source, female, fitFace), out var problem))
			throw new InvalidOperationException(problem);
		try
		{
			References();
			var fitted = Build(source, female, fitFace: fitFace);
			fittedModels[(source, female, fitFace)] = fitted;
			return fitted;
		}
		catch (Exception error)
		{
			var message = "Cannot automatically fit " + source.Name + ": " + error.Message;
			fittingFailures[(source, female, fitFace)] = message;
			throw new InvalidOperationException(message, error);
		}
	}

	static void References()
	{
		if (unifiedMale is not null)
			return;
		using var stream = new MemoryStream(FileSystem.Mounted.ReadAllBytes(ReferencePath).ToArray());
		using var r = new BinaryReader(stream);
		if (r.ReadInt32() != 1)
			throw new InvalidOperationException("Unsupported human fitting reference.");
		originalMale = ReadSurface(r);
		unifiedMale = ReadSurface(r);
		unifiedFemale = new(unifiedMale.Female, unifiedMale.Female, unifiedMale.Triangles, unifiedMale.Names, unifiedMale.Weights);
	}

	static Surface HeadSurface()
	{
		if (headSurface is not null) return headSurface;
		var source = Model.Load(UnifiedHuman.ModelPath);
		var buffers = ReadBuffers(source);
		var points = new List<Vector3>();
		var triangles = new List<int>();
		for (int mi = 0; mi < buffers.Count; mi++)
		{
			var info = source.MeshInfo.Meshes[mi];
			if (!info.Name.StartsWith("head", StringComparison.OrdinalIgnoreCase)) continue;
			var buffer = buffers[mi];
			int callIndex = 0;
			foreach (var call in info.DrawCalls)
			{
				int vb = buffer.Vertices.Length == 1 ? 0 : callIndex;
				int ib = buffer.Indices.Length == 1 ? 0 : callIndex;
				var vertices = Vertices(buffer.Vertices[vb], null, buffer.Bones.Length);
				var indices = Indices(buffer.Indices[ib]);
				int offset = points.Count;
				points.AddRange(vertices.Select(v => v.Position));
				triangles.AddRange(indices.Skip(call.StartIndex).Take(call.Indices).Select(i => i + call.BaseVertex + offset));
				callIndex++;
			}
		}
		var p = points.ToArray();
		return headSurface = new(p, p, triangles.ToArray(), Array.Empty<string>(), Array.Empty<(int, float)[]>());
	}

	static Vector3 V(BinaryReader r) => new(r.ReadSingle(), r.ReadSingle(), r.ReadSingle());
	static Surface ReadSurface(BinaryReader r)
	{
		var names = new string[r.ReadInt32()];
		for (int i = 0; i < names.Length; i++)
			names[i] = Encoding.UTF8.GetString(r.ReadBytes(r.ReadInt32()));
		var count = r.ReadInt32();
		var triangles = new int[r.ReadInt32() * 3];
		var p = new Vector3[count];
		var f = new Vector3[count];
		var weights = new (int, float)[count][];
		for (int i = 0; i < count; i++)
		{
			p[i] = V(r);
			f[i] = V(r);
			weights[i] = new (int, float)[r.ReadInt32()];
			for (int j = 0; j < weights[i].Length; j++)
				weights[i][j] = (r.ReadInt32(), r.ReadSingle());
		}

		for (int i = 0; i < triangles.Length; i++)
			triangles[i] = r.ReadInt32();
		return new(p, f, triangles, names, weights);
	}

	[StructLayout(LayoutKind.Sequential)]
	struct SkinVertex
	{
		[VertexLayout.Position]
		public Vector3 Position;
		[VertexLayout.Normal]
		public Vector3 Normal;
		[VertexLayout.Tangent]
		public Vector4 Tangent;
		[VertexLayout.TexCoord(0)]
		public Vector4 UV;
		[VertexLayout.TexCoord(1)]
		public Vector4 UV2;
		[VertexLayout.Color]
		public Color32 Color;
		[VertexLayout.BlendIndices]
		public Color32 Indices;
		[VertexLayout.BlendWeight]
		public Color32 Weights;
	}

	readonly record struct Attribute(string Name, int Index, int Format, int Offset);
	sealed record Buffer(int Count, int Stride, Attribute[] Attributes, byte[] Data);
	sealed record Buffers(Buffer[] Vertices, Buffer[] Indices, string[] Bones);
	// SaveToVmdl emits unpacked MBUF data, even when the installed resource uses
	// compressed buffers. Reading that public serialization retains authored weights.
	static List<Buffers> ReadBuffers(Model source)
	{
		var data = source.SaveToVmdl();
		using var r = new BinaryReader(new MemoryStream(data));
		r.BaseStream.Position = 8;
		var table = 8 + r.ReadInt32();
		var blocks = r.ReadInt32();
		var result = new List<Buffers>();
		for (int i = 0; i < blocks; i++)
		{
			r.BaseStream.Position = table + i * 12;
			var tag = Encoding.ASCII.GetString(r.ReadBytes(4));
			var anchor = (int)r.BaseStream.Position;
			var offset = r.ReadInt32();
			var size = r.ReadInt32();
			if (tag == "MDAT")
			{
				r.BaseStream.Position = anchor + offset;
				var bones = SerializedMeshBones.Read(r.ReadBytes(size));
				result[^1] = result[^1] with
				{
					Bones = bones
				};
				continue;
			}

			if (tag != "MBUF")
				continue;
			var start = anchor + offset;
			r.BaseStream.Position = start;
			var vo = r.ReadInt32();
			var vc = r.ReadInt32();
			var io = r.ReadInt32();
			var ic = r.ReadInt32();
			result.Add(new(Enumerable.Range(0, vc).Select(j => ReadBuffer(r, start + vo + j * 24)).ToArray(), Enumerable.Range(0, ic).Select(j => ReadBuffer(r, start + 8 + io + j * 24)).ToArray(), null));
		}

		if (result.Count != source.MeshCount)
			throw new InvalidOperationException($"{source.Name}: missing serialized mesh buffers.");
		return result;
	}

	static Buffer ReadBuffer(BinaryReader r, int position)
	{
		r.BaseStream.Position = position;
		var count = r.ReadInt32();
		var stride = r.ReadInt32() & 0x3ffffff;
		var ap = (int)r.BaseStream.Position;
		var ao = r.ReadInt32();
		var ac = r.ReadInt32();
		var dp = (int)r.BaseStream.Position;
		var offset = r.ReadInt32();
		var bytes = r.ReadInt32();
		if (count < 0 || count > 2000000 || stride < 1 || stride > 1024 || bytes != count * stride)
			throw new InvalidOperationException("Unexpected serialized human clothing buffer.");
		var attributes = new Attribute[ac];
		for (int i = 0; i < ac; i++)
		{
			r.BaseStream.Position = ap + ao + i * 56;
			var name = Encoding.ASCII.GetString(r.ReadBytes(32)).TrimEnd('\0');
			attributes[i] = new(name, r.ReadInt32(), r.ReadInt32(), r.ReadInt32());
		}

		r.BaseStream.Position = dp + offset;
		return new(count, stride, attributes, r.ReadBytes(bytes));
	}

	static float F(byte[] b, int p) => BitConverter.ToSingle(b, p);
	static ushort U16(byte[] b, int p) => BitConverter.ToUInt16(b, p);
	static uint U32(byte[] b, int p) => BitConverter.ToUInt32(b, p);
	static Color32 Bytes(byte[] b, int p) => new(b[p], b[p + 1], b[p + 2], b[p + 3]);
	static Vector4 Values(byte[] b, int p, int format)
	{
		float HalfAt(int i) => (float)BitConverter.UInt16BitsToHalf(U16(b, p + i * 2));
		return format switch
		{
			2 => new(F(b, p), F(b, p + 4), F(b, p + 8), F(b, p + 12)),
			6 => new(F(b, p), F(b, p + 4), F(b, p + 8), 0),
			16 => new(F(b, p), F(b, p + 4), 0, 0),
			10 => new(HalfAt(0), HalfAt(1), HalfAt(2), HalfAt(3)),
			34 => new(HalfAt(0), HalfAt(1), 0, 0),
			28 => new(b[p] / 255f, b[p + 1] / 255f, b[p + 2] / 255f, b[p + 3] / 255f),
			_ => throw new InvalidOperationException("Unsupported clothing vertex format " + format)};
	}

	static SkinVertex[] Vertices(Buffer b, Vertex[] native, int boneCount)
	{
		var vertices = new SkinVertex[b.Count];
		for (int i = 0; i < b.Count; i++)
		{
			var v = new SkinVertex
			{
				Color = new(255, 255, 255, 255),
				Weights = new(255, 0, 0, 0)
			};
			foreach (var a in b.Attributes)
			{
				var p = i * b.Stride + a.Offset;
				switch (a.Name)
				{
					case "POSITION":
						var xyz = Values(b.Data, p, a.Format);
						v.Position = new(xyz.x, xyz.y, xyz.z);
						break;
					case "TEXCOORD":
						if (a.Index == 0)
							v.UV = Values(b.Data, p, a.Format);
						else if (a.Index == 1)
							v.UV2 = Values(b.Data, p, a.Format);
						break;
					case "BLENDINDICES":
						if (a.Format == 30)
							v.Indices = Bytes(b.Data, p);
						else if (a.Format == 12)
						{
							var ids = Enumerable.Range(0, 4).Select(j => (int)U16(b.Data, p + j * 2)).ToArray();
							if (ids.Any(j => j >= boneCount || j > 255))
								throw new InvalidOperationException("Clothing bone index exceeds fitting layout.");
							v.Indices = new((byte)ids[0], (byte)ids[1], (byte)ids[2], (byte)ids[3]);
						}
						else
							throw new InvalidOperationException("Unsupported clothing bone-index format " + a.Format);
						break;
					case "BLENDWEIGHT":
						if (a.Format != 28)
							throw new InvalidOperationException("Unsupported clothing weight format " + a.Format);
						v.Weights = Bytes(b.Data, p);
						break;
					case "COLOR":
						if (a.Format == 28)
							v.Color = Bytes(b.Data, p);
						break;
					case "NORMAL":
						if (a.Format == 6)
						{
							var n = Values(b.Data, p, a.Format);
							v.Normal = new(n.x, n.y, n.z);
						}

						break;
					case "TANGENT":
						v.Tangent = Values(b.Data, p, a.Format);
						break;
				}
			}

			if (v.Normal.LengthSquared < .01f && native is not null && native.Length == b.Count && (native[i].Position - v.Position).Length < .001f)
			{
				v.Normal = native[i].Normal;
				v.Tangent = native[i].Tangent;
			}

			if (new[]
			{
				v.Indices.r,
				v.Indices.g,
				v.Indices.b,
				v.Indices.a
			}.Any(j => j >= boneCount))
				throw new InvalidOperationException("Serialized garment bone map does not match its skeleton.");
			vertices[i] = v;
		}

		return vertices;
	}

	static int[] Indices(Buffer b) => Enumerable.Range(0, b.Count).Select(i => b.Stride == 2 ? (int)U16(b.Data, i * 2) : checked((int)U32(b.Data, i * 4))).ToArray();
	static Vector3[] Normals(Vector3[] p, int[] indices)
	{
		var n = new Vector3[p.Length];
		for (int i = 0; i + 2 < indices.Length; i += 3)
		{
			var a = indices[i];
			var b = indices[i + 1];
			var c = indices[i + 2];
			var cross = Vector3.Cross(p[b] - p[a], p[c] - p[a]);
			n[a] += cross;
			n[b] += cross;
			n[c] += cross;
		}

		return n.Select(v => v.Normal).ToArray();
	}

	static Vector4[] Tangents(SkinVertex[] vertices, Vector3[] p, Vector3[] normals, int[] indices)
	{
		var t = new Vector3[p.Length];
		var b = new Vector3[p.Length];
		for (int i = 0; i + 2 < indices.Length; i += 3)
		{
			int a = indices[i], c = indices[i + 1], d = indices[i + 2];
			var e = p[c] - p[a];
			var f = p[d] - p[a];
			var uv = vertices[c].UV - vertices[a].UV;
			var uw = vertices[d].UV - vertices[a].UV;
			float determinant = uv.x * uw.y - uv.y * uw.x;
			if (MathF.Abs(determinant) < 1e-8f)
				continue;
			var tangent = (e * uw.y - f * uv.y) / determinant;
			var bitangent = (f * uv.x - e * uw.x) / determinant;
			foreach (int id in new[]
			{
				a,
				c,
				d
			}

			)
			{
				t[id] += tangent;
				b[id] += bitangent;
			}
		}

		return Enumerable.Range(0, p.Length).Select(i =>
		{
			var n = normals[i];
			var tangent = (t[i] - n * Vector3.Dot(n, t[i])).Normal;
			return new Vector4(tangent, Vector3.Dot(Vector3.Cross(n, tangent), b[i]) < 0 ? -1 : 1);
		}).ToArray();
	}

	static Model Build(Model source, bool female, bool neckOnly = false, bool fitFace = true)
	{
		var buffers = ReadBuffers(source);
		var builder = Model.Builder.WithName("automatic-human-fit:" + source.Name + ":" + (neckOnly ? "replacement-neck" : female ? "female" : "male"));
		for (int i = 0; i < source.BoneCount; i++)
		{
			var transform = source.GetBoneTransform(i);
			var parent = source.GetBoneParent(i);
			builder.AddBone(source.GetBoneName(i), transform.Position, transform.Rotation, parent < 0 ? null : source.GetBoneName(parent));
		}

		var native = source.GetVertices();
		var baseMaterials = (source.MaterialGroupCount > 0 ? source.GetMaterials(0) : source.Materials).ToArray();
		var meshMaterials = new List<Material>();
		var work = new List<Action>();
		var clothMale = new List<Vector3>();
		var clothFemale = new List<Vector3>();
		var clothIndices = new List<int>();
		var head = neckOnly || !fitFace ? null : HeadSurface();
		for (int mi = 0; mi < buffers.Count; mi++)
		{
			var info = source.MeshInfo.Meshes[mi];
			if (neckOnly && !info.Name.StartsWith("head", StringComparison.OrdinalIgnoreCase)) continue;
			var buffer = buffers[mi];
			var boneMap = buffer.Bones.Select(name => Enumerable.Range(0, source.BoneCount).First(i => source.GetBoneName(i) == name)).ToArray();
			int callIndex = 0;
			foreach (var call in info.DrawCalls)
			{
				// SaveToVmdl emits individual buffers in draw-call order when a mesh uses
				// separate streams (its metadata currently labels each stream as zero).
				int vb = buffer.Vertices.Length == 1 ? 0 : callIndex;
				int ib = buffer.Indices.Length == 1 ? 0 : callIndex;
				if (vb >= buffer.Vertices.Length || ib >= buffer.Indices.Length)
					throw new InvalidOperationException("Unexpected garment draw-call buffer mapping.");
				var verts = Vertices(buffer.Vertices[vb], mi == 0 && buffer.Vertices.Length == 1 ? native : null, buffer.Bones.Length);
				var indices = Indices(buffer.Indices[ib]);
				for (int i = 0; i < verts.Length; i++)
				{
					var ids = verts[i].Indices;
					verts[i].Indices = new((byte)boneMap[ids.r], (byte)boneMap[ids.g], (byte)boneMap[ids.b], (byte)boneMap[ids.a]);
				}

				callIndex++;
				var selectedIndices = indices.Skip(call.StartIndex).Take(call.Indices).Select(i => i + call.BaseVertex).ToArray();
				if (neckOnly)
					selectedIndices = Enumerable.Range(0, selectedIndices.Length / 3)
						.Where(t => Enumerable.Range(0, 3).All(j => verts[selectedIndices[t * 3 + j]].Position.z <= 60))
						.SelectMany(t => selectedIndices.Skip(t * 3).Take(3)).ToArray();
				if (selectedIndices.Length == 0) continue;
				var used = selectedIndices.Distinct().ToArray();
				var lookup = used.Select((id, i) => (id, i)).ToDictionary(x => x.id, x => x.i);
				var localIndices = selectedIndices.Select(i => lookup[i]).ToArray();
				var v = used.Select(i => verts[i]).ToArray();
				var original = v.Select(x => x.Position).ToArray();
				var male = new Vector3[v.Length];
				var target = new Vector3[v.Length];
				bool skin = neckOnly || call.Material?.Name.Contains("/bodies/", StringComparison.OrdinalIgnoreCase) == true;
				for (int i = 0; i < v.Length; i++)
				{
					var p = original[i];
					if (!skin && !neckOnly && fitFace && p.z >= 57 && p.z <= 73)
					{
						var hit = head.Nearest(p);
						if (hit.Distance < 1.2f && hit.Distance < originalMale.Nearest(p).Distance)
						{
							var normal = head.FaceNormal(hit.Triangle);
							float clearance = Vector3.Dot(p - hit.Point, normal);
							// Close face coverings need space for the difference in falloff
							// between the skin and fabric under inherited facial deformers.
							float margin = .60f * Math.Clamp((62.2f - p.z) / 1.0f, 0, 1);
							if (margin > 0) p += normal * Math.Clamp(margin - clearance, 0, .85f);
							male[i] = target[i] = p;
							continue;
						}
					}
					if (p.z > 60)
					{
						male[i] = p;
						target[i] = p;
						continue;
					}

					var from = (neckOnly ? unifiedMale : female ? unifiedFemale : originalMale).Nearest(p);
					var mapped = female || neckOnly ? from : unifiedMale.Nearest(from.Point);
					var m = unifiedMale.Interpolate(mapped.Triangle, mapped.Barycentric);
					var f = unifiedFemale.Interpolate(mapped.Triangle, mapped.Barycentric);
					// Preserve authored clearance, including compressed skin under clothing.
					var offset = p - from.Point;
					bool snap = skin && offset.Length < .12f;
					male[i] = m + (snap ? Vector3.Zero : offset);
					target[i] = f + (snap ? Vector3.Zero : offset);
				}

				int clothOffset = clothMale.Count;
				if (!skin)
				{
					clothIndices.AddRange(localIndices.Select(i => i + clothOffset));
					clothMale.AddRange(male);
					clothFemale.AddRange(target);
				}

				work.Add(() =>
				{
					if (!skin)
					{
						for (int i = 0; i < v.Length; i++)
						{
							male[i] = clothMale[clothOffset + i];
							target[i] = clothFemale[clothOffset + i];
						}
					}

					var originalNormals = Normals(original, localIndices);
					var maleNormals = Normals(male, localIndices);
					var femaleNormals = Normals(target, localIndices);
					var tangents = Tangents(v, male, maleNormals, localIndices);
					var morph = new MorphDelta[v.Length];
					for (int i = 0; i < v.Length; i++)
					{
						if (v[i].Normal.LengthSquared < .01f)
							v[i].Normal = originalNormals[i];
						v[i].Normal = (v[i].Normal + maleNormals[i] - originalNormals[i]).Normal;
						v[i].Tangent = tangents[i];
						v[i].Position = male[i];
						morph[i] = new(i, target[i] - male[i], femaleNormals[i] - maleNormals[i]);
					}

					var mesh = new Mesh(call.Material, (MeshPrimitiveType)call.PrimitiveType);
					mesh.UvDensity = call.UvDensity;
					mesh.CreateVertexBuffer(v.Length, v.ToList());
					mesh.CreateIndexBuffer(localIndices.Length, localIndices.ToList());
					mesh.AddMorph(Morph, morph);
					mesh.Bounds = new BBox(male.Concat(target).Aggregate((a, b) => Vector3.Min(a, b)), male.Concat(target).Aggregate((a, b) => Vector3.Max(a, b)));
					var lods = info.LodMask;
					if (source.MeshInfo.LodCount == 0)
						builder.AddMesh(mesh);
					else
						for (int lod = 0; lod < source.MeshInfo.LodCount; lod++)
							if ((lods & (1 << lod)) != 0)
								builder.AddMesh(mesh, lod);
					meshMaterials.Add(call.Material);
				});
			}
		}

		if (clothMale.Count > 0)
		{
			var m = clothMale.ToArray();
			var f = clothFemale.ToArray();
			var t = clothIndices.ToArray();
			ClearTriangles(unifiedMale, m, t);
			ClearTriangles(unifiedFemale, f, t);
			clothMale = m.ToList();
			clothFemale = f.ToList();
		}

		foreach (var build in work)
			build();
		for (int group = 0; group < source.MaterialGroupCount; group++)
		{
			var materials = source.GetMaterials(group).ToArray();
			var g = builder.AddMaterialGroup(source.GetMaterialGroupName(group));
			foreach (var material in meshMaterials)
			{
				var index = Array.IndexOf(baseMaterials, material);
				g.AddMaterial(index >= 0 && index < materials.Length ? materials[index] : material);
			}
		}

		for (int lod = 1; lod < source.MeshInfo.LodSwitchDistances.Length; lod++)
			builder.WithLodDistance(lod, source.MeshInfo.LodSwitchDistances[lod]);
		return builder.Create();
	}

	static void ClearTriangles(Surface body, Vector3[] points, int[] triangles)
	{
		// A smooth spatial clearance field moves folded/layered fabric together.
		// Local per-triangle pushes can pull an inner layer through its outer flap.
		var source = points.ToArray();
		var min = points.Aggregate((a, b) => Vector3.Min(a, b)) - Vector3.One * 1.2f;
		var max = points.Aggregate((a, b) => Vector3.Max(a, b)) + Vector3.One * 1.2f;
		var candidates = body.Points.Where(p => p.x >= min.x && p.x <= max.x && p.y >= min.y && p.y <= max.y && p.z >= min.z && p.z <= max.z).ToArray();
		var pushes = new List<(Vector3 Point, Vector3 Push)>();
		var fabric = new Surface(points, points, triangles, Array.Empty<string>(), Array.Empty<(int, float)[]>());
		foreach (var p in candidates)
		{
			var hit = fabric.Nearest(p);
			if (hit.Distance > 1.2f)
				continue;
			var w = hit.Barycentric;
			if (MathF.Min(w.x, MathF.Min(w.y, w.z)) < .015f)
				continue;
			var n = body.FaceNormal(body.Nearest(p).Triangle);
			if (MathF.Abs(Vector3.Dot(n, fabric.FaceNormal(hit.Triangle))) < .25f)
				continue;
			var penetration = Vector3.Dot(p - hit.Point, n) + .14f;
			if (penetration <= 0 || penetration > 1.2f)
				continue;
			pushes.Add((hit.Point, n * MathF.Min(.7f, penetration + .04f)));
		}

		// Spatial buckets bound the search and apply identical displacement to UV
		// duplicates, overlapping material sections, and neighbouring LOD vertices.
		const float radius = 2.5f;
		(int, int, int) Cell(Vector3 p) => ((int)MathF.Floor(p.x / radius), (int)MathF.Floor(p.y / radius), (int)MathF.Floor(p.z / radius));
		var buckets = pushes.GroupBy(p => Cell(p.Point)).ToDictionary(g => g.Key, g => g.ToArray());
		for (int i = 0; i < points.Length; i++)
		{
			var cell = Cell(source[i]);
			Vector3 best = default;
			for (int x = -1; x <= 1; x++)
				for (int y = -1; y <= 1; y++)
					for (int z = -1; z <= 1; z++)
					{
						if (!buckets.TryGetValue((cell.Item1 + x, cell.Item2 + y, cell.Item3 + z), out var values))
							continue;
						foreach (var p in values)
						{
							var d = (source[i] - p.Point).Length / radius;
							if (d >= 1)
								continue;
							var falloff = 1 - d * d * (3 - 2 * d);
							var push = p.Push * falloff;
							if (push.LengthSquared > best.LengthSquared)
								best = push;
						}
					}

			points[i] += best;
		}
	}

	readonly record struct Hit(int Triangle, Vector3 Barycentric, Vector3 Point, float Distance);
	sealed class Surface
	{
		public readonly Vector3[] Points, Female;
		public readonly int[] Triangles;
		public readonly string[] Names;
		public readonly (int, float)[][] Weights;
		readonly Node tree;
		sealed class Node
		{
			public Vector3 Min, Max;
			public Node Left, Right;
			public int[] Faces;
		}

		public Surface(Vector3[] p, Vector3[] f, int[] t, string[] names, (int, float)[][] weights)
		{
			Points = p;
			Female = f;
			Triangles = t;
			Names = names;
			Weights = weights;
			tree = BuildTree(Enumerable.Range(0, t.Length / 3).ToArray());
		}

		Vector3 A(int t, int corner) => Points[Triangles[t * 3 + corner]];
		Node BuildTree(int[] faces)
		{
			var p = faces.SelectMany(i => new[] { A(i, 0), A(i, 1), A(i, 2) }).ToArray();
			var node = new Node
			{
				Min = p.Aggregate((a, b) => Vector3.Min(a, b)),
				Max = p.Aggregate((a, b) => Vector3.Max(a, b))
			};
			if (faces.Length <= 12)
			{
				node.Faces = faces;
				return node;
			}

			var extent = node.Max - node.Min;
			int axis = extent.x > extent.y ? (extent.x > extent.z ? 0 : 2) : (extent.y > extent.z ? 1 : 2);
			float Coord(Vector3 v) => axis == 0 ? v.x : axis == 1 ? v.y : v.z;
			var sorted = faces.OrderBy(i => Coord(A(i, 0) + A(i, 1) + A(i, 2))).ToArray();
			int half = sorted.Length / 2;
			node.Left = BuildTree(sorted.Take(half).ToArray());
			node.Right = BuildTree(sorted.Skip(half).ToArray());
			return node;
		}

		public Vector3 Interpolate(int t, Vector3 w) => A(t, 0) * w.x + A(t, 1) * w.y + A(t, 2) * w.z;
		public Vector3 FaceNormal(int t) => Vector3.Cross(A(t, 1) - A(t, 0), A(t, 2) - A(t, 0)).Normal;
		static float Distance(Node n, Vector3 p) => (p - Vector3.Max(n.Min, Vector3.Min(n.Max, p))).LengthSquared;
		public Hit Nearest(Vector3 p)
		{
			float best = float.MaxValue;
			int face = 0;
			Vector3 bary = default, point = default;
			void Search(Node node)
			{
				if (Distance(node, p) > best)
					return;
				if (node.Faces is null)
				{
					var a = node.Left;
					var b = node.Right;
					if (Distance(a, p) > Distance(b, p))
						(a, b) = (b, a);
					Search(a);
					Search(b);
					return;
				}

				foreach (var t in node.Faces)
				{
					var w = Closest(p, A(t, 0), A(t, 1), A(t, 2));
					var q = Interpolate(t, w);
					var distance = (p - q).LengthSquared;
					if (distance < best)
					{
						best = distance;
						face = t;
						bary = w;
						point = q;
					}
				}
			}

			Search(tree);
			return new(face, bary, point, MathF.Sqrt(best));
		}

		static Vector3 Closest(Vector3 p, Vector3 a, Vector3 b, Vector3 c)
		{
			var ab = b - a;
			var ac = c - a;
			var ap = p - a;
			float d1 = Vector3.Dot(ab, ap), d2 = Vector3.Dot(ac, ap);
			if (d1 <= 0 && d2 <= 0)
				return new(1, 0, 0);
			var bp = p - b;
			float d3 = Vector3.Dot(ab, bp), d4 = Vector3.Dot(ac, bp);
			if (d3 >= 0 && d4 <= d3)
				return new(0, 1, 0);
			float vc = d1 * d4 - d3 * d2;
			if (vc <= 0 && d1 >= 0 && d3 <= 0)
			{
				var v = d1 / (d1 - d3);
				return new(1 - v, v, 0);
			}

			var cp = p - c;
			float d5 = Vector3.Dot(ab, cp), d6 = Vector3.Dot(ac, cp);
			if (d6 >= 0 && d5 <= d6)
				return new(0, 0, 1);
			float vb = d5 * d2 - d1 * d6;
			if (vb <= 0 && d2 >= 0 && d6 <= 0)
			{
				var w = d2 / (d2 - d6);
				return new(1 - w, 0, w);
			}

			float va = d3 * d6 - d5 * d4;
			if (va <= 0 && d4 - d3 >= 0 && d5 - d6 >= 0)
			{
				var w = (d4 - d3) / (d4 - d3 + d5 - d6);
				return new(0, 1 - w, w);
			}

			float sum = va + vb + vc;
			if (MathF.Abs(sum) < 1e-12f)
				return new(1, 0, 0);
			return new(va / sum, vb / sum, vc / sum);
		}
	}
}// Runtime fitting integration.