3d line-intersection code not working properly

Question

I created this piece of code to get the intersection of two 3d line-segments.

Unfortunately the result of this code is inaccurate, the intersection-point is not always on both lines.

I am confused and unsure what I'm doing wrong.

Here is my code:

--dir = direction
--p1,p2 = represents the line
function GetIntersection(dirStart, dirEnd, p1, p2)
   local s1_x, s1_y, s2_x, s2_y =  dirEnd.x - dirStart.x, dirEnd.z - dirStart.z, p2.x - p1.x, p2.z - p1.z
   local div = (-s2_x * s1_y) + (s1_x * s2_y)

   if div == 0 then return nil end
   local s = (-s1_y * (dirStart.x - p1.x) + s1_x * (dirStart.z - p1.z)) / div
   local t = ( s2_x * (dirStart.z - p1.z) - s2_y * (dirStart.x - p1.x)) / div

   if (s >= 0 and s <= 1 and t >= 0 and t <= 1) and (Vector(dirStart.x + (t * s1_x), 0, dirStart.z + (t * s1_y)) or nil) then
      local v = Vector(dirStart.x + (t * s1_x),0,dirStart.z + (t * s1_y))
      return v
   end
end

Answer 1

This is example of Delphi code to find a distance between two skew lines in 3D. For your purposes it is necessary to check that result if small enough value (intersection does exist), check that s and t parameters are in range 0..1, then calculate point using parameter s

Math of this approach is described in 'the shortest line...' section of Paul Bourke page

VecDiff if vector difference function, Dot id scalar product function

function LineLineDistance(const L0, L1: TLine3D; var s, t: Double): Double;
var
  u: TPoint3D;
  a, b, c, d, e, det, invdet:Double;
begin
  u := VecDiff(L1.Base, L0.Base);
  a := Dot(L0.Direction, L0.Direction);
  b := Dot(L0.Direction, L1.Direction);
  c := Dot(L1.Direction, L1.Direction);
  d := Dot(L0.Direction, u);
  e := Dot(L1.Direction, u);
  det := a * c - b * b;
  if det < eps then   
    Result := -1
  else begin
    invdet := 1 / det;
    s := invdet * (b * e - c * d);
    t := invdet * (a * e - b * d);

    Result := Distance(PointAtParam(L0, s), PointAtParam(L1, t));
  end;
end;

Answer 2

As far as I can tell your code is good. I've implemented this in javascript at https://jsfiddle.net/SalixAlba/kkrc9kcf/

and it seems to work for all the cases I can think of. The only changes I've done is to change things to work in javascript rather than lua. The final condition was commented out

function GetIntersection(dirStart, dirEnd, p1, p2) {
   var s1_x = dirEnd.x - dirStart.x;
   var s1_y = dirEnd.z - dirStart.z;
   var s2_x = p2.x - p1.x;
   var s2_y = p2.z - p1.z;
   var div = (-s2_x * s1_y) + (s1_x * s2_y);

   if (div == 0) 
     return new Vector(0,0);

   var s = (-s1_y * (dirStart.x - p1.x) + s1_x * (dirStart.z - p1.z)) / div;
   var t = ( s2_x * (dirStart.z - p1.z) - s2_y * (dirStart.x - p1.x)) / div;

   if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
      //and (Vector(dirStart.x + (t * s1_x), 0, dirStart.z + (t * s1_y)) or nil) then
      var v = new Vector(
        dirStart.x + (t * s1_x),
        dirStart.z + (t * s1_y));
      return v;
   }
   return new Vector(0,0);
}

Mathmatically it makes sense. If A,B and C,D are your two lines. Let s1 = B-A, s2 = C-D. A point of the line AB is given by A + t s1 and a point on the line CD is given by C + s s2. For an intersection we require

A + t s1 = C + s s2

or

(A-C) + t s1 = s s2

You two formula for s, t are found by taking the 2D cross product with each of the vectors s1 and s2

(A-C)^s1 + t s1^s1 = s s2^s1 (A-C)^s2 + t s1^s2 = s s2^s2

recalling s1^s1=s2^s2=0 and s2^s1= - s1^s2 we get

(A-C)^s1 = s s2^s1 (A-C)^s2 + t s1^s2 = 0

which can be solved to get s and t. This matches your equations.