Computing Distance
Erin Catto
    Blizzard Entertainment
Convex polygons
Closest points
Overlap
Goal
   Compute the distance between convex
    polygons
Keep in mind
   2D
   Code not optimized
Approach
    simple   complex
Geometry
If all else fails …
DEMO!
Outline
1.   Point to line segment
2.   Point to triangle
3.   Point to convex polygon
4.   Convex polygon to convex polygon
Concepts
1.   Voronoi regions
2.   Barycentric coordinates
3.   GJK distance algorithm
4.   Minkowski difference
Section 1
Point to Line Segment
A line segment
          A      B
Query point
          A       B
Closest point
           A        P   B
Projection: region A
           A           B
Projection: region AB
           A            B
Projection: region B
           A           B
Voronoi regions
     region A       region AB       region B
                A               B
Barycentric coordinates
          A        G         B
              G(u,v)=uA+vB
                 u+v=1
Fractional lengths
               v=0.5           u=0.5
           A            G              B
                G(u,v)=uA+vB
                       u+v=1
Fractional lengths
               v=0.25     u=0.75
           A        G              B
                    G(u,v)=uA+vB
                        u+v=1
Fractional lengths
                  u=1.25
        v=-0.25
       G     A                   B
                  G(u,v)=uA+vB
                     u+v=1
Unit vector
              A      n     B
                     B-A
                  n=
                     B-A
(u,v) from G
       A                  G                 B
                   v              u
           v=
               G-A  n
                              u=
                                  B-G n
                B-A                B-A
(u,v) from Q
       A                      G              B
                   v                   u
           v=
               Q-A  n
                                  u=
                                      B-Q  n
                B-A                    B-A
Voronoi region from (u,v)
                v>0        u>0
        A              G         B
     u > 0 and v > 0         region AB
Voronoi region from (u,v)
                 u>0
       v<0
   G         A              B
        v <= 0         region A
Voronoi region from (u,v)
                                u<0
                  v>0
       A                    B         G
     u <= 0             region B
Closet point algorithm
       input: A, B, Q
       compute u and v
       if (u <= 0)
         P = B
       else if (v <= 0)
         P = A
       else
         P = u*A + v*B
Section 2
Point to Triangle
Triangle
           C
Closest feature: vertex
   Q
       P=B
             C
Closest feature: edge
     Q       P
             C
Closest feature: interior
                 P=Q
             C
Voronoi regions
  Region B          B
                               Region AB
                                                   A
                                                       Region A
                            Region ABC
        Region BC
                                           Region CA
                        C
                        Region C
3 line segments
                                                Q
                   uAB ,v AB 
     uBC ,vBC    C
                                  uCA ,vCA 
Vertex regions
  uAB  0
  vBC  0   B
                                A     uCA  0
                                      v AB  0
                uBC  0    Using line segment uv’s
                v CA  0
Edge regions
                   uAB  0
                   v AB  0
         B
                   ?
                                   A
                              Line segment uv’s are
                              not sufficient
Interior region
                      Line segment uv’s don’t
                      help at all
Triangle barycentric coordinates
       B
                        Q  uA  vB  wC
            C
                         u v  w  1
Linear algebra solution
       Ax    Bx   Cx     u  Q x 
       A     By   Cy     v  = Q 
        y                  y
        1   1    1    w   1 
Fractional areas
   B
         C
 The barycenctric coordinates are
 the fractional areas
      B
                    w  Area(ABQ)
                                            A
                    w
                    Q
                u
                        v
u  Area(BCQ)
                            v  Area(CAQ)
                C
Barycentric coordinates
   B
                      w  0           A
                          Q
             u  1
                              v  0
         C
Barycentric coordinates
               area(QBC)
           u 
               area(ABC)
               area(QCA)
           v 
               area(ABC)
               area(QAB)
           w 
               area(ABC)
Barycentric coordinates are
fractional
         line segment : fractional length
         triangles : fractional area
         tetrahedrons : fractional volume
