--------------------------------------------------------------------------------
-- |
-- Module      :  Algorithms.Geometry.SmallestEnclosingBall.Naive
-- Copyright   :  (C) Frank Staals
-- License     :  see the LICENSE file
-- Maintainer  :  Frank Staals
--
-- Naive implementation to compute the smallest enclosing disk of a set of
-- points in \(\mathbb{R}^2\)
--
--------------------------------------------------------------------------------
module Algorithms.Geometry.SmallestEnclosingBall.Naive
  ( smallestEnclosingDisk
  , enclosesAll
  ) where

-- just for the types
import Control.Lens
import Data.Ext
import Algorithms.Geometry.SmallestEnclosingBall.Types
import Data.Geometry.Ball
import Data.Geometry.Point
import Data.List (minimumBy)
import Data.Function (on)
import Data.Maybe (fromMaybe)
import Data.Util(uniquePairs, uniqueTriplets)
import qualified Data.Util as Util
--------------------------------------------------------------------------------

-- | Horrible \( O(n^4) \) implementation that simply tries all disks, checks if they
-- enclose all points, and takes the largest one. Basically, this is only useful
-- to check correctness of the other algorithm(s)
smallestEnclosingDisk          :: (Ord r, Fractional r)
                               => [Point 2 r :+ p]
                               -> DiskResult p r
smallestEnclosingDisk :: [Point 2 r :+ p] -> DiskResult p r
smallestEnclosingDisk pts :: [Point 2 r :+ p]
pts@(Point 2 r :+ p
_:Point 2 r :+ p
_:[Point 2 r :+ p]
_) = [Point 2 r :+ p] -> [DiskResult p r] -> DiskResult p r
forall r p.
(Ord r, Num r) =>
[Point 2 r :+ p] -> [DiskResult p r] -> DiskResult p r
smallestEnclosingDisk' [Point 2 r :+ p]
pts ([DiskResult p r] -> DiskResult p r)
-> [DiskResult p r] -> DiskResult p r
forall a b. (a -> b) -> a -> b
$
                                      [Point 2 r :+ p] -> [DiskResult p r]
forall r p. Fractional r => [Point 2 r :+ p] -> [DiskResult p r]
pairs [Point 2 r :+ p]
pts [DiskResult p r] -> [DiskResult p r] -> [DiskResult p r]
forall a. [a] -> [a] -> [a]
++ [Point 2 r :+ p] -> [DiskResult p r]
forall r p.
(Ord r, Fractional r) =>
[Point 2 r :+ p] -> [DiskResult p r]
triplets [Point 2 r :+ p]
pts
smallestEnclosingDisk [Point 2 r :+ p]
_           = [Char] -> DiskResult p r
forall a. HasCallStack => [Char] -> a
error [Char]
"smallestEnclosingDisk: Too few points"

