-- |

-- Module:      Data.Geo.Jord.Polygon

-- Copyright:   (c) 2020 Cedric Liegeois

-- License:     BSD3

-- Maintainer:  Cedric Liegeois <ofmooseandmen@yahoo.fr>

-- Stability:   experimental

-- Portability: portable

--

-- Types and functions for working with polygons at the surface of a __spherical__ celestial body.

--

-- In order to use this module you should start with the following imports:

--

-- @

-- import qualified Data.Geo.Jord.Geodetic as Geodetic

-- import qualified Data.Geo.Jord.Polygon as Polygon

-- @

module Data.Geo.Jord.Polygon
    (
    -- * The 'Polygon' type

      Polygon
    , vertices
    , edges
    , concave
    -- * smart constructors

    , Error(..)
    , simple
    , circle
    , arc
    -- * calculations

    , contains
    , triangulate
    ) where

import Data.List (find)
import Data.Maybe (isJust, mapMaybe)

import Data.Geo.Jord.Angle (Angle)
import qualified Data.Geo.Jord.Angle as Angle
import Data.Geo.Jord.Ellipsoid (meanRadius)
import Data.Geo.Jord.Geodetic (HorizontalPosition)
import qualified Data.Geo.Jord.Geodetic as Geodetic
import Data.Geo.Jord.GreatCircle (MinorArc)
import qualified Data.Geo.Jord.GreatCircle as GreatCircle
import Data.Geo.Jord.Length (Length)
import qualified Data.Geo.Jord.Length as Length
import qualified Data.Geo.Jord.Math3d as Math3d
import Data.Geo.Jord.Model (Spherical, surface)
import Data.Geo.Jord.Triangle (Triangle)
import qualified Data.Geo.Jord.Triangle as Triangle

-- | A polygon whose vertices are horizontal geodetic positions.

data Polygon a =
    Polygon
        { Polygon a -> [HorizontalPosition a]
vertices :: [HorizontalPosition a] -- ^ vertices of the polygon in __clockwise__ order.

        , Polygon a -> [MinorArc a]
edges :: [MinorArc a] -- ^ edges of the polyon in __clockwise__ order.

        , Polygon a -> Bool
concave :: Bool -- ^ whether the polygon is concave.

        }
    deriving (Polygon a -> Polygon a -> Bool
(Polygon a -> Polygon a -> Bool)
-> (Polygon a -> Polygon a -> Bool) -> Eq (Polygon a)
forall a. Model a => Polygon a -> Polygon a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Polygon a -> Polygon a -> Bool
$c/= :: forall a. Model a => Polygon a -> Polygon a -> Bool
== :: Polygon a -> Polygon a -> Bool
$c== :: forall a. Model a => Polygon a -> Polygon a -> Bool
Eq, Int -> Polygon a -> ShowS
[Polygon a] -> ShowS
Polygon a -> String
(Int -> Polygon a -> ShowS)
-> (Polygon a -> String)
-> ([Polygon a] -> ShowS)
-> Show (Polygon a)
forall a. Spherical a => Int -> Polygon a -> ShowS
forall a. Spherical a => [Polygon a] -> ShowS
forall a. Spherical a => Polygon a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Polygon a] -> ShowS
$cshowList :: forall a. Spherical a => [Polygon a] -> ShowS
show :: Polygon a -> String
$cshow :: forall a. Spherical a => Polygon a -> String
showsPrec :: Int -> Polygon a -> ShowS
$cshowsPrec :: forall a. Spherical a => Int -> Polygon a -> ShowS
Show)

-- | Error returned when attempting to create a polygon from invalid data. 

data Error
    = NotEnoughVertices -- ^ less than 3 vertices were supplied.

    | InvalidEdge -- ^ 2 consecutives vertices are antipodal or equal.

    | InvalidRadius -- ^ radius of circle or arc is <= 0.

    | EmptyArcRange -- ^ arc start angle == end angle.

    | SeflIntersectingEdge -- ^ 2 edges of the polygon intersect.

    deriving (Error -> Error -> Bool
(Error -> Error -> Bool) -> (Error -> Error -> Bool) -> Eq Error
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Error -> Error -> Bool
$c/= :: Error -> Error -> Bool
== :: Error -> Error -> Bool
$c== :: Error -> Error -> Bool
Eq, Int -> Error -> ShowS
[Error] -> ShowS
Error -> String
(Int -> Error -> ShowS)
-> (Error -> String) -> ([Error] -> ShowS) -> Show Error
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Error] -> ShowS
$cshowList :: [Error] -> ShowS
show :: Error -> String
$cshow :: Error -> String
showsPrec :: Int -> Error -> ShowS
$cshowsPrec :: Int -> Error -> ShowS
Show)

