MathGroup Archive 2000

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: winding number

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
> 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 -

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,
    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]];

windingNumber3d[gr_, p_] := Round[-summedSphericalExcesses[gr,

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,
        Polygon[{v2, v4, v3}]}, Lighting -> False];


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. :)



  • Prev by Date: Re: Rounding Numbers in for output in GridBox
  • Next by Date: Bug: Mathematica multipies statements together
  • Previous by thread: winding number
  • Next by thread: Differential operators, Help