[WIP] Connect holes with zero-width corridors in O(kn).

hgeometry/src/Data/Geometry/HalfLine.hs

@@ -1,4 +1,5 @@
-{-# LANGUAGE TemplateHaskell  #-}
+{-# LANGUAGE DeriveAnyClass       #-}
+{-# LANGUAGE TemplateHaskell      #-}
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE DeriveAnyClass #-}
 --------------------------------------------------------------------------------
@@ -14,7 +15,7 @@ module Data.Geometry.HalfLine where
 import           Control.DeepSeq
 import           Control.Lens
 import           Data.Ext
-import qualified Data.Foldable as F
+import qualified Data.Foldable                as F
 import           Data.Geometry.Interval
 import           Data.Geometry.Line
 import           Data.Geometry.LineSegment
@@ -23,9 +24,9 @@ import           Data.Geometry.Properties
 import           Data.Geometry.SubLine
 import           Data.Geometry.Transformation
 import           Data.Geometry.Vector
-import qualified Data.Traversable as T
+import qualified Data.Traversable             as T
 import           Data.UnBounded
-import           GHC.Generics (Generic)
+import           GHC.Generics                 (Generic)
 import           GHC.TypeLits
 
 --------------------------------------------------------------------------------
@@ -70,6 +71,9 @@ instance (Fractional r, Arity d, Arity (d + 1)) => IsTransformable (HalfLine d r
 
 --------------------------------------------------------------------------------
 
+halfLineToLine :: HalfLine d r -> Line d r
+halfLineToLine (HalfLine anchor dir) = Line anchor dir
+
 halfLineToSubLine                :: (Arity d, Num r)
                                  => HalfLine d r -> SubLine d () (UnBounded r) r
 halfLineToSubLine (HalfLine p v) = let l = Line p v
@@ -100,6 +104,8 @@ type instance IntersectionOf (HalfLine 2 r) (LineSegment 2 p r) = [ NoIntersecti
                                                                   , LineSegment 2 () r
                                                                   ]
 
+-- type instance IntersectionOf  (LineSegment 2 p r) (HalfLine 2 r)
+--   = IntersectionOf (HalfLine 2 r) (LineSegment 2 p r)
 
 -- instance (Ord r, Fractional r) => (HalfLine 2 r) `IsIntersectableWith` (Line 2 r) where
   -- hl `intersect` l = match (halfLineToSubLine hl, l)

hgeometry/src/Data/Geometry/Polygon.hs

@@ -22,6 +22,7 @@ module Data.Geometry.Polygon( PolygonType(..)
 
                             , polygonHoles, polygonHoles'
                             , holeList
+                            , connectHoles
 
                             , inPolygon, insidePolygon, onBoundary
 
@@ -45,18 +46,28 @@ module Data.Geometry.Polygon( PolygonType(..)
                             , extremesLinear, cmpExtreme
                             ) where
 
-import           Algorithms.Geometry.LinearProgramming.LP2DRIC
-import           Algorithms.Geometry.LinearProgramming.Types
-import           Control.Lens hiding (Simple)
-import           Control.Monad.Random.Class
+import           Algorithms.Geometry.LinearProgramming.LP2DRIC (solveBoundedLinearProgram)
+import           Algorithms.Geometry.LinearProgramming.Types   (LinearProgram (LinearProgram))
+import           Control.Lens                                  ((^.))
+import           Control.Monad.Random.Class                    (MonadRandom)
+import           Data.Bifunctor
+import           Data.CircularSeq as C
 import           Data.Ext
-import qualified Data.Foldable as F
-import           Data.Geometry.HalfSpace (rightOf)
+import qualified Data.Foldable                                 as F
+import           Data.Function
+import           Data.Geometry.HalfSpace                       (rightOf)
 import           Data.Geometry.Line
-import           Data.Geometry.Point
+import           Data.Geometry.LineSegment
+import           Data.Geometry.Point                           (Point, vector)
 import           Data.Geometry.Polygon.Core
 import           Data.Geometry.Polygon.Extremes
-
+import           Data.Geometry.Vector                          (Vector,
+                                                                pattern Vector2)
+import           Data.Intersection                             (IsIntersectableWith (intersect))
+import           Data.List                                     (sortBy)
+import           Data.Vinyl                                    (Rec (RNil, (:&)))
+import           Data.Vinyl.CoRec                              (Handler (H),
+                                                                match)
 
 --------------------------------------------------------------------------------
 -- * Polygons
@@ -75,3 +86,61 @@ isStarShaped (toClockwiseOrder -> pg) =
     -- the first vertex is the intersection point of the two supporting lines
     -- bounding it, so the first two edges bound the shape in this sirection
     hs = fmap (rightOf . supportingLine) . outerBoundaryEdges $ pg
+
+-- Properties:
+--   \poly -> selfIntersections poly == selfIntersections (connectHoles poly)
+--   \poly -> area poly == area (connectHoles poly)
+connectHoles :: (Ord r, Num r, Fractional r) =>
+  MultiPolygon p r -> SimplePolygon () r
+connectHoles = connectHolesLinear (Vector2 0 1)
+
+-- | Connect holes with the outer hull by linear cuts.
+--
+-- \(O(kn)\)
+connectHolesLinear :: (Ord r, Num r, Fractional r)
+  => Vector 2 r -> MultiPolygon p r -> SimplePolygon () r
+connectHolesLinear vec (MultiPolygon border holes) =
+    connectHolesLinear' vec (trunc $ SimplePolygon border) sorted
+  where
+    trunc = first (const ())
+    extremes = map (snd . extremesLinearSeq vec) (map trunc holes)
+    sorted = sortBy (cmpExtreme vec `on` focus) extremes
+
+connectHolesLinear' :: (Ord r, Num r, Fractional r)
+  => Vector 2 r -> SimplePolygon () r -> [CSeq (Point 2 r :+ ())]
+  -> SimplePolygon () r
+connectHolesLinear' _ border []      = border
+connectHolesLinear' vec border (hole:holes) =
+    connectHolesLinear' vec newBorder holes
+  where
+    newBorder = SimplePolygon (C.fromList (openBorder ++ F.toList hole ++ [focus hole]))
+    line = Line (focus hole ^. core) vec
+    (_pt, openBorder) =
+      F.minimumBy (cmpExtreme vec `on` fst)
+      [ (pt, outer)
+      | (pt, outer) <- cutPolygon border line
+      , cmpExtreme vec (focus hole) pt == LT ]
+
+-- for each hole, find extreme point
+-- pick hole with largest extreme
+-- cut from that hole to outer poly
+-- iterate until no more holes
+
+cutPolygon :: (Fractional r, Ord r) =>
+  SimplePolygon () r -> Line 2 r -> [(Point 2 r :+ (), [Point 2 r :+ ()])]
+cutPolygon poly line = worker [] (listEdges poly)
+  where
+    fromSegment (LineSegment' firstPt _) = firstPt
+    worker _acc [] = []
+    worker acc (segment:segments) =
+      let acc' = fromSegment segment : acc
+          LineSegment' sStart sEnd = segment
+      in match (segment `intersect` line) $
+             (H $ const $ worker acc' segments)
+          :& (H $ \pt ->
+              let toPolyEnd = map fromSegment segments
+                  fromPolyStart = reverse acc ++ [sStart]
+              in (ext pt, ext pt : sEnd : toPolyEnd  ++ fromPolyStart) :
+                  worker acc' segments)
+          :& (H undefined)
+          :& RNil

hgeometry/src/Data/Geometry/Polygon/Extremes.hs

@@ -8,13 +8,18 @@
 -- Finding the Extremal vertex of a polygon in a given direction.
 --
 --------------------------------------------------------------------------------
-module Data.Geometry.Polygon.Extremes( cmpExtreme
-                                     , extremesLinear
-                                     ) where
+module Data.Geometry.Polygon.Extremes
+    ( cmpExtreme
+    , cmpExtreme'
+    , extremesLinear
+    , extremesLinearSeq
+    ) where
 
-import           Control.Lens hiding (Simple)
+import           Control.Lens               hiding (Simple)
+import           Data.CircularSeq
 import           Data.Ext
-import qualified Data.Foldable as F
+import qualified Data.Foldable              as F
+import           Data.Function
 import           Data.Geometry.Point
 import           Data.Geometry.Polygon.Core
 import           Data.Geometry.Vector
@@ -26,9 +31,13 @@ import           Data.Geometry.Vector
 -- the vector u.
 cmpExtreme       :: (Num r, Ord r)
                  => Vector 2 r -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering
-cmpExtreme u p q = u `dot` (p^.core .-. q^.core) `compare` 0
+cmpExtreme u p q = cmpExtreme' u (p^.core) (q^.core)
 
+cmpExtreme'       :: (Num r, Ord r)
+                 => Vector 2 r -> Point 2 r -> Point 2 r -> Ordering
+cmpExtreme' u p q = u `dot` (p .-. q) `compare` 0
 
+-- FIXME: use extremesLinearSeq in extremesLinear? Will the performance be the same?
 -- | Finds the extreme points, minimum and maximum, in a given direction
 --
 -- running time: \(O(n)\)
@@ -37,3 +46,13 @@ extremesLinear     :: (Ord r, Num r) => Vector 2 r -> Polygon t p r
 extremesLinear u p = let vs = p^.outerBoundary
                          f  = cmpExtreme u
                      in (F.minimumBy f vs, F.maximumBy f vs)
+
+-- | Finds the extreme points, minimum and maximum, in a given direction
+--
+-- running time: \(O(n)\)
+extremesLinearSeq     :: (Ord r, Num r) => Vector 2 r -> Polygon t p r
+                   -> (CSeq (Point 2 r :+ p), CSeq (Point 2 r :+ p))
+extremesLinearSeq u p =
+    let vs = allRotations (p^.outerBoundary)
+        f  = cmpExtreme u `on` focus
+    in (F.minimumBy f vs, F.maximumBy f vs)

hgeometry/test/Data/Geometry/PolygonSpec.hs

@@ -84,6 +84,14 @@ spec = do
     forM_ allMultiPolygons' $ \poly -> do
       hasSelfIntersections poly `shouldBe` False
       isCounterClockwise poly `shouldBe` True
+
+  it "has no self-intersections after connectHoles" $
+    property $ \(poly :: MultiPolygon () Rational) ->
+      not (hasSelfIntersections (connectHoles poly))
+
+  it "has the same area after connectHoles" $
+    property $ \(poly :: MultiPolygon () Rational) ->
+      area poly === area (connectHoles poly)
 
 
 testCases    :: FilePath -> Spec