{-# OPTIONS_HADDOCK ignore-exports #-} {-# LANGUAGE ViewPatterns #-} module Algebra.Vector where import Control.Applicative import Control.Arrow ((***)) import Data.List (sortBy) import Diagrams.Coordinates import Diagrams.TwoD.Types import GHC.Float import MyPrelude data Alignment = CW | CCW | CL deriving (Eq) -- |Convert two dimensions such as (xmin, xmax) and (ymin, ymax) -- to proper square coordinates, as in: -- ((xmin, ymin), (xmax, ymax)) dimToSquare :: (Double, Double) -- ^ x dimension -> (Double, Double) -- ^ y dimension -> ((Double, Double), (Double, Double)) -- ^ square describing those dimensions dimToSquare (x1, x2) (y1, y2) = ((x1, y1), (x2, y2)) -- |Checks whether the Point is in a given Square. inRange :: ((Double, Double), (Double, Double)) -- ^ the square: ((xmin, ymin), (xmax, ymax)) -> P2 Double -- ^ Coordinate -> Bool -- ^ result inRange ((xmin, ymin), (xmax, ymax)) (coords -> x :& y) = x >= min xmin xmax && x <= max xmin xmax && y >= min ymin ymax && y <= max ymin ymax -- |Get the angle between two vectors. getAngle :: V2 Double -> V2 Double -> Double getAngle a b = acos . flip (/) (vecLength a * vecLength b) . scalarProd a $ b -- |Get the length of a vector. vecLength :: V2 Double -> Double vecLength v = sqrt (x^(2 :: Int) + y^(2 :: Int)) where (x, y) = unr2 v -- |Compute the scalar product of two vectors. scalarProd :: V2 Double -> V2 Double -> Double scalarProd (V2 a1 a2) (V2 b1 b2) = a1 * b1 + a2 * b2 -- |Multiply a scalar with a vector. scalarMul :: Double -> V2 Double -> V2 Double scalarMul d (V2 a b) = V2 (a * d) (b * d) -- |Construct a vector that points to a point from the origin. pt2Vec :: P2 Double -> V2 Double pt2Vec = r2 . unp2 -- |Give the point which is at the coordinates the vector -- points to from the origin. vec2Pt :: V2 Double -> P2 Double vec2Pt = p2 . unr2 -- |Construct a vector between two points. vp2 :: P2 Double -- ^ vector origin -> P2 Double -- ^ vector points here -> V2 Double vp2 a b = pt2Vec b - pt2Vec a -- |Computes the determinant of 3 points. det :: P2 Double -> P2 Double -> P2 Double -> Double det (coords -> ax :& ay) (coords -> bx :& by) (coords -> cx :& cy) = (bx - ax) * (cy - ay) - (by - ay) * (cx - ax) -- |Get the point where two lines intesect, if any. Excludes the -- case of end-points intersecting. intersectSeg :: (P2 Double, P2 Double) -> (P2 Double, P2 Double) -> Maybe (P2 Double) intersectSeg (a, b) (c, d) = case intersectSegSeg a b c d of Just x -> if x `notElem` [a,b,c,d] then Just a else Nothing Nothing -> Nothing -- |Get the orientation of 3 points which can either be -- * clock-wise -- * counter-clock-wise -- * collinear getOrient :: P2 Double -> P2 Double -> P2 Double -> Alignment getOrient a b c = case compare (det a b c) 0 of LT -> CW GT -> CCW EQ -> CL --- |Checks if 3 points a,b,c do not build a clockwise triangle by --- connecting a-b-c. This is done by computing the determinant and --- checking the algebraic sign. notcw :: P2 Double -> P2 Double -> P2 Double -> Bool notcw a b c = case getOrient a b c of CW -> False _ -> True --- |Checks if 3 points a,b,c do build a clockwise triangle by --- connecting a-b-c. This is done by computing the determinant and --- checking the algebraic sign. cw :: P2 Double -> P2 Double -> P2 Double -> Bool cw a b c = not . notcw a b $ c -- |Sort X and Y coordinates lexicographically. sortedXY :: [P2 Double] -> [P2 Double] sortedXY = fmap p2 . sortLex . fmap unp2 -- |Sort Y and X coordinates lexicographically. sortedYX :: [P2 Double] -> [P2 Double] sortedYX = fmap p2 . sortLexSwapped . fmap unp2 -- |Sort all points according to their X-coordinates only. sortedX :: [P2 Double] -> [P2 Double] sortedX xs = fmap p2 . sortBy (\(a1, _) (a2, _) -> compare a1 a2) $ fmap unp2 xs -- |Sort all points according to their Y-coordinates only. sortedY :: [P2 Double] -> [P2 Double] sortedY xs = fmap p2 . sortBy (\(_, b1) (_, b2) -> compare b1 b2) $ fmap unp2 xs -- |Apply a function on the coordinates of a point. onPT :: ((Double, Double) -> (Double, Double)) -> P2 Double -> P2 Double onPT f = p2 . f . unp2 -- |Compare the y-coordinate of two points. ptCmpY :: P2 Double -> P2 Double -> Ordering ptCmpY (coords -> _ :& y1) (coords -> _ :& y2) = compare y1 y2 -- |Compare the x-coordinate of two points. ptCmpX :: P2 Double -> P2 Double -> Ordering ptCmpX (coords -> x1 :& _) (coords -> x2 :& _) = compare x1 x2 posInfPT :: P2 Double posInfPT = p2 (read "Infinity", read "Infinity") negInfPT :: P2 Double negInfPT = p2 (negate . read $ "Infinity", negate . read $ "Infinity") -- | Given an infinite line which intersects P1 and P2, -- let P4 be the point on the line that is closest to P3. -- -- Return an indication of where on the line P4 is relative to P1 and P2. -- -- @ -- if P4 == P1 then 0 -- if P4 == P2 then 1 -- if P4 is halfway between P1 and P2 then 0.5 -- @ -- -- @ -- | -- P1 -- | -- P4 +---- P3 -- | -- P2 -- | -- @ -- {-# INLINE closestPointOnLineParam #-} closestPointOnLineParam :: P2 Double -- ^ `P1` -> P2 Double -- ^ `P2` -> P2 Double -- ^ `P3` -> Double closestPointOnLineParam p1 p2 p3 = pt2Vec (p3 - p1) `scalarProd` pt2Vec (p2 - p1) / pt2Vec (p2 - p1) `scalarProd` pt2Vec (p2 - p1) -- | Given four points specifying two lines, get the point where the two lines -- cross, if any. Note that the lines extend off to infinity, so the -- intersection point might not line between either of the two pairs of points. -- -- @ -- \\ / -- P1 P4 -- \\ / -- + -- / \\ -- P3 P2 -- / \\ -- @ -- intersectLineLine :: P2 Double -- ^ `P1` -> P2 Double -- ^ `P2` -> P2 Double -- ^ `P3` -> P2 Double -- ^ `P4` -> Maybe (P2 Double) intersectLineLine (coords -> x1 :& y1) (coords -> x2 :& y2) (coords -> x3 :& y3) (coords -> x4 :& y4) = let dx12 = x1 - x2 dx34 = x3 - x4 dy12 = y1 - y2 dy34 = y3 - y4 den = dx12 * dy34 - dy12 * dx34 in if den == 0 then Nothing else let det12 = x1*y2 - y1*x2 det34 = x3*y4 - y3*x4 numx = det12 * dx34 - dx12 * det34 numy = det12 * dy34 - dy12 * det34 in Just $ p2 (numx / den, numy / den) -- | Get the point where a segment @P1-P2@ crosses another segement @P3-P4@, -- if any. intersectSegSeg :: P2 Double -- ^ `P1` -> P2 Double -- ^ `P2` -> P2 Double -- ^ `P3` -> P2 Double -- ^ `P4` -> Maybe (P2 Double) intersectSegSeg p1 p2 p3 p4 -- TODO: merge closest point checks with intersection, reuse subterms. | Just p0 <- intersectLineLine p1 p2 p3 p4 , t12 <- closestPointOnLineParam p1 p2 p0 , t23 <- closestPointOnLineParam p3 p4 p0 , t12 >= 0 && t12 <= 1 , t23 >= 0 && t23 <= 1 = Just p0 | otherwise = Nothing