pairs     :: Fractional r => [Point 2 r :+ p] -> [DiskResult p r]
pairs :: [Point 2 r :+ p] -> [DiskResult p r]
pairs [Point 2 r :+ p]
pts = [ Disk () r -> TwoOrThree (Point 2 r :+ p) -> DiskResult p r
forall p r.
Disk () r -> TwoOrThree (Point 2 r :+ p) -> DiskResult p r
DiskResult (Point 2 r -> Point 2 r -> Disk () r
forall (d :: Nat) r.
(Arity d, Fractional r) =>
Point d r -> Point d r -> Ball d () r
fromDiameter (Point 2 r :+ p
a(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core) (Point 2 r :+ p
b(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core)) ((Point 2 r :+ p) -> (Point 2 r :+ p) -> TwoOrThree (Point 2 r :+ p)
forall a. a -> a -> TwoOrThree a
Two Point 2 r :+ p
a Point 2 r :+ p
b)
            | Util.Two Point 2 r :+ p
a Point 2 r :+ p
b <- [Point 2 r :+ p] -> [Two (Point 2 r :+ p)]
forall a. [a] -> [Two a]
uniquePairs [Point 2 r :+ p]
pts]

triplets     :: (Ord r, Fractional r) => [Point 2 r :+ p] -> [DiskResult p r]
triplets :: [Point 2 r :+ p] -> [DiskResult p r]
triplets [Point 2 r :+ p]
pts = [Disk () r -> TwoOrThree (Point 2 r :+ p) -> DiskResult p r
forall p r.
Disk () r -> TwoOrThree (Point 2 r :+ p) -> DiskResult p r
DiskResult ((Point 2 r :+ p)
-> (Point 2 r :+ p) -> (Point 2 r :+ p) -> Disk () r
forall r p.
(Ord r, Fractional r) =>
(Point 2 r :+ p)
-> (Point 2 r :+ p) -> (Point 2 r :+ p) -> Disk () r
disk' Point 2 r :+ p
a Point 2 r :+ p
b Point 2 r :+ p
c) ((Point 2 r :+ p)
-> (Point 2 r :+ p)
-> (Point 2 r :+ p)
-> TwoOrThree (Point 2 r :+ p)
forall a. a -> a -> a -> TwoOrThree a
Three Point 2 r :+ p
a Point 2 r :+ p
b Point 2 r :+ p
c)
               | Util.Three Point 2 r :+ p
a Point 2 r :+ p
b Point 2 r :+ p
c <- [Point 2 r :+ p] -> [Three (Point 2 r :+ p)]
forall a. [a] -> [Three a]
uniqueTriplets [Point 2 r :+ p]
pts]

{- HLINT ignore disk' -}
disk'       :: (Ord r, Fractional r)
            => Point 2 r :+ p -> Point 2 r :+ p -> Point 2 r :+ p -> Disk () r
disk' :: (Point 2 r :+ p)
-> (Point 2 r :+ p) -> (Point 2 r :+ p) -> Disk () r
disk' Point 2 r :+ p
a Point 2 r :+ p
b Point 2 r :+ p
c = Disk () r -> Maybe (Disk () r) -> Disk () r
forall a. a -> Maybe a -> a
fromMaybe Disk () r
degen (Maybe (Disk () r) -> Disk () r) -> Maybe (Disk () r) -> Disk () r
forall a b. (a -> b) -> a -> b
$ Point 2 r -> Point 2 r -> Point 2 r -> Maybe (Disk () r)
forall r.
(Eq r, Fractional r) =>
Point 2 r -> Point 2 r -> Point 2 r -> Maybe (Disk () r)
disk (Point 2 r :+ p
a(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core) (Point 2 r :+ p
b(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core) (Point 2 r :+ p
c(Point 2 r :+ p)
-> Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r) -> Point 2 r
forall s a. s -> Getting a s a -> a
^.Getting (Point 2 r) (Point 2 r :+ p) (Point 2 r)
forall core extra core'.
Lens (core :+ extra) (core' :+ extra) core core'
core)
  where
    -- if the points are colinear, select the disk by the diametral pair
    degen :: Disk () r
degen = ([Point 2 r :+ p] -> [DiskResult p r] -> DiskResult p r
forall r p.
(Ord r, Num r) =>
[Point 2 r :+ p] -> [DiskResult p r] -> DiskResult p r
smallestEnclosingDisk' [Point 2 r :+ p
a,Point 2 r :+ p
b,Point 2 r :+ p
c] ([DiskResult p r] -> DiskResult p r)
-> [DiskResult p r] -> DiskResult p r
forall a b. (a -> b) -> a -> b
$ [Point 2 r :+ p] -> [DiskResult p r]
forall r p. Fractional r => [Point 2 r :+ p] -> [DiskResult p r]
pairs [Point 2 r :+ p
a,Point 2 r :+ p
b,Point 2 r :+ p
c])DiskResult p r
-> Getting (Disk () r) (DiskResult p r) (Disk () r) -> Disk () r
forall s a. s -> Getting a s a -> a
^.Getting (Disk () r) (DiskResult p r) (Disk () r)
forall p r. Lens' (DiskResult p r) (Disk () r)
enclosingDisk


-- | Given a list of canidate enclosing disks, report the smallest one.
smallestEnclosingDisk'     :: (Ord r, Num r)
                           => [Point 2 r :+ p] -> [DiskResult p r] -> DiskResult p r
smallestEnclosingDisk' :: [Point 2 r :+ p] -> [DiskResult p r] -> DiskResult p r
smallestEnclosingDisk' [Point 2 r :+ p]
pts = (DiskResult p r -> DiskResult p r -> Ordering)
-> [DiskResult p r] -> DiskResult p r
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy (r -> r -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (r -> r -> Ordering)
-> (DiskResult p r -> r)
-> DiskResult p r
-> DiskResult p r
-> Ordering
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` (DiskResult p r -> Getting r (DiskResult p r) r -> r
forall s a. s -> Getting a s a -> a
^.(Disk () r -> Const r (Disk () r))
-> DiskResult p r -> Const r (DiskResult p r)
forall p r. Lens' (DiskResult p r) (Disk () r)
enclosingDisk((Disk () r -> Const r (Disk () r))
 -> DiskResult p r -> Const r (DiskResult p r))
-> ((r -> Const r r) -> Disk () r -> Const r (Disk () r))
-> Getting r (DiskResult p r) r
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(r -> Const r r) -> Disk () r -> Const r (Disk () r)
forall (d :: Nat) p r. Lens' (Ball d p r) r
squaredRadius))
                           ([DiskResult p r] -> DiskResult p r)
-> ([DiskResult p r] -> [DiskResult p r])
-> [DiskResult p r]
-> DiskResult p r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (DiskResult p r -> Bool) -> [DiskResult p r] -> [DiskResult p r]
forall a. (a -> Bool) -> [a] -> [a]
filter (DiskResult p r -> [Point 2 r :+ p] -> Bool
forall r p q.
(Num r, Ord r) =>
DiskResult p r -> [Point 2 r :+ q] -> Bool
`enclosesAll` [Point 2 r :+ p]
pts)

-- | check if a disk encloses all points
enclosesAll   :: (Num r, Ord r) => DiskResult p r -> [Point 2 r :+ q] -> Bool
enclosesAll :: DiskResult p r -> [Point 2 r :+ q] -> Bool
enclosesAll DiskResult p r
d = ((Point 2 r :+ q) -> Bool) -> [Point 2 r :+ q] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\(Point 2 r
p :+ q
_) -> Point 2 r
p Point 2 r -> Ball 2 () r -> Bool
forall (d :: Nat) r p.
(Arity d, Ord r, Num r) =>
Point d r -> Ball d p r -> Bool
`inClosedBall` (DiskResult p r
dDiskResult p r
-> Getting (Ball 2 () r) (DiskResult p r) (Ball 2 () r)
-> Ball 2 () r
forall s a. s -> Getting a s a -> a
^.Getting (Ball 2 () r) (DiskResult p r) (Ball 2 () r)
forall p r. Lens' (DiskResult p r) (Disk () r)
enclosingDisk))