Sim/SkinLayer.cs
using System;
using System.Collections.Generic;

namespace SquishyTerry.Sim;

/// <summary>
/// The rendered skin: a finer copy of the face mesh that rides on smooth PN-triangle patches of the
/// simulation mesh and is allowed to buckle into wrinkles. Port of the original game's GPU cloth
/// pass (JS classes I and g plus attachment/distance/normals shaders), run on the CPU.
///
/// Each skin vertex may drift up to 0.008 * wrinkleStrength from its PN attachment point while
/// neighbour and across-edge distance constraints try to keep its original spacing, so compressed
/// skin has nowhere to go but out: wrinkles.
/// </summary>
public sealed class SkinLayer
{
	const float AttachRadius = 0.008f;
	const int Iterations = 15;

	public Vec3[] Pos;
	public Vec3[] Normals;
	/// <summary>Triangles to draw (hidden base triangles removed).</summary>
	public int[] Indices;
	public float[] Mouthiness;
	public float[] Strength;

	readonly int[] _baseTri;
	readonly float[] _baryW, _baryU, _baryV;
	readonly Vec3[] _attach;
	Vec3[] _temp;
	readonly float[] _px, _py, _pz;

	// Up to 6 neighbours / opposites per vertex, -1 padded.
	readonly int[] _neighbours;
	readonly int[] _opposites;
	readonly float[] _neighbourRest;
	readonly float[] _oppositeRest;
	readonly bool[] _fixed;
	readonly List<int>[] _baseTriangleVertices;
	bool _initialized;

	/// <summary>
	/// Run the wrinkle cloth. When false the skin just follows the smooth PN surface, like the
	/// original's mobile "static skin" path.
	/// </summary>
	public bool Wrinkles = true;

	Vec3[] _prevAttach;
	int _quietFrames;

	public int VertexCount => Pos.Length;

	// Read-only views of the static data, for running the same solve on the GPU.
	public int[] BaseTriangle => _baseTri;
	public float[] BaryW => _baryW;
	public float[] BaryU => _baryU;
	public float[] BaryV => _baryV;
	/// <summary>Ordered one-ring, 6 per vertex, -1 padded.</summary>
	public int[] Neighbours => _neighbours;
	/// <summary>Across-edge vertices, 6 per vertex, -1 padded.</summary>
	public int[] Opposites => _opposites;
	public float[] NeighbourRest => _neighbourRest;
	public float[] OppositeRest => _oppositeRest;
	/// <summary>Skin vertices riding on a base triangle (pinned while that triangle is grabbed).</summary>
	public List<int> VerticesOfBaseTriangle( int triangle ) => _baseTriangleVertices[triangle];

	/// <summary>Compute attachment points and rest lengths without running the cloth.</summary>
	public void InitializeRest( Face face )
	{
		ComputeAttachments( face );
		for ( int v = 0; v < Pos.Length; v++ ) Pos[v] = _attach[v];
		RestLengths( _neighbours, _neighbourRest );
		RestLengths( _opposites, _oppositeRest );
		_prevAttach = Geometry.Copy( _attach );
		_initialized = true;
		ComputeNormals();
	}

	struct Association
	{
		public int Tri;
		public float W, U, V;
	}

	sealed class Vertex
	{
		public readonly List<Association> Associations = new();
		public readonly List<int> Triangles = new();
		public bool OnBorder, BeyondBorder;
	}

	struct Tri
	{
		public int A, B, C;
		public bool Hidden, BeyondBorder;
	}

