{-# OPTIONS -Wall #-}

--------------------------------------------------------------------------------
-- |
-- Module      :  Wumpus.Basic.Geometry.Intersection
-- Copyright   :  (c) Stephen Tetley 2011
-- License     :  BSD3
--
-- Maintainer  :  Stephen Tetley <stephen.tetley@gmail.com>
-- Stability   :  highly unstable
-- Portability :  GHC
--
-- Intersection of line to line and line to plane.
-- 
-- \*\* - WARNING \*\* - this uses quite a high tolerance for 
-- floating point equality. 
--
--------------------------------------------------------------------------------

module Wumpus.Basic.Geometry.Intersection
  ( 

    interLineLine
  , interLinesegLineseg
  , interLinesegLine
  , interCurveLine

  , findIntersect


  ) 
  where

import Wumpus.Basic.Geometry.Base

import Wumpus.Core                              -- package: wumpus-core




-- Potentially lines hould be a new datatype.


-- | '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.
--
interLineLine :: Fractional u 
              => Line u -> Line u -> Maybe (Point2 u)
interLineLine (Line p1 p2) (Line q1 q2) = 
    if det_co == 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



-- | 'interLinesegLineseg' : @ line_segment1 * line_segment2 -> Maybe Point @
-- 
-- Find the intersection of two line segments, if there is one. 
--
-- An answer of @Nothing@ indicates that the line segments 
-- coincide, or that there is no intersection.
--
interLinesegLineseg :: (Fractional u, Ord u, Tolerance u)
                    => LineSegment u -> LineSegment u -> Maybe (Point2 u)
interLinesegLineseg a@(LineSegment p q) b@(LineSegment s t) = 
    interLineLine (Line p q) (Line s t) >>= segcheck
  where
    segcheck = mbCheck (\pt -> withinPoints pt a && withinPoints pt b)
 

-- | 'interLinesegLine' : @ line_segment * line -> Maybe Point @
-- 
-- Find the intersection of a line and a line segment, if there 
-- is one. 
--
-- An answer of @Nothing@ indicates that the the line and line
-- segment coincide, or that there is no intersection.
--
interLinesegLine :: (Fractional u, Ord u, Tolerance u)
                 => LineSegment u -> Line u -> Maybe (Point2 u)
interLinesegLine a@(LineSegment p q) line = 
    interLineLine (Line p q) line >>= segcheck
  where
    segcheck = mbCheck (\pt -> withinPoints pt a)


mbCheck :: (a -> Bool) -> a -> Maybe a
mbCheck test a = if test a then Just a else Nothing

-- | Check the point is \"within\" the span of the line.
--
-- Note - this function is to be used \*after\* an intersection
-- has been found. Hence it is not exported.
--
withinPoints :: (Ord u, Fractional u, Tolerance u) 
             => Point2 u -> LineSegment u -> Bool
withinPoints (P2 x y) (LineSegment (P2 x0 y0) (P2 x1 y1)) =  
    between x (ordpair x0 x1) && between y (ordpair y0 y1)
  where
    ordpair a b     = (min a b, max a b)
    between a (s,t) = (s `tLTE` a) && (a `tLTE` t)



pve :: (Fractional u, Ord u, Tolerance u) => u -> Bool
pve a = a > eq_tolerance

nve :: (Fractional u, Ord u, Tolerance u) => u -> Bool
nve a = a < (negate eq_tolerance)



--------------------------------------------------------------------------------
-- intersection of line and Bezier curve

interCurveLine :: (Floating u , Ord u, Tolerance u)
               => BezierCurve u -> Line u -> Maybe (Point2 u)
interCurveLine c0 (Line p q) = step c0
  where
    eqline  = lineEquation p q
    step c  = case cut c eqline 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
 
-- | 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 make the return type
-- a bit involved.
--
cut :: (Floating u , Ord u, Tolerance u)
    => BezierCurve u -> LineEquation u -> Either (Point2 u) Bool
cut (BezierCurve p0 p1 p2 p3) eqline = 
    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
    d0  = pointLineDistance p0 eqline 
    d1  = pointLineDistance p1 eqline 
    d2  = pointLineDistance p2 eqline 
    d3  = pointLineDistance p3 eqline 



    


--------------------------------------------------------------------------------
-- Intersection of shape boundaries...



-- | 'findIntersect' :: @ radial_origin * theta * [line_segment] -> Maybe Point @
--
-- Find the first intersection of a line through @radial_origin@ 
-- at angle @theta@ and the supplied line segments, if there 
-- is one. 
--
findIntersect :: (Floating u, Real u, Ord u, Tolerance u)
              => Point2 u -> Radian -> [LineSegment u] 
              -> Maybe (Point2 u)
findIntersect radial_ogin ang = step 
  where
    line1       = inclinedLine radial_ogin ang
    step []     = Nothing
    step (x:xs) = case interLinesegLine x line1 of 
                     Just pt | quadrantCheck ang radial_ogin pt -> Just pt
                     _       -> step xs


-- | The tolerance on Radian equality should be acceptable...
--
quadrantCheck :: (Real u, Floating u) 
              => Radian -> Point2 u -> Point2 u -> Bool
quadrantCheck theta ctr pt = theta == lineAngle ctr pt