{-# 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 Data.OrdSeq (OrdSeq)
import qualified Data.OrdSeq as SS
import Data.Semigroup
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)
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
onHole Start = Split
onHole Merge = End
onHole Split = Start
onHole End = Merge
onHole Regular = Regular
classifyVertices' :: (Num r, Ord r)
=> SimplePolygon p r
-> SimplePolygon (p :+ VertexType) r
classifyVertices' (SimplePolygon vs) =
SimplePolygon $ zip3LWith f (rotateL vs) vs (rotateR vs)
where
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
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.OrdSeq (LineSegment 2 Int r))
, _helper :: !(IntMap.IntMap Int)
} deriving (Show)
makeLenses ''StatusStruct
ix' :: Int -> Lens' (V.Vector a) a
ix' i = singular (ix i)
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
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 mempty 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
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
-> OrdSeq (LineSegment 2 q r) -> OrdSeq (LineSegment 2 q r)
insertAt v = SS.insertBy (ordAt $ v^.yCoord)
deleteAt :: (Fractional r, Ord r) => Point 2 r -> LineSegment 2 p r
-> OrdSeq (LineSegment 2 p r) -> OrdSeq (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
tellIfMerge i v iPred
modify $ \ss ->
ss&statusStruct %~ deleteAt v (adj^._1)
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)
getHelper :: Int -> Sweep p r (SP (Point 2 r :+ Int) VertexType)
getHelper i = do Just ui <- gets (^.helper.at i)
STR u _ ut <- asks (^.ix' ui)
pure $ SP (u :+ ui) ut
lookupLE :: (Ord r, Fractional r)
=> Point 2 r -> OrdSeq (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 <> m)
handleSplit :: (Fractional r, Ord r) => Int -> Event r -> Sweep p r ()
handleSplit i (v :+ adj) = do Just ej <- gets $ \ss -> ss^.statusStruct.to (lookupLE v)
let j = ej^.start.extra
SP u _ <- getHelper j
modify $ \(SS t h) ->
SS (insertAt v (adj^._2) t)
(IntMap.insert i i . IntMap.insert j i $ h)
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
tellIfMerge i v ePred
modify $ \ss -> ss&statusStruct %~ deleteAt v (adj^._1)
connectToLeft i v
connectToLeft :: (Fractional r, Ord r) => Int -> Point 2 r -> Sweep p r ()
connectToLeft i v = do Just ej <- gets $ \ss -> ss^.statusStruct.to (lookupLE v)
let j = ej^.start.extra
tellIfMerge i v j
modify $ \ss -> ss&helper %~ IntMap.insert j i
isLeftVertex :: Ord r => Int -> Event r -> Bool
isLeftVertex i (v :+ adj) = case (adj^._1.start) `cmpSweep` (v :+ i) of
GT -> True
_ -> False
handleRegularL :: (Fractional r, Ord r) => Int -> Event r -> Sweep p r ()
handleRegularL i (v :+ adj) = do let ePred = adj^._1.start.extra
tellIfMerge i v ePred
modify $ \ss ->
ss&statusStruct %~ deleteAt v (adj^._1)
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