	public SkinLayer( Face face, int subdivisions )
	{
		int baseTriCount = face.TriangleCount;
		var hiddenBase = new bool[baseTriCount];
		foreach ( int t in face.HiddenTriangles ) hiddenBase[t] = true;

		var verts = new List<Vertex>();
		var tris = new List<Tri>();
		for ( int v = 0; v < face.VertexCount; v++ ) verts.Add( new Vertex() );
		for ( int t = 0; t < baseTriCount; t++ )
		{
			int a = face.Indices[t * 3], b = face.Indices[t * 3 + 1], c = face.Indices[t * 3 + 2];
			tris.Add( new Tri { A = a, B = b, C = c, Hidden = hiddenBase[t] } );
			verts[a].Associations.Add( new Association { Tri = t, W = 1 } );
			verts[b].Associations.Add( new Association { Tri = t, U = 1 } );
			verts[c].Associations.Add( new Association { Tri = t, V = 1 } );
		}

		// --- Midpoint subdivision, tracking barycentric coordinates in the base triangles ---
		var midpoints = new Dictionary<long, int>();
		int Mid( int i, int j )
		{
			long key = i < j ? ((long)i << 32) | (uint)j : ((long)j << 32) | (uint)i;
			if ( midpoints.TryGetValue( key, out int m ) ) return m;
			var nv = new Vertex();
			foreach ( var ai in verts[i].Associations )
			{
				foreach ( var aj in verts[j].Associations )
				{
					if ( ai.Tri != aj.Tri ) continue;
					nv.Associations.Add( new Association
					{
						Tri = ai.Tri,
						W = (ai.W + aj.W) * 0.5f,
						U = (ai.U + aj.U) * 0.5f,
						V = (ai.V + aj.V) * 0.5f,
					} );
				}
			}
			verts.Add( nv );
			midpoints[key] = verts.Count - 1;
			return verts.Count - 1;
		}

		for ( int s = 0; s < subdivisions; s++ )
		{
			var next = new List<Tri>( tris.Count * 4 );
			foreach ( var t in tris )
			{
				int g = Mid( t.A, t.B ), h = Mid( t.B, t.C ), k = Mid( t.C, t.A );
				next.Add( new Tri { A = t.A, B = g, C = k, Hidden = t.Hidden } );
				next.Add( new Tri { A = t.B, B = h, C = g, Hidden = t.Hidden } );
				next.Add( new Tri { A = t.C, B = k, C = h, Hidden = t.Hidden } );
				next.Add( new Tri { A = g, B = h, C = k, Hidden = t.Hidden } );
			}
			tris = next;
		}

		// --- Keep visible vertices plus a one-triangle band past the hidden border ---
		var keep = new bool[verts.Count];
		foreach ( var t in tris )
			if ( !t.Hidden ) { keep[t.A] = true; keep[t.B] = true; keep[t.C] = true; }
		for ( int i = 0; i < tris.Count; i++ )
		{
			var t = tris[i];
			if ( !t.Hidden || !(keep[t.A] || keep[t.B] || keep[t.C]) ) continue;
			foreach ( int v in new[] { t.A, t.B, t.C } )
			{
				if ( keep[v] ) verts[v].OnBorder = true; else verts[v].BeyondBorder = true;
				keep[v] = true;
			}
			t.BeyondBorder = true;
			tris[i] = t;
		}

		var remap = new int[verts.Count];
		var kept = new List<Vertex>();
		for ( int i = 0; i < verts.Count; i++ )
		{
			remap[i] = keep[i] ? kept.Count : -1;
			if ( keep[i] ) kept.Add( verts[i] );
		}
		var allTris = new List<Tri>();
		foreach ( var t in tris )
		{
			if ( t.Hidden && !t.BeyondBorder ) continue;
			allTris.Add( new Tri { A = remap[t.A], B = remap[t.B], C = remap[t.C], Hidden = t.Hidden, BeyondBorder = t.BeyondBorder } );
		}
		for ( int i = 0; i < allTris.Count; i++ )
		{
			kept[allTris[i].A].Triangles.Add( i );
			kept[allTris[i].B].Triangles.Add( i );
			kept[allTris[i].C].Triangles.Add( i );
		}

		int n = kept.Count;
		_neighbours = new int[n * 6];
		_opposites = new int[n * 6];
		Array.Fill( _neighbours, -1 );
		Array.Fill( _opposites, -1 );

		// Ordered one-ring (sortTriangles): walk the fan across each (centre, next) edge.
		for ( int v = 0; v < n; v++ )
		{
			if ( kept[v].BeyondBorder )
			{
				for ( int k = 0; k < 6; k++ ) _neighbours[v * 6 + k] = k;
				continue;
			}
			var fan = kept[v].Triangles;
			int start = fan[0], cur = start, count = 0;
			do
			{
				var t = allTris[cur];
				int after = t.A == v ? t.B : t.B == v ? t.C : t.A;
				if ( count < 6 ) _neighbours[v * 6 + count] = after;
				count++;
				cur = -1;
				foreach ( int ti in fan )
				{
					var o = allTris[ti];
					if ( (o.A == after && o.B == v) || (o.B == after && o.C == v) || (o.C == after && o.A == v) ) { cur = ti; break; }
				}
			}
			while ( cur != -1 && cur != start && count < 12 );
		}

		// Opposites: the two vertices facing each other across every edge.
		var firstOpposite = new Dictionary<long, int>();
		var oppositeCount = new int[n];
		void AddOpposite( int v, int o )
		{
			if ( oppositeCount[v] < 6 ) _opposites[v * 6 + oppositeCount[v]] = o;
			oppositeCount[v]++;
		}
		void Edge( int a, int b, int c )
		{
			long key = a < b ? ((long)a << 32) | (uint)b : ((long)b << 32) | (uint)a;
			if ( firstOpposite.TryGetValue( key, out int other ) )
			{
				AddOpposite( c, other );
				AddOpposite( other, c );
			}
			else firstOpposite[key] = c;
		}
		foreach ( var t in allTris )
		{
			Edge( t.A, t.B, t.C );
			Edge( t.B, t.C, t.A );
			Edge( t.C, t.A, t.B );
		}
		for ( int v = 0; v < n; v++ )
		{
			if ( !kept[v].BeyondBorder ) continue;
			for ( int k = 0; k < 6; k++ ) _opposites[v * 6 + k] = k;
		}

		var draw = new List<int>();
		foreach ( var t in allTris )
		{
			if ( t.BeyondBorder ) continue;
			draw.Add( t.A ); draw.Add( t.B ); draw.Add( t.C );
		}
		Indices = draw.ToArray();

		_baseTri = new int[n];
		_baryW = new float[n];
		_baryU = new float[n];
		_baryV = new float[n];
		Strength = new float[n];
		Mouthiness = new float[n];
		_baseTriangleVertices = new List<int>[baseTriCount];
		for ( int t = 0; t < baseTriCount; t++ ) _baseTriangleVertices[t] = new List<int>();
		for ( int v = 0; v < n; v++ )
		{
			var a0 = kept[v].Associations[0];
			_baseTri[v] = a0.Tri;
			_baryW[v] = a0.W; _baryU[v] = a0.U; _baryV[v] = a0.V;
			int ia = face.Indices[a0.Tri * 3], ib = face.Indices[a0.Tri * 3 + 1], ic = face.Indices[a0.Tri * 3 + 2];
			Strength[v] = a0.W * face.Wrinkle[ia] + a0.U * face.Wrinkle[ib] + a0.V * face.Wrinkle[ic];
			Mouthiness[v] = a0.W * face.Mouthiness[ia] + a0.U * face.Mouthiness[ib] + a0.V * face.Mouthiness[ic];
			if ( kept[v].OnBorder || kept[v].BeyondBorder ) Strength[v] = 0;
			foreach ( var assoc in kept[v].Associations )
				_baseTriangleVertices[assoc.Tri].Add( v );
		}

		Pos = new Vec3[n];
		Normals = new Vec3[n];
		_attach = new Vec3[n];
		_temp = new Vec3[n];
		_px = new float[n];
		_py = new float[n];
		_pz = new float[n];
		_neighbourRest = new float[n * 6];
		_oppositeRest = new float[n * 6];
		_fixed = new bool[n];
	}

