Games103的lab1的学习笔记和代码实现
课程链接: GAMES103-基于物理的计算机动画入门_哔哩哔哩_bilibili
课程主页:GAMES103:基于物理的计算机动画入门 – 计算机图形学与混合现实在线平台
1. 理论基础
1.1 点乘与叉乘
1.2 动量与角动量
转动惯量是刚体绕轴转动时惯性的量度,可以理解为线性动力学中的质量,是一个标量,不需要转换空间。
1.4 速度与角速度
1.3 四元数及计算
为了避免万向节死锁问题(某种情况下降低旋转自由度),引入了四元数的概念。
四元数是一个超复数,其(x,y,z,w)四个值分为三个虚部,一个实部。
与欧拉角的转换可以参考这个文章3D数学基础:四元数与欧拉角之间的转换 - rainbow70626 - 博客园
即xyz三个值是在xyz轴的分量,确定了一个旋转轴的方向,然后w是绕旋转轴旋转的角度。
1.5 坐标空间转换
2. 算法实现
2.1 算法思路
本质上刚体发生变化的就是位置(向量)和旋转(四元数)。
碰撞过程导致了受力变化,从而影响物体速度,进而更新物体的位置坐标和旋转。本算法是当前帧如果发生碰撞,计算出物体局部受到冲量大小(没有发生但是假设发生),从而求出整个物体的速度变化,然后对位置更新(更新当前帧)。优点是可以及时对受力进行反馈。
所以如果发生碰撞,需要对于局部受碰撞点,受到弹力和摩擦力作用,计算出之前的和新的速度大小,得到局部冲量大小。
2.2 刚体碰撞
2.2.1 碰撞检测
检测方法是非常直观的,计算点到表面的距离,大于则在外面,小于在里面。难点是对于复杂网格表面,应该如何计算。
遍历每个点(世界坐标),计算到碰撞体表面的距离(PXi与法线点乘),负数表示碰撞。
2.2.2 碰撞反馈
对于局部受碰撞点:
计算法向和切向速度,法向受弹力作用:速度大小近似表示为法线的反向乘弹力系数
切向受到摩擦力作用:按照a系数减少,但是不能小于0(速度不会反方向回去)
课程表示a系数需要满足库伦定律,具体a的计算如下图所示。
对于受碰撞的物体:
如果我们知道冲量大小,那么就可以求到新的速度和角速度的值。
所以我们需要对局部碰撞点求得冲量大小:
得到计算式如下:
星号表示向量的叉乘可以用一个新的矩阵表示,简化了很多(tql)
位置和旋转更新
知道速度和角速度后我们就可以很快求出新的位置和旋转。
注意问题
1. 碰撞点不止一个,可能有几百上千个,所以我们求的碰撞点是它们的平均坐标
2. 碰撞时候会因为一直有向下的重力加速度,所以弹力也一直存在,导致抖动问题,因此在趋于稳定的时候,将弹力减小,使得速度尽可能控制为0
3. unity里面没有四元数的加法,unity中的四元数表示如下:
具体代码
using UnityEngine;
using System.Collections;
public class Rigid_Bunny : MonoBehaviour
{
bool launched = false;
float dt = 0.015f;
Vector3 v = new Vector3(0, 0, 0); // velocity
Vector3 w = new Vector3(0, 0, 0); // angular velocity
Vector3 gravity = new Vector3(0,-9.8f,0);
float mass; // mass
Matrix4x4 I_ref; // reference inertia
float linear_decay = 0.999f; // for velocity decay
float angular_decay = 0.98f;
float restitution = 0.3f; // for collision
float friction = 0.2f; // 摩擦系数
bool lowCollide = false;
Mesh mesh;
// Use this for initialization
void Start ()
{
transform.rotation = new Quaternion(0, Mathf.Sin(Mathf.Deg2Rad * 30) * 1, 0,Mathf.Sin(Mathf.Deg2Rad * 30));
mesh = GetComponent<MeshFilter>().mesh;
Vector3[] vertices = mesh.vertices;
float m=1;
mass=0;
for (int i=0; i<vertices.Length; i++) // reference inertias
{
mass += m;
float diag=m*vertices[i].sqrMagnitude;
I_ref[0, 0]+=diag;
I_ref[1, 1]+=diag;
I_ref[2, 2]+=diag;
I_ref[0, 0]-=m*vertices[i][0]*vertices[i][0];
I_ref[0, 1]-=m*vertices[i][0]*vertices[i][1];
I_ref[0, 2]-=m*vertices[i][0]*vertices[i][2];
I_ref[1, 0]-=m*vertices[i][1]*vertices[i][0];
I_ref[1, 1]-=m*vertices[i][1]*vertices[i][1];
I_ref[1, 2]-=m*vertices[i][1]*vertices[i][2];
I_ref[2, 0]-=m*vertices[i][2]*vertices[i][0];
I_ref[2, 1]-=m*vertices[i][2]*vertices[i][1];
I_ref[2, 2]-=m*vertices[i][2]*vertices[i][2];
}
I_ref [3, 3] = 1;
}
Matrix4x4 Get_Cross_Matrix(Vector3 a)
{
//Get the cross product matrix of vector a
Matrix4x4 A = Matrix4x4.zero;
A [0, 0] = 0;
A [0, 1] = -a [2];
A [0, 2] = a [1];
A [1, 0] = a [2];
A [1, 1] = 0;
A [1, 2] = -a [0];
A [2, 0] = -a [1];
A [2, 1] = a [0];
A [2, 2] = 0;
A [3, 3] = 1;
return A;
}
Matrix4x4 Matrix_Subtract(Matrix4x4 M1,Matrix4x4 M2)
{
return new Matrix4x4(
M1.GetColumn(0) - M2.GetColumn(0),
M1.GetColumn(1) - M2.GetColumn(1),
M1.GetColumn(2) - M2.GetColumn(2),
M1.GetColumn(3) - M2.GetColumn(3)
);
}
Matrix4x4 MatrixMultiplyFloat(Matrix4x4 M, float f)
{
for (int i = 0;i < 4; i++)
for (int j = 0;j < 4; j++)
{
M[i, j] *= f;
}
return M;
}
// In this function, update v and w by the impulse due to the collision with
// a plane <P, N>
void Collision_Impulse(Vector3 P, Vector3 N)
{
// find vaild vertices collide
Vector3[] vertices = mesh.vertices;
int num = vertices.Length,validNum = 0;
Vector3 x = transform.position;
Vector3 avgPos = Vector3.zero;
Matrix4x4 R = Matrix4x4.Rotate(transform.rotation);// rotate matrix
for(int i = 0; i < num; i++)
{
Vector3 x_i = x + R.MultiplyVector(vertices[i]);
// assume that infinity plane
Vector3 p2xi = x_i - P;
if(Vector3.Dot(p2xi,N) < 0)// inside
{
// is considerable
validNum++;
avgPos += vertices[i];
}
}
// calculate collisional effect
if(validNum > 0)
{
// campute the wanted vi_new
avgPos /= validNum;
Vector3 r_i = avgPos;
Vector3 v_i = v + Vector3.Cross(w, R.MultiplyVector(r_i));
if(Vector3.Dot(v_i,N) > 0)//the velocity is positive
return;
Vector3 x_i = x + R.MultiplyVector(r_i);
float distance = Vector3.Dot(x_i, N);
print("distance" + distance + "posY" + x_i.y);
if (-0.3f < distance)
{
lowCollide = true;
}
Vector3 vi_normal = Vector3.Dot(v_i,N) * N;
Vector3 vi_tangent = v_i - vi_normal;
Vector3 vi_normal_new = -restitution * vi_normal;
float a = Mathf.Max(1- friction * (1+restitution)*(vi_normal.magnitude/vi_tangent.magnitude),0);
Vector3 vi_tangent_new = a * vi_tangent;
Vector3 vi_new = vi_normal_new + vi_tangent_new;
// compute the impulse
//Matrix4x4 I_world = R * I_ref * Matrix4x4.Transpose(R);// inertia is a quantity
Matrix4x4 I_inv = I_ref.inverse;
Matrix4x4 Rc = Get_Cross_Matrix(R * r_i);
Matrix4x4 K = Matrix_Subtract(MatrixMultiplyFloat(Matrix4x4.identity,1/mass) , Rc * I_inv * Rc);
Vector3 j = K.inverse * (vi_new - v_i);
// update v and w
Vector3 dv = j / mass;
Debug.DrawRay(x + R.MultiplyVector(r_i), dv, Color.red);// draw the velocity direction
//Debug.DrawRay(x,v,Color.green);// gobal velocity
v += dv;
w = w + (Vector3)(I_inv * Vector3.Cross(R * r_i , j));
}
}
Quaternion QuaternionAdd(Quaternion q1,Quaternion q2)
{
q1.x += q2.x;
q1.y += q2.y;
q1.z += q2.z;
q1.w += q2.w;
return q1;
}
// Update is called once per frame
void Update ()
{
//Game Control
if(Input.GetKey(KeyCode.R))
{
transform.position = new Vector3 (0, 2f, 0);
restitution = 0.5f;
launched=false;
}
if(Input.GetKey(KeyCode.S))// shoot
{
v = new Vector3 (5, 2, 0);
launched=true;
}
// Part I: Update velocities
v += gravity * dt;
v *= linear_decay;
w *= angular_decay;
// Part II: Collision Impulse
lowCollide = false;
Collision_Impulse(new Vector3(0, 0.01f, 0), new Vector3(0, 1, 0));
Collision_Impulse(new Vector3(2, 0, 0), new Vector3(-1, 0, 0));
if (lowCollide && v.magnitude < 0.13f)
{
restitution *= 0.5f;
if (restitution < 0.01f)
{
restitution = 0;
}
if(v.magnitude < 0.001f)
{
v = Vector3.zero;
}
}
else
{
restitution = 0.5f;
}
v *= linear_decay;
w *= angular_decay;
// Part III: Update position & orientation
//Update linear status
Vector3 x = transform.position;
x += v * dt;
//Update angular status
Quaternion q = transform.rotation;
Quaternion qw = new Quaternion(w.x * dt/2, w.y * dt / 2, w.z * dt / 2, 0);
q = QuaternionAdd(q , qw * q);
// Part IV: Assign to the object
transform.position = x;
transform.rotation = q;
}
}