intersections.lua

function overlapsOBB_OBB(c1, x1, y1, z1, c2, x2, y2, z2)
  local cc = c1 - c2
  local d11, d12, d13 = abs(x1:dot(x2)), abs(x1:dot(y2)), abs(x1:dot(z2))
  local d21, d22, d23 = abs(y1:dot(x2)), abs(y1:dot(y2)), abs(y1:dot(z2))
  local d31, d32, d33 = abs(z1:dot(x2)), abs(z1:dot(y2)), abs(z1:dot(z2))

  return abs(cc:dot(x1))-d11-d12-d13<=x1:squaredLength() and abs(cc:dot(y1))-d21-d22-d23<=y1:squaredLength()
     and abs(cc:dot(z1))-d31-d32-d33<=z1:squaredLength() and abs(cc:dot(x2))-d11-d21-d31<=x2:squaredLength()
     and abs(cc:dot(y2))-d12-d22-d32<=y2:squaredLength() and abs(cc:dot(z2))-d13-d23-d33<=z2:squaredLength()
end

-- untested
function containsOBB_OBB(c1, x1, y1, z1, c2, x2, y2, z2)
  local cc = c1 - c2
  return abs(cc:dot(x1))+abs(x1:dot(x2))+abs(x1:dot(y2))+abs(x1:dot(z2))<=x1:squaredLength()
     and abs(cc:dot(y1))+abs(y1:dot(x2))+abs(y1:dot(y2))+abs(y1:dot(z2))<=y1:squaredLength()
     and abs(cc:dot(z1))+abs(z1:dot(x2))+abs(z1:dot(y2))+abs(z1:dot(z2))<=z1:squaredLength()
end

function overlapsOBB_Sphere(c1, x1, y1, z1, c2, r2)
  local cc = c1 - c2
  local x1len, y1len, z1len = x1:length(), y1:length(), z1:length()
  local ccx, ccy, ccz = abs(cc:dot(x1)), abs(cc:dot(y1)), abs(cc:dot(z1))

  return ccx<=x1len*(x1len+r2) and ccy<=y1len*(y1len+r2) and ccz<=z1len*(z1len+r2)
    and (ccx<=x1len*x1len or ccy<=y1len*y1len or ccz<=z1len*z1len or
      square(ccx/(x1len+1e-30)-x1len)+square(ccy/(y1len+1e-30)-y1len)+square(ccz/(z1len+1e-30)-z1len)<=r2*r2)
end

function overlapsOBB_Plane(c1, x1, y1, z1, plpos, pln)
  return abs((c1 - plpos):dot(pln))<=abs(x1:dot(pln))+abs(y1:dot(pln))+abs(z1:dot(pln))
end

function containsOBB_Sphere(c1, x1, y1, z1, c2, r2)
  local cc = c1 - c2
  local x1len, y1len, z1len = x1:length(), y1:length(), z1:length()
  return abs(cc:dot(x1))<=x1len*(x1len-r2) and abs(cc:dot(y1))<=y1len*(y1len-r2) and abs(cc:dot(z1))<=z1len*(z1len-r2)
end

function containsSphere_OBB(c1, r1, c2, x2, y2, z2)
  local cc = c1 - c2
  local ccx1, cc_x2, y2z2, y2_z2 = cc+x2, cc-x2, y2+z2, y2-z2
  return max((ccx1+y2z2):squaredLength(), (ccx1+y2_z2):squaredLength(), (ccx1-y2_z2):squaredLength(), (ccx1-y2z2):squaredLength(),
      (cc_x2+y2z2):squaredLength(), (cc_x2+y2_z2):squaredLength(), (cc_x2-y2_z2):squaredLength(), (cc_x2-y2z2):squaredLength())<=r1*r1
end

function containsOBB_point(c1, x1, y1, z1, p)
  local cc = c1 - p
  return abs(cc:dot(x1))<=x1:squaredLength() and abs(cc:dot(y1))<=y1:squaredLength() and abs(cc:dot(z1))<=z1:squaredLength()
end

function containsEllipsoid_Point(c1, x1, y1, z1, p)
  local cc = p - c1
  local x, y, z = cc:dot(x1), cc:dot(y1), cc:dot(z1)
  local a2, b2, c2 = x1:squaredLength(), y1:squaredLength(), z1:squaredLength()
  a2, b2, c2 = a2*a2, b2*b2, c2*c2
  local b2c2 = b2*c2
  return x*x*b2c2 + a2*(y*y*c2 + z*z*b2) <= a2*b2c2
end

function constainsCylinder_Point(cposa, cposb, cR, p)
  local xnorm, r2 = p:xnormSquaredDistanceToLineSegment(cposa, cposb)
  return xnorm >=0 and xnorm <= 1 and r2 <= cR*cR
end

