1
$\begingroup$

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]

enter image description here

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

Improve this question
$\endgroup$
3
  • $\begingroup$ 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") $\endgroup$
    – Yves Klett
    Sep 11, 2014 at 11:05
  • $\begingroup$ 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. $\endgroup$
    – george2079
    Sep 11, 2014 at 18:41
  • $\begingroup$ 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) $\endgroup$
    – dime
    Sep 12, 2014 at 0:15

1 Answer 1

Reset to default
1
$\begingroup$ 
(* 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]

Mathematica graphics

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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] $\endgroup$
    – dime
    Sep 11, 2014 at 23:59
  • $\begingroup$ @dime Glad to help. I used the constant 2 to kept the code easier to understand, but you're right $\endgroup$
    – Dr. belisarius
    Sep 12, 2014 at 0:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.