type TMyPoint = record case boolean of false: (X, Y, Z: double); true: (arr: array[0 .. 2] of double); end; function pSUB(const v1, v2: TMyPoint): TMyPoint; begin result.arr[0] := v1.arr[0] - v2.arr[0]; result.arr[1] := v1.arr[1] - v2.arr[1]; result.arr[2] := v1.arr[2] - v2.arr[2]; end; function pCROSS(const v1, v2: TMyPoint): TMyPoint; begin result.arr[0] := v1.arr[1]*v2.arr[2]-v1.arr[2]*v2.arr[1]; result.arr[1] := v1.arr[2]*v2.arr[0]-v1.arr[0]*v2.arr[2]; result.arr[2] := v1.arr[0]*v2.arr[1]-v1.arr[1]*v2.arr[0]; end; function pDOT(const v1, v2: TMyPoint): double; begin result := v1.arr[0]*v2.arr[0]+v1.arr[1]*v2.arr[1]+v1.arr[2]*v2.arr[2]; end; procedure pSORT(var a, b: double); var T: double; begin if a > b then begin T := a; a := b; b := T; end; end; {$define USE_EPSILON_TEST} const EPSILON = 0.000001; function coplanar_tri_tri(N, V0, V1, V2, U0, U1, U2: TMyPoint): boolean; var i0, i1: integer; function EDGE_AGAINST_TRI_EDGES(const V0,V1,U0,U1,U2: TMyPoint): boolean; var Ax,Ay,Bx,By,Cx,Cy,e,d,f: double; function EDGE_EDGE_TEST(const V0,U0,U1: TMyPoint): boolean; begin result := false; Bx:=U0.arr[i0]-U1.arr[i0]; By:=U0.arr[i1]-U1.arr[i1]; Cx:=V0.arr[i0]-U0.arr[i0]; Cy:=V0.arr[i1]-U0.arr[i1]; f:=Ay*Bx-Ax*By; d:=By*Cx-Bx*Cy; if ((f > 0) and (d >= 0) and (d <= f)) or ((f < 0) and (d <= 0) and (d >= f)) then begin e := Ax*Cy-Ay*Cx; if f > 0 then begin if(e >= 0) and (e <= f) then result := true; end else begin if (e <= 0) and (e >= f) then result := true; end; end; end; begin result := false; Ax:=V1.arr[i0]-V0.arr[i0]; Ay:=V1.arr[i1]-V0.arr[i1]; { test edge U0,U1 against V0,V1 } if EDGE_EDGE_TEST(V0,U0,U1) then begin result := true; exit; end; { test edge U1,U2 against V0,V1 } if EDGE_EDGE_TEST(V0,U1,U2) then begin result := true; exit; end; { test edge U2,U1 against V0,V1 } if EDGE_EDGE_TEST(V0,U2,U0) then begin result := true; exit; end; end; function POINT_IN_TRI(const V0,U0,U1,U2: TMyPoint): boolean; var a,b,c,d0,d1,d2: double; begin result := false; { is T1 completly inside T2? check if V0 is inside tri(U0,U1,U2) } a := U1.arr[i1]-U0.arr[i1]; b := -(U1.arr[i0]-U0.arr[i0]); c := -a*U0.arr[i0]-b*U0.arr[i1]; d0 := a*V0.arr[i0]+b*V0.arr[i1]+c; a := U2.arr[i1]-U1.arr[i1]; b := -(U2.arr[i0]-U1.arr[i0]); c := -a*U1.arr[i0]-b*U1.arr[i1]; d1 := a*V0.arr[i0]+b*V0.arr[i1]+c; a := U0.arr[i1]-U2.arr[i1]; b := -(U0.arr[i0]-U2.arr[i0]); c := -a*U2.arr[i0]-b*U2.arr[i1]; d2 := a*V0.arr[i0]+b*V0.arr[i1]+c; if d0*d1 > 0.0 then begin if d0*d2 > 0.0 then begin result := true; exit end; end; end; var A: array[0 .. 2] of double; begin { first project onto an axis-aligned plane, that maximizes the area of the triangles, compute indices: i0,i1. } A[0] := abs(N.arr[0]); A[1] := abs(N.arr[1]); A[2] := abs(N.arr[2]); if A[0] > A[1] then begin if A[0]>A[2] then begin i0:=1; { A[0] is greatest } i1:=2; end else begin i0:=0; { A[2] is greatest } i1:=1; end end else begin { A[0]<=A[1] } if A[2] > A[1] then begin i0:=0; { A[2] is greatest } i1:=1; end else begin i0:=0; { A[1] is greatest } i1:=2; end end; { test all edges of triangle 1 against the edges of triangle 2 } if EDGE_AGAINST_TRI_EDGES(V0,V1,U0,U1,U2) then begin result := true; exit; end; if EDGE_AGAINST_TRI_EDGES(V1,V2,U0,U1,U2) then begin result := true; exit; end; if EDGE_AGAINST_TRI_EDGES(V2,V0,U0,U1,U2) then begin result := true; exit; end; { finally, test if tri1 is totally contained in tri2 or vice versa } if POINT_IN_TRI(V0,U0,U1,U2) then begin result := true; exit; end; if POINT_IN_TRI(U0,V0,V1,V2) then begin result := true; exit; end; result := false; end; procedure ISECT(const VV0,VV1,VV2,D0,D1,D2: double; var isect0,isect1: double); begin isect0 := VV0+(VV1-VV0)*D0/(D0-D1); isect1 := VV0+(VV2-VV0)*D0/(D0-D2); end; function COMPUTE_INTERVALS(const VV0,VV1,VV2,D0,D1,D2,D0D1,D0D2: double; var isect0, isect1: double): boolean; begin result := false; if D0D1 > 0.0 then begin { here we know that D0D2<=0.0 that is D0, D1 are on the same side, D2 on the other or on the plane } ISECT(VV2,VV0,VV1,D2,D0,D1,isect0,isect1); end else if D0D2 > 0.0 then begin { here we know that d0d1<=0.0 } ISECT(VV1,VV0,VV2,D1,D0,D2,isect0,isect1); end else if(D1*D2>0.0) or (D0 <> 0.0) then begin { here we know that d0d1<=0.0 or that D0!=0.0 } ISECT(VV0,VV1,VV2,D0,D1,D2,isect0,isect1); end else if D1 <> 0.0 then begin ISECT(VV1,VV0,VV2,D1,D0,D2,isect0,isect1); end else if D2 <> 0.0 then begin ISECT(VV2,VV0,VV1,D2,D0,D1,isect0,isect1); end else begin { triangles are coplanar } result := true; end; end; function tri_tri_intersect(V0, V1, V2, U0, U1, U2: TMyPoint): boolean; var E1, E2: TMyPoint; N1, N2: TMyPoint; d1, d2: double; du0,du1,du2,dv0,dv1,dv2: double; du0du1,du0du2,dv0dv1,dv0dv2: double; D: TMyPoint; index: integer; vp0,vp1,vp2, up0,up1,up2, b,c,max: double; isect1, isect2: array[0 .. 1] of double; begin { compute plane equation of triangle(V0,V1,V2) } E1 := pSUB(V1,V0); E2 := pSUB(V2,V0); N1 := pCROSS(E1,E2); d1 := -pDOT(N1,V0); { plane equation 1: N1.X+d1=0 } { put U0,U1,U2 into plane equation 1 to compute signed distances to the plane } du0 := pDOT(N1,U0) + d1; du1 := pDOT(N1,U1) + d1; du2 := pDOT(N1,U2) + d1; { coplanarity robustness check } {$ifdef USE_EPSILON_TEST} if abs(du0) < EPSILON then du0 := 0.0; if abs(du1) < EPSILON then du1 := 0.0; if abs(du2) < EPSILON then du2 := 0.0; {$endif} du0du1 := du0*du1; du0du2 := du0*du2; result := false; { same sign on all of them + not equal 0 ? } if (du0du1 > 0.0) and (du0du2 > 0.0) then exit; { no intersection occurs } { compute plane of triangle (U0,U1,U2) } E1 := pSUB(U1,U0); E2 := pSUB(U2,U0); N2 := pCROSS(E1,E2); d2 := -pDOT(N2,U0); { plane equation 2: N2.X+d2=0 } { put V0,V1,V2 into plane equation 2 } dv0 := pDOT(N2,V0) + d2; dv1 := pDOT(N2,V1) + d2; dv2 := pDOT(N2,V2) + d2; {$ifdef USE_EPSILON_TEST} if abs(dv0) < EPSILON then dv0 := 0.0; if abs(dv1) < EPSILON then dv1 := 0.0; if abs(dv2) < EPSILON then dv2 := 0.0; {$endif} dv0dv1 := dv0*dv1; dv0dv2 := dv0*dv2; { same sign on all of them + not equal 0 ? } if (dv0dv1 > 0.0) and (dv0dv2 > 0.0) then exit; { no intersection occurs } { compute direction of intersection line } D := pCROSS(N1,N2); { compute and index to the largest component of D } max := abs(D.arr[0]); index := 0; b := abs(D.arr[1]); c := abs(D.arr[2]); if b > max then begin max := b; index := 1 end; if c > max then begin max := c; index := 2 end; { this is the simplified projection onto L } vp0 := V0.arr[index]; vp1 := V1.arr[index]; vp2 := V2.arr[index]; up0 := U0.arr[index]; up1 := U1.arr[index]; up2 := U2.arr[index]; { compute interval for triangle 1 } if COMPUTE_INTERVALS(vp0,vp1,vp2,dv0,dv1,dv2,dv0dv1,dv0dv2,isect1[0],isect1[1]) then begin result := coplanar_tri_tri(N1,V0,V1,V2,U0,U1,U2); exit; end; { compute interval for triangle 2 } if COMPUTE_INTERVALS(up0,up1,up2,du0,du1,du2,du0du1,du0du2,isect2[0],isect2[1]) then begin result := coplanar_tri_tri(N1,V0,V1,V2,U0,U1,U2); exit; end; pSORT(isect1[0],isect1[1]); pSORT(isect2[0],isect2[1]); if(isect1[1] < isect2[0]) or (isect2[1] < isect1[0]) then exit; result := true; end; begin end.