-- | CG library (minus).
module Data.CG.Minus where

import Data.Complex {- base -}
import Data.Maybe {- base -}
import Text.Printf {- base -}

-- * Types

-- | Two-dimensional point.
--
-- Pt are 'Num', pointwise, ie:
--
-- > Pt 1 2 + Pt 3 4 == Pt 4 6
-- > Pt 1 2 * Pt 3 4 == Pt 3 8
-- > negate (Pt 0 1) == Pt 0 (-1)
-- > abs (Pt (-1) 1) == Pt 1 1
-- > signum (Pt (-1/2) (1/2)) == Pt (-1) 1
data Pt a = Pt {pt_x :: a,pt_y :: a} deriving (Eq,Ord,Show)

instance Num a => Num (Pt a) where
    (+) = pt_binop (+)
    (-) = pt_binop (-)
    (*) = pt_binop (*)
    negate = pt_uop negate
    abs = pt_uop abs
    signum = pt_uop signum
    fromInteger n = let n' = fromInteger n in Pt n' n'

-- | Two-dimensional vector.  Vector are 'Num' in the same manner as
-- 'Pt'.
data Vc a = Vc {vc_x :: a,vc_y :: a} deriving (Eq,Ord,Show)

instance Num a => Num (Vc a) where
    (+) = vc_binop (+)
    (-) = vc_binop (-)
    (*) = vc_binop (*)
    negate = vc_uop negate
    abs = vc_uop abs
    signum = vc_uop signum
    fromInteger n = let n' = fromInteger n in Vc n' n'

-- | Two-dimensional line.
data Ln a = Ln {ln_start :: Pt a,ln_end :: Pt a} deriving (Eq,Ord,Show)

-- | Line segments.
type Ls a = [Pt a]

-- | Window, given by a /lower left/ 'Pt' and an /extent/ 'Vc'.
data Wn a = Wn {wn_ll :: Pt a,wn_ex :: Vc a} deriving (Eq,Show)

-- | Real number, synonym for 'Double'.
type R = Double

-- * R(eal) functions

-- | Epsilon.
epsilon :: Floating n => n
epsilon = 0.000001

-- | Is absolute difference less than 'epsilon'.
(~=) :: (Floating a, Ord a) => a -> a -> Bool
p ~= q = abs (p - q) < epsilon

-- | Degrees to radians.
--
-- > map r_to_radians [-180,-90,0,90,180] == [-pi,-pi/2,0,pi/2,pi]
r_to_radians :: R -> R
r_to_radians x = (x / 180) * pi

-- | Radians to degrees, inverse of 'r_to_radians'.
--
-- > map r_from_radians [-pi,-pi/2,0,pi/2,pi] == [-180,-90,0,90,180]
r_from_radians :: R -> R
r_from_radians x = (x / pi) * 180

-- | 'R' modulo within range.
--
-- > map (r_constrain (3,5)) [2.75,5.25] == [4.75,3.25]
r_constrain :: (R,R) -> R -> R
r_constrain (l,r) =
    let down n i x = if x > i then down n i (x - n) else x
        up n i x = if x < i then up n i (x + n) else x
        both n i j x = up n i (down n j x)
    in both (r - l) l r

-- | Sum of squares.
mag_sq :: Num a => a -> a -> a
mag_sq x y = x * x + y * y

-- | 'sqrt' of 'mag_sq'.
mag :: Floating c => c -> c -> c
mag x = sqrt . mag_sq x

-- * Pt functions

-- | Tuple constructor.
pt' :: (a,a) -> Pt a
pt' (x,y) = Pt x y

-- | Tuple accessor.
pt_xy :: Pt t -> (t, t)
pt_xy (Pt x y) = (x,y)

-- | 'Pt' of (0,0).
--
-- > pt_origin == Pt 0 0
pt_origin :: Num a => Pt a
pt_origin = Pt 0 0

-- | Unary operator at 'Pt', ie. basis for 'Num' instances.
pt_uop :: (a -> b) -> Pt a -> Pt b
pt_uop f (Pt x y) = Pt (f x) (f y)

-- | Binary operator at 'Pt', ie. basis for 'Num' instances.
pt_binop :: (a -> b -> c) -> Pt a -> Pt b -> Pt c
pt_binop f (Pt x1 y1) (Pt x2 y2) = Pt (x1 `f` x2) (y1 `f` y2)

-- | 'Pt' at /(n,n)/.
--
-- > pt_from_scalar 1 == Pt 1 1
pt_from_scalar :: (Num a) => a -> Pt a
pt_from_scalar a = Pt a a

-- | Clip /x/ and /y/ to lie in /(0,n)/.
--
-- > pt_clipu 1 (Pt 0.5 1.5) == Pt 0.5 1
pt_clipu :: (Ord a,Num a) => a -> Pt a -> Pt a
pt_clipu u =
    let f n = if n < 0 then 0 else if n > u then u else n
    in pt_uop f

-- | Swap /x/ and /y/ coordinates at 'Pt'.
--
-- > pt_swap (Pt 1 2) == Pt 2 1
pt_swap :: Pt a -> Pt a
pt_swap (Pt x y) = Pt y x

-- | Negate /y/ element of 'Pt'.
--
-- > pt_negate_y (Pt 1 1) == Pt 1 (-1)
pt_negate_y :: (Num a) => Pt a -> Pt a
pt_negate_y (Pt x y) = Pt x (negate y)

-- | 'Pt' variant of 'r_to_radians'.
--
-- > pt_to_radians (Pt 90 270) == Pt (pi/2) (pi*(3/2))
pt_to_radians :: Pt R -> Pt R
pt_to_radians = pt_uop r_to_radians

