{-# LANGUAGE TypeFamilies #-} {-# OPTIONS -Wall #-} -------------------------------------------------------------------------------- -- | -- Module : Wumpus.Drawing.Paths.Intersection -- Copyright : (c) Stephen Tetley 2011-2012 -- License : BSD3 -- -- Maintainer : Stephen Tetley -- Stability : highly unstable -- Portability : GHC -- -- Intersection of Paths with (infinite) lines. -- -------------------------------------------------------------------------------- module Wumpus.Drawing.Paths.Intersection ( Line(..) , inclinedLine , vectorLine , Ray(..) , inclinedRay , lineLineIntersection , linePathIntersection , linePathSegmentIntersection , rayPathIntersection , rayPathSegmentIntersection , rectangleRadialIntersect , isoscelesTriangleRadialIntersect ) where import Wumpus.Drawing.Basis.BezierCurve import Wumpus.Drawing.Basis.ShapeTrails import Wumpus.Drawing.Paths.Base import Wumpus.Basic.Kernel -- package: wumpus-basic import Wumpus.Core -- package: wumpus-core import Data.AffineSpace -- package: vector-space -------------------------------------------------------------------------------- -- -- Private types - LineEquation, 2x2 matrix and bezier curve. -- -- Although these types are /general/ exposing them just leads -- to a bloated API. -- -- If these types are fund to be more useful - they can go into -- @Drawing.Basis@. -- -- | Line in equational form, i.e. @Ax + By + C = 0@. -- data LineEquation u = LineEquation { _line_eqn_A :: !u , _line_eqn_B :: !u , _line_eqn_C :: !u } deriving (Eq,Show) type instance DUnit (LineEquation u) = u -- | 'lineEquation' : @ point1 * point2 -> LineEquation @ -- -- Construct a line in equational form bisecting the supplied -- points. -- lineEquation :: Num u => Point2 u -> Point2 u -> LineEquation u lineEquation (P2 x1 y1) (P2 x2 y2) = LineEquation a b c where a = y1 - y2 b = x2 - x1 c = (x1*y2) - (x2*y1) -- | 2x2 matrix, considered to be in row-major form. -- -- > (M2'2 a b -- > c d) -- -- data Matrix2'2 u = M2'2 !u !u !u !u deriving (Eq) type instance DUnit (Matrix2'2 u) = u -- | Determinant of a 2x2 matrix. -- det2'2 :: Num u => Matrix2'2 u -> u det2'2 (M2'2 a b c d) = a*d - b*c -------------------------------------------------------------------------------- -- | Infinite line represented by two points. -- data Line u = Line (Point2 u) (Point2 u) deriving (Eq,Show) type instance DUnit (Line u) = u -- | 'inclinedLine' : @ point * ang -> Line @ -- -- | Make an infinite line passing through the supplied point -- inclined by @ang@. -- inclinedLine :: Floating u => Point2 u -> Radian -> Line u inclinedLine radial_ogin ang = Line radial_ogin (radial_ogin .+^ avec ang 100) vectorLine :: Num u => Vec2 u -> Point2 u -> Line u vectorLine v1 p0 = Line p0 (p0 .+^ v1) -- | A 'Ray' extends from the first point, through the second to -- infinity. -- -- ('Line' extends to infinity in both directions. -- data Ray u = Ray (Point2 u) (Point2 u) deriving (Eq,Show) type instance DUnit (Ray u) = u -- | Make an infinite ray starting from the supplied point -- inclined by @ang@. -- inclinedRay :: Floating u => Point2 u -> Radian -> Ray u inclinedRay ray_ogin ang = Ray ray_ogin (ray_ogin .+^ avec ang 100) pointOnLineSeg :: (Real u, Floating u, Ord u, Tolerance u) => Point2 u -> (Point2 u, Point2 u) -> Bool pointOnLineSeg pt (p0,p1) | pt == p0 || pt == p1 = True | otherwise = vdirection v1 == vdirection v0 && vlength v1 `tLTE` vlength v0 where v0 = pvec p0 p1 v1 = pvec p0 pt -- | 'interLineLine' : @ line1 * line2 -> Maybe Point @ -- -- Find the intersection of two lines, if there is one. -- -- Lines are infinite they are represented by points on them, -- they are not line segments. -- -- An answer of @Nothing@ may indicate either the lines coincide -- or the are parallel. -- lineLineIntersection :: (Fractional u, Ord u, Tolerance u) => Line u -> Line u -> Maybe (Point2 u) lineLineIntersection (Line p1 p2) (Line q1 q2) = if det_co `tEQ` 0 then Nothing else Just $ P2 (det_xm / det_co) (det_ym / det_co) where -- Ax + By + C = 0 LineEquation a1 b1 c1 = lineEquation p1 p2 LineEquation a2 b2 c2 = lineEquation q1 q2 coeffM = M2'2 a1 b1 a2 b2 det_co = det2'2 coeffM xM = M2'2 (negate c1) b1 (negate c2) b2 det_xm = det2'2 xM yM = M2'2 a1 (negate c1) a2 (negate c2) det_ym = det2'2 yM linePathIntersection :: (Real u, Floating u, Ord u, Tolerance u) => Line u -> AbsPath u -> Maybe (Point2 u) linePathIntersection ln = step . pathViewL where step EmptyPathL = Nothing step (a :<< bs) = let ans = linePathSegmentIntersection ln a in case ans of Nothing -> step (pathViewL bs) _ -> ans linePathSegmentIntersection :: (Real u, Floating u, Ord u, Tolerance u) => Line u -> PathSegment u -> Maybe (Point2 u) linePathSegmentIntersection ln1 (LineSeg _ p0 p1) = mbWithin p0 p1 $ lineLineIntersection ln1 (Line p0 p1) linePathSegmentIntersection (Line pa pb) (CurveSeg _ p0 p1 p2 p3) = lineEqnCurveIntersection (lineEquation pa pb) (BezierCurve p0 p1 p2 p3) mbWithin :: (Real u, Floating u, Ord u, Tolerance u) => Point2 u -> Point2 u -> Maybe (Point2 u) -> Maybe (Point2 u) mbWithin p0 p1 mb = mb >>= \pt -> if pointOnLineSeg pt (p0,p1) then Just pt else Nothing lineEqnCurveIntersection :: (Floating u, Ord u, Tolerance u) => LineEquation u -> BezierCurve u -> Maybe (Point2 u) lineEqnCurveIntersection eqnline c0 = step c0 where step c = case cut eqnline c of Left pt -> Just pt -- cut at start or end Right False -> Nothing Right True -> let (a,b) = subdivide c in case step a of Just pt -> Just pt Nothing -> step b rayPathIntersection :: (Real u, Floating u, Ord u, Tolerance u) => Ray u -> AbsPath u -> Maybe (Point2 u) rayPathIntersection ry = step . pathViewL where step EmptyPathL = Nothing step (a :<< bs) = let ans = rayPathSegmentIntersection ry a in case ans of Nothing -> step (pathViewL bs) _ -> ans rayPathSegmentIntersection :: (Real u, Floating u, Ord u, Tolerance u) => Ray u -> PathSegment u -> Maybe (Point2 u) rayPathSegmentIntersection (Ray p0 p1) seg = test =<< linePathSegmentIntersection (Line p0 p1) seg where test pt = if vdirection (pvec p0 p1) == vdirection (pvec p0 pt) then Just pt else Nothing -- | Is the curve cut by the line? -- -- The curve might cut at the start or end points - which is good -- as it saves performing a subdivision, but it makes the return -- type a bit involved. -- cut :: (Floating u , Ord u, Tolerance u) => LineEquation u -> BezierCurve u -> Either (Point2 u) Bool cut eqnline (BezierCurve p0 p1 p2 p3) = if d0 `tEQ` 0 then Left p0 else if d3 `tEQ` 0 then Left p3 else let ds = [d0,d1,d2,d3] in Right $ not $ all pve ds || all nve ds where pve = (>= 0) nve = (< 0) d0 = pointLineDistance p0 eqnline d1 = pointLineDistance p1 eqnline d2 = pointLineDistance p2 eqnline d3 = pointLineDistance p3 eqnline -- | 'pointLineDistance' : @ point -> line -> Distance @ -- -- Find the distance from a point to a line in equational form -- using this formula: -- -- > P(u,v) -- > L: Ax + By + C = 0 -- > -- > (A*u) + (B*v) + C -- > ----------------- -- > sqrt $ (A^2) +(B^2) -- -- A positive distance indicates the point is above the line, -- negative indicates below. -- pointLineDistance :: Floating u => Point2 u -> LineEquation u -> u pointLineDistance (P2 u v) (LineEquation a b c) = ((a*u) + (b*v) + c) / base where base = sqrt $ (a^two) + (b^two) two :: Integer two = 2 -------------------------------------------------------------------------------- -- Intersections on common shapes -- | Answer is vector from center. -- rectangleRadialIntersect :: (Real u, Floating u, InterpretUnit u, Tolerance u) => u -> u -> Radian -> Maybe (Vec2 u) rectangleRadialIntersect w h ang = fmap (pvec zeroPt) $ rayPathIntersection (inclinedRay zeroPt ang) rp where rp = anaTrailPath zeroPt $ rectangle_trail w h -- | Answer is vector from centroid. -- isoscelesTriangleRadialIntersect :: (Real u, Floating u , InterpretUnit u, Tolerance u) => u -> u -> Radian -> Maybe (Vec2 u) isoscelesTriangleRadialIntersect bw h ang = fmap (pvec zeroPt) $ rayPathIntersection (inclinedRay zeroPt ang) rp where rp = anaTrailPath zeroPt $ isosceles_triangle_trail bw h