Unity

线性代数之矩阵

矩阵基本计算

Posted by Bob on October 29, 2018

矩阵用于向量从一个坐标系到另一个坐标系变换

坐标系

做算法推算之前一定要明确坐标系类型！

3D坐标系：左手坐标系、右手坐标系,献出双手亲自示范~~~

image

做算法推算之前一定要明确方向！

image

右手法则：大拇指z正轴朝向自己，为旋转轴正朝向，另外四个手指逆时针弯曲的方向为旋转正方向。

世界坐标系(Unity左手坐标系)：又称全局坐标系或者宇宙坐标系，整个游戏世界中，物体没有父节点transform.position可以获得物体世界坐标。

物体坐标系(Unity左手坐标系)：每个物体都有自己特定的坐标系

惯性坐标系：惯性坐标系的原点与物体坐标系的原点重合，惯性坐标系的轴平行于世界坐标系的轴。

image

为什么要引入惯性坐标系？因为从物体坐标系转换到惯性坐标系只需旋转，从惯性坐标系转换到世界坐标系只需平移，把复杂计算的分解成二步简单的。

image

观察坐标系(Unity模型和世界空间用左手坐标系，观察空间用右手坐标系)：以camera为原点，camera红色x+轴是观察空间x+方向，camera绿色y+轴是观察空间y+方向，camera蓝色z+轴是观察空间的蓝色z-方向。观察空间内z轴越小离摄像机越远，depth越大。

屏幕坐标

屏幕坐标是以像素来定义的，它的范围是以左下角为(0,0)，右上角为(Screen.width,Screen.height)定义的这样一个矩形。Z轴是以相机的世界坐标来衡量的。

视口坐标

视口坐标是标准化后的屏幕坐标。它的范围是以左下角为(0,0),右上角为(1,1)定义的这样一个矩形。Z轴是以相机的世界坐标来衡量的。

矩阵

方块矩阵、正方阵(square matrix)：行和列数目相等的矩阵，如unity shader的3x3,4x4的矩阵

对角矩阵(diagonal matrix)：除对角外所有元素都为0，对角线只有一条

单位矩阵(identity matrix)：对角都是1的对角矩阵，任意一个矩阵乘以单位矩阵，都将得到原来的矩阵

image

在unity中Matrix4x4表达

namespace UnityEngine
{

// m00 m01 m02 m03

// m10 m11 m12 m13

// m20 m21 m22 m23

// m30 m31 m32 m33


	public struct Matrix4x4
	{
		public float m00;
		public float m10;
		public float m20;
		public float m30;
		public float m01;
		public float m11;
		public float m21;
		public float m31;
		public float m02;
		public float m12;
		public float m22;
		public float m32;
		public float m03;
		public float m13;
		public float m23;
		public float m33;
      }

// m0   m4   m8	  m12	

// m1   m5   m9   m13		

// m2   m6   m10  m14		

// m3   m7   m11  m15	

	  public float this[int row, int column]
	  {
			get
			{
				return this[row + column * 4];
			}
			set
			{
				this[row + column * 4] = value;
			}
	  }
}

转置矩阵(transposed matrix)

行和列翻转一下得出转置矩阵。一个矩阵转置再转置得到原矩阵

image

运算性质,A和B都是矩阵，a是常数

image

矩阵加减法

相加减的两个矩阵，大小必须一致为mxn

运算性质,A、B、C都是矩阵

A + B = B + A

A + (B + C) = (A + B) + C

image

矩阵乘标量

image

运算性质,A和B都是矩阵，a和b是常数

a(A + B)= aA + aB

a (A * B) = aA * B = A * aB

(a+b)A = aA + bA

a(bA)= (ab)A

一般矩阵乘法(Matmul product)

image

一个mxn的矩阵和nxp的矩阵才能相乘,得到一个mxp的矩阵,下图2行3列矩阵乘以3行2列的矩阵。矩阵乘法不满足交换律，因为交换后不一定能相乘了。这里实际是线性方程。

例如计算结果为第一个矩阵第一行和第二个矩阵的第一列依次相乘求和。 2x0+5x1+7x2 = 19,在经济学中可理解为2块的苹果买了0个+5块的葡萄买了1串+7块的西瓜买了2个=共花了19块钱。

image

矩阵乘法不满足交换律

运算性质,A、B、C都是矩阵

乘法结合律： (AB)C=A(BC)

乘法左分配律：(A+B)C=AC+BC

乘法右分配律：C(A+B)=CA+CB

法线变换

image

法线变换：应该用变换矩阵的逆转置矩阵

T为切线,N为法线。有T = P2-P1,变换后T’= MP2 - MP1 = M(P2-P1) = MT

image

UnityObjectToWorldNormal源码如下:


// Transforms normal from object to world space
inline float3 UnityObjectToWorldNormal( in float3 norm )
{
#ifdef UNITY_ASSUME_UNIFORM_SCALING
    return UnityObjectToWorldDir(norm);
#else
    // mul(IT_M, norm) => mul(norm, I_M) => {dot(norm, I_M.col0), dot(norm, I_M.col1), dot(norm, I_M.col2)}
    return normalize(mul(norm, (float3x3)unity_WorldToObject));
#endif
}
 
// Transforms direction from object to world space
inline float3 UnityObjectToWorldDir( in float3 dir )
{
    return normalize(mul((float3x3)unity_ObjectToWorld, dir));

}

简要分析如下:

  1. 如果模型是等比缩放，则把法线向量转化到世界坐标系上的法线向量时，直接调用 UnityObjectToWorldDir(norm); 这个方法里，unity_ObjectToWorld这个矩阵和物体坐标系下的法线向量相乘，顾名思义，这个矩阵就是物体坐标系转化到世界坐标系的变换矩阵

  2. 如果模型是非等比缩放，如下： mul(norm, (float3x3)unity_WorldToObject) 在等比缩放模式下unity_ObjectToWorld是M矩阵，在非等比缩放模式下unity_WorldToObject则是M的逆矩阵，同样一个向量和矩阵相乘，根据矩阵和向量相乘规则，mul(v,M)和mul(M,v)是在内部做了 一个转置处理。这里的

image

阿达马乘积(Hadamard product)

阿达马乘积记作C=A○B或A*B，常用于shader中颜色的运算

image