-- | Cartesian to polar.
--
-- > pt_to_polar (Pt 0 pi) == Pt pi (pi/2)
pt_to_polar :: Pt R -> Pt R
pt_to_polar (Pt x y) =
    let (x',y') = polar (x :+ y)
    in Pt x' y'

-- | Polar to cartesian, inverse of 'pt_to_polar'.
--
-- > pt_from_polar (Pt pi (pi/2)) ~= Pt 0 pi
pt_from_polar :: Pt R -> Pt R
pt_from_polar (Pt mg ph) =
    let c = mkPolar mg ph
    in Pt (realPart c) (imagPart c)

-- | Scalar 'Pt' '+'.
--
-- > pt_offset 1 pt_origin == Pt 1 1
pt_offset :: Num a => a -> Pt a -> Pt a
pt_offset = pt_uop . (+)

-- | Scalar 'Pt' '*'.
--
-- > pt_scale 2 (Pt 1 2) == Pt 2 4
pt_scale :: Num a => a -> Pt a -> Pt a
pt_scale = pt_uop . (*)

-- | Pointwise 'min'.
pt_min :: (Ord a) => Pt a -> Pt a -> Pt a
pt_min = pt_binop min

-- | Pointwise 'max'.
pt_max :: (Ord a) => Pt a -> Pt a -> Pt a
pt_max = pt_binop max

-- | Apply function to /x/ and /y/ fields of three 'Pt'.
pt_ternary_f :: (a->a->b->b->c->c->d) -> Pt a -> Pt b -> Pt c -> d
pt_ternary_f f (Pt x0 y0) (Pt x1 y1) (Pt x2 y2) = f x0 y0 x1 y1 x2 y2

-- | Given a /(minima,maxima)/ pair, expand so as to include /p/.
--
-- > pt_minmax (Pt 0 0,Pt 1 1) (Pt (-1) 2) == (Pt (-1) 0,Pt 1 2)
pt_minmax :: Ord a => (Pt a,Pt a) -> Pt a -> (Pt a,Pt a)
pt_minmax (p0,p1) p =
    let f x0 y0 x1 y1 x y =
            (Pt (min x x0) (min y y0)
            ,Pt (max x x1) (max y y1))
    in pt_ternary_f f p0 p1 p

-- | 'Pt' variant of 'constrain'.
pt_constrain :: (Pt R,Pt R) -> Pt R -> Pt R
pt_constrain (p0,p1) p =
    let f x0 y0 x1 y1 x y =
            let x' = r_constrain (x0,x1) x
                y' = r_constrain (y0,y1) y
            in Pt x' y'
    in pt_ternary_f f p0 p1 p

-- | Angle to origin.
--
-- > pt_angle_o (Pt 0 1) == pi / 2
pt_angle_o :: Pt R -> R
pt_angle_o (Pt x y) = atan2 y x

-- | Angle from /p/ to /q/.
--
-- > pt_angle (Pt 0 (-1)) (Pt 0 1) == pi/2
-- > pt_angle (Pt 1 0) (Pt 0 1) == pi * 3/4
-- > pt_angle (Pt 0 1) (Pt 0 1) == 0
pt_angle :: Pt R -> Pt R -> R
pt_angle p q = pt_angle_o (q - p)

-- | Pointwise '+'.
--
-- pt_translate (Pt 0 0) (vc 1 1) == pt 1 1
pt_translate :: (Num a,Eq a) => Pt a -> Vc a -> Pt a
pt_translate (Pt x y) (Vc dx dy) = Pt (x + dx) (y + dy)

-- | 'pt_uop' 'fromIntegral'.
pt_from_i :: (Integral i,Num a) => Pt i -> Pt a
pt_from_i = pt_uop fromIntegral

-- | 'mag_sq' of /x/ /y/ of 'Pt'.
pt_mag_sq :: Num a => Pt a -> a
pt_mag_sq (Pt x y) = mag_sq x y

-- | 'mag' of /x/ /y/ of 'Pt'.
pt_mag :: Floating a => Pt a -> a
pt_mag (Pt x y) = mag x y

-- | Distance from 'Pt' /p/ to 'Pt' /q/.
--
-- > pt_distance (Pt 0 0) (Pt 0 1) == 1
-- > pt_distance (Pt 0 0) (Pt 1 1) == sqrt 2
pt_distance :: (Floating a,Eq a) => Pt a -> Pt a -> a
pt_distance (Pt x1 y1) (Pt x2 y2) = pt_mag (Pt (x2 - x1) (y2 - y1))

-- | Are /x/ and /y/ of 'Pt' /p/ in range (0,1).
--
-- > map pt_is_normal [Pt 0 0,Pt 1 1,Pt 2 2] == [True,True,False]
pt_is_normal :: (Ord a,Num a) => Pt a -> Bool
pt_is_normal (Pt x y) = x >= 0 && x <= 1 && y >= 0 && y <= 1

-- | Rotate 'Pt' /n/ radians.
--
-- > pt_rotate pi (Pt 1 0) ~= Pt (-1) 0
pt_rotate :: Floating a => a -> Pt a -> Pt a
pt_rotate a (Pt x y) =
    let s = sin a
        c = cos a
    in Pt (x * c - y * s) (y * c + x * s)

-- * Vc functions

-- | Unary operator at 'Vc', ie. basis for 'Num' instances.
vc_uop :: (a -> b) -> Vc a -> Vc b
vc_uop f (Vc x y) = Vc (f x) (f y)