Computing Area
                C
      A                             B
                         1
          signed area=     cross B-A,C-A 
                         2
Q outside the triangle
    B
         C
Q outside the triangle
    B
                                 A
          w0
Q
             v0
                         v+w>1
         C
Q outside the triangle
    B
Q
    u0
          C
Voronoi versus Barycentric
   Voronoi regions != barycentric coordinate
    regions
   The barycentric regions are still useful
Barycentric regions of a triangle
              C
Interior
               C       u  0, v  0, w  0
Negative u
                  A
    Q
        u0
              C
Negative v
             C   v0
Negative w
             Q
       B
                     w0
                           A
                 C
The uv regions are not exclusive
            Q
Finding the Voronoi region
   Use the barycentric coordinates to identify
    the Voronoi region
   Coordinates for the 3 line segments and the
    triangle
   Regions must be considered in the correct
    order
First: vertex regions
  uAB  0
  vBC  0   B
                           A   v AB  0
                               uCA  0
                uBC  0
                v CA  0
Second: edge regions
                 uAB  0
         B       v AB  0
                 ?
                            A
             C
Second: edge regions solved
                 uAB  0
                 v AB  0
         B
                 w ABC  0
                             A
             C
Third: interior region
                              A
                  uABC > 0
                  v ABC > 0
                  w ABC > 0
              C
Closest point
   Find the Voronoi region for point Q
   Use the barycentric coordinates to compute
    the closest point Q
Example 1
 Q
     B
            C
Example 1
 Q
     B
            C
                    uAB <= 0
Example 1
 Q
     B
            C
                uAB <= 0 and vBC <= 0
Example 1
 Q
     P=B
            C       Conclusion:
                    P=B
Example 2
                Q
     B
            C
Example 2
                Q
     B
            C
                    Q is not in any vertex region
Example 2
                Q
     B
                        uAB > 0
Example 2
                Q
     B
            C
                    uAB > 0 and vAB > 0
Example 2
                Q
     B
            C
                    uAB > 0 and vAB > 0
                    and wABC <= 0
Example 2
                    Q
     B
                        A
                P
            C
                        Conclusion:
                        P = uAB*A + vAB*B
Implementation
      input: A, B, C, Q
      compute uAB, vAB, uBC, vBC, uCA, vCA
      compute uABC, vABC, wABC
      // Test vertex regions
      …
      // Test edge regions
      …
      // Else interior region
      …
Testing the vertex regions
       // Region A
       if (vAB <= 0 && uCA <= 0)
         P = A
         return
       // Similar tests for Region B and C
Testing the edge regions
       // Region AB
       if (uAB > 0 && vAB > 0 && wABC <= 0)
         P = uAB * A + vAB * B
         return
       // Similar for Regions BC and CA
Testing the interior region
      // Region ABC
      assert(uABC > 0 && vABC > 0 && wABC > 0)
      P = Q
      return
Section 3
Point to Convex Polygon
Convex polygon
          B
  D
              E
Polygon structure
       struct Polygon
       {
         Vec2* points;
         int count;
       };
Convex polygon: closest point
            B
                         Query point Q
        C
                    A
Q
    D
                E
Convex polygon: closest point
                B
                            Closest point Q
        C
                        A
Q
            P
    D
                    E
How do we compute P?
What do we know?
                   Closest point to point
                   Closest point to line segment
                   Closest point to triangle
Simplex
 0-simplex   1-simplex   2-simplex
Idea: inscribe a simplex
            B
                    A
Q
    D
                E
Idea: closest point on simplex
          B
    P=C
                  A
Q
    D
              E
Idea: evolve the simplex
            B
                    A
Q
    D
                E
Simplex vertex
       struct SimplexVertex
       {
           Vec2 point;
           int index;
           float u;
       };
Simplex
      struct Simplex
      {
        SimplexVertex vertexA;
        SimplexVertex vertexB;
        SimplexVertex vertexC;
        int count;
      };
