-- SPDX-License-Identifier: MPL-2.0

{- | Planar hulls, polygon measurements and Boolean arrangements. Hulls return
unique counterclockwise extreme vertices, omitting collinear interior points;
0/1/2-point hulls are valid. Boolean inputs are simple polygons with at most
128 vertices each. They use exact Rational edge splitting and local side
classification, and return oriented rings: counterclockwise exterior rings,
clockwise holes. Rings may share isolated boundary vertices. Coordinates are
rounded to Double only at the output. This is a bounded reference arrangement,
not a sweep-line implementation.
-}
module HMath.Geometry.Polygon (
    HullAlgorithm (..),
    convexHull,
    polygonArea,
    polygonCentroid,
    shoelace,
    BooleanOperation (..),
    polygonBoolean,
    clipPolygon,
    minimumBoundingCircle,
    minimumBoundingRectangle,
    minkowskiSum,
) where

import Control.Monad (foldM)
import Data.List (delete, maximumBy, minimumBy, nub, sort)
import Data.Ord (comparing)
import HMath.Error (MathError (..))
import HMath.Geometry.Internal
import HMath.Geometry.Primitives (circumcenter)

-- | Four distinct planar convex-hull construction algorithms.
data HullAlgorithm = GrahamScan | MonotoneChain | Quickhull | GiftWrapping deriving (HullAlgorithm -> HullAlgorithm -> Bool
(HullAlgorithm -> HullAlgorithm -> Bool)
-> (HullAlgorithm -> HullAlgorithm -> Bool) -> Eq HullAlgorithm
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: HullAlgorithm -> HullAlgorithm -> Bool
== :: HullAlgorithm -> HullAlgorithm -> Bool
$c/= :: HullAlgorithm -> HullAlgorithm -> Bool
/= :: HullAlgorithm -> HullAlgorithm -> Bool
Eq, Int -> HullAlgorithm -> ShowS
[HullAlgorithm] -> ShowS
HullAlgorithm -> String
(Int -> HullAlgorithm -> ShowS)
-> (HullAlgorithm -> String)
-> ([HullAlgorithm] -> ShowS)
-> Show HullAlgorithm
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> HullAlgorithm -> ShowS
showsPrec :: Int -> HullAlgorithm -> ShowS
$cshow :: HullAlgorithm -> String
show :: HullAlgorithm -> String
$cshowList :: [HullAlgorithm] -> ShowS
showList :: [HullAlgorithm] -> ShowS
Show)

-- | Convex hull with at most 4096 input points. Exact orientation breaks degeneracies.
convexHull :: HullAlgorithm -> [Point2] -> Either MathError [Point2]
convexHull :: HullAlgorithm -> [Point2] -> Either MathError [Point2]
convexHull HullAlgorithm
algorithm [Point2]
input = do
    Bool -> String -> Either MathError ()
require ([Point2] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Point2]
input Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
4096) String
"hull reference budget exceeded"
    (Point2 -> Either MathError ()) -> [Point2] -> Either MathError ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ Point2 -> Either MathError ()
check2 [Point2]
input
    let ps :: [Point2]
ps = [Point2] -> [Point2]
forall a. Ord a => [a] -> [a]
sort ([Point2] -> [Point2]
forall a. Eq a => [a] -> [a]
nub [Point2]
input)
    [Point2] -> Either MathError [Point2]
forall a b. b -> Either a b
Right ([Point2] -> Either MathError [Point2])
-> [Point2] -> Either MathError [Point2]
forall a b. (a -> b) -> a -> b
$
        if [Point2] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Point2]
ps Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
3
            then [Point2]
ps
            else [Point2] -> [Point2]
canonical ([Point2] -> [Point2]) -> [Point2] -> [Point2]
forall a b. (a -> b) -> a -> b
$ case HullAlgorithm
algorithm of
                HullAlgorithm
MonotoneChain -> [Point2] -> [Point2]
chain [Point2]
ps
                HullAlgorithm
GrahamScan -> [Point2] -> [Point2]
graham [Point2]
ps
                HullAlgorithm
Quickhull -> [Point2] -> [Point2]
quick [Point2]
ps
                HullAlgorithm
GiftWrapping -> [Point2] -> [Point2]
gift [Point2]
ps
  where
    canonical :: [Point2] -> [Point2]
canonical [Point2]
xs = let ys :: [Point2]
ys = [Point2] -> [Point2]
ccw [Point2]
xs; p :: Point2
p = [Point2] -> Point2
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Point2]
ys; ([Point2]
a, [Point2]
b) = (Point2 -> Bool) -> [Point2] -> ([Point2], [Point2])
forall a. (a -> Bool) -> [a] -> ([a], [a])
break (Point2 -> Point2 -> Bool
forall a. Eq a => a -> a -> Bool
== Point2
p) [Point2]
ys in [Point2]
b [Point2] -> [Point2] -> [Point2]
forall a. [a] -> [a] -> [a]
++ [Point2]
a
    push :: [Point2] -> Point2 -> [Point2]
push [Point2]
stack Point2
p = case [Point2]
stack of
        Point2
b : Point2
a : [Point2]
rest | Point2 -> Point2 -> Point2 -> Ordering
orient Point2
a Point2
b Point2
p Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
/= Ordering
GT -> [Point2] -> Point2 -> [Point2]
push (Point2
a Point2 -> [Point2] -> [Point2]
forall a. a -> [a] -> [a]
: [Point2]
rest) Point2
p
        [Point2]
_ -> Point2
p Point2 -> [Point2] -> [Point2]
forall a. a -> [a] -> [a]
: [Point2]
stack
    half :: [Point2] -> [Point2]
half = [Point2] -> [Point2]
forall a. [a] -> [a]
reverse ([Point2] -> [Point2])
-> ([Point2] -> [Point2]) -> [Point2] -> [Point2]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Point2] -> Point2 -> [Point2])
-> [Point2] -> [Point2] -> [Point2]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl [Point2] -> Point2 -> [Point2]
push []
    chain :: [Point2] -> [Point2]
chain [Point2]
ps = [Point2] -> [Point2]
forall a. HasCallStack => [a] -> [a]
init ([Point2] -> [Point2]
half [Point2]
ps) [Point2] -> [Point2] -> [Point2]
forall a. [a] -> [a] -> [a]
++ [Point2] -> [Point2]
forall a. HasCallStack => [a] -> [a]
init ([Point2] -> [Point2]
half ([Point2] -> [Point2]
forall a. [a] -> [a]
reverse [Point2]
ps))
    graham :: [Point2] -> [Point2]
