--------------------------------------------------------------------------------
-- |
-- Module      :  Algorithms.Geometry.InPolygon
-- Copyright   :  (C) Frank Staals
-- License     :  see the LICENSE file
-- Maintainer  :  Frank Staals
--
-- Testing if a point lies in a polygon
--
--------------------------------------------------------------------------------
module Algorithms.Geometry.InPolygon
  ( inPolygon
  , insidePolygon
  , onBoundary
  ) where

import           Control.Lens

import           Data.Ext
import qualified Data.Foldable as F
import           Data.Geometry.Boundary
import           Data.Geometry.Line
import           Data.Geometry.LineSegment
import           Data.Geometry.Point
import           Data.Geometry.Polygon.Core
import           Data.Geometry.Properties

import qualified Data.List.Util as List
import           Data.Maybe (mapMaybe)
import           Data.Vinyl.CoRec (asA)
--------------------------------------------------------------------------------

{- $setup
>>> import Data.RealNumber.Rational
>>> import Data.Foldable
>>> import Control.Lens.Extras
>>> :{
-- import qualified Data.Vector.Circular as CV
let simplePoly :: SimplePolygon () (RealNumber 10)
    simplePoly = fromPoints . map ext $
      [ Point2 0 0
      , Point2 10 0
      , Point2 10 10
      , Point2 5 15
      , Point2 1 11
      ]
    simpleTriangle :: SimplePolygon () (RealNumber 10)
    simpleTriangle = fromPoints  . map ext $
      [ Point2 0 0, Point2 2 0, Point2 1 1]
    multiPoly :: MultiPolygon () (RealNumber 10)
    multiPoly = MultiPolygon
      (fromPoints . map ext $ [Point2 (-1) (-1), Point2 3 (-1), Point2 2 2])
      [simpleTriangle]
:} -}


-- | \( O(n) \) Test if q lies on the boundary of the polygon.
--
-- >>> Point2 1 1 `onBoundary` simplePoly
-- False
-- >>> Point2 0 0 `onBoundary` simplePoly
-- True
-- >>> Point2 10 0 `onBoundary` simplePoly
-- True
-- >>> Point2 5 13 `onBoundary` simplePoly
-- False
-- >>> Point2 5 10 `onBoundary` simplePoly
-- False
-- >>> Point2 10 5 `onBoundary` simplePoly
-- True
-- >>> Point2 20 5 `onBoundary` simplePoly
-- False
--
-- TODO: testcases multipolygon
onBoundary        :: (Num r, Ord r) => Point 2 r -> Polygon t p r -> Bool
Point 2 r
q onBoundary :: Point 2 r -> Polygon t p r -> Bool
`onBoundary` Polygon t p r
pg = (LineSegment 2 p r -> Bool) -> [LineSegment 2 p r] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Point 2 r
q Point 2 r -> LineSegment 2 p r -> Bool
forall g h. IsIntersectableWith g h => g -> h -> Bool
`intersects`) [LineSegment 2 p r]
es
  where
    out :: SimplePolygon p r
out = Polygon t p r
pgPolygon t p r
-> Getting (SimplePolygon p r) (Polygon t p r) (SimplePolygon p r)
-> SimplePolygon p r
forall s a. s -> Getting a s a -> a
^.Getting (SimplePolygon p r) (Polygon t p r) (SimplePolygon p r)
forall (t :: PolygonType) p r.
Lens' (Polygon t p r) (SimplePolygon p r)
outerBoundary
    es :: [LineSegment 2 p r]
