Sunday, October 12, 2008

Points inside a polygon and the Residue Theorem

A usual way to determine whether a 2D point p lies inside a simple polygon P consists in tracing a ray from p to infinity and computing how many times this ray intersects P: p is inside P iff the number of intersections is odd. We investigate a different technique for solving the problem based on complex analysis, only to finally see that this new method is really equivalent to the classical ray algorithm.

We consider p and P to be represented in the complex plane and assume, without loss of generality, that p = 0. The Residue Theorem tells us that

P (1/z) dz

is null if the point 0 (the only pole of 1/z) is outside P, and

2πi Resz=0(1/z) = 2πi

if 0 lies inside P. So we need only evaluate the integral P (1/z) dz; as P is a loop sequence of segments going counterclockwise like this

P = (p1,p2), (p2,p3), ... , (pN−1,pN), (pN,p1),

and ln z is an antiderivative of 1/z, we have that

P (1/z) dz = (pi,pi+1) (1/z) dz =
= (ln pi+1 − ln pi).

But now, this sum seems to be unconditionally null:

(ln pi+1 − ln pi) =
= ln p2 − ln p1 + ln p3 − ln p2 + ln p4 − ln p3 + ... + ln p1 − ln pN = 0.

Each term is cancelled out by some other with different sign, so the total sum is, or seems to be, zero! Obviously we have made some mistake during the process, and in fact the wrong assumption is this: the equality

(pi,pi+1) (1/z) dz = ln pi+1 − ln pi

only holds if the derivative of ln z is 1/z along all the segment (pi,pi+1) (more precisely, at a region containing the segment): but ln z, or, strictly speaking, the principal branch of ln z, has a cut at the negative real axis, so the equality above fails to hold when (pi,pi+1) intersects this cut. When this is the case, we can evaluate the integral by using a different branch of the complex logarithm so that (pi,pi+1) does not intersect the new cut. It is not hard to see that doing so yields

(pi,pi+1) (1/z) dz = ln pi+1 − ln pi + σ2πi,

where ln denotes again the principal branch of the complex logarithm and σ is +1 or −1 depending on whether (pi,pi+1) crosses the negative real axis southwards or northwards. Having the situation with the negative real axis covered, we can easily deduce that

P (1/z) dz =(sn)2πi,

where s is the number of segments of P crossing the negative real axis southwards, and n corresponds to those crossing northwards. And, in the end, this is but a simple reformulation of the old ray-based algorithm.

No comments:

Post a Comment