{-# LANGUAGE TypeFamilies               #-}
{-# OPTIONS -Wall #-}

--------------------------------------------------------------------------------
-- |
-- Module      :  Wumpus.Drawing.Paths.Intersection
-- Copyright   :  (c) Stephen Tetley 2011-2012
-- License     :  BSD3
--
-- Maintainer  :  Stephen Tetley <stephen.tetley@gmail.com>
-- 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