graham [Point2]
ps =
        let pivot :: Point2
pivot = (Point2 -> Point2 -> Ordering) -> [Point2] -> Point2
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy ((Point2 -> (Double, Double)) -> Point2 -> Point2 -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (\(P2 Double
x Double
y) -> (Double
y, Double
x))) [Point2]
ps
            polar :: Point2 -> Point2 -> Ordering
polar Point2
a Point2
b = case Point2 -> Point2 -> Point2 -> Ordering
orient Point2
pivot Point2
a Point2
b of
                Ordering
GT -> Ordering
LT
                Ordering
LT -> Ordering
GT
                Ordering
EQ -> Double -> Double -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Point2 -> Point2 -> Double
dot2 (Point2 -> Point2 -> Point2
sub2 Point2
a Point2
pivot) (Point2 -> Point2 -> Point2
sub2 Point2
a Point2
pivot)) (Point2 -> Point2 -> Double
dot2 (Point2 -> Point2 -> Point2
sub2 Point2
b Point2
pivot) (Point2 -> Point2 -> Point2
sub2 Point2
b Point2
pivot))
            ordered :: [Point2]
ordered = (Point2 -> Point2 -> Ordering) -> [Point2] -> [Point2]
forall {t :: * -> *} {t}.
Foldable t =>
(t -> t -> Ordering) -> t t -> [t]
sortByLocal Point2 -> Point2 -> Ordering
polar ((Point2 -> Bool) -> [Point2] -> [Point2]
forall a. (a -> Bool) -> [a] -> [a]
filter (Point2 -> Point2 -> Bool
forall a. Eq a => a -> a -> Bool
/= Point2
pivot) [Point2]
ps)
         in [Point2] -> [Point2]
forall a. [a] -> [a]
reverse (([Point2] -> Point2 -> [Point2])
-> [Point2] -> [Point2] -> [Point2]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl [Point2] -> Point2 -> [Point2]
push [Point2
pivot] [Point2]
ordered)
    quick :: [Point2] -> [Point2]
quick [Point2]
ps =
        let a :: Point2
a = Point2 -> [Point2] -> Point2
forall a. a -> [a] -> a
firstOr (Double -> Double -> Point2
P2 Double
0 Double
0) [Point2]
ps
            b :: Point2
b = [Point2] -> Point2
forall a. HasCallStack => [a] -> a
last [Point2]
ps
            side :: Point2 -> Point2 -> [Point2] -> [Point2]
side Point2
u Point2
v [Point2]
qs = case (Point2 -> Bool) -> [Point2] -> [Point2]
forall a. (a -> Bool) -> [a] -> [a]
filter (\Point2
p -> Point2 -> Point2 -> Point2 -> Ordering
orient Point2
u Point2
v Point2
p Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
GT) [Point2]
qs of
                [] -> []
                [Point2]
inside ->
                    let p :: Point2
p = (Point2 -> Point2 -> Ordering) -> [Point2] -> Point2
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
maximumBy ((Point2 -> Rational) -> Point2 -> Point2 -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (\Point2
q -> Point2 -> Point2 -> Point2 -> Rational
orientExact Point2
u Point2
v Point2
q)) [Point2]
inside
                     in Point2 -> Point2 -> [Point2] -> [Point2]
side Point2
u Point2
p [Point2]
inside [Point2] -> [Point2] -> [Point2]
forall a. [a] -> [a] -> [a]
++ [Point2
p] [Point2] -> [Point2] -> [Point2]
forall a. [a] -> [a] -> [a]
++ Point2 -> Point2 -> [Point2] -> [Point2]
side Point2
p Point2
v [Point2]
inside
         in [Point2
a] [Point2] -> [Point2] -> [Point2]
forall a. [a] -> [a] -> [a]
++ Point2 -> Point2 -> [Point2] -> [Point2]
side Point2
a Point2
b [Point2]
ps [Point2] -> [Point2] -> [Point2]
forall a. [a] -> [a] -> [a]
++ [Point2
b] [Point2] -> [Point2] -> [Point2]
forall a. [a] -> [a] -> [a]
++ Point2 -> Point2 -> [Point2] -> [Point2]
side Point2
b Point2
a [Point2]
ps
    gift :: [Point2] -> [Point2]
gift [Point2]
ps =
        let start :: Point2
start = Point2 -> [Point2] -> Point2
forall a. a -> [a] -> a
firstOr (Double -> Double -> Point2
P2 Double
0 Double
0) [Point2]
ps
            walk :: Point2 -> [Point2] -> [Point2]
walk Point2
p [Point2]
acc =
                let candidates :: [Point2]
candidates = (Point2 -> Bool) -> [Point2] -> [Point2]
forall a. (a -> Bool) -> [a] -> [a]
filter (Point2 -> Point2 -> Bool
forall a. Eq a => a -> a -> Bool
/= Point2
p) [Point2]
ps
                    better :: Point2 -> Point2 -> Point2
better Point2
a Point2
b = case Point2 -> Point2 -> Point2 -> Ordering
orient Point2
p Point2
a Point2
b of Ordering
LT -> Point2
b; Ordering
GT -> Point2
a; Ordering
EQ -> if Point2 -> Double
norm2 (Point2 -> Point2 -> Point2
sub2 Point2
b Point2
p) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Point2 -> Double
norm2 (Point2 -> Point2 -> Point2
sub2 Point2
a Point2
p) then Point2
b else Point2
a
                    next :: Point2
next = (Point2 -> Point2 -> Point2) -> [Point2] -> Point2
forall a. (a -> a -> a) -> [a] -> a
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 Point2 -> Point2 -> Point2
better [Point2]
candidates
                 in if Point2
next Point2 -> Point2 -> Bool
forall a. Eq a => a -> a -> Bool
== Point2
start then [Point2]
acc [Point2] -> [Point2] -> [Point2]
forall a. [a] -> [a] -> [a]
++ [Point2
p] else Point2 -> [Point2] -> [Point2]
walk Point2
next ([Point2]
acc [Point2] -> [Point2] -> [Point2]
forall a. [a] -> [a] -> [a]
++ [Point2
p])
         in Point2 -> [Point2] -> [Point2]
walk Point2
start []
    sortByLocal :: (t -> t -> Ordering) -> t t -> [t]