-- | Simple polygon (outer ring only and not self-intersecting) from given vertices. Returns an error ('Left') if:

--

--     * less than 3 vertices are given.

--     * the given vertices defines self-intersecting edges.

--     * the given vertices contains duplicated positions or antipodal positions.

simple :: (Spherical a) => [HorizontalPosition a] -> Either Error (Polygon a)
simple :: [HorizontalPosition a] -> Either Error (Polygon a)
simple [HorizontalPosition a]
vs
    | [HorizontalPosition a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [HorizontalPosition a]
vs = Error -> Either Error (Polygon a)
forall a b. a -> Either a b
Left Error
NotEnoughVertices
    | [HorizontalPosition a] -> HorizontalPosition a
forall a. [a] -> a
head [HorizontalPosition a]
vs HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
== [HorizontalPosition a] -> HorizontalPosition a
forall a. [a] -> a
last [HorizontalPosition a]
vs = [HorizontalPosition a] -> Either Error (Polygon a)
forall a.
Spherical a =>
[HorizontalPosition a] -> Either Error (Polygon a)
simple ([HorizontalPosition a] -> [HorizontalPosition a]
forall a. [a] -> [a]
init [HorizontalPosition a]
vs)
    | [HorizontalPosition a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [HorizontalPosition a]
vs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
3 = Error -> Either Error (Polygon a)
forall a b. a -> Either a b
Left Error
NotEnoughVertices
    | Bool
otherwise = [HorizontalPosition a] -> Either Error (Polygon a)
forall a.
Spherical a =>
[HorizontalPosition a] -> Either Error (Polygon a)
simple' [HorizontalPosition a]
vs

-- | Circle from given centre and radius. The resulting polygon contains @nb@ vertices equally distant from one

-- another. Returns an error ('Left') if:

--

--     * given radius is 0

--     * given number of positions is less than 3

circle :: (Spherical a) => HorizontalPosition a -> Length -> Int -> Either Error (Polygon a)
circle :: HorizontalPosition a -> Length -> Int -> Either Error (Polygon a)
circle HorizontalPosition a
c Length
r Int
nb
    | Length
r Length -> Length -> Bool
forall a. Ord a => a -> a -> Bool
<= Length
Length.zero = Error -> Either Error (Polygon a)
forall a b. a -> Either a b
Left Error
InvalidRadius
    | Int
nb Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
3 = Error -> Either Error (Polygon a)
forall a b. a -> Either a b
Left Error
NotEnoughVertices
    | Bool
otherwise = Polygon a -> Either Error (Polygon a)
forall a b. b -> Either a b
Right (HorizontalPosition a -> Length -> [Double] -> Polygon a
forall a.
Spherical a =>
HorizontalPosition a -> Length -> [Double] -> Polygon a
discretiseArc HorizontalPosition a
c Length
r [Double]
as)
  where
    n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
nb :: Double
    as :: [Double]
as = Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
nb ((Double -> Double) -> Double -> [Double]
forall a. (a -> a) -> a -> [a]
iterate (\Double
x -> Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
2.0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
n) Double
0.0)

-- | Arc from given centre, radius and given start/end angles. The resulting polygon contains @nb@ vertices equally

-- distant from one another. Returns an error ('Left') if:

--

--     * given radius is 0

--     * given number of positions is less than 3

--     * difference between start and end angle is 0

arc :: (Spherical a)
    => HorizontalPosition a
    -> Length
    -> Angle
    -> Angle
    -> Int
    -> Either Error (Polygon a)
arc :: HorizontalPosition a
-> Length -> Angle -> Angle -> Int -> Either Error (Polygon a)
arc HorizontalPosition a
c Length
r Angle
sa Angle
ea Int
nb
    | Length
r Length -> Length -> Bool
forall a. Ord a => a -> a -> Bool
<= Length
Length.zero = Error -> Either Error (Polygon a)
forall a b. a -> Either a b
Left Error
InvalidRadius
    | Int
nb Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
3 = Error -> Either Error (Polygon a)
forall a b. a -> Either a b
Left Error
NotEnoughVertices
    | Angle
range Angle -> Angle -> Bool
forall a. Eq a => a -> a -> Bool
== Angle
Angle.zero = Error -> Either Error (Polygon a)
forall a b. a -> Either a b
Left Error
EmptyArcRange
    | Bool
otherwise = Polygon a -> Either Error (Polygon a)
forall a b. b -> Either a b
Right (HorizontalPosition a -> Length -> [Double] -> Polygon a
forall a.
Spherical a =>
HorizontalPosition a -> Length -> [Double] -> Polygon a
discretiseArc HorizontalPosition a
c Length
r [Double]
as)
  where
    range :: Angle
range = Angle -> Angle -> Angle
Angle.clockwiseDifference Angle
sa Angle
ea
    n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
nb :: Double
    inc :: Double
inc = Angle -> Double
Angle.toRadians Angle
range Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1.0)
    r0 :: Double
r0 = Angle -> Double
Angle.toRadians Angle
sa
    as :: [Double]
as = Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
nb ((Double -> Double) -> Double -> [Double]
forall a. (a -> a) -> a -> [a]
iterate (Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
inc) Double
r0)

-- | @contains poly p@ returns 'True' if position @p@ is enclosed by the vertices of polygon

-- @poly@ - see 'GreatCircle.enclosedBy'.

contains :: (Spherical a) => Polygon a -> HorizontalPosition a -> Bool
contains :: Polygon a -> HorizontalPosition a -> Bool
contains Polygon a
poly HorizontalPosition a
p = HorizontalPosition a -> [HorizontalPosition a] -> Bool
forall a.
Spherical a =>
HorizontalPosition a -> [HorizontalPosition a] -> Bool
GreatCircle.enclosedBy HorizontalPosition a
p (Polygon a -> [HorizontalPosition a]
forall a. Polygon a -> [HorizontalPosition a]
vertices Polygon a
poly)

-- | Triangulates the given polygon using the ear clipping method.

--

-- May return an empty list if the algorithm fails to find an ear (which probably indicates a bug in the implementation).

triangulate :: (Spherical a) => Polygon a -> [Triangle a]
triangulate :: Polygon a -> [Triangle a]
triangulate Polygon a
p
    | [HorizontalPosition a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [HorizontalPosition a]
vs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
3 = [[HorizontalPosition a] -> Triangle a
forall a. Spherical a => [HorizontalPosition a] -> Triangle a
triangle [HorizontalPosition a]
vs]
    | Bool
otherwise = [HorizontalPosition a] -> [Triangle a] -> [Triangle a]
forall a.
Spherical a =>
[HorizontalPosition a] -> [Triangle a] -> [Triangle a]
earClipping [HorizontalPosition a]
vs []
  where
    vs :: [HorizontalPosition a]
vs = Polygon a -> [HorizontalPosition a]
forall a. Polygon a -> [HorizontalPosition a]
vertices Polygon a
p

-- private

triangle :: (Spherical a) => [HorizontalPosition a] -> Triangle a
triangle :: [HorizontalPosition a] -> Triangle a
triangle [HorizontalPosition a]
vs = HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Triangle a
forall a.
Spherical a =>
HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Triangle a
Triangle.unsafeMake ([HorizontalPosition a] -> HorizontalPosition a
forall a. [a] -> a
head [HorizontalPosition a]
vs) ([HorizontalPosition a]
vs [HorizontalPosition a] -> Int -> HorizontalPosition a
forall a. [a] -> Int -> a
!! Int
1) ([HorizontalPosition a]
vs [HorizontalPosition a] -> Int -> HorizontalPosition a
forall a. [a] -> Int -> a
!! Int
2)

earClipping :: (Spherical a) => [HorizontalPosition a] -> [Triangle a] -> [Triangle a]
earClipping :: [HorizontalPosition a] -> [Triangle a] -> [Triangle a]
earClipping [HorizontalPosition a]
vs [Triangle a]
ts
    | [HorizontalPosition a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [HorizontalPosition a]
vs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
3 = [Triangle a] -> [Triangle a]
forall a. [a] -> [a]
reverse ([HorizontalPosition a] -> Triangle a
forall a. Spherical a => [HorizontalPosition a] -> Triangle a
triangle [HorizontalPosition a]
vs Triangle a -> [Triangle a] -> [Triangle a]
forall a. a -> [a] -> [a]
: [Triangle a]
ts)
    | Bool
otherwise =
        case [HorizontalPosition a]
-> Maybe
     (HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
forall a.
Spherical a =>
[HorizontalPosition a]
-> Maybe
     (HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
findEar [HorizontalPosition a]
vs of
            Maybe
  (HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
Nothing -> []
            (Just (HorizontalPosition a
p, HorizontalPosition a
e, HorizontalPosition a
n)) -> [HorizontalPosition a] -> [Triangle a] -> [Triangle a]
forall a.
Spherical a =>
[HorizontalPosition a] -> [Triangle a] -> [Triangle a]
earClipping [HorizontalPosition a]
vs' [Triangle a]
ts'
                where ts' :: [Triangle a]
ts' = HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Triangle a
forall a.
Spherical a =>
HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Triangle a
Triangle.unsafeMake HorizontalPosition a
p HorizontalPosition a
e HorizontalPosition a
n Triangle a -> [Triangle a] -> [Triangle a]
forall a. a -> [a] -> [a]
: [Triangle a]
ts
                      vs' :: [HorizontalPosition a]
vs' = (HorizontalPosition a -> Bool)
-> [HorizontalPosition a] -> [HorizontalPosition a]
forall a. (a -> Bool) -> [a] -> [a]
filter (HorizontalPosition a -> HorizontalPosition a -> Bool
forall a. Eq a => a -> a -> Bool
/= HorizontalPosition a
e) [HorizontalPosition a]
vs

findEar ::
       (Spherical a)
    => [HorizontalPosition a]
    -> Maybe (HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
findEar :: [HorizontalPosition a]
-> Maybe
     (HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
findEar [HorizontalPosition a]
ps = ((HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
 -> Bool)
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
-> Maybe
     (HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
find ((HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
-> [HorizontalPosition a] -> Bool
forall a.
Spherical a =>
(HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
-> [HorizontalPosition a] -> Bool
`isEar` [HorizontalPosition a]
rs) [(HorizontalPosition a, HorizontalPosition a,
  HorizontalPosition a)]
convex
  where
    rs :: [HorizontalPosition a]
rs = [HorizontalPosition a] -> [HorizontalPosition a]
forall a.
Spherical a =>
[HorizontalPosition a] -> [HorizontalPosition a]
reflices [HorizontalPosition a]
ps
    t3 :: [(HorizontalPosition a, HorizontalPosition a,
  HorizontalPosition a)]
t3 = [HorizontalPosition a]
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
forall a.
Spherical a =>
[HorizontalPosition a]
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
tuple3 [HorizontalPosition a]
ps
    convex :: [(HorizontalPosition a, HorizontalPosition a,
  HorizontalPosition a)]
convex = ((HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
 -> Bool)
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
forall a. (a -> Bool) -> [a] -> [a]
filter (\(HorizontalPosition a
_, HorizontalPosition a
v, HorizontalPosition a
_) -> HorizontalPosition a
v HorizontalPosition a -> [HorizontalPosition a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [HorizontalPosition a]
rs) [(HorizontalPosition a, HorizontalPosition a,
  HorizontalPosition a)]
t3

-- | a convex vertex @c@ is an ear if triangle (prev, c, next) contains no reflex.

isEar ::
       (Spherical a)
    => (HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
    -> [HorizontalPosition a]
    -> Bool
isEar :: (HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
-> [HorizontalPosition a] -> Bool
isEar (HorizontalPosition a
p, HorizontalPosition a
c, HorizontalPosition a
n) = (HorizontalPosition a -> Bool) -> [HorizontalPosition a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\HorizontalPosition a
r -> Bool -> Bool
not (HorizontalPosition a -> [HorizontalPosition a] -> Bool
forall a.
Spherical a =>
HorizontalPosition a -> [HorizontalPosition a] -> Bool
GreatCircle.enclosedBy HorizontalPosition a
r [HorizontalPosition a]
vs))
  where
    vs :: [HorizontalPosition a]
vs = [HorizontalPosition a
p, HorizontalPosition a
c, HorizontalPosition a
n]

-- | A reflex is a vertex where the polygon is concave.

-- a vertex is a reflex if previous vertex is left (assuming a clockwise polygon), otherwise it is a convex vertex

reflices :: (Spherical a) => [HorizontalPosition a] -> [HorizontalPosition a]
reflices :: [HorizontalPosition a] -> [HorizontalPosition a]
reflices [HorizontalPosition a]
ps = ((HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
 -> HorizontalPosition a)
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
-> [HorizontalPosition a]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\(HorizontalPosition a
_, HorizontalPosition a
c, HorizontalPosition a
_) -> HorizontalPosition a
c) [(HorizontalPosition a, HorizontalPosition a,
  HorizontalPosition a)]
rs
  where
    t3 :: [(HorizontalPosition a, HorizontalPosition a,
  HorizontalPosition a)]
t3 = [HorizontalPosition a]
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
forall a.
Spherical a =>
[HorizontalPosition a]
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
tuple3 [HorizontalPosition a]
ps
    rs :: [(HorizontalPosition a, HorizontalPosition a,
  HorizontalPosition a)]
rs = ((HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
 -> Bool)
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
forall a. (a -> Bool) -> [a] -> [a]
filter (\(HorizontalPosition a
p, HorizontalPosition a
c, HorizontalPosition a
n) -> HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Side
forall a.
Spherical a =>
HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Side
GreatCircle.side HorizontalPosition a
p HorizontalPosition a
c HorizontalPosition a
n Side -> Side -> Bool
forall a. Eq a => a -> a -> Bool
== Side
GreatCircle.LeftOf) [(HorizontalPosition a, HorizontalPosition a,
  HorizontalPosition a)]
t3

-- | [mAB, mBC, mCD, mDE, mEA]

-- no intersections:

-- mAB vs [mCD, mDE]

-- mBC vs [mDE, mEA]

-- mCD vs [mEA]

selfIntersects :: (Spherical a) => [MinorArc a] -> Bool
selfIntersects :: [MinorArc a] -> Bool
selfIntersects [MinorArc a]
ps
    | [MinorArc a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [MinorArc a]
ps Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
4 = Bool
False
    | Bool
otherwise = ((MinorArc a, [MinorArc a]) -> Bool)
-> [(MinorArc a, [MinorArc a])] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (MinorArc a, [MinorArc a]) -> Bool
forall a. Spherical a => (MinorArc a, [MinorArc a]) -> Bool
intersects [(MinorArc a, [MinorArc a])]
pairs
  where
    ([MinorArc a]
_, [(MinorArc a, [MinorArc a])]
pairs) = ([MinorArc a], [(MinorArc a, [MinorArc a])])
-> ([MinorArc a], [(MinorArc a, [MinorArc a])])
forall a.
Spherical a =>
([MinorArc a], [(MinorArc a, [MinorArc a])])
-> ([MinorArc a], [(MinorArc a, [MinorArc a])])
makePairs' ([MinorArc a]
ps, [])

intersects :: (Spherical a) => (MinorArc a, [MinorArc a]) -> Bool
intersects :: (MinorArc a, [MinorArc a]) -> Bool
intersects (MinorArc a
ma, [MinorArc a]
mas) = (MinorArc a -> Bool) -> [MinorArc a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Maybe (HorizontalPosition a) -> Bool
forall a. Maybe a -> Bool
isJust (Maybe (HorizontalPosition a) -> Bool)
-> (MinorArc a -> Maybe (HorizontalPosition a))
-> MinorArc a
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. MinorArc a -> MinorArc a -> Maybe (HorizontalPosition a)
forall a.
Spherical a =>
MinorArc a -> MinorArc a -> Maybe (HorizontalPosition a)
GreatCircle.intersection MinorArc a
ma) [MinorArc a]
mas

makePairs' ::
       (Spherical a)
    => ([MinorArc a], [(MinorArc a, [MinorArc a])])
    -> ([MinorArc a], [(MinorArc a, [MinorArc a])])
makePairs' :: ([MinorArc a], [(MinorArc a, [MinorArc a])])
-> ([MinorArc a], [(MinorArc a, [MinorArc a])])
makePairs' ([MinorArc a]
xs, [(MinorArc a, [MinorArc a])]
ps)
    | [MinorArc a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [MinorArc a]
xs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
3 = ([MinorArc a]
xs, [(MinorArc a, [MinorArc a])]
ps)
    | Bool
otherwise = ([MinorArc a], [(MinorArc a, [MinorArc a])])
-> ([MinorArc a], [(MinorArc a, [MinorArc a])])
forall a.
Spherical a =>
([MinorArc a], [(MinorArc a, [MinorArc a])])
-> ([MinorArc a], [(MinorArc a, [MinorArc a])])
makePairs' ([MinorArc a]
nxs, (MinorArc a, [MinorArc a])
np (MinorArc a, [MinorArc a])
-> [(MinorArc a, [MinorArc a])] -> [(MinorArc a, [MinorArc a])]
forall a. a -> [a] -> [a]
: [(MinorArc a, [MinorArc a])]
ps)
  where
    nxs :: [MinorArc a]
nxs = [MinorArc a] -> [MinorArc a]
forall a. [a] -> [a]
tail [MinorArc a]
xs
    -- if ps is empty (first call), drop last minor arc as it connect to first minor arc

    versus :: [MinorArc a]
versus =
        if [(MinorArc a, [MinorArc a])] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(MinorArc a, [MinorArc a])]
ps
            then [MinorArc a] -> [MinorArc a]
forall a. [a] -> [a]
init ([MinorArc a] -> [MinorArc a])
-> ([MinorArc a] -> [MinorArc a]) -> [MinorArc a] -> [MinorArc a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [MinorArc a] -> [MinorArc a]
forall a. [a] -> [a]
tail ([MinorArc a] -> [MinorArc a])
-> ([MinorArc a] -> [MinorArc a]) -> [MinorArc a] -> [MinorArc a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [MinorArc a] -> [MinorArc a]
forall a. [a] -> [a]
tail ([MinorArc a] -> [MinorArc a]) -> [MinorArc a] -> [MinorArc a]
forall a b. (a -> b) -> a -> b
$ [MinorArc a]
xs
            else [MinorArc a] -> [MinorArc a]
forall a. [a] -> [a]
tail ([MinorArc a] -> [MinorArc a])
-> ([MinorArc a] -> [MinorArc a]) -> [MinorArc a] -> [MinorArc a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [MinorArc a] -> [MinorArc a]
forall a. [a] -> [a]
tail ([MinorArc a] -> [MinorArc a]) -> [MinorArc a] -> [MinorArc a]
forall a b. (a -> b) -> a -> b
$ [MinorArc a]
xs
    np :: (MinorArc a, [MinorArc a])
np = ([MinorArc a] -> MinorArc a
forall a. [a] -> a
head [MinorArc a]
xs, [MinorArc a]
versus)

simple' :: (Spherical a) => [HorizontalPosition a] -> Either Error (Polygon a)
simple' :: [HorizontalPosition a] -> Either Error (Polygon a)
simple' [HorizontalPosition a]
vs
    | [MinorArc a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [MinorArc a]
es Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [HorizontalPosition a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [HorizontalPosition a]
vs = Error -> Either Error (Polygon a)
forall a b. a -> Either a b
Left Error
InvalidEdge
    | Bool
si = Error -> Either Error (Polygon a)
forall a b. a -> Either a b
Left Error
SeflIntersectingEdge
    | Bool
otherwise = Polygon a -> Either Error (Polygon a)
forall a b. b -> Either a b
Right ([HorizontalPosition a] -> [MinorArc a] -> Bool -> Polygon a
forall a.
[HorizontalPosition a] -> [MinorArc a] -> Bool -> Polygon a
Polygon [HorizontalPosition a]
os [MinorArc a]
es Bool
isConcave)
  where
    zs :: [(HorizontalPosition a, HorizontalPosition a,
  HorizontalPosition a)]
zs = [HorizontalPosition a]
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
forall a.
Spherical a =>
[HorizontalPosition a]
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
tuple3 [HorizontalPosition a]
vs
    clockwise :: Bool
clockwise = [Double] -> Double
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (((HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
 -> Double)
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
-> [Double]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\(HorizontalPosition a
a, HorizontalPosition a
b, HorizontalPosition a
c) -> Angle -> Double
Angle.toRadians (HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Angle
forall a.
Spherical a =>
HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Angle
GreatCircle.turn HorizontalPosition a
a HorizontalPosition a
b HorizontalPosition a
c)) [(HorizontalPosition a, HorizontalPosition a,
  HorizontalPosition a)]
zs) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.0
    os :: [HorizontalPosition a]
os =
        if Bool
clockwise
            then [HorizontalPosition a]
vs
            else [HorizontalPosition a] -> [HorizontalPosition a]
forall a. [a] -> [a]
reverse [HorizontalPosition a]
vs
    es :: [MinorArc a]
es = [HorizontalPosition a] -> [MinorArc a]
forall a. Spherical a => [HorizontalPosition a] -> [MinorArc a]
mkEdges [HorizontalPosition a]
os
    zzs :: [(HorizontalPosition a, HorizontalPosition a,
  HorizontalPosition a)]
zzs = [HorizontalPosition a]
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
forall a.
Spherical a =>
[HorizontalPosition a]
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
tuple3 [HorizontalPosition a]
os
    isConcave :: Bool
isConcave =
        [HorizontalPosition a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [HorizontalPosition a]
vs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
3 Bool -> Bool -> Bool
&& ((HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)
 -> Bool)
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
-> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\(HorizontalPosition a
a, HorizontalPosition a
b, HorizontalPosition a
c) -> HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Side
forall a.
Spherical a =>
HorizontalPosition a
-> HorizontalPosition a -> HorizontalPosition a -> Side
GreatCircle.side HorizontalPosition a
a HorizontalPosition a
b HorizontalPosition a
c Side -> Side -> Bool
forall a. Eq a => a -> a -> Bool
== Side
GreatCircle.LeftOf) [(HorizontalPosition a, HorizontalPosition a,
  HorizontalPosition a)]
zzs
    si :: Bool
si = Bool
isConcave Bool -> Bool -> Bool
&& [MinorArc a] -> Bool
forall a. Spherical a => [MinorArc a] -> Bool
selfIntersects [MinorArc a]
es

tuple3 ::
       (Spherical a)
    => [HorizontalPosition a]
    -> [(HorizontalPosition a, HorizontalPosition a, HorizontalPosition a)]
tuple3 :: [HorizontalPosition a]
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
tuple3 [HorizontalPosition a]
ps = [HorizontalPosition a]
-> [HorizontalPosition a]
-> [HorizontalPosition a]
-> [(HorizontalPosition a, HorizontalPosition a,
     HorizontalPosition a)]
forall a b c. [a] -> [b] -> [c] -> [(a, b, c)]
zip3 [HorizontalPosition a]
l1 [HorizontalPosition a]
l2 [HorizontalPosition a]
l3
  where
    l1 :: [HorizontalPosition a]
l1 = [HorizontalPosition a] -> HorizontalPosition a
forall a. [a] -> a
last [HorizontalPosition a]
ps HorizontalPosition a
-> [HorizontalPosition a] -> [HorizontalPosition a]
forall a. a -> [a] -> [a]
: [HorizontalPosition a] -> [HorizontalPosition a]
forall a. [a] -> [a]
init [HorizontalPosition a]
ps
    l2 :: [HorizontalPosition a]
l2 = [HorizontalPosition a]
ps
    l3 :: [HorizontalPosition a]
l3 = [HorizontalPosition a] -> [HorizontalPosition a]
forall a. [a] -> [a]
tail [HorizontalPosition a]
ps [HorizontalPosition a]
-> [HorizontalPosition a] -> [HorizontalPosition a]
forall a. [a] -> [a] -> [a]
++ [[HorizontalPosition a] -> HorizontalPosition a
forall a. [a] -> a
head [HorizontalPosition a]
ps]

mkEdges :: (Spherical a) => [HorizontalPosition a] -> [MinorArc a]
mkEdges :: [HorizontalPosition a] -> [MinorArc a]
mkEdges [HorizontalPosition a]
ps = ((HorizontalPosition a, HorizontalPosition a)
 -> Maybe (MinorArc a))
-> [(HorizontalPosition a, HorizontalPosition a)] -> [MinorArc a]
forall a b. (a -> Maybe b) -> [a] -> [b]
mapMaybe ((HorizontalPosition a
 -> HorizontalPosition a -> Maybe (MinorArc a))
-> (HorizontalPosition a, HorizontalPosition a)
-> Maybe (MinorArc a)
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry HorizontalPosition a -> HorizontalPosition a -> Maybe (MinorArc a)
forall a.
Spherical a =>
HorizontalPosition a -> HorizontalPosition a -> Maybe (MinorArc a)
GreatCircle.minorArc) ([HorizontalPosition a]
-> [HorizontalPosition a]
-> [(HorizontalPosition a, HorizontalPosition a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [HorizontalPosition a]
ps ([HorizontalPosition a] -> [HorizontalPosition a]
forall a. [a] -> [a]
tail [HorizontalPosition a]
ps [HorizontalPosition a]
-> [HorizontalPosition a] -> [HorizontalPosition a]
forall a. [a] -> [a] -> [a]
++ [[HorizontalPosition a] -> HorizontalPosition a
forall a. [a] -> a
head [HorizontalPosition a]
ps]))

discretiseArc :: (Spherical a) => HorizontalPosition a -> Length -> [Double] -> Polygon a
discretiseArc :: HorizontalPosition a -> Length -> [Double] -> Polygon a
discretiseArc HorizontalPosition a
c Length
r [Double]
as = [HorizontalPosition a] -> [MinorArc a] -> Bool -> Polygon a
forall a.
[HorizontalPosition a] -> [MinorArc a] -> Bool -> Polygon a
Polygon [HorizontalPosition a]
ps ([HorizontalPosition a] -> [MinorArc a]
forall a. Spherical a => [HorizontalPosition a] -> [MinorArc a]
mkEdges [HorizontalPosition a]
ps) Bool
False
  where
    m :: a
m = HorizontalPosition a -> a
forall a. Model a => HorizontalPosition a -> a
Geodetic.model HorizontalPosition a
c
    lat :: Angle
lat = HorizontalPosition a -> Angle
forall a. HasCoordinates a => a -> Angle
Geodetic.latitude HorizontalPosition a
c
    lon :: Angle
lon = HorizontalPosition a -> Angle
forall a. HasCoordinates a => a -> Angle
Geodetic.longitude HorizontalPosition a
c
    erm :: Double
erm = Length -> Double
Length.toMetres (Length -> Double) -> (a -> Length) -> a -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Ellipsoid -> Length
meanRadius (Ellipsoid -> Length) -> (a -> Ellipsoid) -> a -> Length
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> Ellipsoid
forall a. Model a => a -> Ellipsoid
surface (a -> Double) -> a -> Double
forall a b. (a -> b) -> a -> b
$ a
m
    rm :: Double
rm = Double
erm Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin (Length -> Double
Length.toMetres Length
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
erm)
    z :: Double
z = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
erm Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
erm Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
rm Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rm)
    rya :: Double
rya = Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Angle -> Double
Angle.toRadians Angle
lat
    cy :: Double
cy = Double -> Double
forall a. Floating a => a -> a
cos Double
rya
    sy :: Double
sy = Double -> Double
forall a. Floating a => a -> a
sin Double
rya
    ry :: [V3]
ry = [Double -> Double -> Double -> V3
Math3d.vec3 Double
cy Double
0 Double
sy, Double -> Double -> Double -> V3
Math3d.vec3 Double
0 Double
1 Double
0, Double -> Double -> Double -> V3
Math3d.vec3 (-Double
sy) Double
0 Double
cy]
    rza :: Double
rza = Angle -> Double
Angle.toRadians Angle
lon
    cz :: Double
cz = Double -> Double
forall a. Floating a => a -> a
cos Double
rza
    sz :: Double
sz = Double -> Double
forall a. Floating a => a -> a
sin Double
rza
    rz :: [V3]
rz = [Double -> Double -> Double -> V3
Math3d.vec3 Double
cz (-Double
sz) Double
0, Double -> Double -> Double -> V3
Math3d.vec3 Double
sz Double
cz Double
0, Double -> Double -> Double -> V3
Math3d.vec3 Double
0 Double
0 Double
1]
    anp :: [V3]
anp = (Double -> V3) -> [Double] -> [V3]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\Double
a -> Double -> Double -> Double -> V3
Math3d.vec3 ((-Double
rm) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
a) (Double
rm Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
a) Double
z) [Double]
as -- arc at north pole

    rot :: [V3]
rot = (V3 -> V3) -> [V3] -> [V3]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\V3
v -> V3 -> [V3] -> V3
Math3d.multM (V3 -> [V3] -> V3
Math3d.multM V3
v [V3]
ry) [V3]
rz) [V3]
anp -- rotate each point to arc centre

    ps :: [HorizontalPosition a]
ps = (V3 -> HorizontalPosition a) -> [V3] -> [HorizontalPosition a]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (V3 -> a -> HorizontalPosition a
forall a. Model a => V3 -> a -> HorizontalPosition a
`Geodetic.nvectorPos'` a
m) [V3]
rot