GJK implementation

Basic use cases

Given two sets of points in nd (possibly lifted from a 2d space):

1. separation problem: are they linearly separable ?
2. distance problem: they are separable, but how far they are from each other ?

Sketch of the algorithm

Let X1 and X2, two sets of points in nd, of size m1 and m2. The set difference is denoted by X = X1 - X2. Then: dist(X1, X2) = dist(X,o), where o is the origin. X1 and X2 are separable iff dist(X,o) > 0. To check if dist(X,o) > 0, the idea is to generate a sequence of polytope Vk contained in conv(X) (having at most n+1 affinely independent points), such that dist(Vk,o) converges towards (and reaches) dist(X,o):

for problem 2.

1. k = 0, Vk has v vertices, which are affinely independent points of X (1 <= v <= n+1, but v = 1 is the easiest init).
2. if conv(Vk) contains o, then return '0', else go to 3.
3. find Vk* in Vk, the closest polytope to o (of at most n affinely independent points) and p* in Vk*, the closest point to o.
4. set m to farthest_X( normal( Vk* ) ).
5. if m belongs to Vk* , then return 'dist(p* ,o)', else set Vk+1 to Vk* U m, k++, go to 2.

NB: in 4., farthest_X( normal( Vk* ) ) returns one of the points in X, which are the farthest along the direction normal to Vk* and towards o. Note that this can be computed in O(m1+m2) from points of X1 and X2 by inner products.

NB: 2. (resp. 3.) requires the computation of the projection of o onto Vk (resp. Vk* ). This is computed by solving a linear system ( p.M = 0, where matrix M is given by Vk (resp. Vk* ) )

required operations

• inner products for 4.
• linear system solver for 2. and 3.

problem

• round-offs errors or fractions with great denominators, due to the computation of the projection in 2. and 3.

for problem 1.

1. k = 0, Vk has n+1 vertices.
2. if conv(Vk) contains o, then return 'not separable', else go to 3.
3. find Vk* in Vk, one of the polytopes (of n affinely independent points) that separates X from o.
4. set m to farthest_X( normal( Vk* ) ).
5. if m belongs to Vk* , then return 'separable', else set Vk+1 to Vk* U m, k++, go to 2.

NB: 4. can be done with determinant computations (even with determinant sign evaluations) using Vk* instead of its normal.

NB: 2. can be done with a simple orientation predicate (determinant sign evaluations).

NB: 3. can be performed by applying recursively the whole algorithm for X = Vk.

required operations

• determinant sign evaluation.

problems

Determinant sign evaluation fails in degenerate cases (less than n affinely independant points). As a consequence:

• Vk must be initialized with n+1 affinely.
• The final polytope Vk* has n affinely independent points, even if the closest solution is a polytope of less than n affinely independant points.
  Morevoer, the orientation of the n points of Vk* must be always the same.

on-line version

linear algebra libs