-- | Binary operator at 'Vc', ie. basis for 'Num' instances.
vc_binop :: (a -> b -> c) -> Vc a -> Vc b -> Vc c
vc_binop f (Vc x1 y1) (Vc x2 y2) = Vc (x1 `f` x2) (y1 `f` y2)

-- | 'mag_sq' of 'Vc'.
vc_mag_sq :: Floating c => Vc c -> c
vc_mag_sq (Vc dx dy) = mag_sq dx dy

-- | 'mag' of 'Vc'.
vc_mag :: Floating c => Vc c -> c
vc_mag (Vc dx dy) = mag dx dy

-- | Multiply 'Vc' pointwise by scalar.
--
-- > vc_scale 2 (Vc 3 4) == Vc 6 8
vc_scale :: Num a => a -> Vc a -> Vc a
vc_scale n (Vc x y) = Vc (x * n) (y * n)

-- | 'Vc' dot product.
--
-- > vc_dot (Vc 1 2) (Vc 3 4) == 11
vc_dot :: Num a => Vc a -> Vc a -> a
vc_dot (Vc x y) (Vc x' y') = (x * x') + (y * y')

-- | Scale 'Vc' to have unit magnitude (to within tolerance).
--
-- > vc_unit (Vc 1 1) ~= let x = (sqrt 2) / 2 in Vc x x
vc_unit :: (Ord a, Floating a) => Vc a -> Vc a
vc_unit v
    | abs (vc_mag_sq v - 1) < epsilon = v
    | vc_mag_sq v == 0 = v
    | otherwise = let Vc x y = v
                      m = mag x y
                  in Vc (x / m) (y / m)

-- | The angle between two vectors on a plane. The angle is from v1 to
-- v2, positive anticlockwise.  The result is in (-pi,pi)
vc_angle :: Vc R -> Vc R -> R
vc_angle (Vc x1 y1) (Vc x2 y2) =
    let t1 = atan2 y1 x1
        t2 = atan2 y2 x2
    in r_constrain (-pi,pi) (t2 - t1)

-- * Line functions

-- | Variant on 'Ln' which takes 'Pt' co-ordinates as duples.
--
-- > ln' (0,0) (1,1) == Ln (Pt 0 0) (Pt 1 1)
-- > ln_start (Ln (Pt 0 0) (Pt 1 1)) == Pt 0 0
-- > ln_end (Ln (Pt 0 0) (Pt 1 1)) == Pt 1 1
ln' :: (Num a,Eq a) => (a,a) -> (a,a) -> Ln a
ln' (x1,y1) (x2,y2) = Ln (Pt x1 y1) (Pt x2 y2)

-- | 'Vc' that 'pt_translate's start 'Pt' to end 'Pt' of 'Ln'.
--
-- > let l = Ln (Pt 0 0) (Pt 1 1)
-- > in ln_start l `pt_translate` ln_vc l == Pt 1 1
ln_vc :: (Num a,Eq a) => Ln a -> Vc a
ln_vc (Ln p q) = let Pt x y = q - p in Vc x y

-- | 'Pt' UOp at 'Ln'.
ln_uop :: (Pt a -> Pt b) -> Ln a -> Ln b
ln_uop f (Ln l r) = Ln (f l) (f r)

-- | 'pt_scale' at 'Ln'.
ln_scale :: Num b => b -> Ln b -> Ln b
ln_scale m = ln_uop (pt_scale m)

