Re: Re: winding number
- To: mathgroup at smc.vnet.net
- Subject: [mg25436] Re: [mg25405] Re: winding number
- From: Martin Kraus <Martin.Kraus at informatik.uni-stuttgart.de>
- Date: Sun, 1 Oct 2000 02:44:20 -0400 (EDT)
- 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