	/// <summary>Runs one frame of the wrinkle cloth against the current face (g.prototype.update).</summary>
	public void Update( Face face, bool dragging, int dragTriangle )
	{
		ComputeAttachments( face );

		if ( !_initialized )
		{
			for ( int v = 0; v < Pos.Length; v++ ) Pos[v] = _attach[v];
			RestLengths( _neighbours, _neighbourRest );
			RestLengths( _opposites, _oppositeRest );
			_prevAttach = Geometry.Copy( _attach );
			_initialized = true;
		}

		if ( !Wrinkles )
		{
			for ( int v = 0; v < Pos.Length; v++ ) Pos[v] = _attach[v];
			ComputeNormals();
			return;
		}

		// Once the face has been still for a while the cloth has converged; skip the solve.
		float maxMove = 0;
		for ( int v = 0; v < _attach.Length; v++ )
		{
			maxMove = MathF.Max( maxMove, (_attach[v] - _prevAttach[v]).LengthSquared );
			_prevAttach[v] = _attach[v];
		}
		_quietFrames = dragging || maxMove > 1e-10f ? 0 : _quietFrames + 1;
		if ( _quietFrames > 30 ) return;

		// The original only rewrites the fixed set while dragging, so the last grabbed patch stays
		// pinned (except on the final iteration) until the next grab.
		if ( dragging )
		{
			Array.Clear( _fixed );
			foreach ( int v in _baseTriangleVertices[dragTriangle] ) _fixed[v] = true;
		}

		for ( int it = 0; it < Iterations; it++ )
		{
			if ( it == 0 )
			{
				for ( int v = 0; v < Pos.Length; v++ )
				{
					float radius = AttachRadius * Strength[v];
					var d = _attach[v] - Pos[v];
					float dist = d.Length;
					if ( dist > radius ) Pos[v] += d * ((dist - radius) / dist);
				}
			}

			bool respectFixed = it < Iterations - 1;
			Relax( _neighbours, _neighbourRest, respectFixed );
			Relax( _opposites, _oppositeRest, respectFixed );
		}

		ComputeNormals();
	}

	void ComputeAttachments( Face face )
	{
		var p = face.Pos;
		var nrm = face.Normals;
		var idx = face.Indices;
		for ( int v = 0; v < _attach.Length; v++ )
		{
			int t = _baseTri[v];
			var p1 = p[idx[t * 3]]; var p2 = p[idx[t * 3 + 1]]; var p3 = p[idx[t * 3 + 2]];
			var n1 = nrm[idx[t * 3]]; var n2 = nrm[idx[t * 3 + 1]]; var n3 = nrm[idx[t * 3 + 2]];
			_attach[v] = PointNormalTriangle( p1, p2, p3, n1, n2, n3, _baryW[v], _baryU[v], _baryV[v] );
		}
	}