function altitudeOBB_Plane(c1, x1, y1, z1, plpos, pln)
  return (c1 - plpos):dot(pln)+abs(x1:dot(pln))+abs(y1:dot(pln))+abs(z1:dot(pln))
end

-- returns signed distance of plane on the ray
function intersectsRay_Plane(rpos, rdir, plpos, pln)
  return min((plpos - rpos):dot(pln) / rdir:dot(pln), math.huge)
end

-- hit: minhit < maxhit, inside: minhit < 0
function intersectsRay_Sphere(rpos, rdir, cpos, cr)
  local rcpos = cpos - rpos
  local dcr = rdir:dot(rcpos)
  local s = dcr*dcr - rcpos:squaredLength() + cr*cr
  if s < 0 then return math.huge, math.huge end
  s = sqrt(s)
  return dcr - s, dcr + s
end

function intersectsRay_Ellipsoid(rpos, rdir, c1, x1, y1, z1)
  local invx1, invy1, invz1 = 1 / (x1:squaredLength() + 1e-30), 1 / (y1:squaredLength() + 1e-30), 1 / (z1:squaredLength() + 1e-30)
  local cc = rpos - c1
  local pM = vec3(cc:dot(x1)*invx1, cc:dot(y1)*invy1, cc:dot(z1)*invz1)
  local dirM = vec3(rdir:dot(x1)*invx1, rdir:dot(y1)*invy1, rdir:dot(z1)*invz1)

  local a, b, c = dirM:squaredLength(), 2*pM:dot(dirM), pM:squaredLength() - 1
  local d = b*b - 4*a*c
  if d < 0 then return math.huge, math.huge end
  d = -b -sign(b)*sqrt(d)
  local r1, r2 = 0.5*d / a, 2*c / d
  return min(r1, r2), max(r1,r2)
end

function intersectsRay_Cylinder(rpos, rdir, cposa, cposb, cR)
  local rca, cba = cposa - rpos, cposb - cposa
  local cpnorm = cba:normalized()
  local cp = rca:projectToOriginPlane(cpnorm)
  local rdp = rdir:projectToOriginPlane(cpnorm)
  local minhit, maxhit = intersectsRay_Sphere(vec3(0,0,0), rdp:normalized(), cp, cR)
  local invrdplen = 1 / (rdp:length() + 1e-30)
  minhit, maxhit = minhit * invrdplen, maxhit * invrdplen
  local plhita, plhitb = intersectsRay_Plane(rpos, rdir, cposa, cpnorm), intersectsRay_Plane(rpos, rdir, cposb, cpnorm)
  minhit, maxhit = max(minhit, min(plhita, plhitb)), min(maxhit, max(plhita, plhitb))
  return (minhit <= maxhit and minhit or math.huge), maxhit
end

function intersectsRay_OBB(rpos, rdir, c1, x1, y1, z1)
  local rposc1 = c1 - rpos
  local rposc1x1, x1sq, invrdirx1 = rposc1:dot(x1), x1:squaredLength(), 1 / rdir:dot(x1)
  local dx1, dx2 = (rposc1x1 - x1sq) * invrdirx1, min((rposc1x1 + x1sq) * invrdirx1, math.huge)
  local rposc1y1, y1sq, invrdiry1 = rposc1:dot(y1), y1:squaredLength(), 1 / rdir:dot(y1)
  local dy1, dy2 = (rposc1y1 - y1sq) * invrdiry1, min((rposc1y1 + y1sq) * invrdiry1, math.huge)
  local rposc1z1, z1sq, invrdirz1 = rposc1:dot(z1), z1:squaredLength(), 1 / rdir:dot(z1)
  local dz1, dz2 = (rposc1z1 - z1sq) * invrdirz1, min((rposc1z1 + z1sq) * invrdirz1, math.huge)

  local minhit, maxhit = max(min(dx1, dx2), min(dy1, dy2), min(dz1, dz2)), min(max(dx1, dx2), max(dy1, dy2), max(dz1, dz2))
  return (minhit <= maxhit and minhit or math.huge), maxhit
end

-- returns hit distance, barycentric x, y
function intersectsRay_triangle(rpos, rdir, a, b, c)
  local ca, bc = c - a, b - c
  local norm = ca:cross(bc)
  local rposc = rpos - c
  local pOnTri = rposc:dot(norm) / rdir:dot(norm)
  if pOnTri <= 0 then
    local pacnorm = (rposc - rdir * pOnTri):cross(norm)
    local bx, by = bc:dot(pacnorm), ca:dot(pacnorm)
    local normSq = norm:squaredLength() + 1e-30
    if min(bx, by) >= 0 and bx + by <= normSq then
      return -pOnTri, bx / normSq, by / normSq
    end
  end
  return math.huge, -1, -1
end