--------------------------------------------------------------------------------
-- |
-- Module      :  Data.Geometry.Polygon.Monotone
-- Copyright   :  (C) 1ndy
-- License     :  see the LICENSE file
-- Maintainer  :  David Himmelstrup
--
-- A polygon is monotone in a certain direction if rays orthogonal to that
-- direction intersects the polygon at most twice. See
-- <https://en.wikipedia.org/wiki/Monotone_polygon>
--
--------------------------------------------------------------------------------
module Data.Geometry.Polygon.Monotone
  ( isMonotone
  , randomMonotone
  , randomMonotoneDirected
  , monotoneFrom
  , randomNonZeroVector
  ) where

import           Control.Monad.Random
import           Data.Ext
import qualified Data.Foldable                  as F
import           Data.Geometry.Line             (Line (..))
import           Data.Geometry.LineSegment
import           Data.Geometry.Point
import           Data.Geometry.Polygon.Core
import           Data.Geometry.Polygon.Extremes
import           Data.Geometry.Vector
import           Data.Intersection
import           Data.List
import           Data.Vinyl
import           Data.Vinyl.CoRec
import           Prelude                        hiding (max, min)

-- | \( O(n \log n) \)
--   A polygon is monotone if a straight line in a given direction
--   cannot have more than two intersections.
isMonotone :: (Fractional r, Ord r) => Vector 2 r -> SimplePolygon p r -> Bool
-- Check for each vertex that the number of intersections with the
-- line starting at the vertex and going out in the given direction
-- intersects with the edges of the polygon no more than 2 times.
isMonotone :: Vector 2 r -> SimplePolygon p r -> Bool
isMonotone Vector 2 r
direction SimplePolygon p r
p = (Point 2 r -> Bool) -> [Point 2 r] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Point 2 r -> Bool
isMonotoneAt (((Point 2 r :+ p) -> Point 2 r) -> [Point 2 r :+ p] -> [Point 2 r]
forall a b. (a -> b) -> [a] -> [b]
map (Point 2 r :+ p) -> Point 2 r
forall core extra. (core :+ extra) -> core
_core ([Point 2 r :+ p] -> [Point 2 r])
-> [Point 2 r :+ p] -> [Point 2 r]
forall a b. (a -> b) -> a -> b
$ SimplePolygon p r -> [Point 2 r :+ p]
forall (t :: PolygonType) p r. Polygon t p r -> [Point 2 r :+ p]
toPoints SimplePolygon p r
p)
  where
    isMonotoneAt :: Point 2 r -> Bool
isMonotoneAt Point 2 r
pt =
      [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((LineSegment 2 p r -> Integer) -> [LineSegment 2 p r] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Point 2 r -> LineSegment 2 p r -> Integer
intersectionsThrough Point 2 r
pt) (CircularVector (LineSegment 2 p r) -> [LineSegment 2 p r]
forall (t :: * -> *) a. Foldable t => t a -> [a]
F.toList (CircularVector (LineSegment 2 p r) -> [LineSegment 2 p r])
-> CircularVector (LineSegment 2 p r) -> [LineSegment 2 p r]
forall a b. (a -> b) -> a -> b
$ SimplePolygon p r -> CircularVector (LineSegment 2 p r)
forall (t :: PolygonType) p r.
Polygon t p r -> CircularVector (LineSegment 2 p r)
outerBoundaryEdges SimplePolygon p r
p)) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
2
    intersectionsThrough :: Point 2 r -> LineSegment 2 p r -> Integer
intersectionsThrough Point 2 r
pt LineSegment 2 p r
edge =
      CoRec Identity '[NoIntersection, Point 2 r, LineSegment 2 p r]
-> Handlers '[NoIntersection, Point 2 r, LineSegment 2 p r] Integer
-> Integer
forall (ts :: [*]) b. CoRec Identity ts -> Handlers ts b -> b
match (LineSegment 2 p r
-> Line 2 r -> Intersection (LineSegment 2 p r) (Line 2 r)
forall g h. IsIntersectableWith g h => g -> h -> Intersection g h
Data.Intersection.intersect LineSegment 2 p r
edge Line 2 r
line) (Handlers '[NoIntersection, Point 2 r, LineSegment 2 p r] Integer
 -> Integer)
