{-# 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))