ที่ผ่านมาเราพูดถึงแค่วัตถุชิ้นเดียวที่เคลื่อนที่ในพื้นที่: ตำแหน่งของมัน การหมุนของมัน เวกเตอร์ที่ใช้อธิบายมัน บทนี้จะพูดถึงสิ่งที่เกิดขึ้นเมื่อวัตถุพวกนี้หลายชิ้นอยู่ในพื้นที่เดียวกันแล้วเริ่มชนกัน — ลูกบอลตกลงพื้น กล่องไถลไปชนกำแพง กองกล่องที่วางซ้อนกันแล้วไม่ทะลุกันไปเอง นี่คือ rigid body physics และมันตั้งอยู่บนสองอย่างที่เรารู้จักอยู่แล้วทั้งหมด คือ semi-implicit Euler stepping และ fixed timestep จากบทแคลคูลัส กับเวกเตอร์ dot product และ AABB (axis-aligned bounding box) จากบทเรขาคณิต จริง ๆ แล้วกล่อง tip ในบทเรขาคณิตเรื่องวงกลมกับกล่องสี่เหลี่ยมเคยพูดคำว่า "broad phase" กับ "narrow phase" ไปแล้วโดยไม่ได้อธิบาย บทนี้คือคำสัญญานั้นที่จะจ่ายคืนให้ครบ
ทุกไอเดียด้านล่างนี้จะเดินตามรูปแบบเดิม: โปรแกรม C++ เล็ก ๆ ที่รันได้จริง เอาต์พุตจริงที่มันพิมพ์ออกมา แล้วค่อยอธิบายง่าย ๆ ว่าเกิดอะไรขึ้น หัวข้อ 1 ถึง 11 จะสร้าง physics engine เล็ก ๆ ขึ้นมาเองตั้งแต่ต้นด้วย C++ โดยใช้ struct Vec2 กับฟังก์ชัน add/sub/scale/dot/length จากบทเวกเตอร์ซ้ำ เพื่อให้เห็นชัด ๆ ว่าข้างใต้มันทำงานยังไง หัวข้อ 12 จะสลับไปใช้ C# แล้วโชว์ว่า component Rigidbody กับ Collider ของ Unity ทำงานแบบเดียวกันนี้ให้เราโดยอัตโนมัติยังไง เพื่อให้เรามองออกว่าแต่ละคอนเซปต์ที่เรียนไปอยู่ตรงไหนในเครื่องมือที่จะใช้จริงทุกวัน
Rigid body คือวัตถุที่รูปร่างไม่เปลี่ยนแปลงเลย — ระยะห่างระหว่างจุดสองจุดใด ๆ บนตัวมันจะคงที่เสมอ ไม่ว่ามันจะเคลื่อนที่หรือหมุนยังไงก็ตาม กล่อง ลูกบอล ตัวถังรถ ล้วนเป็น rigid body ส่วนผ้า เชือก หรือผิวหนังที่ยุบตัวได้ ไม่ใช่ — พวกนั้นเรียกว่า soft body ซึ่งเป็นหัวข้อที่ต่างออกไปและยากกว่า บทนี้จะไม่พูดถึง การจำกัดตัวเองให้อยู่กับรูปร่างที่ไม่เปลี่ยนรูปเลยนี่แหละ ที่ทำให้ rigid body physics คำนวณได้จริงในทางปฏิบัติ: แทนที่จะต้องติดตามทุกจุดบนวัตถุแยกกัน เราแค่ต้องติดตามตัวเลขไม่กี่ตัวสำหรับวัตถุทั้งชิ้นในคราวเดียว
ตัวเลขพวกนี้คือ state ของ rigid body:
นี่คือ state นั้นในรูปแบบ struct ของ C++ ซึ่งเป็นฐานที่ทุกหัวข้อถัดไปจะต่อยอดขึ้นมา:
#include <iostream>
struct Vec2 { float x, y; };
struct RigidBody {
Vec2 position;
float rotation; // 2D: a single angle, in radians
Vec2 velocity; // linear velocity, units per second
float angularVelocity; // 2D: radians per second
float mass;
float invMass; // 1 / mass; 0 means "infinite mass", i.e. never moves
};
RigidBody makeBody(Vec2 pos, float mass) {
RigidBody b;
b.position = pos;
b.rotation = 0.0f;
b.velocity = {0.0f, 0.0f};
b.angularVelocity = 0.0f;
b.mass = mass;
b.invMass = (mass > 0.0f) ? 1.0f / mass : 0.0f;
return b;
}
int main() {
RigidBody box = makeBody({0.0f, 5.0f}, 2.0f);
RigidBody floor = makeBody({0.0f, 0.0f}, 0.0f); // mass 0 = static, never moves
std::cout << "box: mass=" << box.mass << " invMass=" << box.invMass << "\n";
std::cout << "floor: mass=" << floor.mass << " invMass=" << floor.invMass << "\n";
}ผลลัพธ์:
box: mass=2 invMass=0.5
floor: mass=0 invMass=0สังเกตว่าเราเก็บ invMass (inverse mass หรือ 1/mass) แทนที่จะเก็บแค่ mass เฉย ๆ เพราะหัวข้อถัดไปทุกหัวข้อจะต้องหารด้วย "วัตถุสองชิ้นนี้รวมกันต้านการเคลื่อนที่มากแค่ไหน" — และกำแพงหรือพื้นไม่ควรขยับเลยไม่ว่าจะโดนกระแทกแรงแค่ไหน การเก็บ "mass เป็นอนันต์" ตรง ๆ ทำไม่ได้ (หารด้วยอนันต์ไม่ได้) แต่การเก็บ "inverse mass เป็นศูนย์" ทำได้ง่าย ๆ แถมคูณอะไรกับศูนย์ก็ได้ศูนย์แบบไม่มีปัญหา เคล็ดลับตัวนี้แหละที่ทำให้ invMass ไม่ใช่ mass โผล่มาในทุกสูตรตลอดบทนี้ วัตถุที่ invMass = 0 เรียกว่า static ส่วนที่เหลือเรียกว่า dynamic
rotation/angularVelocity เป็นตัวเลขธรรมดาเพื่อความง่าย การแปลงแรงผลักจาก contact ให้กลายเป็นการหมุน — คิด torque ใช้ moment of inertia — เป็นคณิตศาสตร์เพิ่มเติมจริง ๆ ที่ต่อยอดจากทุกอย่างในบทนี้ และเป็นหัวข้อถัดไปที่เหมาะจะเรียนพอไอเดียในบทนี้เข้าใจแล้ว ทุกคอนเซปต์ด้านล่าง (integration, broad phase, narrow phase, contacts, impulses) ทำงานแบบเดียวกันใน 3D เต็มรูปแบบ แค่การจัดการเรื่อง rotation จะซับซ้อนขึ้นย้อนกลับไปที่บทแคลคูลัส: physics engine ไม่เคยแก้สมการการเคลื่อนที่แบบตรง ๆ ทั้งหมด มันจะประมาณค่าไปทีละก้าวเล็ก ๆ คงที่ ด้วยวิธี semi-implicit (symplectic) Euler — อัปเดต velocity ก่อนจาก acceleration ปัจจุบัน แล้วค่อยอัปเดต position โดยใช้ velocity ตัวใหม่ ไม่ใช่ตัวเก่า แค่ลำดับบรรทัดเดียวนี้แหละที่ทำให้พลังงานไม่รั่วหายหรือระเบิดออกไป และนี่คือสิ่งที่ physics engine จริง ๆ ทุกตัวทำในทุกสเต็ป มาดูกันว่าเอามาใช้กับ RigidBody แล้วสเต็ปภายใต้แรงโน้มถ่วงเป็นยังไง:
#include <iostream>
struct Vec2 { float x, y; };
struct RigidBody {
Vec2 position;
Vec2 velocity;
float invMass;
};
void integrate(RigidBody& body, Vec2 gravity, float dt) {
if (body.invMass <= 0.0f) return; // static bodies never move
// 1. velocity FIRST, using the current acceleration (here: just gravity)
body.velocity.x += gravity.x * dt;
body.velocity.y += gravity.y * dt;
// 2. position using the NEW velocity -- this is what makes it semi-implicit
body.position.x += body.velocity.x * dt;
body.position.y += body.velocity.y * dt;
}
int main() {
RigidBody ball;
ball.position = {0.0f, 0.0f};
ball.velocity = {3.0f, 0.0f};
ball.invMass = 1.0f;
Vec2 gravity = {0.0f, -10.0f};
float dt = 0.1f;
for (int tick = 1; tick <= 4; tick++) {
integrate(ball, gravity, dt);
std::cout << "tick " << tick << ": pos=(" << ball.position.x << ", " << ball.position.y
<< ") vel=(" << ball.velocity.x << ", " << ball.velocity.y << ")\n";
}
}ผลลัพธ์:
tick 1: pos=(0.3, -0.1) vel=(3, -1)
tick 2: pos=(0.6, -0.3) vel=(3, -2)
tick 3: pos=(0.9, -0.6) vel=(3, -3)
tick 4: pos=(1.2, -1) vel=(3, -4)ถ้าสเต็ป 2 ใช้ velocity ตัวเก่าในการขยับ position แทน (นี่คือ explicit Euler ตัวที่ไม่เสถียรซึ่งบทแคลคูลัสเคยเตือนไว้) ตำแหน่ง y หลังจากสี่ tick เดิมจะกลายเป็น 0, -0.1, -0.3, -0.6 แทนที่จะเป็น -0.1, -0.3, -0.6, -1.0 — ช้าไปหนึ่งสเต็ปเต็ม ๆ เพราะมันขยับ position ด้วย velocity ที่ยังไม่ "เห็น" แรงโน้มถ่วงของสเต็ปนี้เลย จำนวนการคำนวณเท่ากัน ตัวเลขอินพุตเท่ากัน แต่ trajectory ต่างกันอย่างมีนัยสำคัญ (แถมเสถียรน้อยกว่าด้วย) นี่คือลำดับบรรทัดเดียวกันเป๊ะกับที่ตัวอย่าง FixedUpdate ในบท game loop เคยใช้ไปแล้วโดยไม่ได้เอ่ยชื่อมัน — อัปเดต velocityY ก่อนเอาไปขยับ transform.position
ทำไมต้องยืนกรานลำดับ "velocity ก่อน"? เพราะลำดับอีกแบบ — explicit (forward) Euler ที่ขยับ position ด้วย velocity ตัวเก่าแล้วค่อยอัปเดต velocity ทีหลัง — ไม่ได้แค่ช้าไปหนึ่งสเต็ปเฉย ๆ แต่กับอะไรก็ตามที่แกว่งไปมา มันจะเพิ่มพลังงานเข้าไปทุกสเต็ปอย่างเงียบ ๆ จนกระทั่ง simulation ระเบิด ที่ที่เห็นชัดที่สุดคือสปริง: มวลที่ถูกดึงกลับเข้าหา origin ด้วย acceleration a = -k*x ซึ่งในโลกจริงจะแกว่งไปมาตลอดกาลด้วย amplitude คงที่ ลองรัน integrator ทั้งสองตัวกับสปริงตัวเดียวกันนี้ แล้วพิมพ์พลังงานรวม (0.5*(v*v + k*x*x) ซึ่งทางฟิสิกส์ต้องคงที่) ออกมา:
#include <cstdio>
#include <cmath>
int main() {
float k = 1.0f, dt = 0.2f; // spring constant, fixed step
float ex = 1.0f, ev = 0.0f; // explicit Euler: position, velocity
float sx = 1.0f, sv = 0.0f; // semi-implicit Euler: position, velocity
printf("step | explicit x energy | semi-impl x energy\n");
for (int step = 0; step <= 40; step++) {
if (step % 8 == 0) {
float eE = 0.5f * (ev * ev + k * ex * ex);
float sE = 0.5f * (sv * sv + k * sx * sx);
printf("%3d | %9.4f %.4f | %9.4f %.4f\n",
step, ex, eE, sx, sE);
}
// explicit (forward) Euler: move x with OLD v, then v with OLD x
float exNew = ex + ev * dt;
float evNew = ev - k * ex * dt;
ex = exNew; ev = evNew;
// semi-implicit: v FIRST (current x), then x with the NEW v
sv = sv - k * sx * dt;
sx = sx + sv * dt;
}
}ผลลัพธ์:
step | explicit x energy | semi-impl x energy
0 | 1.0000 0.5000 | 1.0000 0.5000
8 | -0.0098 0.6843 | -0.1323 0.5133
16 | -1.3684 0.9365 | -0.9916 0.4936
24 | 0.0402 1.2817 | 0.1955 0.5196
32 | 1.8719 1.7540 | 0.9791 0.4875
40 | -0.0917 2.4005 | -0.2580 0.5256สปริงตัวนี้ควรจะรักษาพลังงานไว้ที่ 0.5 ตลอดไป semi-implicit Euler (สองคอลัมน์ขวา) ทำแบบนั้นเป๊ะ — พลังงานของมันสั่นอยู่รอบ ๆ 0.5 นิดเดียว (0.49 ถึง 0.53) และไม่เคยเลื่อนหนีไปไหน มวลก็แกว่งอยู่ระหว่างประมาณ -1 กับ +1 นานเท่าที่เรารันเลย ส่วน explicit Euler (สองคอลัมน์ซ้าย) กลับไต่ขึ้นไม่หยุด: 0.50, 0.68, 0.94, 1.28, 1.75, 2.40 — เกือบห้าเท่าของพลังงานเริ่มต้นในแค่ 40 สเต็ป และยังโตต่อ ช่วงการแกว่งก็กว้างขึ้นทุกรอบ ปล่อยไปอีกไม่กี่ร้อยสเต็ปตัวเลขก็หนีไปเป็นอนันต์ นี่ไม่ใช่ error จากการปัดเศษที่ย่อให้เล็กลงได้ด้วยการใช้ double แต่เป็น bias เชิงทิศทางที่ฝังอยู่ในลำดับการอัปเดต นี่แหละคือเหตุผลที่ไม่มี engine ตัวไหนใช้ explicit Euler กับ dynamics และเป็นเหตุผลว่าทำไมการสลับลำดับบรรทัดเดียวจากบทแคลคูลัสถึงสำคัญมาก
มี integrator ตัวที่สามที่ควรรู้จักไว้ เพราะเอฟเฟกต์ทั้งตระกูล — ผ้า เชือก เส้นผม และ position-based dynamics — สร้างขึ้นบนมัน นั่นคือ Verlet integration เคล็ดลับของมันคือไม่เก็บ velocity เลย แต่จะเก็บ position ก่อนหน้าไว้แทน แล้วอนุมานการเคลื่อนที่จากช่องว่างระหว่าง "เมื่อกี้ฉันอยู่ไหน" กับ "ตอนนี้ฉันอยู่ไหน":
#include <cstdio>
int main() {
float g = -10.0f, dt = 0.1f; // gravity, fixed step
float y = 0.0f; // current position (start at rest at 0)
float yPrev = 0.0f; // previous position; at rest it equals y
printf("step | y (Verlet) | implied v = (y - yPrev)/dt\n");
for (int step = 1; step <= 5; step++) {
// position Verlet: next = 2*current - previous + a*dt*dt
float yNext = 2.0f * y - yPrev + g * dt * dt;
float impliedV = (yNext - y) / dt; // velocity is never stored -- it is inferred
yPrev = y;
y = yNext;
printf("%3d | %8.2f | %8.1f\n", step, y, impliedV);
}
}ผลลัพธ์:
step | y (Verlet) | implied v = (y - yPrev)/dt
1 | -0.10 | -1.0
2 | -0.30 | -2.0
3 | -0.60 | -3.0
4 | -1.00 | -4.0
5 | -1.50 | -5.0ค่า y พวกนั้น — -0.10, -0.30, -0.60, -1.00, -1.50 — คือการตกแบบเดียวกันเป๊ะกับที่ semi-implicit Euler สร้างไว้เมื่อสองลิสต์ก่อน แต่ได้มาโดยไม่เคยเก็บตัวแปร velocity เลย Verlet เป็น symplectic เหมือนกัน (มันจะไม่ทำให้สปริงระเบิด) และมีคุณสมบัติหนึ่งที่ทำให้มันเป็นตัวเลือกหลักสำหรับผ้ากับเชือก: เพราะไม่มี velocity แยกมาขัดกับ position เราเลยจับจุดแล้วย้ายมันได้เลย — ตรึงมุมธง ดึงเชือกกลับมาที่ความยาวสูงสุด — แล้ว (current - previous) ในสเต็ปถัดไปจะแปลงการเปลี่ยนแปลงนั้นให้เป็น velocity ที่สมเหตุสมผลให้เอง การแก้ position ตรง ๆ แล้วปล่อยให้ integrator กู้ velocity คืนมาฟรี ๆ นี่แหละคือไอเดียทั้งหมดของ position-based dynamics
บท game loop เคยโชว์ให้เห็นว่า deltaTime ของ Update สั่นไปมาในแต่ละเฟรม และบทแคลคูลัสก็สร้าง fixed-timestep-plus-accumulator loop เต็มรูปแบบขึ้นมาแก้ปัญหานี้ collision detection นี่แหละคือเหตุผลที่ physics ข้ามกลไกนั้นไปไม่ได้ มีเหตุผลจริง ๆ สามข้อที่ซ้อนทับกันอยู่:
dt ที่ใหญ่และไม่สม่ำเสมอทำให้ error ของ integrator ตระกูล Euler โตเร็วขึ้น lag spike แค่ครั้งเดียวก็ทำให้หนึ่งสเต็ปผิดเพี้ยนไปไกลได้dt คงที่ตัวเดียวกันทุกครั้ง จะได้ผลลัพธ์เดียวกันเป๊ะเสมอ — สำคัญมากสำหรับ replay และสำหรับให้ simulation ของผู้เล่นสองคนใน networked game ตรงกันspeed * dt ถ้า dt ใหญ่และคาดเดาไม่ได้ การกระโดดก็จะใหญ่และคาดเดาไม่ได้ตามไปด้วย — และหัวข้อ 10 ของบทนี้เองจะโชว์ให้เห็นว่าการกระโดดที่ใหญ่พอสามารถกระโดดข้ามกำแพงบาง ๆ ไปได้เลยโดยไม่มีใครจับได้ยังไง fixed step ทำให้ขนาดการกระโดดนั้นคงที่และรู้ล่วงหน้าได้ ซึ่งเป็นสิ่งเดียวที่ทำให้วิธีแก้ในหัวข้อ 10 เป็นไปได้ตั้งแต่แรกดังนั้นวัตถุทุกชิ้นในโลกจะถูก integrate ด้วย dt คงที่ตัวเดียวกันในทุก physics step — ซึ่งเซ็ตอัพส่วนที่เหลือของบทนี้ได้พอดี: พอรู้ตำแหน่งใหม่ของวัตถุทุกชิ้นในสเต็ปนี้แล้ว engine ก็ต้องหาว่าตัวไหนบ้างที่ตอนนี้ทับกันอยู่
"ใช้ dt คงที่" ไปเจอกำแพงเข้าอย่างจัง: เฟรมจริง ๆ ไม่ได้มาตามตารางคงที่ เฟรมหนึ่งใช้ 12 ms เฟรมถัดไป 22 ms เฟรมถัดไปอีก 9 ms ขึ้นอยู่กับว่า GPU และส่วนที่เหลือของเกมกำลังทำอะไรอยู่ตอนนั้น ถ้าเราส่งระยะเวลาจริงของแต่ละเฟรมเข้า integrate ตรง ๆ physics จะสเต็ปด้วย dt ต่างกันทุกเฟรม — และคำตอบของ Euler ขึ้นอยู่กับ dt การเคลื่อนที่แบบเดียวกันเลยจะไปตกคนละที่ ขึ้นอยู่กับว่าเฟรมมันบังเอิญตกลงมายังไง:
#include <iostream>
// integrate a falling body for a total of 1.0s, split into `steps` equal parts
float fallFor(float totalTime, int steps) {
float y = 0.0f, v = 0.0f, g = -10.0f;
float dt = totalTime / steps;
for (int i = 0; i < steps; i++) {
v += g * dt; // semi-implicit Euler
y += v * dt;
}
return y;
}
int main() {
std::cout << "1.0s of gravity, stepped different ways:\n";
std::cout << " 1 step of 1.00s : y = " << fallFor(1.0f, 1) << "\n";
std::cout << " 4 steps of 0.25s : y = " << fallFor(1.0f, 4) << "\n";
std::cout << " 10 steps of 0.10s: y = " << fallFor(1.0f, 10) << "\n";
}ผลลัพธ์:
1.0s of gravity, stepped different ways:
1 step of 1.00s : y = -10
4 steps of 0.25s : y = -6.25
10 steps of 0.10s: y = -5.5แรงโน้มถ่วงตัวเดียวกันหนึ่งวินาที ทำให้วัตถุไปตกที่ -10 หรือ -6.25 หรือ -5.5 ขึ้นอยู่กับว่าวินาทีนั้นถูกซอยยังไงล้วน ๆ (สเต็ปเล็กกว่าจะใกล้คำตอบจริงที่ -5.0 มากกว่า แต่นั่นไม่ใช่ประเด็นตรงนี้) ความต่างนั้นคือความไม่ determinism: ส่ง dt อะไรก็ตามที่เฟรมบังเอิญใช้เข้าไปให้ physics แล้วผลลัพธ์ก็ขึ้นอยู่กับจังหวะเฟรม — ทำ replay ซ้ำไม่ได้ และรับประกันว่าจะทำให้ผู้เล่นสองคนใน networked game เลื่อนหนีจากกัน วิธีแก้คือ accumulator: บวกเวลาจริงที่ผ่านไปของแต่ละเฟรมเข้าไปในยอดสะสม แล้วค่อยใช้ยอดนั้นเป็นก้อนขนาดคงที่ พร้อมยกเศษที่เหลือไปเฟรมถัดไป
#include <iostream>
int main() {
const float FIXED = 0.25f; // physics ALWAYS steps by this, never anything else
float acc = 0.0f;
float frameTimes[] = {0.30f, 0.10f, 0.28f, 0.40f}; // jittery real frame durations
for (int f = 0; f < 4; f++) {
acc += frameTimes[f];
int stepsThisFrame = 0;
while (acc >= FIXED) {
acc -= FIXED; // integrate(FIXED) would run here
stepsThisFrame++;
}
std::cout << " frame dt=" << frameTimes[f]
<< " -> ran " << stepsThisFrame
<< " fixed step(s), leftover=" << acc << "\n";
}
}ผลลัพธ์:
frame dt=0.3 -> ran 1 fixed step(s), leftover=0.05
frame dt=0.1 -> ran 0 fixed step(s), leftover=0.15
frame dt=0.28 -> ran 1 fixed step(s), leftover=0.18
frame dt=0.4 -> ran 2 fixed step(s), leftover=0.08เฟรมไทม์สั่นสะบัด — 0.30, 0.10, 0.28, 0.40 — แต่ทุก physics step เป็น 0.25 เป๊ะ เฟรมยาว (0.40) รันสองสเต็ปเพื่อไล่ให้ทัน ส่วนเฟรมสั้น (0.10) รันศูนย์สเต็ปแล้วแค่ฝากเวลาไว้ใน accumulator ไว้ใช้ทีหลัง physics ไม่เคยเห็น dt ที่เปลี่ยนไปมาเลยสักครั้ง มันเลยยังคง deterministic เต็มที่ไม่ว่าเฟรมจะกระตุกยังไง นี่คือสิ่งที่ FixedUpdate ของ Unity ทำเป๊ะ ๆ: มันถูกเรียกศูนย์ หนึ่ง หรือหลายครั้งต่อเฟรมที่เรนเดอร์ เพื่อให้ physics เดินหน้าเป็นก้อนขนาด Time.fixedDeltaTime เสมอ มีอีกเรื่องหนึ่งที่ควรรู้: เฟรมที่เรนเดอร์มักตกลงระหว่างสองสเต็ปคงที่ (เศษที่ค้างอยู่ใน accumulator นั่นแหละ) engine ที่ลื่น ๆ เลยจะ interpolate ตำแหน่งภาพระหว่าง state ของ physics สองตัวล่าสุด — ไม่งั้นวัตถุที่เร็ว ๆ จะกระตุกให้เห็นทั้ง ๆ ที่ simulation ข้างใต้สมบูรณ์แบบ
ถ้ามีวัตถุ n ชิ้นในโลก จะมีคู่ที่อาจจะแตะกันได้ทั้งหมด n * (n-1) / 2 คู่ การเทสต์ทุกคู่ด้วยการเช็กรูปร่างแบบละเอียดเป็นการเสียเวลาโดยใช่เหตุ: วัตถุ 100 ชิ้นหมายถึงต้องเช็ก 4,950 คู่ในทุก physics step ทั้ง ๆ ที่เกือบทั้งหมดอยู่ห่างกันคนละที่ งานเดียวของ broad phase collision detection คือคัดคู่ที่เห็นชัด ๆ ว่าไม่มีทางแตะกันทิ้งไปแบบไว ๆ และถูก ๆ โดยใช้ AABB จากบทเรขาคณิต เหลือแค่ลิสต์คู่ candidate เล็ก ๆ ไว้ให้เทสต์แบบละเอียดต่อไป
นี่คือวงกลมสี่วง โดยใช้ struct Circle/AABB จากบทเรขาคณิต พร้อม broad phase แบบง่าย ๆ ที่เช็ก AABB ทับกันของทุกคู่:
#include <iostream>
struct Vec2 { float x, y; };
struct AABB { Vec2 min, max; };
AABB circleAABB(Vec2 center, float radius) {
return { {center.x - radius, center.y - radius},
{center.x + radius, center.y + radius} };
}
bool aabbOverlap(AABB a, AABB b) {
if (a.max.x < b.min.x || a.min.x > b.max.x) return false;
if (a.max.y < b.min.y || a.min.y > b.max.y) return false;
return true;
}
int main() {
struct Named { const char* name; Vec2 center; float radius; };
Named bodies[] = {
{"A", {0.0f, 0.0f}, 1.0f},
{"B", {1.7f, 0.0f}, 1.0f},
{"C", {5.0f, 5.0f}, 1.0f},
{"D", {0.0f, 3.0f}, 1.0f},
};
int pairsChecked = 0, candidates = 0;
for (int i = 0; i < 4; i++) {
for (int j = i + 1; j < 4; j++) {
AABB boxI = circleAABB(bodies[i].center, bodies[i].radius);
AABB boxJ = circleAABB(bodies[j].center, bodies[j].radius);
pairsChecked++;
if (aabbOverlap(boxI, boxJ)) {
std::cout << bodies[i].name << "-" << bodies[j].name << ": CANDIDATE\n";
candidates++;
}
}
}
std::cout << "checked " << pairsChecked << " pairs, found " << candidates << " candidate(s)\n";
}ผลลัพธ์:
A-B: CANDIDATE
checked 6 pairs, found 1 candidate(s)เช็กไปหกคู่ เหลือรอดมาคู่เดียว คู่นั้นคือ A กับ B ซึ่งเป็นคู่เดียวที่ narrow phase (หัวข้อ 4) ต้องไปยุ่งด้วย เวอร์ชันง่าย ๆ นี้ยังนับว่าเป็น O(n^2) ในแง่จำนวนคู่ที่ต้อง "เช็ก" อยู่ดี (แค่ทำให้แต่ละการเช็กถูกลง) ซึ่งก็โอเคถ้ามีวัตถุแค่ไม่กี่ชิ้น แต่พอฉากมีวัตถุเป็นร้อย ๆ ตัว มันจะเริ่มสะสมจนช้าลง มีวิธีทั่วไปสองแบบที่หลีกเลี่ยงการเช็กทุกคู่ไปเลย:
แบ่งโลกออกเป็นเซลล์สี่เหลี่ยมขนาดคงที่ แล้วโยนแต่ละวัตถุลงในเซลล์ที่จุดศูนย์กลางของมันตกอยู่ (ใช้ floor(position / cellSize)) จะมีแค่วัตถุที่อยู่เซลล์เดียวกัน (หรือเซลล์ข้างเคียง สำหรับวัตถุที่อยู่ใกล้ขอบเซลล์) เท่านั้นที่จะถูกเอามาเทียบกัน วัตถุที่อยู่เซลล์ห่างกันมากจะไม่ถูกพิจารณาเลย:
#include <iostream>
#include <cmath>
struct Vec2 { float x, y; };
int cellOf(float coord, float cellSize) {
return (int)std::floor(coord / cellSize);
}
int main() {
struct Named { const char* name; Vec2 center; };
Named bodies[] = {
{"A", {0.0f, 0.0f}},
{"B", {1.7f, 0.0f}},
{"C", {5.0f, 5.0f}},
{"D", {0.0f, 3.0f}},
};
float cellSize = 2.0f;
for (auto& b : bodies) {
int cx = cellOf(b.center.x, cellSize);
int cy = cellOf(b.center.y, cellSize);
std::cout << b.name << ": cell (" << cx << ", " << cy << ")\n";
}
}ผลลัพธ์:
A: cell (0, 0)
B: cell (0, 0)
C: cell (2, 2)
D: cell (0, 1)A กับ B ตกเซลล์เดียวกัน ส่วน C กับ D ต่างก็อยู่เซลล์ของตัวเองคนเดียว การ implement จริงจะจัดกลุ่มวัตถุตามเซลล์ (ปกติใช้ hash map ที่ key เป็นพิกัดเซลล์) แล้วเทสต์แค่คู่ที่อยู่เซลล์เดียวกันหรือเซลล์ข้างเคียง — ได้คู่ candidate เดียวกับที่ AABB sweep เจอ แค่ไปถึงด้วยวิธีคนละแบบ grid นั้นเรียบง่ายและไวเมื่อวัตถุกระจายตัวสม่ำเสมอ แต่จะเปลืองหน่วยความจำไปกับเซลล์ที่ว่างเปล่าเป็นส่วนใหญ่ ถ้าวัตถุไปกระจุกกันแน่น ๆ ในพื้นที่เดียว หรือกระจายอยู่ในโลกที่ใหญ่และโปร่งมาก
Sweep and prune (บางทีก็เรียกว่าวิธี "sort and sweep") มองปัญหาจากมุมที่ต่างออกไป: เรียง AABB ของแต่ละวัตถุตามพิกัด x ต่ำสุด แล้วไล่เดินไปตามลิสต์ที่เรียงแล้วจากซ้ายไปขวา คอยจดว่า AABB ตัวไหน "เปิด" อยู่ตอนนี้บ้าง (ช่วงของมันเริ่มแล้วแต่ยังไม่จบ) AABB สองตัวจะมีโอกาสทับกันได้ก็ต่อเมื่อทั้งคู่เปิดอยู่พร้อมกันเท่านั้น:
engine อย่าง Box2D ชอบใช้ sweep and prune เพราะเพื่อนบ้านของวัตถุแต่ละตัวแทบไม่เปลี่ยนจากเฟรมหนึ่งไปอีกเฟรม — ลำดับที่เรียงไว้จากเฟรมก่อนแทบจะถูกอยู่แล้ว การเรียงใหม่เลยแทบไม่เสียอะไรเลย (สลับกันแค่ไม่กี่ตัว ไม่ใช่เรียงใหม่ทั้งหมด) ในเกือบทุกสเต็ป ไม่ว่า broad phase จะใช้วิธีไหน งานของมันก็เหมือนเดิมเสมอ: เปลี่ยน "คู่ที่เป็นไปได้ทั้งหมด" ให้กลายเป็น "คู่ candidate" อย่างไว เพื่อให้การเทสต์แบบละเอียด (ที่แพงกว่ามาก) ในหัวข้อถัดไปต้องรันแค่ไม่กี่ครั้ง ไม่ใช่เป็นพัน
O(n^2) exact shape test ทุก physics step ทันที แล้ว frame rate ก็จะร่วงลงเหวเลย สร้างนิสัยใช้ broad phase ไว้ตั้งแต่เนิ่น ๆ ก่อนที่มันจะกลายเป็นไฟที่ต้องมาดับทีหลังNarrow phase collision detection เอาลิสต์คู่ candidate เล็ก ๆ ที่ broad phase สร้างมา แล้วรันการเทสต์แบบละเอียดเฉพาะรูปร่างนั้น ๆ กับแต่ละคู่: รูปร่างสองอันนี้ทับกันจริงไหม แล้วถ้าใช่ ทับกันแค่ไหน สำหรับวงกลมสองวงเรื่องนี้ง่าย — เป็นการเทสต์แบบละเอียดที่ต่อยอดตรง ๆ จากไอเดีย point-in-circle ในบทเรขาคณิต:
#include <iostream>
#include <cmath>
struct Vec2 { float x, y; };
Vec2 sub(Vec2 a, Vec2 b) { return { a.x - b.x, a.y - b.y }; }
float length(Vec2 v) { return std::sqrt(v.x * v.x + v.y * v.y); }
bool circlesOverlap(Vec2 centerA, float rA, Vec2 centerB, float rB) {
float dist = length(sub(centerB, centerA));
return dist < (rA + rB);
}
int main() {
Vec2 A = {0.0f, 0.0f}, B = {1.7f, 0.0f};
std::cout << "A-B exact test: " << (circlesOverlap(A, 1.0f, B, 1.0f) ? "COLLIDING" : "clear") << "\n";
}ผลลัพธ์:
A-B exact test: COLLIDINGbroad phase บอกเราไปแล้วว่า A กับ B เป็นคู่เดียวที่คุ้มค่าจะเช็ก narrow phase เลยมีแค่การเทสต์เดียวที่ต้องรัน ไม่ใช่หกครั้งแบบที่การเช็กทุกคู่ตรง ๆ ต้องทำ "จุดศูนย์กลางสองจุดอยู่ใกล้กันน้อยกว่าผลรวมของรัศมี" คือการเทสต์แบบละเอียดทั้งหมดสำหรับวงกลม นี่คือเหตุผลที่วงกลม (และทรงกลมใน 3D) เป็นรูปร่างที่ถูกที่สุดที่ physics engine จะเอามาชนกันได้
ระหว่าง "วงกลมสองวง" กับ "polygon สองอัน" มีเคสหนึ่งที่พบบ่อยพอจะมีทางลัดของตัวเอง คือวงกลมกับกล่องแบบ axis-aligned การเทสต์แบบละเอียดเป็นทริกสวย ๆ — หาจุดบนกล่องที่ใกล้ศูนย์กลางวงกลมที่สุด ด้วยการ clamp ศูนย์กลางเข้าไปในช่วงของกล่องในแต่ละแกน จากนั้นก็เหลือแค่การเช็กวงกลมกับจุดเท่านั้น:
#include <iostream>
#include <cmath>
struct Vec2 { float x, y; };
struct AABB { Vec2 min, max; };
float clampf(float v, float lo, float hi) {
return v < lo ? lo : (v > hi ? hi : v);
}
int main() {
AABB box = { {0.0f, 0.0f}, {3.0f, 2.0f} };
Vec2 center = {3.5f, 1.0f};
float radius = 1.0f;
// closest point on the box to the circle center = clamp center into the box
Vec2 closest = { clampf(center.x, box.min.x, box.max.x),
clampf(center.y, box.min.y, box.max.y) };
Vec2 d = { center.x - closest.x, center.y - closest.y };
float dist = std::sqrt(d.x * d.x + d.y * d.y);
std::cout << "closest point on box = (" << closest.x << ", " << closest.y << ")\n";
std::cout << "distance to center = " << dist << "\n";
if (dist < radius) {
Vec2 normal = { d.x / dist, d.y / dist }; // from box toward circle
float penetration = radius - dist;
std::cout << "COLLIDING normal=(" << normal.x << ", " << normal.y
<< ") penetration=" << penetration << "\n";
} else {
std::cout << "clear\n";
}
}ผลลัพธ์:
closest point on box = (3, 1)
distance to center = 0.5
COLLIDING normal=(1, 0) penetration=0.5วงกลมอยู่ที่ x=3.5 เลยขอบขวาของกล่องที่ x=3 ไปนิดเดียว การ clamp ศูนย์กลางเข้าช่วง [0,3] x [0,2] ของกล่องจะดึงมันไปที่ (3, 1) — จุดที่ใกล้ที่สุดบนกล่อง — ซึ่งอยู่ห่างไป 0.5 อยู่ในรัศมี 1 มันเลยชนกัน 0.5 โดยมี normal ชี้ตรงออกไปตามแนว +x เคสเดียวที่เวอร์ชันง่าย ๆ นี้จัดการไม่ได้คือตอนที่ศูนย์กลางอยู่ข้างในกล่อง (ตอนนั้น clamp จะคืนศูนย์กลางกลับมาเอง dist เป็น 0 แล้วเราหารเพื่อหา normal ไม่ได้) การ implement เต็ม ๆ จะตรวจจับเคสนั้นแล้วผลักออกไปตามหน้าที่ใกล้ที่สุดแทน engine เก็บทางลัดเขียนมือแบบนี้ไว้สำหรับ circle-vs-circle, circle-vs-box และ box-vs-box ก็เพราะมันถูกขนาดนี้แหละ
วงกลมเป็นเคสง่าย สำหรับกล่องและ convex polygon อื่น ๆ (รูปร่างจะเป็น convex ถ้าเส้นตรงระหว่างจุดสองจุดใด ๆ ข้างในมันไม่มีทางออกนอกรูปร่างเลย) การเทสต์แบบละเอียดคลาสสิกคือ Separating Axis Theorem (SAT): รูปทรง convex สองอันจะไม่ชนกันก็ต่อเมื่อเราหาแกนได้อย่างน้อยหนึ่งแกนที่พอฉายรูปทรงทั้งสองลงไปแล้วมีช่องว่างระหว่างกัน ลองทุกแกนที่เป็นไปได้ (สำหรับ polygon หมายถึง normal ที่ตั้งฉากกับทุกขอบ) ถ้าทุกแกนแสดงว่าทับกันหมด แสดงว่ารูปทรงชนกันจริง และแกนที่มีช่วงทับซ้อนน้อยที่สุดจะบอก collision normal ให้เรา
นี่คือทฤษฎีนั้นในรูปแบบการเทสต์ที่รันได้จริง สำหรับแต่ละแกนที่เป็นไปได้ (edge normal ทุกอันของรูปทรงทั้งสอง) ให้ ฉาย polygon ทั้งสองลงบนมัน — ย่อแต่ละอันให้เหลือแค่เงา [min, max] ที่มันทอดลงบนแกนนั้น — แล้วเช็กว่าเงาทับกันไหม พอมีแกนใดแกนหนึ่งแสดงช่องว่าง ให้หยุดเลย: รูปทรงแยกกันอยู่ ถ้าทุกแกนทับกันหมด แสดงว่าชนกัน และแกนที่มีช่วงทับน้อยที่สุดคือ collision normal ส่วนช่วงทับนั้นคือ penetration depth:
#include <iostream>
#include <cmath>
#include <vector>
struct Vec2 { float x, y; };
float dot(Vec2 a, Vec2 b) { return a.x * b.x + a.y * b.y; }
// project every vertex of `poly` onto `axis`, return the [min,max] shadow
void projectOnto(const std::vector<Vec2>& poly, Vec2 axis, float& mn, float& mx) {
mn = mx = dot(poly[0], axis);
for (const Vec2& v : poly) {
float p = dot(v, axis);
if (p < mn) mn = p;
if (p > mx) mx = p;
}
}
bool satOverlap(const std::vector<Vec2>& A, const std::vector<Vec2>& B,
Vec2& normal, float& depth) {
depth = 1e30f;
for (int poly = 0; poly < 2; poly++) { // axes from BOTH shapes
const std::vector<Vec2>& P = (poly == 0) ? A : B;
for (size_t i = 0; i < P.size(); i++) {
Vec2 a = P[i], b = P[(i + 1) % P.size()];
Vec2 edge = { b.x - a.x, b.y - a.y };
Vec2 axis = { -edge.y, edge.x }; // edge normal (perpendicular)
float len = std::sqrt(axis.x * axis.x + axis.y * axis.y);
axis = { axis.x / len, axis.y / len };
float minA, maxA, minB, maxB;
projectOnto(A, axis, minA, maxA);
projectOnto(B, axis, minB, maxB);
float overlap = std::min(maxA, maxB) - std::max(minA, minB);
if (overlap <= 0.0f) return false; // a gap -> separated, done
if (overlap < depth) { depth = overlap; normal = axis; }
}
}
return true;
}
void test(const char* name, std::vector<Vec2> A, std::vector<Vec2> B) {
Vec2 n; float d;
if (satOverlap(A, B, n, d))
std::cout << name << ": COLLIDING normal=(" << n.x << ", " << n.y
<< ") depth=" << d << "\n";
else
std::cout << name << ": separated (a gap axis was found)\n";
}
int main() {
std::vector<Vec2> A = {{0,0},{2,0},{2,2},{0,2}}; // axis-aligned square
std::vector<Vec2> Bhit = {{2.8f,0},{3.8f,1},{2.8f,2},{1.8f,1}}; // diamond, overlaps A
std::vector<Vec2> Bmiss = {{4.2f,0},{5.2f,1},{4.2f,2},{3.2f,1}}; // diamond, clears A
test("overlapping", A, Bhit);
test("separated ", A, Bmiss);
}ผลลัพธ์:
overlapping: COLLIDING normal=(-1, 0) depth=0.2
separated : separated (a gap axis was found)กล่องใบที่สองเป็นสี่เหลี่ยมข้าวหลามตัด (สี่เหลี่ยมจัตุรัสหมุน 45 องศา) นี่เลยเป็น SAT ของจริง ไม่ใช่การเช็ก AABB ที่ปลอมตัวมา — การฉายลงบน edge normal เอียง ๆ ของข้าวหลามตัดคืองานจริง ๆ ในเคสที่ทับกัน สี่เหลี่ยม A กินช่วง x เป็น [0, 2] ส่วนปลายซ้ายสุดของข้าวหลามตัดอยู่ที่ x=1.8 เงาของมันบนแกน x เลยทับกันแค่ 0.2 — น้อยกว่าช่วงทับบนแกนอื่น ๆ ที่ลองมาทั้งหมด แกนนั้นเลยชนะเป็น collision normal โดยมี penetration 0.2 (เครื่องหมายของ normal ขึ้นอยู่กับว่าขอบไหนสร้างมันขึ้นมา engine จริง ๆ จะกลับด้านให้ชี้จาก A ไปทาง B เสมอ ที่นี่มันเลยจะรายงานเป็น (1, 0)) ส่วนในเคสที่แยกกัน ข้าวหลามตัดถูกเลื่อนไปทางขวาจนปลายของมันพ้นขอบ A แล้วแกนแรกที่แสดงช่องว่างก็ทำให้ satOverlap คืนค่าออกไปเลยทันที — การออกก่อนกำหนดนี่แหละที่ทำให้ SAT ไวในเคสทั่วไปที่ไม่ได้แตะกัน
SAT ใช้ได้ดีกับ polygon เพราะ "ลองทุก edge normal" เป็นลิสต์แกนที่สั้นและรู้ล่วงหน้าได้ แต่สำหรับรูปทรง convex ทั่วไป — capsule โค้งมน หรือ convex mesh แบบใดก็ได้ — ไม่มีลิสต์แกนตายตัวสั้น ๆ ให้ลอง คำตอบมาตรฐานคือ GJK algorithm (Gilbert-Johnson-Keerthi ตั้งชื่อตามผู้เขียนทั้งสามคน) ซึ่งเลี่ยงปัญหาการค้นหาแกนไปเลย GJK ทำงานบน Minkowski difference ของรูปทรงสองอัน (พูดง่าย ๆ คือ "รูปทรง A ลบรูปทรง B" สร้างขึ้นโดยเอาทุกจุดของ B ไปลบออกจากทุกจุดของ A) — ข้อเท็จจริงที่สวยงามของโครงสร้างนี้คือ รูปทรงต้นฉบับทั้งสองจะทับกันก็ต่อเมื่อรูปทรงที่รวมกันนี้มีจุด origin (0,0) อยู่ข้างในเท่านั้น
GJK ไม่จำเป็นต้องสร้างรูปทรง Minkowski difference เต็มรูปแบบเลย (ซึ่งอาจใหญ่มาก) แต่จะใช้ support function แทน — "ขอจุดของรูปทรงนี้ที่อยู่ไกลที่สุดในทิศทาง D" — เพื่อค่อย ๆ ขยายสามเหลี่ยมเล็ก ๆ (เรียกว่า simplex) ไปทีละขั้น แต่ละครั้งจะเลือกทิศทางใหม่ที่ทำให้ simplex เข้าใกล้การล้อมรอบจุด origin มากขึ้น ถ้า simplex ล้อมรอบ origin ได้เมื่อไหร่ แสดงว่ารูปทรงทับกัน แต่ถ้าทิศทางการค้นหาชี้ออกจาก origin โดยไม่มีทางเข้าใกล้ได้อีก แสดงว่าไม่ทับกัน อัลกอริทึมเต็มรูปแบบนั้น implement ให้ถูกต้องยากจริง ๆ ซึ่งเป็นเหตุผลที่แทบไม่มีใคร hand-write มันเอง — PhysX backend ของ Unity, Chaos physics ของ Unreal และ physics engine ทั่วไปแทบทุกตัวมากับ GJK ที่ผ่านการพิสูจน์มาแล้ว (มักจะจับคู่กับอัลกอริทึมที่เกี่ยวข้องอย่าง EPA เพื่อหา penetration depth หลังจาก GJK ยืนยันว่าทับกัน) สิ่งที่สำคัญสำหรับเราคือจำชื่อนี้ให้ได้และรู้ว่ามันแก้ปัญหาอะไร: การเทสต์การทับกันแบบละเอียดสำหรับรูปทรง convex ทั่วไป ในกรณีที่ไม่มีลิสต์แกนสั้น ๆ แบบ SAT ให้ใช้
circlesOverlap ในหัวข้อ 4) แทนที่จะส่งผ่าน GJK — มันไวกว่าถ้าข้ามกลไกทั่วไปไปเมื่อรู้ทางลัดอยู่แล้ว GJK คือตัวที่จัดการกับทุกอย่างที่เหลือ: capsule, convex hull, convex mesh แบบใดก็ได้คำตอบแบบ yes/no ว่า "ชนกันไหม" ยังไม่พอที่จะแก้ปัญหาการชนได้จริง — engine ต้องรู้ด้วยว่ารูปทรงแตะกันตรงไหนแน่ ๆ และต้องผลักแยกกันไปทางไหน เอาต์พุตจริงของ narrow phase หลังจากยืนยันว่าทับกันแล้ว คือ struct เล็ก ๆ ที่เรียกว่า contact (หรือ manifold เมื่อคู่เดียวสร้าง contact point ได้มากกว่าหนึ่งจุด เช่น กล่องที่วางราบอยู่บนพื้น):
สำหรับวงกลมสองวง ทั้งสามอย่างนี้ได้มาตรงจากการเทสต์แบบละเอียดในหัวข้อ 4 เลย — เส้นจากศูนย์กลางถึงศูนย์กลางก็คือทิศทาง normal อยู่แล้ว:
#include <iostream>
#include <cmath>
struct Vec2 { float x, y; };
Vec2 add(Vec2 a, Vec2 b) { return { a.x + b.x, a.y + b.y }; }
Vec2 sub(Vec2 a, Vec2 b) { return { a.x - b.x, a.y - b.y }; }
Vec2 scale(Vec2 a, float s) { return { a.x * s, a.y * s }; }
float length(Vec2 v) { return std::sqrt(v.x * v.x + v.y * v.y); }
struct Contact {
Vec2 point;
Vec2 normal; // points from A toward B
float penetration;
};
Contact makeCircleContact(Vec2 posA, float rA, Vec2 posB, float rB) {
Vec2 delta = sub(posB, posA);
float dist = length(delta);
Vec2 normal = scale(delta, 1.0f / dist);
Vec2 pointOnA = add(posA, scale(normal, rA));
Vec2 pointOnB = sub(posB, scale(normal, rB));
Vec2 midpoint = scale(add(pointOnA, pointOnB), 0.5f);
float penetration = (rA + rB) - dist;
return { midpoint, normal, penetration };
}
int main() {
Contact c = makeCircleContact({0.0f, 0.0f}, 1.0f, {1.7f, 0.0f}, 1.0f);
std::cout << "point=(" << c.point.x << ", " << c.point.y << ")\n";
std::cout << "normal=(" << c.normal.x << ", " << c.normal.y << ")\n";
std::cout << "penetration=" << c.penetration << "\n";
}ผลลัพธ์:
point=(0.85, 0)
normal=(1, 0)
penetration=0.3A (ศูนย์กลางอยู่ที่ origin) กับ B (ศูนย์กลางอยู่ทางขวา 1.7 หน่วย) ทับกันอยู่ 0.3 หน่วย contact จะอยู่ตรงกึ่งกลางของช่วงที่ทับกันนั้น และ normal ก็ชี้ตรงจาก A ไปทาง B ทุกหัวข้อที่เหลือในบทนี้ทำงานอยู่บน struct เล็ก ๆ ตัวนี้ทั้งหมด — มันคือจุดส่งต่อระหว่าง "ชนกันไหม" (หัวข้อ 1-4) กับ "แล้วเราจะทำยังไงกับมัน" (หัวข้อ 6 เป็นต้นไป)
ถ้าปล่อยไว้เฉย ๆ วัตถุสองชิ้นที่ทับกันก็จะติดค้างอยู่ในกันและกันแบบนั้น Positional correction จะขยับวัตถุสองชิ้นให้แยกออกจากกันตรง ๆ ตามแนว contact normal ให้พอดีกับที่ทำให้ penetration หายไป ว่าวัตถุไหนจะขยับมากแค่ไหนก็ใช้ตรรกะ invMass เดิมจากหัวข้อ 1: วัตถุที่มี inverse mass มากกว่า (เบากว่า ผลักง่ายกว่า) จะถูกผลักไปไกลกว่า ส่วนวัตถุ static (invMass = 0) จะไม่ขยับเลย และวัตถุอีกตัวจะรับ correction ไปเต็ม ๆ:
#include <iostream>
struct Vec2 { float x, y; };
struct RigidBody { Vec2 position; float invMass; };
struct Contact { Vec2 point; Vec2 normal; float penetration; };
Vec2 add(Vec2 a, Vec2 b) { return { a.x + b.x, a.y + b.y }; }
Vec2 sub(Vec2 a, Vec2 b) { return { a.x - b.x, a.y - b.y }; }
Vec2 scale(Vec2 a, float s) { return { a.x * s, a.y * s }; }
void correctPenetration(RigidBody& a, RigidBody& b, Contact c) {
float totalInvMass = a.invMass + b.invMass;
if (totalInvMass <= 0.0f) return; // both static -- nothing can move
float shareA = a.invMass / totalInvMass;
float shareB = b.invMass / totalInvMass;
Vec2 correction = scale(c.normal, c.penetration);
a.position = sub(a.position, scale(correction, shareA));
b.position = add(b.position, scale(correction, shareB));
}
int main() {
RigidBody a = { {0.0f, 0.0f}, 0.5f }; // mass 2 -> invMass 0.5
RigidBody b = { {1.7f, 0.0f}, 1.0f }; // mass 1 -> invMass 1.0
Contact c = { {0.85f, 0.0f}, {1.0f, 0.0f}, 0.3f };
correctPenetration(a, b, c);
std::cout << "a.position=(" << a.position.x << ", " << a.position.y << ")\n";
std::cout << "b.position=(" << b.position.x << ", " << b.position.y << ")\n";
}ผลลัพธ์:
a.position=(-0.1, 0)
b.position=(1.9, 0)วัตถุ A (หนักเป็นสองเท่า inverse mass ครึ่งเดียว) ขยับไปทางซ้ายแค่ 0.1 หน่วย ส่วนวัตถุ B (เบากว่า inverse mass เต็ม) ขยับไปทางขวา 0.2 หน่วย — รวมกันแล้วครอบคลุม penetration เต็ม 0.3 พอดี แบ่งกันตามสัดส่วนว่าแต่ละตัว "ผลักง่าย" แค่ไหน ลองเช็กระยะห่างใหม่ระหว่างสองตัว: 1.9 - (-0.1) = 2.0 ซึ่งเท่ากับผลรวมของรัศมีพอดีเป๊ะ — การทับกันหายไปหมดแล้ว
การแก้ position ตรง ๆ แบบข้างบนใช้ได้ แต่มันไปสู้กับ velocity solver ที่กำลังจะมาในหัวข้อ 7 — เราขยับวัตถุไปทางหนึ่ง แล้วสเต็ป impulse ก็ขยับมันไปอีกทาง engine ส่วนใหญ่เลยเลือกรวมการแก้ penetration เข้าไปในการแก้ velocity ด้วยวิธี Baumgarte stabilization: เปลี่ยนการทับที่เหลือให้เป็นความเร็วแยกจากกันเพิ่มขึ้นนิดหน่อย เรียกว่า bias แล้วปล่อยให้ normal impulse เป็นตัวส่งมันไป bias คือ (beta / dt) * (penetration - slop) โดยที่ beta เป็นเศษส่วนเล็ก ๆ (ปกติ 0.1 ถึง 0.2) ที่กำหนดว่าจะปิดช่องว่างดุแค่ไหน:
#include <iostream>
int main() {
float penetration = 0.30f; // how deep the bodies overlap
float slop = 0.01f; // permitted overlap left uncorrected
float beta = 0.20f; // Baumgarte factor (fraction closed per step)
float dt = 1.0f / 60.0f;
float correctable = penetration - slop;
if (correctable < 0.0f) correctable = 0.0f;
float biasVel = (beta / dt) * correctable; // added to the target separating speed
float pushThisStep = biasVel * dt; // how far it moves them apart this step
std::cout << "penetration = " << penetration << "\n";
std::cout << "bias velocity = " << biasVel << " units/s\n";
std::cout << "push-out this dt = " << pushThisStep << "\n";
std::cout << "= beta*(pen-slop)= " << beta * correctable << " (same number)\n";
}ผลลัพธ์:
penetration = 0.3
bias velocity = 3.48 units/s
push-out this dt = 0.058
= beta*(pen-slop)= 0.058 (same number)จุดสำคัญอยู่ในสองบรรทัดสุดท้าย: การดันวัตถุแยกออกจากกันด้วย bias velocity 3.48 เป็นเวลาหนึ่งสเต็ป ขยับมันไป 0.058 — ซึ่งเท่ากับ beta * (penetration - slop) = 0.2 * 0.29 พอดี ตัว /dt ใน bias กับ *dt ของหนึ่งสเต็ป integration ตัดกันพอดี beta เลยเป็นแค่ "จะกำจัดการทับกันกี่เศษส่วนต่อสเต็ป" — ไอเดีย 20% ต่อสเต็ปเดียวกับ split positional correction ข้างบน แค่ส่งมาในรูปของความเร็วที่ impulse solver เดิมรู้จักวิธีใช้อยู่แล้ว นั่นคือเหตุผลที่มันเป็นที่นิยม: ไม่ต้องมีสเต็ปแก้ position แยก solver ตัวเดียวทำงานทั้งสองอย่าง
beta สูงเกินไป (เช่น 0.8 เพื่อ "ปิดช่องว่างให้ไว") เพราะ bias เป็นความเร็วจริงที่เพิ่มเข้าไปในวัตถุ ส่วนที่เหลือของมันหลังจากรูปทรงแยกกันแล้วจะกลายเป็น momentum จริงที่วัตถุเก็บไว้ — กองวัตถุจะได้พลังงานเพิ่มแล้วเด้งหรือสั่น ซึ่งเป็นอาการ pop ที่วิธีนี้ควรจะป้องกันพอดี การตั้ง beta ให้ต่ำ (0.1 ถึง 0.2) คือวิธีแก้ทั่วไป ส่วนวิธีที่ครบถ้วนกว่าที่ Box2D และตัวอื่น ๆ ใช้คือ split impulse (หรือ position projection): ใช้ bias บน "pseudo-velocity" ที่ทิ้งได้ต่างหาก ซึ่งดันรูปทรงให้แยกกันเพื่อจัดตำแหน่ง แต่ถูกโยนทิ้งก่อนที่มันจะป้อนกลับเข้าไปใน momentum จริงการแก้ position ทำให้วัตถุหยุดทับกัน แต่ไม่ได้บอกอะไรเลยว่าหลังชนแล้วมันควรจะเคลื่อนที่ยังไง — ลูกบอลควรจะเด้ง ไม่ใช่แค่เทเลพอร์ตออกจากพื้นแล้วตกทะลุมันต่อในสเต็ปถัดไป นี่คืองานของ impulse: การเปลี่ยน velocity แบบทันทีทันใด ใช้ตรง ๆ (ไม่ค่อย ๆ สะสมแบบที่แรงที่กระทำตลอดช่วงเวลาทำ) ตามแนว contact normal โดยเฉพาะ
ส่วนประกอบสำคัญคือ relative velocity ตามแนว normal — วัตถุสองชิ้นเข้าใกล้กันเร็วแค่ไหน วัดตามทิศทางที่มันกำลังจะแยกออกจากกัน ถ้าค่านี้เป็นศูนย์หรือบวกอยู่แล้ว แสดงว่ามันกำลังแยกออกจากกัน (หรือแค่แตะกันพอดี) เองอยู่แล้ว ไม่ต้องใช้ impulse เลย แต่ถ้าเป็นลบ แสดงว่ามันกำลังเข้าใกล้กัน impulse ต้องยกเลิกการเข้าใกล้นั้น และขึ้นอยู่กับ restitution (ค่าสัมประสิทธิ์ความเด้ง ตั้งแต่ 0 = วัตถุติดกัน "inelastic เต็มที่" ไปจนถึง 1 = เด้งแบบ elastic สมบูรณ์ไม่มีพลังงานสูญเสีย) ส่งมันแยกออกจากกันอีกครั้ง:
สูตรขนาดของ impulse (เป็นผลลัพธ์มาตรฐาน ได้มาจากการกำหนดให้ momentum ต้องอนุรักษ์ และความเร็วที่แยกออกจากกันหลังชนต้องเท่ากับ restitution * approaching speed) คือ:
#include <iostream>
struct Vec2 { float x, y; };
struct RigidBody { Vec2 velocity; float invMass; };
struct Contact { Vec2 normal; };
Vec2 add(Vec2 a, Vec2 b) { return { a.x + b.x, a.y + b.y }; }
Vec2 sub(Vec2 a, Vec2 b) { return { a.x - b.x, a.y - b.y }; }
Vec2 scale(Vec2 a, float s) { return { a.x * s, a.y * s }; }
float dot(Vec2 a, Vec2 b) { return a.x * b.x + a.y * b.y; }
float resolveVelocity(RigidBody& a, RigidBody& b, Contact c, float restitution) {
Vec2 relVel = sub(b.velocity, a.velocity);
float velAlongNormal = dot(relVel, c.normal);
if (velAlongNormal >= 0.0f) return 0.0f; // already separating -- nothing to do
float totalInvMass = a.invMass + b.invMass;
float j = -(1.0f + restitution) * velAlongNormal / totalInvMass;
Vec2 impulse = scale(c.normal, j);
a.velocity = sub(a.velocity, scale(impulse, a.invMass));
b.velocity = add(b.velocity, scale(impulse, b.invMass));
return j;
}
int main() {
RigidBody a = { {4.0f, 0.0f}, 0.5f }; // mass 2
RigidBody b = { {-2.0f, 0.0f}, 1.0f }; // mass 1
Contact c = { {1.0f, 0.0f} };
float momentumBefore = (1.0f / a.invMass) * a.velocity.x + (1.0f / b.invMass) * b.velocity.x;
float j = resolveVelocity(a, b, c, 0.5f);
float momentumAfter = (1.0f / a.invMass) * a.velocity.x + (1.0f / b.invMass) * b.velocity.x;
std::cout << "impulse j=" << j << "\n";
std::cout << "a.velocity=(" << a.velocity.x << ", " << a.velocity.y << ")\n";
std::cout << "b.velocity=(" << b.velocity.x << ", " << b.velocity.y << ")\n";
std::cout << "momentum before=" << momentumBefore << " after=" << momentumAfter << "\n";
}ผลลัพธ์:
impulse j=6
a.velocity=(1, 0)
b.velocity=(4, 0)
momentum before=6 after=6วัตถุ A (มวล 2) กำลังเคลื่อนที่ไปทางขวาด้วยความเร็ว 4 วัตถุ B (มวล 1) กำลังเคลื่อนที่ไปทางซ้ายด้วยความเร็ว 2 — มันกำลังปิดช่องว่างกันด้วยความเร็ว 6 หน่วยต่อวินาที หลังจาก impulse แล้ว A ช้าลงเหลือ 1 และ B เร็วขึ้นเป็น 4 ตอนนี้มันแยกออกจากกันด้วยความเร็ว 3 หน่วยต่อวินาที ซึ่งเท่ากับครึ่งหนึ่งของความเร็วที่เข้าใกล้กันพอดี ตรงกับ restitution 0.5 momentum (mass * velocity รวมทั้งสองวัตถุ) อ่านได้ 6 ทั้งก่อนและหลัง — เป็นกฎฟิสิกส์จริง ๆ ที่ได้มาฟรี ๆ จากสูตร ไม่ใช่สิ่งที่โค้ดบังคับไว้ตรง ๆ ซึ่งเป็นสัญญาณที่ดีว่าสูตรถูกต้อง
impulse ในหัวข้อ 7 แตะแค่ velocity ตามแนว normal เท่านั้น — ไม่ได้พูดถึงการไถลด้านข้างเลย Friction คือ impulse ตัวที่สอง ใช้ตามแนวพื้นผิว (ทิศทาง tangent ซึ่งตั้งฉากกับ normal) ที่ต้านการไถลด้านข้างนั้น friction จริง ๆ เดินตามกฎประมาณที่เรียกว่า Coulomb friction: แรงด้านข้างสูงสุดที่มีได้จะแปรผันตรงกับแรงที่พื้นผิวถูกกดเข้าหากัน จำกัดด้วย friction coefficient (มักเขียนเป็น mu ตัวอักษรกรีก) พูดให้เป็นรูปธรรม: คำนวณ impulse ที่จะยกเลิกความเร็วสัมพัทธ์ด้านข้างได้ เต็มที่ แล้ว clamp ขนาดของมันไว้ที่ mu * normalImpulse — ถ้าการหยุดเต็มที่ต้องการมากกว่านั้น พื้นผิวจะไถล (kinetic friction) ถ้าไม่ ก็จะยึดติดกันสนิท:
#include <iostream>
#include <cmath>
struct Vec2 { float x, y; };
struct RigidBody { Vec2 velocity; float invMass; };
struct Contact { Vec2 normal; };
Vec2 add(Vec2 a, Vec2 b) { return { a.x + b.x, a.y + b.y }; }
Vec2 sub(Vec2 a, Vec2 b) { return { a.x - b.x, a.y - b.y }; }
Vec2 scale(Vec2 a, float s) { return { a.x * s, a.y * s }; }
float dot(Vec2 a, Vec2 b) { return a.x * b.x + a.y * b.y; }
float length(Vec2 v) { return std::sqrt(v.x * v.x + v.y * v.y); }
void applyFriction(RigidBody& a, RigidBody& b, Contact c, float normalImpulse, float friction) {
Vec2 relVel = sub(b.velocity, a.velocity);
// build the tangent: the part of relVel that is NOT along the normal
float velAlongNormal = dot(relVel, c.normal);
Vec2 tangentRaw = sub(relVel, scale(c.normal, velAlongNormal));
float tangentLen = length(tangentRaw);
if (tangentLen < 0.0001f) return; // no sideways motion to resist
Vec2 tangent = scale(tangentRaw, 1.0f / tangentLen);
float totalInvMass = a.invMass + b.invMass;
float velAlongTangent = dot(relVel, tangent);
float jt = -velAlongTangent / totalInvMass; // impulse to fully stop sliding
float maxFriction = friction * normalImpulse; // Coulomb's clamp
if (jt > maxFriction) jt = maxFriction;
if (jt < -maxFriction) jt = -maxFriction;
Vec2 frictionImpulse = scale(tangent, jt);
a.velocity = sub(a.velocity, scale(frictionImpulse, a.invMass));
b.velocity = add(b.velocity, scale(frictionImpulse, b.invMass));
}
int main() {
RigidBody floorBody = { {0.0f, 0.0f}, 0.0f }; // static
RigidBody box = { {3.0f, -2.0f}, 1.0f }; // sliding in and down
Contact c = { {0.0f, 1.0f} }; // normal points up
// step 1: resolve the normal impulse (section 7's formula, restitution 0.2)
Vec2 relVel = sub(box.velocity, floorBody.velocity);
float velAlongNormal = dot(relVel, c.normal);
float totalInvMass = floorBody.invMass + box.invMass;
float j = -(1.0f + 0.2f) * velAlongNormal / totalInvMass;
box.velocity = add(box.velocity, scale(c.normal, j * box.invMass));
std::cout << "after normal impulse: box.velocity=(" << box.velocity.x << ", " << box.velocity.y << ")\n";
// step 2: friction, clamped by mu = 0.5
applyFriction(floorBody, box, c, j, 0.5f);
std::cout << "after friction: box.velocity=(" << box.velocity.x << ", " << box.velocity.y << ")\n";
}ผลลัพธ์:
after normal impulse: box.velocity=(3, 0.4)
after friction: box.velocity=(1.8, 0.4)กล่องกระแทกพื้นด้วยความเร็วไปทางขวา 3 และลง 2 normal impulse (restitution 0.2) หยุดการเคลื่อนที่ลงและทำให้เด้งขึ้นเล็กน้อย 0.4 แต่ไม่แตะความเร็วด้านข้างเต็ม ๆ 3 เลย จากนั้น friction จะพยายามยกเลิกความเร็วด้านข้างนั้นให้เต็มที่ — ซึ่งต้องการ impulse ขนาด 3 — แต่ Coulomb clamp อนุญาตแค่ 0.5 * 2.4 = 1.2 เท่านั้น กล่องเลยยังไถลต่อ แค่ช้าลง: จาก 3 เหลือ 1.8 ถ้า friction coefficient สูงพอที่จะครอบคลุมทั้ง 3 กล่องก็จะหยุดไถลสนิทในสเต็ปเดียวแทน — clamp นี่แหละคือเส้นแบ่งพอดีระหว่าง "ยึดติด" กับ "ไถล"
staticFriction (ต้านการเริ่มไถลจากท่านิ่ง) แยกจาก dynamicFriction (ต้านการไถลที่กำลังเกิดขึ้นอยู่แล้ว) เพราะในความเป็นจริงมักต้องใช้แรงมากกว่าในการเริ่มให้อะไรไถล เทียบกับการทำให้มันไถลต่อไป เวอร์ชันที่ทำให้ง่ายด้วย mu ตัวเดียวข้างบนนี้ก็คือไอเดียหลักตัวเดียวกัน แค่ตัดความแตกต่างนั้นออกเพื่อความง่ายในการอธิบายcontact เดียวแก้ได้เกลี้ยงในรอบเดียว ตามที่หัวข้อ 6 กับ 7 เพิ่งโชว์ไป แต่ฉากจริง ๆ แทบไม่มีแค่ contact เดียวหรอก — กล่องที่วางอยู่บนพื้นแล้วมีอีกกล่องวางซ้อนอยู่ข้างบน มี contact สองอันที่ใช้วัตถุร่วมกันตัวหนึ่ง (กล่องตรงกลาง) และการแก้ contact หนึ่งอาจทำให้อีก contact หนึ่งพังไปแบบเงียบ ๆ ได้ ลองดูให้เห็นภาพ: พื้น static, box2 วางอยู่บนมัน (จมลงไป 0.3), box3 วางอยู่บน box2 (จมลงไป 0.3 เหมือนกัน) แต่ละกล่องสูงพอดี 1 หน่วย:
#include <iostream>
struct Box { float bottom; float invMass; }; // 1D vertical stack; height is always 1
void resolveContact(Box& lower, Box& upper) {
float lowerTop = lower.bottom + 1.0f;
float penetration = lowerTop - upper.bottom;
if (penetration <= 0.0f) return;
float totalInv = lower.invMass + upper.invMass;
if (totalInv <= 0.0f) return;
float shareLower = lower.invMass / totalInv;
float shareUpper = upper.invMass / totalInv;
lower.bottom -= penetration * shareLower;
upper.bottom += penetration * shareUpper;
}
int main() {
Box floorBox = { 0.0f, 0.0f }; // static: top always at 1.0
Box box2 = { 0.7f, 1.0f }; // sunk 0.3 into the floor
Box box3 = { 1.4f, 1.0f }; // sunk 0.3 into box2
for (int iter = 1; iter <= 3; iter++) {
resolveContact(floorBox, box2); // contact A: floor - box2
resolveContact(box2, box3); // contact B: box2 - box3
float penA = (floorBox.bottom + 1.0f) - box2.bottom;
float penB = (box2.bottom + 1.0f) - box3.bottom;
std::cout << "iteration " << iter << ": penetrationA=" << penA
<< " penetrationB=" << penB << "\n";
}
}ผลลัพธ์:
iteration 1: penetrationA=0.3 penetrationB=0
iteration 2: penetrationA=0.15 penetrationB=0
iteration 3: penetrationA=0.075 penetrationB=0ลองไล่ดูว่าเกิดอะไรขึ้นใน iteration 1: การแก้ contact A (floor-box2) ผลัก box2 ขึ้นเต็ม ๆ 0.3 ซึ่งดัน box2 ให้จมลึกเข้าไปใน box3 มากขึ้น — penetration ระหว่าง box2-box3 กระโดดจาก 0.3 เป็น 0.6 จากนั้นการแก้ contact B จะแบ่งช่องว่าง 0.6 นั้น 50/50 ระหว่าง box2 กับ box3 ซึ่งดึง box2 กลับลงมา 0.3 — กลับไปที่จุดเริ่มต้นเป๊ะ ทำให้ contact A เปิดขึ้นมาใหม่ หลังจากผ่านไปหนึ่งรอบเต็ม ๆ contact A ก็พังพอ ๆ กับตอนที่สเต็ปเริ่มต้น ต้องประมวลผลทั้งสอง contact อีกครั้ง (iteration 2) contact A ถึงจะดีขึ้น — penetration ลดลงครึ่งหนึ่งเหลือ 0.15 แล้วลดลงครึ่งอีกทีเหลือ 0.075 ใน iteration 3 และจะลดลงครึ่งไปเรื่อย ๆ แบบนี้ตลอดไปโดยไม่มีทางถึงศูนย์เป๊ะในจำนวนสเต็ปจำกัด
นี่คือเหตุผลที่ physics engine ถูกเรียกว่า iterative solver: การไล่ผ่านทุก contact แค่รอบเดียว (แพทเทิร์นนี้เรียกว่า sequential impulse เมื่อใช้กับ velocity และในทางคณิตศาสตร์มันคือรูปแบบหนึ่งของวิธี Gauss-Seidel) แทบไม่เคยทำให้ทุก contact พอใจพร้อมกันได้เลย ในตอนที่มี contact สองอันขึ้นไปใช้วัตถุร่วมกัน การรันหลาย ๆ รอบ — สัก 4 ถึง 10 iteration โดยทั่วไป — จะทำให้ทั้งระบบลู่เข้าใกล้สถานะที่ "ทุก contact พอใจพร้อมกัน" มากพอที่จะดูและรู้สึกแน่นหนา ยิ่ง iteration เยอะ กองวัตถุก็จะแน่นและแม่นยำขึ้นแต่แลกมาด้วย CPU ที่มากขึ้นต่อสเต็ป iteration น้อยจะถูกกว่าแต่ทำให้กองวัตถุจมลงเห็นได้ชัด สั่น หรือรู้สึกยวบ ๆ นิดหน่อย นี่คือ trade-off จริง ๆ ที่ปรับได้ ไม่ใช่บั๊กที่แก้ครั้งเดียวแล้วจบ — ซึ่งเป็นเหตุผลที่ engine อย่าง Unity เปิดให้ปรับจำนวน iteration เป็นค่าตั้งค่าแทนที่จะฝังตายตัวไว้ในโค้ด
ทุกการเทสต์ที่ผ่านมาเป็นแบบ discrete ทั้งหมด: มันดูแค่ว่าวัตถุอยู่ตรงไหนตอนต้นกับตอนจบของสเต็ป แล้วเช็กว่าสองสแนปช็อตนั้นทับกันไหม วิธีนี้ใช้ได้ดีกับอะไรก็ตามที่เคลื่อนที่ช้ากว่าขนาดตัวมันเองต่อสเต็ป แต่จะพังสนิทกับอะไรก็ตามที่เร็วและบาง — กระสุนปืน มีดที่ถูกขว้าง รถที่วิ่งเร็วชนรั้วลวดตาข่าย — เพราะวัตถุอาจอยู่ฝั่งหนึ่งของกำแพงบาง ๆ เต็ม ๆ ตอนต้นสเต็ป แล้วไปอยู่อีกฝั่งเต็ม ๆ ตอนจบสเต็ป โดยทะลุผ่านตรงกลางไปโดยไม่มีสแนปช็อตไหนจับมันได้ระหว่างทางเลย ความผิดพลาดแบบนี้เรียกว่า tunneling
#include <iostream>
struct Wall { float xStart, xEnd; };
bool discreteHitsWall(float x, Wall w) {
return x >= w.xStart && x <= w.xEnd;
}
int main() {
Wall wall = { 5.0f, 5.1f }; // a thin wall, only 0.1 units thick
float x0 = 4.5f; // position at the START of this step
float speed = 80.0f;
float dt = 0.02f;
float x1 = x0 + speed * dt; // position at the END of this step
std::cout << "x0=" << x0 << " x1=" << x1 << "\n";
std::cout << "discrete check at x0: " << (discreteHitsWall(x0, wall) ? "HIT" : "clear") << "\n";
std::cout << "discrete check at x1: " << (discreteHitsWall(x1, wall) ? "HIT" : "clear") << "\n";
}ผลลัพธ์:
x0=4.5 x1=6.1
discrete check at x0: clear
discrete check at x1: clearลูกบอลเคลื่อนที่ไป 1.6 หน่วยในสเต็ปนี้ (80 * 0.02) แต่กำแพงหนาแค่ 0.1 หน่วย — ลูกบอลกระโดดข้ามความกว้างทั้งหมดของกำแพงไปเลยระหว่างแซมเปิลหนึ่งกับอีกแซมเปิลหนึ่ง แล้วการเช็กแบบ discrete ทั้งสองครั้งก็รายงานว่า "clear" ลูกบอล tunnel ทะลุไปโดยไม่มีใครสังเกตเห็นเลย
Continuous collision detection (CCD) แก้ปัญหานี้ด้วยการเทสต์เส้นทางทั้งหมดที่วัตถุกวาดผ่านไปในสเต็ปนี้ ไม่ใช่แค่จุดปลายสองจุด — เรียกว่า swept test เวอร์ชันง่ายที่สุดจะมองการเคลื่อนที่เป็นเส้นตรงจากตำแหน่งเริ่มต้นไปยังตำแหน่งสุดท้าย แล้วหาว่าเส้นนั้นตัดกำแพงเป็นครั้งแรกตรงไหน:
#include <iostream>
struct Wall { float xStart, xEnd; };
bool sweptHitsWall(float x0, float x1, Wall w, float& hitT, float& hitX) {
if (x0 >= w.xStart) return false; // already at/past the wall's near face
if (x1 < w.xStart) return false; // never reached the wall at all
hitT = (w.xStart - x0) / (x1 - x0);
hitX = x0 + hitT * (x1 - x0);
return true;
}
int main() {
Wall wall = { 5.0f, 5.1f };
float x0 = 4.5f, x1 = 6.1f;
float hitT, hitX;
if (sweptHitsWall(x0, x1, wall, hitT, hitX)) {
std::cout << "swept test: HIT at t=" << hitT << ", x=" << hitX << "\n";
} else {
std::cout << "swept test: clear\n";
}
}ผลลัพธ์:
swept test: HIT at t=0.3125, x=5swept test จับสิ่งที่สแนปช็อต discrete ทั้งสองตัวพลาดไปได้พอดี: ลูกบอลตัดผ่านหน้ากำแพงด้านใกล้ที่ 31.25% ของระยะการเคลื่อนที่ในสเต็ปนี้ ที่ตำแหน่ง x = 5.0 engine จริง ๆ จะตอบสนองด้วยการหยุดวัตถุตรงนั้นเลย (หรือแก้ collision เต็มรูปแบบตรงจุดนั้น) แทนที่จะปล่อยให้มันไปต่อถึง x = 6.1 การ implement CCD จริง ๆ ไปไกลกว่าการเทสต์เส้นตรงง่าย ๆ แบบนี้อีก — conservative advancement จะขยาย swept segment ด้วยรัศมีของรูปทรงที่เคลื่อนที่เอง เพื่อให้ทรงกลมความเร็วสูง ไม่ใช่แค่จุดความเร็วสูง ถูกจับได้อย่างถูกต้อง และบาง engine ก็ใช้ speculative contacts แทน ซึ่งจะทำนาย collision ที่กำลังจะเกิดล่วงหน้าหนึ่งสเต็ป แล้วเริ่มแก้มันตั้งแต่ก่อนที่รูปทรงจะแตะกันจริง ๆ ไม่ว่าจะแบบไหน ไอเดียหลักก็เหมือนเดิม: มองเส้นทางทั้งหมดที่กวาดผ่านในสเต็ปนี้ ไม่ใช่แค่จุดปลายสองจุด
กล่องที่หล่นลงพื้นควรจะ — หลังจากเด้งหนึ่งหรือสองครั้ง — วางนิ่งอยู่ตรงนั้น นิ่งสนิท ไม่กินทรัพยากรอะไรเลย การทำให้วัตถุมาหยุดนิ่งจริง ๆเป็นปัญหาเล็ก ๆ ของตัวมันเอง เพราะ loop แบบง่าย ๆ จากหัวข้อก่อน ๆ ไม่เคยยอมให้มันนิ่งได้จริง แรงโน้มถ่วงดึงกล่องที่วางนิ่งลงมานิดเดียวทุกสเต็ป contact ก็รับมันไว้ แล้ว restitution ในหัวข้อ 7 ก็เด้งเจ้านิดเดียวนั้นกลับขึ้นไปอย่างซื่อสัตย์ ผลลัพธ์คือวัตถุที่ไม่เคยหยุดกระตุกเลย:
#include <iostream>
// one step of a ball above a floor at y=0: semi-implicit Euler, then a bounce.
// `threshold`: approach speeds slower than this get restitution 0 (no bounce).
void step(float& y, float& v, float g, float dt, float e, float threshold) {
v += g * dt;
y += v * dt;
if (y < 0.0f) { // hit the floor this step
y = 0.0f;
if (v < 0.0f) {
float useE = (-v < threshold) ? 0.0f : e; // slow contact -> no bounce
v = -useE * v;
}
}
}
void run(const char* label, float threshold) {
float y = 0.0f, v = 0.0f, g = -10.0f, dt = 1.0f / 60.0f, e = 0.5f;
std::cout << label << "\n";
for (int i = 1; i <= 6; i++) {
step(y, v, g, dt, e, threshold);
std::cout << " step " << i << ": y=" << y << " v=" << v << "\n";
}
}
int main() {
run("no threshold (restitution always on) -> never settles:", 0.0f);
run("with threshold 1.0 -> settles to rest:", 1.0f);
}ผลลัพธ์:
no threshold (restitution always on) -> never settles:
step 1: y=0 v=0.0833333
step 2: y=0 v=0.0416667
step 3: y=0 v=0.0625
step 4: y=0 v=0.0520833
step 5: y=0 v=0.0572917
step 6: y=0 v=0.0546875
with threshold 1.0 -> settles to rest:
step 1: y=0 v=0
step 2: y=0 v=0
step 3: y=0 v=0
step 4: y=0 v=0
step 5: y=0 v=0
step 6: y=0 v=0เมื่อ restitution เปิดอยู่ตลอด (ด้านบน) velocity ไม่เคยถึงศูนย์ — มันล็อกอยู่ที่การสั่นถาวรรอบ ๆ +0.05 เพราะทุกสเต็ปแรงโน้มถ่วงป้อนความเร็วลงเข้ามานิดหน่อย แล้ว restitution ก็เตะมันกลับขึ้นไปทันที ในฉากจริง ๆ ความเร็วขึ้นที่เหลืออยู่นั้นจะยกกล่องลอยพ้นพื้นนิดหนึ่งทุกสเต็ป: นี่คือ resting jitter ที่เห็นได้ชัดของกองวัตถุที่ไม่ยอมนิ่ง วิธีแก้ (ด้านล่าง) คือ restitution threshold (เรียกอีกอย่างว่า bounce threshold): ถ้าความเร็วที่เข้าใกล้กันต่ำกว่าค่าเล็ก ๆ ค่าหนึ่ง ให้ถือว่า restitution เป็นศูนย์ contact ที่ช้าเลยแค่หยุดแทนที่จะเด้ง บรรทัดเดียว แล้วกล่องก็นิ่งตั้งแต่สเต็ปแรกและอยู่นิ่งเลย ทุก engine ทำแบบนี้ Unity เปิดให้ปรับผ่าน Physics.bounceThreshold
threshold ทำให้วัตถุมานิ่งได้ แต่วัตถุที่นิ่งก็ยังกิน CPU อยู่ดี: ทุกสเต็ปมันถูก integrate ถูกเทสต์การชน และถูกส่งผ่าน solver แค่เพื่อจะไปจบที่เดิมเป๊ะ ในฉากที่มีวัตถุนิ่งอยู่เป็นร้อย ๆ ชิ้น — โกดังที่เต็มไปด้วยลัง กองเศษซากปรักหักพัง — นั่นคือความสิ้นเปลืองมหาศาล และแย่กว่านั้น numerical noise เล็ก ๆ ที่เหลืออยู่ใน solver ก็คอยสะกิดวัตถุพวกนั้นจนมันไม่เคยนิ่งสนิทเลย คำตอบคือ sleeping: พอความเร็วของวัตถุอยู่ต่ำกว่า threshold นานพอ ให้ทำเครื่องหมายว่ามันหลับ แล้วหยุด simulate มันไปเลย จนกว่าจะมีอะไรมาแตะแล้วปลุกมันขึ้นมา
#include <iostream>
#include <cmath>
struct Body {
float x, v;
float idleTime; // how long we have been below the sleep threshold
bool asleep;
};
int main() {
Body b = {0.0f, 5.0f, 0.0f, false};
float dt = 1.0f / 60.0f;
float sleepSpeed = 0.10f; // below this counts as "still"
float sleepDelay = 0.05f; // must stay still this long to fall asleep
for (int i = 1; i <= 14; i++) {
if (!b.asleep) {
b.v *= 0.6f; // friction bleeds speed off each step
b.x += b.v * dt;
if (std::fabs(b.v) < sleepSpeed) b.idleTime += dt;
else b.idleTime = 0.0f;
if (b.idleTime >= sleepDelay) b.asleep = true; // go to sleep
}
// an asleep body skips integration entirely -- zero CPU until woken
std::cout << "step " << i << ": x=" << b.x << " v=" << b.v
<< (b.asleep ? " [ASLEEP]" : "") << "\n";
}
}ผลลัพธ์:
step 1: x=0.05 v=3
step 2: x=0.08 v=1.8
step 3: x=0.098 v=1.08
step 4: x=0.1088 v=0.648
step 5: x=0.11528 v=0.3888
step 6: x=0.119168 v=0.23328
step 7: x=0.121501 v=0.139968
step 8: x=0.1229 v=0.0839808
step 9: x=0.12374 v=0.0503885
step 10: x=0.124244 v=0.0302331 [ASLEEP]
step 11: x=0.124244 v=0.0302331 [ASLEEP]
step 12: x=0.124244 v=0.0302331 [ASLEEP]
step 13: x=0.124244 v=0.0302331 [ASLEEP]
step 14: x=0.124244 v=0.0302331 [ASLEEP]วัตถุค่อย ๆ ช้าลงด้วย friction จนกระทั่งที่สเต็ป 10 ความเร็วของมันอยู่ต่ำกว่า 0.1 ครบ 0.05 วินาทีตามที่กำหนด มันก็หลับ แล้วจากนั้นตำแหน่งก็ถูกแช่แข็งและ loop ก็ไม่ทำงานอะไรกับมันอีกเลย engine จริง ๆ ไปไกลกว่านั้นอีกขั้นด้วยการทำให้ทั้ง island หลับไปพร้อมกัน — กองกล่องที่วางพิงกันอยู่ต้องหลับและตื่นเป็นกลุ่ม เพราะการปลุกกล่องล่างสุดต้องปลุกกล่องที่อยู่ข้างบนมันด้วย นี่ยังเป็นแหล่งกำเนิด "บั๊กที่ไม่ใช่บั๊ก" ที่โด่งดัง: ลังที่วางอยู่บนแป้นกดที่ดีไซเนอร์คาดว่าจะทำงาน หรือวัตถุที่ควรตอบสนองต่อการดันช้า ๆ แต่กลับไม่ทำอะไรเลย — เพราะมันหลับอยู่ ใน Unity ปุ่มปรับคือ Rigidbody.sleepThreshold คำสั่งเช็กคือ rb.IsSleeping() และทางออกฉุกเฉินคือ rb.WakeUp() เพื่อบังคับปลุกวัตถุตอนที่โค้ดของเราเองเปลี่ยนอะไรบางอย่างที่ physics engine มองไม่เห็น
ทุกไอเดียในหัวข้อ 1 ถึง 11 คือสิ่งที่ physics ในตัวของ Unity (ขับเคลื่อนด้วย PhysX) ทำให้เราอยู่แล้วในทุก fixed step โดยที่เราไม่ต้องเขียนโค้ด broad phase, narrow phase หรือ impulse-resolution แม้แต่บรรทัดเดียว การมองออกว่าคอนเซปต์ไหนซ่อนอยู่ใน component ตัวไหนคือประเด็นทั้งหมดของหัวข้อนี้ — เราจะใช้เวลาไปกับการตั้งค่า Rigidbody กับ Collider มากกว่าการเขียน solver เองเยอะมาก
using UnityEngine;
public class BouncyBall : MonoBehaviour
{
void Start()
{
Rigidbody rb = GetComponent<Rigidbody>();
rb.mass = 2f; // section 1's RigidBody.mass
rb.linearVelocity = new Vector3(3f, 0f, 0f); // section 1's RigidBody.velocity
rb.useGravity = true; // section 2's integrate(), running every FixedUpdate
// section 10: switch on CCD only for objects at real risk of tunneling
rb.collisionDetectionMode = CollisionDetectionMode.ContinuousDynamic;
}
}Rigidbody ก็คือ struct RigidBody ของหัวข้อ 1 นั่นแหละ ตรงกันทุก field: mass, linearVelocity, angularVelocity, position/rotation (ตัวหลังเก็บเป็น quaternion ตามที่บทเรื่อง rotation เคยอธิบายไว้) ทุก tick ของ FixedUpdate — fixed timestep ตัวเดียวกับที่บท game loop พูดถึง — Unity จะรัน semi-implicit Euler step จากหัวข้อ 2 ให้เองภายใน โดยใช้ Time.fixedDeltaTime ด้วยเหตุผลเรื่อง stability กับ determinism แบบเดียวกับที่หัวข้อ 2 อธิบายไปเป๊ะ ๆ
A Collider (BoxCollider, SphereCollider, CapsuleCollider, MeshCollider สำหรับ convex mesh) คือรูปทรงที่ broad phase กับ narrow phase ของ PhysX ใช้ทำงานจริง ๆ — PhysX รัน broad phase ของตัวเองบน world-space bounds ของทุก collider (งานของหัวข้อ 3) และ narrow phase ของตัวเอง (การเทสต์แบบ SAT สำหรับรูปทรงง่าย ๆ และอัลกอริทึมตระกูล GJK สำหรับ convex mesh ตรงกับหัวข้อ 4 เป๊ะ) เพื่อหา contact ทั้งหมดนี้เกิดขึ้นเบื้องหลังโดยที่โค้ดของเราไม่เห็นเลย:
using UnityEngine;
public class ContactLogger : MonoBehaviour
{
void OnCollisionEnter(Collision collision)
{
// section 5's Contact struct, handed to you directly by Unity
foreach (ContactPoint contact in collision.contacts)
{
Debug.Log("point=" + contact.point + " normal=" + contact.normal
+ " separation=" + contact.separation); // separation < 0 means penetrating
}
}
}ทุก field ของ struct Contact จากหัวข้อ 5 มีอยู่ใน ContactPoint ครบเลย: point, normal, และ separation (ชื่อที่ Unity ใช้เรียก penetration depth แต่กลับเครื่องหมาย — ค่าติดลบหมายถึงทับกันอยู่) ส่วน restitution กับ friction — หัวข้อ 7 กับ 8 — อยู่ใน Physic Material asset ที่ผูกกับ collider พร้อม combine mode (average, minimum, maximum, หรือ multiply) ที่คุมว่า material ของวัตถุสองตัวที่ต่างกันจะผสมกันยังไงตอนแตะกัน:
using UnityEngine;
public class BouncyMaterialSetup : MonoBehaviour
{
void Start()
{
PhysicsMaterial mat = new PhysicsMaterial();
mat.bounciness = 0.5f; // section 7's restitution
mat.dynamicFriction = 0.5f; // section 8's friction (sliding)
mat.staticFriction = 0.6f; // section 8's friction (starting from rest)
mat.frictionCombine = PhysicsMaterialCombine.Average;
GetComponent<Collider>().material = mat;
}
}จำนวน iteration จากหัวข้อ 9 ก็เป็นค่าตั้งค่าจริง ๆ ที่เปิดให้ปรับได้เหมือนกัน — Project Settings > Physics > Solver Iterations (position iterations) และ Solver Velocity Iterations คุมพอดีว่า solver ภายในของ PhysX จะรันกี่รอบต่อสเต็ป เป็นไอเดียเดียวกับการลู่เข้าด้วยการทำซ้ำจากตัวอย่างกล่องซ้อนกันในหัวข้อ 9 แค่รันกับทุก contact ในฉากพร้อมกัน แทนที่จะทำแค่สอง contact ด้วยมือ
invMass = 0) ไม่ขยับเลย เช่น พื้นหรือกำแพง1/mass ใช้เพื่อให้ "mass เป็นอนันต์" (ไม่ขยับเลย) แทนได้ด้วยค่า 0 แทนที่จะเป็นการหารที่ไม่มีนิยามdt เดียวกันทุกสเต็ป เพื่อความเสถียร determinism และความเสี่ยง tunneling ที่คาดเดาได้dt เดียวกันเสมอไม่ว่า frame rate จะสั่นแค่ไหนbeta/dt * (penetration - slop)) เข้าไปใน velocity solver แทนที่จะขยับ position ตรง ๆAABBs:
E: min=(1, 1) max=(3, 3)
F: min=(2.2, 1) max=(4.2, 3)
G: min=(9, 9) max=(11, 11)
E-F: overlap (x: 1<=4.2 and 3>=2.2; y: identical range) -> CANDIDATE
E-G: no overlap (E.max.x=3 < G.min.x=9)
F-G: no overlap (F.max.x=4.2 < G.min.x=9)
narrow phase E-F:
distance = 1.2, sum of radii = 2 -> 1.2 < 2, colliding
penetration = 2 - 1.2 = 0.8
normal = (1, 0)
point on E = (3, 2), point on F = (2.2, 2), contact = midpoint = (2.6, 2)มีแค่ E-F เท่านั้นที่ AABB ทับกัน มันเลยเป็นคู่เดียวที่ narrow phase ต้องเทสต์ — G ไม่เคยถูกรัน exact test กับใครเลย ซึ่งนั่นคือประเด็นทั้งหมดของ broad phase การเทสต์แบบละเอียดยืนยันว่าชนกันจริง โดยมี penetration ค่อนข้างมาก (0.8 ซึ่งเกือบเท่ารัศมีของแต่ละวงกลมเลย) มี normal ที่ชี้ตรงจาก E ไปทาง F ตามแนวแกน x (เพราะวงกลมทั้งสองอยู่ที่ y เดียวกัน) และมี contact point อยู่ตรงกึ่งกลางของช่วงที่ทับกัน
j และ velocity ของวัตถุทั้งสองหลัง impulse แล้วตรวจสอบว่า momentum อนุรักษ์อยู่จริงinvMassA = 1/3, invMassB = 1
relVel = velB - velA = (-4,0) - (2,0) = (-6, 0)
velocityAlongNormal = -6
j = -(1 + 1) * (-6) / (1/3 + 1) = 12 / (4/3) = 9
impulse = (9, 0)
velA_new = (2,0) - (9,0)*(1/3) = (2-3, 0) = (-1, 0)
velB_new = (-4,0) + (9,0)*1 = (-4+9, 0) = (5, 0)
momentum before = 3*2 + 1*(-4) = 6 - 4 = 2
momentum after = 3*(-1) + 1*5 = -3 + 5 = 2 -- matchesด้วย restitution 1 วัตถุจะแยกออกจากกันด้วยความเร็วเท่ากับตอนที่เข้าใกล้กันพอดี: มันกำลังปิดช่องว่างกันด้วยความเร็ว 6 หน่วยต่อวินาที และหลัง impulse มันแยกออกจากกันด้วย velB_new - velA_new = 5 - (-1) = 6 หน่วยต่อวินาที — คุณสมบัติที่นิยามการชนแบบ elastic สมบูรณ์เลย momentum ที่อ่านได้ 2 พอดีทั้งก่อนและหลังคือการเช็กความถูกต้องแบบฟรี ๆ เดียวกับที่หัวข้อ 7 ใช้ และมันก็ยังใช้ได้ที่นี่เหมือนกัน
x = 9.0 เคลื่อนที่ด้วยความเร็ว 150 หน่วยต่อวินาที ด้วย fixed timestep 0.02 วินาที กำแพงบาง ๆ อยู่ระหว่าง x = 10.0 ถึง x = 10.2 การเช็ก collision แบบ discrete (แค่ต้น/ปลาย) จะจับการชนในสเต็ปนี้ได้ไหม ถ้าจับไม่ได้ ให้ใช้สูตร swept-test จากหัวข้อ 10 หาว่าวัตถุไปถึงกำแพงเป็นครั้งแรกตรงไหนของการเคลื่อนที่travel this step = 150 * 0.02 = 3.0
x0 = 9.0, x1 = 9.0 + 3.0 = 12.0
discrete check at x0=9.0: 10.0 <= 9.0 <= 10.2 ? NO
discrete check at x1=12.0: 10.0 <= 12.0 <= 10.2 ? NO
-> tunneling: the discrete check misses the wall completely
swept test:
hitT = (10.0 - 9.0) / (12.0 - 9.0) = 1.0 / 3.0 = 0.3333
hitX = 9.0 + 0.3333 * 3.0 = 10.0ระยะที่วัตถุเคลื่อนที่ต่อสเต็ป (3.0 หน่วย) มากกว่าความหนาของกำแพง (0.2 หน่วย) ถึงสิบห้าเท่า สแนปช็อต discrete ทั้งสองเลยตกลงไปคนละฝั่งของกำแพงได้ง่าย ๆ โดยไม่เคยทับมันเลย — ตรงกับอาการ tunneling จากหัวข้อ 10 เป๊ะ swept test พบว่าวัตถุไปถึงหน้ากำแพงด้านใกล้จริง ๆ ที่หนึ่งในสามของการเคลื่อนที่ในสเต็ปนี้ ที่ x = 10.0 engine จริง ๆ จะหยุดหรือแก้ collision ตรงนั้น แทนที่จะปล่อยให้วัตถุไปถึง x = 12.0 เหมือนไม่มีอะไรขวางอยู่เลย
closest point = clamp center (6,4) into [2,5] x [1,3]
= (5, 3) -- a corner of the box
d = center - closest = (6-5, 4-3) = (1, 1)
dist = sqrt(1*1 + 1*1) = sqrt(2) = 1.41421
1.41421 < 1.5 -> COLLIDING
penetration = 1.5 - 1.41421 = 0.08579
normal = d / dist = (0.7071, 0.7071)ศูนย์กลางวงกลม clamp ไปที่ (5, 3) ซึ่งเป็นมุมของกล่อง ฟีเจอร์ที่ใกล้ที่สุดเลยเป็นจุดเดียว ไม่ใช่หน้า — และนั่นคือเหตุผลที่ normal ออกมาเป็นแนวทแยง (0.7071, 0.7071) แทนที่จะขนานแกน วงกลมทับมุมนั้นอยู่แค่ตื้น ๆ 0.086 และทิศทางการผลักออกก็ชี้ออกจากมุมอย่างสวยงามตามแนวเส้น 45 องศา เคสมุมนี้แหละคือเหตุผลทั้งหมดที่ closest-point test ชนะการเช็กแบบ "อยู่ฝั่งไหน" ที่ไร้เดียงสา: มันจัดการทั้งหน้าและมุมด้วย clamp ครั้งเดียวเหมือนกัน
นี่คือ pipeline เต็มรูปแบบที่ physics engine รันทุก fixed step: integrate position กับ rotation ของวัตถุทุกชิ้นไปข้างหน้าด้วย semi-implicit Euler ย่อโลกที่เต็มไปด้วยคู่ที่เป็นไปได้ให้เหลือแค่ candidate จริง ๆ ไม่กี่คู่ด้วย broad phase ยืนยันและวัดแต่ละคู่อย่างละเอียดด้วย narrow phase แปลงการทับกันที่ยืนยันแล้วให้กลายเป็น contact แล้วแก้ contact นั้นสองรอบ — รอบหนึ่งแก้ position (ดัน penetration ออก) และอีกรอบแก้ velocity (impulse, restitution, และ friction) — ทำซ้ำขั้นตอนแก้ทั้งหมดนี้สักไม่กี่รอบ เพื่อให้กองและกลุ่มของ contact ลงตัวเป็นสถานะที่สอดคล้องกัน แทนที่จะสู้กันเอง และสุดท้ายปล่อยให้วัตถุที่นิ่งสนิทแล้ว sleep เพื่อให้ฉากที่สงบแทบไม่กินทรัพยากรเลย Continuous collision detection คือสิ่งเดียวที่ต้องเพิ่มเข้ามาเพื่อไม่ให้ pipeline นี้พลาดอะไรก็ตามที่เร็วและบางไป Rigidbody กับ Collider ของ Unity รันทุกขั้นตอนนี้ให้เราเอง ซึ่งเป็นเหตุผลที่การเข้าใจมันในบทนี้เปลี่ยน "physics มันพังว่ะ" ให้กลายเป็น checklist สั้น ๆ เจาะจง แทนที่จะเป็นการเดาสุ่ม