Triangle Intersection in 3D Space

Could you explain more ?

• The 2 triangles are in an horizontal plane or not ?
• If it is 3D objects, intersection is probably a plan, not a line.
• What do you mean by their "vertecies" ?

A small drawing, even a hand sketch would really help.

Some suggestion.

Your triangles T1 (A1-B1-C1) and T2 (A2-B2-C2) are in Planes (P1 and P2).

In 3D, those planes are defined by equations:

• (P1) : a1 x + b1 y + c1 z + d1 = 0
• (P2) : a2 x + b2 y + c2 z + d2 = 0

You find a, b, c, d by resolving the system with 3 unknown, considering that each vertices of a triangle belongs to the plan.

For instance A1 is in P1: a1 xA1 + b1 yA1 + c1 zA1+ d1 = 0

You will have to consider 2 cases:

• d = 0 (the origin is in the plan : you completely resolve the system of 3 equations
• d ≠ 0 Origin is out of plan : divide consider d = 1 to resolve

Then the intersection of T1 and T2 is in fact intersection of P1 and P2, which is usually a line (may be empty if parallel plans or plan if plans are the same)

The intersection is given by:

• a1 x + b1 y + c1 z + d1 = 0
• a2 x + b2 y + c2 z + d2 = 0

solution may be expressed as:

• x = a t + xA
• y = b t + yA
• z = c t + zA

Hope that helps.

First, you didn't answer my question on what you are looking for exactly !

Back to your new question:

But how do I solve for those 4 vars.(a,b,c,d) with only one equation

You don't have one equation, you have 3.

Triangle T1 (with vertices A1, B1, C1) in plane P1:

• A1 is (xA1, yA1, zA1)
• B1 is (xB1, yB1, zB1)
• C1 is (xC1, yC1, zC1)

A1, B1, C1 are in plane P1; hence it satisfies 3 equations:

• a1 xA1 + b1 yA1 + c1 zA1+ d1 = 0
• a1 xB1 + b1 yB1 + c1 zB1+ d1 = 0
• a1 xC1 + b1 yC1 + c1 zC1+ d1 = 0

if Origine (0, 0, 0) is in plane P1, we have d1 = 0, hence 3 equations and 3 unknown:

• a1 xA1 + b1 yA1 + c1 zA1 = 0
• a1 xB1 + b1 yB1 + c1 zB1 = 0
• a1 xC1 + b1 yC1 + c1 zC1 = 0

if Origine (0, 0, 0) is not in plane P1, we have d1 ≠ 0, hence we can divide everything by d1, which is equivalent to d1 = 1 ; we have again 3 equations and 3 unknown:

• a1 xA1 + b1 yA1 + c1 zA1 + 1 = 0
• a1 xB1 + b1 yB1 + c1 zB1 + 1 = 0
• a1 xC1 + b1 yC1 + c1 zC1 + 1 = 0

There are tutorials to tell you how to solve this problem.

And Apple's example: https://developer.apple.com/documentation/accelerate/solving_systems_of_linear_equations_with_lapack

Once you have computed these coefficients, you have the equations of P1 and P2.

The intersection is given by satisfying both:

• a1 x + b1 y + c1 z + d1 = 0
• a2 x + b2 y + c2 z + d2 = 0

solution may be expressed as:

• x = a t + xA
• y = b t + yA
• z = c t + zA

For example, if we have

• (P1) : 2 x + 3 y - z + 2 = 0
• (P2) : x + y - 2z + 5 = 0

We get

• (P1) : z = 2 x + 3 y + 2
• (P2) : x + y - 2 (2 x + 3 y + 2) + 5 = 0 <=> -3x - 5y + 1 = 0

By multiplying first by 3 second by 2:

• (P1) : 3z = 6x + 9 y + 6
• (P2) : -6x - 10y + 2 = 0

summing the 2, we eliminate x

• (P1) : 6x + 9 y + 6 = 3z
• (P2) : -6x - 10y + 2 = 0

To get:

• (P1) : 6x + 9 y + 6 = 3z
• (P1) + (P2) : -y + 8 = 3z. Hence y = 8 - 3z

Thus

• (P1) : z = 2 x + 3 y + 2 => z = 2 x + 3 (8 - 3z) + 2 <=> 10 z = 2 x + 26 <=> x = 5 z - 13
• (P1) + (P2) : y = 8 - 3z

So finally, the parametric equation of the intersect:

• x = 5 z - 13
• y = - 3z + 8
• z = z

You can now transform this in code step by step (it is a bit tedious though).

But may be your problem is a bit simpler, depending what your original triangles are. For instance if one in in the plane z = 0, that will simplify a lot.

Note: this is to get the equation of intersection.

If you just want to test if a point in on the intersection, it is much easier ; check if it satisfies (P1) and (P2)

Please explain.

 If I am not wrong two triangles intersect at a line.

Are the 2 triangles in the same plane (you speak of 3D, so I guess they may be on different planes).

Could you post a simple drawing to illustrate what you are looking for.

This is what I would like to achive.

There are a couple papers I found that discuss this, I would try implementing one of those

https://stackoverflow.com/questions/1496215/triangle-triangle-intersection-in-3d-space

https://web.stanford.edu/class/cs277/resources/papers/Moller1997b.pdf

https://web.mst.edu/~chaman/home/pubs/2015WimoTriangleTrianglePublished.pdf

Thanks, that is much more clear.

Note that the 2 triangles could well not intersect.

I wrote a basic algorithm to find those intersection points which are inside the triangle.

Play with it by changing triangles (ABC and DEF) definition and tell me how it works.

I added a lot of comments to help you understand what's going on.

// Intersection of 2 triangles

struct Point3D {
    var x: Double
    var y: Double
    var z: Double
    
    static func vector(from a: Point3D, to b: Point3D) -> Vector3D {
        Vector3D(x: b.x - a.x, y: b.y - a.y, z: b.z - a.z)
    }
}

struct Vector3D {    // Same structure as Point3D
    var x: Double
    var y: Double
    var z: Double
    
    static func add(v1 u: Vector3D, v2 v: Vector3D) -> Vector3D {
        Vector3D(x: u.x + v.x, y: u.y + v.y, z: u.z + v.z)
    }
    
    static func prodScal(v1 u: Vector3D, v2 v: Vector3D) -> Double {
        u.x * v.x + u.y * v.y + u.z * v.z
    }
    
    static func prodVector(v1 u: Vector3D, v2 v: Vector3D) -> Vector3D {
        let wx = u.y*v.z - u.z*v.y
        let wy = u.z*v.x - u.x*v.z
        let wz = u.x*v.y - u.y*v.x
        return Vector3D(x: wx, y: wy, z: wz)
    }
    
    func multi(by d: Double) -> Vector3D {
        return Vector3D(x: d*self.x, y: d*self.y, z: d*self.z)
    }
}

/* Not needed so far
struct Triangle3D {
    var a : Point3D
    var b : Point3D
    var c : Point3D
}
*/

let pA = Point3D(x: 0, y: 5, z: 0)      // For simplicity in this test, ABC is in the plane (x, y) ; but put anything you want
let pB = Point3D(x: 5, y: 1, z: 0)
let pC = Point3D(x: 3, y: 1, z: 0)     

// let tr1 = Triangle3D(a: pA, b: pB, c: pC)

let pD = Point3D(x: 1, y: 2, z: 4)
let pE = Point3D(x: 3, y: 1, z: 6)
let pF = Point3D(x: 2, y: 4, z: -2)

// let tr2 = Triangle3D(a: pD, b: pE, c: pF)

// Compute "vertical" (perpendicular) vector AZ to ABC at point A

// AZ is the vector product : AB ^ AC, given by formula:
//     u =  AB          ^       v =  AC              w =  result
//  ux = [pB.x - pA.x]     vx = [pC.x - pA.x]      [uy vz - uz vy]
//  uy = [pB.y - pA.y]     vy = [pC.y - pA.y]      [uz vx - ux vz]
//  uz = [pB.z - pA.z]     vz = [pC.z - pA.z]      [ux vy - uy vx]
 
let u = Point3D.vector(from: pA, to: pB)     // vector AB
let v = Point3D.vector(from: pA, to: pC)     // vector AC

let w = Vector3D.prodVector(v1: u, v2: v)     // vector perpendicular to plane ABC

// a point inside ABC triangle is defined as:
// OM = α OA  β OB  γ OC (vectors)
//  with α + β + γ = 1
//  and α > 0  β > 0  γ > 0
// a point inside DEF triangle is defined as:
// OM = α' OD  β' OE  γ' OF (vectors)
//  with α' + β' + γ' = 1
//  and α' ≥ 0  β' ≥ 0  γ' ≥ 0
// Let's use D as O point:
// DM = α' DD +  β' DE +  γ' DF with α' + β' + γ' = 1 and α' ≥ 0  β' ≥ 0  γ' ≥ 0 ; α' no more used;  β' ≥ 0  γ' ≥ 0  β' + γ' ≤ 1
// AM = AD + DM

let ab = Point3D.vector(from: pA, to: pB)   // To be used later
let ac = Point3D.vector(from: pA, to: pC)   // To be used later
let ad = Point3D.vector(from: pA, to: pD)
let de = Point3D.vector(from: pD, to: pE)
let df = Point3D.vector(from: pD, to: pF)
// With our above definitions and calling vector AM as t : t = ad +  β de +  γ df

// We can now test all points in triangle DEF by looping on β and γ with a small increment
// And see when they are in the plane of ABC : we test that scalar product AM * w is zero (in fact less than a small epsilon value
// SIMD should be used here, but I do it rapidly

// For each point M which is in ABC plan (AM * w is zero), we need to check if it is inside the ABC triangle.
// That will be the case if AM = k AB + l AC with k ≥ 0, l ≥ 0 and k + l ≤ 1
// We have thus to resolve 3 equations with 2 unknown
// Noting AM is vector (t.x, t.y, t.z)
// (1) t.x = k ab.x + l ac.x
// (2) t.y = k ab.y + l ac.y
// (3) t.z = k ab.z + l ac.z
// However, because M is in ABC plane, those equations are redundant : if 2 are satisfied, the third will be.
// So we just have to solve by finding k and l :
// (1) t.x = k ab.x + l ac.x
// (2) t.y = k ab.y + l ac.y
// Which is simple: (delta = ab.x ac.y - ac.x ab.y is called the determinant)
//  k = (t.x ac.y - t.y ac.x) / (ab.x ac.y - ac.x ab.y)
//  l = (t.y ab.x - t.x ab.y) / (ab.x ac.y - ac.x ab.y)
// We shall test whether k and l are in correct bounds

let epsilon = 0.001

for beta in stride(from: 0.0, through: 1.0, by: 0.01) {
    for gamma in stride(from: 0.0, through: 1.0 - beta, by: 0.01) {
        let t1 = Vector3D.add(v1: de.multi(by: beta), v2: df.multi(by: gamma))
        let t = Vector3D.add(v1: ad, v2: t1)   // AM: t = ad +  β de +  γ df ; for a point M inside DEF
        let prod = Vector3D.prodScal(v1: t, v2: w)
        if abs(prod) < epsilon {    // M is in plane ABC
            // Is it inside triangle ABC ?
            let delta = ab.x * ac.y - ac.x * ab.y
            if delta == 0 {
                print("Determinant is zero ; use equations (2) and (3) instead")
                continue // skip this point, to not crash
            }
            let k = (t.x * ac.y - t.y * ac.x) / delta
            let l = (t.y * ab.x - t.x * ab.y) / delta
            print("k", k, "l", l)
            if k >= 0 && l >= 0 && k + l <= 1 {
                print("Point is inside ABC")
            }
        }
        
    }
}

Now that you can compute the intersection points, you can draw the line (between points defined by 2 extreme sets of valid values for k and l.

Or do anything else you need. Good continuation.

Sorry that I could not respond to your comment earlier. I had to learn for an exam. If someone finds some wrong math or wrong results please let me know.

import Foundation
import simd



// Dfinition of the two planes (Q,R)

let Q = (SIMD3<Double>(0,5,0.00000000000000001),SIMD3<Double>(5,1,0),SIMD3<Double>(3, 1, 0)) // first triangle

let R = (SIMD3<Double>(1,2,4),SIMD3<Double>(3,1,6),SIMD3<Double>(2,4,-2))

// That are the same points as in the answer above // In this example you can not set P.z = 0 because the line is paralle to xy



// normals

// normalized to get a better comperiosion

let NQ = normalize(cross((Q.1 - Q.0),(Q.2 - Q.0)))

let NR = normalize(cross((R.1 - R.0),(R.2 - R.0)))



// Plane equtions

// x -> a , y -> b, z -> c // That is the normal to plane equtions variable convertion ((a,b,c) are used normally)

// D of first plane

let DQ = -(NQ.x*Q.0.x + NQ.y*Q.0.y + NQ.z*Q.0.z)

// D of second plane

let DR = -(NR.x*R.0.x + NR.y*R.0.y + NR.z*R.0.z)



// Find P on Vektor V

// P.z = 0

let x = (NQ.y*DR - DQ*NR.x)/(NQ.x*NR.y - NQ.y*NR.x) // took to plane equtions and removed z and solved for x and y

let y = -(NR.x*x+DR) / NR.y // computes y of the Point P

let P = SIMD3<Double>(x,y,0) // That is a point in the Intesectionline of the two planes



//Linevector

let VI = cross(NQ,NR) // VI is normal on both palne normlas



// Lines on the Vektor of triangle Q

let IntersectionLine = [(t:computesides(P,VI,Q.0,Q.1),first: true),(t:computesides(P,VI,Q.0,Q.2),first: true),(t:computesides(P,VI,Q.2,Q.1),first: true),(t:computesides(P,VI,R.0,R.1),first: false),(t:computesides(P,VI,R.0,R.2),first: false),(t:computesides(P,VI,R.2,R.1),first: false)]

var UnwrapedLine = IntersectionLine.filter { $0.t != nil} // No nil ones left

if UnwrapedLine.count != 4 { print("No Intersection")}

UnwrapedLine.sort { $0.t! < $1.t!} // There should be no nils left

if UnwrapedLine[0].first == UnwrapedLine[1].first {

   print("Linesegments do not intersect")

} else {

   print("Line: \(P+VI*UnwrapedLine[0].t!) to \(P+VI*UnwrapedLine[1].t!) )")

}



func computesides(_ P: SIMD3<Double>,_ IV: SIMD3<Double>,_ VS: SIMD3<Double>,_ VE: SIMD3<Double>) -> Double? { // This func returns the distance to the Point P in the form of a Vector

   let v = VE - VS // v for Vector

   let s = (P.x * IV.y-P.y*IV.x+IV.x*VS.y-IV.y*VS.x)/(IV.y*v.x-IV.x*v.y) // Scalar shapeedge of v

   if s > 1 || s < 0 {

       print("No Intersection with the plane")

       return nil

   }

   let t = (VS.x+v.x*s-P.x)/IV.x // it can be solved with x,y or z//P+IV*t = VS+v*s // (VS+v*s-P)/IV

return t

}

Thanks for your code above.

If someone finds some wrong math or wrong results please let me know.

¡¿¡¿ You should be the one to tell if that works !?!?

Not very fair you didn't even test or feedback to the code I wrote for your question…