I have segments of lines and ellipses (NOT planes and ellipsoids) transformed into 3d space and need to calculate whether two given segments intersect and where.
I'm using Java but haven't yet found a library which solves my problem, nor came across some algorithms that I could use for my own implementation.
For line-line intersection test there are several ways to solve. The classic way is using linear algebra (i.e., solving a linear matrix system) but from software development point of view I like more the Geometric-Algebra way, in the form of Plucker Coordinates, which only requires to implement vector algebra operations (i.e., cross-product and dot-product) which are simpler to code than matrix operations for solving linear systems.
I will show both for the sake of comparison then you decide:
Linear Algebra Way
Given line segment P limited by points P1 and P2 and line segment Q limited by points Q1 and Q2.
The parametric form of the lines is given by:
P(t) = P1 + t (P2 - P1)
Q(t) = Q1 + t (Q2 - Q1)
Where t is a real number in the interval [0 1].
If two lines intersect then the following equation holds:
P(t0) = Q(t1)
Provided that the two unknown numbers t0 and t1 exist. Expanding the above equation we get:
t0 (P2 - P1) - t1 (Q2 - Q1) = Q1 - P1
We can solve for t0 and t1 by expressing the above equation in matrix algebra:
A x = B
Where A is a 3x2 matrix with coordinates of vector (P2 - P1) in the first column and coordinates of the vector (Q2 - Q1) in the second column; x is a 2x1 column vector of unknowns t0 and t1 and B is a 3x1column vector with coordinates of vector (Q1 - P1).
Classically the system can be solved calculating the pseudo-inverse of matrix A, denoted by A^+:
A^+ = (A^T A)^-1 A^T
See: https://en.m.wikipedia.org/wiki/Generalized_inverse
Fortunately any matrix package in Java should be able to compute the above calculations very easily and perhaps very efficiently too.
If multiplying A with its pseudo-inverse A^+ is equal to the identity matrix I, i.e., (A A^+) == I, then there IS a unique answer (intersection) and you can get it computing the following product:
x = A^+ B
Of course if you cannot compute the pseudo-inverse in the first place, e.g., because (A^T A) is singular (i.e. determinant is zero), then no intersection exists.
Since we are dealing with segments, the intersection point is at point P(x0) or Q(x1) iff x0 and x1 are both in the interval [0 1].
OTHER METHOD OF SOLUTION
To avoid dealing with matrix algebra we can try to solve the system using vector algebra and substitution method:
t0 (P2 - P1) - t1 (Q2 - Q1) = Q1 - P1
t0 = a + t1 b
t1 = C • (Q1 - (1 - a) P1 - a P2) / |C|^2
Where:
a = (P2 - P1) • (Q1 - P1) / |P2 - P1|^2
b = (P2 - P1) • (Q2 - Q1) / |P2 - P1|^2
C = b (P2 - P1) - (Q2 - Q1)
I cannot provide a geometric intuition of the above results yet.
Plucker Coordinates way
Given line segment P limited by points P1 and P2 and line segment Q limited by points Q1 and Q2.
The Plucker Coordinates of line P is given by a pair of 3d vectors (Pd, Pm):
Pd = P2 - P1
Pm = P1 x P2
Where Pm is the cross-product of P1 and P2. The Plucker Coordinates (Qd, Qm) of line Q can be calculated in exactly the same way.
The lines P and Q only can intersect if they are coplanar. Thr lines P and Q are coplnanar iif:
Pd • Qm + Qd • Pm = 0
Where (•) is the dot-product. Since machines have finite precision a robust test should check not for zero but for a small number. If Pd x Qd = 0 then lines are parallel (here 0 is the zero vector). Again a robust test should be for instamce that the squared length of (Pd x Qd) is small.
If the lines are not parallel then they are coplanar and their intersection (called "meet" in Plucker's jargon) will be the point:
x = ((Pm • N) Qd - (Qm • N) Pd - (Pm • Qd) N) / (Pd x Qd) • N
Where N is any coordinate axis vector (i.e., (1,0,0) or (0,1,0) or (0,0,1)), such that (Pd x Qd) • N is non-zero.
If the neither P nor Q pass through the origin, then their Plucker coordinate Pm and Qm respectively will be non-zero and the following sinpler formula will work
x = Pm x Qm / Pd • Qm
For an introduction to Plucker Coordinates see:
https://en.m.wikipedia.org/wiki/Plücker_coordinates
http://www.realtimerendering.com/resources/RTNews/html/rtnv11n1.html#art3
For a general intersection formula see "Corollary 6" of:
http://web.cs.iastate.edu/~cs577/handouts/plucker-coordinates.pdf
Transforming Ellipses to Circles Forth and Back
We can always transform an ellipse into a circle. An ellipse has two "radius", called semi-axes, which you can visualize in your mind as two orthogonal vectors, one big called major semi-axes and one small called minor semi-axes. You can apply a non-uniform scaling transformation to both semi-axes vectors for making them of equal size, so you get a circle.
We define an ellipse "E" by its center O, which is a 3d point and its two semi-axes A1 and A2, which are also 3d vectors. A normal vector N to the ellipse's plane can be computed by the cross product of its semi-axes N = A1 x A2 and then normalizing it to have unit length.
For now suppose there is a linear function M that you can use to transform (scale) your ellipse E into a circle C, with radius equal to the minor semi-axes, by applying it to your ellipse's semi-axes A1 and A2 and to the ellipse's center O.
Notice that the eliipse's center O and ellipse's normal vector N are not changed by M. So M(N) = N and M(O) = O. That means that the circle is in the same plane and has the same position C than the ellipse. The linear function M has a corresponding inverse function M^-1 so we can transform back the vectors of the circle to get the original ellipse E.
Of course we can transform the endpoints of lines P and Q also using function M for sending them to the "circle space" and we can send them back to "ellipse space" using M^-1.
Using M we can compute the intersection of the lines P and Q with the ellipse E in the circle space. So now we can focus on the line-circle intersection.
Line-Plane Intersection
Given a plane with normal vector N and distance D such that:
N • x + D = 0
For every point x in the plane. Then the intersection with line P with Plucker Coordinates (Pd, Pm) is given by:
x = (N x Pm - D Pd) / N • Pd
This works only if the line P is not in the plane i.e.,:
(N • P1 + D) != 0 and (N • P2 + D) != 0
And for our ellipse E we have:
N = (A1 x A2)/|A1 x A2|
D = -N • O
Line-Circle and Point-Circle Intersection
Now having x, the point-in-circle check is easy:
|O - x| <= |A2|
Equality holds only when x is in circle boundary.
If line P is in circle's plane then the following sinple check will give you the answer:
https://stackoverflow.com/a/1079478/9147444
How to Compute the Linear Function M
A simple way to compute M is the following. Use the Ellipse's normal vector N and semi-axes A1 and A2 to create a 3x3 matrix U. Such that U has the vectors A1, A2 and N as columns.
Then scale the major semi-axes A1 to have the same length to the minor semi-axes A2. Then create the matrix V auch that V has the new vector A1 and A2 and N as columns.
Then we define the linear matrix system:
R U = V
Where R is a 3x3 (non-uniform-)scaling-rotation matrix.
We can solve for R by multiplying both sides of the equation by the inverse of U which is denoted by U^-1
R U U^-1 = V U^-1
Since U U^-1 is the identity matrix we get:
R = V U^+
Using R we define the affine transformation
M(x) = R (x - O) + O
We can use M to transform points to circle space, such as O, P1, P2, Q1 and Q2. But if we need to transform vectors such as A1, A2, N, Pd and Qd. We need to use the simpler M:
M(x) = R x
Since A1, A2 and N are eigenvectors of R then R is not singular and has an inverse. We define the inverse M as:
M^-1(x) = R^-1 (x - O) + O
And for vectors:
M^-1(x) = R^-1 x
Update: Ellipse-Ellipse intersection
Two intersecting non-coplanar 3d-ellipses have their intersection points on the line formed by the intersection between their planes. So you first find the line formed by the intersection of the planes containig the ellipses (if planes do not intersect i.e., they are parallel, then neither the ellipses intersect).
The line of intersection belong to both planes, so it is perpendicular to both normals. The direction vector V is the cross product of the plane normals:
V = N1 × N2
To fully determine the line we also need to find a point on the line. We can do that solving the linear equations of the planes. Given the 2x3 matrix N = [N1^T N2^T] with the normals N1 and N2 as rows, and the 2x1 column vector b = [N1 • C1, N2 • C2], where C1 and C2 are the centers of the ellipses, we can form the linear matrix system:
N X = b
Where X is some point satifying both plane equations. The system is underdetermined since there are infinite number of points in the line satifying the system. We can still find a particular solution closer to the origin by using the pseudo-inverse of matrix N, denoted by N^+.
X = N^+ b
The line equation is
L(t) = X + t V
For some scalar t.
Then you can apply the method described earlier to test the line-ellipse intersection i.e., scaling the ellipse to a circle and intersect with the coplanar line. If both ellipses intersect the line in the same points then they intersect.
Now, the case in which the ellipses actually lie on the same plane is more complex. You can ceck the approach taken by Dr Eberly in his excellent book "Geometric Tools" (also available online):
https://www.geometrictools.com/Documentation/IntersectionOfEllipses.pdf
And also you can check the C++ source code which is open source:
https://www.geometrictools.com/GTEngine/Include/Mathematics/GteIntrEllipse2Ellipse2.h
