hgeometry-0.14: src/Algorithms/Geometry/SmallestEnclosingBall/RIC.hs
--------------------------------------------------------------------------------
-- |
-- Module : Algorithms.Geometry.SmallestEnclosingBall.RIC
-- Copyright : (C) Frank Staals
-- License : see the LICENSE file
-- Maintainer : Frank Staals
--
-- An randomized algorithm to compute the smallest enclosing disk of a set of
-- \(n\) points in \(\mathbb{R}^2\). The expected running time is \(O(n)\).
--
--------------------------------------------------------------------------------
module Algorithms.Geometry.SmallestEnclosingBall.RIC(
smallestEnclosingDisk'
, smallestEnclosingDisk
, smallestEnclosingDiskWithPoint
, smallestEnclosingDiskWithPoints
) 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 List
import Data.List.NonEmpty(NonEmpty(..))
import Data.Maybe (fromMaybe, mapMaybe, catMaybes)
import Data.Ord (comparing)
import System.Random.Shuffle (shuffle)
-- import Data.RealNumber.Rational
-- import Debug.Trace
--------------------------------------------------------------------------------
-- | Compute the smallest enclosing disk of a set of points,
-- implemented using randomized incremental construction.
--
-- pre: the input has at least two points.
--
-- running time: expected \(O(n)\) time, where \(n\) is the number of input points.
smallestEnclosingDisk :: (Ord r, Fractional r, MonadRandom m
-- , Show r, Show p
)
=> [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"
-- | Smallest enclosing disk.
smallestEnclosingDisk' :: (Ord r, Fractional r
-- , Show r, Show p
)
=> Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p]
-> DiskResult p r
smallestEnclosingDisk' a b = foldr addPoint (initial a b) . List.tails
where
-- The empty case occurs only initially
addPoint [] br = br
addPoint (p:pts) br@(DiskResult d _)
| (p^.core) `inClosedBall` d = br
| otherwise = fromJust' $ smallestEnclosingDiskWithPoint p (a :| (b : pts))
fromJust' = fromMaybe (error "smallestEncosingDisk' : fromJust, absurd")
-- | Smallest enclosing disk, given that p should be on it.
smallestEnclosingDiskWithPoint :: (Ord r, Fractional r
-- , Show r, Show p
)
=> Point 2 r :+ p -> NonEmpty (Point 2 r :+ p)
-> Maybe (DiskResult p r)
smallestEnclosingDiskWithPoint p (a :| pts) = foldr addPoint (Just $ initial p a) $ List.tails pts
where
addPoint [] br = br
addPoint (q:pts') br@(Just (DiskResult d _))
| (q^.core) `inClosedBall` d = br
| otherwise = smallestEnclosingDiskWithPoints p q (a:pts')
addPoint _ br = br
-- | Smallest enclosing disk, given that p and q should be on it
--
-- running time: \(O(n)\)
smallestEnclosingDiskWithPoints :: (Ord r, Fractional r
-- , Show r, Show p
)
=> Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p]
-> Maybe (DiskResult p r)
smallestEnclosingDiskWithPoints p q ps = minimumOn (^.enclosingDisk.squaredRadius)
$ catMaybes [mkEnclosingDisk dl, mkEnclosingDisk dr, mdc]
where
centers = mapMaybe disk' ps
-- generate a disk with p q and r
disk' r = (r:+) <$> disk (p^.core) (q^.core) (r^.core)
-- partition the points in to those on the left and those on the
-- right. Note that centers still contains only those points (and
-- disks) for which the three points are not colinear. So the
-- points are either on the left or on the right.
(leftCenters,rightCenters) = List.partition (\(r :+ _) -> ccw' p q r == CCW) centers
-- note that we consider 'leftmost' with respect to going from p
-- to q. This does not really have a global meaning.
-- we need to find the leftmost and rightmost center on the
-- bisector. In case there are left-centers, this means that among
-- the left centers we want to find the point that is furthest way
-- from p (or q). If there are no left-centers, we with to find
-- the closest one among the right-centers.
leftDist z = let c = z^.extra.center
s = if ccw' p q c == CCW then 1 else -1
in s * squaredEuclideanDist (p^.core) (c^.core)
dl = maximumOn leftDist leftCenters -- disk that has the "leftmost" center
dr = minimumOn leftDist rightCenters -- disk that has the "rightmost" center
-- diameteral disk
dd = fromDiameter (p^.core) (q^.core)
mdc | isEnclosingDisk dd ps = Just $ DiskResult dd (Two p q)
| otherwise = Nothing
-- test if d is an enclosing disk.
mkEnclosingDisk md = md >>= mkEnclosingDisk'
mkEnclosingDisk' (r :+ d) | isEnclosingDisk d ps = Just (DiskResult d (Three p q r))
| otherwise = Nothing
isEnclosingDisk :: (Foldable t, Ord r, Num r)
=> Disk p r -> t (Point 2 r :+ extra) -> Bool
isEnclosingDisk d = all (\s -> (s^.core) `inClosedBall` d)
-- | Constructs the initial 'DiskResult' from two points
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)
maximumOn :: Ord b => (a -> b) -> [a] -> Maybe a
maximumOn f = \case
[] -> Nothing
xs -> Just $ List.maximumBy (comparing f) xs
minimumOn :: Ord b => (a -> b) -> [a] -> Maybe a
minimumOn f = \case
[] -> Nothing
xs -> Just $ List.minimumBy (comparing f) xs
--------------------------------------------------------------------------------
-- test :: Maybe (DiskResult () Rational)
-- test = smallestEnclosingDiskWithPoints p q myPts
-- where
-- p = ext $ Point2 0 (-6)
-- q = ext $ Point2 0 6
-- myPts = map ext [Point2 5 1, Point2 3 3, Point2 (-2) 2, Point2 (-4) 5]
-- disk'' r = (r:+) <$> disk (p^.core) (q^.core) (r^.core)
-- where
-- p = ext $ Point2 0 (-6)
-- q = ext $ Point2 0 6
-- maartenBug :: DiskResult () Double
-- maartenBug = let (p:q:rest) = maartenBug'
-- in smallestEnclosingDisk' p q rest
-- maartenBug' :: [Point 2 Double :+ ()]
-- maartenBug' = [ Point2 (7.2784424e-3) (249.23) :+ ()
-- , Point2 (-5.188493 ) (249.23) :+ ()
-- , Point2 (-10.382694 ) (249.23) :+ ()
-- , Point2 (-15.575621 ) (249.23) :+ ()
-- , Point2 (0.0 ) (249.23) :+ ()
-- , Point2 (0.0 ) (239.9031) :+ ()
-- , Point2 (0.0 ) (230.37791) :+ ()
-- , Point2 (0.0 ) (220.67882) :+ ()
-- ]