Location and time of impact for a moving point and a moving line segment (Continuous Collision Detection)

TLDR

I did read "...like 5 minutes to evaluate..."

No way too long, this is a real-time solution for many lines and points.

Sorry this is not a complete answer (I did not rationalize and simplify the equation) that will find the point of intercept, that I leave to you.

Also I can see several approaches to the solution as it revolves around a triangle (see image) that when flat is the solution. The approach bellow finds the point in time when the long side of the triangle is equal to the sum of the shorter two.

Solving for u (time)

This can be done as a simple quadratic with the coefficients derived from the 3 starting points, the vector over unit time of each point. Solving for u

The image below give more details.

• The point P is the start pos of point
• The points L1, L2 are the start points of line ends.
• The vector V1 is for the point, over unit time (along green line).
• The vectors V2,V3 are for the line ends over unit time.
• u is the unit time
• A is the point (blue), and B and C are the line end points (red)

There is (may) a point in time u where A is on the line B,C. At this point in time the length of the lines AB (as a) and AC (as c) sum to equal the length of line BC (as b) (orange line).

That means that when b - (a + c) == 0 the point is on the line. In the image the points are squared as this simplifies it a little. b² - (a² + c²) == 0

At the bottom of image is the equation (quadratic) in terms of u, P, L1, L2, V1, V2, V3.

That equation needs to be rearranged such that you get (???)u² + (???)u + (???) = 0