-> Handlers '[NoIntersection, Point 2 r, LineSegment 2 p r] Integer
-> Integer
forall a b. (a -> b) -> a -> b
$
           (NoIntersection -> Integer) -> Handler Integer NoIntersection
forall b a. (a -> b) -> Handler b a
H (\NoIntersection
NoIntersection -> Integer
0)
        Handler Integer NoIntersection
-> Rec (Handler Integer) '[Point 2 r, LineSegment 2 p r]
-> Handlers '[NoIntersection, Point 2 r, LineSegment 2 p r] Integer
forall u (a :: u -> *) (r :: u) (rs :: [u]).
a r -> Rec a rs -> Rec a (r : rs)
:& (Point 2 r -> Integer) -> Handler Integer (Point 2 r)
forall b a. (a -> b) -> Handler b a
H (\Point{} -> Integer
1)
        -- This happens when an edge is parallel with the given direction.
        -- I think it's correct to count it as a single intersection.
        Handler Integer (Point 2 r)
-> Rec (Handler Integer) '[LineSegment 2 p r]
-> Rec (Handler Integer) '[Point 2 r, LineSegment 2 p r]
forall u (a :: u -> *) (r :: u) (rs :: [u]).
a r -> Rec a rs -> Rec a (r : rs)
:& (LineSegment 2 p r -> Integer)
-> Handler Integer (LineSegment 2 p r)
forall b a. (a -> b) -> Handler b a
H (\LineSegment{} -> Integer
1)
        Handler Integer (LineSegment 2 p r)
-> Rec (Handler Integer) '[]
-> Rec (Handler Integer) '[LineSegment 2 p r]
forall u (a :: u -> *) (r :: u) (rs :: [u]).
a r -> Rec a rs -> Rec a (r : rs)
:& Rec (Handler Integer) '[]
forall u (a :: u -> *). Rec a '[]
RNil
      where
        line :: Line 2 r
line = Point 2 r -> Vector 2 r -> Line 2 r
forall (d :: Nat) r. Point d r -> Vector d r -> Line d r
Line Point 2 r
pt (Vector 2 r -> Vector 2 r
forall r. Num r => Vector 2 r -> Vector 2 r
rot90 Vector 2 r
direction)
        rot90 :: Vector 2 r -> Vector 2 r
rot90 (Vector2 r
x r
y) = r -> r -> Vector 2 r
forall r. r -> r -> Vector 2 r
Vector2 (-r
y) r
x

{- Algorithm overview:

  1. Create N `Point 2 Rational` (N >= 3)
  2. Create a random `Vector 2 Rational`
  3. Find the extremes (min and max) of the points when sorted in the direction of the vector.
      We already have code for this. See `maximumBy (cmpExtreme vector)` and
      `minimumBy (cmpExtreme vector)`.
  4. Take out the two extremal points from the set.
  5. Partition the remaining points according to whether they're on the left side or right side
    of the imaginary line between the two extremal points.
  6. Sort the two partitioned sets, one in the direction of the vector and one in the opposite
    direction.
  7. Connect the points, starting from the minimal extreme point, going through the set of points
    that are increasing in the direction of the vector, then to the maximal point, and finally
    down through the points that are decreasing in the direction of the vector.
-}
-- | \( O(n \log n) \)
--   Generate a random N-sided polygon that is monotone in a random direction.
randomMonotone :: (RandomGen g, Random r, Ord r, Num r) => Int -> Rand g (SimplePolygon () r)
randomMonotone :: Int -> Rand g (SimplePolygon () r)
randomMonotone Int
nVertices = Int -> Vector 2 r -> Rand g (SimplePolygon () r)
forall g r.
(RandomGen g, Random r, Ord r, Num r) =>
Int -> Vector 2 r -> Rand g (SimplePolygon () r)
randomMonotoneDirected Int
nVertices (Vector 2 r -> Rand g (SimplePolygon () r))
-> RandT g Identity (Vector 2 r) -> Rand g (SimplePolygon () r)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< RandT g Identity (Vector 2 r)
forall g r.
(RandomGen g, Random r, Eq r, Num r) =>
Rand g (Vector 2 r)
randomNonZeroVector

