← Back to index

Collision Detection Using GJK and a Modified EPA

Epicrelyt
#23914
Hello, I am a highschool student and have been working on a 2D game with a friend. I have done most of the programming from scratch using only 3 STB libraries. Recently I went to rewrite the physics engine for the game as it had a few bugs and didn't scale well—I had initially written it from following Handmade Hero. Initially I was planning for the new system to use discrete collision detection with GJK and EPA. After implenting the two, I found that GJK worked fine, however I encountered a problem with EPA.

The problem was that EPA doesn't take into account the change in position of the colliding objects and only returns the closest edge. I quickly found out that this does not work if you have multiple tiles constituting a platform or floor for the player to run on. EPA would end up returning the side edge of the tile, whereas the correct edge would be the top one, this unfotunately resulted in the player falling through the floor. I read about a few solutions to this, however I thought that it should be possible to modify EPA to take into account the dP of the object.

Rather simply this leads to a modified EPA where the polygon is not expanded based on the closest edge, but instead based on whichever edge is intersecting with the dP of the object. Like EPA the algorithm would iterate on the simplex returned by GJK, meaning that the polygon must contain the origin so if the dP is treated as a ray it must intersect with an edge. A few cases were possible:

1: The dP doesn't intersect with the edge at all, this edge isn't considered
2: The dP is colinear with the edge, this edge is kept unless the dP intersects another edge
3: The dP intersects with the edge, this edge is kept and the loop searching for the edge is broken
3.5. If the dP were a ray it would intersect with the edge, this edge is kept unless the dP is colinear with another edge or intersects another

The rest of the algorithm basically does the same thing as EPA. It uses the intersecting edge to expand the polygon in the direction given by the normal of the edge. If it is able to expand the polygon then the algorithm would loop again, otherwise you then have the intersecting edge. Using this edge you can calculate the intersection point. The time of impact could be calculated by either 1-IntersectionPoint.X/dP.X or 1-IntersectionPoint.Y/dP.Y. The time of impact would be clamped to be within the range of 0 to 1. The normal would then be the inverse of the normal used to expand the polygon. If the intersection point is further along the direction of dP than dP itself (case 3.5—and the time of impact would <0), then the object was colliding before the start of the frame. This allows the algorithm to do a bit of penetration correction.

The GJK support functions would be modified to return the object swept along its dP, in order to prevent tunnelling. In addition to account for multiple moving objects you would use the difference of the dPs of the two objects, then reduce it to one moving object with the new dP and one static one.

I have implented this algorithm and support functions for squares, wedges, and circles. Thankfully I haven't encountered any major flaws in the way this algorithm works. One flaw I have observed is that the moving object can clip the corner of a wedge, which isn't too desirable but can be solved by adding a step height. It seems similar to Casey's idea of searching for a valid position instead of searching for a valid time. I however think that it likely has been documented somewhere else in a paper or something, however I am not entirely sure how to find it if it actually has or not. So far I have been unable to think of any obvious major drawbacks of using this algorithm as it currently works just fine for what I have done so far. I am also unsure how well this would scale into 3 dimensions, as it will be more complicated to calculate the intersecting triangle face of the polytope. Hopefully the explanation I have given isn't too bad.
Edited by Epicrelyt on Reason: Clarification
Simon Anciaux
#23915
Welcome to the forums.

I'm not sure if you're asking a question or just sharing you experience with GJK and EPA, but I would appreciate if you could share some code. I implemented GJK and EPA in the past and I remember having issues at the time, part of it was floating point precision issues. I never used it much as I decided to focus on other projects, but I'd like to see how you implemented it.
Epicrelyt
#23917
Here is the code I have currently in the game. However I have not yet implented collision between multiple moving objects, so currently this code only supports a single moving object but multiple static objects. I have also not included my GJK implentation as it just an ordinary GJK implementation.

Support functions: 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
struct physics_object {
    v2 P, dP, ddP;
    f32 Mass;
    collision_boundary *Boundaries;
    u8 BoundaryCount;
};

