This is a follow-up for convex hull collisions in 3d.
This article extends the separating axis collision implementation from that previous article to calculate the contact manifold using polygon clipping on the collision plane normal to the separating axis.
Read on or jump to the final result: hull-manifold.js
Or jump ahead to the interactive demo.
The contact manifold is the set of points shared between two colliding polyhedra once the penetration along the minimum translation vector is applied to separate them.
The contact manifold is necessary to generate physically stable results. In the previous implementation, the implementation produced only a single contact point even when two objects may have a line or polygon area in common. If only one of those points is used to apply a force to each object, the result will not look physically plausible.
Luckily, for the case of convex polyhedra colliding, the convex manifold will also be convex. We also do not need the edges of the polygon, only the points which comprise a 2d convex hull. This allows for a simpler and faster implementation compared to a more general polygon clipping algorithm.
The following sections will be additional or ammended methods in the Hull class from
the previous article.
First we need to update the collide() function to write multiple points to the output.
The collision type E represents an edge/edge collision.
For these, we need only use the existing lineLine() method which results in a single
collision point.
We also need to ensure that this edge/edge intersection has a result with mu0
mu1 in the bounds of 0 to 1 to use this
candidate separation.
The collision types A and B are face/face collisions and we will run our
clipping algorithm defined later in the clip() method on these collision types.
For face/face collision types A and B, we will also need to calculate the
incident face in a later findIncident() method.
static TYPE_A = 0; static TYPE_B = 1; static TYPE_E = 2
static #qA = { type: this.TYPE_A, face: -1, support: -1, separation: -Infinity }
static #qB = { type: this.TYPE_B, face: -1, support: -1, separation: -Infinity }
static #qE = { type: this.TYPE_E, index: [-1,-1], separation: -Infinity, normal: [0,0,0] }
static #mu = [0,0]
static collide(out, A, B) {
out.length = 0
if (!this.queryFaceDirections(this.#qA, A, B)) return null
if (!this.queryFaceDirections(this.#qB, B, A)) return null
if (!this.queryEdgeDirections(this.#qE, A, B)) return null
out.info = this.#qA
if (this.#qB.separation > out.info.separation) {
out.info = this.#qB
}
if (this.#qE.separation > out.info.separation) {
let ea = A.edges[this.#qE.index[0]]
let eb = B.edges[this.#qE.index[1]]
let a0 = A.positions[ea[0]]
let a1 = A.positions[ea[1]]
let b0 = B.positions[eb[0]]
let b1 = B.positions[eb[1]]
vec3.scale(out.mtv, this.#qE.normal, -this.#qE.separation)
this.#lineLine(Hull.#mu, Hull.#tv0, Hull.#tv1, a0, a1, b0, b1)
if (Math.max(Math.abs(Hull.#mu[0]*2.0-1.0), Math.abs(Hull.#mu[1]*2.0-1.0)) <= 1.0) {
out.info = this.#qE
}
}
if (out.info.type === this.TYPE_A) {
let faceA = this.#qA.face
let N = A.normals[faceA]
let faceB = this.#findIncident(B, N)
vec3.scale(out.mtv, N, -this.#qA.separation)
Hull.#clip(out, N, A, B, faceA, faceB, 1e-8)
} else if (out.info.type === this.TYPE_B) {
let faceB = this.#qB.face
let N = B.normals[faceB]
let faceA = this.#findIncident(A, N)
vec3.scale(out.mtv, N, this.#qB.separation)
Hull.#clip(out, N, B, A, faceB, faceA, 1e-8)
} else if (out.info.type === this.TYPE_E) {
vec3.copy(out.points[out.length++], Hull.#tv0)
}
return out
}
In a face/face collision, the reference face is perpendicular to the separating axis while the incident face does not necessarily have this characteristic. See the diagram on page 91 of Robust Contact Creation for Physics Simulations.
The incident face is the "most anti-parallel face on the other hull". This face's normal will have the minimum dot product with separating axis, which we will return the index of in our implementation:
static #findIncident(X, N) {
let min = Infinity, minI = -1
for (let i = 0; i < X.normals.length; i++) {
let fN = X.normals[i]
let d = vec3.dot(fN, N)
if (d < min) {
min = d
minI = i
}
}
return minI
}
The clip() implementation also needs a point-in-polygon test function.
This point-in-polygon implementation tests if point P projected onto the separating
axis plane defined by the point P itself and the separating axis (face normal)
N is contained inside but not on the edge of face, an array of indices
into the positions array.
Instead of projecting the clip geometry into 2d and potentially allocating more memory to do a more conventional ray-casting point-in-polygon test, here we use the property that the faces are necessarily convex because the colliding hulls are convex to do a different technique.
Each edge in face is crossed with the separating axis N, which produces a
vector that points either towards or away from the interior of the face. So long as the face is an
adjacent set of points along the exterior of the shape, all of these vectors will point inward or
outward. We can describe an edge plane using one of the edge points and each inward/outward vector.
The point P is inside the face when P is on the same side of
all edge planes. We can use the formula
from the previous article here to check which side of the edge plane the point P is on.
The side of the plane that point P is on when the plane is described by a point on the
plane Q with a plane normal of N is given by:
vec3.dot(N, P) - vec3.dot(N, Q)and so our implementation becomes:
static #fN = [0,0,0]
static #pip(P, N, positions, face, EPSILON=1e-3) {
// assumes face is convex with an adjacent list of points
let sides = 0
for (let i = 0; i < face.length; i++) {
let x0 = positions[face[i]]
let x1 = positions[face[(i+1)%face.length]]
vec3.subtract(this.#fN, x0, x1)
vec3.cross(this.#fN, this.#fN, N)
vec3.normalize(this.#fN, this.#fN)
// (fN,x0) represents an edge plane to check for sides
let side = vec3.dot(this.#fN, x0) - vec3.dot(this.#fN, P)
sides += Math.abs(side) < EPSILON ? 0 : Math.sign(side)
}
// checking for +/- N makes this check winding-insensitive
return Math.abs(sides) === face.length
}
Now finally we can put everything together to implement clip().
This static method adds points from the clipped geometry to the contact manifold. This clipped
geometry is itself convex and is the intersection of the incident and reference faces, when certain
conditions are met.
The are 3 cases when we add points which are shown by the inline comments below:
static #clip(out, aN, A, B, iFA, iFB, EPSILON=1e-3) {
out.length = 0
let bN = B.normals[iFB]
let faceA = A.faces[iFA]
let faceB = B.faces[iFB]
let aQ = A.positions[faceA[0]]
let dNQ = vec3.dot(aN, aQ)
// add points from incident face (on B) which are on the separation plane
// and inside of the reference face (on A)
for (let i = 0; i < faceB.length; i++) {
let p = B.positions[faceB[i]]
let d = vec3.dot(aN, p) - dNQ
if (Math.abs(d - out.info.separation) < EPSILON) {
if (this.#pip(p, aN, A.positions, faceA, EPSILON)) {
vec3.copy(out.points[out.length++], p)
}
}
}
// add points from reference face (on A) which are inside the incident face (on B)
// only when face normals are nearly exactly flipped
if (vec3.dot(aN, bN) < EPSILON - 1.0) {
for (let i = 0; i < faceA.length; i++) {
let p = A.positions[faceA[i]]
if (this.#pip(p, aN, B.positions, faceB, EPSILON)) {
vec3.subtract(out.points[out.length++], p, out.mtv)
}
}
}
// add all intersections of the edges between the incident (on B) and reference (on B) faces
// which are on the separation plane
for (let i = 0; i < faceA.length; i++) {
let A0 = A.positions[faceA[i]]
let A1 = A.positions[faceA[(i+1)%faceA.length]]
for (let j = 0; j < faceB.length; j++) {
let B0 = B.positions[faceB[j]]
let B1 = B.positions[faceB[(j+1)%faceB.length]]
if (!this.#lineLine(this.#mu, this.#tv0, this.#tv1, A0, A1, B0, B1)) continue
if (Math.abs(this.#mu[0]*2.0-1.0) > 1.0) continue
if (Math.abs(this.#mu[1]*2.0-1.0) > 1.0) continue
let d = vec3.dot(aN, this.#tv1) - dNQ
if (Math.abs(d - out.info.separation) < EPSILON) {
this.#insertPoint(out, this.#tv1, EPSILON)
}
}
}
}
// plus a little helper to avoid adding duplicate points:
static #insertPoint(out, p, EPSILON=1e-3) {
for (let i = 0; i < out.length; i++) {
if (vec3.distance(out.points[i], p) < EPSILON) return
}
vec3.copy(out.points[out.length++], p)
}
Only the last case adding points to the manifold appears to generate duplicate points in practice so
insertPoint() is skipped for the earlier cases.
Check out the interactive demo which features an interactive editor and viewer along with preset geometries.