sortByLocal t -> t -> Ordering
cmp = (t -> [t] -> [t]) -> [t] -> t t -> [t]
forall a b. (a -> b -> b) -> b -> t a -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\t
x [t]
xs -> let ([t]
a, [t]
b) = (t -> Bool) -> [t] -> ([t], [t])
forall a. (a -> Bool) -> [a] -> ([a], [a])
span (\t
y -> t -> t -> Ordering
cmp t
x t
y Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
GT) [t]
xs in [t]
a [t] -> [t] -> [t]
forall a. [a] -> [a] -> [a]
++ [t
x] [t] -> [t] -> [t]
forall a. [a] -> [a] -> [a]
++ [t]
b) []

-- | Signed shoelace area; clockwise polygons have negative area.
shoelace :: Polygon -> Either MathError Double
shoelace :: Polygon -> Either MathError Double
shoelace Polygon
poly = [Point2] -> Double
signedArea ([Point2] -> Double)
-> Either MathError [Point2] -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Polygon -> Either MathError [Point2]
checkPolygon Polygon
poly

-- | Unsigned simple polygon area.
polygonArea :: Polygon -> Either MathError Double
polygonArea :: Polygon -> Either MathError Double
polygonArea = (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall a b. (a -> b) -> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Double -> Double
forall a. Num a => a -> a
abs (Either MathError Double -> Either MathError Double)
-> (Polygon -> Either MathError Double)
-> Polygon
-> Either MathError Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Polygon -> Either MathError Double
shoelace

-- | Area-weighted centroid; works for nonconvex simple polygons.
polygonCentroid :: Polygon -> Either MathError Point2
polygonCentroid :: Polygon -> Either MathError Point2
polygonCentroid Polygon
poly = do
    ps <- Polygon -> Either MathError [Point2]
checkPolygon Polygon
poly
    let rational (P2 Double
x Double
y) = (Double -> Rational
forall a. Real a => a -> Rational
toRational Double
x, Double -> Rational
forall a. Real a => a -> Rational
toRational Double
y)
        es = [(Rational, Rational)]
-> [((Rational, Rational), (Rational, Rational))]
forall a. [a] -> [(a, a)]
edges ((Point2 -> (Rational, Rational))
-> [Point2] -> [(Rational, Rational)]
forall a b. (a -> b) -> [a] -> [b]
map Point2 -> (Rational, Rational)
rational [Point2]
ps)
        terms = [(Rational
ax Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
by Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
ay Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
bx, Rational
ax Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
bx, Rational
ay Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
by) | ((Rational
ax, Rational
ay), (Rational
bx, Rational
by)) <- [((Rational, Rational), (Rational, Rational))]
es]
        area6 = Rational
3 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* [Rational] -> Rational
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Rational
c | (Rational
c, Rational
_, Rational
_) <- [(Rational, Rational, Rational)]
terms]
    point2
        (fromRational (sum [c * x | (c, x, _) <- terms] / area6))
        (fromRational (sum [c * y | (c, _, y) <- terms] / area6))

-- | Boolean set operation on closed polygon regions; ring boundaries have zero area.
data BooleanOperation = Union | Intersection | Difference deriving (BooleanOperation -> BooleanOperation -> Bool
(BooleanOperation -> BooleanOperation -> Bool)
-> (BooleanOperation -> BooleanOperation -> Bool)
-> Eq BooleanOperation
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: BooleanOperation -> BooleanOperation -> Bool
== :: BooleanOperation -> BooleanOperation -> Bool
$c/= :: BooleanOperation -> BooleanOperation -> Bool
/= :: BooleanOperation -> BooleanOperation -> Bool
Eq, Int -> BooleanOperation -> ShowS
[BooleanOperation] -> ShowS
BooleanOperation -> String
(Int -> BooleanOperation -> ShowS)
-> (BooleanOperation -> String)
-> ([BooleanOperation] -> ShowS)
-> Show BooleanOperation
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> BooleanOperation -> ShowS
showsPrec :: Int -> BooleanOperation -> ShowS
$cshow :: BooleanOperation -> String
show :: BooleanOperation -> String
$cshowList :: [BooleanOperation] -> ShowS
showList :: [BooleanOperation] -> ShowS
Show)

