Skew Lines Midpoint from 4 points?

Question

I'm trying to understand the basics underlying a piece of source code I was given to use. It works, this is proven. I'm just trying to wrap my head around the why of it well enough that I could do it myself, or possibly extend/expand upon it.

The code in question finds the midpoint of the shortest line between two skew lines in 3D space. This paper is the closest I've come to finding something that matches, but I'm still missing some conceptual steps (and my linear algebra skills are decades out of use)

In this application, P1, P2, P3, and P4 are 3D (X,Y,Z) points in space. The lines we're concerned with are P1-P2 and P3-P4.

The language this system runs on doesn't include a Determinant function, hence why the original programmer wrote their own. VectMagn is a system function that simply returns the Norm of a 3D point value (ie, SQRT(X^2 + Y^2 + Z^2)). Pow(i,j) is just what it looks like, returning i^j.

Where I'm most getting stuck is the large formula for t. Based on my research so far, I would expect that both lines would need to be converted into unit-vector lines, then processed as per the first PDF. But the t formula appears to be doing all of this in one jump, and I'm missing the intermediate steps. It's obviously creating two 2x2 matrices from various matrix math on Ps 1-4, then dividing the Determinant of one matrix by the other.

If I'm understanding this correctly, t is r1 and r2 from the PDF, depending on which order the Points were passed to iv3DSkewLinePoint. But I haven't yet found any papers or formulae that explain why/how this particular algorithm works. So far, everything I've found starts with unit-vector lines and moves on from there.

! Returns the derterminant of a matrix
LOCAL FUNC num Det(num a,num b,num c,num d)
  RETURN (a*d-c*b);
ENDFUNC
       
! Returns a point on a line (P1-P2) closest to a point on a skewed line (P3-P4)
FUNC pos iv3DSkewLinePoint(pos P1,pos P2,pos P3,pos P4)
  VAR num x;
  VAR num y;
  VAR num z;
  VAR num t;

  t:=Det(DotProd(P3-P1,P2-P1),DotProd(P4-P3,P2-P1),DotProd(P3-P1,P4-P3),Pow(VectMagn(P4-P3),2))/Det(Pow(VectMagn(P2-P1),2),DotProd(P4-P3,P2-P1),DotProd(P2-P1,P4-P3),Pow(VectMagn(P4-P3),2));
  x:=P1.x+(P2.x-P1.x)*t;
  y:=P1.y+(P2.y-P1.y)*t;
  z:=P1.z+(P2.z-P1.z)*t;

  RETURN [x,y,z];
ENDFUNC

! Returns the closest point to two skewed lines in space
FUNC pos iv3DSkewLineMidpoint(pos P1,pos P2,pos P3,pos P4)
  RETURN 0.5*(iv3DSkewLinePoint(P1,P2,P3,P4)+iv3DSkewLinePoint(P3,P4,P1,P2));
ENDFUNC

user_844416 · Answer

The shortest segment between two skew lines in 3D must be perpendicular to both (perpendicular projection is the shortest one).

So we have to get two points A and B that fulfill the next conditions:

A lies on P1..P2 line, so in parametric form using vector notation:

A = P1 + u*(P2-P1)

B lies on P3..P4 line, so in parametric form

B = P3 + v*(P4-P3)

vector AB is perpendicular to P1P2, so dot product is zero

(B-A).dot.(P2-P1) = 0

vector AB is perpendicular to P3P4, so dot product is zero

(B-A).dot.(P4-P3) = 0

Other stuff is vector algebra calculation to find u and v parameters (t in iv3DSkewLinePoint in your code).

I think that first expressions in the paper with dot products are simpler to calculate, than expressions with many vector product at the end - note wiki approach requires the only vector product calculation.

Seems yout long formula for t represents r1 formula from paper, but it is rather hard readable.

in addition:
Paul Bourke short article with codes
My code based on "Geometric Tools for Computer Graphics" book of D.Eberly (geometrictools.com)

Skew Lines Midpoint from 4 points?

1 Answers