es = (SimplePolygon p r -> [LineSegment 2 p r])
-> [SimplePolygon p r] -> [LineSegment 2 p r]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (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])
-> (SimplePolygon p r -> CircularVector (LineSegment 2 p r))
-> SimplePolygon p r
-> [LineSegment 2 p r]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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] -> [LineSegment 2 p r])
-> [SimplePolygon p r] -> [LineSegment 2 p r]
forall a b. (a -> b) -> a -> b
$ SimplePolygon p r
out SimplePolygon p r -> [SimplePolygon p r] -> [SimplePolygon p r]
forall a. a -> [a] -> [a]
: Polygon t p r -> [SimplePolygon p r]
forall (t :: PolygonType) p r.
Polygon t p r -> [Polygon 'Simple p r]
holeList Polygon t p r
pg

-- | Check if a point lies inside a polygon, on the boundary, or outside of the polygon.
-- Running time: O(n).
--
-- >>> Point2 1 1 `inPolygon` simplePoly
-- Inside
-- >>> Point2 0 0 `inPolygon` simplePoly
-- OnBoundary
-- >>> Point2 10 0 `inPolygon` simplePoly
-- OnBoundary
-- >>> Point2 5 13 `inPolygon` simplePoly
-- Inside
-- >>> Point2 5 10 `inPolygon` simplePoly
-- Inside
-- >>> Point2 10 5 `inPolygon` simplePoly
-- OnBoundary
-- >>> Point2 20 5 `inPolygon` simplePoly
-- Outside
--
-- TODO: Add some testcases with multiPolygons
-- TODO: Add some more onBoundary testcases
inPolygon             :: forall t p r. (Fractional r, Ord r)
                      => Point 2 r -> Polygon t p r -> PointLocationResult
Point 2 r
q inPolygon :: Point 2 r -> Polygon t p r -> PointLocationResult
`inPolygon` Polygon t p r
pg
  | Point 2 r
q Point 2 r -> Polygon t p r -> Bool
forall r (t :: PolygonType) p.
(Num r, Ord r) =>
Point 2 r -> Polygon t p r -> Bool
`onBoundary` Polygon t p r
pg = PointLocationResult
OnBoundary
  | Bool
inHole            = PointLocationResult
Outside
  | Bool
otherwise         = Point 2 r
q Point 2 r -> SimplePolygon p r -> PointLocationResult
forall p r.
(Fractional r, Ord r) =>
Point 2 r -> SimplePolygon p r -> PointLocationResult
`inPolygon'` (Polygon t p r
pgPolygon t p r
-> Getting (SimplePolygon p r) (Polygon t p r) (SimplePolygon p r)
-> SimplePolygon p r
forall s a. s -> Getting a s a -> a
^.Getting (SimplePolygon p r) (Polygon t p r) (SimplePolygon p r)
forall (t :: PolygonType) p r.
Lens' (Polygon t p r) (SimplePolygon p r)
outerBoundary)
  where
    inHole :: Bool
inHole = (SimplePolygon p r -> Bool) -> [SimplePolygon p r] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Point 2 r
q Point 2 r -> SimplePolygon p r -> Bool
forall r (t :: PolygonType) p.
(Fractional r, Ord r) =>
Point 2 r -> Polygon t p r -> Bool
`insidePolygon`) ([SimplePolygon p r] -> Bool) -> [SimplePolygon p r] -> Bool
forall a b. (a -> b) -> a -> b
$ Polygon t p r -> [SimplePolygon p r]
forall (t :: PolygonType) p r.
Polygon t p r -> [Polygon 'Simple p r]
holeList Polygon t p r
pg

-- | Returns true if the point lies in the polygon
-- pre: point lies inside or outside the polygon, not on its boundary.
inPolygon'        :: forall p r. (Fractional r, Ord r)
                  => Point 2 r -> SimplePolygon p r
                  -> PointLocationResult
Point 2 r
q inPolygon' :: Point 2 r -> SimplePolygon p r -> PointLocationResult
`inPolygon'` SimplePolygon p r
pg = if Int -> Bool
forall a. Integral a => a -> Bool
odd (Int -> Bool)
-> ([LineSegment 2 p r] -> Int) -> [LineSegment 2 p r] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Point 2 r] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Point 2 r] -> Int)
-> ([LineSegment 2 p r] -> [Point 2 r])
-> [LineSegment 2 p r]
-> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (LineSegment 2 p r -> Maybe (Point 2 r))
-> [LineSegment 2 p r] -> [Point 2 r]
forall a b. (a -> Maybe b) -> [a] -> [b]
mapMaybe LineSegment 2 p r -> Maybe (Point 2 r)
intersectionPoint ([LineSegment 2 p r] -> Bool) -> [LineSegment 2 p r] -> Bool
forall a b. (a -> b) -> a -> b
$ [LineSegment 2 p r]
ups [LineSegment 2 p r] -> [LineSegment 2 p r] -> [LineSegment 2 p r]
forall a. Semigroup a => a -> a -> a
<> [LineSegment 2 p r]
downs
                    then PointLocationResult
