{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE TemplateHaskell #-}
module Algorithms.Geometry.SmallestEnclosingBall.RandomizedIncrementalConstruction where
import Algorithms.Geometry.SmallestEnclosingBall.Types
import Control.Lens
import Control.Monad.Random.Class
import Data.Ext
import qualified Data.Foldable as F
import Data.Geometry
import Data.Geometry.Ball
import qualified Data.List as L
import Data.List.NonEmpty
import Data.Maybe (fromMaybe)
import System.Random.Shuffle (shuffle)
smallestEnclosingDisk :: (Ord r, Fractional r, MonadRandom m)
=> [Point 2 r :+ p]
-> m (DiskResult p r)
smallestEnclosingDisk pts@(_:_:_) = ((\(p:q:pts') -> smallestEnclosingDisk' p q pts')
. F.toList) <$> shuffle pts
smallestEnclosingDisk _ = error "smallestEnclosingDisk: Too few points"
smallestEnclosingDisk' :: (Ord r, Fractional r)
=> Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p]
-> DiskResult p r
smallestEnclosingDisk' a b = foldr addPoint (initial a b) . L.tails
where
addPoint [] br = br
addPoint (p:pts) br@(DiskResult d _)
| (p^.core) `inClosedBall` d = br
| otherwise = smallestEnclosingDiskWithPoint p (a :| (b : pts))
smallestEnclosingDiskWithPoint :: (Ord r, Fractional r)
=> Point 2 r :+ p -> NonEmpty (Point 2 r :+ p)
-> DiskResult p r
smallestEnclosingDiskWithPoint p (a :| pts) = foldr addPoint (initial p a) $ L.tails pts
where
addPoint [] br = br
addPoint (q:pts') br@(DiskResult d _)
| (q^.core) `inClosedBall` d = br
| otherwise = smallestEnclosingDiskWithPoints p q (a:pts')
smallestEnclosingDiskWithPoints :: (Ord r, Fractional r)
=> Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p]
-> DiskResult p r
smallestEnclosingDiskWithPoints p q = foldr addPoint (initial p q)
where
addPoint r br@(DiskResult d _)
| (r^.core) `inClosedBall` d = br
| otherwise = DiskResult (circle' r) (Three p q r)
circle' r = fromMaybe degen $ disk (p^.core) (q^.core) (r^.core)
degen = error "smallestEnclosingDisk: Unhandled degeneracy, three points on a line"
initial :: Fractional r => Point 2 r :+ p -> Point 2 r :+ p -> DiskResult p r
initial p q = DiskResult (fromDiameter (p^.core) (q^.core)) (Two p q)