type RPoint = (Rational, Rational)
rp :: Point2 -> RPoint
rp :: Point2 -> (Rational, Rational)
rp (P2 Double
x Double
y) = (Double -> Rational
forall a. Real a => a -> Rational
toRational Double
x, Double -> Rational
forall a. Real a => a -> Rational
toRational Double
y)
radd :: RPoint -> RPoint -> RPoint
radd :: (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
radd (Rational
a, Rational
b) (Rational
x, Rational
y) = (Rational
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
x, Rational
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
y)
rsub :: RPoint -> RPoint -> RPoint
rsub :: (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational
a, Rational
b) (Rational
x, Rational
y) = (Rational
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
x, Rational
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
y)
rscale :: Rational -> RPoint -> RPoint
rscale :: Rational -> (Rational, Rational) -> (Rational, Rational)
rscale Rational
t (Rational
x, Rational
y) = (Rational
t Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
x, Rational
t Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
y)
rcross :: RPoint -> RPoint -> Rational
rcross :: (Rational, Rational) -> (Rational, Rational) -> Rational
rcross (Rational
a, Rational
b) (Rational
x, Rational
y) = Rational
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
y Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
x
rdot :: RPoint -> RPoint -> Rational
rdot :: (Rational, Rational) -> (Rational, Rational) -> Rational
rdot (Rational
a, Rational
b) (Rational
x, Rational
y) = Rational
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
x Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
y
rOn :: RPoint -> RPoint -> RPoint -> Bool
rOn :: (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational) -> Bool
rOn (Rational, Rational)
a (Rational, Rational)
b (Rational, Rational)
p = (Rational, Rational) -> (Rational, Rational) -> Rational
rcross ((Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
b (Rational, Rational)
a) ((Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
p (Rational, Rational)
a) Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0 Bool -> Bool -> Bool
&& (Rational, Rational) -> (Rational, Rational) -> Rational
rdot ((Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
p (Rational, Rational)
a) ((Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
p (Rational, Rational)
b) Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Rational
0
rInside :: [RPoint] -> RPoint -> Bool
rInside :: [(Rational, Rational)] -> (Rational, Rational) -> Bool
rInside [(Rational, Rational)]
ps (Rational
x, Rational
y) =
    Int -> Bool
forall a. Integral a => a -> Bool
odd
        ( [()] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length
            [ ()
            | ((Rational
ax, Rational
ay), (Rational
bx, Rational
by)) <- [(Rational, Rational)]
-> [((Rational, Rational), (Rational, Rational))]
forall a. [a] -> [(a, a)]
edges [(Rational, Rational)]
ps
            , (Rational
ay Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
y) Bool -> Bool -> Bool
forall a. Eq a => a -> a -> Bool
/= (Rational
by Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
y)
            , Rational
x Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< (Rational
bx Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
ax) Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* (Rational
y Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
ay) Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ (Rational
by Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
ay) Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
ax
            ]
        )

{- | Exact edge-arrangement Boolean. Empty output is allowed, including an
intersection that consists only of points/segments (zero-area components).
-}
polygonBoolean :: BooleanOperation -> Polygon -> Polygon -> Either MathError [Polygon]
polygonBoolean :: BooleanOperation
-> Polygon -> Polygon -> Either MathError [Polygon]
polygonBoolean BooleanOperation
operation Polygon
a Polygon
b = do
    pa <- Polygon -> Either MathError [Point2]
checkPolygon Polygon
a
    pb <- checkPolygon b
    require (length pa <= 128 && length pb <= 128) "polygon Boolean budget is 128 vertices per input"
    let aa = (Point2 -> (Rational, Rational))
-> [Point2] -> [(Rational, Rational)]
forall a b. (a -> b) -> [a] -> [b]
map Point2 -> (Rational, Rational)
rp ([Point2] -> [Point2]
ccw [Point2]
pa)
        bb = (Point2 -> (Rational, Rational))
-> [Point2] -> [(Rational, Rational)]
forall a b. (a -> b) -> [a] -> [b]
map Point2 -> (Rational, Rational)
rp ([Point2] -> [Point2]
ccw [Point2]
pb)
        allEdges = [(Rational, Rational)]
-> [((Rational, Rational), (Rational, Rational))]
forall a. [a] -> [(a, a)]
edges [(Rational, Rational)]
aa [((Rational, Rational), (Rational, Rational))]
-> [((Rational, Rational), (Rational, Rational))]
-> [((Rational, Rational), (Rational, Rational))]
forall a. [a] -> [a] -> [a]
++ [(Rational, Rational)]
-> [((Rational, Rational), (Rational, Rational))]
forall a. [a] -> [(a, a)]
edges [(Rational, Rational)]
bb
        splitEdge ((Rational, Rational)
u, (Rational, Rational)
v) =
            let dir :: (Rational, Rational)
dir = (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
v (Rational, Rational)
u
                denom :: Rational
denom = (Rational, Rational) -> (Rational, Rational) -> Rational
rdot (Rational, Rational)
dir (Rational, Rational)
dir
                params :: [Rational]
params = [(Rational, Rational) -> (Rational, Rational) -> Rational
rdot ((Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
p (Rational, Rational)
u) (Rational, Rational)
dir Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
denom | ((Rational, Rational)
c, (Rational, Rational)
d) <- [((Rational, Rational), (Rational, Rational))]
allEdges, (Rational, Rational)
p <- ((Rational, Rational), (Rational, Rational))
-> ((Rational, Rational), (Rational, Rational))
-> [(Rational, Rational)]
intersections ((Rational, Rational)
u, (Rational, Rational)
v) ((Rational, Rational)
c, (Rational, Rational)
d), (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational) -> Bool
rOn (Rational, Rational)
u (Rational, Rational)
v (Rational, Rational)
p]
                ts :: [Rational]
ts = [Rational] -> [Rational]
forall a. Ord a => [a] -> [a]
sort ([Rational] -> [Rational]
forall a. Eq a => [a] -> [a]
nub (Rational
0 Rational -> [Rational] -> [Rational]
forall a. a -> [a] -> [a]
: Rational
1 Rational -> [Rational] -> [Rational]
forall a. a -> [a] -> [a]
: [Rational]
params))
                at :: Rational -> (Rational, Rational)
at Rational
t = (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
radd (Rational, Rational)
u (Rational -> (Rational, Rational) -> (Rational, Rational)
rscale Rational
t (Rational, Rational)
dir)
             in [(Rational, Rational)]
-> [(Rational, Rational)]
-> [((Rational, Rational), (Rational, Rational))]
forall a b. [a] -> [b] -> [(a, b)]
zip ((Rational -> (Rational, Rational))
-> [Rational] -> [(Rational, Rational)]
forall a b. (a -> b) -> [a] -> [b]
map Rational -> (Rational, Rational)
at ([Rational] -> [Rational]
forall a. HasCallStack => [a] -> [a]
init [Rational]
ts)) ((Rational -> (Rational, Rational))
-> [Rational] -> [(Rational, Rational)]
forall a b. (a -> b) -> [a] -> [b]
map Rational -> (Rational, Rational)
at (Int -> [Rational] -> [Rational]
forall a. Int -> [a] -> [a]
drop Int
1 [Rational]
ts))
        pieces = [((Rational, Rational), (Rational, Rational))]
-> [((Rational, Rational), (Rational, Rational))]
forall a. Eq a => [a] -> [a]
nub [if (Rational, Rational)
u (Rational, Rational) -> (Rational, Rational) -> Bool
forall a. Ord a => a -> a -> Bool
< (Rational, Rational)
v then ((Rational, Rational)
u, (Rational, Rational)
v) else ((Rational, Rational)
v, (Rational, Rational)
u) | ((Rational, Rational), (Rational, Rational))
e <- [((Rational, Rational), (Rational, Rational))]
allEdges, ((Rational, Rational)
u, (Rational, Rational)
v) <- ((Rational, Rational), (Rational, Rational))
-> [((Rational, Rational), (Rational, Rational))]
splitEdge ((Rational, Rational), (Rational, Rational))
e, (Rational, Rational)
u (Rational, Rational) -> (Rational, Rational) -> Bool
forall a. Eq a => a -> a -> Bool
/= (Rational, Rational)
v]
        truth Bool
x Bool
y = case BooleanOperation
operation of BooleanOperation
Union -> Bool
x Bool -> Bool -> Bool
|| Bool
y; BooleanOperation
Intersection -> Bool
x Bool -> Bool -> Bool
&& Bool
y; BooleanOperation
Difference -> Bool
x Bool -> Bool -> Bool
&& Bool -> Bool
not Bool
y
        boundary ((Rational, Rational)
u, (Rational, Rational)
v) =
            let midpoint :: (Rational, Rational)
midpoint = Rational -> (Rational, Rational) -> (Rational, Rational)
rscale (Rational
1 Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
2) ((Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
radd (Rational, Rational)
u (Rational, Rational)
v)
                (Rational
dx, Rational
dy) = (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
v (Rational, Rational)
u
                normal :: (Rational, Rational)
normal = (-Rational
dy, Rational
dx)
                cuts :: [Rational]
cuts =
                    [ Rational -> Rational
forall a. Num a => a -> a
abs Rational
t
                    | ((Rational, Rational)
c, (Rational, Rational)
d) <- [((Rational, Rational), (Rational, Rational))]
allEdges
                    , let edge :: (Rational, Rational)
edge = (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
d (Rational, Rational)
c
                    , let den :: Rational
den = (Rational, Rational) -> (Rational, Rational) -> Rational
rcross (Rational, Rational)
normal (Rational, Rational)
edge
                    , Rational
den Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
0
                    , let t :: Rational
t = (Rational, Rational) -> (Rational, Rational) -> Rational
rcross ((Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
c (Rational, Rational)
midpoint) (Rational, Rational)
edge Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
den
                    , let s :: Rational
s = (Rational, Rational) -> (Rational, Rational) -> Rational
rcross ((Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
c (Rational, Rational)
midpoint) (Rational, Rational)
normal Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
den
                    , Rational
t Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
0
                    , Rational
s Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
>= Rational
0
                    , Rational
s Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Rational
1
                    ]
                epsilon :: Rational
epsilon = [Rational] -> Rational
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum (Rational
1 Rational -> [Rational] -> [Rational]
forall a. a -> [a] -> [a]
: [Rational]
cuts) Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
2
                l :: (Rational, Rational)
l = (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
radd (Rational, Rational)
midpoint (Rational -> (Rational, Rational) -> (Rational, Rational)
rscale Rational
epsilon (Rational, Rational)
normal)
                r :: (Rational, Rational)
r = (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
midpoint (Rational -> (Rational, Rational) -> (Rational, Rational)
rscale Rational
epsilon (Rational, Rational)
normal)
                left :: Bool
left = Bool -> Bool -> Bool
truth ([(Rational, Rational)] -> (Rational, Rational) -> Bool
rInside [(Rational, Rational)]
aa (Rational, Rational)
l) ([(Rational, Rational)] -> (Rational, Rational) -> Bool
rInside [(Rational, Rational)]
bb (Rational, Rational)
l)
                right :: Bool
right = Bool -> Bool -> Bool
truth ([(Rational, Rational)] -> (Rational, Rational) -> Bool
rInside [(Rational, Rational)]
aa (Rational, Rational)
r) ([(Rational, Rational)] -> (Rational, Rational) -> Bool
rInside [(Rational, Rational)]
bb (Rational, Rational)
r)
             in if Bool
left Bool -> Bool -> Bool
forall a. Eq a => a -> a -> Bool
== Bool
right then [] else if Bool
left then [((Rational, Rational)
u, (Rational, Rational)
v)] else [((Rational, Rational)
v, (Rational, Rational)
u)]
        segments = (((Rational, Rational), (Rational, Rational))
 -> [((Rational, Rational), (Rational, Rational))])
-> [((Rational, Rational), (Rational, Rational))]
-> [((Rational, Rational), (Rational, Rational))]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap ((Rational, Rational), (Rational, Rational))
-> [((Rational, Rational), (Rational, Rational))]
boundary [((Rational, Rational), (Rational, Rational))]
pieces
    loops <- stitch segments
    mapM (\[(Rational, Rational)]
ring -> [Point2] -> Polygon
Polygon ([Point2] -> Polygon)
-> Either MathError [Point2] -> Either MathError Polygon
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> ((Rational, Rational) -> Either MathError Point2)
-> [(Rational, Rational)] -> Either MathError [Point2]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM (\(Rational
x, Rational
y) -> Double -> Double -> Either MathError Point2
point2 (Rational -> Double
forall a. Fractional a => Rational -> a
fromRational Rational
x) (Rational -> Double
forall a. Fractional a => Rational -> a
fromRational Rational
y)) [(Rational, Rational)]
ring) loops
  where
    intersections :: ((Rational, Rational), (Rational, Rational))
-> ((Rational, Rational), (Rational, Rational))
-> [(Rational, Rational)]
intersections ((Rational, Rational)
a0, (Rational, Rational)
b0) ((Rational, Rational)
c0, (Rational, Rational)
d0) =
        let u :: (Rational, Rational)
u = (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
b0 (Rational, Rational)
a0; v :: (Rational, Rational)
v = (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
d0 (Rational, Rational)
c0; den :: Rational
den = (Rational, Rational) -> (Rational, Rational) -> Rational
rcross (Rational, Rational)
u (Rational, Rational)
v; w :: (Rational, Rational)
w = (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
c0 (Rational, Rational)
a0
         in if Rational
den Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0
                then if (Rational, Rational) -> (Rational, Rational) -> Rational
rcross (Rational, Rational)
w (Rational, Rational)
u Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0 then ((Rational, Rational) -> Bool)
-> [(Rational, Rational)] -> [(Rational, Rational)]
forall a. (a -> Bool) -> [a] -> [a]
filter (\(Rational, Rational)
p -> (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational) -> Bool
rOn (Rational, Rational)
a0 (Rational, Rational)
b0 (Rational, Rational)
p Bool -> Bool -> Bool
&& (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational) -> Bool
rOn (Rational, Rational)
c0 (Rational, Rational)
d0 (Rational, Rational)
p) [(Rational, Rational)
a0, (Rational, Rational)
b0, (Rational, Rational)
c0, (Rational, Rational)
d0] else []
                else
                    let t :: Rational
t = (Rational, Rational) -> (Rational, Rational) -> Rational
rcross (Rational, Rational)
w (Rational, Rational)
v Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
den; s :: Rational
s = (Rational, Rational) -> (Rational, Rational) -> Rational
rcross (Rational, Rational)
w (Rational, Rational)
u Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
den
                     in if Rational
t Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
>= Rational
0 Bool -> Bool -> Bool
&& Rational
t Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Rational
1 Bool -> Bool -> Bool
&& Rational
s Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
>= Rational
0 Bool -> Bool -> Bool
&& Rational
s Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Rational
1 then [(Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
radd (Rational, Rational)
a0 (Rational -> (Rational, Rational) -> (Rational, Rational)
rscale Rational
t (Rational, Rational)
u)] else []
    stitch :: [((Rational, Rational), (Rational, Rational))]
-> Either MathError [[(Rational, Rational)]]
stitch [] = [[(Rational, Rational)]]
-> Either MathError [[(Rational, Rational)]]
forall a b. b -> Either a b
Right []
    stitch (((Rational, Rational)
start, (Rational, Rational)
next) : [((Rational, Rational), (Rational, Rational))]
remaining) = do
        (ring, rest) <- (Rational, Rational)
-> (Rational, Rational)
-> (Rational, Rational)
-> [(Rational, Rational)]
-> [((Rational, Rational), (Rational, Rational))]
-> Int
-> Either
     MathError
     ([(Rational, Rational)],
      [((Rational, Rational), (Rational, Rational))])
forall {t}.
(Ord t, Num t) =>
(Rational, Rational)
-> (Rational, Rational)
-> (Rational, Rational)
-> [(Rational, Rational)]
-> [((Rational, Rational), (Rational, Rational))]
-> t
-> Either
     MathError
     ([(Rational, Rational)],
      [((Rational, Rational), (Rational, Rational))])
walk (Rational, Rational)
start (Rational, Rational)
start (Rational, Rational)
next [(Rational, Rational)
start] [((Rational, Rational), (Rational, Rational))]
remaining ([((Rational, Rational), (Rational, Rational))] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [((Rational, Rational), (Rational, Rational))]
remaining Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)
        others <- stitch rest
        Right (if length ring >= 3 then ring : others else others)
    walk :: (Rational, Rational)
-> (Rational, Rational)
-> (Rational, Rational)
-> [(Rational, Rational)]
-> [((Rational, Rational), (Rational, Rational))]
-> t
-> Either
     MathError
     ([(Rational, Rational)],
      [((Rational, Rational), (Rational, Rational))])
walk (Rational, Rational)
start (Rational, Rational)
before (Rational, Rational)
current [(Rational, Rational)]
acc [((Rational, Rational), (Rational, Rational))]
remaining t
budget
        | (Rational, Rational)
current (Rational, Rational) -> (Rational, Rational) -> Bool
forall a. Eq a => a -> a -> Bool
== (Rational, Rational)
start = ([(Rational, Rational)],
 [((Rational, Rational), (Rational, Rational))])
-> Either
     MathError
     ([(Rational, Rational)],
      [((Rational, Rational), (Rational, Rational))])
forall a b. b -> Either a b
Right ([(Rational, Rational)]
acc, [((Rational, Rational), (Rational, Rational))]
remaining)
        | t
budget t -> t -> Bool
forall a. Ord a => a -> a -> Bool
<= t
0 = MathError
-> Either
     MathError
     ([(Rational, Rational)],
      [((Rational, Rational), (Rational, Rational))])
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"polygon ring stitching" Int
0)
        | Bool
otherwise = case (((Rational, Rational), (Rational, Rational)) -> Bool)
-> [((Rational, Rational), (Rational, Rational))]
-> [((Rational, Rational), (Rational, Rational))]
forall a. (a -> Bool) -> [a] -> [a]
filter (((Rational, Rational) -> (Rational, Rational) -> Bool
forall a. Eq a => a -> a -> Bool
== (Rational, Rational)
current) ((Rational, Rational) -> Bool)
-> (((Rational, Rational), (Rational, Rational))
    -> (Rational, Rational))
-> ((Rational, Rational), (Rational, Rational))
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Rational, Rational), (Rational, Rational))
-> (Rational, Rational)
forall a b. (a, b) -> a
fst) [((Rational, Rational), (Rational, Rational))]
remaining of
            [] -> MathError
-> Either
     MathError
     ([(Rational, Rational)],
      [((Rational, Rational), (Rational, Rational))])
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"polygon arrangement has an open boundary")
            [((Rational, Rational), (Rational, Rational))]
candidates ->
                let incoming :: (Rational, Rational)
incoming = (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
current (Rational, Rational)
before
                    turn :: (a, (Rational, Rational)) -> Double
turn (a
_, (Rational, Rational)
p) =
                        let (Rational
x, Rational
y) = (Rational, Rational)
-> (Rational, Rational) -> (Rational, Rational)
rsub (Rational, Rational)
p (Rational, Rational)
current; (Rational
u, Rational
v) = (Rational, Rational)
incoming
                         in Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Rational -> Double
forall a. Fractional a => Rational -> a
fromRational (Rational
u Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
y Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
v Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
x) :: Double) (Rational -> Double
forall a. Fractional a => Rational -> a
fromRational (Rational
u Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
x Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
v Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
y))
                    chosen :: ((Rational, Rational), (Rational, Rational))
chosen = (((Rational, Rational), (Rational, Rational))
 -> ((Rational, Rational), (Rational, Rational)) -> Ordering)
-> [((Rational, Rational), (Rational, Rational))]
-> ((Rational, Rational), (Rational, Rational))
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
maximumBy ((((Rational, Rational), (Rational, Rational)) -> Double)
-> ((Rational, Rational), (Rational, Rational))
-> ((Rational, Rational), (Rational, Rational))
-> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing ((Rational, Rational), (Rational, Rational)) -> Double
forall {a}. (a, (Rational, Rational)) -> Double
turn) [((Rational, Rational), (Rational, Rational))]
candidates
                 in (Rational, Rational)
-> (Rational, Rational)
-> (Rational, Rational)
-> [(Rational, Rational)]
-> [((Rational, Rational), (Rational, Rational))]
-> t
-> Either
     MathError
     ([(Rational, Rational)],
      [((Rational, Rational), (Rational, Rational))])
walk (Rational, Rational)
start (Rational, Rational)
current (((Rational, Rational), (Rational, Rational))
-> (Rational, Rational)
forall a b. (a, b) -> b
snd ((Rational, Rational), (Rational, Rational))
chosen) ([(Rational, Rational)]
acc [(Rational, Rational)]
-> [(Rational, Rational)] -> [(Rational, Rational)]
forall a. [a] -> [a] -> [a]
++ [(Rational, Rational)
current]) (((Rational, Rational), (Rational, Rational))
-> [((Rational, Rational), (Rational, Rational))]
-> [((Rational, Rational), (Rational, Rational))]
forall a. Eq a => a -> [a] -> [a]
delete ((Rational, Rational), (Rational, Rational))
chosen [((Rational, Rational), (Rational, Rational))]
remaining) (t
budget t -> t -> t
forall a. Num a => a -> a -> a
- t
1)

{- | Sutherland-Hodgman clipping against a convex CCW/clockwise clip polygon.
Output may be empty; nonconvex clippers fail instead of silently corrupting output.
-}
clipPolygon :: Polygon -> Polygon -> Either MathError [Point2]
clipPolygon :: Polygon -> Polygon -> Either MathError [Point2]
clipPolygon Polygon
subject Polygon
clipper = do
    ps <- Polygon -> Either MathError [Point2]
checkPolygon Polygon
subject
    cs <- ccw <$> checkPolygon clipper
    require (and [orient a b c /= LT | (a, b, c) <- triples cs]) "clip polygon is not convex"
    foldM clip ps (edges cs)
  where
    triples :: [c] -> [(c, c, c)]
triples [c]
xs = [c] -> [c] -> [c] -> [(c, c, c)]
forall a b c. [a] -> [b] -> [c] -> [(a, b, c)]
zip3 [c]
xs (Int -> [c] -> [c]
forall a. Int -> [a] -> [a]
drop Int
1 [c]
xs [c] -> [c] -> [c]
forall a. [a] -> [a] -> [a]
++ Int -> [c] -> [c]
forall a. Int -> [a] -> [a]
take Int
1 [c]
xs) (Int -> [c] -> [c]
forall a. Int -> [a] -> [a]
drop Int
2 [c]
xs [c] -> [c] -> [c]
forall a. [a] -> [a] -> [a]
++ Int -> [c] -> [c]
forall a. Int -> [a] -> [a]
take Int
2 [c]
xs)
    clip :: [Point2] -> (Point2, Point2) -> Either MathError [Point2]
clip [] (Point2, Point2)
_ = [Point2] -> Either MathError [Point2]
forall a b. b -> Either a b
Right []
    clip [Point2]
input (Point2
a, Point2
b) =
        ([[Point2]] -> [Point2])
-> Either MathError [[Point2]] -> Either MathError [Point2]
forall a b. (a -> b) -> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap [[Point2]] -> [Point2]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat (Either MathError [[Point2]] -> Either MathError [Point2])
-> Either MathError [[Point2]] -> Either MathError [Point2]
forall a b. (a -> b) -> a -> b
$
            ((Point2, Point2) -> Either MathError [Point2])
-> [(Point2, Point2)] -> Either MathError [[Point2]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM
                ( \(Point2
u, Point2
v) ->
                    let iu :: Bool
iu = Point2 -> Point2 -> Point2 -> Ordering
orient Point2
a Point2
b Point2
u Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
/= Ordering
LT
                        iv :: Bool
iv = Point2 -> Point2 -> Point2 -> Ordering
orient Point2
a Point2
b Point2
v Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
/= Ordering
LT
                        hit :: Either MathError [Point2]
hit = do
                            result <- Point2
-> Point2 -> Point2 -> Point2 -> Either MathError (Maybe Point2)
lineIntersectionLocal Point2
u (Point2 -> Point2 -> Point2
sub2 Point2
v Point2
u) Point2
a (Point2 -> Point2 -> Point2
sub2 Point2
b Point2
a)
                            case result of { Just Point2
p -> [Point2] -> Either MathError [Point2]
forall a b. b -> Either a b
Right [Point2
p]; Maybe Point2
Nothing -> [Point2] -> Either MathError [Point2]
forall a b. b -> Either a b
Right [] }
                     in case (Bool
iu, Bool
iv) of
                            (Bool
True, Bool
True) -> [Point2] -> Either MathError [Point2]
forall a b. b -> Either a b
Right [Point2
v]
                            (Bool
True, Bool
False) -> Either MathError [Point2]
hit
                            (Bool
False, Bool
True) -> ([Point2] -> [Point2] -> [Point2]
forall a. [a] -> [a] -> [a]
++ [Point2
v]) ([Point2] -> [Point2])
-> Either MathError [Point2] -> Either MathError [Point2]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Either MathError [Point2]
hit
                            (Bool, Bool)
_ -> [Point2] -> Either MathError [Point2]
forall a b. b -> Either a b
Right []
                )
                ([Point2] -> [(Point2, Point2)]
forall a. [a] -> [(a, a)]
edges [Point2]
input)
    lineIntersectionLocal :: Point2
-> Point2 -> Point2 -> Point2 -> Either MathError (Maybe Point2)
lineIntersectionLocal Point2
a Point2
u Point2
b Point2
v =
        let den :: Rational
den = Point2 -> Point2 -> Point2 -> Rational
orientExact (Double -> Double -> Point2
P2 Double
0 Double
0) Point2
u Point2
v
         in if Rational
den Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0
                then Maybe Point2 -> Either MathError (Maybe Point2)
forall a b. b -> Either a b
Right Maybe Point2
forall a. Maybe a
Nothing
                else
                    let t :: Double
t = Rational -> Double
forall a. Fractional a => Rational -> a
fromRational (Point2 -> Point2 -> Point2 -> Rational
orientExact (Double -> Double -> Point2
P2 Double
0 Double
0) (Point2 -> Point2 -> Point2
sub2 Point2
b Point2
a) Point2
v Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
den); P2 Double
x Double
y = Point2 -> Point2 -> Point2
add2 Point2
a (Double -> Point2 -> Point2
scale2 Double
t Point2
u)
                     in Point2 -> Maybe Point2
forall a. a -> Maybe a
Just (Point2 -> Maybe Point2)
-> Either MathError Point2 -> Either MathError (Maybe Point2)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Double -> Double -> Either MathError Point2
point2 Double
x Double
y

{- | Deterministic incremental minimum enclosing circle, O(n^3), at most 256
points. Each violated boundary introduces a 1-, 2-, or 3-point support.
-}
minimumBoundingCircle :: [Point2] -> Either MathError Circle
minimumBoundingCircle :: [Point2] -> Either MathError Circle
minimumBoundingCircle [Point2]
ps = do
    Bool -> String -> Either MathError ()
require (Bool -> Bool
not ([Point2] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Point2]
ps) Bool -> Bool -> Bool
&& [Point2] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Point2]
ps Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
256) String
"minimum circle requires 1..256 points"
    (Point2 -> Either MathError ()) -> [Point2] -> Either MathError ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ Point2 -> Either MathError ()
check2 [Point2]
ps
    (Circle -> (Int, Point2) -> Either MathError Circle)
-> Circle -> [(Int, Point2)] -> Either MathError Circle
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM Circle -> (Int, Point2) -> Either MathError Circle
outer (Point2 -> Double -> Circle
Circle (Point2 -> [Point2] -> Point2
forall a. a -> [a] -> a
firstOr (Double -> Double -> Point2
P2 Double
0 Double
0) [Point2]
ps) Double
0) ([Int] -> [Point2] -> [(Int, Point2)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [Point2]
ps)
  where
    contains :: Circle -> Point2 -> Bool
contains (Circle Point2
c Double
r) Point2
p = Point2 -> Double
norm2 (Point2 -> Point2 -> Point2
sub2 Point2
p Point2
c) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1e-12)
    diameter :: Point2 -> Point2 -> Circle
diameter Point2
a Point2
b = Point2 -> Double -> Circle
Circle (Double -> Point2 -> Point2
scale2 Double
0.5 (Point2 -> Point2 -> Point2
add2 Point2
a Point2
b)) (Point2 -> Double
norm2 (Point2 -> Point2 -> Point2
sub2 Point2
a Point2
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
    outer :: Circle -> (Int, Point2) -> Either MathError Circle
outer Circle
circle (Int
i, Point2
p)
        | Circle -> Point2 -> Bool
contains Circle
circle Point2
p = Circle -> Either MathError Circle
forall a b. b -> Either a b
Right Circle
circle
        | Bool
otherwise = (Circle -> (Int, Point2) -> Either MathError Circle)
-> Circle -> [(Int, Point2)] -> Either MathError Circle
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM (Point2 -> Circle -> (Int, Point2) -> Either MathError Circle
middle Point2
p) (Point2 -> Double -> Circle
Circle Point2
p Double
0) ([Int] -> [Point2] -> [(Int, Point2)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] (Int -> [Point2] -> [Point2]
forall a. Int -> [a] -> [a]
take Int
i [Point2]
ps))
    middle :: Point2 -> Circle -> (Int, Point2) -> Either MathError Circle
middle Point2
p Circle
circle (Int
j, Point2
q)
        | Circle -> Point2 -> Bool
contains Circle
circle Point2
q = Circle -> Either MathError Circle
forall a b. b -> Either a b
Right Circle
circle
        | Bool
otherwise = (Circle -> Point2 -> Either MathError Circle)
-> Circle -> [Point2] -> Either MathError Circle
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM (Point2 -> Point2 -> Circle -> Point2 -> Either MathError Circle
inner Point2
p Point2
q) (Point2 -> Point2 -> Circle
diameter Point2
p Point2
q) (Int -> [Point2] -> [Point2]
forall a. Int -> [a] -> [a]
take Int
j [Point2]
ps)
    inner :: Point2 -> Point2 -> Circle -> Point2 -> Either MathError Circle
inner Point2
p Point2
q Circle
circle Point2
r
        | Circle -> Point2 -> Bool
contains Circle
circle Point2
r = Circle -> Either MathError Circle
forall a b. b -> Either a b
Right Circle
circle
        | Point2 -> Point2 -> Point2 -> Ordering
orient Point2
p Point2
q Point2
r Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
EQ =
            Circle -> Either MathError Circle
forall a b. b -> Either a b
Right
                ((Circle -> Circle -> Ordering) -> [Circle] -> Circle
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
maximumBy ((Circle -> Double) -> Circle -> Circle -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (\(Circle Point2
_ Double
radius) -> Double
radius)) [Point2 -> Point2 -> Circle
diameter Point2
p Point2
q, Point2 -> Point2 -> Circle
diameter Point2
p Point2
r, Point2 -> Point2 -> Circle
diameter Point2
q Point2
r])
        | Bool
otherwise = do center <- Triangle -> Either MathError Point2
circumcenter (Point2 -> Point2 -> Point2 -> Triangle
Triangle Point2
p Point2
q Point2
r); Right (Circle center (norm2 (sub2 p center)))

{- | Minimum-area oriented enclosing rectangle. Testing every hull edge uses
the rotating-calipers extremal-orientation theorem; the reference scans all
vertices for each orientation (O(h^2)), at most 4096 points.
-}
minimumBoundingRectangle :: [Point2] -> Either MathError [Point2]
minimumBoundingRectangle :: [Point2] -> Either MathError [Point2]
minimumBoundingRectangle [Point2]
ps = do
    hull <- HullAlgorithm -> [Point2] -> Either MathError [Point2]
convexHull HullAlgorithm
MonotoneChain [Point2]
ps
    require (not (null hull)) "empty bounding rectangle"
    Right $
        if length hull == 1
            then replicate 4 (firstOr (P2 0 0) hull)
            else snd (minimumBy (comparing fst) [candidate a b | (a, b) <- edges hull])
  where
    candidate :: Point2 -> Point2 -> (Double, [Point2])
candidate Point2
a Point2
b =
        let d :: Point2
d = Point2 -> Point2 -> Point2
sub2 Point2
b Point2
a
            P2 Double
dx Double
dy = Point2
d
            size :: Double
size = Point2 -> Double
norm2 Point2
d
            ux :: Double
ux = Double
dx Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
size
            uy :: Double
uy = Double
dy Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
size
            u :: Point2
u = Double -> Double -> Point2
P2 Double
ux Double
uy
            v :: Point2
v = Double -> Double -> Point2
P2 (-Double
uy) Double
ux
            xs :: [Double]
xs = (Point2 -> Double) -> [Point2] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Point2 -> Point2 -> Double
dot2 Point2
u) [Point2]
ps
            ys :: [Double]
ys = (Point2 -> Double) -> [Point2] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Point2 -> Point2 -> Double
dot2 Point2
v) [Point2]
ps
            lo :: Double
lo = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
xs
            hi :: Double
hi = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
xs
            bottom :: Double
bottom = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
ys
            top :: Double
top = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
ys
         in ( (Double
hi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lo) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
top Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
bottom)
            , [Point2 -> Point2 -> Point2
add2 (Double -> Point2 -> Point2
scale2 Double
x Point2
u) (Double -> Point2 -> Point2
scale2 Double
y Point2
v) | (Double
x, Double
y) <- [(Double
lo, Double
bottom), (Double
hi, Double
bottom), (Double
hi, Double
top), (Double
lo, Double
top)]]
            )

{- | Convex Minkowski sum, via pairwise sums followed by a hull. Nonconvex
inputs are deliberately convexified; at most 64 vertices per operand.
-}
minkowskiSum :: [Point2] -> [Point2] -> Either MathError [Point2]
minkowskiSum :: [Point2] -> [Point2] -> Either MathError [Point2]
minkowskiSum [Point2]
a [Point2]
b = do
    Bool -> String -> Either MathError ()
require
        (Bool -> Bool
not ([Point2] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Point2]
a Bool -> Bool -> Bool
|| [Point2] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Point2]
b) Bool -> Bool -> Bool
&& [Point2] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Point2]
a Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
64 Bool -> Bool -> Bool
&& [Point2] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Point2]
b Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
64)
        String
"Minkowski sum requires 1..64 points per operand"
    (Point2 -> Either MathError ()) -> [Point2] -> Either MathError ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ Point2 -> Either MathError ()
check2 ([Point2]
a [Point2] -> [Point2] -> [Point2]
forall a. [a] -> [a] -> [a]
++ [Point2]
b)
    HullAlgorithm -> [Point2] -> Either MathError [Point2]
convexHull HullAlgorithm
MonotoneChain [Point2 -> Point2 -> Point2
add2 Point2
x Point2
y | Point2
x <- [Point2]
a, Point2
y <- [Point2]
b]