Inside else PointLocationResult
Outside
  where
    -- we don't care about horizontal edges
    ([LineSegment 2 p r]
ups',[LineSegment 2 p r]
_horizontals,[LineSegment 2 p r]
downs') = [LineSegment 2 p r]
-> ([LineSegment 2 p r], [LineSegment 2 p r], [LineSegment 2 p r])
partitionEdges ([LineSegment 2 p r]
 -> ([LineSegment 2 p r], [LineSegment 2 p r], [LineSegment 2 p r]))
-> (SimplePolygon p r -> [LineSegment 2 p r])
-> SimplePolygon p r
-> ([LineSegment 2 p r], [LineSegment 2 p r], [LineSegment 2 p r])
forall b c a. (b -> c) -> (a -> b) -> a -> c
. SimplePolygon p r -> [LineSegment 2 p r]
forall (t :: PolygonType) p r. Polygon t p r -> [LineSegment 2 p r]
listEdges (SimplePolygon p r
 -> ([LineSegment 2 p r], [LineSegment 2 p r], [LineSegment 2 p r]))
-> SimplePolygon p r
-> ([LineSegment 2 p r], [LineSegment 2 p r], [LineSegment 2 p r])
forall a b. (a -> b) -> a -> b
$ SimplePolygon p r
pg
    partitionEdges :: [LineSegment 2 p r]
-> ([LineSegment 2 p r], [LineSegment 2 p r], [LineSegment 2 p r])
partitionEdges = (LineSegment 2 p r -> Ordering)
-> [LineSegment 2 p r]
-> ([LineSegment 2 p r], [LineSegment 2 p r], [LineSegment 2 p r])
forall (f :: * -> *) a.
Foldable f =>
(a -> Ordering) -> f a -> ([a], [a], [a])
List.partition3 ((LineSegment 2 p r -> Ordering)
 -> [LineSegment 2 p r]
 -> ([LineSegment 2 p r], [LineSegment 2 p r], [LineSegment 2 p r]))
-> (LineSegment 2 p r -> Ordering)
-> [LineSegment 2 p r]
-> ([LineSegment 2 p r], [LineSegment 2 p r], [LineSegment 2 p r])
forall a b. (a -> b) -> a -> b
$ \LineSegment 2 p r
s -> (LineSegment 2 p r
sLineSegment 2 p r -> Getting r (LineSegment 2 p r) r -> r
forall s a. s -> Getting a s a -> a
^.((Point 2 r :+ EndExtra (LineSegment 2 p r))
 -> Const r (Point 2 r :+ EndExtra (LineSegment 2 p r)))
-> LineSegment 2 p r -> Const r (LineSegment 2 p r)
forall t. HasEnd t => Lens' t (EndCore t :+ EndExtra t)
end(((Point 2 r :+ EndExtra (LineSegment 2 p r))
  -> Const r (Point 2 r :+ EndExtra (LineSegment 2 p r)))
 -> LineSegment 2 p r -> Const r (LineSegment 2 p r))
-> ((r -> Const r r)
    -> (Point 2 r :+ EndExtra (LineSegment 2 p r))
    -> Const r (Point 2 r :+ EndExtra (LineSegment 2 p r)))
-> Getting r (LineSegment 2 p r) r
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(Point 2 r -> Const r (Point 2 r))
-> (Point 2 r :+ EndExtra (LineSegment 2 p r))
-> Const r (Point 2 r :+ EndExtra (LineSegment 2 p r))
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core((Point 2 r -> Const r (Point 2 r))
 -> (Point 2 r :+ EndExtra (LineSegment 2 p r))
 -> Const r (Point 2 r :+ EndExtra (LineSegment 2 p r)))
-> ((r -> Const r r) -> Point 2 r -> Const r (Point 2 r))
-> (r -> Const r r)
-> (Point 2 r :+ EndExtra (LineSegment 2 p r))
-> Const r (Point 2 r :+ EndExtra (LineSegment 2 p r))
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(r -> Const r r) -> Point 2 r -> Const r (Point 2 r)
forall (d :: Nat) (point :: Nat -> * -> *) r.
(2 <= d, Arity d, AsAPoint point) =>
Lens' (point d r) r
yCoord) r -> r -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` (LineSegment 2 p r
sLineSegment 2 p r -> Getting r (LineSegment 2 p r) r -> r
forall s a. s -> Getting a s a -> a
^.((Point 2 r :+ StartExtra (LineSegment 2 p r))
 -> Const r (Point 2 r :+ StartExtra (LineSegment 2 p r)))
-> LineSegment 2 p r -> Const r (LineSegment 2 p r)
forall t. HasStart t => Lens' t (StartCore t :+ StartExtra t)
start(((Point 2 r :+ StartExtra (LineSegment 2 p r))
  -> Const r (Point 2 r :+ StartExtra (LineSegment 2 p r)))
 -> LineSegment 2 p r -> Const r (LineSegment 2 p r))
-> ((r -> Const r r)
    -> (Point 2 r :+ StartExtra (LineSegment 2 p r))
    -> Const r (Point 2 r :+ StartExtra (LineSegment 2 p r)))
-> Getting r (LineSegment 2 p r) r
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(Point 2 r -> Const r (Point 2 r))
-> (Point 2 r :+ StartExtra (LineSegment 2 p r))
-> Const r (Point 2 r :+ StartExtra (LineSegment 2 p r))
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core((Point 2 r -> Const r (Point 2 r))
 -> (Point 2 r :+ StartExtra (LineSegment 2 p r))
 -> Const r (Point 2 r :+ StartExtra (LineSegment 2 p r)))
-> ((r -> Const r r) -> Point 2 r -> Const r (Point 2 r))
-> (r -> Const r r)
-> (Point 2 r :+ StartExtra (LineSegment 2 p r))
-> Const r (Point 2 r :+ StartExtra (LineSegment 2 p r))
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(r -> Const r r) -> Point 2 r -> Const r (Point 2 r)
forall (d :: Nat) (point :: Nat -> * -> *) r.
(2 <= d, Arity d, AsAPoint point) =>
Lens' (point d r) r
yCoord)

    -- upward edges include start, exclude end
    ups :: [LineSegment 2 p r]
ups   = (LineSegment 2 p r -> LineSegment 2 p r)
-> [LineSegment 2 p r] -> [LineSegment 2 p r]
forall a b. (a -> b) -> [a] -> [b]
map (\(LineSegment' Point 2 r :+ p
a Point 2 r :+ p
b) -> EndPoint (Point 2 r :+ p)
-> EndPoint (Point 2 r :+ p) -> LineSegment 2 p r
forall (d :: Nat) r p.
EndPoint (Point d r :+ p)
-> EndPoint (Point d r :+ p) -> LineSegment d p r
LineSegment ((Point 2 r :+ p) -> EndPoint (Point 2 r :+ p)
forall a. a -> EndPoint a
Closed Point 2 r :+ p
a) ((Point 2 r :+ p) -> EndPoint (Point 2 r :+ p)
forall a. a -> EndPoint a
Open Point 2 r :+ p
b)) [LineSegment 2 p r]
ups'
    -- downward edges exclude start, include end
    downs :: [LineSegment 2 p r]
downs = (LineSegment 2 p r -> LineSegment 2 p r)
-> [LineSegment 2 p r] -> [LineSegment 2 p r]
forall a b. (a -> b) -> [a] -> [b]
map (\(LineSegment' Point 2 r :+ p
a Point 2 r :+ p
b) -> EndPoint (Point 2 r :+ p)
-> EndPoint (Point 2 r :+ p) -> LineSegment 2 p r
forall (d :: Nat) r p.
EndPoint (Point d r :+ p)
-> EndPoint (Point d r :+ p) -> LineSegment d p r
LineSegment ((Point 2 r :+ p) -> EndPoint (Point 2 r :+ p)
forall a. a -> EndPoint a
Open Point 2 r :+ p
a) ((Point 2 r :+ p) -> EndPoint (Point 2 r :+ p)
forall a. a -> EndPoint a
Closed Point 2 r :+ p
b)) [LineSegment 2 p r]
downs'

    -- Given an edge, compute the intersection point (if a point) with
    -- the line through the query point, and test if it lies strictly
    -- right of q.
    --
    -- See http://geomalgorithms.com/a03-_inclusion.html for more information.
    intersectionPoint :: LineSegment 2 p r -> Maybe (Point 2 r)
intersectionPoint =  (Point 2 r -> Bool) -> Maybe (Point 2 r) -> Maybe (Point 2 r)
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
F.find (\Point 2 r
p -> Point 2 r
pPoint 2 r
-> ((r -> Const r r) -> Point 2 r -> Const r (Point 2 r)) -> r
forall s a. s -> Getting a s a -> a
^.(r -> Const r r) -> Point 2 r -> Const r (Point 2 r)
forall (d :: Nat) (point :: Nat -> * -> *) r.
(1 <= d, Arity d, AsAPoint point) =>
Lens' (point d r) r
xCoord r -> r -> Bool
forall a. Ord a => a -> a -> Bool
> Point 2 r
qPoint 2 r
-> ((r -> Const r r) -> Point 2 r -> Const r (Point 2 r)) -> r
forall s a. s -> Getting a s a -> a
^.(r -> Const r r) -> Point 2 r -> Const r (Point 2 r)
forall (d :: Nat) (point :: Nat -> * -> *) r.
(1 <= d, Arity d, AsAPoint point) =>
Lens' (point d r) r
xCoord) (Maybe (Point 2 r) -> Maybe (Point 2 r))
-> (LineSegment 2 p r -> Maybe (Point 2 r))
-> LineSegment 2 p r
-> Maybe (Point 2 r)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (ts :: [*]).
NatToInt (RIndex (Point 2 r) ts) =>
CoRec Identity ts -> Maybe (Point 2 r)
forall t (ts :: [*]).
NatToInt (RIndex t ts) =>
CoRec Identity ts -> Maybe t
asA @(Point 2 r) (CoRec Identity '[NoIntersection, Point 2 r, LineSegment 2 p r]
 -> Maybe (Point 2 r))
-> (LineSegment 2 p r
    -> CoRec Identity '[NoIntersection, Point 2 r, LineSegment 2 p r])
-> LineSegment 2 p r
-> Maybe (Point 2 r)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (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
`intersect` Line 2 r
l)
    l :: Line 2 r
l = r -> Line 2 r
forall r. Num r => r -> Line 2 r
horizontalLine (r -> Line 2 r) -> r -> Line 2 r
forall a b. (a -> b) -> a -> b
$ Point 2 r
qPoint 2 r
-> ((r -> Const r r) -> Point 2 r -> Const r (Point 2 r)) -> r
forall s a. s -> Getting a s a -> a
^.(r -> Const r r) -> Point 2 r -> Const r (Point 2 r)
forall (d :: Nat) (point :: Nat -> * -> *) r.
(2 <= d, Arity d, AsAPoint point) =>
Lens' (point d r) r
yCoord


-- | Test if a point lies strictly inside the polgyon.
insidePolygon        :: (Fractional r, Ord r) => Point 2 r -> Polygon t p r -> Bool
Point 2 r
q insidePolygon :: Point 2 r -> Polygon t p r -> Bool
`insidePolygon` Polygon t p r
pg = Point 2 r
q Point 2 r -> Polygon t p r -> PointLocationResult
forall (t :: PolygonType) p r.
(Fractional r, Ord r) =>
Point 2 r -> Polygon t p r -> PointLocationResult
`inPolygon` Polygon t p r
pg PointLocationResult -> PointLocationResult -> Bool
forall a. Eq a => a -> a -> Bool
== PointLocationResult
Inside


-- testQ = map (`inPolygon` testPoly) [ Point2 1 1    -- Inside
--                                    , Point2 0 0    -- OnBoundary
--                                    , Point2 5 14   -- Inside
--                                    , Point2 5 10   -- Inside
--                                    , Point2 10 5   -- OnBoundary
--                                    , Point2 20 5   -- Outside
--                                    ]

-- testPoly :: SimplePolygon () Rational
-- testPoly = fromPoints . map ext $ [ Point2 0 0
--                                                   , Point2 10 0
--                                                   , Point2 10 10
--                                                   , Point2 5 15
--                                                   , Point2 1 11
--                                                   ]