'--------------------------------------------------------------------------------------------------------- '************************************ Routine SphereSegmentIntersect ************************************* '--------------------------------------------------------------------------------------------------------- 'VBA implementation of the theory provided by Paul Bourke (http://paulbourke.net/geometry/circlesphere/) 'author: ing. Giuseppe Iaria - rev. 12/08/2014 '--------------------------------------------------------------------------------------------------------- 'Routine determines the intersection point of a segment and a sphere 'The segment is defined by its end points P1 and P2, the sphere is defined by its centre C and radius R 'If exsist, actual points P of intersection are returned in a 2-dimensions array P() 'The P() array is implemented assuming a Base 1 condition, so the "Option Base 1" statement is required 'The first dimension of the array P() refers to point number (maybe 1 or 2) 'The second dimension of the array P() refers to point axis coordinates (1 = x, 2 = y, 3 = z) 'Possible points of intersection coordinates are defined as (P(1,1);P(1,2);P(1,3)) (P(2,1);P(2,2);P(2,3)) 'Both conditions of line and segment intersection are valued and returned 'Routine arguments are defined as follows: '- end points of segment: P1 (x1 ; y1), P2 (x2 ; y2) '- sphere: C (xC ; yC), R '- condition of line intersection: LineIntersection '- condition of segment intersection: SegmentIntersection '- array of points of intersection: P() '--------------------------------------------------------------------------------------------------------- Option Base 1 Public Sub SphereSegmentIntersect(ByVal x1 As Double, ByVal y1 As Double, ByVal z1 As Double, _ ByVal x2 As Double, ByVal y2 As Double, ByVal z2 As Double, _ ByVal xC As Double, ByVal yC As Double, zC As Double, _ ByVal R As Double, ByRef P() As Double, _ ByRef LineIntersection As Boolean, ByRef SegmentIntersection As Boolean) Dim a As Double, b As Double, c As Double Dim bbm4ac As Double, u1 As Double, u2 As Double, u As Double Const EPS As Single = 0.0001 a = (x2 - x1) ^ 2 + (y2 - y1) ^ 2 + (z2 - z1) ^ 2 b = 2 * ((x2 - x1) * (x1 - xC) + (y2 - y1) * (y1 - yC) + (z2 - z1) * (z1 - zC)) c = xC ^ 2 + yC ^ 2 + zC ^ 2 + x1 ^ 2 + y1 ^ 2 + z1 ^ 2 - 2 * (xC * x1 + yC * y1 + zC * z1) - R ^ 2 bbm4ac = b ^ 2 - 4 * a * c 'Verify if line-sphere intersection occurs If (bbm4ac < 0 Or Abs(a) < EPS) Then 'NO intersection line-sphere LineIntersection = False SegmentIntersection = False Exit Sub Else 'line-sphere intersection occurs LineIntersection = True u1 = (-b - Sqr(bbm4ac)) / 2 / a u2 = (-b + Sqr(bbm4ac)) / 2 / a End If 'segment-sphere condition n. 1: segment is outside sphere, no point of intersection If ((u1 < 0 And u2 < 0) Or (u1 > 1 And u2 > 1)) Then SegmentIntersection = False Exit Sub End If 'segment-sphere condition n. 2: segment is inside sphere, no point of intersection If ((u1 < 0 And u2 > 1) Or (u2 < 0 And u1 > 1)) Then SegmentIntersection = False Exit Sub End If 'segment-sphere condition n. 3: one point of segment is inside sphere, 1 point of intersection If (((u1 >= 0 And u1 <= 1) And (u2 < 0 Or u2 > 1)) Or ((u2 >= 0 And u2 <= 1) And (u1 < 0 Or u1 > 1))) Then SegmentIntersection = True If (u2 >= 0 And u2 <= 1) Then u = u2 Else u = u1 ReDim P(1, 3) P(1, 1) = x1 + u * (x2 - x1) P(1, 2) = y1 + u * (y2 - y1) P(1, 3) = z1 + u * (z2 - z1) Exit Sub End If 'segment-sphere condition n. 4: segment is tangent sphere, 1 point of intersection If ((Abs(u1 - u2) < EPS) And (u1 >= 0 And u1 <= 1)) Then SegmentIntersection = True ReDim P(1, 3) P(1, 1) = x1 + u1 * (x2 - x1) P(1, 2) = y1 + u1 * (y2 - y1) P(1, 3) = z1 + u1 * (z2 - z1) Exit Sub End If 'segment-sphere condition n. 5: segment is outside sphere, 2 point of intersection If ((u1 >= 0 And u1 <= 1) And (u2 >= 0 And u2 <= 1)) Then SegmentIntersection = True ReDim P(2, 3) P(1, 1) = x1 + u1 * (x2 - x1) P(1, 2) = y1 + u1 * (y2 - y1) P(1, 3) = z1 + u1 * (z2 - z1) P(2, 1) = x1 + u2 * (x2 - x1) P(2, 2) = y1 + u2 * (y2 - y1) P(2, 3) = z1 + u2 * (z2 - z1) Exit Sub End If End Sub