--------------------------------------------------------------------------------
-- |
-- Module      :  HGeometry.Ball.BoundaryPoints
-- Copyright   :  (C) Frank Staals
-- License     :  see the LICENSE file
-- Maintainer  :  Frank Staals
--
-- A ball in \(d\) dimensional space represented by \(d+1\) points on the boundary.
--
--------------------------------------------------------------------------------
{-# LANGUAGE UndecidableInstances #-}
module HGeometry.Ball.BoundaryPoints
  ( BallByPoints, BallByPoints'(..)
  , diskFromPoints
  ) where

import Control.Lens
import Data.Foldable1
import GHC.TypeLits
import HGeometry.Ball.Class
import HGeometry.Boundary
import HGeometry.Intersection
import HGeometry.Matrix
import HGeometry.Point
import HGeometry.Properties
import HGeometry.Vector
import Linear.Matrix (det33)
import Linear.V3 (V3(..))

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

-- | A ball specified by d+1 points which form a CCW simplex.
type BallByPoints point = BallByPoints' (1 + Dimension point) point

-- | A ball specified by k points which form a CCW simplex.
--
-- This is just a helper type to implement the BallByPoints type.
newtype BallByPoints' k point = BoundaryPoints (Vector k point)


deriving stock   instance Show (Vector k point) => Show (BallByPoints' k point)
deriving newtype instance Eq   (Vector k point) => Eq (BallByPoints' k point)
deriving newtype instance Ord  (Vector k point) => Ord (BallByPoints' k point)

deriving newtype instance Functor   (Vector k) => Functor   (BallByPoints' k)
deriving newtype instance Foldable  (Vector k) => Foldable  (BallByPoints' k)
deriving newtype instance Foldable1 (Vector k) => Foldable1 (BallByPoints' k)

instance Traversable (Vector k) => Traversable (BallByPoints' k) where
  traverse :: forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> BallByPoints' k a -> f (BallByPoints' k b)
traverse a -> f b
f (BoundaryPoints Vector k a
v) = Vector k b -> BallByPoints' k b
forall (k :: Nat) point. Vector k point -> BallByPoints' k point
BoundaryPoints (Vector k b -> BallByPoints' k b)
-> f (Vector k b) -> f (BallByPoints' k b)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (a -> f b) -> Vector k a -> f (Vector k b)
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> Vector k a -> f (Vector k b)
traverse a -> f b
f Vector k a
v
instance Traversable1 (Vector k) => Traversable1 (BallByPoints' k) where
  traverse1 :: forall (f :: * -> *) a b.
Apply f =>
(a -> f b) -> BallByPoints' k a -> f (BallByPoints' k b)
traverse1 a -> f b
f (BoundaryPoints Vector k a
v) = Vector k b -> BallByPoints' k b
forall (k :: Nat) point. Vector k point -> BallByPoints' k point
BoundaryPoints (Vector k b -> BallByPoints' k b)
-> f (Vector k b) -> f (BallByPoints' k b)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (a -> f b) -> Vector k a -> f (Vector k b)
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable1 t, Apply f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Apply f =>
(a -> f b) -> Vector k a -> f (Vector k b)
traverse1 a -> f b
f Vector k a
v


-- | smart constructor that makes sure the points are given in the right orientation may
-- fail if the points are colinear
diskFromPoints       :: (Point_ point 2 r, Ord r, Num r)
                     => point -> point -> point -> Maybe (BallByPoints point)
diskFromPoints :: forall point r.
(Point_ point 2 r, Ord r, Num r) =>
point -> point -> point -> Maybe (BallByPoints point)
diskFromPoints point
a point
b point
c = case point -> point -> point -> CCW
forall point r point' point''.
(Point_ point 2 r, Point_ point' 2 r, Point_ point'' 2 r, Num r,
 Ord r) =>
point -> point' -> point'' -> CCW
ccw point
a point
b point
c of
                         CCW
CCW      -> BallByPoints' 3 point -> Maybe (BallByPoints' 3 point)
BallByPoints' 3 point -> Maybe (BallByPoints point)
forall a. a -> Maybe a
Just (BallByPoints' 3 point -> Maybe (BallByPoints point))
-> (Vector 3 point -> BallByPoints' 3 point)
-> Vector 3 point
-> Maybe (BallByPoints point)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector 3 point -> BallByPoints' 3 point
forall (k :: Nat) point. Vector k point -> BallByPoints' k point
BoundaryPoints (Vector 3 point -> Maybe (BallByPoints point))
-> Vector 3 point -> Maybe (BallByPoints point)
forall a b. (a -> b) -> a -> b
$ point -> point -> point -> Vector 3 point
forall r. r -> r -> r -> Vector 3 r
Vector3 point
c point
b point
a
                         CCW
CW       -> BallByPoints' 3 point -> Maybe (BallByPoints' 3 point)
BallByPoints' 3 point -> Maybe (BallByPoints point)
forall a. a -> Maybe a
Just (BallByPoints' 3 point -> Maybe (BallByPoints point))
-> (Vector 3 point -> BallByPoints' 3 point)
-> Vector 3 point
-> Maybe (BallByPoints point)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector 3 point -> BallByPoints' 3 point
forall (k :: Nat) point. Vector k point -> BallByPoints' k point
BoundaryPoints (Vector 3 point -> Maybe (BallByPoints point))
-> Vector 3 point -> Maybe (BallByPoints point)
forall a b. (a -> b) -> a -> b
$ point -> point -> point -> Vector 3 point
forall r. r -> r -> r -> Vector 3 r
Vector3 point
a point
b point
c
                         CCW
CoLinear -> Maybe (BallByPoints' 3 point)
Maybe (BallByPoints point)
forall a. Maybe a
Nothing

type instance Dimension (BallByPoints' k point) = Dimension point
type instance NumType   (BallByPoints' k point) = NumType   point

instance ( Point_ point d r, Fractional r, 1 <= k, Has_ Vector_ k point, Has_ Metric_ d r
         , HasCenter (BallByPoints' k point) (Point d r)
         ) => Ball_ (BallByPoints' k point) (Point d r) where
  squaredRadius :: Getter (BallByPoints' k point) (NumType (BallByPoints' k point))
squaredRadius = (BallByPoints' k point -> NumType (BallByPoints' k point))
-> (NumType (BallByPoints' k point)
    -> f (NumType (BallByPoints' k point)))
-> BallByPoints' k point
-> f (BallByPoints' k point)
forall (p :: * -> * -> *) (f :: * -> *) s a.
(Profunctor p, Contravariant f) =>
(s -> a) -> Optic' p f s a
to ((BallByPoints' k point -> NumType (BallByPoints' k point))
 -> (NumType (BallByPoints' k point)
     -> f (NumType (BallByPoints' k point)))
 -> BallByPoints' k point
 -> f (BallByPoints' k point))
-> (BallByPoints' k point -> NumType (BallByPoints' k point))
-> (NumType (BallByPoints' k point)
    -> f (NumType (BallByPoints' k point)))
-> BallByPoints' k point
-> f (BallByPoints' k point)
forall a b. (a -> b) -> a -> b
$ \ball :: BallByPoints' k point
ball@(BoundaryPoints Vector k point
v) ->
                         let p :: Point d r
p = Vector k point
vVector k point
-> Getting (Point d r) (Vector k point) (Point d r) -> Point d r
forall s a. s -> Getting a s a -> a
^.(point -> Const (Point d r) point)
-> Vector k point -> Const (Point d r) (Vector k point)
forall vector (d :: Nat) r.
(Vector_ vector d r, 1 <= d) =>
IndexedLens' Int vector r
IndexedLens' Int (Vector k point) point
xComponent((point -> Const (Point d r) point)
 -> Vector k point -> Const (Point d r) (Vector k point))
-> ((Point d r -> Const (Point d r) (Point d r))
    -> point -> Const (Point d r) point)
-> Getting (Point d r) (Vector k point) (Point d r)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(Point d r -> Const (Point d r) (Point d r))
-> point -> Const (Point d r) point
forall point (d :: Nat) r.
Point_ point d r =>
Lens' point (Point d r)
Lens' point (Point d r)
asPoint
                         in Point d r -> Point d r -> r
forall r point (d :: Nat) point'.
(Num r, Point_ point d r, Point_ point' d r,
 Metric_ (Vector d r) d r) =>
point -> point' -> r
squaredEuclideanDist (BallByPoints' k point
ballBallByPoints' k point
-> Getting (Point d r) (BallByPoints' k point) (Point d r)
-> Point d r
forall s a. s -> Getting a s a -> a
^.Getting (Point d r) (BallByPoints' k point) (Point d r)
forall geom point. HasCenter geom point => Lens' geom point
Lens' (BallByPoints' k point) (Point d r)
center) Point d r
p

--------------------------------------------------------------------------------
-- * In ball

instance Point_ point 2 r => HasInBall (BallByPoints' 3 point) where
  inBall :: forall point (d :: Nat) r.
(Point_ point d r, Ord r, Num r,
 NumType (BallByPoints' 3 point) ~ r,
 Dimension (BallByPoints' 3 point) ~ d) =>
point -> BallByPoints' 3 point -> PointLocationResult
inBall point
q (BoundaryPoints (Vector3 point
a point
b point
c)) = let v :: point -> Vector 4 (NumType point)
v (Point2_ NumType point
x NumType point
y) = NumType point
-> NumType point
-> NumType point
-> NumType point
-> Vector 4 (NumType point)
forall r. r -> r -> r -> r -> Vector 4 r
Vector4 NumType point
x NumType point
y (NumType point
xNumType point -> NumType point -> NumType point
forall a. Num a => a -> a -> a
*NumType point
xNumType point -> NumType point -> NumType point
forall a. Num a => a -> a -> a
+NumType point
yNumType point -> NumType point -> NumType point
forall a. Num a => a -> a -> a
*NumType point
y) NumType point
1
                                                  m :: Matrix 4 4 r
m = Vector 4 (Vector 4 r) -> Matrix 4 4 r
forall (n :: Nat) (m :: Nat) r.
Vector n (Vector m r) -> Matrix n m r
Matrix (Vector 4 (Vector 4 r) -> Matrix 4 4 r)
-> Vector 4 (Vector 4 r) -> Matrix 4 4 r
forall a b. (a -> b) -> a -> b
$ Vector 4 r
-> Vector 4 r -> Vector 4 r -> Vector 4 r -> Vector 4 (Vector 4 r)
forall r. r -> r -> r -> r -> Vector 4 r
Vector4 (point -> Vector 4 (NumType point)
forall {point}.
(Dimension point ~ 2, Point_ point 2 (NumType point),
 Num (NumType point)) =>
point -> Vector 4 (NumType point)
v point
a)
                                                                       (point -> Vector 4 (NumType point)
forall {point}.
(Dimension point ~ 2, Point_ point 2 (NumType point),
 Num (NumType point)) =>
point -> Vector 4 (NumType point)
v point
b)
                                                                       (point -> Vector 4 (NumType point)
forall {point}.
(Dimension point ~ 2, Point_ point 2 (NumType point),
 Num (NumType point)) =>
point -> Vector 4 (NumType point)
v point
c)
                                                                       (point -> Vector 4 (NumType point)
forall {point}.
(Dimension point ~ 2, Point_ point 2 (NumType point),
 Num (NumType point)) =>
point -> Vector 4 (NumType point)
v point
q)
                                              in case Matrix 4 4 r -> r
forall (d :: Nat) r matrix.
(HasDeterminant d, Num r, Matrix_ matrix d d r) =>
matrix -> r
forall r matrix. (Num r, Matrix_ matrix 4 4 r) => matrix -> r
det Matrix 4 4 r
m r -> r -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` r
0 of
                                                   Ordering
LT -> PointLocationResult
Inside
                                                   Ordering
EQ -> PointLocationResult
OnBoundary
                                                   Ordering
GT -> PointLocationResult
Outside
    -- see:
    -- Lecture Notes on Geometric Robustness
    -- Jonathan Richard Shewchuk jrs@cs.berkeley.edu
    -- April 15, 2013


type instance Intersection (Point d r) (BallByPoints' k point) = Maybe (Point d r)

instance ( Point_ point d r, Ord r, Num r, HasInBall (BallByPoints' k point)
         ) => (Point d r) `HasIntersectionWith` (BallByPoints' k point) where
  intersects :: Point d r -> BallByPoints' k point -> Bool
intersects Point d r
q BallByPoints' k point
b = Point d r
q Point d r -> BallByPoints' k point -> PointLocationResult
forall ball point (d :: Nat) r.
(HasInBall ball, Point_ point d r, Ord r, Num r, NumType ball ~ r,
 Dimension ball ~ d) =>
point -> ball -> PointLocationResult
forall point (d :: Nat) r.
(Point_ point d r, Ord r, Num r,
 NumType (BallByPoints' k point) ~ r,
 Dimension (BallByPoints' k point) ~ d) =>
point -> BallByPoints' k point -> PointLocationResult
`inBall` BallByPoints' k point
b PointLocationResult -> PointLocationResult -> Bool
forall a. Eq a => a -> a -> Bool
/= PointLocationResult
Outside

instance ( Point_ point d r, Ord r, Num r, HasInBall (BallByPoints' k point)
         ) => (Point d r) `IsIntersectableWith` (BallByPoints' k point) where
  intersect :: Point d r
-> BallByPoints' k point
-> Intersection (Point d r) (BallByPoints' k point)
intersect Point d r
q BallByPoints' k point
b | Point d r
q Point d r -> BallByPoints' k point -> Bool
forall g h. HasIntersectionWith g h => g -> h -> Bool
`intersects` BallByPoints' k point
b = Point d r -> Maybe (Point d r)
forall a. a -> Maybe a
Just Point d r
q
                | Bool
otherwise        = Maybe (Point d r)
Intersection (Point d r) (BallByPoints' k point)
forall a. Maybe a
Nothing

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

instance (Point_ point 2 r, Fractional r) => HasCenter  (BallByPoints' 3 point) (Point 2 r) where
  center :: Lens' (BallByPoints' 3 point) (Point 2 r)
center = (BallByPoints' 3 point -> Point 2 r)
-> (BallByPoints' 3 point -> Point 2 r -> BallByPoints' 3 point)
-> Lens' (BallByPoints' 3 point) (Point 2 r)
forall s a b t. (s -> a) -> (s -> b -> t) -> Lens s t a b
lens BallByPoints' 3 point -> Point 2 r
BallByPoints point -> Point 2 r
forall point r.
(Point_ point 2 r, Fractional r) =>
BallByPoints point -> Point 2 r
computeCenter (\BallByPoints' 3 point
ball Point 2 r
c' -> let c :: Point 2 r
c = BallByPoints point -> Point 2 r
forall point r.
(Point_ point 2 r, Fractional r) =>
BallByPoints point -> Point 2 r
computeCenter BallByPoints' 3 point
BallByPoints point
ball
                                               v :: Vector 2 r
v = Point 2 r
c' Point 2 r -> Point 2 r -> Vector 2 r
forall point (d :: Nat) r.
(Affine_ point d r, Num r) =>
point -> point -> Vector d r
.-. Point 2 r
c
                                           in (point -> Vector 2 r -> point
forall point (d :: Nat) r.
(Affine_ point d r, Num r) =>
point -> Vector d r -> point
.+^ Vector 2 r
v) (point -> point) -> BallByPoints' 3 point -> BallByPoints' 3 point
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> BallByPoints' 3 point
ball
                              ) -- shifts the center to the new position

-- | computes the center of a disk
computeCenter :: (Point_ point 2 r, Fractional r) => BallByPoints point -> Point 2 r
computeCenter :: forall point r.
(Point_ point 2 r, Fractional r) =>
BallByPoints point -> Point 2 r
computeCenter (BoundaryPoints (Vector3 (Point2_ r
px r
py) (Point2_ r
qx r
qy) (Point2_ r
sx r
sy))) = Point 2 r
c
  where
    f :: a -> a -> a
f  a
x a
y = a
xa -> Integer -> a
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 a -> a -> a
forall a. Num a => a -> a -> a
+ a
ya -> Integer -> a
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2
    fx :: a -> a -> V3 a
fx a
x a
y = a -> a -> a -> V3 a
forall a. a -> a -> a -> V3 a
V3 (a -> a -> a
forall a. Num a => a -> a -> a
f a
x a
y) a
y       a
1
    fy :: a -> a -> V3 a
fy a
x a
y = a -> a -> a -> V3 a
forall a. a -> a -> a -> V3 a
V3 a
x       (a -> a -> a
forall a. Num a => a -> a -> a
f a
x a
y) a
1

    xnom :: r
xnom   = M33 r -> r
forall a. Num a => M33 a -> a
det33 (M33 r -> r) -> M33 r -> r
forall a b. (a -> b) -> a -> b
$ V3 r -> V3 r -> V3 r -> M33 r
forall a. a -> a -> a -> V3 a
V3 (r -> r -> V3 r
forall {a}. Num a => a -> a -> V3 a
fx r
px r
py) (r -> r -> V3 r
forall {a}. Num a => a -> a -> V3 a
fx r
qx r
qy) (r -> r -> V3 r
forall {a}. Num a => a -> a -> V3 a
fx r
sx r
sy)
    ynom :: r
ynom   = M33 r -> r
forall a. Num a => M33 a -> a
det33 (M33 r -> r) -> M33 r -> r
forall a b. (a -> b) -> a -> b
$ V3 r -> V3 r -> V3 r -> M33 r
forall a. a -> a -> a -> V3 a
V3 (r -> r -> V3 r
forall {a}. Num a => a -> a -> V3 a
fy r
px r
py) (r -> r -> V3 r
forall {a}. Num a => a -> a -> V3 a
fy r
qx r
qy) (r -> r -> V3 r
forall {a}. Num a => a -> a -> V3 a
fy r
sx r
sy)

    denom :: r
denom  = (r
2 r -> r -> r
forall a. Num a => a -> a -> a
*) (r -> r) -> (M33 r -> r) -> M33 r -> r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. M33 r -> r
forall a. Num a => M33 a -> a
det33 (M33 r -> r) -> M33 r -> r
forall a b. (a -> b) -> a -> b
$ V3 r -> V3 r -> V3 r -> M33 r
forall a. a -> a -> a -> V3 a
V3 (r -> r -> r -> V3 r
forall a. a -> a -> a -> V3 a
V3 r
px r
py r
1) (r -> r -> r -> V3 r
forall a. a -> a -> a -> V3 a
V3 r
qx r
qy r
1) (r -> r -> r -> V3 r
forall a. a -> a -> a -> V3 a
V3 r
sx r
sy r
1)
    c :: Point 2 r
c      = r -> r -> Point 2 r
forall r. r -> r -> Point 2 r
Point2 (r
xnom r -> r -> r
forall a. Fractional a => a -> a -> a
/ r
denom) (r
ynom r -> r -> r
forall a. Fractional a => a -> a -> a
/ r
denom)


-- instance Ball_ (BallByPoints point) point where
--   squaredRadius