hcg-minus-0.12: Data/CG/Minus.hs
-- | CG library (minus).
module Data.CG.Minus where
import Control.Applicative
import Data.Complex
import Data.Maybe
import qualified Data.SG as G {- SG -}
import Text.Printf
-- * Types
-- | Two-dimensional point.
type Pt = G.Point2'
-- type Pt a = Pt {pt_x :: a, pt_y :: a}
-- | Two-dimensional vector.
type Vc = G.Rel2'
-- type Vc a = Vc {vc_x :: a, vc_y :: a}
-- | Two-dimensional line.
type Ln = G.Line2'
-- type Ln a = Ln (Pt a) (Pt a)
-- | Line segments.
type Ls a = [Pt a]
-- | Window, given by a /lower left/ 'Pt' and an /extent/ 'Vc'.
type Wn a = (Pt a,Vc a)
-- | Real number, synonym for 'Double'.
type R = Double
-- * R(eal) functions
-- | 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
-- * Pt functions
-- | 'Pt' constructor.
--
-- > pt_xy (pt 0 pi) == (0,pi)
pt :: a -> a -> Pt a
pt = curry G.Point2
-- pt = Pt
-- | Variant 'Pt' constructor, ie. 'uncurry' 'pt'.
--
-- > pt_xy (pt' (0,pi)) == (0,pi)
pt' :: (a,a) -> Pt a
pt' = G.Point2
--pt' = uncurry pt
-- | /x/ field of 'Pt'.
--
-- > pt_x (pt 0 pi) == 0
pt_x :: Pt t -> t
pt_x = G.getX
-- | /y/ field of 'Pt'.
--
-- > pt_y (pt 0 pi) == pi
pt_y :: Pt t -> t
pt_y = G.getY
-- | /x/ and /y/ fields of 'Pt'.
--
-- > pt_xy (pt 0 pi) == (0,pi)
pt_xy :: Pt a -> (a,a)
pt_xy p = (pt_x p,pt_y p)
-- | 'Pt' of (0,0).
--
-- > pt_origin == pt 0 0
pt_origin :: Num a => Pt a
pt_origin = pt 0 0
-- | Binary operator at 'Pt'. Given the 'Applicative' instance for
-- 'Pt' this is a synonym for 'liftA2'.
--
-- > pt_binary_op (+) (pt 1 1) (pt 2 2) == pt 3 3
pt_binary_op_ :: (a -> b -> c) -> Pt a -> Pt b -> Pt c
pt_binary_op_ f p1 p2 = pt (pt_x p1 `f` pt_x p2) (pt_y p1 `f` pt_y p2)
-- | Variant applicative definition as 'liftA2'
pt_binary_op :: (a -> b -> c) -> Pt a -> Pt b -> Pt c
pt_binary_op = liftA2
-- | Pointwise '+'.
--
-- > pt_add (pt 1 2) (pt 3 4) == pt 4 6
pt_add :: (Num a) => Pt a -> Pt a -> Pt a
pt_add = pt_binary_op (+)
-- | Pointwise '+' (applicative definition).
--
-- > pt_add_ (pt 1 2) (pt 3 4) == pt 4 6
-- > (liftA2 (+)) (pt 1 2) (pt 3 4) == pt 4 6
-- > (pure (+) <*> pt 1 2 <*> pt 3 4) == pt 4 6
pt_add_ :: (Num a) => Pt a -> Pt a -> Pt a
pt_add_ = pt_binary_op_ (+)
-- | Pointwise '-'.
pt_sub :: (Num a) => Pt a -> Pt a -> Pt a
pt_sub = pt_binary_op (-)
-- | Pointwise '*'.
pt_mul :: (Num a) => Pt a -> Pt a -> Pt a
pt_mul = pt_binary_op (*)
-- | Unary operator at 'Pt'. Given 'Applicative' instance for 'Pt'
-- this is a synonym for 'liftA'.
--
-- > pt_unary_op negate (pt 0 1) == pt 0 (-1)
-- > pt_unary_op_ negate (pt 0 1) == pt 0 (-1)
-- > (liftA negate) (pt 0 1) == pt 0 (-1)
-- > (pure negate <*> pt 0 1) == pt 0 (-1)
pt_unary_op_ :: (a -> b) -> Pt a -> Pt b
pt_unary_op_ f p = pt (f (pt_x p)) (f (pt_y p))
-- | Variant applicative definition as 'liftA'.
pt_unary_op :: (a -> b) -> Pt a -> Pt b
pt_unary_op = liftA
-- | Pointwise 'negate'.
--
-- > pt_negate (pt 0 1) == pt 0 (-1)
pt_negate :: (Num a) => Pt a -> Pt a
pt_negate = pt_unary_op negate
-- | Pointwise 'abs'.
--
-- > pt_abs (pt (-1) 1) == pt 1 1
pt_abs :: (Num a) => Pt a -> Pt a
pt_abs = pt_unary_op abs
-- | Pointwise 'signum'.
--
-- > pt_signum (pt (-1/2) (1/2)) == pt (-1) 1
pt_signum :: (Num a) => Pt a -> Pt a
pt_signum = pt_unary_op signum
-- | '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_unary_op f
-- | Swap /x/ and /y/ coordinates at 'Pt'.
--
-- > pt_swap (pt 1 2) == pt 2 1
pt_swap :: Pt a -> Pt a
pt_swap p = pt (pt_y p) (pt_x p)
-- | 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 p = pt (pt_x p) (negate (pt_y p))
-- | '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_unary_op 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 p = pt' (polar (pt_x p :+ pt_y p))
-- | 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 p =
let (mg,ph) = pt_xy p
c = mkPolar mg ph
x = realPart c
y = imagPart c
in pt x y
-- | Scalar 'Pt' '+'.
--
-- > pt_offset 1 pt_origin == pt 1 1
pt_offset :: Num a => a -> Pt a -> Pt a
pt_offset = pt_unary_op . (+)
-- | Pointwise 'min'.
pt_min :: (Ord a) => Pt a -> Pt a -> Pt a
pt_min = pt_binary_op min
-- | Pointwise 'max'.
pt_max :: (Ord a) => Pt a -> Pt a -> Pt a
pt_max = pt_binary_op 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 p0 p1 p2 =
let (x0,y0) = pt_xy p0
(x1,y1) = pt_xy p1
(x2,y2) = pt_xy p2
in 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 p = atan2 (pt_y p) (pt_x p)
-- | 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 `pt_sub` p)
-- | Pointwise '+'.
pt_translate :: (Num a,Eq a) => Pt a -> Vc a -> Pt a
pt_translate = G.plusDir
-- pt_translate p v = pt (pt_x p + vc_x v) (pt_y p + vc_y v)
-- pt_translate (Pt x y) (Vc dx dy) = Pt (x + dx) (y + dy)
-- | Alternate implementation of 'pt_translate'.
pt_translate_ :: Num a => Pt a -> Vc a -> Pt a
pt_translate_ p v =
let (dx,dy) = vc_xy v
(x,y) = pt_xy p
in pt (x+dx) (y+dy)
-- | 'pt_unary_op' 'fromIntegral'.
pt_from_i :: (Integral i,Num a) => Pt i -> Pt a
pt_from_i = pt_unary_op fromIntegral
-- | 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 = G.distFrom
-- | 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 p =
let (x,y) = pt_xy p
in 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 p =
let (x,y) = pt_xy p
s = sin a
c = cos a
in pt (x * c - y * s) (y * c + x * s)
-- * Vc functions
-- | Construct 'Vc'.
vc :: Num a => a -> a -> Vc a
vc = curry G.makeRel2
-- vc = Vc
-- | Alernate constructor for 'Vc'.
--
-- > (vc 0 0 == vc' (0,0)) == True
vc' :: Num a => (a,a) -> Vc a
vc' = G.makeRel2
-- vc' = uncurry vc
-- | 'Vc' /x/ field.
vc_x :: Vc t -> t
vc_x = G.getX
-- | 'Vc' /y/ field.
vc_y :: Vc t -> t
vc_y = G.getY
-- | 'Vc' /x/ and /y/ fields.
vc_xy :: Vc a -> (a,a)
vc_xy v = (vc_x v,vc_y v)
-- | 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 = G.scaleRel
-- vc_scale n v = vc (vc_x v * n) (vc_y v * n)
-- | 'Vc' dot product.
--
-- > vc_dot (vc 1 2) (vc 3 4) == 11
vc_dot :: Num a => Vc a -> Vc a -> a
vc_dot = G.dotProduct
-- vc_dot p q = (vc_x p * vc_x q) + (vc_y p * vc_y q)
-- | 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 = G.unitVector
-- | 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 v1 v2 =
let (x1,y1) = vc_xy v1
(x2,y2) = vc_xy v2
t1 = atan2 y1 x1
t2 = atan2 y2 x2
in r_constrain (-pi,pi) (t2 - t1)
-- * Line functions
-- | 'Ln' constructor.
--
-- > ln_start (ln (pt 0 0) (pt 1 1)) == pt 0 0
ln :: (Num a,Eq a) => Pt a -> Pt a -> Ln a
ln = G.lineTo
-- | Variant constructor.
--
-- > ln_start (ln_ (pt 1 1) (pt 0 0)) == pt 1 1
ln_ :: (Num a,Eq a) => Pt a -> Pt a -> Ln a
ln_ p q = G.Line2 p (G.fromPt q p)
-- | Variant on 'ln' which takes 'Pt' co-ordinates as duples.
--
-- > ln' (0,0) (1,1) == ln (pt 0 0) (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)
-- | Initial 'Pt' of 'Ln'.
--
-- > ln_start (ln (pt 0 0) (pt 1 1)) == pt 0 0
ln_start :: (Num a,Eq a) => Ln a -> Pt a
ln_start = G.getLineStart
-- | Alternate implementation of 'ln_start' (without 'Num' constraint).
ln_start_ :: Ln a -> Pt a
ln_start_ (G.Line2 p _) = p
-- | End 'Pt' of 'Ln'.
--
-- > ln_end (ln (pt 0 0) (pt 1 1)) == pt 1 1
ln_end :: (Num a,Eq a) => Ln a -> Pt a
ln_end = G.getLineEnd
-- | Alternate implementation of 'ln_end'.
ln_end_ :: (Num a,Eq a) => Ln a -> Pt a
ln_end_ (G.Line2 p v) = p `pt_translate` v
-- | '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 = G.getLineDir
-- | Alternate implementation of 'ln_vc', without 'Num' constraint.
ln_vc_ :: Ln a -> Vc a
ln_vc_ (G.Line2 _ v) = v
-- | 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 = G.toAngle . ln_vc
-- | 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 l = (ln_start l,ln_end l)
-- | 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' l =
let (p1,p2) = ln_pt l
in (pt_xy p1,pt_xy p2)
-- | 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 l =
let (p1,p2) = ln_pt l
x = (pt_x p1 + pt_x p2) / 2
y = (pt_y p1 + pt_y p2) / 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 = G.mag . ln_vc
-- | Variant definition of 'ln_magnitude'.
ln_magnitude_ :: Ln R -> R
ln_magnitude_ l =
let ((x1,y1),(x2,y2)) = ln_pt' l
x = x2 - x1
y = y2 - y1
in sqrt (x * x + y * y)
-- | 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 l =
let (p,q) = ln_pt l
in case compare (pt_x p) (pt_x q) of
LT -> l
EQ -> if pt_y p <= pt_y q then l 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 = G.makeLength
-- | 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
-- | Does 'Pt' /p/ lie on 'Ln' (inclusive).
--
-- > let f = pt_on_line (ln (pt 0 0) (pt 1 1))
-- > in map f [pt 0.5 0.5,pt 2 2,pt (-1) (-1),pt 0 0] == [True,False,False,True]
pt_on_line :: Ln R -> Pt R -> Bool
pt_on_line l r =
let (p,q) = ln_pt l
(i,j) = pt_xy (pt_to_polar (q `pt_sub` p))
(i',j') = pt_xy (pt_to_polar (r `pt_sub` p))
in r == p || r == q || (j == j' && i' <= i)
-- | Variant definition of 'pt_on_line', exclusive of starting point.
--
-- > let f = pt_on_line_ (ln (pt 0 0) (pt 1 1))
-- > in map f [pt 0.5 0.5,pt 2 2,pt (-1) (-1),pt 0 0] == [True,False,False,False]
pt_on_line_ :: Ln R -> Pt R -> Bool
pt_on_line_ l p =
case G.distAlongLine p l of
Nothing -> False
Just d -> d >= 0 && d <= 1
-- | Calculate the point that extends a line by length 'n'.
--
-- > pt_linear_extension (sqrt 2) (ln (pt 0 0) (pt 1 1)) ~= pt 2 2
pt_linear_extension :: R -> Ln R -> Pt R
pt_linear_extension n l =
let (p0,p1) = ln_pt l
(mg,ph) = pt_xy (pt_to_polar (p1 `pt_sub` p0))
in pt_from_polar (pt (mg+n) ph) `pt_add` p0
-- * Intersection
-- | 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
-- | 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 G.intersectLines2 l0 l1 of
Nothing -> Nothing
Just (i,j) -> if i >= 0 && i <= 1 && j >= 0 && j <= 1
then Just (G.alongLine 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 (ln' (1,1) (3,8)) (ln' (0.5,2) (4,7)) == True
-- > ln_intersect (ln' (3.5,9) (3.5,0.5)) (ln' (3,1) (9,1)) == True
ln_intersect :: (Ord a, Fractional a) => Ln a -> Ln a -> Bool
ln_intersect 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)
-- | 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 `G.sameDirection` 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) `G.sameDirection` 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
-- * 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 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
-- | 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
-- | 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 pt' [(0,0),(1,1),(3,3)])) == [2,1]
ls_separate :: (Ord a,Num a) => Vc a -> Ls a -> [Ls a]
ls_separate v =
let (dx,dy) = vc_xy v
f p0 p1 = let (x0,y0) = pt_xy p0
(x1,y1) = pt_xy p1
in abs (x1 - x0) < dx &&
abs (y1 - y0) < dy
in segment_f f
-- | 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
-- | 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 v =
let (x,y) = vc_xy v
too_close p0 p1 =
let (x0,y0) = pt_xy p0
(x1,y1) = pt_xy p1
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'
-- | 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)
_ -> []
-- | 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 p =
let (x,y) = pt_xy p
in case s of
[] -> undefined
l0:l -> let xs = pairs ((l0:l)++[l0])
f (p1,p2) =
let (x1,y1) = pt_xy p1
(x2,y2) = pt_xy p2
in 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
-- * Window
-- | 'Wn' constructor.
wn :: Pt a -> Vc a -> Wn a
wn p v = (p,v)
-- | Variant 'Wn' constructor.
wn' :: Num a => (a,a) -> (a,a) -> Wn a
wn' (x,y) (i,j) = (pt x y,vc i j)
-- | Lower-left 'Pt' of 'Wn'.
wn_ll :: Wn a -> Pt a
wn_ll (p,_) = p
-- | Extent 'Vc' of 'Wn'.
wn_ex :: Wn a -> Vc a
wn_ex (_,v) = v
-- | 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 (p,v) = (pt_xy p,vc_xy v)
-- | 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 w =
let ((x0,y0),(dx,dy)) = wn_extract w
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 w p =
let ((lx,ly),(dx,dy)) = wn_extract w
(x,y) = pt_xy p
(ux,uy) = (lx+dx,ly+dy)
in x > lx && x < ux && y > ly && y < uy
-- | '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 l =
let (p0,p1) = ls_minmax l
(x0,y0) = pt_xy p0
(x1,y1) = pt_xy p1
in (pt x0 y0,vc (x1-x0) (y1-y0))
-- | A 'Wn' that encompasses both input 'Wn's.
wn_join :: (Num a,Ord a) => Wn a -> Wn a -> Wn a
wn_join w0 w1 =
let ((x0,y0),(dx0,dy0)) = wn_extract w0
((x1,y1),(dx1,dy1)) = wn_extract w1
x = min x0 x1
y = min y0 y1
dx = max (x0+dx0) (x1+dx1) - x
dy = max (y0+dy0) (y1+dy1) - y
in (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 ((x0,y0),(dx0,dy0)) = wn_extract w0
((x1,y1),(dx1,dy1)) = wn_extract 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
-- | Given a 'Wn' for a 'Ls', normalise the 'Ls' to lie within (0,1).
ls_normalise_w :: Wn R -> Ls R -> Ls R
ls_normalise_w w =
let ((x0,y0),(dx,dy)) = wn_extract w
z = max dx dy
f p = let (x,y) = pt_xy p
in pt ((x - x0) / z) ((y - y0) / z)
in map f
-- | Shift lower left 'Pt' of 'Wn' by indicated 'Pt' (ie. 'pt_add').
pt_shift_w :: Num a => Pt a -> Wn a -> Wn a
pt_shift_w p (dp,ex) = (p `pt_add` dp,ex)
-- | Negate /y/ field of lower left 'Pt' of 'Wn'.
wn_negate_y :: Num a => Wn a -> Wn a
wn_negate_y (p,v) = (pt_negate_y p,v)