We are onto a winner!
The GJK distance algorithm
   Computes the closest point on a convex
    polygon
   Invented by Gilbert, Johnson, and Keerthi
The GJK distance algorithm
   Inscribed simplexes
   Simplex evolution
Starting simplex
            B
                        Start with arbitrary
        C               vertex. Pick E.
                    A   This is our starting
Q                       simplex.
    D
                E
Closest point on simplex
            B
                           P is the closest
        C                  point.
                      A
Q
    D
                P=E
Search vector
                B
                              Draw a vector
        C                     from P to Q.
                          A   Call this vector d.
Q
            d
    D
                    P=E
Find the support point
                B
                              Find the vertex on
        C                     polygon furthest in
                              direction d.
                          A
Q                             This is the support
            d                 point.
    D
                    P=E
Support point code
    int Support(const Polygon& poly, const Vec2& d)
    {
      int index = 0;
      float maxValue = Dot(d, poly.points[index]);
      for (int i = 1; i < poly.count; ++i)
      {
        float value = Dot(d, poly.points[i]);
        if (value > maxValue)
        {
          index = i;
          maxValue = value;
        }
      }
      return index;
    }
Support point found
                B
                            C is the support
        C                   point.
                        A
Q
            d
    D
                    E
Evolve the simplex
            B
                        Create a line
        C               segment CE.
                    A   We now have a
Q                       1-simplex.
    D
                E
Repeat the process
                B
                            Find closest point
        C                   P on CE.
            P
                        A
Q
    D
                    E
New search direction
                B
                            Build d as a line
        C                   pointing from P to
                            Q.
    d       P
                        A
Q
    D
                    E
New support point
            B
                        D is the support
        C               point.
    d               A
Q
    D
                E
Evolve the simplex
            B
                        Create triangle CDE.
        C
                        This is a
                    A   2-simplex.
Q
    D
                E
Closest point
                B
                            Compute closest
        C                   point on CDE to Q.
                        A
Q
            P
    D
                    E
E is worthless
                B
                            Closest point is on
        C                   CD.
                        A   E does not
Q                           contribute.
            P
    D
                    E
Reduced simplex
                B
                            We dropped E,
        C                   so we now have
                            a 1-simplex.
                        A
Q
            P
    D
                    E
Termination
                B
                            Compute support
        C                   point in direction d.
            d           A   We find either C or
Q                           D. Since this is a
            P               repeat, we are
                            done.
    D
                    E
GJK algorithm
Input: polygon and point Q
pick arbitrary initial simplex S
loop
  compute closest point P on S
  cull non-contributing vertices from S
  build vector d pointing from P to Q
  compute support point in direction d
  add support point to S
end
DEMO!!!
Numerical Issues
   Search direction
   Termination
   Poorly formed polygons
A bad direction
                B
                            d can be built from
        C                   PQ.
            P
    d                   A   Due to round-off:
Q
                            dot(Q-P, C-E) != 0
    D
                    E
A real example in single precision
         Line Segment
         A = [0.021119118, 79.584320]
         B = [0.020964622, -31.515678]
         Query Point
         Q = [0.0 0.0]
         Barycentric Coordinates
         (u, v) = (0.28366947, 0.71633047)
         Search Direction
         d = Q – P = [-0.021008447, 0.0]
         dot(d, B – A) = 3.2457051e-006
Small errors matter
              correct
              support
                                   wrong
                                   support
    exact direction     numerical direction
An accurate search direction
            B
                            Directly compute a
        C                   vector
                            perpendicular to
                        A   CE.
                d
Q
                            d = cross(C-E,z)
    D                       Where z is normal
                    E       to the plane.
The dot product is exactly zero
   edge direction:     e  x y
   search direction:   d   -y x 
   dot product:        e d  -xy+yx=0
Fixing the sign
                B
                            Flip the sign of d so
        C                   that:
                        A   dot(d, Q – C) > 0
Q
                            Perk: no divides
            d
    D
                    E
Termination conditions
                }
Case 1: repeated support point
                B
            d           A
Q
            P
    D
                    E