克罗内克积(Kronecker product)

Kronecker积是两个任意大小矩阵间的运算，表示为 A ⊗ B。如果A是一个 m x n 的矩阵，而B是一个 p x q 的矩阵，克罗内克积则是一个 mp x nq 的矩阵。克罗内克积也称为直积或张量积，以德国数学家利奥波德·克罗内克命名。

image

矩阵行列式(determinant of a matrix)

一阶行列式

丨A丨= A

二阶行列式

image

三阶伴随矩阵和行列式

image

A为n阶矩阵，丨A丨是A的行列式(注意不是绝对值,也不是向量求模)，又记为det(A),取值为一个标量。当n=1时，D=丨a11丨= a11;当n>=2时，D=a11A11+a12A12+…+a1nA1n;其几何意义参见向量叉乘使用三阶行列式。

A*是指矩阵A的伴随矩阵(adjoint matrix),记作adj(A)，是由A的元素的代数余子式按照交换行列标的顺序构成的同级矩阵。

伴随矩阵的一些基本性质如下：

  • A可逆当且仅当A*可逆；
  • 如果A可逆，则image,当丨A丨 = 0 时，A的逆矩阵不存在。
四阶行列式
double Matrix4x4f::GetDeterminant () const
{
	double m00 = Get(0, 0);  double m01 = Get(0, 1);  double m02 = Get(0, 2);  double m03 = Get(0, 3);
	double m10 = Get(1, 0);  double m11 = Get(1, 1);  double m12 = Get(1, 2);  double m13 = Get(1, 3);
	double m20 = Get(2, 0);  double m21 = Get(2, 1);  double m22 = Get(2, 2);  double m23 = Get(2, 3);
	double m30 = Get(3, 0);  double m31 = Get(3, 1);  double m32 = Get(3, 2);  double m33 = Get(3, 3);

	double result =
		m03 * m12 * m21 * m30 - m02 * m13 * m21 * m30 - m03 * m11 * m22 * m30 + m01 * m13 * m22 * m30 +
		m02 * m11 * m23 * m30 - m01 * m12 * m23 * m30 - m03 * m12 * m20 * m31 + m02 * m13 * m20 * m31 +
		m03 * m10 * m22 * m31 - m00 * m13 * m22 * m31 - m02 * m10 * m23 * m31 + m00 * m12 * m23 * m31 +
		m03 * m11 * m20 * m32 - m01 * m13 * m20 * m32 - m03 * m10 * m21 * m32 + m00 * m13 * m21 * m32 +
		m01 * m10 * m23 * m32 - m00 * m11 * m23 * m32 - m02 * m11 * m20 * m33 + m01 * m12 * m20 * m33 +
		m02 * m10 * m21 * m33 - m00 * m12 * m21 * m33 - m01 * m10 * m22 * m33 + m00 * m11 * m22 * m33;
	return result;
}

逆矩阵(inverse matrix)

逆矩阵必须是方块矩阵，不是所有矩阵都有逆矩阵

设A为n阶矩阵，若存在n阶矩阵B使得:AB=BA=I(单位矩阵)，则称A是可逆的且矩阵B是矩阵A的逆矩阵。并称A是一个非奇异矩阵，如果不存在B矩阵，则A是一个奇异矩阵。

image

当矩阵A可逆，则有下图(其中I是单位矩阵)

image

根据矩阵乘法的定义，单位矩阵I的重要性质为：AI = A 和 IB = B

求逆的具体算法参见intel库中引用方法invert_matrix_general或Trapas3D matrix,Unity的c++内部实现也是来自该算法(作者Jacques Leroy)。

正交矩阵(orthogonal matrix)

一个方正矩阵M和它的转置矩阵的乘积是单位矩阵，这个方正矩阵M就是正交的。并且这个它的转置矩阵和逆矩阵是相等的。

image

M矩阵各行是单位向量(点积是1)且两两正交(点积是0)

M矩阵各列是单位向量且两两正交

齐次坐标(Homogeneous Coordinate)

在欧式几何原理中二条平行线不能相交。但是在投影空间中不是这样，物体二条平行边在远处会相交，变成一个无穷远的点。这样的点在欧几里德空间中无法处理。

August Ferdinand Möbius引入齐次坐标可以在投影空间中计算图形和几何，齐次坐标是用N+1个数表示N维坐标的一种方式。

2d齐次坐标中，引入一个分量w，由此笛卡尔坐标中的点（X，Y）在齐次坐标中为（x,y,w）,那么有X = x/w,Y = y/w,例如笛卡尔坐标中的点（1，2）在齐次坐标中为（1,2,1），如果点移动到无穷远在笛卡尔坐标中为(∞,∞) ，在齐次坐标中她为(1,2,0), 由此（1 / 0,2 / 0）≈（∞，∞）。我们可以在不使用“∞”的情况下表达无穷远处的点。笛卡尔坐标系就是齐次坐标系中w=1的那个平面，(x,y,1)是齐次坐标(kx,ky,k)表示的点在w=1上的映射。