	/// <summary>Cubic PN-triangle position (Vlachos et al. 2001), as in attachmentpositions.frag.</summary>
	static Vec3 PointNormalTriangle( Vec3 p1, Vec3 p2, Vec3 p3, Vec3 n1, Vec3 n2, Vec3 n3, float w, float u, float v )
	{
		const float third = 1.0f / 3.0f;
		var b210 = (p1 * 2 + p2 - n1 * Vec3.Dot( p2 - p1, n1 )) * third;
		var b120 = (p2 * 2 + p1 - n2 * Vec3.Dot( p1 - p2, n2 )) * third;
		var b021 = (p2 * 2 + p3 - n2 * Vec3.Dot( p3 - p2, n2 )) * third;
		var b012 = (p3 * 2 + p2 - n3 * Vec3.Dot( p2 - p3, n3 )) * third;
		var b102 = (p3 * 2 + p1 - n3 * Vec3.Dot( p1 - p3, n3 )) * third;
		var b201 = (p1 * 2 + p3 - n1 * Vec3.Dot( p3 - p1, n1 )) * third;
		var e = (b210 + b120 + b021 + b012 + b102 + b201) / 6.0f;
		var vc = (p1 + p2 + p3) * third;
		var b111 = e + (e - vc) * 0.5f;

		float u2 = u * u, v2 = v * v, w2 = w * w;
		return p1 * (w2 * w) + p2 * (u2 * u) + p3 * (v2 * v)
			+ b210 * (3 * w2 * u) + b120 * (3 * w * u2) + b201 * (3 * w2 * v)
			+ b021 * (3 * u2 * v) + b102 * (3 * w * v2) + b012 * (3 * u * v2)
			+ b111 * (6 * w * u * v);
	}

	void RestLengths( int[] links, float[] rest )
	{
		for ( int v = 0; v < Pos.Length; v++ )
		{
			for ( int k = 0; k < 6; k++ )
			{
				int o = links[v * 6 + k];
				rest[v * 6 + k] = o < 0 ? 0 : (_attach[v] - _attach[o]).Length;
			}
		}
	}

	/// <summary>One Jacobi pass of distanceconstraint.frag (over-relaxed by 1.5).</summary>
	void Relax( int[] links, float[] rest, bool respectFixed )
	{
		// Flat float copies keep the hot loop free of struct shuffling.
		int n = Pos.Length;
		var px = _px; var py = _py; var pz = _pz;
		for ( int v = 0; v < n; v++ ) { px[v] = Pos[v].X; py[v] = Pos[v].Y; pz[v] = Pos[v].Z; }

		for ( int v = 0; v < n; v++ )
		{
			float s = Strength[v];
			if ( s == 0 || (respectFixed && _fixed[v]) ) { _temp[v] = Pos[v]; continue; }

			float cx = px[v], cy = py[v], cz = pz[v];
			float ox = 0, oy = 0, oz = 0;
			int count = 0;
			int baseIdx = v * 6;
			for ( int k = 0; k < 6; k++ )
			{
				int o = links[baseIdx + k];
				if ( o < 0 ) continue;
				count++;
				float dx = px[o] - cx, dy = py[o] - cy, dz = pz[o] - cz;
				float dist = MathF.Sqrt( dx * dx + dy * dy + dz * dz );
				if ( dist < 1e-9f ) continue;
				float f = 1.5f * (dist - rest[baseIdx + k]) / dist;
				ox += dx * f; oy += dy * f; oz += dz * f;
			}
			float scale = count > 0 ? s / count : 0;
			_temp[v] = new Vec3( cx + ox * scale, cy + oy * scale, cz + oz * scale );
		}
		(Pos, _temp) = (_temp, Pos);
	}

	/// <summary>Normal from the ordered one-ring (normals.frag).</summary>
	void ComputeNormals()
	{
		for ( int v = 0; v < Pos.Length; v++ )
		{
			var c = Pos[v];
			var sum = Vec3.Zero;
			int count = 0;
			for ( int k = 0; k < 6; k++ ) if ( _neighbours[v * 6 + k] >= 0 ) count++;
			for ( int k = 0; k < count; k++ )
			{
				var a = Pos[_neighbours[v * 6 + k]] - c;
				var b = Pos[_neighbours[v * 6 + (k + 1) % count]] - c;
				sum += Vec3.Cross( a, b );
			}
			float len = sum.Length;
			Normals[v] = len > 0 ? -sum / len : Vec3.Zero;
		}
	}
}