ESSA.hs

{-# LANGUAGE UnicodeSyntax, ForeignFunctionInterface #-}
{-# CFILES Numeric/Geometric/Predicates/ESSA/ESSAPrimitives.c #-}

{- | Hardware based, exact computation using the ESSA algorithm in double precision (1) 

     * We're using Float inputs on Double precision ESSA at the moment. Hopefully later we can add support for Double inputs on Quadruple ESSA

     * Line intersection is based on the algorithm presented in (4)

     * ccw and incircle based on (5)

     * See (2) and (3) for more information on the splitDouble operation

     * We assume realToFrac is broken and that CFloat == Float and CDouble == Double

     (1) Helmut Ratschek, Jon Rokne. \"Exact computation of the sign of a finite sum\". Applied Mathematics and Computation, Volume 99, Issue 2-3, Pages: 99-127, ISSN:0096-3003, 1999.   

     (2) Siegfried M. Rump. \"High precision evaluation of nonlinear functions\" 

     (3) T.J. Dekker. \"A Floating-Point Technique for Extending the Available Precision\". Numerische Mathematik, 18:224-242, 1971.

     (4) Marina Gavrilova, Jon Rokne. \"Reliable line segment intersection testing\"

     (5) Marina Gavrilova, Helmut Ratschek and Jon Rokne. \"Exact computation of Voronoi diagram and Delaunay triangulation\" 
-}

module Numeric.Geometric.Predicates.ESSA (cinttESSA, intersectESSA_SS2D, ccwESSA, incircleESSA, essa, splitDouble) where
import Numeric.Geometric.Primitives

import Data.Foldable (toList,Foldable)
import Control.Applicative
import Foreign.Ptr
import Foreign.Marshal.Alloc
import Foreign.Storable
import Foreign.Marshal.Array
import Foreign.C.Types
import System.IO.Unsafe


foreign import ccall unsafe ccw_essa ∷ Float → Float → Float → Float → Float → Float → IO CInt 
foreign import ccall unsafe incircle_essa ∷ Float → Float → Float → Float → Float → Float → Float → Float → IO CInt 
foreign import ccall unsafe intersect2D_essa ∷ Float → Float → Float → Float → Float → Float → Float → Float → Ptr CInt → Ptr CInt → IO CInt 
foreign import ccall unsafe cintt_essa ∷ Float → Float → Float → IO CInt 

foreign import ccall unsafe split_double ∷ Double → Ptr Double → Ptr Double → IO () 
foreign import ccall unsafe essa_double ∷ Ptr Double → CSize → IO CInt 

-- | Test if p3 is within the closed interval specified by [p1,p2]

cinttESSA ∷ Float → Float → Float → Bool
cinttESSA lo hi p = (unsafePerformIO $ cintt_essa (cast lo) (cast hi) (cast p)) /= 0

-- | Intersect two line segments

intersectESSA_SS2D ∷ LineSegment (Vector2 Float) → LineSegment (Vector2 Float) → LineIntersection 
intersectESSA_SS2D (a,b) (c,d) = unsafePerformIO $ 
         alloca (\ip1p → 
         alloca (\ip2p → do
                         x ← intersect2D_essa xi yi xj yj xk yk xl yl ip1p ip2p
                         case x of
                           0 → return NINP
                           1 → return Coincident
                           2 → return Parallel
                           3 → do
                               ip1 ← ip <$> peek ip1p
                               ip2 ← ip <$> peek ip2p
                               return (Intersecting (ip1,ip2))

                           _ → error "intersectESSA_SS2D: unexpected result from FFI"
                ))


    where 
      ip 0 = Endpoint0
      ip 1 = Endpoint1
      ip 2 = Between
      ip _ = error "intersectESSA_SS2D: unexpected intersection result from FFI"

      (xi,yi) = castVector a 
      (xj,yj) = castVector b
      (xk,yk) = castVector c
      (xl,yl) = castVector d      



{- | Counter-clockwise orientation test. Classifies p3 in relation to the line formed by p1 and p2. 

     Result: LT=Right, GT=Left, EQ=Coincident 
-}

ccwESSA ∷ Vector2 Float → Vector2 Float → Vector2 Float → Ordering
ccwESSA p1 p2 p3 = compare (unsafePerformIO $ ccw_essa x1 y1 x2 y2 x3 y3) 0
    where
      (x1,y1) = castVector p1  
      (x2,y2) = castVector p2
      (x3,y3) = castVector p3



{- | Test the relation of a point to the circle formed by (p1..p3). (p1..p3) must be in counterclockwise order. 

     Result: GT=inside, EQ=border, LT=outside.

     Note: this is the sum of 192 multiplications.
-}

incircleESSA ∷ (Vector2 Float, Vector2 Float, Vector2 Float) → Vector2 Float → Ordering
incircleESSA (a,b,c) d = compare (unsafePerformIO (incircle_essa xi yi xj yj xk yk xl yl)) 0

    where 
      (xi,yi) = castVector a 
      (xj,yj) = castVector b
      (xk,yk) = castVector c
      (xl,yl) = castVector d


{- | Compute the exact sign of the sum of the input sequence.   
     It is the caller's responsibility to ensure that the inputs have not suffered a loss of precision.
-}

essa ∷ (Functor t, Foldable t) => t Double → Ordering
essa = doubleESSA . fmap realToFrac -- maybe support 128bit quad later



{- | Split a 53 bit double into two 26 bit halves so that: @ let (lo,hi) = splitDouble x in x == lo + hi @

     The trick is that the sign is used as the additional bit.

     Note that the multiplication of two 26-bit floating point numbers is exact in double precision.

     If you're new to this function you may want to read paper (5), as using this function properly may be trickier than it seems.
-}

splitDouble ∷ Double → (Double,Double)
splitDouble a = unsafePerformIO (alloca (\xp → 
                                 alloca (\yp → do
                                               split_double a xp yp
                                               x ← peek xp
                                               y ← peek yp
                                               return (x,y))))




{--

foreign import ccall unsafe slope_essa ∷ Double → Double → Double → Double → IO CInt 

slopeESSA ∷ LineSegment (Vector2 Float) → Slope
slopeESSA (a,b) = case unsafePerformIO $ slope_essa x1 y1 x2 y2 of
                                0 → ZeroSlope
                                1 → UndefinedSlope
                                2 → NegativeSlope
                                3 → PositiveSlope

    where
      (x1,y1) = castVector a
      (x2,y2) = castVector b

--}

----------------------------------------------------

-- At the time realToFrac couldn't be trusted to do the right thing 
cast ∷ Float → Float
cast = id

castVector ∷ Vector2 Float → Vector2 Float
castVector = id

doubleESSA ∷ (Functor a, Foldable a) => a Double → Ordering
doubleESSA xs = compare v 0

    where
      v = unsafePerformIO $ withArrayLen (toList xs) f
      f = (\i p → essa_double p (fromIntegral i))