module Geometry.SmallestEnclosingCircle where

import Data.Function
import qualified Data.List as L
import qualified Data.Vector.Generic as G
import Geometry
import Geometry.Circle
import System.Random.Utils

smallestEnclosingCircle ::
  (Floating a, Ord a, G.Vector v (Point a)) => v (Point a) -> Circle a
smallestEnclosingCircle :: forall a (v :: * -> *).
(Floating a, Ord a, Vector v (Point a)) =>
v (Point a) -> Circle a
smallestEnclosingCircle =
  [Point a] -> [Point a] -> Circle a
forall {a}.
(Floating a, Ord a) =>
[Point a] -> [Point a] -> Circle a
welzl []
    ([Point a] -> Circle a)
-> (v (Point a) -> [Point a]) -> v (Point a) -> Circle a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. v (Point a) -> [Point a]
forall (v :: * -> *) a. Vector v a => v a -> [a]
G.toList
    (v (Point a) -> [Point a])
-> (v (Point a) -> v (Point a)) -> v (Point a) -> [Point a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (forall s. Mutable v s (Point a) -> ST s ())
-> v (Point a) -> v (Point a)
forall (v :: * -> *) a.
Vector v a =>
(forall s. Mutable v s a -> ST s ()) -> v a -> v a
G.modify
      ( \Mutable v s (Point a)
mps -> do
          StdGen
rng <- ST s StdGen
forall (m :: * -> *). PrimMonad m => m StdGen
newStdGenPrim
          StdGen -> Mutable v (PrimState (ST s)) (Point a) -> ST s ()
forall (m :: * -> *) (mv :: * -> * -> *) a g.
(PrimMonad m, MVector mv a, RandomGen g) =>
g -> mv (PrimState m) a -> m ()
shuffle StdGen
rng Mutable v s (Point a)
Mutable v (PrimState (ST s)) (Point a)
mps
      )
  where
    welzl :: [Point a] -> [Point a] -> Circle a
welzl [Point a
p0, Point a
p1, Point a
p2] [Point a]
_ =
      -- not equal to the smallest enclosing circle
      let c :: Point a
c = Point a -> Point a -> Point a -> Point a
forall a. Fractional a => Point a -> Point a -> Point a -> Point a
circumcenter Point a
p0 Point a
p1 Point a
p2
       in Point a -> a -> Circle a
forall a. Point a -> a -> Circle a
Circle Point a
c (Point a -> a
forall a. Floating a => Point a -> a
norm (Point a
p0 Point a -> Point a -> Point a
forall a. Num a => a -> a -> a
- Point a
c))
    welzl [] [] = Point a -> a -> Circle a
forall a. Point a -> a -> Circle a
Circle (a -> a -> Point a
forall a. a -> a -> Point a
P a
0.0 a
0.0) a
0.0
    welzl [Point a
p] [] = Point a -> a -> Circle a
forall a. Point a -> a -> Circle a
Circle Point a
p a
0.0
    welzl [Point a
p0, Point a
p1] [] = Point a -> a -> Circle a
forall a. Point a -> a -> Circle a
Circle ((a -> a) -> Point a -> Point a
forall a b. (a -> b) -> Point a -> Point b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (a
0.5 *) (Point a
p0 Point a -> Point a -> Point a
forall a. Num a => a -> a -> a
+ Point a
p1)) (Point a -> a
forall a. Floating a => Point a -> a
norm (Point a
p1 Point a -> Point a -> Point a
forall a. Num a => a -> a -> a
- Point a
p0) a -> a -> a
forall a. Num a => a -> a -> a
* a
0.5)
    welzl [Point a]
rs (Point a
p : [Point a]
ps)
      | Point a -> Circle a -> Bool
forall a. (Num a, Ord a) => Point a -> Circle a -> Bool
inCircle Point a
p Circle a
c = Circle a
c
      | Bool
otherwise = [Point a] -> [Point a] -> Circle a
welzl (Point a
p Point a -> [Point a] -> [Point a]
forall a. a -> [a] -> [a]
: [Point a]
rs) [Point a]
ps
      where
        !c :: Circle a
c = [Point a] -> [Point a] -> Circle a
welzl [Point a]
rs [Point a]
ps
    welzl [Point a]
_ [] = [Char] -> Circle a
forall a. HasCallStack => [Char] -> a
error [Char]
"unreachable"

naiveSmallestEnclosingCircle :: (Floating a, Ord a) => [Point a] -> Circle a
naiveSmallestEnclosingCircle :: forall a. (Floating a, Ord a) => [Point a] -> Circle a
naiveSmallestEnclosingCircle [] = Point a -> a -> Circle a
forall a. Point a -> a -> Circle a
Circle (a -> a -> Point a
forall a. a -> a -> Point a
P a
0.0 a
0.0) a
0.0
naiveSmallestEnclosingCircle [Point a
p] = Point a -> a -> Circle a
forall a. Point a -> a -> Circle a
Circle Point a
p a
0.0
naiveSmallestEnclosingCircle [Point a
p0, Point a
p1] =
  Point a -> a -> Circle a
forall a. Point a -> a -> Circle a
Circle ((a -> a) -> Point a -> Point a
forall a b. (a -> b) -> Point a -> Point b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (a
0.5 *) (Point a
p0 Point a -> Point a -> Point a
forall a. Num a => a -> a -> a
+ Point a
p1)) (Point a -> a
forall a. Floating a => Point a -> a
norm (Point a
p1 Point a -> Point a -> Point a
forall a. Num a => a -> a -> a
- Point a
p0) a -> a -> a
forall a. Num a => a -> a -> a
* a
0.5)
naiveSmallestEnclosingCircle [Point a]
ps =
  (a, Circle a) -> Circle a
forall a b. (a, b) -> b
snd
    ((a, Circle a) -> Circle a) -> (a, Circle a) -> Circle a
forall a b. (a -> b) -> a -> b
$ ((a, Circle a) -> (a, Circle a) -> Ordering)
-> [(a, Circle a)] -> (a, Circle a)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
L.minimumBy
      (a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (a -> a -> Ordering)
-> ((a, Circle a) -> a)
-> (a, Circle a)
-> (a, Circle a)
-> Ordering
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` (a, Circle a) -> a
forall a b. (a, b) -> a
fst)
      [(a
r, Point a -> a -> Circle a
forall a. Point a -> a -> Circle a
Circle Point a
c a
r) | Point a
c <- [Point a]
centers, let !r :: a
r = Point a -> a
radius Point a
c]
  where
    centers :: [Point a]
centers =
      [(a -> a) -> Point a -> Point a
forall a b. (a -> b) -> Point a -> Point b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (a
0.5 *) (Point a
p0 Point a -> Point a -> Point a
forall a. Num a => a -> a -> a
+ Point a
p1) | Point a
p0 <- [Point a]
ps, Point a
p1 <- [Point a]
ps]
        [Point a] -> [Point a] -> [Point a]
forall a. [a] -> [a] -> [a]
++ [ Point a -> Point a -> Point a -> Point a
forall a. Fractional a => Point a -> Point a -> Point a -> Point a
circumcenter Point a
p0 Point a
p1 Point a
p2
           | Point a
p0 <- [Point a]
ps
           , Point a
p1 <- [Point a]
ps
           , Point a
p2 <- [Point a]
ps
           , a -> a
forall a. Num a => a -> a
abs (Point a -> Point a -> Point a -> a
forall a. Num a => Point a -> Point a -> Point a -> a
area Point a
p0 Point a
p1 Point a
p2) a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0
           ]
    radius :: Point a -> a
radius Point a
c = [a] -> a
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Point a -> a
forall a. Floating a => Point a -> a
norm (Point a
p Point a -> Point a -> Point a
forall a. Num a => a -> a -> a
- Point a
c) | Point a
p <- [Point a]
ps]