Add function for cutting polygons

src/Geometry/Polygon.hs

@@ -202,3 +202,41 @@ shrinkPolyOnEdges _ _ = error "too few vertices in polygon"
 shrinkVert :: Float -> [Point2] -> Point2
 shrinkVert d [x, y, z] = x +.+ (d *.* normalizeV (x -.- y)) +.+ (d *.* normalizeV (z -.- y))
 shrinkVert _ _ = error "wrong number of vertices"
+
+-- divide a polygon into two along a line
+-- assumes the line intersects the polygon exactly twice, nocolinearity
+-- this may duplicate points
+-- this should be tested: there are many possible points of failure...
+cutPoly :: Point2 -> Point2 -> [Point2] -> ([Point2],[Point2])
+cutPoly a b (p:ps)
+    | not (isLHS a b p) && not (isRHS a b p) = cutPolyL a b ps p p [p] [p]
+    | isLHS a b p = cutPolyL a b ps p p [p] []
+    | otherwise = cutPolyR a b ps p p [] [p]
+cutPoly _ _ [] = error "cutPoly empty poly"
+
+cutPolyL :: Point2 -> Point2 -> [Point2] -> Point2 -> Point2 -> [Point2] -> [Point2] -> ([Point2], [Point2])
+cutPolyL a b (p:ps) e x ls rs
+    | isLHS a b p = cutPolyL a b ps e p (p\:ls) rs
+    | otherwise = case intersectLineLine x p a b of
+        Nothing -> error "cutPolyL nonintersecting lines"
+        Just p' -> cutPolyR a b ps e p (p'\:ls) (p\:p'\:rs)
+cutPolyL a b [] e x ls rs = case intersectSegLine x e a b of
+    Nothing -> (reverse ls,reverse rs)
+    Just p -> (reverse (p\:ls), reverse (p\:rs))
+
+cutPolyR :: Point2 -> Point2 -> [Point2] -> Point2 -> Point2 -> [Point2] -> [Point2] -> ([Point2], [Point2])
+cutPolyR a b (p:ps) e x ls rs
+    | isRHS a b p = cutPolyR a b ps e p ls (p\:rs)
+    | otherwise = case intersectLineLine x p a b of
+        Nothing -> error "cutPolyR nonintersecting lines"
+        Just p' -> cutPolyL a b ps e p (p\:p'\:ls) (p'\:rs)
+cutPolyR a b [] e x ls rs = case intersectSegLine x e a b of
+    Nothing -> (reverse ls,reverse rs)
+    Just p -> (reverse (p\:ls), reverse (p\:rs))
+
+infixr 5 \:
+(\:) :: Eq a => a -> [a] -> [a]
+(\:) x (y:ys) 
+    | x /= y = x:y:ys
+    | otherwise = y:ys 
+(\:) x [] = [x]