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.