Sorry doing that manually is very tedious and very prone to mistakes. I don`t have the tools at hand to do that nor do I use python so the math lib you are using is unknown to me. However it should be able to help you find how to calculate the coefficients for (???)u² + (???)u + (???) = 0

Update

Ignore most of the above as I made a mistake. b - (a + c) == 0 is not the same as b² - (a² + c²) == 0. The first one is the one needed and that is a problem when dealing with radicals (Note that there could still be a solution using a + bi == sqrt(a^2 + b^2) where i is the imaginary number).

Another solution

So I explored the other options.

The simplest has a slight flaw. It will return the time of intercept. However that must be validated as it will also return the time for intercepts when it intercepts the line, rather than the line segment BC

Thus when a result is found you then test it by dividing the dot product of the found point and line segment with the square of the line segments length. See function isPointOnLine in test snippet.

To solve I use the fact that the cross product of the line BC and the vector from B to A will be 0 when the point is on the line.

Some renaming

Using the image above I renamed the variables so that it is easier for me to do all the fiddly bits.

/*
point P is  {a,b}
point L1 is  {c,d}
point L2 is  {e,f}
vector V1 is {g,h}
vector V2 is {i,j}
vector V3 is {k,l}

Thus for points A,B,C over time u    */
Ax = (a+g*u)
Ay = (b+h*u)
Bx = (c+i*u)
By = (d+j*u)
Cx = (e+k*u)
Cy = (f+l*u)

/* Vectors BA and BC at u */
Vbax = ((a+g*u)-(c+i*u))
Vbay = ((b+h*u)-(d+j*u))
Vbcx = ((e+k*u)-(c+i*u))
Vbcy = ((f+l*u)-(d+j*u))

/*
   thus Vbax * Vbcy - Vbay * Vbcx == 0 at intercept 
*/

This gives the quadratic

0 = ((a+g*u)-(c+i*u)) * ((f+l*u)-(d+j*u)) - ((b+h*u)-(d+j*u)) * ((e+k*u)-(c+i*u))

Rearranging we get

0 = -((i*l)-(h*k)+g*l+i*h+(i+k)*j-(g+i)*j)*u* u -(d*g-c*l-k*b-h*e+l*a+g*f+i*b+c*h+(i+k)*d+(c+e)*j-((f+d)*i)-((a+c)*j))*u +(c+e)*d-((a+c)*d)+a*f-(c*f)-(b*e)+c*b

The coefficients are thus

 A = -((i*l)-(h*k)+g*l+i*h+(i+k)*j-(g+i)*j)
 B = -(d*g-c*l-k*b-h*e+l*a+g*f+i*b+c*h+(i+k)*d+(c+e)*j-((f+d)*i)-((a+c)*j))
 C = (c+e)*d-((a+c)*d)+a*f-(c*f)-(b*e)+c*b

We can solve using the quadratic formula (see image top right).

Note that there could be two solutions. In the example I ignored the second solution. However as the first may not be on the line segment you need to keep the second solution if within the range 0 <= u <= 1 just in case the first fails. You also need to validate that result.

Testing

To avoid errors I had to test the solution

Below is a snippet that generates a random random pair of lines and then generate random lines until an intercept is found.

The functions of interest are

• movingLineVPoint which return the unit time of first intercept if any.
• isPointOnLine to validate the result.

const ctx = canvas.getContext("2d");
canvas.addEventListener("click",test);
const W = 256, H = W, D = (W ** 2 * 2) ** 0.5;
canvas.width  = W; canvas.height = H;
const rand = (m, M) => Math.random() * (M - m) + m;
const Tests = 300;
var line1, line2, path, count = 0; 
setTimeout(test, 0);

// creating P point L line
const P = (x,y) => ({x,y,get arr() {return [this.x, this.y]}}); 
const L = (l1, l2) => ({l1,l2,vec: P(l2.x - l1.x, l2.y - l1.y), get arr() {return [this.l1, this.l2]}}); 
const randLine = () => L(P(rand(0, W), rand(0, H)), P(rand(0, W), rand(0, H)));
const isPointOnLine = (p, l) =>  {
    const x = p.x - l.l1.x;
    const y = p.y - l.l1.y;
    const u = (l.vec.x * x + l.vec.y * y) / (l.vec.x * l.vec.x + l.vec.y * l.vec.y);
    return u >= 0 && u <= 1;
}
// See answer illustration for names
// arguments in order Px,Py,L1x,l1y,l2x,l2y,V1x,V1y,V2x,V2y,V3x,V3y
function movingLineVPoint(a,b, c,d, e,f, g,h, i,j, k,l) {
    var A = -(i*l)-(h*k)+g*l+i*h+(i+k)*j-(g+i)*j;
    var B = -d*g-c*l-k*b-h*e+l*a+g*f+i*b+c*h+(i+k)*d+(c+e)*j-((f+d)*i)-((a+c)*j)
    var C = +(c+e)*d-((a+c)*d)+a*f-(c*f)-(b*e)+c*b

    // Find roots if any. Could be up to 2
    // Using the smallest root >= 0 and <= 1
    var u, D, u1, u2;
    // if A is tiny we can ignore
    if (Math.abs(A) < 1e-6) { 
        if (B !== 0) {
            u = -C / B;
            if (u < 0 || u > 1) { return }  // !!!!  no solution  !!!!
        } else { return }                   // !!!!  no solution  !!!!
    } else {
        B /= A;
        D = B * B - 4 * (C / A);
        if (D > 0) {
            D **= 0.5;
            u1 = 0.5 * (-B + D);
            u2 = 0.5 * (-B - D);
            if ((u1 < 0 || u1 > 1) && (u2 < 0 || u2 > 1))  { return }  // !!!!  no solution  !!!!
            if (u1 < 0 || u1 > 1) { u = u2 }        // is first out of range
            else if (u2 < 0 || u2 > 1) { u = u1 }   // is second out of range
            else if (u1 < u2) { u = u1 }            // first is smallest
            else { u = u2 }
        } else if (D === 0) {
            u = 0.5 * -B;
            if (u < 0 || u > 1)  { return }  // !!!!  no solution  !!!!            
        } else { return }                    // !!!!  no solution  !!!! 
    }    
    return u;
}

function test() {
   if (count> 0) { return }
   line1 = randLine();
   line2 = randLine();
   count = Tests
   subTest();
}
function subTest() {
   path = randLine()
   ctx.clearRect(0,0,W,H);
   drawLines();
   const u = movingLineVPoint(
       path.l1.x, path.l1.y,
       line1.l1.x, line1.l1.y,
       line2.l1.x, line2.l1.y,
       path.vec.x, path.vec.y,
       line1.vec.x, line1.vec.y,
       line2.vec.x, line2.vec.y
   );
   
   if (u !== undefined) { // intercept found maybe
      pointAt = P(path.l1.x + path.vec.x * u, path.l1.y + path.vec.y * u);
      lineAt = L(
          P(line1.l1.x + line1.vec.x * u, line1.l1.y + line1.vec.y * u),
          P(line2.l1.x + line2.vec.x * u, line2.l1.y + line2.vec.y * u)
      );
      const isOn = isPointOnLine(pointAt, lineAt);
      if (isOn) {
          drawResult(pointAt, lineAt);
          count = 0;
          info.textContent = "Found at: u= " + u.toFixed(4) + ". Click for another";
          return;
      }
   }
   setTimeout((--count < 0 ? test : subTest), 18);
}   








function drawLine(line, col = "#000", lw = 1) {
    ctx.lineWidth = lw;
    ctx.strokeStyle = col;
    ctx.beginPath();
    ctx.lineTo(...line.l1.arr);
    ctx.lineTo(...line.l2.arr);
    ctx.stroke();
}
function markPoint(p, size = 3, col = "#000", lw = 1) {
    ctx.lineWidth = lw;
    ctx.strokeStyle = col;
    ctx.beginPath();
    ctx.arc(...p.arr, size, 0, Math.PI * 2);
    ctx.stroke();
}
function drawLines() {
   drawLine(line1);
   drawLine(line2);
   markPoint(line1.l1);
   markPoint(line2.l1);
   drawLine(path, "#0B0", 1);
   markPoint(path.l1, 2, "#0B0", 2);
}
function drawResult(pointAt, lineAt) {
   ctx.clearRect(0,0,W,H);
   drawLines();
   markPoint(lineAt.l1, 2, "red", 1.5);
   markPoint(lineAt.l2, 2, "red", 1.5);
   markPoint(pointAt, 2, "blue", 3);
   drawLine(lineAt, "#BA0", 2);
}

div {position: absolute; top: 10px; left: 12px}
canvas {border: 2px solid black}

<canvas id="canvas" width="1024" height="1024"></canvas>
    <div><span id="info">Click to start</span></div>

There are two parts of @Blindman67's solution I don't understand:

• Solving for b^2 - (a^2 + c^2) = 0 instead of sqrt(b^2)-(sqrt(a^2)+sqrt(b^2)) = 0
• The returned timestamp being clamped in the range [0,1]

Maybe I'm missing something obvious, but in any case, I designed a solution that addresses these concerns:

• All quadratic terms are solved for, not just one
• The returned time stamp has no limits
• sqrt(b^2)-(sqrt(a^2)+sqrt(b^2)) = 0 is solved for, instead of b^2 - (a^2 + c^2) = 0

Feel free to recommend ways this could be optimized:

# pnt, crt_1, and crt_2 are points, each with x,y and dx,dy attributes
# returns a list of timestamps for which pnt is on the segment
# whose endpoints are crt_1 and crt_2
def colinear_points_collision(pnt, crt_1, crt_2):
    a, b, c, d = pnt.x, pnt.y, pnt.dx, pnt.dy
    e, f, g, h = crt_1.x, crt_1.y, crt_1.dx, crt_1.dy
    i, j, k, l = crt_2.x, crt_2.y, crt_2.dx, crt_2.dy

    m = a - e
    n = c - g
    o = b - f
    p = d - h
    q = a - i
    r = c - k
    s = b - j
    u = d - l
    v = e - i
    w = g - k
    x = f - j
    y = h - l

    # Left-hand expansion
    r1 = n * n + p * p
    r2 = 2 * o * p + 2 * m * n
    r3 = m * m + o * o
    r4 = r * r + u * u
    r5 = 2 * q * r + 2 * s * u
    r6 = q * q + s * s

    coef_a = 4 * r1 * r4  # t^4 coefficient
    coef_b = 4 * (r1 * r5 + r2 * r4)  # t^3 coefficient
    coef_c = 4 * (r1 * r6 + r2 * r5 + r3 * r4)  # t^2 coefficient
    coef_d = 4 * (r2 * r6 + r3 * r5)  # t coefficient
    coef_e = 4 * r3 * r6  # constant

    # Right-hand expansion
    q1 = (w * w + y * y - n * n - p * p - r * r - u * u)
    q2 = 2 * (v * w + x * y - m * n - o * p - q * r - s * u)
    q3 = v * v + x * x - m * m - o * o - q * q - s * s

    coef1 = q1 * q1  # t^4 coefficient
    coef2 = 2 * q1 * q2  # t^3 coefficient
    coef3 = 2 * q1 * q3 + q2 * q2  # t^2 coefficient
    coef4 = 2 * q2 * q3  # t coefficient
    coef5 = q3 * q3  # constant

    # Moves all the coefficients onto one side of the equation to get
    # at^4 + bt^3 + ct^2 + dt + e
    # solve for possible values of t
    p = np.array([coef1 - coef_a, coef2 - coef_b, coef3 - coef_c, coef4 - coef_d, coef5 - coef_e])

    def fun(x):
        return p[0] * x**4 + p[1] * x**3 + p[2] * x**2 + p[3] * x + p[4]

    # could use np.root, but I found this to be more numerically stable
    sol = optimize.root(fun, [0, 0], tol=0.002)
    r = sol.x

    uniques = np.unique(np.round(np.real(r[np.isreal(r)]), 4))
    final = []

    for r in uniques[uniques > 0]:
        if point_between(e + g * r, f + h * r, i + k * r, j + l * r, a + c * r, b + d * r):
            final.append(r)

    return np.array(final)

# Returns true if the point (px,py) is between the endpoints
# of the line segment whose endpoints lay at (ax,ay) and (bx,by)
def point_between(ax, ay, bx, by, px, py):
    # colinear already checked above, this checks between the other two.
    return (min(ax, bx) <= px <= max(ax, bx) or abs(ax - bx) < 0.001) and (min(ay, by) <= py <= max(ay, by) or abs(ay - by) < 0.001)

An example (L1 and L2 are endpoints of line):
P = (0,0) with velocity (0, +1)
L1 = (-1,2) with velocity (0, -1)
L2 = (1,2) with velocity (0, -1)

The returned result would be t=1, because after 1 time step, P will be one unit higher, and both endpoints of the line will each be one unit lower, therefore, the point intersects the segment at t=1.