hgeometry-0.9.0.0: src/Algorithms/Geometry/PolygonTriangulation/MakeMonotone.hs
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Algorithms.Geometry.PolygonTriangulation.MakeMonotone where
import Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann ( xCoordAt
, ordAt)
import Algorithms.Geometry.PolygonTriangulation.Types
import Control.Lens
import Control.Monad (forM_, when)
import Control.Monad.Reader
import Control.Monad.State.Strict
import Control.Monad.Writer (WriterT, execWriterT,tell)
import Data.Bifunctor
import Data.CircularSeq (rotateL, rotateR, zip3LWith)
import qualified Data.DList as DList
import Data.Ext
import qualified Data.Foldable as F
import Data.Geometry.LineSegment
import Data.Geometry.PlanarSubdivision.Basic
import Data.Geometry.Point
import Data.Geometry.Polygon
import qualified Data.IntMap as IntMap
import qualified Data.List.NonEmpty as NonEmpty
import Data.Ord (comparing, Down(..))
import qualified Data.Set as SS
import qualified Data.Set.Util as SS
import Data.Util
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as MV
----------------------------------------------------------------------------------
data VertexType = Start | Merge | Split | End | Regular deriving (Show,Read,Eq)
-- How about the hole vertices?
-- | assigns a vertex type to each vertex
--
-- pre: the polygon is given in CCW order
--
-- running time: \(O(n)\).
classifyVertices :: (Num r, Ord r)
=> Polygon t p r
-> Polygon t (p :+ VertexType) r
classifyVertices p@(SimplePolygon _) = classifyVertices' p
classifyVertices (MultiPolygon vs h) = MultiPolygon vs' h'
where
(SimplePolygon vs') = classifyVertices' $ SimplePolygon vs
h' = map (first (&extra %~ onHole) . classifyVertices') h
-- the roles on hole vertices are slightly different
onHole Start = Split
onHole Merge = End
onHole Split = Start
onHole End = Merge
onHole Regular = Regular
-- | assigns a vertex type to each vertex
--
-- pre: the polygon is given in CCW order
--
-- running time: \(O(n)\).
classifyVertices' :: (Num r, Ord r)
=> SimplePolygon p r
-> SimplePolygon (p :+ VertexType) r
classifyVertices' (SimplePolygon vs) =
SimplePolygon $ zip3LWith f (rotateL vs) vs (rotateR vs)
where
-- is the angle larger than > 180 degrees
largeInteriorAngle p c n = case ccw (p^.core) (c^.core) (n^.core) of
CCW -> False
CW -> True
_ -> error "classifyVertices -> largeInteriorAngle: colinear points"
f p c n = c&extra %~ (:+ vt)
where
vt = case (p `cmpSweep` c, n `cmpSweep` c, largeInteriorAngle p c n) of
(LT, LT, False) -> Start
(LT, LT, True) -> Split
(GT, GT, False) -> End
(GT, GT, True) -> Merge
_ -> Regular
-- | p < q = p.y < q.y || p.y == q.y && p.x > q.y
cmpSweep :: Ord r => Point 2 r :+ e -> Point 2 r :+ e -> Ordering
p `cmpSweep` q =
comparing (^.core.yCoord) p q <> comparing (Down . (^.core.xCoord)) p q
--------------------------------------------------------------------------------
type Event r = Point 2 r :+ (Two (LineSegment 2 Int r))
data StatusStruct r = SS { _statusStruct :: !(SS.Set (LineSegment 2 Int r))
, _helper :: !(IntMap.IntMap Int)
-- ^ for every e_i, the id of the helper vertex
} deriving (Show)
makeLenses ''StatusStruct
ix' :: Int -> Lens' (V.Vector a) a
ix' i = singular (ix i)
-- | Given a polygon, find a set of non-intersecting diagonals that partition
-- the polygon into y-monotone pieces.
--
-- running time: \(O(n\log n)\)
computeDiagonals :: forall t r p. (Fractional r, Ord r)
=> Polygon t p r -> [LineSegment 2 p r]
computeDiagonals p' = map f . sweep
. NonEmpty.sortBy (flip cmpSweep)
. polygonVertices . withIncidentEdges
. first (^._1) $ pg
where
-- remaps to get the p value rather than the vertexId
f = first (\i -> vertexInfo^.ix' i._2)
pg :: Polygon t (SP Int (p :+ VertexType)) r
pg = numberVertices . classifyVertices . toCounterClockWiseOrder $ p'
vertexInfo :: V.Vector (STR (Point 2 r) p VertexType)
vertexInfo = let vs = polygonVertices pg
n = F.length vs
in V.create $ do
v <- MV.new n
forM_ vs $ \(pt :+ SP i (p :+ vt)) ->
MV.write v i (STR pt p vt)
return v
initialSS = SS SS.empty mempty
sweep es = flip runReader vertexInfo $ evalStateT (sweep' es) initialSS
sweep' es = DList.toList <$> execWriterT (sweep'' es)
sweep'' :: NonEmpty.NonEmpty (Event r) -> Sweep p r ()
sweep'' = mapM_ handle
-- | Computes a set of diagionals that decompose the polygon into y-monotone
-- pieces.
--
-- running time: \(O(n\log n)\)
makeMonotone :: (Fractional r, Ord r)
=> proxy s -> Polygon t p r
-> PlanarSubdivision s p PolygonEdgeType PolygonFaceData r
makeMonotone px pg = let (e:es) = listEdges pg
in constructSubdivision px e es (computeDiagonals pg)
type Sweep p r = WriterT (DList.DList (LineSegment 2 Int r))
(StateT (StatusStruct r)
(Reader (V.Vector (VertexInfo p r))))
type VertexInfo p r = STR (Point 2 r) p VertexType
tell' :: LineSegment 2 Int r -> Sweep p r ()
tell' = tell . DList.singleton
getIdx :: Event r -> Int
getIdx = view (extra._1.end.extra)
getVertexType :: Int -> Sweep p r VertexType
getVertexType v = asks (^.ix' v._3)
getEventType :: Event r -> Sweep p r VertexType
getEventType = getVertexType . getIdx
handle :: (Fractional r, Ord r) => Event r -> Sweep p r ()
handle e = let i = getIdx e in getEventType e >>= \case
Start -> handleStart i e
End -> handleEnd i e
Split -> handleSplit i e
Merge -> handleMerge i e
Regular | isLeftVertex i e -> handleRegularL i e
| otherwise -> handleRegularR i e
insertAt :: (Ord r, Fractional r) => Point 2 r -> LineSegment 2 q r
-> SS.Set (LineSegment 2 q r) -> SS.Set (LineSegment 2 q r)
insertAt v = SS.insertBy (ordAt $ v^.yCoord)
deleteAt :: (Fractional r, Ord r) => Point 2 r -> LineSegment 2 p r
-> SS.Set (LineSegment 2 p r) -> SS.Set (LineSegment 2 p r)
deleteAt v = SS.deleteAllBy (ordAt $ v^.yCoord)
handleStart :: (Fractional r, Ord r)
=> Int -> Event r -> Sweep p r ()
handleStart i (v :+ adj) = modify $ \(SS t h) ->
SS (insertAt v (adj^._2) t)
(IntMap.insert i i h)
handleEnd :: (Fractional r, Ord r)
=> Int -> Event r -> Sweep p r ()
handleEnd i (v :+ adj) = do let iPred = adj^._1.start.extra -- i-1
-- lookup p's helper; if it is a merge vertex
-- we insert a new segment
tellIfMerge i v iPred
-- delete e_{i-1} from the status struct
modify $ \ss ->
ss&statusStruct %~ deleteAt v (adj^._1)
-- | Adds edge (i,j) if e_j's helper is a merge vertex
tellIfMerge :: Int -> Point 2 r -> Int -> Sweep p r ()
tellIfMerge i v j = do SP u ut <- getHelper j
when (ut == Merge) (tell' $ ClosedLineSegment (v :+ i) u)
-- | Get the helper of edge i, and its vertex type
getHelper :: Int -> Sweep p r (SP (Point 2 r :+ Int) VertexType)
getHelper i = do ui <- gets (^?!helper.ix i)
STR u _ ut <- asks (^.ix' ui)
pure $ SP (u :+ ui) ut
lookupLE :: (Ord r, Fractional r)
=> Point 2 r -> SS.Set (LineSegment 2 Int r)
-> Maybe (LineSegment 2 Int r)
lookupLE v s = let (l,m,_) = SS.splitOn (xCoordAt $ v^.yCoord) (v^.xCoord) s
in SS.lookupMax (l `SS.join` m)
handleSplit :: (Fractional r, Ord r) => Int -> Event r -> Sweep p r ()
handleSplit i (v :+ adj) = do ej <- gets $ \ss -> ss^?!statusStruct.to (lookupLE v)._Just
let j = ej^.start.extra
SP u _ <- getHelper j
-- update the status struct:
-- insert the new edge into the status Struct and
-- set the helper of e_j to be v_i
modify $ \(SS t h) ->
SS (insertAt v (adj^._2) t)
(IntMap.insert i i . IntMap.insert j i $ h)
-- return the diagonal
tell' $ ClosedLineSegment (v :+ i) u
handleMerge :: (Fractional r, Ord r) => Int -> Event r -> Sweep p r ()
handleMerge i (v :+ adj) = do let ePred = adj^._1.start.extra -- i-1
tellIfMerge i v ePred
-- delete e_{i-1} from the status struct
modify $ \ss -> ss&statusStruct %~ deleteAt v (adj^._1)
connectToLeft i v
-- | finds the edge j to the left of v_i, and connect v_i to it if the helper
-- of j is a merge vertex
connectToLeft :: (Fractional r, Ord r) => Int -> Point 2 r -> Sweep p r ()
connectToLeft i v = do ej <- gets $ \ss -> ss^?!statusStruct.to (lookupLE v)._Just
let j = ej^.start.extra
tellIfMerge i v j
modify $ \ss -> ss&helper %~ IntMap.insert j i
-- | returns True if v the interior of the polygon is to the right of v
isLeftVertex :: Ord r => Int -> Event r -> Bool
isLeftVertex i (v :+ adj) = case (adj^._1.start) `cmpSweep` (v :+ i) of
GT -> True
_ -> False
-- if the predecessor occurs before the sweep, this must be a left vertex
handleRegularL :: (Fractional r, Ord r) => Int -> Event r -> Sweep p r ()
handleRegularL i (v :+ adj) = do let ePred = adj^._1.start.extra -- i-1
tellIfMerge i v ePred
-- delete e_{i-1} from the status struct
modify $ \ss ->
ss&statusStruct %~ deleteAt v (adj^._1)
-- insert a e_i in the status struct, and set its helper
-- to be v_i
modify $ \(SS t h) ->
SS (insertAt v (adj^._2) t)
(IntMap.insert i i h)
handleRegularR :: (Fractional r, Ord r) => Int -> Event r -> Sweep p r ()
handleRegularR i (v :+ _) = connectToLeft i v
--------------------------------------------------------------------------------
-- testPolygon :: SimplePolygon Int Rational
-- testPolygon = fromPoints [ point2 20 20 :+ 1
-- , point2 18 19 :+ 2
-- , point2 16 25 :+ 3
-- , point2 13 23 :+ 4
-- , point2 10 24 :+ 5
-- , point2 6 22 :+ 6
-- , point2 8 21 :+ 7
-- , point2 7 18 :+ 8
-- , point2 2 19 :+ 9
-- , point2 1 10 :+ 10
-- , point2 3 5 :+ 11
-- , point2 11 7 :+ 12
-- , point2 15 1 :+ 13
-- , point2 12 15 :+ 14
-- , point2 15 12 :+ 15
-- ]
-- vertexTypes = [Start,Merge,Start,Merge,Start,Regular,Regular,Merge,Start,Regular,End,Split,End,Split,End]
-- loadT = do pgs <- readAllFrom "/Users/frank/tmp/testPoly.ipe"
-- :: IO [SimplePolygon () Rational :+ IpeAttributes Path Rational]
-- mapM_ print pgs
-- let diags = map (computeDiagonals . (^.core)) pgs
-- f = asIpeGroup . map (asIpeObject' mempty)
-- out = [ asIpeGroup $ map (\(pg :+ a) -> asIpeObject pg a) pgs
-- , asIpeGroup $ map f diags
-- ]
-- outFile = "/Users/frank/tmp/out.ipe"
-- writeIpeFile outFile . singlePageFromContent $ out