为什么叫齐次？齐次坐标（x,y,w）转化为笛卡尔坐标(x/w,y/w），我们可以简单地将x和y除以w。

（1,2,3） => (1/3,2/3)

（2,4,6） => (1/3,2/3)

（4,8,12）=> (1/3,2/3)

齐次坐标中（1,2,3）、2,4,6）、（4,8,12）都对应笛卡尔坐标中(1/3,2/3)，这些点具有尺度不变性（Scale Invariant），是“齐性的”(同族的)，所以称之为齐次坐标。

齐次坐标表示是计算机图形学的重要手段之一，它既能够用来明确区分向量和点，同时也更易用于进行仿射(线性)几何变换。”—— F.S. Hill, JR。

a.点是三维空间中的某个坐标，是绝对的，它的值是参照原点的。

b.向量用于表示力和速度等具有方向和大小的量， 通常用具有长度和方向的线段来表示

image

在普通坐标(Ordinary Coordinate)和齐次坐标(Homogeneous Coordinate)之间进行转换：

(1)从普通坐标转换成齐次坐标时,齐次坐标就是将一个原本是n维的向量用一个n+1维向量来表示

一个普通坐标的点P=(x, y,z)，有对应的多个不唯一齐次坐标(wx, wy, wz, w)，其中w不等于零。比如当w取不同值时，P(1, 2, 3)的齐次坐标有(1, 2, 3, 1)、(1, 4, 6, 2)、(-0.1, -0.2, -0.3, -0.1)等。当w=1时产生的齐次坐标叫“规格化坐标”。

如果(x,y,z)是个点，则变为(x,y,z,1);

如果(x,y,z)是个向量，则变为(x,y,z,0)

(2)从齐次坐标转换成普通坐标时,一个三维坐标的三个分量x，y，z用齐次坐标表示为变为x，y，z，w的四维空间，变换成三维坐标是方式是x/w,y/w,z/w，当w为0时，在数学上代表无穷远点，即并非一个具体的坐标位置，而是一个具有大小和方向的向量。从而，通过w我们就可以用同一系统表示两种不同的量。

image

如果是(x,y,z,1)，则知道它是个点，变成(x,y,z);

如果是(x,y,z,0)，则知道它是个向量，仍然变成(x,y,z)

标量集合(scalar aggregates)、向量空间(Vector space)、仿射空间(affine space)

  • 标量集合 中的任何两个标量都可以经过加减乘除法这两种运算得到另一个标量，标量有实数等。

  • 向量空间，标量乘以向量得另一个向量，向量加法得另一个向量等。

方程L为线性变换(linear transformation)满足如下性质：

L(p+q) = L(p) + L(q)

aL(p) = L(ap)

p=(px,py,pz)和q=(qx,qy,qz)是任意3d向量，a为一个标量。

a.线性变换变换前是直线的，变换后依然是直线

b.线性变换直线比例保持不变

c.线性变换变换前是原点的，变换后依然是原点

例如：定义函数L(x,y,z)=(x^2+y^2+z^2),有L(1,2,3)=(1,4,9),函数为非线性函数，因为当a=2时，p=(1,2,3), L(ap)=L(2,4,6)=(4,16,36),而2L(p)=2(1,4,9)=(2,8,18).

平移变换不是线性变换

  • 仿射空间中，点与点之间做差可以得到向量，点与向量做加法将得到另一个点，但是点与点之间不可以做加法。

image

求三点构成平面上的一点P，a是标量，用仿射组合(Affine Combination)表示一个点。

image

将多个标量a，简化为一个标量t，

image

得到曲线公式

image

转化为矩阵计算

image

a. 一个任意的仿射变换都能表示为向量u乘以一个矩阵A(线性变换)接着再 加上一个向量b(平移).

image

b.仿射变形，其特征就是一切变形都不会破坏线条的线性。变形后水平和垂直方向上的长度比例可以发生变化。但直线永远不会变成曲线。坐标系内各点的变换都是均匀的，不存在局部扭曲和象限的塌缩。一对平行线，无论经过多少次仿射变形，都将保持平行，不会有交集。

c.仿射变形不考虑作为参照的原点.

变换

2d点

image

2d平移、缩放、旋转

平移

image

缩放

image

复合变换

image

旋转

px逆时针旋转到px’,py旋转旋转到py’

image

得出下式，矩阵R左乘旋转列向量

image

被动旋转和主动旋转

被动旋转(Alias)：坐标系旋转，点相对坐标系不变。

主动旋转(Alibi)：点自己主动旋转，坐标系不变.

上面的R(a)矩阵就表示逆时针旋转角度a，或坐标系旋转相同角度但相反方向(即顺时针方向)。

image

下列Unity Shader示例实现了一个UV顺时针旋转动画。

image

Shader "Unlit/UVRotation"
{
	Properties{_MainTex ("Texture", 2D) = "white" {}}
	SubShader
	{
		Tags { "RenderType"="Opaque" }
		LOD 100
		Blend SrcAlpha OneMinusSrcAlpha

		Pass
		{
			CGPROGRAM
			#pragma vertex vert
			#pragma fragment frag
			#include "UnityCG.cginc"

			struct appdata
			{
				float4 vertex : POSITION;
				float2 uv : TEXCOORD0;
			};

			struct v2f
			{
				float2 uv : TEXCOORD0;
				float4 vertex : SV_POSITION;
			};

			sampler2D _MainTex;
			float4 _MainTex_ST;
			
			v2f vert (appdata v)
			{
				v2f o;
				o.vertex = UnityObjectToClipPos(v.vertex);
				o.uv = TRANSFORM_TEX(v.uv, _MainTex);
				return o;
			}
			
			fixed4 frag (v2f i) : SV_Target
			{
				float2 defaultUV = i.uv;
				float2 offset = float2(0.5,0.5);				
				defaultUV.xy -= offset;
				if (length(defaultUV) > 0.5) {
					return fixed4(0,0,0,0);
				}
				float angle = _Time.y;

				float x = defaultUV.x*cos(angle) - defaultUV.y*sin(angle);
				float y = defaultUV.x*sin(angle) + defaultUV.y*cos(angle);
				i.uv.x = x;
				i.uv.y = y;
				i.uv.xy += offset;

				//float2x2 rotationMatrix = float2x2(cos(angle),-sin(angle),sin(angle),cos(angle));
				
				//i.uv.xy = mul(rotationMatrix,defaultUV) + offset;				
				
				fixed4 col = tex2D(_MainTex, i.uv);
				return col;
			}
			ENDCG
		}
	}
}

3d平移、缩放、旋转

点平移

根据矩阵乘法可得点平移。右手坐标系有OpenGL API，在空间中，右手拇指x正轴朝右，食指y正轴朝上，中指z正轴朝自己，则称为右手坐标系。 左手坐标系有Direct X API、Unity、OSG。

矩阵乘法的顺序与坐标系是左手系还是右手系有关系么？根本没啥关系！

DirectX、OSG、中使用的是行向量。


//row-major    position 行向量  矩阵在右，行向量在左   
 
mul(position, matrix)

Unity Shader、OpenGL、CG 中使用的是列向量。


//column-major  position 列向量 矩阵在左，列向量在右

o.pos = mul(UNITY_MATRIX_MVP, v.position);

HLSL：　存储方式和DirectX相反(column-major)

NVIDIA公司的CG(C for Graphic)或OpenG的GLSL(OpenGL Shading Language)转Direct3D 可以用转置实现：

mul(matrix, position) --->  mul(position, transpose(matrix))

矩阵T和点p(列矩阵)相乘记做Tp，列矩阵从右向左读，读列向量p右乘矩阵T得列向量,或矩阵T左乘p。

image

在unity中表达Matrix4x4


		public static Matrix4x4 Translate(Vector3 v)
		{
			return new Matrix4x4
			{
				m00 = 1f,m01 = 0f,m02 = 0f,m03 = v.x,
				m10 = 0f,m11 = 1f,m12 = 0f,m13 = v.y,
				m20 = 0f,m21 = 0f,m22 = 1f,m23 = v.z,
				m30 = 0f,m31 = 0f,m32 = 0f,m33 = 1f
			};
		}

		
		public Vector3 MultiplyPoint(Vector3 v)
		{
			Vector3 result;
			result.x = this.m00 * v.x + this.m01 * v.y + this.m02 * v.z + this.m03;
			result.y = this.m10 * v.x + this.m11 * v.y + this.m12 * v.z + this.m13;
			result.z = this.m20 * v.x + this.m21 * v.y + this.m22 * v.z + this.m23;
			float num = this.m30 * v.x + this.m31 * v.y + this.m32 * v.z + this.m33;
			num = 1f / num;
			result.x *= num;
			result.y *= num;
			result.z *= num;
			return result;
		}

		public Vector3 MultiplyPoint3x4(Vector3 v)
		{
			Vector3 result;
			result.x = this.m00 * v.x + this.m01 * v.y + this.m02 * v.z + this.m03;
			result.y = this.m10 * v.x + this.m11 * v.y + this.m12 * v.z + this.m13;
			result.z = this.m20 * v.x + this.m21 * v.y + this.m22 * v.z + this.m23;
			return result;
		}

		public Vector3 MultiplyVector(Vector3 v)
		{
			Vector3 result;
			result.x = this.m00 * v.x + this.m01 * v.y + this.m02 * v.z;
			result.y = this.m10 * v.x + this.m11 * v.y + this.m12 * v.z;
			result.z = this.m20 * v.x + this.m21 * v.y + this.m22 * v.z;
			return result;
		}

p(行矩阵)和矩阵相乘记做pT。行矩阵从左向右读，读行向量p左乘矩阵T得行向量。

image

旋转和缩放对于向量和点都有意义。平移变换只对于点才有意义，因为普通向量(例如法向量)没有位置概念，只有大小和方向,而平移变换是个加法运算(点r=[rx,ry,rz] ，平移t=[tx,ty,tz] ，r+t=[rx+tx,ry+ty,rz+tz])，不方便以后和旋转、缩放做乘法混合计算,所以引入齐次坐标区分，写成两个矩阵乘积形式的。使用齐次坐标使得仿射变换可以以统一的矩阵形式进行表示。

image

引入齐次坐标系表达 p̃ =(x,y,z,1)，（尺度不变性，实际上在高一维的空间映射到 w=1 平面, 这样计算后结果直接可导出到欧式空间）

点缩放

image

在unity中表达Matrix4x4


		public static Matrix4x4 Scale(Vector3 v)
		{
			return new Matrix4x4
			{
				m00 = v.x,m01 = 0f,m02 = 0f,m03 = 0f,
				m10 = 0f,m11 = v.y,m12 = 0f,m13 = 0f,
				m20 = 0f,m21 = 0f,m22 = v.z,m23 = 0f,
				m30 = 0f,m31 = 0f,m32 = 0f,m33 = 1f
			};
		}
在右手系点绕x轴逆时针旋转的矩阵

image

在右手系点绕y轴逆时针旋转的矩阵

image

在右手系点绕z轴逆时针旋转的矩阵

image

红点p逆时针旋转到绿点p’，可推导上图公式

image

已知p(x,y),z值保持不变，求旋转β度后的p’(x’,y’)，右手逆时针四指弯曲为旋转正方向，拇指为z轴正方向。

x = rcosα,y = rsinα

x’ = rcos(β + α),y’ = rsin(β + α)

根据二角和差公式

image

可得

x’ = rcosβcosα - rsinβsinα = xcosβ - ysinβ

y’ = rsinβcosα + rcosβsinα = xsinβ + ycosβ

这个线性方程就可以变换成下图的列向量右乘矩阵。

image

在左手系点绕x轴顺时针旋转的矩阵

image

在左手系点绕y轴顺时针旋转的矩阵

image

在左手系点绕z轴顺时针旋转的矩阵

image

推导:如果点P是绕z轴顺时针(z+向z-看)旋转a到P’.

image

x’ = rcos(θ+h),y’ = rsin(θ+h)

x = rcosh,y = rsinh

x’ = rcosθcosh - rsinθsinh = xcosθ - ysinθ

y’ = rsinθcosh + rcosθsinh = xsinθ + ycosθ

矩阵计算如下图，对比上面右手坐标系绕z轴逆时针旋转矩阵，实际推导的方式是一致的，只是这里用的行向量左乘的旋转矩阵，规定的正方向不同而已。

image

旋转的逆变换只需反方向旋转相同角度,任何旋转矩阵的逆都是其转置矩阵，因为旋转矩阵是一种正交矩阵。

image

正方向定义
  左手坐标系 左手坐标系 右手坐标系 右手坐标系
从哪里看 正旋转 负旋转 正旋转 负旋转
轴的负端，朝向轴的正端 逆时针 顺时针 顺时针 逆时针
轴的正端，朝向轴的负端 顺时针 逆时针 逆时针 顺时针

下列Unity Shader示例(左手法则，左手坐标系)中从左向右依次为,绕x轴顺时针旋转，绕y轴顺时针旋转，绕z轴顺时针旋转;观察方向均为轴正端向原点观测。

image

Shader "Unlit/Rotation"
{
	Properties
	{
	 _Switch("Switch",Vector) = (0,0,0,0)
	}
	SubShader
	{

		Pass
		{
			CGPROGRAM
			#pragma vertex vert
			#pragma fragment frag
			#include "UnityCG.cginc"

			struct appdata{float4 vertex : POSITION;};

			struct v2f{float4 vertex : SV_POSITION;};

			float4 _Switch;

			v2f vert (appdata v)
			{
				v2f o;

				float angle = _Time.y;

				if (_Switch.x == 1) {
					float y = v.vertex.y*cos(angle) - v.vertex.z*sin(angle);
					float z = v.vertex.y*sin(angle) + v.vertex.z*cos(angle);
					v.vertex.y = y;
					v.vertex.z = z;
				};

				if (_Switch.y == 1) {
					float x = v.vertex.x*cos(angle) + v.vertex.z*sin(angle);
					float z = -v.vertex.x*sin(angle) + v.vertex.z*cos(angle);
					v.vertex.x = x;
					v.vertex.z = z;
				};

				if (_Switch.z == 1) {
					float x = v.vertex.x*cos(angle) - v.vertex.y*sin(angle);
					float y = v.vertex.x*sin(angle) + v.vertex.y*cos(angle);
					v.vertex.x = x;
					v.vertex.y = y;
				};

				
				o.vertex = UnityObjectToClipPos(v.vertex);
				return o;
			}
			
			fixed4 frag (v2f i) : SV_Target{return fixed4(1,1,0,1);}
			ENDCG
		}
	}
}
绕任意轴的三维旋转

在右手坐标系下点H =(x，y，z)绕任意轴V旋转θ得到L =(x’，y’，z’)可通过以下步骤实现

(1)平移空间，使旋转轴穿过原点 (T)

(2)围绕y轴旋转空间，使旋转轴位于zy平面内 (Ry)

(3)绕x轴旋转空间，使旋转轴沿z轴 (Rx)

(4)通过θ绕z轴进行所需的旋转 (Rz)

(5)应用步骤(3)的逆 (Rx逆)

(6)应用步骤(2)的逆 (Ry逆)

(7)应用步骤(1)的逆 (T逆)

image

上图是2,3,4步骤

罗德里格(Rodrigues)旋转公式

向量V绕任意轴向量N顺时针旋转θ到向量V’,设‖N‖ = 1 image

上图可知：

向量Proj_v为向量V的在向量N方向上投影。

Proj_v = (N·V) * N / ‖N‖² 推导见前文点乘的几何意义：投影

V⊥ = V - Proj_v

image

‖N × V‖ = ‖N‖‖V‖sina = ‖v‖sina = ‖ V⊥ ‖

V⊥’ = cosθV⊥ + sinθ(N × V)

则：V’ = V⊥’ + Proj_v

V’ = cosθV⊥ + sinθ(N × V) + (N·V) * N

V’ = cosθ(V - (N·V) * N) + sinθ(N × V) + (N·V) * N

V’ = cosθV + (1-cosθ)(N·V)*N + sinθ(N × V)

设V为(x,y,z),N为(Nx,Ny,Nz) 带入得如下:

V’ = (Vx,Vy,Vz) cosθ + (NyVz - NzVy, NzVx - NxVz, NxVy - NyVx) sinθ + (Nx, Ny, Nz)(VxNx + VyNy + VzNz)(1 - cosθ)

V’.x = Vx cosθ + (NyVz - NzVy) sinθ + Nx (VxNx + VyNy + VzNz) ( 1- cosθ)

V’.y = Vy cosθ + (NzVx - NxVz) sinθ + Ny (VxNx + VyNy + VzNz) ( 1- cosθ)

V’.z = Vz cosθ + (NxVy - NyVx) sinθ + Nz (VxNx + VyNy + VzNz) ( 1- cosθ)

整理 V’.x = Vx(cosθ + NxNx(1 -cosθ)) + Vy(Nzsinθ + NxNy(1 -cosθ)) + Vz(Nysinθ + NxNz(1 -cosθ)) ，依次整理V’.y和V’.z 可得旋转矩阵表示如下：

image

当以x轴旋转N为=(1,0,0)，y轴旋转N为=(0,1,0)，z轴旋转N为=(0,0,1)可得上面推导的旋转矩阵

image

欧拉角

任何一个旋转可以表示为依次绕着三个旋转轴旋三个角度的组合。这三个角度称为欧拉角(pitch,heading,bank)。

三个轴可以指固定的世界坐标系轴，也可以指被旋转的物体坐标系的轴。三个旋转轴次序不同，会导致结果不同。

旋转矩阵转欧拉角

下图为欧拉旋转公式，采用上面右手坐标系中旋转矩阵(列向量右乘方式)，旋转顺序是zxy，计算顺序从左向右。

image

sinp = m32 ,角度p = asin(m32);

tan h = y/x ,tan h = sin h/cos h ,h = atan2(y，x)求的是y/x的反正切;

h = atan2(-m31/cosp ,m33/cosp) ,简化为角度h = atan2(-m31, m33)

b = atan2(-m12/cosp ,m22/cosp) ,简化为角度b = atan2(-m12, m22)

Unity atan2方法定义

采用上面左手坐标系中旋转矩阵(行向量左乘方式)，旋转顺序是zxy，计算顺序从左向右。

image

-sinp = m32 ,角度p = asin(-m32);

h = atan2(m31/cosp ,m33/cosp) ,简化为角度h = atan2(m31, m33)

b = atan2(m12/cosp ,m22/cosp) ,简化为角度b = atan2(m12, m22)

上面讨论都是当cosp为分母，cosp ！= 0的情况。如果p(X)=±90度时,cosp = 0，sinp=±1，那么m31,m33,m12,m22均为0,Y和Z共面。

对比如上二种计算方式，只是符号相反而已。

欧拉角分类

欧拉角按旋转的坐标系分为内旋(intrinsic rotation)和外旋(extrinsic rotation)。

  • 内旋(intrinsic rotation)- 动态: 绕物体自身的坐标系object-space 旋转，每一次旋转都会改变下一次旋转的轴，旋转的轴是动态。

  • 外旋(extrinsic rotation)- 静态:绕惯性坐标系upright-space 旋转(upright space指基向量平行于world-space或parent-space，原点与object-space的原点重合的空间)。 无论是三步旋转中的哪一步，轴都是固定的，是不会动的。

欧拉角顺规

按旋转轴分为经典欧拉角(Proper Euler Angle)和泰特布莱恩角(Tait–Bryan angles)

  • 经典欧拉角(Proper Euler Angle：按(z-x-z, x-y-x, y-z-y, z-y-z, x-z-x, y-x-y)轴序列旋转，即第一个旋转轴和最后一个旋转轴相同。

  • 泰特布莱恩角(Tait–Bryan angles)：按(x-y-z, y-z-x, z-x-y, x-z-y, z-y-x, y-x-z)轴序列旋转，即三个不同的轴

在unity中transform.Rotate(Vector3 eularAngles Space relativeTo) 如果relativeTo为Space.Self，采用yxz顺规（heading - pitch - bank），是指从惯性坐标系到物体坐标系，依据左手法则，+x向右，+y向上，+z向前

  1. 此时物体坐标系和惯性坐标系重合，heading为绕y轴的旋转量，向右旋转为正（如果从上面看，旋转正方向就是顺时针方向），
  2. 经过heading旋转后，pitch为绕x轴的旋转量，注意是物体坐标系的x轴，不是原惯性坐标系的x轴。依然遵守左手法则，向下旋转为正
  3. 最后，经过了heading和pitch，bank为绕z轴的旋转量。再次提醒，是物体坐标系的z轴，不是原惯性坐标系的z轴。依据左手法则，从原点向+z看，逆时针方向为正

由此可见，按yxz顺序旋转的话，每一步都是在上一步变化后的基础上进行旋转，即每一步都是相对于当前物体坐标系进行旋转，所以此系统是“yxz物体空间旋转系统”，即正是unity所采用的欧拉角系统。

image

Unity欧拉角(Transform.eulerAngles)使用zxy(roll - pitch - yaw)的顺规，是指向量从物体坐标系到惯性坐标系，和(heading - pitch - bank)顺序相反

image

  1. 首先绕z轴进行旋转，由于z轴是终端节点，所以x轴和y轴都不会发生变化。

  2. 再绕x轴进行旋转，由于x轴在上一步旋转中没有发生变化，所以就等价于绕惯性空间的x轴旋转。由于y轴是x轴父节点，所以y轴不会发生变化。

  3. 再绕y轴进行旋转，由于y轴在上一步旋转中没有发生变化，所以就等价于绕惯性空间的y轴旋转。

由此可见，按zxy顺序旋转的话，每一步都等价于相对于惯性坐标系进行旋转，所以此系统是“zxy惯性空间旋转系统”，即正是unity所采用的欧拉角系统。

万向锁(Gimble lock)问题

  • 使用动态欧拉角会出现万向锁现象；静态欧拉角不存在万向锁的问题。

  • 四元数不存在万向锁的问题。

  • 最简单的万向锁现象：一个模型先把X旋转改为90(蘑菇pitch俯仰低头)，模型的Z轴和世界的Y共面了，那么现在可以试着调整一下模型的Y旋转，表现为模型local坐标轴Z轴的旋转。如果你调整模型的Z旋转就会发现这Y和Z两个调整的是同一个轴的旋转(yaw和roll)，(90,90,0)和(90,0,-90)是同一个旋转结果。

image

万向锁是指物体的两个旋转轴指向同一个方向。实际上，当两个旋转轴平行时，我们就说万向节锁现象发生了，换句话说，绕一个轴旋转可能会覆盖住另一个轴的旋转，从而失去一维自由度。

四元数(Quaternion)

Quaternion = (xi + yj + zk + w ) = (x,y,z,w)

四元数是最简单的超复数。 复数是由实数加上元素 i 组成，其中i^2 = -1。 相似地，四元数都是由实数加上三个元素 i、j、k 组成，而且它们有如下的关系： i^2 = j^2 = k^2 = -1 ,ij=k、ji=-k、jk=i、kj=-i、ki=j、ik=-j, 每个四元数都是 1、i、j 和 k 的线性组合，即是四元数一般可表示为w + xi + yj + zk，其中x、y、z 、w是实数

在unity中表示


namespace UnityEngine
{
	public struct Quaternion
	{
		public float x;
		public float y;
		public float z;
		public float w;
  }
}

四元数的加法

定义二个四元数，由标量a与向量v构成。

image

跟复数、向量和矩阵一样，两个四元数之和需要将不同的元素加起来。

image

四元数的乘法

两个四元数之间的非可换乘积通常被格拉斯曼(Hermann Grassmann)称为积

image

四元数乘法的非可换性，$q_{1}q_{2}$并不等于$q_{2}q_{1}$

四元数乘法可观察出向量点积 “·” 和叉积”×” 。

image

四元数的模

image

模为1.称为单位四元数

四元数的共轭和逆

四元数和它的共轭(conjugate)q*代表相反的角位移，因为相当于旋转轴反向。四元数的逆$q^{-1}$。

image

如果$\vert q \vert = 1$,q视为单位四元数,则$q^{-1} = q$，所以$q^{-1}$和$q$是相同的旋转。

四元数描述旋转

如图所示，u为旋转轴，旋转角度为σ，向量v旋转到w处.

image

image

如果我们取每个术语的平方根，我们得到四元数的部分

image

四元数转旋转矩阵

v旋转后的坐标v’为：v’= qvq^-1,当q为单位四元数则v’= qvq*,需求出一个矩阵M,有v’=Mv,假设v=0+xi+yi+zk,q=w+qxi+qyj+qzk,带入式有：

image

全部展开后如下式：

image

将x,y和z项分组并将它们放在矩阵中给出：

image

Unity中四元数转矩阵的c++实现如下：


//Quaternion.cpp


// m0,0 m0,1 m0,2 m0,3

// m1,0 m1,1 m1,2 m1,3

// m2,0 m2,1 m2,2 m2,3

// m3,0 m3,1 m3,2 m3,3


// The floats are laid out:

// m0   m4   m8	  m12	

// m1   m5   m9   m13		

// m2   m6   m10  m14		

// m3   m7   m11  m15		

void QuaternionToMatrix (const Quaternionf& q, Matrix4x4f& m)
{
	// If q is guaranteed to be a unit quaternion, s will always
	// be 1.  In that case, this calculation can be optimized out.

	#if DEBUGMODE
	if (!CompareApproximately (SqrMagnitude (q), 1.0F, Vector3f::epsilon))
	{
		AssertString(Format("Quaternion To Matrix conversion failed because input Quaternion is invalid %f, %f, %f, %f l=%f", q.x, q.y, q.z, q.w, SqrMagnitude(q)));		
	}
	#endif
 
	//float norm = GetNorm (q);
	//float s = (norm > 0.0) ? 2.0/norm : 0;
 
	// Precalculate coordinate products

	float x = q.x * 2.0F;
	float y = q.y * 2.0F;
	float z = q.z * 2.0F;
	float xx = q.x * x;
	float yy = q.y * y;
	float zz = q.z * z;
	float xy = q.x * y;
	float xz = q.x * z;
	float yz = q.y * z;
	float wx = q.w * x;
	float wy = q.w * y;
	float wz = q.w * z;
 
	// Calculate 3x3 matrix from orthonormal basis

	m.m_Data[0] = 1.0f - (yy + zz);
	m.m_Data[1] = xy + wz;
	m.m_Data[2] = xz - wy;
	m.m_Data[3] = 0.0F;
 
	m.m_Data[4] = xy - wz;
	m.m_Data[5] = 1.0f - (xx + zz);
	m.m_Data[6] = yz + wx;
	m.m_Data[7] = 0.0F;
 
	m.m_Data[8]  = xz + wy;
	m.m_Data[9]  = yz - wx;
	m.m_Data[10] = 1.0f - (xx + yy);
	m.m_Data[11] = 0.0F;
 
	m.m_Data[12] = 0.0F;
	m.m_Data[13] = 0.0F;
	m.m_Data[14] = 0.0F;
	m.m_Data[15] = 1.0F;
}

旋转矩阵转四元数

image

image

在线性代数中，一个n×n矩阵M的主对角线（从左上方至右下方的对角线）上各个元素的总和被称为矩阵A的迹（或迹数），一般记作tr(M)。

image

再次观察该矩阵有：

image

四元数转欧拉角

下面旋转矩阵使用泰特布莱恩角(Tait–Bryan angles)的yxz顺规,矩阵代表惯性到物体的旋转.

image

分析如上二式有

image

在unity中c++实现


//Quaternion.cpp

Vector3f QuaternionToEuler (const Quaternionf& quat)
{
	Matrix3x3f m;
	Vector3f rot;
	QuaternionToMatrix (quat, m);
	MatrixToEuler (m, rot);
	return rot;
}


void QuaternionToMatrix (const Quaternionf& q, Matrix3x3f& m)
{
	// If q is guaranteed to be a unit quaternion, s will always
	// be 1.  In that case, this calculation can be optimized out.
	#if DEBUGMODE
	if (!CompareApproximately (SqrMagnitude (q), 1.0F, Vector3f::epsilon))
	{
		AssertString(Format("Quaternion To Matrix conversion failed because input Quaternion is invalid %f, %f, %f, %f l=%f", q.x, q.y, q.z, q.w, SqrMagnitude(q)));		
	}
	#endif
	//float norm = GetNorm (q);
	//float s = (norm > 0.0) ? 2.0/norm : 0;

	// Precalculate coordinate products
	float x = q.x * 2.0F;
	float y = q.y * 2.0F;
	float z = q.z * 2.0F;
	float xx = q.x * x;
	float yy = q.y * y;
	float zz = q.z * z;
	float xy = q.x * y;
	float xz = q.x * z;
	float yz = q.y * z;
	float wx = q.w * x;
	float wy = q.w * y;
	float wz = q.w * z;

	// Calculate 3x3 matrix from orthonormal basis
	m.m_Data[0] = 1.0f - (yy + zz);
	m.m_Data[1] = xy + wz;
	m.m_Data[2] = xz - wy;

	m.m_Data[3] = xy - wz;
	m.m_Data[4] = 1.0f - (xx + zz);
	m.m_Data[5] = yz + wx;
	
	m.m_Data[6]  = xz + wy;
	m.m_Data[7]  = yz - wx;
	m.m_Data[8] = 1.0f - (xx + yy);
}

///Matrix3x3.cpp

// m0,0 m0,1 m0,2

// m1,0 m1,1 m1,2

// m2,0 m2,1 m2,2


// The floats are laid out:

// m0   m3   m6	

// m1   m4   m7	

// m2   m5   m8	

const float& Get (int row, int column)const 	{ return m_Data[row + (column * 3)]; }

/// This is YXZ euler conversion

bool MatrixToEuler (const Matrix3x3f& matrix, Vector3f& v)
{
	// from http://www.geometrictools.com/Documentation/EulerAngles.pdf
	// YXZ order
	if ( matrix.Get(1,2) < 0.999F ) // some fudge for imprecision
	{
		if ( matrix.Get(1,2) > -0.999F ) // some fudge for imprecision
		{
			v.x = asin(-matrix.Get(1,2));
			v.y = atan2(matrix.Get(0,2), matrix.Get(2,2));
			v.z = atan2(matrix.Get(1,0), matrix.Get(1,1));
			SanitizeEuler (v);
            return true;
        }
        else
        {
            // WARNING.  Not unique.  YA - ZA = atan2(r01,r00)
            v.x = kPI * 0.5F;
            v.y = atan2(matrix.Get (0,1), matrix.Get(0,0));
            v.z = 0.0F;
			SanitizeEuler (v);
            
            return false;
        }
    }
    else
    {
        // WARNING.  Not unique.  YA + ZA = atan2(-r01,r00)
        v.x = -kPI * 0.5F;
        v.y = atan2(-matrix.Get(0,1),matrix.Get(0,0));
        v.z = 0.0F;
 		SanitizeEuler (v);
        return false;
    }
}

inline void SanitizeEuler (Vector3f& euler)
{
	MakePositive (euler);
}

inline void MakePositive (Vector3f& euler)
{
	const float negativeFlip = -0.0001F;
	const float positiveFlip = (kPI * 2.0F) - 0.0001F;
	
	if (euler.x < negativeFlip)
		euler.x += 2.0 * kPI;
	else if (euler.x > positiveFlip)
		euler.x -= 2.0 * kPI;
    
	if (euler.y < negativeFlip)
		euler.y += 2.0 * kPI;
	else if (euler.y > positiveFlip)
		euler.y -= 2.0 * kPI;
    
	if (euler.z < negativeFlip)
		euler.z += 2.0 * kPI;
	else if (euler.z > positiveFlip)
		euler.z -= 2.0 * kPI;
}

欧拉角转四元数

image


inline friend Quaternionf operator * (const Quaternionf& lhs, const Quaternionf& rhs)
{
	return Quaternionf (
				lhs.w*rhs.x + lhs.x*rhs.w + lhs.y*rhs.z - lhs.z*rhs.y,
				lhs.w*rhs.y + lhs.y*rhs.w + lhs.z*rhs.x - lhs.x*rhs.z,
				lhs.w*rhs.z + lhs.z*rhs.w + lhs.x*rhs.y - lhs.y*rhs.x,
				lhs.w*rhs.w - lhs.x*rhs.x - lhs.y*rhs.y - lhs.z*rhs.z);
}

image

在unity中c++实现


//Quaternion.cpp

Quaternionf EulerToQuaternion (const Vector3f& someEulerAngles)
{
	float cX (cos (someEulerAngles.x / 2.0f));
	float sX (sin (someEulerAngles.x / 2.0f));

	float cY (cos (someEulerAngles.y / 2.0f));
	float sY (sin (someEulerAngles.y / 2.0f));

	float cZ (cos (someEulerAngles.z / 2.0f));
	float sZ (sin (someEulerAngles.z / 2.0f));
	
	Quaternionf qX (sX, 0.0F, 0.0F, cX);
	Quaternionf qY (0.0F, sY, 0.0F, cY);
	Quaternionf qZ (0.0F, 0.0F, sZ, cZ);
	
	Quaternionf q = (qY * qX) * qZ;
	AssertIf (!CompareApproximately (SqrMagnitude (q), 1.0F));
	return q;
}

inline friend Quaternionf operator * (const Quaternionf& lhs, const Quaternionf& rhs)
{
	return Quaternionf (
				lhs.w*rhs.x + lhs.x*rhs.w + lhs.y*rhs.z - lhs.z*rhs.y,
				lhs.w*rhs.y + lhs.y*rhs.w + lhs.z*rhs.x - lhs.x*rhs.z,
				lhs.w*rhs.z + lhs.z*rhs.w + lhs.x*rhs.y - lhs.y*rhs.x,
				lhs.w*rhs.w - lhs.x*rhs.x - lhs.y*rhs.y - lhs.z*rhs.z);
}

Unity四元数

四元数可提供平滑差值，没有Euler旋转的万向锁。


var rotation = Quaternion.Euler(0, 30, 0);

//返回一个旋转角度，绕z轴旋转z度，绕x轴旋转x度，绕y轴旋转y=30度(像这样的顺序)


Quaternion q3 = new Quaternion();

q3.eulerAngles = new Vector3(10, 30, 20);

Quaternion qx3 = Quaternion.AngleAxis(10,Vector3.right);

//绕y轴旋转30度

Quaternion qy3 = Quaternion.AngleAxis(30,Vector3.up);

Quaternion qz3 = Quaternion.AngleAxis(20,Vector3.forward);
		
Quaternion qxyz3 = qz3*qy3*qx3;

//q3和qxyz3值一样

Unity复合变换

有了以上平移、缩放、旋转矩阵后，我们就可以通过矩阵乘法求得点P任意变化后坐标P’，注意下图是P是列矩阵，用的右乘。

image

矩阵乘法不满足交换律，先后乘的顺序会导致结果是不一致的。Unity里约定变换 1.先缩放2.再旋转3.最后平移

为什么这样约定?如果在空间放一个cube，先移动(10,20),再自身旋转60度，最后放大5五倍。你会发现cube绕自己本身原点以外位置的原点的旋转和缩放，这并不我们想要的效果。而先放大5五倍，再自身旋转60度，最后移动(10,20)比较符合使用习惯。


//c#

public static Matrix4x4 TRS(Vector3 pos, Quaternion q, Vector3 s)
{
	return Matrix4x4.INTERNAL_CALL_TRS(ref pos, ref q, ref s);
}

//c++

void Matrix4x4f::SetTRS (const Vector3f& pos, const Quaternionf& q, const Vector3f& s)
{
	QuaternionToMatrix (q, *this);
 
	m_Data[0] *= s[0];
	m_Data[1] *= s[0];
	m_Data[2] *= s[0];
 
	m_Data[4] *= s[1];
	m_Data[5] *= s[1];
	m_Data[6] *= s[1];
 
	m_Data[8] *= s[2];
	m_Data[9] *= s[2];
	m_Data[10] *= s[2];
 
	m_Data[12] = pos[0];
	m_Data[13] = pos[1];
	m_Data[14] = pos[2];

}

Unity旋转转换

image

旋转对比

各方法比较 Matrix4x4 Euler Quaternion
在坐标系间(物体和惯性)旋转点 能 不能(必须转换到 矩阵)   不能(必须转换到矩阵) 
连接或增量旋转 能,但经常比四元数慢,小心矩阵蠕变的情况 不能 能,比矩阵快
插值 基本上不能   能,但可能遭遇万 向锁或其他问题  Slerp提供了平滑插值
易用程度 难 易 难
在内存或文件中存储 9个数 3个数 4个数
对给定方位的表达方式是否唯一 是 不是,对同一方位 有无数多种方法 不是,有两种方法,它们相互为负
可能导致非法 矩阵蠕变 任意三个数都能构成合法的欧拉角 可能会出现误差积累,从而产生非法的四元数

渲染管线

image

在unity shader中一个模型顶点到裁剪空间，通过mul(UNITY_MATRIX_MVP, v.vertex)就可计算。分解计算步骤如下：

1.模型到世界用mul(UNITY_MATRIX_M, v.vertex) ,c#中矩阵为Matrix4x4 world = transform.localToWorldMatrix;

2.世界到观察用mul(UNITY_MATRIX_V, v.vertex),c#中矩阵为Matrix4x4 view = camera.worldToCameraMatrix;

3.观察到裁剪用mul(UNITY_MATRIX_P, v.vertex),c#中矩阵为Matrix4x4 proj = camera.projectionMatrix;

4.裁剪到屏幕是Unity系统完成的。

如果直接就可以用UnityObjectToClipPos(v.vertex)来处理效率更高一些，

到此角色描边顶点变换理论基础总结完毕。

参考资料:

Youtube旋转矩阵推导

绕任意轴旋转

万向锁

四元数

四元数转化

罗德里格旋转公式

四元数和欧拉角转换

欧拉角

3D Rotation Converter计算器

Nasa推导

CATALOG
    FEATURED TAGS
    Lua Unity HardWare C++ C# Unreal
    FRIENDS