internal inline v2
GJKSupport(collision_boundary *Boundary, 
           physics_object *Object, v2 Delta,
           v2 Direction){
    v2 Result = {};
    
    switch(Boundary->Type){
        case BoundaryType_Rect: {
            v2 Min = Boundary->Bounds.Min;
            v2 Max = Boundary->Bounds.Max;
            v2 Points[4] = {
                Min,
                V2(Min.X, Max.Y),
                Max,
                V2(Max.X, Min.Y),
            };
            for(u32 I=0; I < 4; I++){
                if(Dot(Points[I], Direction) > Dot(Result, Direction)){
                    Result = Points[I];
                }
            }
        }break;
        case BoundaryType_Circle: {
            f32 DirectionLength = Length(Direction);
            v2 UnitDirection = Direction/DirectionLength;
            Result = UnitDirection * Boundary->Radius;
        }break;
        case BoundaryType_Wedge: {
            for(u32 I=0; I < 3; I++){
                v2 Point = Boundary->Points[I];
                if(Dot(Point, Direction) > Dot(Result, Direction)){
                    Result = Point;
                }
            }
        }break;
        default: INVALID_CODE_PATH;
    }
    
    Result += Object->P;
    if(Dot(Delta, Direction) > 0.0f){
        Result += Delta;
    }
    
    return(Result);
}

internal inline v2 
CalculateSupport(physics_object *ObjectA, physics_object *ObjectB, v2 Delta, v2 Direction){
    v2 Result = GJKSupport(ObjectA->Boundaries, ObjectA, Delta, Direction) - GJKSupport(ObjectB->Boundaries, ObjectB, V20, -Direction);
    return(Result);
}


The modified EPA: 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
struct epa_result {
    f32 ActualTimeOfImpact;
    f32 WorkingTimeOfImpact;
    v2 Normal;
};

