diff --git a/Data/CG/Minus.hs b/Data/CG/Minus.hs --- a/Data/CG/Minus.hs +++ b/Data/CG/Minus.hs @@ -1,37 +1,67 @@ -- | 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 +import Data.Complex {- base -} +import Data.Maybe {- base -} +import Text.Printf {- base -} -- * Types -- | Two-dimensional point. -type Pt = G.Point2' --- type Pt a = Pt {pt_x :: a, pt_y :: a} +-- +-- 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) --- | Two-dimensional vector. -type Vc = G.Rel2' --- type Vc a = Vc {vc_x :: a, vc_y :: a} +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. -type Ln = G.Line2' --- type Ln a = Ln (Pt a) (Pt a) +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'. -type Wn a = (Pt a,Vc a) +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] @@ -54,190 +84,118 @@ 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 +-- | Sum of squares. +mag_sq :: Num a => a -> a -> a +mag_sq x y = x * x + y * y --- | Variant 'Pt' constructor, ie. 'uncurry' 'pt'. --- --- > pt_xy (pt' (0,pi)) == (0,pi) -pt' :: (a,a) -> Pt a -pt' = G.Point2 ---pt' = uncurry pt +-- | 'sqrt' of 'mag_sq'. +mag :: Floating c => c -> c -> c +mag x = sqrt . mag_sq x --- | /x/ field of 'Pt'. --- --- > pt_x (pt 0 pi) == 0 -pt_x :: Pt t -> t -pt_x = G.getX +-- * Pt functions --- | /y/ field of 'Pt'. --- --- > pt_y (pt 0 pi) == pi -pt_y :: Pt t -> t -pt_y = G.getY +-- | Tuple constructor. +pt' :: (a,a) -> Pt a +pt' (x,y) = Pt x y --- | /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) +-- | 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 == 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 +pt_origin = Pt 0 0 --- | Pointwise 'abs'. --- --- > pt_abs (pt (-1) 1) == pt 1 1 -pt_abs :: (Num a) => Pt a -> Pt a -pt_abs = pt_unary_op abs +-- | 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) --- | 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 +-- | 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 1 == Pt 1 1 pt_from_scalar :: (Num a) => a -> Pt a -pt_from_scalar a = pt a 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 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 + in pt_uop f -- | Swap /x/ and /y/ coordinates at 'Pt'. -- --- > pt_swap (pt 1 2) == pt 2 1 +-- > pt_swap (Pt 1 2) == Pt 2 1 pt_swap :: Pt a -> Pt a -pt_swap p = pt (pt_y p) (pt_x p) +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 (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_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 90 270) == Pt (pi/2) (pi*(3/2)) pt_to_radians :: Pt R -> Pt R -pt_to_radians = pt_unary_op r_to_radians +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 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)) +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 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 +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 1 pt_origin == Pt 1 1 pt_offset :: Num a => a -> Pt a -> Pt a -pt_offset = pt_unary_op . (+) +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_binary_op min +pt_min = pt_binop min -- | Pointwise 'max'. pt_max :: (Ord a) => Pt a -> Pt a -> Pt a -pt_max = pt_binary_op max +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 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 +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 (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)) + (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'. @@ -246,172 +204,137 @@ 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 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 0 1) == pi / 2 pt_angle_o :: Pt R -> R -pt_angle_o p = atan2 (pt_y p) (pt_x p) +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 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) +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 = 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_translate (Pt x y) (Vc dx dy) = Pt (x + dx) (y + dy) --- | 'pt_unary_op' 'fromIntegral'. +-- | 'pt_uop' 'fromIntegral'. pt_from_i :: (Integral i,Num a) => Pt i -> Pt a -pt_from_i = pt_unary_op fromIntegral +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 (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 +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] +-- > 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 +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 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 +pt_rotate a (Pt x y) = + let s = sin a c = cos a - in pt (x * c - y * s) (y * c + x * s) + 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 +-- | 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) --- | 'Vc' /x/ field. -vc_x :: Vc t -> t -vc_x = G.getX +-- | 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) --- | 'Vc' /y/ field. -vc_y :: Vc t -> t -vc_y = G.getY +-- | 'mag_sq' of 'Vc'. +vc_mag_sq :: Floating c => Vc c -> c +vc_mag_sq (Vc dx dy) = mag_sq dx dy --- | 'Vc' /x/ and /y/ fields. -vc_xy :: Vc a -> (a,a) -vc_xy v = (vc_x v,vc_y v) +-- | '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 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_scale n (Vc x y) = Vc (x * n) (y * n) -- | 'Vc' dot product. -- --- > vc_dot (vc 1 2) (vc 3 4) == 11 +-- > 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) +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 (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 +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 v1 v2 = - let (x1,y1) = vc_xy v1 - (x2,y2) = vc_xy v2 - t1 = atan2 y1 x1 +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 --- | '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. +-- | Variant on 'Ln' which takes 'Pt' co-ordinates as duples. -- --- > ln' (0,0) (1,1) == ln (pt 0 0) (pt 1 1) +-- > 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) - --- | 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 +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 +-- > 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 +ln_vc (Ln p q) = let Pt x y = q - p in Vc x y --- | Alternate implementation of 'ln_vc', without 'Num' constraint. -ln_vc_ :: Ln a -> Vc a -ln_vc_ (G.Line2 _ v) = v +-- | '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 @@ -419,151 +342,149 @@ -- > 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 +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 l = (ln_start l,ln_end l) +-- > 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' (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) +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 (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 +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 (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 + (Nothing,Nothing) -> Pt 0 0 (Just p,Nothing) -> p (Nothing, Just q) -> q - (Just p, Just q) -> ln_midpoint (ln p 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 (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) +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 (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 +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 (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 +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 (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) +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_ (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)) --- > in map f [pt 0.5 0.5,pt 2 2,pt (-1) (-1),pt 0 0] == [True,False,False,True] +-- > 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 - (i,j) = pt_xy (pt_to_polar (q `pt_sub` p)) - (i',j') = pt_xy (pt_to_polar (r `pt_sub` p)) + 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) --- | 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 +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' (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 (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 + case ln_intersect l0 l1 of Nothing -> Nothing Just (i,j) -> if i >= 0 && i <= 1 && j >= 0 && j <= 1 - then Just (G.alongLine i l0) + then Just (ln_pt_along i l0) else Nothing -- | Variant definition of 'ln_intersection', using algorithm at -- . -- -- > 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_ (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 @@ -574,21 +495,21 @@ in if d == 0 then Nothing else if ua' == 0 && ub' == 0 - then Just (pt x1 y1) + 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) + 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 +-- > 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 @@ -620,18 +541,26 @@ 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 `G.sameDirection` ln_vc q +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) `G.sameDirection` ln_vc (ln_sort q) +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. -- @@ -646,6 +575,18 @@ 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. @@ -654,7 +595,7 @@ -- | Variant 'Ls' constructor from 'Pt' co-ordinates as duples. ls' :: [(a,a)] -> Ls a -ls' = map pt' +ls' = map (uncurry Pt) -- | Negate /y/ elements. ls_negate_y :: (Num a) => Ls a -> Ls a @@ -667,64 +608,22 @@ [] -> 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] +-- > 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 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 +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 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) +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 @@ -736,38 +635,26 @@ 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' [(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))) +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' [(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 @@ -778,75 +665,64 @@ [] -> True p:z -> pt_is_normal p && ls_check_normalised z --- * Window +-- | 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]) --- | 'Wn' constructor. -wn :: Pt a -> Vc a -> Wn a -wn p v = (p,v) +-- * Window -- | 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 +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 (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) +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 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 +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] +-- > 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) +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 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)) +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 w0 w1 = - let ((x0,y0),(dx0,dy0)) = wn_extract w0 - ((x1,y1),(dx1,dy1)) = wn_extract w1 - x = min x0 x1 +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 (pt x y,vc dx dy) + 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 ((x0,y0),(dx0,dy0)) = wn_extract w0 - ((x1,y1),(dx1,dy1)) = wn_extract w1 +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'. @@ -869,19 +745,202 @@ 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). +-- | 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 = - 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 +ls_normalise_w w = map (wn_normalise_f w) --- | Shift lower left 'Pt' of 'Wn' by indicated 'Pt' (ie. 'pt_add'). +-- | 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 (dp,ex) = (p `pt_add` dp,ex) +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 (p,v) = (pt_negate_y p,v) +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) + _ -> [] diff --git a/Data/CG/Minus/Arrow.hs b/Data/CG/Minus/Arrow.hs --- a/Data/CG/Minus/Arrow.hs +++ b/Data/CG/Minus/Arrow.hs @@ -6,7 +6,7 @@ -- | Given the arrow body 'Ln' and the arrow length and arrow angle -- (in radians) 'R' calculate the 'Pt' of each arrow tip. -- --- > arrow_coord (ln' (0,0) (1,1)) 0.1 (pi/9) +-- > arrow_coord (Ln (Pt 0 0) (Pt 1 1)) 0.1 (pi/9) arrow_coord :: Ln R -> R -> R -> (Pt R,Pt R) arrow_coord l n a = let ((x0,y0),(x1,y1)) = ln_pt' l @@ -15,4 +15,4 @@ y2 = y1 + n * sin (a' - a) x3 = x1 + n * cos (a' + a) y3 = y1 + n * sin (a' + a) - in (pt x2 y2,pt x3 y3) + in (Pt x2 y2,Pt x3 y3) diff --git a/Data/CG/Minus/Bearing.hs b/Data/CG/Minus/Bearing.hs --- a/Data/CG/Minus/Bearing.hs +++ b/Data/CG/Minus/Bearing.hs @@ -12,7 +12,7 @@ -- | Bearing from 'Pt' /p/ to /q/. -- --- > let f (x,y) = bearing (pt 0 0) (pt x y) +-- > let f (x,y) = bearing (Pt 0 0) (Pt x y) -- > map f [(0,1),(1,1),(1,0),(1,-1)] == [N,NE,E,SE] -- > map f [(0,-1),(-1,-1),(-1,0),(-1,1)] == [S,SW,W,NW] -- > map f [(1/4,1),(1,1/4),(1,-1/4),(1/4,-1)] == [NNE,ENE,ESE,SSE] @@ -25,7 +25,7 @@ -- | Bearing to nearest eight point compass bearing -- --- > let f (x,y) = bearing_8 (pt 0 0) (pt x y) +-- > let f (x,y) = bearing_8 (Pt 0 0) (Pt x y) -- > map f [(1/4,1),(1,1/4),(1,-1/4),(1/4,-1)] == [N,E,E,S] bearing_8 :: Pt R -> Pt R -> Bearing bearing_8 p q = diff --git a/Data/CG/Minus/Colour.hs b/Data/CG/Minus/Colour.hs --- a/Data/CG/Minus/Colour.hs +++ b/Data/CG/Minus/Colour.hs @@ -2,8 +2,8 @@ module Data.CG.Minus.Colour where import Data.Colour {- colour -} -import Data.Colour.SRGB -import Data.Colour.Names as N +import qualified Data.Colour.SRGB as S {- colour -} +import qualified Data.Colour.Names as N {- colour -} -- | Opaque colour. type C = Colour Double @@ -13,12 +13,12 @@ -- | Grey 'Colour'. mk_grey :: (Ord a,Floating a) => a -> Colour a -mk_grey x = sRGB x x x +mk_grey x = S.sRGB x x x -- | Reduce 'Colour' to grey. Constants are @0.3@, @0.59@ and @0.11@. to_greyscale :: (Ord a,Floating a) => Colour a -> a to_greyscale c = - let (RGB r g b) = toSRGB c + let (S.RGB r g b) = S.toSRGB c in r * 0.3 + g * 0.59 + b * 0.11 -- | 'mk_grey' '.' 'to_greyscale'. @@ -37,13 +37,13 @@ -- | Tuple to 'C', inverse of 'unC'. toC :: (Double,Double,Double) -> C -toC (r,g,b) = sRGB r g b +toC (r,g,b) = S.sRGB r g b -- | 'C' to /(red,green,blue)/ tuple. unC :: C -> (Double,Double,Double) unC x = - let x' = toSRGB x - in (channelRed x', channelGreen x', channelBlue x') + let x' = S.toSRGB x + in (S.channelRed x',S.channelGreen x',S.channelBlue x') -- | Tuple to 'Ca', inverse of 'unCa'. toCa :: (Double,Double,Double,Double) -> Ca @@ -52,42 +52,45 @@ -- | 'Ca' to /(red,green,blue,alpha)/ tuple unCa :: Ca -> (Double,Double,Double,Double) unCa x = - let x' = toSRGB (pureColour x) - in (channelRed x', channelGreen x', channelBlue x', alphaChannel x) + let x' = S.toSRGB (pureColour x) + in (S.channelRed x' + ,S.channelGreen x' + ,S.channelBlue x' + ,alphaChannel x) -- * Constants -- | Venetian red (@#c80815@). venetianRed :: C -venetianRed = sRGB24read "#c80815" +venetianRed = S.sRGB24read "#c80815" -- | Swedish azure blue (@#005b99@). swedishAzureBlue :: C -swedishAzureBlue = sRGB24read "#005b99" +swedishAzureBlue = S.sRGB24read "#005b99" -- | Safety orange (@#ff6600@). safetyOrange :: C -safetyOrange = sRGB24read "#ff6600" +safetyOrange = S.sRGB24read "#ff6600" -- | Dye magenta (@#ca1f7b@). dyeMagenta :: C -dyeMagenta = sRGB24read "#ca1f7b" +dyeMagenta = S.sRGB24read "#ca1f7b" -- | Candlelight yellow (@#fcd116@). candlelightYellow :: C -candlelightYellow = sRGB24read "#fcd116" +candlelightYellow = S.sRGB24read "#fcd116" -- | Subtractive primary cyan (@#00B7EB@). subtractivePrimaryCyan :: C -subtractivePrimaryCyan = sRGB24read "#00B7EB" +subtractivePrimaryCyan = S.sRGB24read "#00B7EB" -- | Fern green (@#009246@). fernGreen :: C -fernGreen = sRGB24read "#009246" +fernGreen = S.sRGB24read "#009246" -- | Sepia brown (@#704214@). sepiaBrown :: C -sepiaBrown = sRGB24read "#704214" +sepiaBrown = S.sRGB24read "#704214" -- | The set of named colours defined in this module. non_svg_colour_set :: [C] @@ -107,150 +110,150 @@ -- alphabetical order). svg_colour_set :: [C] svg_colour_set = - [aliceblue - ,antiquewhite - ,aqua - ,aquamarine - ,azure - ,beige - ,bisque - ,black - ,blanchedalmond - ,blue - ,blueviolet - ,brown - ,burlywood - ,cadetblue - ,chartreuse - ,chocolate - ,coral - ,cornflowerblue - ,cornsilk - ,crimson - ,cyan - ,darkblue - ,darkcyan - ,darkgoldenrod - ,darkgray - ,darkgreen - ,darkgrey - ,darkkhaki - ,darkmagenta - ,darkolivegreen - ,darkorange - ,darkorchid - ,darkred - ,darksalmon - ,darkseagreen - ,darkslateblue - ,darkslategray - ,darkslategrey - ,darkturquoise - ,darkviolet - ,deeppink - ,deepskyblue - ,dimgray - ,dimgrey - ,dodgerblue - ,firebrick - ,floralwhite - ,forestgreen - ,fuchsia - ,gainsboro - ,ghostwhite - ,gold - ,goldenrod - ,gray - ,grey - ,green - ,greenyellow - ,honeydew - ,hotpink - ,indianred - ,indigo - ,ivory - ,khaki - ,lavender - ,lavenderblush - ,lawngreen - ,lemonchiffon - ,lightblue - ,lightcoral - ,lightcyan - ,lightgoldenrodyellow - ,lightgray - ,lightgreen - ,lightgrey - ,lightpink - ,lightsalmon - ,lightseagreen - ,lightskyblue - ,lightslategray - ,lightslategrey - ,lightsteelblue - ,lightyellow - ,lime - ,limegreen - ,linen - ,magenta - ,maroon - ,mediumaquamarine - ,mediumblue - ,mediumorchid - ,mediumpurple - ,mediumseagreen - ,mediumslateblue - ,mediumspringgreen - ,mediumturquoise - ,mediumvioletred - ,midnightblue - ,mintcream - ,mistyrose - ,moccasin - ,navajowhite - ,navy - ,oldlace - ,olive - ,olivedrab - ,orange - ,orangered - ,orchid - ,palegoldenrod - ,palegreen - ,paleturquoise - ,palevioletred - ,papayawhip - ,peachpuff - ,peru - ,pink - ,plum - ,powderblue - ,purple - ,red - ,rosybrown - ,royalblue - ,saddlebrown - ,salmon - ,sandybrown - ,seagreen - ,seashell - ,sienna - ,silver - ,skyblue - ,slateblue - ,slategray - ,slategrey - ,snow - ,springgreen - ,steelblue + [N.aliceblue + ,N.antiquewhite + ,N.aqua + ,N.aquamarine + ,N.azure + ,N.beige + ,N.bisque + ,N.black + ,N.blanchedalmond + ,N.blue + ,N.blueviolet + ,N.brown + ,N.burlywood + ,N.cadetblue + ,N.chartreuse + ,N.chocolate + ,N.coral + ,N.cornflowerblue + ,N.cornsilk + ,N.crimson + ,N.cyan + ,N.darkblue + ,N.darkcyan + ,N.darkgoldenrod + ,N.darkgray + ,N.darkgreen + ,N.darkgrey + ,N.darkkhaki + ,N.darkmagenta + ,N.darkolivegreen + ,N.darkorange + ,N.darkorchid + ,N.darkred + ,N.darksalmon + ,N.darkseagreen + ,N.darkslateblue + ,N.darkslategray + ,N.darkslategrey + ,N.darkturquoise + ,N.darkviolet + ,N.deeppink + ,N.deepskyblue + ,N.dimgray + ,N.dimgrey + ,N.dodgerblue + ,N.firebrick + ,N.floralwhite + ,N.forestgreen + ,N.fuchsia + ,N.gainsboro + ,N.ghostwhite + ,N.gold + ,N.goldenrod + ,N.gray + ,N.grey + ,N.green + ,N.greenyellow + ,N.honeydew + ,N.hotpink + ,N.indianred + ,N.indigo + ,N.ivory + ,N.khaki + ,N.lavender + ,N.lavenderblush + ,N.lawngreen + ,N.lemonchiffon + ,N.lightblue + ,N.lightcoral + ,N.lightcyan + ,N.lightgoldenrodyellow + ,N.lightgray + ,N.lightgreen + ,N.lightgrey + ,N.lightpink + ,N.lightsalmon + ,N.lightseagreen + ,N.lightskyblue + ,N.lightslategray + ,N.lightslategrey + ,N.lightsteelblue + ,N.lightyellow + ,N.lime + ,N.limegreen + ,N.linen + ,N.magenta + ,N.maroon + ,N.mediumaquamarine + ,N.mediumblue + ,N.mediumorchid + ,N.mediumpurple + ,N.mediumseagreen + ,N.mediumslateblue + ,N.mediumspringgreen + ,N.mediumturquoise + ,N.mediumvioletred + ,N.midnightblue + ,N.mintcream + ,N.mistyrose + ,N.moccasin + ,N.navajowhite + ,N.navy + ,N.oldlace + ,N.olive + ,N.olivedrab + ,N.orange + ,N.orangered + ,N.orchid + ,N.palegoldenrod + ,N.palegreen + ,N.paleturquoise + ,N.palevioletred + ,N.papayawhip + ,N.peachpuff + ,N.peru + ,N.pink + ,N.plum + ,N.powderblue + ,N.purple + ,N.red + ,N.rosybrown + ,N.royalblue + ,N.saddlebrown + ,N.salmon + ,N.sandybrown + ,N.seagreen + ,N.seashell + ,N.sienna + ,N.silver + ,N.skyblue + ,N.slateblue + ,N.slategray + ,N.slategrey + ,N.snow + ,N.springgreen + ,N.steelblue ,N.tan - ,teal - ,thistle - ,tomato - ,turquoise - ,violet - ,wheat - ,white - ,whitesmoke - ,yellow - ,yellowgreen] + ,N.teal + ,N.thistle + ,N.tomato + ,N.turquoise + ,N.violet + ,N.wheat + ,N.white + ,N.whitesmoke + ,N.yellow + ,N.yellowgreen] diff --git a/Data/CG/Minus/Colour/Planck.hs b/Data/CG/Minus/Colour/Planck.hs --- a/Data/CG/Minus/Colour/Planck.hs +++ b/Data/CG/Minus/Colour/Planck.hs @@ -1,14 +1,13 @@ -- | Planck radiation equation. module Data.CG.Minus.Colour.Planck where -import Data.CG.Minus (R) import Data.CG.Minus.Colour -- | Given wavelength (in microns) and temperature (in degrees Kelvin) -- solve Planck's radiation equation. -- -- > planck_rad_eq 0.7 2600 == 8.22656629154115e7 -planck_rad_eq :: R -> R -> R +planck_rad_eq :: Floating a => a -> a -> a planck_rad_eq l t = let k0 = 3.7403e10 k1 = -5.0 @@ -28,7 +27,7 @@ -- -- > let h (r,g,b) = let f = floor . (*) 255 in (f r,f g,f b) -- > in map (h . k_to_rgb) [2600,28000] == [(255,95,22),(49,118,254)] -k_to_rgb :: R -> (R,R,R) +k_to_rgb :: (Floating t, Ord t) => t -> (t,t,t) k_to_rgb k = let r = planck_rad_eq 0.7000 k g = planck_rad_eq 0.5461 k @@ -37,5 +36,5 @@ in (r * s,g * s,b * s) -- | 'toC' '.' 'k_to_rgb'. -k_to_colour :: R -> C +k_to_colour :: Double -> C k_to_colour = toC . k_to_rgb diff --git a/README b/README --- a/README +++ b/README @@ -5,7 +5,7 @@ [hs]: http://haskell.org/ -© [rohan drape][rd], 2009-2012, [gpl] +© [rohan drape][rd], 2009-2013, [gpl] [rd]: http://rd.slavepianos.org/ [gpl]: http://gnu.org/copyleft/ diff --git a/hcg-minus.cabal b/hcg-minus.cabal --- a/hcg-minus.cabal +++ b/hcg-minus.cabal @@ -1,10 +1,10 @@ name: hcg-minus -version: 0.12 +version: 0.14 synopsis: haskell cg (minus) description: cg (minus) library license: BSD3 category: Math -copyright: (c) rohan drape, 2011-2012 +copyright: (c) rohan drape, 2011-2013 author: Rohan Drape maintainer: rd@slavepianos.org stability: Experimental @@ -16,7 +16,7 @@ data-files: README library - build-depends: base==4.*,colour,SG + build-depends: base==4.*,colour ghc-options: -Wall -fwarn-tabs exposed-modules: Data.CG.Minus Data.CG.Minus.Arrow