Intersection of two lines in 3D by QR / SVD

2022-01

Finding the intersection of two lines in 3D in a less ad-hoc manner using surface constraints and QR / SVD factorization.

Doing TA for vision class this quarter. The formulas in lecture 3 slides 14 are neat: cross-product of 2 points gives the common line; cross-product of 2 lines gives the intersecting point. But all these happen in 2D. What to do in 3D? If you look it up online, some suggestions like this, this, this look ad-hoc. Not as memorable as cross-products.

The essence of cross-product in the 2D case is to find the orthogonal constraints that define a line. In

\R^2

, a single linear constraint

\va \cdot \vx = b

(or in homogeneous coordinates

[\va, -b] \cdot [\vx, 1] = 0

) defines a line. We can generalize this notion of orthogonality to 3D and beyond. In

\R^3

, one linear constraint cuts a plane, and two linear constraints cut a line. Orthogonality is easy to implement using QR or SVD factorization.

Therefore the generalized way to handle intersection in 3D is that QR/SVD of 2 points give the common line; QR/SVD of 2 lines give the common point. Here is a simple test.

import numpy as np
from numpy.linalg import svd, qr


def line_constraint_repr_by_qr(p_1, p_2):
    A = np.block([
        [p_1, 1],
        [p_2, 1]
    ])
    Q, R = qr(A.T, mode="complete")
    constraints = Q[:, 2:].T
    return constraints


def line_constraint_repr_by_svd(p_1, p_2):
    A = np.block([
        [p_1, 1],
        [p_2, 1]
    ])
    U, s, V_t = svd(A)
    constraints = V_t[2:, :]
    return constraints


def main():
    """
    In 3D:
        qr/svd of 2 points gives the line
        qr/svd of 2 lines gives the point
    """
    p_intersect = np.array([1, 2, 3])
    p_a = np.array([9, 4, 5])
    p_b = np.array([-7, 3, 16])

    # random scaling to avoid directly using the shared point
    line1_p1, line1_p2 = p_a, p_a + (p_intersect - p_a) * 0.56
    line2_p1, line2_p2 = p_b, p_b + (p_intersect - p_b) * 0.88

    line1_cons = line_constraint_repr_by_qr(line1_p1, line1_p2)
    line2_cons = line_constraint_repr_by_qr(line2_p1, line2_p2)

    cons = np.block([
        [line1_cons],
        [line2_cons]
    ])

    _, s, V_t = svd(cons)  # or QR.

    # last singular value being 0 means the lines intersect exactly
    assert np.allclose(s[-1], 0)
    p = V_t[-1]
    p = p / p[-1]
    assert np.allclose(p[:3], p_intersect)


if __name__ == "__main__":
    main()

Last updated on 2024-12-31. Design inspired by distill.