Re: winding number

*To*: mathgroup at smc.vnet.net*Subject*: [mg25405] Re: winding number*From*: Martin Kraus <Martin.Kraus at informatik.uni-stuttgart.de>*Date*: Fri, 29 Sep 2000 01:06:48 -0400 (EDT)*Organization*: Institut fuer Informatik, Universitaet Stuttgart*References*: <8qho3l$j8s@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Will Self wrote: > > There has been some discussion of finding out whether a point in the > plane is enclosed in a given polygon. Here is some code that will > calculate the winding number of a polygonal closed curve about a given > point. The curve should be given as a list of vertices (the first > vertex does not have to be repeated at the end of the list). The list > of vertices is called vList in the code. This is really from Cauchy's > theorem with the integration along segments already done. > > Will Self > > --------------------------------------------------- > > windingNumber[point_, vList_] := Module[ > {g, angle}, > g[t_] := Mod[t + Pi, 2Pi] - Pi; > angle[p_, a1_, a2_] := g[ArcTan @@ (a2 - p) - ArcTan @@ (a1 - p)]; > (Plus @@ MapThread [ angle[point, #1, #2] & , {vList, RotateLeft[vList]} > ])/2/Pi] > > Sent via Deja.com http://www.deja.com/ > Before you buy. Hi MathGroup! I have to admit that this is probably the "best" way of testing whether a point is inside a polygon or not (in Mathematica). However, it is not easily generalized to the 3d case of a point inside a polyhedron! (At least I think it is not easy.) Here is my solution: The basic idea is to take a list of well oriented triangles (ordering of vertices anticlockwise as seen from outside the polyhedron) and calculating the solid angle (or spherical excess) of each triangle with respect to the point in question after adding signs (depending whether the front or back face is visible). If the point is inside the polyhedron, the sum of the excesses should be a 4 pi or (-4 pi depending on our choice of signs); otherwise it should be 0. An implementation is this: tangentVector::usage = "tangentVector[a, b, p] gives the tangent vector in a parallel to the \ line from a to b on a unit sphere with center p."; tangentVector[a_, b_, p_] := Module[{na, nb, t}, na = N[(a - p)/Sqrt[(a - p).(a - p)]]; nb = N[(b - p)/Sqrt[(b - p).(b - p)]]; t = (nb - na) - na(na.(nb - na)); N[t/Sqrt[t.t]]] sphericalExcess::usage = "sphericalExcess[a, b, c, p] gives the spherical exess of the triangle a, \ b, c with respect to a sphere with center p."; sphericalExcess[a_, b_, c_, p_] := Module[{tab, tac, tba, tbc, tca, tcb, aa, ab, ac}, tab = tangentVector[a, b, p]; tac = tangentVector[a, c, p]; tba = tangentVector[b, a, p]; tbc = tangentVector[b, c, p]; tca = tangentVector[c, a, p]; tcb = tangentVector[c, b, p]; aa = ArcCos[tab.tac]; ab = ArcCos[tba.tbc]; ac = ArcCos[tca.tcb]; aa + ab + ac - Pi] signedSphericalExcess::usage := "signedSphericalExcess[a, b, c, p] gives the \ spherical excess of the triangle a, b, c with respect to a sphere with center \ p if a, b, c is in anti-clockwise order as seen from p, and the negative \ value if a, b, c is in clockwise order." signedSphericalExcess[a_, b_, c_, p_] := Module[{n, e}, n = Cross[b - a, c - a]; e = sphericalExcess[a, b, c, p]; If[(a - p).n < 0, e, -e]] summedSphericalExcesses::usage := "summedSphericalExcesses[gr,p] gives the \ sum of all signed spherical excesses of all triangles in gr with respect to \ point p." summedSphericalExcesses[gr_, p_] := Module[{polysa, a, b, c, sum}, polys = Cases[gr, Polygon[{a_, b_, c_}], Infinity]; sum = Apply[Plus, polys /. Polygon[{a_, b_, c_}] :> signedSphericalExcess[a, b, c, p]]; sum]; windingNumber3d[gr_, p_] := Round[-summedSphericalExcesses[gr, p]/4./Pi]; Here is an example: (Note that the order of the vertices in the polygons does matter! The vertices have to be anticlock-wise as seen from outside.) v1 = {0, 0, 0}; v2 = {1, 1, 0}; v3 = {1, 0, 1}; v4 = {0, 1, 1}; g = Graphics3D[{FaceForm[RGBColor[1, 0, 0], RGBColor[0, 1, 0]], Polygon[{v1, v2, v3}], Polygon[{v1, v4, v2}], Polygon[{v1, v3, v4}], Polygon[{v2, v4, v3}]}, Lighting -> False]; Show[g]; windingNumber3d[g, {0.5, 0.5, 0.5}] windingNumber3d[g, {0.5, 0.5, -0.5}] The problem for many applications is the requirement of well oriented triangles. (I once programmed a hack for MathWorld to do this but never refined it for efficiency and robustness. Does anyone know about a good implementation in Mathematica?) Anyway, I am happy to know that the winding number solution can be generalized so straightforwardly to three dimensions. :) Cheers Martin