-- Pick a random vector and then call 'randomMonotone'.
-- | \( O(n \log n) \)
--   Generate a random N-sided polygon that is monotone in the given direction.
randomMonotoneDirected :: (RandomGen g, Random r, Ord r, Num r)
  => Int -> Vector 2 r -> Rand g (SimplePolygon () r)
randomMonotoneDirected :: Int -> Vector 2 r -> Rand g (SimplePolygon () r)
randomMonotoneDirected Int
nVertices Vector 2 r
direction = do
    [Point 2 r]
points <- Int -> RandT g Identity (Point 2 r) -> RandT g Identity [Point 2 r]
forall (m :: * -> *) a. Applicative m => Int -> m a -> m [a]
replicateM Int
nVertices RandT g Identity (Point 2 r)
forall (m :: * -> *) a. (MonadRandom m, Random a) => m a
getRandom
    SimplePolygon () r -> Rand g (SimplePolygon () r)
forall (m :: * -> *) a. Monad m => a -> m a
return (Vector 2 r -> [Point 2 r] -> SimplePolygon () r
forall r.
(Ord r, Num r) =>
Vector 2 r -> [Point 2 r] -> SimplePolygon () r
monotoneFrom Vector 2 r
direction [Point 2 r]
points)

-- | \( O(n \log n) \)
--   Assemble a given set of points in a polygon that is monotone in the given direction.
monotoneFrom :: (Ord r, Num r) => Vector 2 r -> [Point 2 r] -> SimplePolygon () r
monotoneFrom :: Vector 2 r -> [Point 2 r] -> SimplePolygon () r
monotoneFrom Vector 2 r
direction [Point 2 r]
vertices = [Point 2 r :+ ()] -> SimplePolygon () r
forall p r. (Eq r, Num r) => [Point 2 r :+ p] -> SimplePolygon p r
fromPoints ([Point 2 r :+ ()
min] [Point 2 r :+ ()] -> [Point 2 r :+ ()] -> [Point 2 r :+ ()]
forall a. [a] -> [a] -> [a]
++ [Point 2 r :+ ()]
rightHalf [Point 2 r :+ ()] -> [Point 2 r :+ ()] -> [Point 2 r :+ ()]
forall a. [a] -> [a] -> [a]
++ [Point 2 r :+ ()
max] [Point 2 r :+ ()] -> [Point 2 r :+ ()] -> [Point 2 r :+ ()]
forall a. [a] -> [a] -> [a]
++ [Point 2 r :+ ()]
leftHalf)
    where
        specialPoints :: [Point 2 r :+ ()]
specialPoints = (Point 2 r -> Point 2 r :+ ()) -> [Point 2 r] -> [Point 2 r :+ ()]
forall a b. (a -> b) -> [a] -> [b]
map (\Point 2 r
x -> Point 2 r
x Point 2 r -> () -> Point 2 r :+ ()
forall core extra. core -> extra -> core :+ extra
:+ ()) [Point 2 r]
vertices
        min :: Point 2 r :+ ()
