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;
}
}
}