# how to find the vector intersection in a collection of triangulated 3D points?

I'm trying to calculate the 3D point coordinates for the intersection of a vector with a plane. I have been using the following code to do this:

PointPlaneIntersection[{{x1_, y1_, z1_}, {x2_, y2_, z2_}, {x3_, y3_,
z3_}}, {{x4_, y4_, z4_}, {x5_, y5_, z5_}}] :=  Module[{t = -Det[{{1, 1, 1, 1}, {x1, x2, x3, x4}, {y1, y2, y3,
y4}, {z1, z2, z3, z4}}]/
Det[{{1, 1, 1, 0}, {x1, x2, x3, x5 - x4}, {y1, y2, y3,
y5 - y4}, {z1, z2, z3, z5 - z4}}]},  {x4 + t (x5 - x4), y4 + t (y5 - y4), z4 + t (z5 - z4)}]


which I found at: http://mathworld.wolfram.com/Line-PlaneIntersection.html

I have the vector and a point of origin for this vector:

vec={1.11593, 0.0607521, 0.125472};pt={-45.0446, 0.0900457, -1.42599}


However, rather than having a know plane defined by 3 points, as the PointPlaneIntersection code needs, I have a group of planes/triangles defined by 3D points:

triangles={{{-46.4755, -0.0930855, -2.55178}, {-46.6708, -1.22548, -2.25304},{-46.9169, -0.664019, -0.916458}}, {{-46.4736,  0.544871, -1.66483}, {-46.4755, 0.0930855, -2.55178}, {-46.9169,-0.664019, -0.916458}}, {{-46.0117, 0.605676, -3.13173}, {-46.4755, -0.0930855, -2.55178}, {-46.4736, 0.544871, -1.66483}}, {{-45.6652, 1.74667, -1.92923}, {-46.0117, 0.605676, -3.13173}, {-46.4736, 0.544871, -1.66483}}, {{-45.2907, 1.39215, -3.80111}, {-46.0117, 0.605676, -3.13173}, {-45.6652, 1.74667, -1.92923}}, {{-46.4736, 0.544871, -1.66483}, {-46.9169, -0.664019, -0.916458}, {-46.8575, -0.142585, 0.141953}}, {{-46.2786, 1.32328, -0.512178}, {-46.4736, 0.544871, -1.66483}, {-46.8575, -0.142585, 0.141953}}, {{-45.6652, 1.74667, -1.92923}, {-46.4736, 0.544871, -1.66483}, {-46.2786, 1.32328, -0.512178}}};


This can be visualised:

Show[Graphics3D[{Yellow, Specularity[White, 20], Directive[Opacity[0.75]], EdgeForm[{Gray, Thin}], Polygon[triangles], Red, AbsolutePointSize[8], Point[pt], Blue, Thick, Arrowheads[0.02], Arrow[{pt, (pt + (vec))}]}], PlotRange -> All, Lighting -> "Neutral", Boxed -> False, Axes -> False, Ticks -> None, AxesLabel -> {"X", "Y", "Z"}, AxesStyle -> Gray, BoxRatios -> Automatic, AspectRatio -> Automatic, ImageSize -> 650, ViewPoint -> vec*10]


Also, note that in the graphic produced the vector direction is away from the collection of triangles. I'm looking for a solution where it doesn't matter whether the vector is towards or away from the triangles...ie the code will effectively flip the vector round.

Can anyone suggest a solution for selecting which triangle/plane the vector passes through?

Or alternatively, could anyone suggest a good way to fit a surface (eg a bezier surface) to the points, then find the intersection of the vector with the surface?

thank you for reading my problem,

Dime

• Do you have Version 10? This has some inbuilt functionality. 2D examples here: mathematica.stackexchange.com/q/56299/131. Docs: reference.wolfram.com/language/ref/Solve.html (see "Scope / Geometric Regions") Sep 11, 2014 at 11:05
• note this can be done in closed form: use your PointPlaneIntersection on each triangle, then the resulting 2d insideness check has an analytic formula. This probably is a tad faster than belisarius approach. Sep 11, 2014 at 18:41
• thanks for your comments. At the moment I'm having trouble getting version 10 (site licence issue)... @george2079. This is an interesting approach that I hadn't thought of. Do you think it would be quicker than belisarius' solution if there were a lot of triangles? would you approach this by converting something like the following code for 3D? testpoint2[poly_, pt_] := GraphicsMeshInPolygonQ[poly, pt] testpoint2[{{-1, 0}, {0, 1}, {1, 0}}, {1/3, 1/3}] (True)
– dime
Sep 12, 2014 at 0:15

(* equation for your vector*)
lin[a_, vec_, pt_] := a vec + pt

(* Param eq. for the surface of a triangle *)
c[s_, t_, v_] := v[[1]] (1 - t) s + v[[2]] t s + v[[3]] (1 - s)

(*which triangle and the parameters*)
s1 = NSolve[{lin[a, vec, pt] == c[s, t, #], 0 < s < 1, 0 < t < 1}, {a, s,  t}] & /@ triangles

(* {{}, {{a -> -1.42384, s -> 0.640281,  t -> 0.368053}}, {}, {}, {}, {}, {}, {}} *)

So it's in the second triangle, and the point of intersection is:

int = lin[a, vec, pt] /. s1[[2]]
(* {{-46.6335, 0.00354414, -1.60464}} *)

Graphics3D[{Opacity[.3], Polygon@ triangles,
Green,  Polygon[triangles[[2]]],
Blue,   Arrow[{pt, lin[a, vec, pt]}] /. s1[[2, 1]]},
Opacity[1], Red, PointSize[Large], Point@int,
Boxed -> False]


• thank you Belisarius. That does exactly what I needed! I also added the following so that the triangle number (ie number 2) doesn't have to be manually input: triagNum = Position[s1, a][[1, 1]]; int = Flatten[lin[a, vec, pt] /. s1[[triagNum]], 1]; which gives the same, Graphics3D[{Opacity[.3], Polygon[triangles], Green, Polygon[triangles[[triagNum]]], Blue, Arrow[{pt, lin[a, vec, pt]}] /. s1[[2, 1]], Opacity[1], Red, PointSize[Large], Point[int]}, Boxed -> False]
– dime
Sep 11, 2014 at 23:59
• @dime Glad to help. I used the constant 2 to kept the code easier to understand, but you're right Sep 12, 2014 at 0:11