阅读记录How to Create a Custom 2D Physics Engine - Rangy Gaul（1）
本文为物理相关文章的阅读记录，包括了翻译、个人勘误以及一些个人理解，原文链接https://gamedevelopment.tutsplus.com/series/how-to-create-a-custom-physics-engine--gamedev-12715
A custom physics engine can tackle any sort of technical effect the creator has the skill to create.自定义物理引擎可以处理创作者有能力创造的任何技术效果。
The Basics and Impulse Resolution
Axis Aligned Bounding Boxes
基础讲了AABB,拓展有AABB的intersecting检测
AABB盒为复杂几何体做用于一个简单测试
AABB的数据结构：
struct AABB
{
Vec2 min;
Vec2 max;
};
In order to tell whether two AABB shapes are intersecting you will need to have a basic understanding of the Separating Axis Theorem (SAT).
Here's a quick test taken from Real-Time Collision Detection by Christer Ericson, which makes use of the SAT:实时碰撞检测中AABB测试的SAT的算法大致如下
bool AABBvsAABB( AABB a, AABB b )
{
// Exit with no intersection if found separated along an axis
if(a.max.x < b.min.x or a.min.x > b.max.x) return false
if(a.max.y < b.min.y or a.min.y > b.max.y) return false
// No separating axis found, therefor there is at least one overlapping axis
return true
}
Circle
数据结构
struct Circle
{
float radius
Vec position
};
测试两个圆是否相交很简单:取两个圆的半径并将它们相加，然后检查这个和是否大于两个圆之间的距离。
计算距离时一般使用平方而不开根号。
冲量方法
是一种碰撞解决方法
这里给出的是2维的例子
一般来说，物理引擎中的物体有三个主要的自由度(在二维空间中):在xy平面上的运动（2个）和垂直平面轴的旋转（1个）。在本文中，我们隐式地限制旋转，只使用aabb和圆，所以我们真正需要考虑的唯一自由度是沿着xy平面的运动。
碰撞处理的目的是我们要让这两个参与碰撞的物体不再Intersecting。
冲量方法的核心就是使用一个速度突变值来分离物体。为了做到这一点，必须以某种方式考虑每个物体的质量、位置和速度:我们希望大的物体与小的物体碰撞时，在碰撞过程中移动一点，并让小的物体飞离。我们还希望质量无穷大的物体完全不运动。
我们将使用刚体和一些数学知识。刚体的定义是受到外力后不可变形的物体。如果可变形我们物理计算公式的复杂度将提升一个档次。
碰撞检测得到碰撞法线+渗透深度
为了给两个物体施加脉冲使它们分开，我们需要知道推动它们的方向和力度。碰撞法线是冲击作用的方向。穿透深度(以及其他一些东西)决定了脉冲的大小。这意味着唯一需要解的值就是脉冲的大小。
Equation 1
从A指向B的向量，我们计算时的公式是：末位点 减去 开始点
n用于表示碰撞的法向量
Equation 2
Equation 3 点乘式子
下一步是引入所谓的恢复系数e。恢复是指弹性或弹性，可能中学物理学的弹性碰撞的那个系数就是指这个东西。你应该总是使用碰撞中最低恢复系数。
// Given two objects A and B
e = min( A.restitution, B.restitution )
结合牛顿恢复定律，用恢复系数 计算得到一个碰撞后的速度
Equation 4
这就是说，碰撞后的速度等于碰撞前的速度，乘以某个常数。这个常数表示“反弹因子”。
整合到Equation 5(和原文有改动，我加了 ' )
在牛顿恢复定律中，反弹后的矢量，实际上是指向V的反方向，那么在推导过程中，我们如何表示反方向呢?引入一个负号。
现在我们需要能够表达在脉冲影响下的速度。下面是一个简单的方程，用于沿特定方向用脉冲标量j沿着特定方向n修改速度矢量
Equation 6
其实就是把一个想来你加到另一个向量上，j其实就表示这个冲量的大小
在形式上，冲量被定义为动量的变化。动量是质量*速度。
Equation 7
当物体a和物体b发生碰撞时，a会被推向b相反的方向:
Equation 8 (这里的正负号也有修改，我是觉得碰撞A的话用减更好，符号和碰撞向量的定义有关，取正负我觉得都可以，这里取负主要是为了和后面的推导过程不冲突)
现在需要做的就是合并方程8和5，下式减去上式然后左边代入公式5，再左右移项。我们得到的方程是这样的:（这里也和原文有改动，估计是原文公式写错了）
Equation 9
这里也有改动，j应该是提出来了
Equation 10
n是方向向量，单位向量，自乘就等于1
这意味着我们可以编写几行代码来求解脉冲标量(j)。代码比数学符号可读性强多了!
void ResolveCollision( Object A, Object B )
{
// Calculate relative velocity
Vec2 rv = B.velocity - A.velocity
// Calculate relative velocity in terms of the normal direction
float velAlongNormal = DotProduct( rv, normal )
// Do not resolve if velocities are separating
if(velAlongNormal > 0)
return;
// Calculate restitution
float e = min( A.restitution, B.restitution)
// Calculate impulse scalar
float j = -(1 + e) * velAlongNormal
j /= 1 / A.mass + 1 / B.mass
// Apply impulse
Vec2 impulse = j * normal
A.velocity -= 1 / A.mass * impulse
B.velocity += 1 / B.mass * impulse
}
在上面的代码示例中，有几件关键的事情需要注意。首先是第10行检查，if(VelAlongNormal > 0).这个检查非常重要;它确保只在物体相互移动时才解决碰撞。
通常会记住质量的导数
A.inv_mass = 1 / A.mass
许多物理引擎实际上并不存储原始质量。物理引擎通常只存储质量倒数。碰巧大多数涉及质量的数学都是使用质量倒数。
最后要注意的是，我们智能地将脉冲标量(j)分配到两个对象上。我们希望小的物体的（j）比较大，大的物体的速度只被一小部分的(j)所改变。
float mass_sum = A.mass + B.mass
float ratio = A.mass / mass_sum
A.velocity -= ratio * impulse
ratio = B.mass / mass_sum
B.velocity += ratio * impulse
Sinking Object问题
如果我们继续使用到目前为止的代码，对象将会相互碰撞并反弹。这很好，不过如果一个物体的质量是无穷大呢?我们需要一个好方法来表示无限质量在我们的模拟中。
我建议用零作为无限质量——尽管如果我们试图计算一个物体的质量为零的逆，我们将会被零除法。解决这个问题的方法是在计算逆质量时做下面的事情:
if(A.mass == 0)
A.inv_mass = 0
else
A.inv_mass = 1 / A.mass
简单的说就是设置特殊值
这还是可以的。Sinking Objects问题出现在物体由于重力开始下沉到另一个物体时。也许是某种恢复力很低的东西撞上了质量无限的墙，然后开始下沉。
这种下沉是由于浮点错误。在每次浮点运算过程中，由于硬件的原因，都会引入一个小的浮点误差。随着时间的推移，这个错误累积在位置错误中，导致对象陷入另一个对象中。
为了纠正这个错误，必须对它进行解释。为了纠正这个位置误差，我将向你展示一种叫做线性投影的方法。线性投影使两个物体的穿透率降低了一个小百分比，这是在施加脉冲后进行的。位置校正非常简单:按穿透深度的百分比沿碰撞法(n)移动每个物体:
void PositionalCorrection( Object A, Object B )
{
const float percent = 0.2 // usually 20% to 80%
Vec2 correction = penetrationDepth / (A.inv_mass + B.inv_mass)) * percent * n
A.position -= A.inv_mass * correction
B.position += B.inv_mass * correction
这就是位置级别的纠正，有类似PBD的思想
请注意，我们用系统的总质量来衡量穿透深度。这将给出一个与我们所处理的质量成比例的位置修正。小的物体比重的物体推得快。
这个实现有一个小问题:如果我们总是解决我们的位置错误，那么对象将在相互依赖时来回抖动。为了防止这种情况的发生，必须采取一些宽松措施。我们只在穿透超过某个被称为“slop”的任意阈值时执行位置校正:
void PositionalCorrection( Object A, Object B )
{
const float percent = 0.2 // usually 20% to 80%
const float slop = 0.01 // usually 0.01 to 0.1
Vec2 correction = max( penetration - k_slop, 0.0f ) / (A.inv_mass + B.inv_mass)) * percent * n
A.position -= A.inv_mass * correction
B.position += B.inv_mass * correction
}
这使得物体可以在没有位置校正的情况下略微穿透。
简单流型生成
本文要讨论的最后一个主题是简单流形生成。在数学术语中，流形是“在空间中代表一个区域的点的集合”。然而，当我提到术语流形时，我指的是一个包含关于两个对象之间碰撞信息的小对象。
其实理解成碰撞信息就好了，在冲量方法中就是渗透深度和法向量
struct Manifold
{
Object *A;
Object *B;
float penetration;
Vec2 normal;
};
在碰撞检测过程中，既要计算穿透度，也要计算碰撞法线。为了找到这些信息，必须对本文开头的原始碰撞检测算法进行扩展。这里比较长就建议看原文了。
Core Engine
积分
F = ma
Erin Catto's Integration Demo from GDC 2009 and Hannu's addition to symplectic Euler for more stability in low FPS environments.
// Explicit Euler
x += v * dt
v += (1/m * F) * dt
dt表示delta time
// Symplectic Euler
v += (1/m * F) * dt
x += v * dt
为何需要定步长？
我们不希望游戏在不同电脑上的物理的运行速度和电脑的性能有关
我们希望我们物理引擎每次只按定时间步长运作
Using the exact same dt value in your code everywhere will actually make your physics engine deterministic, and is known as a fixed timestep
在你的代码中使用完全相同的dt值将使你的物理引擎具有确定性，这被称为固定的时间步长。使用固定时间步长的其他理由：verlet积分的使用对步长有要求
在守望先锋中服务器和客户端为了两端模拟的确定性，也采用了定时间步长发包的方式。
假设给出了相同的输入，确定性物理引擎总是在每次运行时都做同样的事情。
const float fps = 100
const float dt = 1 / fps
float accumulator = 0
// In units of seconds
float frameStart = GetCurrentTime( )
// main loop
while(true)
const float currentTime = GetCurrentTime( )
// Store the time elapsed since the last frame began
accumulator += currentTime - frameStart( )
// Record the starting of this frame
frameStart = currentTime
while(accumulator > dt)
UpdatePhysics( dt )
accumulator -= dt
RenderGame( )
记录过去的时间，当其满足了一次物理迭代时就进行物理引擎更新，或许我们应该把这个放到别的线程，主线程和物理线程用同步的方法去做
Chunks of dt are removed from the accumulator until the accumulator is smaller than a dt chunk.
这个循环是有问题的，物理引擎的更新也作为当前帧的时间计算了，如果物理算的太久，就一直卡物理了，这就是死亡螺旋
为了解决这个问题，如果累加器变得过高，引擎就需要运行更少的物理更新。一个简单的方法是将累加器固定在某个任意值以下。
// Avoid spiral of death and clamp dt, thus clamping
// how many times the UpdatePhysics can be called in
// a single game loop.
if(accumulator > 0.2f)
accumulator = 0.2f
物理机制不会死亡，但这样还是会运行的慢一些，而且这样还会有物理跳跃现象
const float fps = 100
const float dt = 1 / fps
float accumulator = 0
// In units seconds
float frameStart = GetCurrentTime( )
// main loop
while(true)
const float currentTime = GetCurrentTime( )
// Store the time elapsed since the last frame began
accumulator += currentTime - frameStart( )
// Record the starting of this frame
frameStart = currentTime
// Avoid spiral of death and clamp dt, thus clamping
// how many times the UpdatePhysics can be called in
// a single game loop.
if(accumulator > 0.2f)
accumulator = 0.2f
while(accumulator > dt)
UpdatePhysics( dt )
accumulator -= dt
const float alpha = accumulator / dt;
RenderGame( alpha )
void RenderGame( float alpha )
for shape in game do
// calculate an interpolated transform for rendering
Transform i = shape.previous * alpha + shape.current * (1.0f - alpha)
shape.previous = shape.current
shape.Render( i )
令渲染落后物理一帧，然后插值平滑所有的累加器余量。这样的实际表现就会显得很平滑
为什么不去让渲染快一点呢，去插值到当前物理的位置+时间累加器预测的位置？
1、这样需要预测
2、物理引擎中突变的情况很多，对于这种突变的预测产生的错误预测往往会导致一些瞬移
模块设计
Bodies 物理主体是一个包含了关于某个给定物理对象的所有信息的对象。它将存储对象所代表的形状、大量数据、变换(位置、旋转)、速度、扭矩等。
struct body
{
Shape *shape;
Transform tx;
Material material;
MassData mass_data;
Vec2 velocity;
Vec2 force;
real gravityScale;
};
body应该有多种Shape，以组合的方式提供，可以添加多种shape
struct MassData
{
float mass;
float inv_mass;
// For rotations (not covered in this article)
float inertia;
float inverse_inertia;
};
对旋转就有转动惯量
力，矢量，力在每次迭代计算，每次迭代都不同，算完之后可以归零
宽相位计算
我们使用宽相位碰撞检测来确定哪些对物体可能发生碰撞，然后用窄相位碰撞检测来检查它们是否真的发生碰撞
struct Pair
{
body *A;
body *B;
};
广义相应该收集一系列可能的碰撞，并将它们全部存储在Pair结构中。这些对可以传递到引擎的另一部分(狭窄阶段)，然后解决。
这是一个O（n^2）的方法，使用空间数据结构可以将它降到O（n）
// Generates the pair list.
// All previous pairs are cleared when this function is called.
void BroadPhase::GeneratePairs( void )
{
pairs.clear( )
// Cache space for AABBs to be used in computation
// of each shape's bounding box
AABB A_aabb
AABB B_aabb
for(i = bodies.begin( ); i != bodies.end( ); i = i->next)
{
for(j = bodies.begin( ); j != bodies.end( ); j = j->next)
{
Body *A = &i->GetData( )
Body *B = &j->GetData( )
// Skip check with self
if(A == B)
continue
A->ComputeAABB( &A_aabb )
B->ComputeAABB( &B_aabb )
if(AABBtoAABB( A_aabb, B_aabb ))
pairs.push_back( A, B )
}
}
}
这里也会有很多重复的对子，我们需要将他们排序，删除重复的
// Sort pairs to expose duplicates
sort( pairs, pairs.end( ), SortPairs );
// Queue manifolds for solving
{
int i = 0;
while(i < pairs.size( ))
{
Pair *pair = pairs.begin( ) + i;
uniquePairs.push_front( pair );
++i;
// Skip duplicate pairs by iterating i until we find a unique pair
while(i < pairs.size( ))
{
Pair *potential_dup = pairs + i;
if(pair->A != potential_dup->B || pair->B != potential_dup->A)
break;
++i;
}
}
}
bool SortPairs( Pair lhs, Pair rhs )
{
if(lhs.A < rhs.A)
return true;
if(lhs.A == rhs.A)
return lhs.B < rhs.B;
return false;
}
术语lhs和rhs可以理解为“左手边”和“右手边”。这些术语通常用来指函数的参数，在这些参数中，事物在逻辑上可以被视为某个方程或算法的左右两边。
分层
不同层的body不会碰撞
// Generates the pair list.
// All previous pairs are cleared when this function is called.
void BroadPhase::GeneratePairs( void )
{
pairs.clear( )
// Cache space for AABBs to be used in computation
// of each shape's bounding box
AABB A_aabb
AABB B_aabb
for(i = bodies.begin( ); i != bodies.end( ); i = i->next)
{
for(j = bodies.begin( ); j != bodies.end( ); j = j->next)
{
Body *A = &i->GetData( )
Body *B = &j->GetData( )
// Skip check with self
if(A == B)
continue
// Only matching layers will be considered
if(!(A->layers & B->layers))
continue;
A->ComputeAABB( &A_aabb )
B->ComputeAABB( &B_aabb )
if(AABBtoAABB( A_aabb, B_aabb ))
pairs.push_back( A, B )
}
}
}
Halfspace Intersection
二维中，半空间表示一条直线的一面
判断在那一边，把点代入线的方程，然后看符号，结果为0表示该点在直线上，正/负表示直线的不同边。
法向量不一定是标准化的（也就是说，它的长度不一定是 1）。
2021.8.25