internal epa_result
DoVelocityEPA(physics_object *ObjectA, physics_object *ObjectB, v2 Delta, v2 Simplex[3]){
    TIMED_FUNCTION();
    const f32 Epsilon = 0.0001f;
    
    dynamic_array<v2> Polytope; 
    DynamicArrayInitialize(&Polytope, 20, &TransientStorageArena);
    DynamicArrayPushBack(&Polytope, Simplex[0]);
    DynamicArrayPushBack(&Polytope, Simplex[1]);
    DynamicArrayPushBack(&Polytope, Simplex[2]);
    
    v2 DeltaNormal = Normalize(Clockwise90(Delta));
    
    epa_result Result = {};
    while(true){
        b8 FoundEdge = false;
        b8 IsColinear = false;
        u32 EdgeIndex = 0;
        f32 EdgeDistance = 0;
        v2 InverseNormal = {};
        v2 IntersectionPoint = {};
        
        for(u32 I=0; I < Polytope.Count; I++){
            v2 A = Polytope[I];
            v2 B = Polytope[(I+1)%Polytope.Count];
            
            f32 AAlongDeltaNormal = Dot(DeltaNormal, A);
            f32 BAlongDeltaNormal = Dot(DeltaNormal, B);
            
            if(((-Epsilon <= AAlongDeltaNormal) && (AAlongDeltaNormal <= Epsilon)) && // Collinear
               ((-Epsilon <= BAlongDeltaNormal) && (BAlongDeltaNormal <= Epsilon))){
                FoundEdge = true;
                IsColinear = true;
                EdgeIndex = I;
                InverseNormal = Normalize(TripleProduct(B-A, -A));
                EdgeDistance = Dot(InverseNormal, -A);
                IntersectionPoint = V20;
            }else if(((AAlongDeltaNormal >= 0) && (BAlongDeltaNormal <= 0)) || 
                     ((AAlongDeltaNormal <= 0) && (BAlongDeltaNormal >= 0))){ // The delta line  intersects
                AAlongDeltaNormal = AbsoluteValue(AAlongDeltaNormal);
                BAlongDeltaNormal = AbsoluteValue(BAlongDeltaNormal);
                f32 ABPercent = AAlongDeltaNormal/(AAlongDeltaNormal + BAlongDeltaNormal);
                v2 Point = A + ABPercent*(B-A);
                
                v2 DeltaDirection = Normalize(Delta);
                f32 PointAlongDeltaDirection = Dot(DeltaDirection, Point);
                f32 DeltaLength = Dot(DeltaDirection, Delta);
                
                if((-Epsilon <= PointAlongDeltaDirection) && 
                   (PointAlongDeltaDirection <= DeltaLength + Epsilon)){ // The delta intersects
                    FoundEdge = true;
                    EdgeIndex = I;
                    InverseNormal = Normalize(TripleProduct(B-A, -A));
                    EdgeDistance = Dot(InverseNormal, -A);
                    IntersectionPoint = Point;
                    
                    break;
                }else if(-Epsilon <= PointAlongDeltaDirection){ // The  delta ray intersects
                    if(!IsColinear){
                        FoundEdge = true;
                        EdgeIndex = I;
                        InverseNormal = Normalize(TripleProduct(B-A, -A));
                        EdgeDistance = Dot(InverseNormal, -A);
                        IntersectionPoint = Point;
                    }
                }
            }
        }
        Assert(FoundEdge); // Catch any potential bugs
        
        f32 Distance;
        v2 NewPoint = V20;
        if((InverseNormal.X == 0.0f) && (InverseNormal.Y == 0.0f)){
            Distance = EdgeDistance;
        }else{
            v2 Normal = -InverseNormal;
            NewPoint = CalculateSupport(ObjectA, ObjectB, Delta, Normal);
            Distance = Dot(NewPoint, Normal);
        }

        if((-Epsilon <= (Distance-EdgeDistance)) && ((Distance-EdgeDistance) <= Epsilon)){
            f32 Percent = IntersectionPoint.X/Delta.X;
            if((Delta.X == 0.0f) && (Delta.Y == 0.0f)){
                Percent = 0.0f;
            }else if(Delta.X == 0.0f){
                Percent = IntersectionPoint.Y/Delta.Y;
            }
            
            
            Result.WorkingTimeOfImpact = 1.0f - Percent;
            if(Result.WorkingTimeOfImpact > 1.0f){
                Result.WorkingTimeOfImpact = 1.0f;
            }else if(Result.WorkingTimeOfImpact < 0.0f){
                Result.WorkingTimeOfImpact = 0.0f;
            }
            Result.ActualTimeOfImpact = 1.0f - Percent; // This is for penetration correction
            
            Result.Normal = InverseNormal;
            
            break;
        }else{
            DynamicArrayInsertNewArrayItem(&Polytope, EdgeIndex+1, NewPoint);
        }
    }
    
    return(Result);
}


Using the result from the modified EPA 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
    epa_result Colliding = {};
    Colliding.ActualTimeOfImpact = F32_POSITIVE_INFINITY;
    f32 CollidingObjectDistanceAlongDelta = F32_POSITIVE_INFINITY;
                
    for(u32 I = 0; I < CollidingObjects.Count; I++){ // Colliding objects is just an array of objects that GJK has deemed to be colliding
        physics_object *ObjectB = CollidingObjects[I];
        v2 *Simplex = &CollidingSimplices[I*3];
        epa_result EPAResult = DoVelocityEPA(ObjectA, ObjectB, Delta, Simplex);
        if(EPAResult.ActualTimeOfImpact < Colliding.ActualTimeOfImpact){
            Colliding = EPAResult;
            CollidingObjectDistanceAlongDelta = Dot(Delta, ObjectB->P-ObjectA->P);
        }else if(EPAResult.ActualTimeOfImpact == Colliding.ActualTimeOfImpact){
            f32 ObjectBDistanceAlongDelta = Dot(Delta, ObjectB->P-ObjectA->P);
            if(ObjectBDistanceAlongDelta < CollidingObjectDistanceAlongDelta){
                Colliding = EPAResult;
                CollidingObjectDistanceAlongDelta = ObjectBDistanceAlongDelta;
            }
        }
    }
    
    // I still need to properly implement collision response, I have just finished collision detection
    f32 COR = 1.0f;
                
    FrameTimeRemaining -= FrameTimeRemaining*Colliding.WorkingTimeOfImpact;
    ObjectA->P += Colliding.ActualTimeOfImpact*Delta;
    ObjectA->dP = ObjectA->dP - COR*Dot(Colliding.Normal, ObjectA->dP)*Colliding.Normal;
    Delta = Delta - COR*Dot(Colliding.Normal, Delta)*Colliding.Normal;


