Intersections
The objective of this section is to cover geometric constructions using points and vectors. In particular, this section focuses on computing intersections between geometric entities.
Line Plane
Provided with a line ( a, b ), defined by two points a and b; and plane ( o, n ), defined by its origin point o and its unit normal vector n, determine the point of their intersection.
""" Inputs
"""
a, b = Point3d( ax, ay, az ), Point3d( bx, by, bz )
o, n = Point3d( ox, oy, oz ), Vector3d( nx, ny, nz )
Solution
The point of intersection p must be on the line, therefore it can be expressed as p = a * ( 1 - t ) + b * t, or by using vector u = b - a as p = a + u * t. Meanwhile, the point is also on the plane, therefore its projected distance to the plane must be zero. This is expressed as n · ( p - o ) = 0, see the implicit form of the plane.
n · ( p - o ) = 0
or with substituting
n · ( a + u * t - o ) = 0
or with expanding
n · ( a - o ) + n · u * t = 0
or with solving
t = n · ( o - a ) / ( n · u )
Computing the point of intersection requires constructing a couple of vectors and obtaining their dot products with the plane's normal vector.
""" Line Plane Intersection
"""
def LinePlaneFast( a, b, o, n ):
u = b - a
v = o - a
return ( n * v ) / ( n * u )
t = LinePlaneFast( a, b, o, n )
print( 't: {:>.5f}'.format( t ) )
p = a * ( 1 - t ) + b * t
Note however that the denominator dot product n · u can become zero if the plane's normal and the line's direction are orthogonal. This is because the line is parallel to the plane and thus they do not intersect.
""" Line Plane Intersection
"""
def LinePlaneSafe( a, b, o, n, epsilon = 1e-5 ):
u = b - a
v = o - a
d = n * u
if( abs( d ) <= epsilon ):
return False, 0.0
return True, ( n * v ) / d
success, t = LinePlaneSafe( a, b, o, n )
if( success ):
print( 't: {:>.5f}'.format( t ) )
p = a * ( 1 - t ) + b * t
else:
print( 'Parallel' )
Rhino's Intersection class provides the LinePlane( ) method which computes the point of intersection, if it exists, or returns a failure status otherwise.
success, t = Intersection.LinePlane(
Line( a, b ), Plane( o, n ) )
if( success ):
print( 't: {:>.5f}'.format( t ) )
p = a * ( 1 - t ) + b * t
else:
print( 'Parallel' )
Plane Plane Plane
Provided with three planes ( a, u ), ( b, v ) and ( c, n ), defined by their origin points a, b, and c; and their unit normal vectors u, v and n, respectively, determine the point of their intersection.
""" Inputs
"""
a, u = Point3d( ax, ay, az ), Vector3d( ux, uy, uz )
b, v = Point3d( bx, by, bz ), Vector3d( vx, vy, vz )
c, n = Point3d( cx, cy, cz ), Vector3d( nx, ny, nz )
Solution
This problem can be approached using standard linear algebra methods, for instance Cramer's rule, since there are three plane equations and three unknown coordinates, of the point of intersection. The conversion between the implicit origin-normal and explicit equation form of a plane is described in the plane section.
u.X * x + u.Y * y + u.Z * z = a · u
v.X * x + n.Y * y + n.Z * z = b · v
n.X * x + v.Y * y + v.Z * z = c · n
or with matrix notation
[ u.X u.Y u.Z ] [ x ] [ a · u ]
[ v.X v.Y v.Z ] × [ y ] = [ b · v ]
[ n.X n.Y n.Z ] [ z ] [ c · n ]
Note that the intersection will fail if the determinant of the basis formed by the plane normals is zero. Geometrically, this the result of two or more planes being parallel with one another.
''' Determinant of 3x3 Matrix
'''
def Det3x3( xx, yx, zx,
xy, yy, zy,
xz, yz, zz ):
return ( xx * ( yy * zz - zy * yz ) -
yx * ( xy * zz - zy * xz ) +
zx * ( xy * yz - yy * xz ) )
''' Triple Plane Intersection
'''
def TriplePlane( a, u, b, v, c, n ):
''' Plane Distances
'''
d = Vector3d(
u * Vector3d( a ),
v * Vector3d( b ),
n * Vector3d( c ) )
''' Basis Determinant
'''
D = Det3x3( u.X, u.Y, u.Z,
v.X, v.Y, v.Z,
n.X, n.Y, n.Z )
if( D == 0.0 ): return False, D
''' Coordinate Determinants
'''
X = Det3x3( d.X, u.Y, u.Z,
d.Y, v.Y, v.Z,
d.Z, n.Y, n.Z )
Y = Det3x3( u.X, d.X, u.Z,
v.X, d.Y, v.Z,
n.X, d.Z, n.Z )
Z = Det3x3( u.X, u.Y, d.X,
v.X, v.Y, d.Y,
n.X, n.Y, d.Z )
return True, Point3d( X / D, Y / D, Z / D )
success, point = TriplePlane( a, u, b, v, c, n )
Rhino's Intersection class provides the PlanePlanePlane( ) method for computing the triple-plane intersection, using exactly the same semantics.
''' Triple Plane Intersection
'''
success, point = Intersection.PlanePlanePlane(
Plane( a, u ), Plane( b, v ), Plane( c, n ) )
Plane Plane
Provided with two planes ( a, u ) and ( b, v ), defined by their origin points a and b; and their unit normal vectors u and v, respectively, determine the point of their intersection.
""" Inputs
"""
a, u = Point3d( ax, ay, az ), Vector3d( ux, uy, uz )
b, v = Point3d( bx, by, bz ), Vector3d( vx, vy, vz )
Solution
The intersection between two planes, if they are not parallel, produces a line expressing their common points. The direction of the line must be orthogonal to both plane normals, so that the line is in both planes. Therefore, it can be trivially computed by the normals' cross product n = u × v.
The simplest way to find an origin point o for the line of intersection is by reusing the triple-plane method. The third plane's origin c can be defined by either of the two planes' origins a or b, or their mid-point ( a + b ) / 2. The direction of the line can be used as the third plane's normal because if the two other planes are not parallel, then the cross product of their normals ensures an orthogonal direction which will not cause problems with the determinant.
''' Plane Plane Intersection
'''
def PlanePlane( a, u, b, v, epsilon = 1e-10 ):
n = Vector3d.CrossProduct( u, v )
if( n * n <= epsilon ):
return False, a, n
n.Unitize( )
c = ( a + b ) / 2
success, o = TriplePlane( a, u, b, v, c, n )
return success, o, n
success, o, n = PlanePlane( a, u, b, v )
Alternatively it is possible to derive the intersection from first principles. The system of equations used in the triple-plane intersection captures the relationships the line of intersection should obey, namely its projected distance to both planes must be zero, its direction is the cross product, and its points are free along the line. Replacing the term c · n with a parameter t is equivalent of allowing the point to slide along the line. If we multiply both sides with the inverse of the normal vectors' matrix and expand the result into vector products produces the equation of line.
[ u.X u.Y u.Z ] [ x ] [ a · u ]
[ v.X v.Y v.Z ] × [ y ] = [ b · v ]
[ n.X n.Y n.Z ] [ z ] [ t ]
or with multiplying the inverse of A
[ x ] [ m.X o.X w.X ] [ a · u ]
[ y ] = [ m.Y o.Y w.Y ] × [ b · v ]
[ z ] [ m.Z o.Z w.Z ] [ t ]
or with expanding the vector products
p = m * ( a · u ) +
o * ( b · v ) +
w * ( t )
Note that we do not care about the term w because we already know that the direction of the line is n. All we need to resolve is the origin point m * a · u + o * b · v. The Cramer's rule can be used here to invert the matrix, but instead of the repeating the same 3x3 determinant, let's instead replace them with regular vector arithmetic.
d = ( u.X * ( v.Y * n.Z - v.Z * n.Y ) -
u.Y * ( v.X * n.Z - v.Z * n.X ) +
u.Z * ( v.X * n.Y - v.Y * n.X ) )
or with vector algebra
d = u * ( v × n )
Moreover, the inverse of the matrix, in terms of its columns, can be reduced to three cross product vectors scaled by the inverse of the determinant. Note that the column vector w is indeed the line's direction n = u × v.
mxx = ( v.Y * n.Z - v.Z * n.Y ) / d
mxy = ( v.Z * n.X - v.X * n.Z ) / d
mxz = ( v.X * n.Y - v.Y * n.X ) / d
oyx = ( u.Z * n.Y - u.Y * n.Z ) / d
oyy = ( u.X * n.Z - u.Z * n.X ) / d
oyz = ( u.Y * n.X - u.X * n.Y ) / d
wzx = ( u.Y * v.Z - u.Z * v.Y ) / d
wzy = ( u.Z * v.X - u.X * v.Z ) / d
wzz = ( u.X * v.Y - u.Y * v.X ) / d
or with vector algebra
m = ( v × n ) / d
o = ( n × u ) / d
w = ( u × v ) / d
The python code below implements the plane-plane intersection logic using only simple vector algebra operations.
''' Plane Plane Intersection
'''
def PlanePlane( a, u, b, v, epsilon = 1e-10 ):
''' Line's Direction
'''
n = Vector3d.CrossProduct( u, v )
if( n * n <= epsilon ):
return False, a, n
n.Unitize( )
''' Compute Determinant
'''
vn = Vector3d.CrossProduct( v, n )
nu = Vector3d.CrossProduct( n, u )
d = u * vn
''' Origin of Line
'''
au = Vector3d( a ) * u
bv = Vector3d( b ) * v
o = Point3d( au * vn / d +
bv * nu / d )
return True, o, n
Rhino's Intersection class provides the PlanePlane( ) method for computing the intersection between two planes, with the same semantics.
''' Plane Plane Intersection
'''
success, line = Intersection.PlanePlane(
Plane( a, u ),
Plane( b, v ) )
Line Line (Fast)
Provided with two lines ( a, b ) and ( c, d ) defined by two points each, determine the point of their intersection. Assume that the intersection exists.
""" Inputs
"""
a, b = Point3d( ax, ay, az ), Point3d( bx, by, bz )
c, d = Point3d( cx, cy, cz ), Point3d( dx, dy, dz )
Solution
Constructive Geometry
Intersection between two lines can be approached in a constructive geometry sense without need for algebraic transformations. This is performed by using the concept of mapping between global to local coordinates: (1) The line a->b and the line a->c are used to establish an orthonormal basis ( a, u, v, n ). (2) The line c->d can be expressed in 2D local coordinates C->D in the ( a, u, v ) plane. (3) Finally, the point of intersection between the u-axis and C->D line is computed using the similar triangles ratios.
''' Constructive Intersection
'''
def LineLineFast( a, b, c, d ):
''' Orthonormal ( a, b, c ) Basis
'''
u = b - a
n = Vector3d.CrossProduct( u, c - a )
v = Vector3d.CrossProduct( n, u )
u.Unitize( )
v.Unitize( )
n.Unitize( )
''' World to Local Mapping
'''
C = Point3d( u * ( c - a ), v * ( c - a ), 0 )
D = Point3d( u * ( d - a ), v * ( d - a ), 0 )
''' Similar Triangles
'''
t = ( D.Y * C.X - D.X * C.Y ) / ( D.Y - C.Y )
return t / a.DistanceTo( b )
t = LineLineFast( a, b, c, d )
print( 't: {:>8.5f}'.format( t ) )
p = a * ( 1 - t ) + b * t
Symbolic Perspective
Symbolically there are quite a few different ways to compute the point of intersection. We consider two points at parameters s and t on the lines a->b and c->d, respectively. The two points are p = a * ( 1 - s ) + b * s and q = c * ( 1 - t ) + d * t must be equal p = q for the lines to intersect. Substituting the point equations yields a * ( 1 - s ) + b * s = c * ( 1 - t ) + d * t or simply a + ( b - a ) * s = c + ( d - c ) * t.
Since the points a, b, c and d are known, we can define the vectors u = b - a, v = d - c and w = c - a to simplify the expression to u * s - v * t = w. Therefore, we have one equation with two unknowns. The points p and q can be situated anywhere on the lines so we need to constraint them based on some geometric principle. To solve this we either (1) need at least one more equation to solve a 2x2 system, (2) apply a distance constraint such as |p - q| = 0, or (3) perform some symbolic trick to get rid of one unknown, such as either u * s or v * t.
Geometric Constraints
One way of thinking about this is by requiring that the line p->q is perpendicular to both a->b and c->d. If this is the case then p and q are the closest points between the lines. Therefore, if the lines do intersect, then p and q they must be coincident. This can be expressed using the dot product as u · ( q - p ) = 0 and at the same time v · ( q - p ) = 0. But q - p = w + v * t - u * s, thus expanding the above yields two equations with two unknown. Note that the coefficients of s and t are easily computed dot products.
u · ( q - p ) = 0
v · ( q - p ) = 0
or with substituting
u · ( w + v * t - u * s ) = 0
v · ( w + v * t - u * s ) = 0
or with expanding
uw + uv * t - uu * s = 0
vw + vv * t - vu * s = 0
or with solving
s = ( uv * vw - vv * uw ) / ( uv * uv - vv * uu )
t = ( uu * vw - uv * uw ) / ( uv * uv - vv * uu )
This approach is based on Daniel Sunday's derivation. The python code below reflects this process minimally without considering special cases and potential geometric and numerical issues involved.
def LineLineFast( a, b, c, d ):
u = b - a
v = d - c
w = c - a
uu = u * u
uv = u * v
vv = v * v
uw = u * w
vw = v * w
s = ( uv * vw - vv * uw ) / ( uv * uv - vv * uu )
t = ( uu * vw - uv * uw ) / ( uv * uv - vv * uu )
return s, t
Calculous / Minimization
Requiring that the distance between the points p and q is zero, yields an unpleasant relationship |q - p| = 0 involving square roots. However, we can approach this as a minimization problem instead. This requires setting the partial derivatives θs( |q - p| ** 2 ) = 0 and θt( |q - p| ** 2 ) = 0. The square of the distance between points is equivalent to the dot product of the vectors ( q - p ) · ( q - p ). Performing the requisite algebra yields a quadratic expression seen below.
( q - p ) · ( q - p )
or with substituting
( w + v * t - u * s ) · ( w + v * t - u * s )
or with expanding
w · w + w · v * t - w · u * s +
w · v * t + v · v * t * t - u · v * s * t -
w · u * s - u · v * s * t + u · u * s * s
or with combining
u · u * s * s +
v · v * t * t -
u · v * s * t * 2 +
u · w * s * 2 +
v · w * t * 2 -
w · w
The partial derivatives of this expression with respect to s and t, drop a degree and create a system of two linear equations with two unknowns. Note that this is exactly the same linear system seen earlier.
θs = 0
θt = 0
or with differentiating
u · u * s * 2 - u · v * t * 2 - u · w * 2 = 0
v · v * t * 2 - u · v * s * 2 + v · w * 2 = 0
or with compacting
uu * s - uv * t - uw = 0
vv * t - uv * s + vw = 0
or with solving
s = ( uv * vw - vv * uw ) / ( uv * uv - vv * uu )
t = ( uu * vw - uv * uw ) / ( uv * uv - vv * uu )
Algebraic Perspective
Alternatively, we can get rid of one unknown in the equation u * s - v * t = w by multiplying, in the cross product sense, with u or v. This is because the cross product with a vector and itself is the zero vector. Therefore, we arrive at two vector-scalar product expressions as seen below.
u * s - v * t = w
or with crossing
v × ( u * s - v * t ) = v × w
u × ( u * s - v * t ) = u × w
or with expanding
v × u * s - v × v * t = v × w
u × u * s - u × v * t = u × w
or with eliminating
v × u * s = v × w
- u × v * t = u × w
or with reordering
u × v * s = w × v
u × v * t = w × u
Finally, to drop from an equation between vectors and scalars we multiply, in the dot product sense, with u × v. Note that the vector being dotted must not be parallel to the terms involved, otherwise we are introducing zeros. The normal vector n = u × v is safe because it is orthogonal to u, v and w.
u × v * s = w × v
u × v * t = w × u
or with n = u × v
n * s = w × v
n * t = w × u
or with dotting
n · n * s = ( w × v ) · n
n · n * t = ( w × u ) · n
or with solving
s = ( ( w × v ) · n ) / ( n · n )
t = ( ( w × u ) · n ) / ( n · n )
This approach is based on Steven Jenke's derivation. The python code below reflects this process minimally without considering special cases and potential geometric and numerical issues involved.
''' Line and Line Intersection
'''
def LineLineFast( a, b, c, d ):
u = b - a
v = d - c
w = c - a
n = Vector3d.CrossProduct( u, v )
s = Vector3d.CrossProduct( w, v ) * n / ( n * n )
t = Vector3d.CrossProduct( w, u ) * n / ( n * n )
return s, t
s, t = LineLineFast( a, b, c, d )
print( 's: {:>8.5f}'.format( s ) )
print( 't: {:>8.5f}'.format( t ) )
Rhino Intersection
Rhino's Intersection class supports computing the intersection between two lines using the LineLine( ) static method. The method returns a success / failure status code, as well as the two parameters of lines where the lines intersect, if they are on the same plane.
""" Outputs
"""
from Rhino.Geometry.Intersect import Intersection
success, s, t = Intersection.LineLine(
Line( a, b ),
Line( c, d ) )
if( success ):
print( 'Success' )
p = a * ( 1 - s ) + b * s
q = c * ( 1 - t ) + d * t
else:
print( 'Failure' )
print( 's: {:>8.5f}'.format( s ) )
print( 't: {:>8.5f}'.format( t ) )
Line Line (Safe)
Provided with two lines ( a, b ) and ( c, d ) defined by two points each, determine the point of their intersection.
""" Inputs
"""
a, b = Point3d( ax, ay, az ), Point3d( bx, by, bz )
c, d = Point3d( cx, cy, cz ), Point3d( dx, dy, dz )
Solution
This problem has the same solution as seen above but asks for considering the special cases that can cause either of the methods to fail. Geometrically, there are several scenarios seen below. The objective is to account for them.
- If either line is degenerate, in the sense that the distance between its two points is zero, then this will generate zero-length vectors, which will eventually raise division by zero exceptions.
- If all points are colinear, then it is even more troublesome because it can either mean that there are infinite solutions represented by their overlap, or none if we assume the lines as bounded segments.
- If the lines are parallel, then the cross product of their directions will become the zero vector causing some denominator to also become zero, resulting again in a division by zero exceptions.
- If all points of both lines are not on the same plane, then there is no intersection point, because they are disjoint.
''' Line and Line Intersection
'''
def LineLineSafe( a, b, c, d, epsilon = 1e-5 ):
''' Degenerate Cases
'''
u = b - a
if( u * u <= epsilon ):
return 'Degenerate', 1.0, 0.0
v = d - c
if( v * v <= epsilon ):
return 'Degenerate', 0.0, 1.0
''' Parallel & Colinear Cases
'''
w = c - a
n = Vector3d.CrossProduct( u, v )
if( n * n <= epsilon ):
''' Compute Distance Between Lines
'''
t = ( u * w ) / ( u * u )
d = w - u * t
if( d * d <= epsilon ):
return 'Colinear', 0.0, 0.0
else:
return 'Parallel', d, 0.0
''' Intersection & Disjoint Cases
'''
s, t = LineLineFast( a, b, c, d )
if( abs( w * n ) <= epsilon ):
return 'Success', s, t
else:
return 'Disjoint', s, t
r, s, t = LineLineSafe( a, b, c, d )
print( 'Result: '.format( r ) )
print( 's: {:>8.5f}'.format( s ) )
print( 't: {:>8.5f}'.format( t ) )