Case 2: containment
                B
                            We find a 2-simplex
      C                     and all vertices
                            contribute.
                        A
          P=Q
  D
                    E
 Case 3a: vertex overlap
            B
                           We will compute
        C                  d=Q-P as zero.
                    A      So we terminate if
                           d=0.
P=Q=D           E
Case 3b: edge overlap
                  B
                              d will have an
          C                   arbitrary sign.
              d
                          A
P=Q
      D
                      E
Case 3b: d points left
            B
                         If we search left, we
        C                get a duplicate
d                        support point.
                    A    In this case we
                         terminate.
P=Q
    D
                E
Case 3b: d points right
              B
                          If we search right,
      C                   we get a new
                          support point (A).
          d
                      A
Q=P
  D
                  E
Case 3b: d points right
              B
                          But then we get
      C                   back the same P,
                          and then the same
          d           A   d.
P=Q                       Soon, we detect a
                          repeated support
  D                       point or detect
                  E       containment.
Case 4: interior edge
              B
                          d will have an
      C                   arbitrary sign.
          d
                      A
              P=Q
  D
                  E
Case 4: interior edge
              B
                          Similar to Case 3b
      C
          d
                      A
              P=Q
  D
                  E
Termination in 3D
   May require new/different conditions
   Check for distance progression
Non-convex polygon
                      Vertex B is non-
      C               convex
          B       A
  D
              E
Non-convex polygon
                      B is never a
      C               support point
          B       A
  D
              E
Collinear vertices
                         B, C, and D are
      B                  collinear
                     A
  C
  D
              E
Collinear vertices
                         2-simplex BCD
      B
                     A
Q C
  D
              E
Collinear vertices
                         area(BCD) = 0
      B
                     A
Q C
  D
              E
Section 4
Convex Polygon to Convex Polygon
Closest point between convex
polygons
          X        Y
What do we know?
             GJK
What do we need to know?
              ???
Idea
   Convert polygon to polygon into point to
    polygon
   Use GJK to solve point to polygon
Minkowski difference
        Y       Y-X
  X                    Z
Minkowski difference definition
      Z = y j - xi : xi  X, y j  Y
Building the Minkowski difference
Input: polygon X and Y
array points
for all xi in X
  for all yj in Y
    points.push_back(yj – xi)
  end
end
polygon Z = ConvexHull(points)
Example point cloud
           Y             Y-X
  X
      Compute Yi – Xj for i = 1 to 4 and j = 1 to 3
Building the convex hull
     Compute the convex hull by shrink wrapping the points.
Building the convex hull
     Compute the convex hull by shrink wrapping the points.
Building the convex hull
     Compute the convex hull by shrink wrapping the points.
Building the convex hull
     Compute the convex hull by shrink wrapping the points.
Building the convex hull
     Compute the convex hull by shrink wrapping the points.
Building the convex hull
     Compute the convex hull by shrink wrapping the points.
Building the convex hull
     Compute the convex hull by shrink wrapping the points.
Building the convex hull
     Compute the convex hull by shrink wrapping the points.
Building the convex hull
     Compute the convex hull by shrink wrapping the points.
The final polygon
                    Z
Property 1: distances are equal
  X      Y                        O          Z
         distance(X,Y) == distance(O, Y-X)
Property 2: support points
  -d            d                                d
  X            Y                     O          Z
       support(Z, d) = support(Y, d) – support(X, -d)
Convex Hull?
Modifying GJK
   Change the support function
   Simplex vertices hold two indices
Closest point on polygons
   Use the barycentric coordinates to compute
    the closest points on X and Y
   See the demo code for details
    DEMO!!!
Download: box2d.org
Further reading
   Collision Detection in Interactive 3D
    Environments, Gino van den Bergen, 2004
   Real-Time Collision Detection, Christer
    Ericson, 2005
   Implementing GJK:
    http://mollyrocket.com/849, Casey Muratori.
Box2D
   An open source 2D physics engine
   http://www.box2d.org
   Written in C++