Exact (integer arithmetic only) capsule intersection test

Posted on June 6, 2015 by Aaron Rotenberg
Post tagged: Haskell, Mathematics

The capsule (also called a “pill”) is one of a few common types of shape used for collision detection in games and other applications. A 3-D capsule is a cylinder capped by hemispheres. The 2-D equivalent, a rectangle capped by semicircles, is technically called a stadium, but they are commonly just called capsules as well.

Capsules used for collision detection in Super Smash Bros. Melee. From here.

Capsules are nice because they can form both spherical and elongated shapes in any direction. The animation above shows how Super Smash Bros. Melee uses spherical hitboxes that are “stretched” across frames into capsules to prevent fast-moving attacks from going through opponents without hitting them. (Marvel vs. Capcom 3 uses the same trick.) What really makes capsules useful is that they have a very simple mathematical description: a capsule is the set of all points less than a certain radius from a line segment. This means you can check whether two capsules intersect each other by just finding the shortest distance between the two line segments and checking whether it is less than the sum of the radii.

Calculating the distance between two line segments is a well-known problem. This StackOverflow answer gives the code to do that with floating-point arithmetic. Sometimes, though, approximating the correct answer with floating-point isn’t good enough. What if we want an exact intersection test for capsules using only integer arithmetic?

I’ll be giving code examples in Haskell. The code will be for 2-D capsules, but the 3-D case is not too different. Let’s start with some basic definitions using the vector-space package. Since we’re using integer arithmetic, all of our vectors and radii should have integer values only.

{-# LANGUAGE TypeFamilies #-}
import Data.VectorSpace

-- Use arbitrary-size integers to avoid overflow in later calculations.
-- If you are using very small values only, this may not be necessary.
type GeomInt = Integer

data Vec = Vec { vecX, vecY :: !GeomInt }
    deriving (Show)
instance AdditiveGroup Vec where
    zeroV = Vec 0 0
    Vec x1 y1 ^+^ Vec x2 y2 = Vec (x1 + x2) (y1 + y2)
    negateV (Vec x y) = Vec (-x) (-y)
instance VectorSpace Vec where
    type Scalar Vec = GeomInt
    s *^ Vec x y = Vec (s * x) (s * y)
instance InnerSpace Vec where
    Vec x1 y1 <.> Vec x2 y2 = x1 * x2 + y1 * y2

-- Represents a *closed* 2-D line segment. Zero-length segments are allowed.
data Segment = Segment { segmentEnd1, segmentEnd2 :: !Vec }
    deriving (Show)

-- Represents an *open* 2-D stadium (disk-capped rectangle). It is required
-- that capsuleRadius > 0.
data Capsule = Capsule {
    capsuleSegment :: !Segment,
    capsuleRadius :: !GeomInt
  } deriving (Show)

The first step of the line segment distance computation is to test whether the line segments intersect. The test shown in the StackOverflow answer doesn’t work for our purposes because it uses floating-point division and because it treats parallel segments as never intersecting. Instead, we can use the exact test from this page.

segmentsIntersect :: Segment -> Segment -> Bool
segmentsIntersect (Segment p1 q1) (Segment p2 q2) =
    (o1 /= o2 && o3 /= o4)
    || (o1 == Collinear && onSegment p1 p2 q1)
    || (o2 == Collinear && onSegment p1 q2 q1)
    || (o3 == Collinear && onSegment p2 p1 q2)
    || (o4 == Collinear && onSegment p2 q1 q2)
  where o1 = orientation p1 q1 p2
        o2 = orientation p1 q1 q2
        o3 = orientation p1 q1 p1
        o4 = orientation p1 q1 q1

data Orientation = Collinear | Clockwise | Counterclockwise
    deriving (Show, Eq)

orientation :: Vec -> Vec -> Vec -> Orientation
orientation (Vec px py) (Vec qx qy) (Vec rx ry) =
    case compare val 0 of
        LT -> Counterclockwise
        EQ -> Collinear
        GT -> Clockwise
  where val = (qy - py) * (rx - qx) - (qx - px) * (ry - qy)

-- onSegment p q r checks if q lies on the segment pr assuming that
-- p, q, and r are collinear.
onSegment :: Vec -> Vec -> Vec -> Bool
onSegment (Vec px py) (Vec qx qy) (Vec rx ry) =
    qx <= max px rx && qx >= min px rx && qy <= max py ry && qy >= min py ry

Now here’s the tricky bit. If the segments do not intersect, we can’t simply find the distance between them to check against the radii, because the shortest distance may not be an integer. The standard trick of doing all comparisons on squared distance values to avoid square root operations doesn’t completely solve the problem either, because the closest point on one segment to the other may not even have integer coordinates.

If we were using imprecise floating-point arithmetic, the test would look like this:

capsulesIntersect :: Capsule -> Capsule -> Bool
capsulesIntersect (Capsule s1@(Segment p1 q1) r1)
        (Capsule s2@(Segment p2 q2) r2) =
    segmentsIntersect s1 s2 ||
        check p1 s2 || check q1 s2 || check p2 s1 || check q2 s1
  where thresholdSq = (r1 + r2)^2
        
        check :: Vec -> Segment -> Bool
        check p (Segment e1 e2)
            | segLenSq == 0 || t <= 0 =
                magnitudeSq (p ^-^ e1) < thresholdSq
            | t >= 1 = magnitudeSq (p ^-^ e2) < thresholdSq
            | otherwise = magnitudeSq (p ^-^ near) < thresholdSq
          where d = e2 ^-^ e1
                segLenSq = magnitudeSq d
                near = e1 ^+^ t *^ d
                t = ((p ^-^ e1) <.> d) / segLenSq

Since we’re using integer arithmetic, though, the (/) operator is banned. The trick to pulling this off with integers only is to scale both sides of the third inequality by segLenSq. This will cancel the denominator so that we don’t have to do any division. We can use primed variables to denote “multiplied by a factor of segLenSq”. The exact capsule intersection is then:

capsulesIntersect :: Capsule -> Capsule -> Bool
capsulesIntersect (Capsule s1@(Segment p1 q1) r1)
        (Capsule s2@(Segment p2 q2) r2) =
    segmentsIntersect s1 s2 ||
        check p1 s2 || check q1 s2 || check p2 s1 || check q2 s1
  where thresholdSq = (r1 + r2)^2
        
        check :: Vec -> Segment -> Bool
        check p (Segment e1 e2)
            | t' <= 0 = magnitudeSq (p ^-^ e1) < thresholdSq
            | t' >= segLenSq = magnitudeSq (p ^-^ e2) < thresholdSq
            | otherwise = magnitudeSq (p' ^-^ near') < thresholdSq''
          where d = e2 ^-^ e1
                segLenSq = magnitudeSq d
                
                thresholdSq'' = segLenSq^2 * thresholdSq
                p' = segLenSq *^ p
                near' = segLenSq *^ e1 ^+^ t' *^ d
                t' = (p ^-^ e1) <.> d

Notice also that we don’t have to check segLenSq == 0 anymore, because the t’ <= 0 case implicitly covers that.

Yay, math!