GAMES 103 動畫基礎作業1 Shape Matching 淺淺解析

簡介

作業1簡單實現了一個以一定初始速度和角速度的模型和牆壁碰撞的效果.
總共講解了三種算法

1. impulse (脈沖法)

2. Shape Matching(基於形狀保持的算法, 不包含物理特性)

3. Penalty methods

Shape Matching 可以說是最簡單的方法之一

因為完全不涉及角速度. 基礎邏輯理論就是, 當模型和牆壁發生碰撞的時候.這些碰撞的粒子會產生相反的形變但是. 我們強制將其整理保持形狀的一致性. 就是先形變. 再變回來.
由於形變的步驟. 是中間步驟. 不會展示出來. 所以效果還是挺好的.

Image

TIPS

code

using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class Rigid_Bunny_by_Shape_Matching : MonoBehaviour
{
	public bool launched = false;
	Vector3[] X; // world position
	Vector3[] Y; // temp world position
	Vector3[] Q; // Local coordinates
	Vector3[] V; // speed
	Vector3[] OldX;
	Matrix4x4 QQt = Matrix4x4.zero;
	Vector3 G = new Vector3(0.0f, -9.8f, 0.0f);
	float linear_decay = 0.999f;
	Vector3 ground = new Vector3(0, 0.01f, 0);
	Vector3 groundNormal = new Vector3(0, 1, 0);
	Vector3 wall = new Vector3(2.01f, 0, 0);
	Vector3 wallNormal = new Vector3(-1, 0, 0);
	float mu_T = 0.5f;  // μ_T may be coefficient of air resistance
	float mu_N = 5.0f; // μ_N may be Coefficient of Restitution
	float m_timer = 0;


	// Start is called before the first frame update
	void Start()
    {
    	Mesh mesh = GetComponent<MeshFilter>().mesh;
        V = new Vector3[mesh.vertices.Length];
        X = mesh.vertices;
		Y = mesh.vertices;
		OldX = mesh.vertices;
        Q = mesh.vertices;

        //Centerizing Q.
        Vector3 c=Vector3.zero;
        for(int i=0; i<Q.Length; i++)
        	c+=Q[i];
        c/=Q.Length;
		Debug.Log("C: " + c);
        for(int i=0; i<Q.Length; i++)
        	Q[i]-=c;

        //Get QQ^t ready.
		for(int i=0; i<Q.Length; i++)
		{
			QQt[0, 0]+=Q[i][0]*Q[i][0];
			QQt[0, 1]+=Q[i][0]*Q[i][1];
			QQt[0, 2]+=Q[i][0]*Q[i][2];
			QQt[1, 0]+=Q[i][1]*Q[i][0];
			QQt[1, 1]+=Q[i][1]*Q[i][1];
			QQt[1, 2]+=Q[i][1]*Q[i][2];
			QQt[2, 0]+=Q[i][2]*Q[i][0];
			QQt[2, 1]+=Q[i][2]*Q[i][1];
			QQt[2, 2]+=Q[i][2]*Q[i][2];
		}
		QQt[3, 3]=1;

		for(int i=0; i<X.Length; i++)
			V[i][0]=4.0f;

		Update_Mesh(transform.position, Matrix4x4.Rotate(transform.rotation), 0);
		transform.position=Vector3.zero;
		transform.rotation=Quaternion.identity;
   	}

	Matrix4x4 vector3x1dotvector1x3(Vector3 A, Vector3 B)
    {
		Matrix4x4 rlt = Matrix4x4.zero;
		rlt[3, 3] = 1.0f;

		rlt[0, 0] = A[0] * B[0];
		rlt[0, 1] = A[0] * B[1];
		rlt[0, 2] = A[0] * B[2];

		rlt[1, 0] = A[1] * B[0];
		rlt[1, 1] = A[1] * B[1];
		rlt[1, 2] = A[1] * B[2];

		rlt[2, 0] = A[2] * B[0];
		rlt[2, 1] = A[2] * B[1];
		rlt[2, 2] = A[2] * B[2];

		return rlt;
    }

   	// Polar Decomposition that returns the rotation from F.
   	Matrix4x4 Get_Rotation(Matrix4x4 F)
	{
		Matrix4x4 C = Matrix4x4.zero;
	    for(int ii=0; ii<3; ii++)
	    for(int jj=0; jj<3; jj++)
	    for(int kk=0; kk<3; kk++)
	        C[ii,jj]+=F[kk,ii]*F[kk,jj];
	   
	   	Matrix4x4 C2 = Matrix4x4.zero;
		for(int ii=0; ii<3; ii++)
	    for(int jj=0; jj<3; jj++)
	    for(int kk=0; kk<3; kk++)
	        C2[ii,jj]+=C[ii,kk]*C[jj,kk];
	    
	    float det    =  F[0,0]*F[1,1]*F[2,2]+
	                    F[0,1]*F[1,2]*F[2,0]+
	                    F[1,0]*F[2,1]*F[0,2]-
	                    F[0,2]*F[1,1]*F[2,0]-
	                    F[0,1]*F[1,0]*F[2,2]-
	                    F[0,0]*F[1,2]*F[2,1];
	    
	    float I_c    =   C[0,0]+C[1,1]+C[2,2];
	    float I_c2   =   I_c*I_c;
	    float II_c   =   0.5f*(I_c2-C2[0,0]-C2[1,1]-C2[2,2]);
	    float III_c  =   det*det;
	    float k      =   I_c2-3*II_c;
	    
	    Matrix4x4 inv_U = Matrix4x4.zero;
	    if(k<1e-10f)
	    {
	        float inv_lambda=1/Mathf.Sqrt(I_c/3);
	        inv_U[0,0]=inv_lambda;
	        inv_U[1,1]=inv_lambda;
	        inv_U[2,2]=inv_lambda;
	    }
	    else
	    {
	        float l = I_c*(I_c*I_c-4.5f*II_c)+13.5f*III_c;
	        float k_root = Mathf.Sqrt(k);
	        float value=l/(k*k_root);
	        if(value<-1.0f) value=-1.0f;
	        if(value> 1.0f) value= 1.0f;
	        float phi = Mathf.Acos(value);
	        float lambda2=(I_c+2*k_root*Mathf.Cos(phi/3))/3.0f;
	        float lambda=Mathf.Sqrt(lambda2);
	        
	        float III_u = Mathf.Sqrt(III_c);
	        if(det<0)   III_u=-III_u;
	        float I_u = lambda + Mathf.Sqrt(-lambda2 + I_c + 2*III_u/lambda);
	        float II_u=(I_u*I_u-I_c)*0.5f;
	        
	        
	        float inv_rate, factor;
	        inv_rate=1/(I_u*II_u-III_u);
	        factor=I_u*III_u*inv_rate;
	        
	       	Matrix4x4 U = Matrix4x4.zero;
			U[0,0]=factor;
	        U[1,1]=factor;
	        U[2,2]=factor;
	        
	        factor=(I_u*I_u-II_u)*inv_rate;
	        for(int i=0; i<3; i++)
	        for(int j=0; j<3; j++)
	            U[i,j]+=factor*C[i,j]-inv_rate*C2[i,j];
	        
	        inv_rate=1/III_u;
	        factor=II_u*inv_rate;
	        inv_U[0,0]=factor;
	        inv_U[1,1]=factor;
	        inv_U[2,2]=factor;
	        
	        factor=-I_u*inv_rate;
	        for(int i=0; i<3; i++)
	        for(int j=0; j<3; j++)
	            inv_U[i,j]+=factor*U[i,j]+inv_rate*C[i,j];
	    }
	    
	    Matrix4x4 R=Matrix4x4.zero;
	    for(int ii=0; ii<3; ii++)
	    for(int jj=0; jj<3; jj++)
	    for(int kk=0; kk<3; kk++)
	        R[ii,jj]+=F[ii,kk]*inv_U[kk,jj];
	    R[3,3]=1;
	    return R;
	}

	// Update the mesh vertices according to translation c and rotation R.
	// It also updates the velocity.
	void Update_Mesh(Vector3 c, Matrix4x4 R, float inv_dt)
   	{
   		for(int i=0; i<Q.Length; i++)
		{
			Vector3 x=(Vector3)(R*Q[i])+c;

			V[i] = (x-X[i])*inv_dt;
			X[i] = x;
		}	
		Mesh mesh = GetComponent<MeshFilter>().mesh;
		mesh.vertices=X;
   	}

	void Collision(float inv_dt)
	{
		for(int i=0; i<Q.Length; i++)
        {
			if(Vector3.Dot(X[i] - ground, groundNormal) < 0 && Vector3.Dot(V[i], groundNormal) < 0)// collision with ground
            {
				Vector3 VN = Vector3.Dot(V[i], groundNormal) * groundNormal;
				Vector3 VT = V[i] - VN;
				float a = Mathf.Max(0, 1.0f - mu_T * (1.0f + mu_N)) * Vector3.Magnitude(VN) / Vector3.Magnitude(VT);
				V[i] = -1.0f * mu_N * VN + 2.0f * a * VT;
            }
			else if(Vector3.Dot(X[i] - wall, wallNormal) < 0 && Vector3.Dot(V[i], wallNormal) < 0) // collision with wall
            {
				Vector3 VN = Vector3.Dot(V[i], wallNormal) * wallNormal;
				Vector3 VT = V[i] - VN;
				float a = Mathf.Max(0, 1.0f - mu_T * (1.0f + mu_N)) * Vector3.Magnitude(VN) / Vector3.Magnitude(VT);
				V[i] = -1.0f * mu_N * VN + 2.0f * a * VT;
			}
        }
	}

    // Update is called once per frame
    void Update()
    {

		if (Input.GetKey("l"))
        {
			launched = true;
			for(int i=0; i<V.Length; i++)
            {
				V[i] = new Vector3(5.0f, 2.0f, 0.0f);
			}
        }
		if (Input.GetKey("r"))
		{
			launched = false;
			for (int i = 0; i < V.Length; i++)
			{
				V[i] = new Vector3(4.0f, 0.0f, 0.0f);
			}

			Update_Mesh(new Vector3(0, 0.6f, 0), Matrix4x4.Rotate(transform.rotation), 0);
		}
		if(!launched)
        {
			return;
        }
		//m_timer += Time.time;
		//if (m_timer <= 500)
		//{
		//	return;
		//}
		//else
		//{
		//	m_timer = 0;
		//}


		float dt = 0.015f;

  		//Step 1: run a simple particle system.
        for(int i=0; i<V.Length; i++)
        {
			V[i] = V[i] + G * dt;
			V[i] *= linear_decay;
        }

        //Step 2: Perform simple particle collision.
		Collision(1/dt);

		// Step 3: Use shape matching to get new translation c and 
		// new rotation R. Update the mesh by c and R.
        //Shape Matching (translation)
		for(int i=0; i<V.Length; i++)
        {
			Y[i] = X[i] + V[i] * dt;
        }

		// calc c
		Vector3 c = new Vector3(0, 0, 0); 
		for(int i=0; i<V.Length; i++)
        {
			c += Y[i];
        }
		c = c / V.Length;

		// calc A
		Matrix4x4 A = Matrix4x4.zero;
		A[3, 3] = 1.0f;
		for (int i=0; i<V.Length; i++)
        {
			Matrix4x4 o = vector3x1dotvector1x3(Y[i] - c, Q[i]);
			A[0, 0] += o[0, 0];
			A[0, 1] += o[0, 1];
			A[0, 2] += o[0, 2];
			A[1, 0] += o[1, 0];
			A[1, 1] += o[1, 1];
			A[1, 2] += o[1, 2];
			A[2, 0] += o[2, 0];
			A[2, 1] += o[2, 1];
            A[2, 2] += o[2, 2];
		}
		A = A * QQt.inverse;

		//Shape Matching (rotation)
		// calc R
		Matrix4x4 R = Matrix4x4.zero;
		R = Get_Rotation(A);

		Update_Mesh(c, R, 1/dt);

	}
}

TRICK

PIPELINE

計算 C (模型質心坐標)

計算 A (通過A計算得到R, 即 旋轉矩陣) == 從A得到R. 老師已經提供了.

更新頂點坐標

Trick

在處理碰撞的函數中. 里面的參數是自己調整的參數. 比如 5.0f 邏輯上 μ_T μ_N 都應該是一個小於1的數. 但是, 測試后感覺效果不是特別好. 如果你知道為什么的話, 請留言.

請看完bilibli 4個視屏后開始自己的作業.

參考鏈接

配置VS 作為 Unity 的配置環境
https://blog.csdn.net/qq_34405576/article/details/105572069

源碼參考師弟寫出來的, 大部分是抄的嘿嘿~~ 侵權刪除.