```{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE IncoherentInstances #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE NoMonomorphismRestriction #-}
module Math.Geometry
( module Math.Geometry
,module Data.VectorSpace
)
where

import Data.VectorSpace
import Data.List (foldl')

type Vertex2 a = (a, a)
vec2 :: Double -> Double -> Vertex2 Double
vec2 x y = (x, y)

vx :: VectorSpace (v, v) => (v, v) -> v
vx (x, _) = x

vy :: VectorSpace (v, v) => (v, v) -> v
vy (_, y) = y

type Edge2 a = (Vertex2 a, Vertex2 a)
type LineSegment = Edge2 Double

-- "intersects" is a commutative binary predicate on two shapes.
class Intersectable s0 s1 where
intersects :: s0 -> s1 -> Bool

--instance Intersectable s0 s1 => Intersectable s1 s0 where
--   intersects s1 s0 = intersects s0 s1

intersections e es = filter (intersects e) es

-- A Boundary describes an inclusive lower bound with a corner point and an exclusive upper bound
-- with a size.
data Boundary = Boundary
{
boundary_corner :: Vertex2 Double,
boundary_size   :: Double
}
deriving (Show)

-- Boundaries b0 and b1 intersect if the min extent of the intersection of b1 with (the plane +x
-- including b0.p unioned with the plane +y including b0.p) is within b0.
instance Intersectable Boundary Boundary where
intersects b0 b1 =
let c = union_min_extent_planes_of b0
in if b1 `intersects` c
then let (Boundary p _) = intersection b1 c
in intersects b0 p
else False

-- The union of the min extent planes of a boundary is represented by the min extent
newtype MinExtentPlanes = MinExtentPlanes (Vertex2 Double)
union_min_extent_planes_of (Boundary p _) = MinExtentPlanes p

-- A boundary intersects the min extent planes if the far extent extent of the boundary is within
-- the range defined by the min extent planes.
-- The comparison is > and not >= since the far extent is the point just beyond the boundary. Which
-- needs to be just inside the planes in order for the boundary to be inside the planes.
instance Intersectable Boundary MinExtentPlanes where
intersects b (MinExtentPlanes (min_x, min_y)) =
let (_, (b_max_x, b_max_y)) = boundary_extents b
in (b_max_x > min_x) && (b_max_y > min_y)

intersection :: Boundary -> MinExtentPlanes -> Boundary
intersection (Boundary p size) (MinExtentPlanes min_p) = Boundary (ext_max min_p p) size

instance Intersectable Boundary LineSegment where
intersects b l@(p0, p1) =
-- If any point of the line segment is contained in the boundary then the line segment intersects the
-- element.
intersects b p0 || intersects b p1
-- If niether point is in the element the line segment could still intersect the boundary. The line
-- segment must, in this case, intersect an edge of the boundary.
|| any (intersects l) (boundary_edges b)

--The equations for line intersection are pulled from
--  http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/
-- without much thought.

instance Intersectable LineSegment LineSegment where
intersects (p0a, p0b) (p1a, p1b) =
let x1 = vx p0a
y1 = vy p0a
x2 = vx p0b
y2 = vy p0b
x3 = vx p1a
y3 = vy p1a
x4 = vx p1b
y4 = vy p1b
div = (y4 - y3)*(x2 - x1) - (x4 - x3)*(y2 - y1)
in if div < 1e-9
then False
else
let t0n = (x4 - x3)*(y1 - y3) - (y4 - y3)*(x1 - x3)
t0 = t0n / div
t1n = (x2 - x1)*(y1 - y3) - (y2 - y1)*(x1 - x3)
t1 = t1n / div
in t0 > 0.0 && t0 < 1.0 && t1 > 0.0 && t1 < 1.0

union_boundaries :: Boundary -> Boundary -> Boundary
union_boundaries b0 b1 =
let (min0, max0) = boundary_extents b0
(min1, max1) = boundary_extents b1
p = ext_min min0 min1
ext = ext_max max0 max1
(w,h) = ext ^-^ p
size = max w h
in Boundary p size

ext_min (x0,y0) (x1,y1) = (min x0 x1, min y0 y1)
ext_max (x0,y0) (x1,y1) = (max x0 x1, max y0 y1)

--instance Show Boundary where
--    show (Boundary p size) = show p ++ " -> " ++ show size

instance Intersectable Boundary (Vertex2 Double) where
intersects bounds (px, py) =
let (x, y) = boundary_corner bounds
s = boundary_size bounds
in px < (x + s) && px >= x && py < (y + s) && py >= y

class HasBoundary s where
boundary_points :: s -> [Vertex2 Double]
boundary_edges :: s -> [Edge2 Double]
boundary_edges s =
let ps@(p0 : ps') = boundary_points s
in zip ps (ps' ++ [p0])
boundary_extents :: s -> (Vertex2 Double, Vertex2 Double)
boundary_extents s =
let (p0 : ps) = boundary_points s
initial_min_extent = p0
initial_max_extent = p0
union_extents ((min_x, min_y), (max_x,max_y)) (x, y) =
let min_x' = min min_x x
min_y' = min min_y y
max_x' = max max_x x
max_y' = max max_y y
in ((min_x', min_y'), (max_x', max_y'))
in foldl' union_extents (initial_min_extent, initial_max_extent) ps
boundary_square :: s -> Boundary
boundary_square s =
let (min_extent, max_extent) = boundary_extents s
width  = fst max_extent - fst min_extent
height = snd max_extent - snd min_extent
size = max width height
in Boundary (fst min_extent, snd min_extent) size

-- A boundary cleary has itself as it's boundary.
instance HasBoundary Boundary where
boundary_points (Boundary p s) =
[ p
, p ^+^ (0, s)
, p ^+^ (s, s)
, p ^+^ (s, 0)
]
boundary_extents (Boundary p s) = (p, p ^+^ (s,s))
boundary_square b = b

```