There are a couple of ways to do this (the basic problem being how to deal with concave polygons). One way is to convert the polygon to a convex hull and a list of triangles representing subtractions from the convex hull, then check containment within the hull and non-containment within the triangles. However, I present a simpler way: pick a point guaranteed not to be in the polygon, form the ray between that and the point in question, and count the number of times that ray crosses the polygon sides. If it’s even, the point is outside the polygon. This method also works for self-crossing polygons (provided that one considers any internal spaces formed as outside the polygon: consider a pentagram for example).
raysCross :: Ray -> Ray -> Bool
raysCross (a,b) (x,y) = let h = isLeftOf x (a,b) && isRightOf y (a,b)
i = isLeftOf y (a,b) && isRightOf x (a,b)
j = isLeftOf a (x,y) && isRightOf b (x,y)
k = isLeftOf b (x,y) && isRightOf a (x,y)
in (h || i) && (j || k)
countCrossings :: Ray -> [Ray] -> Int
countCrossings r rs = foldl (+) 0
(map (\x -> if raysCross r x then 1 else 0) rs)
guaranteedOutside :: [Vertex] -> Vertex
guaranteedOutside vs
= (foldl maxop 0 xs + 1, foldl maxop 0 ys + 1)
where (xs, ys) = unzip vs
maxop a b = if a > b then a else b
containsS' :: Shape -> Coordinate -> Bool
(Polygon ps) `containsS'` p
= let poutside = guaranteedOutside ps
in (countCrossings (p,poutside)
(zip ps (tail ps ++ [head ps]))) `mod` 2 == 1
See Exercise 8.4 for definitions of isLeftOf and isRightOf. I am also left feeling that there must be a more elegant way to count the number of list elements that fulfil a predicate (which is what countCrossings is really doing).
3 responses to “Exercise 8.12”
Yes, there is a more elegant way to count the number of list elements that fulfil a predicate. I used:
crossings = length (filter (raysCross ray) sides)
There’s a more elegant method than maxop, which is to use the Standard Prelude maximum function.
Unfortunately, however, the method shown does not work in this case:
let poly = Polygon ([(-1,-1),(1,-1),(1,1),(-1,1)]) p = (0,0) in poly `containsS'` pgives a result of False, even though the point (0,0) clearly lies inside the square.
The reason is that the point outside that is chosen is (2,2), and the ray (0,0) to (2,2) passes through the vertex (1,1), and so lies on two edges.
This problem is therefore far subtler than it first appears.
One approach is to use the given idea but to take account of the fixed ray (from chosen point to outside point) passing through vertices, by having each one counting for 0.5. However, this will also fail if the ray just touches a vertex (if the shape is concave). A fix for this is to use signed crossings, taking account of the direction in which the rays cross. But this is getting fairly awkward and difficult.
I’m currently working on an alternative approach which uses winding numbers – we’ll see how it works out ….
Julian
Here goes ….
We include a little test program at the end.
{- A version of containsS for polygons which handles non-convex polygons. We compute the winding number of the polygon around the point as follows. For the point p and consecutive vertices a and b of the polygon (in that order), we have to turn from the vector PA = a-p to the vector PB = b-p. We wish to know the angle between these vectors, measured in an anticlockwise direction, and such the absolute angle is 0 <= angle <= pi. If the angle is pi or either PA or PB is zero, then P lies on the line segment AB, possibly at an endpoint. In this case, the point is definitely contained within the polygon, and we need go no further. Otherwise, we proceed as follows. We consider the plane as the complex plane, and write PA = a-p = z1 = x1 + i*y1, PB = b-p = z2 = x2 + i*y2. Then the rotation+enlargement from PA to PB is given by z2/z1 = (x2 + i*y2) / (x1 + i*y1) = ( (x2 + i*y2) * (x1 - i*y1) ) / (x1^2 + y1^2) = ( (x1*x2 + y1*y2) + i*(x1*y2 - x2*y1) ) / (x1^2 + y1^2) Since we are only interested in the angle and not the enlargement, we can find this as the angle of the vector (x1*x2 + y1*y2, x1*y2 - x2*y1) from the origin. This also save us from having to perform complicated calculations to determine whether either of PA or PB is zero; PA or PB is zero if both the real part (x1*x2 + y1*y2) and the imaginary part (x1*y2 - x2*y1) are zero, and the angle is pi if the imaginary part is zero and the real part is negative. We thus compute this for each pair of adjacent vertices, and return a pair (hit, angle), with hit==True if the point p lies on the side AB and False otherwise, and angle is the required angle. We use the Standard Prelude function atan2 y x which returns the angle from the origin to the point (x,y), measured anticlockwise from the x-axis and lying in the range -pi <= atan2 y x <= pi. Note the order of the arguments of atan2!! Finally, if any hit is True, the point is contained within the polygon. Otherwise, we add up all the angles, divide the answer by 2*pi, and round to the nearest integer (to overcome rounding errors). If the final answer is 0, then the point is outside the polygon, otherwise it is inside. (If the final answer is 2, it means that the polygon winds twice around this point, and the polygon is self-intersecting.) -} winding :: Coordinate -> Ray -> (Bool,Float) winding (px,py) ((ax,ay),(bx,by)) = let x1 = ax - px y1 = ay - py x2 = bx - px y2 = by - py re = x1*x2 + y1*y2 im = x1*y2 - x2*y1 in aux re im where aux r i | r <= 0 && i == 0 = (True,0) | otherwise = (False, atan2 i r) containsS :: Shape -> Coordinate -> Bool (Rectangle s1 s2) `containsS` (x,y) = let t1 = s1/2; t2 = s2/2 in (-t1<=x) && (x<=t1) && (-t2<=y) && (y<=t2) (Ellipse r1 r2) `containsS` (x,y) = (x/r1)^2 + (y/r2)^2 <= 1 (Polygon pts) `containsS` p = let windings = map (winding p) (zip pts (tail pts ++ [head pts])) (hits, angles) = unzip windings in or hits || round (sum angles / (2 * pi)) /= 0 (RtTriangle s1 s2) `containsS` p = (Polygon [(0,0),(s1,0),(0,s2)]) `containsS` p shs :: [Shape] shs = [ Polygon [(-1,-1),(-1,1),(1,1),(1,-1)], Polygon [(-1,-1),(1,-1),(1,1),(-1,1)], Polygon [(-1,-1),(0,-1),(-1,1),(0,1)], Polygon [(-1,-1),(-1,1),(1,1),(1,-1),(-1,-1)], Polygon [(-1,-1),(-1,1),(1,-1)], Polygon [(-1,-1),(-1,1),(0,0)], Polygon [(-1,-1),(-1,1),(0,1)] ] shcontains = map (\x -> x `containsS` (0,0)) shs main :: IO() main = do putStrLn (show shcontains)