-- | The angle, in /radians/, anti-clockwise from the /x/-axis.
--
-- > ln_angle (ln' (0,0) (0,0)) == 0
-- > ln_angle (ln' (0,0) (1,1)) == pi/4
-- > ln_angle (ln' (0,0) (0,1)) == pi/2
-- > ln_angle (ln' (0,0) (-1,1)) == pi * 3/4
ln_angle :: Ln R -> R
ln_angle ln =
    let Vc dx dy = ln_vc ln
    in if dx == 0 && dy == 0 then 0 else atan2 dy dx

-- | Start and end points of 'Ln'.
--
-- > ln_pt (Ln (Pt 1 0) (Pt 0 0)) == (Pt 1 0,Pt 0 0)
ln_pt :: (Num a,Eq a) => Ln a -> (Pt a,Pt a)
ln_pt (Ln s e) = (s,e)

-- | Variant of 'ln_pt' giving co-ordinates as duples.
--
-- > ln_pt' (Ln (Pt 1 0) (Pt 0 0)) == ((1,0),(0,0))
ln_pt' :: (Num a,Eq a) => Ln a -> ((a,a),(a,a))
ln_pt' (Ln (Pt x1 y1) (Pt x2 y2)) = ((x1,y1),(x2,y2))

-- | Midpoint of a 'Ln'.
--
-- > ln_midpoint (Ln (Pt 0 0) (Pt 2 1)) == Pt 1 (1/2)
ln_midpoint :: (Fractional a,Eq a) => Ln a -> Pt a
ln_midpoint (Ln (Pt x1 y1) (Pt x2 y2)) =
    let x = (x1 + x2) / 2
        y = (y1 + y2) / 2
    in Pt x y

-- | Variant on 'ln_midpoint'.
--
-- > cc_midpoint (Just (Pt 0 0),Nothing) == Pt 0 0
-- > cc_midpoint (Nothing,Just (Pt 2 1)) == Pt 2 1
-- > cc_midpoint (Just (Pt 0 0),Just (Pt 2 1)) == Pt 1 (1/2)
cc_midpoint :: (Maybe (Pt R), Maybe (Pt R)) -> Pt R
cc_midpoint cc =
    case cc of
      (Nothing,Nothing) -> Pt 0 0
      (Just p,Nothing) -> p
      (Nothing, Just q) -> q
      (Just p, Just q) -> ln_midpoint (Ln p q)

-- | Magnitude of 'Ln', ie. length of line.
--
-- > ln_magnitude (Ln (Pt 0 0) (Pt 1 1)) == sqrt 2
-- > pt_x (pt_to_polar (Pt 1 1)) == sqrt 2
ln_magnitude :: Ln R -> R
ln_magnitude = vc_mag . ln_vc

-- | Order 'Pt' at 'Ln' so that /p/ is to the left of /q/.  If /x/
-- fields are equal, sort on /y/.
--
-- > ln_sort (Ln (Pt 1 0) (Pt 0 0)) == Ln (Pt 0 0) (Pt 1 0)
-- > ln_sort (Ln (Pt 0 1) (Pt 0 0)) == Ln (Pt 0 0) (Pt 0 1)
ln_sort :: (Num a,Ord a) => Ln a -> Ln a
ln_sort ln =
    let Ln p q = ln
        Pt x1 y1 = p
        Pt x2 y2 = q
    in case compare x1 x2 of
         LT -> ln
         EQ -> if y1 <= y2 then ln else Ln q p
         GT -> Ln q p

-- | Adjust 'Ln' to have equal starting 'Pt' but magnitude 'R'.
--
-- > ln_adjust (sqrt 2) (Ln (Pt 0 0) (Pt 2 2)) == Ln (Pt 0 0) (Pt 1 1)
ln_adjust :: (Floating a, Ord a) => a -> Ln a -> Ln a
ln_adjust z ln =
    let Ln p _ = ln
        v = vc_scale z (vc_unit (ln_vc ln))
    in Ln p (pt_translate p v)

-- | Extend 'Ln' by 'R', ie. 'ln_adjust' with /n/ added to
-- 'ln_magnitude'.
--
-- > ln_extend (sqrt 2) (Ln (Pt 0 0) (Pt 1 1)) ~= Ln (Pt 0 0) (Pt 2 2)
ln_extend :: R -> Ln R -> Ln R
ln_extend n l = Ln (ln_start l) (pt_linear_extension n l)

-- | Variant definition of 'ln_extend'.
--
-- > ln_extend_ (sqrt 2) (Ln (Pt 0 0) (Pt 1 1)) == Ln (Pt 0 0) (Pt 2 2)
ln_extend_ :: R -> Ln R -> Ln R
ln_extend_ n l = ln_adjust (n + ln_magnitude l) l

-- | Calculate the point that extends a line by length 'n'.
--
-- > pt_linear_extension (sqrt 2) (Ln (Pt 1 1) (Pt 2 2)) ~= Pt 3 3
-- > pt_linear_extension 1 (Ln (Pt 1 1) (Pt 1 2)) ~= Pt 1 3
pt_linear_extension :: R -> Ln R -> Pt R
pt_linear_extension n (Ln p q) =
    let Pt mg ph = pt_to_polar (q - p)
    in pt_from_polar (Pt (mg + n) ph) + p

-- | Does 'Pt' /p/ lie on 'Ln' (inclusive).
--
-- > let {f = pt_on_line (Ln (Pt 0 0) (Pt 1 1))
-- >     ;r = [True,False,False,True]}
-- > in map f [Pt 0.5 0.5,Pt 2 2,Pt (-1) (-1),Pt 0 0] == r
pt_on_line :: Ln R -> Pt R -> Bool
pt_on_line l r =
    let (p,q) = ln_pt l
        Pt i j = pt_to_polar (q - p)
        Pt i' j' = pt_to_polar (r - p)
    in r == p || r == q || (j == j' && i' <= i)

-- * Intersection

ln_intersect :: (Eq t, Fractional t) => Ln t -> Ln t -> Maybe (t,t)
ln_intersect l1 l2 =
    let Ln (Pt x1 y1) _ = l1
        Vc dx1 dy1 = ln_vc l1
        Ln (Pt x2 y2) _ = l2
        Vc dx2 dy2 = ln_vc l2
        a = (dx2 * dy1) - (dx1 * dy2)
        t' = ((dx1 * (y2 - y1)) - (dy1 * (x2 - x1))) / a
        t = ((dx2 * (y1 - y2)) - (dy2 * (x1 - x2))) / (negate a)
    in if a == 0 then Nothing else Just (t,t')

ln_pt_along :: (Eq a, Num a) => a -> Ln a -> Pt a
ln_pt_along z ln =
    let v = vc_scale z (ln_vc ln)
        Ln p _ = ln
    in pt_translate p v

-- | Do two 'Ln's intersect, and if so at which 'Pt'.
--
-- > ln_intersection (ln' (0,0) (5,5)) (ln' (5,0) (0,5)) == Just (Pt 2.5 2.5)
-- > ln_intersection (ln' (1,3) (9,3)) (ln' (0,1) (2,1)) == Nothing
-- > ln_intersection (ln' (1,5) (6,8)) (ln' (0.5,3) (6,4)) == Nothing
-- > ln_intersection (ln' (1,2) (3,6)) (ln' (2,4) (4,8)) == Nothing
-- > ln_intersection (ln' (2,3) (7,9)) (ln' (1,2) (5,7)) == Nothing
-- > ln_intersection (ln' (0,0) (1,1)) (ln' (0,0) (1,0)) == Just (Pt 0 0)
ln_intersection :: (Ord a,Fractional a) => Ln a -> Ln a -> Maybe (Pt a)
ln_intersection l0 l1 =
    case ln_intersect l0 l1 of
      Nothing -> Nothing
      Just (i,j) -> if i >= 0 && i <= 1 && j >= 0 && j <= 1
                    then Just (ln_pt_along i l0)
                    else Nothing

-- | Variant definition of 'ln_intersection', using algorithm at
-- <http://paulbourke.net/geometry/lineline2d/>.
--
-- > ln_intersection_ (ln' (1,2) (3,6)) (ln' (2,4) (4,8)) == Nothing
-- > ln_intersection_ (ln' (0,0) (1,1)) (ln' (0,0) (1,0)) == Just (Pt 0 0)
ln_intersection_ :: (Ord a,Fractional a) => Ln a -> Ln a -> Maybe (Pt a)
ln_intersection_ l0 l1 =
    let ((x1,y1),(x2,y2)) = ln_pt' l0
        ((x3,y3),(x4,y4)) = ln_pt' l1
        d = (y4 - y3) * (x2 - x1) - (x4 - x3) * (y2 - y1)
        ua' = (x4 - x3) * (y1 - y3) - (y4 - y3) * (x1 - x3)
        ub' = (x2 - x1) * (y1 - y3) - (y2 - y1) * (x1 - x3)
    in if d == 0
       then Nothing
       else if ua' == 0 && ub' == 0
            then Just (Pt x1 y1)
            else let ua = ua' / d
                     ub = ub' / d
                 in if in_range 0 1 ua && in_range 0 1 ub
                    then let x = x1 + ua * (x2 - x1)
                             y = y1 + ua * (y2 - y1)
                         in Just (Pt x y)
                    else Nothing

-- | Predicate variant of 'ln_intersection'.
--
-- > ln_intersect_p (ln' (1,1) (3,8)) (ln' (0.5,2) (4,7)) == True
-- > ln_intersect_p (ln' (3.5,9) (3.5,0.5)) (ln' (3,1) (9,1)) == True
ln_intersect_p :: (Ord a, Fractional a) => Ln a -> Ln a -> Bool
ln_intersect_p l = isJust . ln_intersection l

-- * Line slope

-- | Slope of 'Ln' or 'Nothing' if /vertical/.
--
-- > let l = zipWith ln' (repeat (0,0)) [(1,0),(2,1),(1,1),(0,1),(-1,1)]
-- > in map ln_slope l == [Just 0,Just (1/2),Just 1,Nothing,Just (-1)]
ln_slope :: (Fractional a,Eq a) => Ln a -> Maybe a
ln_slope l =
    let ((x1,y1),(x2,y2)) = ln_pt' l
    in case x2 - x1 of
         0 -> Nothing
         dx -> Just ((y2 - y1) / dx)

-- | Are 'Ln's parallel, ie. have equal 'ln_slope'.  Note that the
-- direction of the 'Ln' is not relevant, ie. this is not equal to
-- 'ln_same_direction'.
--
-- > ln_parallel (ln' (0,0) (1,1)) (ln' (2,2) (1,1)) == True
-- > ln_parallel (ln' (0,0) (1,1)) (ln' (2,0) (1,1)) == False
-- > ln_parallel (ln' (1,2) (3,6)) (ln' (2,4) (4,8)) == True
-- > map ln_slope [ln' (2,2) (1,1),ln' (2,0) (1,1)] == [Just 1,Just (-1)]
ln_parallel :: (Ord a,Fractional a) => Ln a -> Ln a -> Bool
ln_parallel p q = ln_slope p == ln_slope q

-- | Are 'Ln's parallel, ie. have equal 'ln_angle'.
--
-- > ln_parallel_ (ln' (0,0) (1,1)) (ln' (2,2) (1,1)) == True
ln_parallel_ :: Ln R -> Ln R -> Bool
ln_parallel_ p q = ln_angle (ln_sort p) == ln_angle (ln_sort q)

-- | Are two vectors are in the same direction (to within a small
-- tolerance).
vc_same_direction :: (Ord a, Floating a) => Vc a -> Vc a -> Bool
vc_same_direction v w =
    let Vc dx1 dy1 = vc_unit v
        Vc dx2 dy2 = vc_unit w
    in abs (dx2 - dx1) < epsilon && abs (dy2 - dy1) < epsilon

-- | Do 'Ln's have same direction (within tolerance).
--
-- > ln_same_direction (ln' (0,0) (1,1)) (ln' (0,0) (2,2)) == True
-- > ln_same_direction (ln' (0,0) (1,1)) (ln' (2,2) (0,0)) == False
ln_same_direction :: (Ord a, Floating a) => Ln a -> Ln a -> Bool
ln_same_direction p q = ln_vc p `vc_same_direction` ln_vc q

-- | Are 'Ln's parallel, ie. does 'ln_vc' of each equal 'ln_same_direction'.
--
-- > ln_parallel__ (ln' (0,0) (1,1)) (ln' (2,2) (1,1)) == True
ln_parallel__ :: Ln R -> Ln R -> Bool
ln_parallel__ p q = ln_vc (ln_sort p) `vc_same_direction` ln_vc (ln_sort q)

-- | Is 'Ln' horizontal, ie. is 'ln_slope' zero.
--
-- > ln_horizontal (ln' (0,0) (1,0)) == True
-- > ln_horizontal (ln' (1,0) (0,0)) == True
ln_horizontal :: (Fractional a,Eq a) => Ln a -> Bool
ln_horizontal = (== Just 0) . ln_slope

-- | Is 'Ln' vertical, ie. is 'ln_slope' 'Nothing'.
--
-- > ln_vertical (ln' (0,0) (0,1)) == True
ln_vertical :: (Fractional a,Eq a) => Ln a -> Bool
ln_vertical = (== Nothing) . ln_slope

-- * Ln sets

-- | 'pt_minmax' for set of 'Ln'.
lns_minmax :: [Ln R] -> (Pt R,Pt R)
lns_minmax = ls_minmax . concatMap (\(Ln l r) -> [l,r])

-- | Normalise to (0,m).
lns_normalise :: R -> [Ln R] -> [Ln R]
lns_normalise m l =
    let w = wn_from_extent (lns_minmax l)
    in map (ln_scale m . ln_normalise_w w) l

-- * L(ine) s(egment) functions

-- | 'Ls' constructor.
ls :: [Pt a] -> Ls a
ls = id

-- | Variant 'Ls' constructor from 'Pt' co-ordinates as duples.
ls' :: [(a,a)] -> Ls a
ls' = map (uncurry Pt)

-- | Negate /y/ elements.
ls_negate_y :: (Num a) => Ls a -> Ls a
ls_negate_y = map pt_negate_y

-- | Generate /minima/ and /maxima/ 'Point's from 'Ls'.
ls_minmax :: Ord a => Ls a -> (Pt a,Pt a)
ls_minmax s =
    case s of
      [] -> undefined
      p:ps -> foldl pt_minmax (p,p) ps

-- | Separate 'Ls' at points where the 'Vc' from one element to the
-- next exceeds the indicated distance.
--
-- > map length (ls_separate (Vc 2 2) (map (uncurry Pt) [(0,0),(1,1),(3,3)])) == [2,1]
ls_separate :: (Ord a,Num a) => Vc a -> Ls a -> [Ls a]
ls_separate (Vc dx dy) =
    let f (Pt x0 y0) (Pt x1 y1) = abs (x1 - x0) < dx &&
                                  abs (y1 - y0) < dy
    in segment_f f

-- | Delete 'Pt' from 'Ls' so that no two 'Pt' are within a tolerance
-- given by 'Vc'.
ls_tolerate :: (Ord a,Num a) => Vc a -> Ls a -> Ls a
ls_tolerate (Vc x y) =
    let too_close (Pt x0 y0) (Pt x1 y1) =
            let dx = abs (x1 - x0)
                dy = abs (y1 - y0)
            in dx < x && dy < y
    in delete_f too_close

-- | Variant of 'ls_tolerate' where 'Vc' is optional, and 'Nothing' gives 'id'.
ls_tolerate' :: (Ord a,Num a) => Maybe (Vc a) -> Ls a -> Ls a
ls_tolerate' i =
    case i of
      Nothing -> id
      Just i' -> ls_tolerate i'

-- | Test if point 'Pt' lies inside polygon 'Ls'.
--
-- > ls_pt_inside (ls' [(0,0),(1,0),(1,1),(0,1)]) (Pt 0.5 0.5) == True
ls_pt_inside :: Ls R -> Pt R -> Bool
ls_pt_inside s (Pt x y) =
    case s of
      [] -> undefined
      l0:l -> let xs = pairs ((l0:l)++[l0])
                  f (Pt x1 y1,Pt x2 y2) =
                      and [y > min y1 y2
                          ,y <= max y1 y2
                          ,x <= max x1 x2
                          ,y1 /= y2
                          ,x1 == x2 ||
                           x <= (y-y1)*(x2-x1)/(y2-y1)+x1]
              in odd (length (filter id (map f xs)))

-- | Variant that counts points at vertices as inside.
--
-- > ls_pt_inside' (ls' [(0,0),(1,0),(1,1),(0,1)]) (Pt 0 1) == True
ls_pt_inside' :: Ls R -> Pt R -> Bool
ls_pt_inside' l p = p `elem` l || ls_pt_inside l p

-- | Check all 'Pt' at 'Ls' are 'pt_is_normal'.
ls_check_normalised :: (Ord a,Num a) => Ls a -> Bool
ls_check_normalised s =
    case s of
      [] -> True
      p:z -> pt_is_normal p && ls_check_normalised z

-- | Line co-ordinates as /x/,/y/ list.
--
-- > ls_xy [Pt 0 0,Pt 1 1] == [0,0,1,1]
ls_xy :: Ls a -> [a]
ls_xy = concatMap (\(Pt x y) -> [x,y])

-- * Window

-- | Variant 'Wn' constructor.
wn' :: Num a => (a,a) -> (a,a) -> Wn a
wn' (x,y) (i,j) = Wn (Pt x y) (Vc i j)

-- | Extract /(x,y)/ and /(dx,dy)/ pairs.
--
-- > wn_extract (Wn (Pt 0 0) (Vc 1 1)) == ((0,0),(1,1))
wn_extract :: Wn a -> ((a,a),(a,a))
wn_extract (Wn (Pt x y) (Vc dx dy)) = ((x,y),(dx,dy))

-- | Show function for window with fixed precision of 'n'.
--
-- > wn_show 1 (Wn (Pt 0 0) (Vc 1 1)) == "((0.0,0.0),(1.0,1.0))"
wn_show :: Int -> Wn R -> String
wn_show n (Wn (Pt x0 y0) (Vc dx dy)) =
    let fs = printf "((%%.%df,%%.%df),(%%.%df,%%.%df))" n n n n
    in printf fs x0 y0 dx dy

-- | Is 'Pt' within 'Wn' exclusive of edge.
--
-- > map (pt_in_window (wn' (0,0) (1,1))) [Pt 0.5 0.5,Pt 1 1] == [True,False]
pt_in_window :: (Ord a,Num a) => Wn a -> Pt a -> Bool
pt_in_window (Wn (Pt lx ly) (Vc dx dy)) (Pt x y) =
    let (ux,uy) = (lx+dx,ly+dy)
    in x > lx && x < ux && y > ly && y < uy

-- | 'Wn' from /(lower-left,upper-right)/ extent.
wn_from_extent :: (Num a,Ord a) => (Pt a,Pt a) -> Wn a
wn_from_extent (Pt x0 y0,Pt x1 y1) = Wn (Pt x0 y0) (Vc (x1-x0) (y1-y0))

-- | 'Wn' containing 'Ls'.
--
-- > ls_window (ls' [(0,0),(1,1),(2,0)]) == wn' (0,0) (2,1)
ls_window :: (Num a,Ord a) => Ls a -> Wn a
ls_window = wn_from_extent . ls_minmax

-- | A 'Wn' that encompasses both input 'Wn's.
wn_join :: (Num a,Ord a) => Wn a -> Wn a -> Wn a
wn_join (Wn (Pt x0 y0) (Vc dx0 dy0)) (Wn (Pt x1 y1) (Vc dx1 dy1)) =
    let x = min x0 x1
        y = min y0 y1
        dx = max (x0+dx0) (x1+dx1) - x
        dy = max (y0+dy0) (y1+dy1) - y
    in Wn (Pt x y) (Vc dx dy)

-- | Predictate to determine if two 'Wn's intersect.
wn_intersect :: (Num a,Ord a) => Wn a -> Wn a -> Bool
wn_intersect w0 w1  =
    let (Wn (Pt x0 y0) (Vc dx0 dy0)) = w0
        (Wn (Pt x1 y1) (Vc dx1 dy1)) = w1
    in not (x0 > x1+dx1 || x1 > x0+dx0 || y0 > y1+dy1 || y1 > y0+dy0)

-- | Are all points at 'Ls' within the 'Wn'.
ls_in_window :: Wn R -> Ls R -> Bool
ls_in_window w = all (pt_in_window w)

-- | Are any points at 'Ls' within the window 'Wn'.
ls_enters_window :: Wn R -> Ls R -> Bool
ls_enters_window w = any (pt_in_window w)

-- | Are all points at 'Ls' outside the 'Wn'.
ls_not_in_window :: Wn R -> Ls R -> Bool
ls_not_in_window w = all (not . pt_in_window w)

-- | Break 'Ls' into segments that are entirely within the 'Wn'.
ls_segment_window :: Wn R -> Ls R -> [Ls R]
ls_segment_window w =
    let g [] = []
        g xs = let (i,xs') = span (pt_in_window w) xs
               in i : g (dropWhile (not . pt_in_window w) xs')
    in filter (not . null) . g

-- | Normalisation function for 'Wn', ie. map 'Pt' to lie within (0,1).
wn_normalise_f :: Wn R -> Pt R -> Pt R
wn_normalise_f (Wn (Pt x0 y0) (Vc dx dy)) (Pt x y) =
    let z = max dx dy
    in Pt ((x - x0) / z) ((y - y0) / z)

-- | Given 'Wn' normalise the 'Ls'.
ls_normalise_w :: Wn R -> Ls R -> Ls R
ls_normalise_w w = map (wn_normalise_f w)

-- | Given 'Wn' normalise 'Ln'.
ln_normalise_w :: Wn R -> Ln R -> Ln R
ln_normalise_w w (Ln p q) =
    let f = wn_normalise_f w
    in Ln (f p) (f q)

-- | Shift lower left 'Pt' of 'Wn' by indicated 'Pt'.
pt_shift_w :: Num a => Pt a -> Wn a -> Wn a
pt_shift_w p (Wn dp ex) = Wn (p + dp) ex

-- | Negate /y/ field of lower left 'Pt' of 'Wn'.
wn_negate_y :: Num a => Wn a -> Wn a
wn_negate_y (Wn p v) = Wn (pt_negate_y p) v

-- * Matrix

-- | Transformation matrix data type.
data Matrix n = Matrix n n n n n n deriving (Eq,Show)

-- | Enumeration of 'Matrix' indices.
data Matrix_Index = I0 | I1 | I2

mx_row :: Num n => Matrix n -> Matrix_Index -> (n,n,n)
mx_row (Matrix a b c d e f) i =
    case i of
      I0 -> (a,b,0)
      I1 -> (c,d,0)
      I2 -> (e,f,1)

mx_col :: Num n => Matrix n -> Matrix_Index -> (n,n,n)
mx_col (Matrix a b c d e f) i =
    case i of
      I0 -> (a,c,e)
      I1 -> (b,d,f)
      I2 -> (0,0,1)

mx_multiply :: Num n => Matrix n -> Matrix n -> Matrix n
mx_multiply a b =
    let f i j = let (r1,r2,r3) = mx_row a i
                    (c1,c2,c3) = mx_col b j
                in r1 * c1 + r2 * c2 + r3 * c3
    in Matrix (f I0 I0) (f I0 I1) (f I1 I0) (f I1 I1) (f I2 I0) (f I2 I1)

-- | Pointwise unary operator.
mx_uop :: (n -> n) -> Matrix n -> Matrix n
mx_uop g (Matrix a b c d e f) =
    Matrix (g a) (g b) (g c) (g d) (g e) (g f)

-- | Pointwise binary operator.
mx_binop :: (n -> n -> n) -> Matrix n -> Matrix n -> Matrix n
mx_binop g (Matrix a b c d e f) (Matrix a' b' c' d' e' f') =
    Matrix (g a a') (g b b') (g c c') (g d d') (g e e') (g f f')

instance Num n => Num (Matrix n) where
    (*) = mx_multiply
    (+) = mx_binop (+)
    (-) = mx_binop (-)
    abs = mx_uop abs
    signum = mx_uop signum
    fromInteger n = let n' = fromInteger n
                    in Matrix n' 0 0 n' 0 0

-- | A translation matrix with independent x and y offsets.
mx_translation :: Num n => n -> n -> Matrix n
mx_translation = Matrix 1 0 0 1

-- | A scaling matrix with independent x and y scalars.
mx_scaling :: Num n => n -> n -> Matrix n
mx_scaling x y = Matrix x 0 0 y 0 0

-- | A rotation matrix through the indicated angle (in radians).
mx_rotation :: Floating n => n -> Matrix n
mx_rotation a =
    let c = cos a
        s = sin a
        t = negate s
    in Matrix c s t c 0 0

-- | The identity matrix.
mx_identity :: Num n => Matrix n
mx_identity = Matrix 1 0 0 1 0 0

mx_translate :: Num n => n -> n -> Matrix n -> Matrix n
mx_translate x y m = m * (mx_translation x y)

mx_scale :: Num n => n -> n -> Matrix n -> Matrix n
mx_scale x y m = m * (mx_scaling x y)

mx_rotate :: Floating n => n -> Matrix n -> Matrix n
mx_rotate r m = m * (mx_rotation r)

mx_scalar_multiply :: Num n => n -> Matrix n -> Matrix n
mx_scalar_multiply scalar = mx_uop (* scalar)

mx_adjoint :: Num n => Matrix n -> Matrix n
mx_adjoint (Matrix a b c d x y) =
    Matrix d (-b) (-c) a (c * y - d * x) (b * x - a * y)

mx_invert :: Fractional n => Matrix n -> Matrix n
mx_invert m =
    let Matrix xx yx xy yy _ _ = m
        d = xx*yy - yx*xy
    in mx_scalar_multiply (recip d) (mx_adjoint m)

mx_list :: Matrix n -> [n]
mx_list (Matrix a b c d e f) = [a,b,c,d,e,f]

-- | Apply a transformation matrix to a point.
pt_transform :: Num n => Matrix n -> Pt n -> Pt n
pt_transform (Matrix a b c d e f) (Pt x y) =
    let x' = x * a + y * c + e
        y' = x * b + y * d + f
    in Pt x' y'

-- * Bezier functions.

bezier3 :: Num n => Pt n -> Pt n -> Pt n -> n -> Pt n
bezier3 (Pt x1 y1) (Pt x2 y2) (Pt x3 y3) mu = (Pt x y)
    where a = mu*mu
          b = 1 - mu
          c = b*b
          x = x1*c + 2*x2*b*mu + x3*a
          y = y1*c + 2*y2*b*mu + y3*a

-- | Four-point bezier curve interpolation.  The index /mu/ is
--   in the range zero to one.
bezier4 :: Num n => Pt n -> Pt n -> Pt n -> Pt n -> n -> Pt n
bezier4 (Pt x1 y1) (Pt x2 y2) (Pt x3 y3) (Pt x4 y4) mu =
    let a = 1 - mu
        b = a*a*a
        c = mu*mu*mu
        x = b*x1 + 3*mu*a*a*x2 + 3*mu*mu*a*x3 + c*x4
        y = b*y1 + 3*mu*a*a*y2 + 3*mu*mu*a*y3 + c*y4
    in Pt x y

-- * Ord

-- | Given /left/ and /right/, is /x/ in range (inclusive).
--
-- > map (in_range 0 1) [-1,0,1,2] == [False,True,True,False]
in_range :: Ord a => a -> a -> a -> Bool
in_range l r x = l <= x && x <= r

-- * List

-- | Split list at element where predicate /f/ over adjacent elements
-- first holds.
--
-- > split_f (\p q -> q - p < 3) [1,2,4,7,11] == ([1,2,4],[7,11])
split_f :: (a -> a -> Bool) -> [a] -> ([a],[a])
split_f f =
    let go i [] = (reverse i,[])
        go i [p] = (reverse (p:i), [])
        go i (p:q:r) =
            if f p q
            then go (p:i) (q:r)
            else (reverse (p:i),q:r)
    in go []

-- | Variant on 'split_f' that segments input.
--
-- > segment_f (\p q -> abs (q - p) < 3) [1,3,7,9,15] == [[1,3],[7,9],[15]]
segment_f :: (a -> a -> Bool) -> [a] -> [[a]]
segment_f f xs =
    let (p,q) = split_f f xs
    in if null q
       then [p]
       else p : segment_f f q

-- | Delete elements of a list using a predicate over the
-- previous and current elements.
delete_f :: (a -> a -> Bool) -> [a] -> [a]
delete_f f =
    let go [] = []
        go [p] = [p]
        go (p:q:r) =
            if f p q
            then go (p:r)
            else p : go (q:r)
    in go

-- | All adjacent pairs of a list.
--
-- > pairs [1..5] == [(1,2),(2,3),(3,4),(4,5)]
pairs :: [x] -> [(x,x)]
pairs l =
    case l of
      x:y:z -> (x,y) : pairs (y:z)
      _ -> []