min = ((Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Ordering)
-> [Point 2 r :+ ()] -> Point 2 r :+ ()
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
Data.List.minimumBy (Vector 2 r -> (Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Ordering
forall r p q.
(Num r, Ord r) =>
Vector 2 r -> (Point 2 r :+ p) -> (Point 2 r :+ q) -> Ordering
cmpExtreme Vector 2 r
direction) [Point 2 r :+ ()]
specialPoints
        max :: Point 2 r :+ ()
max = ((Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Ordering)
-> [Point 2 r :+ ()] -> Point 2 r :+ ()
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
Data.List.maximumBy (Vector 2 r -> (Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Ordering
forall r p q.
(Num r, Ord r) =>
Vector 2 r -> (Point 2 r :+ p) -> (Point 2 r :+ q) -> Ordering
cmpExtreme Vector 2 r
direction) [Point 2 r :+ ()]
specialPoints
        -- 4
        pointsWithoutExtremes :: [Point 2 r :+ ()]
pointsWithoutExtremes = ((Point 2 r :+ ()) -> Bool)
-> [Point 2 r :+ ()] -> [Point 2 r :+ ()]
forall a. (a -> Bool) -> [a] -> [a]
filter (\Point 2 r :+ ()
x -> Point 2 r :+ ()
x (Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Bool
forall a. Eq a => a -> a -> Bool
/= Point 2 r :+ ()
min Bool -> Bool -> Bool
&& Point 2 r :+ ()
x (Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Bool
forall a. Eq a => a -> a -> Bool
/= Point 2 r :+ ()
max) [Point 2 r :+ ()]
specialPoints
        -- 5, 6
        ([Point 2 r :+ ()]
leftHalfUnsorted,[Point 2 r :+ ()]
rightHalfUnsorted) = ((Point 2 r :+ ()) -> Bool)
-> [Point 2 r :+ ()] -> ([Point 2 r :+ ()], [Point 2 r :+ ()])
forall a. (a -> Bool) -> [a] -> ([a], [a])
Data.List.partition ((Point 2 r :+ ()) -> (Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Bool
forall r.
(Ord r, Num r) =>
(Point 2 r :+ ()) -> (Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Bool
toTheLeft Point 2 r :+ ()
min Point 2 r :+ ()
max) [Point 2 r :+ ()]
pointsWithoutExtremes
        leftHalf :: [Point 2 r :+ ()]
leftHalf = ((Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Ordering)
-> [Point 2 r :+ ()] -> [Point 2 r :+ ()]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy (((Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Ordering)
-> (Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Ordering
forall a b c. (a -> b -> c) -> b -> a -> c
flip (((Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Ordering)
 -> (Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Ordering)
-> ((Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Ordering)
-> (Point 2 r :+ ())
-> (Point 2 r :+ ())
-> Ordering
forall a b. (a -> b) -> a -> b
$ Vector 2 r -> (Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Ordering
forall r p q.
(Num r, Ord r) =>
Vector 2 r -> (Point 2 r :+ p) -> (Point 2 r :+ q) -> Ordering
cmpExtreme Vector 2 r
direction) [Point 2 r :+ ()]
leftHalfUnsorted
        rightHalf :: [Point 2 r :+ ()]
rightHalf = ((Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Ordering)
-> [Point 2 r :+ ()] -> [Point 2 r :+ ()]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy (Vector 2 r -> (Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Ordering
forall r p q.
(Num r, Ord r) =>
Vector 2 r -> (Point 2 r :+ p) -> (Point 2 r :+ q) -> Ordering
cmpExtreme Vector 2 r
direction) [Point 2 r :+ ()]
rightHalfUnsorted

-------------------------------------------------------------------------------------------------
-- helper functions

-- for partitioning points
toTheLeft :: (Ord r, Num r) => Point 2 r :+ () -> Point 2 r :+ () -> Point 2 r :+ () -> Bool
toTheLeft :: (Point 2 r :+ ()) -> (Point 2 r :+ ()) -> (Point 2 r :+ ()) -> Bool
toTheLeft Point 2 r :+ ()
min Point 2 r :+ ()
max Point 2 r :+ ()
x = (Point 2 r :+ ()) -> (Point 2 r :+ ()) -> (Point 2 r :+ ()) -> CCW
forall r a b c.
(Ord r, Num r) =>
(Point 2 r :+ a) -> (Point 2 r :+ b) -> (Point 2 r :+ c) -> CCW
ccw' Point 2 r :+ ()
min Point 2 r :+ ()
max Point 2 r :+ ()
x CCW -> CCW -> Bool
forall a. Eq a => a -> a -> Bool
== CCW
CCW

-- | \( O(1) \)
--   Create a random 2D vector which has a non-zero magnitude.
randomNonZeroVector :: (RandomGen g, Random r, Eq r, Num r) => Rand g (Vector 2 r)
randomNonZeroVector :: Rand g (Vector 2 r)
randomNonZeroVector = do
    Vector 2 r
v <- Rand g (Vector 2 r)
forall (m :: * -> *) a. (MonadRandom m, Random a) => m a
getRandom
    if (Vector 2 r -> r
forall (f :: * -> *) a. (Metric f, Num a) => f a -> a
quadrance Vector 2 r
vr -> r -> Bool
forall a. Eq a => a -> a -> Bool
==r
0)
      then Rand g (Vector 2 r)
forall g r.
(RandomGen g, Random r, Eq r, Num r) =>
Rand g (Vector 2 r)
randomNonZeroVector
      else Vector 2 r -> Rand g (Vector 2 r)
forall (f :: * -> *) a. Applicative f => a -> f a
pure Vector 2 r
v