Beware, this is the current implementation—minus the debug code—I have fixed all the bugs I have encountered, but more likely remain. It is also possible that there are certain cases which I have not considered but might discover later. I have also just implemented VelocityEPA and finished debugging it, so it could probably be cleaned up significantly. To clarify as well, I do have actually a few questions:

Firstly, from the little bit of googling I did trying to find what this modification of EPA, I haven't found anything. I figure it likely might exist somewhere, but I don't know how I would begin finding if it truly does. It might not matter if it already exists, but even if it doesn't I am still curious. How would one go about finding if a certain algorithm has been described somewhere else?

Secondly, I believe my reasoning is sound in that the delta - IntersectionPoint, is actually the distance the object can move along delta. As it is basically what EPA returns, just along a certain vector. Is that correct or have I not realized something? How would one know if an algorithm works correctly? or realize major drawbacks? I don't entirely understand this process.

Thirdly, as might be evident by the code, I don't fully understand when floating point precision becomes a problem. The epsilons I have might not actually be necessary if I did something different. So when should epsilons be used? How small should they be?
Simon Anciaux
#23918
Thanks for the code, I'll try to have a look when I have some time.

Epicrelyt
How would one go about finding if a certain algorithm has been described somewhere else?


I have no idea. The best way I think would be to find someone knowledgeable about the subject and ask them. Unfortunately twitter becomes the "easiest" way to contact people and that's not great to express complex ideas. You can ask in the handmade.network forums, irc and discord, but we might not know the answer either.

Sometimes reading wikipedia pages and sources at the bottom can help you find things you didn't knew.

For your second question, I believe you're correct (from what I remember about EPA, which isn't much). If your game run correctly, than the algorithm is working, at least most of the time. If at times there is something unexpected happening, than there might be an issue somewhere. Identifying the cause of the issue will help you know if the problem is with the algorithm, the way you use it or something else.

You can build some test cases (making a square collide with a square, a circle...) in a controlled environment where you control the stepping process (e.g.: you have to press a button to advance the simulation; you can set an arbitrary time step).

Adding some visualization is almost always worth the time (for example being able to visualize current point tested in GJK, visualize each step, not just the final result of the physics step...).

A good thing to have but sometimes hard to do is to record the state of the application (or part of it) to be able to replay and debug it.

For your third question: I had similar question s.

I use epsilon when I would expect a computation to be exact, but after testing, the result is not exact and it matters (it's noticeable) in the result of what I'm doing. That doesn't help much but I don't think there is a single answer. Knowing more about how floating point number works helps though.

How small should epsilon be depends on each case (e.g. if 1.0f is a meter, and the smallest thing you represent is 1cm, than having an epsilon of 0.001f is enough, at least for position). And you can have different epsilon for different things in a single algorithm (position, angle...).

Epicrelyt
#23920
Thank you for your advice. I had thought about making some sort of visual debugging thing for collisions. I didn't think it would be worth the effort, though it would have been highly useful. I expect there will likely be more collision bugs that come out of hiding. So I'll follow your advice and write a visual debugging thing for collisions, for even if it isn't much use it will probably also be a good thing to learn for the future.
Edited by Epicrelyt on
Reply to Thread
← Back to index