-- | Fast, high-quality Poisson disc sampling, based on [Bridson’s algorithm](https://dl.acm.org/doi/10.1145/1278780.1278807),
-- then [improved by Martin Robers](https://observablehq.com/@techsparx/an-improvement-on-bridsons-algorithm-for-poisson-disc-samp/2),
-- then improved an \(\varepsilon\) bit more here (4 cells in 'neighbouringPoints' can be skipped).
module Geometry.Algorithms.Sampling.PoissonDisc (
      poissonDisc
    , poissonDiscForest
) where



import           Control.Monad.Primitive
import           Data.List
import           Data.Map                (Map)
import qualified Data.Map                as M
import           Data.Maybe
import           Data.Set                (Set)
import qualified Data.Set                as S
import           Data.Vector             (Vector)
import qualified Data.Vector.Extended    as V
import           System.Random.MWC

import Geometry



-- $setup
-- >>> import           Draw
-- >>> import           Control.Monad.ST
-- >>> import           Numerics.Interpolation
-- >>> import           Graphics.Rendering.Cairo as C
-- >>> import qualified System.Random.MWC        as MWC

-- | Newtype safety wrapper
newtype CellSize = CellSize Double
    deriving (CellSize -> CellSize -> Bool
(CellSize -> CellSize -> Bool)
-> (CellSize -> CellSize -> Bool) -> Eq CellSize
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: CellSize -> CellSize -> Bool
== :: CellSize -> CellSize -> Bool
$c/= :: CellSize -> CellSize -> Bool
/= :: CellSize -> CellSize -> Bool
Eq, Eq CellSize
Eq CellSize
-> (CellSize -> CellSize -> Ordering)
-> (CellSize -> CellSize -> Bool)
-> (CellSize -> CellSize -> Bool)
-> (CellSize -> CellSize -> Bool)
-> (CellSize -> CellSize -> Bool)
-> (CellSize -> CellSize -> CellSize)
-> (CellSize -> CellSize -> CellSize)
-> Ord CellSize
CellSize -> CellSize -> Bool
CellSize -> CellSize -> Ordering
CellSize -> CellSize -> CellSize
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: CellSize -> CellSize -> Ordering
compare :: CellSize -> CellSize -> Ordering
$c< :: CellSize -> CellSize -> Bool
< :: CellSize -> CellSize -> Bool
$c<= :: CellSize -> CellSize -> Bool
<= :: CellSize -> CellSize -> Bool
$c> :: CellSize -> CellSize -> Bool
> :: CellSize -> CellSize -> Bool
$c>= :: CellSize -> CellSize -> Bool
>= :: CellSize -> CellSize -> Bool
$cmax :: CellSize -> CellSize -> CellSize
max :: CellSize -> CellSize -> CellSize
$cmin :: CellSize -> CellSize -> CellSize
min :: CellSize -> CellSize -> CellSize
Ord, Int -> CellSize -> ShowS
[CellSize] -> ShowS
CellSize -> String
(Int -> CellSize -> ShowS)
-> (CellSize -> String) -> ([CellSize] -> ShowS) -> Show CellSize
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> CellSize -> ShowS
showsPrec :: Int -> CellSize -> ShowS
$cshow :: CellSize -> String
show :: CellSize -> String
$cshowList :: [CellSize] -> ShowS
showList :: [CellSize] -> ShowS
Show)

-- | Sample points using the Poisson Disc algorithm, which yields a visually
-- uniform distribution. This is opposed to uniformly distributed points yield
-- clumps and empty areas, which is often undesirable for generative art.
--
-- <<docs/haddock/Geometry/Algorithms/Sampling/PoissonDisc/poisson_disc.svg>>
--
-- === __(image code)__
-- >>> :{
-- haddockRender "Geometry/Algorithms/Sampling/PoissonDisc/poisson_disc.svg" 300 300 $ \_ -> do
--     let points = runST $ do
--             gen <- MWC.create
--             poissonDisc gen (shrinkBoundingBox 30 [zero, Vec2 300 300]) 10 4
--     for_ (zip [0..] points) $ \(i,p) -> do
--         setColor (mma i)
--         sketch (Circle p 2)
--         fill
-- :}
-- Generated file: size 115KB, crc32: 0x43236d3e
poissonDisc
    :: (PrimMonad m, HasBoundingBox boundingBox)
    => Gen (PrimState m) -- ^ RNG from mwc-random. 'create' yields the default (static) RNG.
    -> boundingBox -- ^ Region to generate points in
    -> Double      -- ^ Radius around each point no other points are genereted. Smaller values yield more points.
    -> Int         -- ^ \(k\) parameter: per point, how many attempts should be made to find an empty spot?
                   --   Typical values: 3-10. Higher values are slower, but increase result quality.
    -> m [Vec2]
poissonDisc :: forall (m :: * -> *) boundingBox.
(PrimMonad m, HasBoundingBox boundingBox) =>
Gen (PrimState m) -> boundingBox -> Double -> Int -> m [Vec2]
poissonDisc Gen (PrimState m)
gen boundingBox
bb' Double
radius Int
k = do
    let bb :: BoundingBox
bb@(BoundingBox Vec2
minV Vec2
maxV) = boundingBox -> BoundingBox
forall a. HasBoundingBox a => a -> BoundingBox
boundingBox boundingBox
bb'
        cellSize :: CellSize
cellSize = Double -> CellSize
CellSize (Double
radiusDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double -> Double
forall a. Floating a => a -> a
sqrt Double
2)
    Vec2
initialPoint <- (Vec2, Vec2) -> Gen (PrimState m) -> m Vec2
forall a g (m :: * -> *).
(UniformRange a, StatefulGen g m) =>
(a, a) -> g -> m a
forall g (m :: * -> *).
StatefulGen g m =>
(Vec2, Vec2) -> g -> m Vec2
uniformRM (Vec2
minV, Vec2
maxV) Gen (PrimState m)
gen
    let initialGrid :: Map (Int, Int) Vec2
initialGrid = (Int, Int) -> Vec2 -> Map (Int, Int) Vec2
forall k a. k -> a -> Map k a
M.singleton (CellSize -> Vec2 -> (Int, Int)
gridCell CellSize
cellSize Vec2
initialPoint) Vec2
initialPoint
        initialActive :: Set Vec2
initialActive = Vec2 -> Set Vec2
forall a. a -> Set a
S.singleton Vec2
initialPoint

        sampleLoop :: Map (Int, Int) Vec2 -> Set Vec2 -> f [Vec2]
sampleLoop Map (Int, Int) Vec2
_grid Set Vec2
active
            | Set Vec2 -> Bool
forall a. Set a -> Bool
S.null Set Vec2
active = [Vec2] -> f [Vec2]
forall a. a -> f a
forall (f :: * -> *) a. Applicative f => a -> f a
pure []
        sampleLoop Map (Int, Int) Vec2
grid Set Vec2
active = do
            Vec2
activeSample <- Gen (PrimState f) -> Set Vec2 -> f Vec2
forall (m :: * -> *) a.
PrimMonad m =>
Gen (PrimState m) -> Set a -> m a
randomSetElement Gen (PrimState m)
Gen (PrimState f)
gen Set Vec2
active
            Vector Vec2
candidates <- Gen (PrimState f)
-> Int -> BoundingBox -> Double -> Vec2 -> f (Vector Vec2)
forall (m :: * -> *).
PrimMonad m =>
Gen (PrimState m)
-> Int -> BoundingBox -> Double -> Vec2 -> m (Vector Vec2)
candidatesAroundSample Gen (PrimState m)
Gen (PrimState f)
gen Int
k BoundingBox
bb Double
radius Vec2
activeSample

            let validPoint :: Vec2 -> Bool
validPoint Vec2
candidate =
                    let neighbours :: [Vec2]
neighbours = CellSize -> Map (Int, Int) Vec2 -> Vec2 -> [Vec2]
neighbouringPoints CellSize
cellSize Map (Int, Int) Vec2
grid Vec2
candidate
                        tooClose :: Vec2 -> Bool
tooClose Vec2
neighbour = Vec2 -> Double
normSquare (Vec2
candidate Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
neighbour) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
radiusDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2
                    in Bool -> Bool
not ((Vec2 -> Bool) -> [Vec2] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any Vec2 -> Bool
tooClose [Vec2]
neighbours)

            case (Vec2 -> Bool) -> Vector Vec2 -> Maybe Vec2
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
find Vec2 -> Bool
validPoint Vector Vec2
candidates of
                Maybe Vec2
Nothing -> do
                    [Vec2]
rest <- Map (Int, Int) Vec2 -> Set Vec2 -> f [Vec2]
sampleLoop Map (Int, Int) Vec2
grid (Vec2 -> Set Vec2 -> Set Vec2
forall a. Ord a => a -> Set a -> Set a
S.delete Vec2
activeSample Set Vec2
active)
                    [Vec2] -> f [Vec2]
forall a. a -> f a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Vec2
activeSample Vec2 -> [Vec2] -> [Vec2]
forall a. a -> [a] -> [a]
: [Vec2]
rest)
                Just Vec2
sample -> Map (Int, Int) Vec2 -> Set Vec2 -> f [Vec2]
sampleLoop
                    ((Int, Int) -> Vec2 -> Map (Int, Int) Vec2 -> Map (Int, Int) Vec2
forall k a. Ord k => k -> a -> Map k a -> Map k a
M.insert (CellSize -> Vec2 -> (Int, Int)
gridCell CellSize
cellSize Vec2
sample) Vec2
sample Map (Int, Int) Vec2
grid)
                    (Vec2 -> Set Vec2 -> Set Vec2
forall a. Ord a => a -> Set a -> Set a
S.insert Vec2
sample Set Vec2
active)

    Map (Int, Int) Vec2 -> Set Vec2 -> m [Vec2]
forall {f :: * -> *}.
(PrimState f ~ PrimState m, PrimMonad f) =>
Map (Int, Int) Vec2 -> Set Vec2 -> f [Vec2]
sampleLoop Map (Int, Int) Vec2
initialGrid Set Vec2
initialActive

-- | 'poissonDisc', but keeps track of which parent spawned which children. While
-- this algorithm does a bit more processing to reassemble the trees, it allows
-- specifying the starting points explicitly. ('poissonDisc' would also allow this,
-- but since it’s mostly used for simply sampling, this additional config option
-- would just clutter the API.)
--
-- <<docs/haddock/Geometry/Algorithms/Sampling/PoissonDisc/poisson_disc_forest.svg>>
--
-- === __(image code)__
-- >>> :{
-- haddockRender "Geometry/Algorithms/Sampling/PoissonDisc/poisson_disc_forest.svg" 300 300 $ \_ -> do
--     let initialPoints = [Vec2 50 50, Vec2 150 150, Vec2 250 250]
--         forest = runST $ do
--             gen <- MWC.create
--             poissonDiscForest gen (shrinkBoundingBox 30 [zero, Vec2 300 300]) 7 10 initialPoints
--     let paint color i parent = do
--             let parentRadius = 2.5
--             cairoScope $ do
--                 setColor (color (sin (lerpID (0,20) (0,pi) i)^2))
--                 sketch (Circle parent parentRadius)
--                 fill
--             case M.lookup parent forest of
--                 Nothing -> pure () -- Should never happen: only the roots have no parents
--                 Just children -> for_ children $ \child -> do
--                     cairoScope $ do
--                         setColor black
--                         sketch (resizeLineSymmetric (\l -> l-2*parentRadius) (Line parent child))
--                         stroke
--                     paint color (i+1) child
--     for_ (zip [flare, crest, flare] initialPoints) $ \(color, p0) -> paint color 0 p0
-- :}
-- Generated file: size 474KB, crc32: 0x28fe4a90
poissonDiscForest
    :: (PrimMonad m, HasBoundingBox boundingBox)
    => Gen (PrimState m)       -- ^ RNG from mwc-random. 'create' yields the default (static) RNG.
    -> boundingBox             -- ^ Region to generate points in
    -> Double                  -- ^ Radius around each point no other points are genereted. Smaller values yield more points.
    -> Int                     -- ^ \(k\) parameter: per point, how many attempts should be made to find an empty spot?
                               --   Typical values: 3-10. Higher values are slower, but increase result quality.
    -> [Vec2]                  -- ^ Initial points to start sampling from.
    -> m (Map Vec2 (Set Vec2)) -- ^ Map of parent to children.
poissonDiscForest :: forall (m :: * -> *) boundingBox.
(PrimMonad m, HasBoundingBox boundingBox) =>
Gen (PrimState m)
-> boundingBox
-> Double
-> Int
-> [Vec2]
-> m (Map Vec2 (Set Vec2))
poissonDiscForest Gen (PrimState m)
gen boundingBox
bb' Double
radius Int
k [Vec2]
initialPoints = do
    let bb :: BoundingBox
bb = boundingBox -> BoundingBox
forall a. HasBoundingBox a => a -> BoundingBox
boundingBox boundingBox
bb'
        cellSize :: CellSize
cellSize = Double -> CellSize
CellSize (Double
radiusDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double -> Double
forall a. Floating a => a -> a
sqrt Double
2)
    let initialGrid :: Map (Int, Int) Vec2
initialGrid = [((Int, Int), Vec2)] -> Map (Int, Int) Vec2
forall k a. Ord k => [(k, a)] -> Map k a
M.fromList [(CellSize -> Vec2 -> (Int, Int)
gridCell CellSize
cellSize Vec2
initialPoint, Vec2
initialPoint) | Vec2
initialPoint <- [Vec2]
initialPoints]
        initialActive :: Set ChildParent
initialActive = [ChildParent] -> Set ChildParent
forall a. Ord a => [a] -> Set a
S.fromList [Vec2 -> Maybe Vec2 -> ChildParent
ChildParent Vec2
initialPoint Maybe Vec2
forall a. Maybe a
Nothing | Vec2
initialPoint <- [Vec2]
initialPoints]

        sampleLoop :: Map (Int, Int) Vec2 -> Set ChildParent -> f [ChildParent]
sampleLoop Map (Int, Int) Vec2
_grid Set ChildParent
active
            | Set ChildParent -> Bool
forall a. Set a -> Bool
S.null Set ChildParent
active = [ChildParent] -> f [ChildParent]
forall a. a -> f a
forall (f :: * -> *) a. Applicative f => a -> f a
pure []
        sampleLoop Map (Int, Int) Vec2
grid Set ChildParent
active = do
            ChildParent Vec2
activeSample Maybe Vec2
parent <- Gen (PrimState f) -> Set ChildParent -> f ChildParent
forall (m :: * -> *) a.
PrimMonad m =>
Gen (PrimState m) -> Set a -> m a
randomSetElement Gen (PrimState m)
Gen (PrimState f)
gen Set ChildParent
active
            Vector Vec2
candidates <- Gen (PrimState f)
-> Int -> BoundingBox -> Double -> Vec2 -> f (Vector Vec2)
forall (m :: * -> *).
PrimMonad m =>
Gen (PrimState m)
-> Int -> BoundingBox -> Double -> Vec2 -> m (Vector Vec2)
candidatesAroundSample Gen (PrimState m)
Gen (PrimState f)
gen Int
k BoundingBox
bb Double
radius Vec2
activeSample

            let validPoint :: Vec2 -> Bool
validPoint Vec2
candidate =
                    let neighbours :: [Vec2]
neighbours = CellSize -> Map (Int, Int) Vec2 -> Vec2 -> [Vec2]
neighbouringPoints CellSize
cellSize Map (Int, Int) Vec2
grid Vec2
candidate
                        tooClose :: Vec2 -> Bool
tooClose Vec2
neighbour = Vec2 -> Double
normSquare (Vec2
candidate Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vec2
neighbour) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
radiusDouble -> Integer -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2
                    in Bool -> Bool
not ((Vec2 -> Bool) -> [Vec2] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any Vec2 -> Bool
tooClose [Vec2]
neighbours)

            case (Vec2 -> Bool) -> Vector Vec2 -> Maybe Vec2
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
find Vec2 -> Bool
validPoint Vector Vec2
candidates of
                Maybe Vec2
Nothing -> do
                    [ChildParent]
rest <- Map (Int, Int) Vec2 -> Set ChildParent -> f [ChildParent]
sampleLoop Map (Int, Int) Vec2
grid (ChildParent -> Set ChildParent -> Set ChildParent
forall a. Ord a => a -> Set a -> Set a
S.delete (Vec2 -> Maybe Vec2 -> ChildParent
ChildParent Vec2
activeSample Maybe Vec2
parent) Set ChildParent
active)
                    [ChildParent] -> f [ChildParent]
forall a. a -> f a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Vec2 -> Maybe Vec2 -> ChildParent
ChildParent Vec2
activeSample Maybe Vec2
parent ChildParent -> [ChildParent] -> [ChildParent]
forall a. a -> [a] -> [a]
: [ChildParent]
rest)
                Just Vec2
sample -> Map (Int, Int) Vec2 -> Set ChildParent -> f [ChildParent]
sampleLoop
                    ((Int, Int) -> Vec2 -> Map (Int, Int) Vec2 -> Map (Int, Int) Vec2
forall k a. Ord k => k -> a -> Map k a -> Map k a
M.insert (CellSize -> Vec2 -> (Int, Int)
gridCell CellSize
cellSize Vec2
sample) Vec2
sample Map (Int, Int) Vec2
grid)
                    (ChildParent -> Set ChildParent -> Set ChildParent
forall a. Ord a => a -> Set a -> Set a
S.insert (Vec2 -> Maybe Vec2 -> ChildParent
ChildParent Vec2
sample (Vec2 -> Maybe Vec2
forall a. a -> Maybe a
Just Vec2
activeSample)) Set ChildParent
active)

    [ChildParent]
edges <- Map (Int, Int) Vec2 -> Set ChildParent -> m [ChildParent]
forall {f :: * -> *}.
(PrimState f ~ PrimState m, PrimMonad f) =>
Map (Int, Int) Vec2 -> Set ChildParent -> f [ChildParent]
sampleLoop Map (Int, Int) Vec2
initialGrid Set ChildParent
initialActive
    let families :: [(Vec2, Set Vec2)]
families = (ChildParent -> Maybe (Vec2, Set Vec2))
-> [ChildParent] -> [(Vec2, Set Vec2)]
forall a b. (a -> Maybe b) -> [a] -> [b]
mapMaybe
            (\(ChildParent Vec2
child Maybe Vec2
maybeParent) -> case Maybe Vec2
maybeParent of
                Maybe Vec2
Nothing -> Maybe (Vec2, Set Vec2)
forall a. Maybe a
Nothing
                Just Vec2
parent -> (Vec2, Set Vec2) -> Maybe (Vec2, Set Vec2)
forall a. a -> Maybe a
Just (Vec2
parent, Vec2 -> Set Vec2
forall a. a -> Set a
S.singleton Vec2
child))
            [ChildParent]
edges
    Map Vec2 (Set Vec2) -> m (Map Vec2 (Set Vec2))
forall a. a -> m a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ((Set Vec2 -> Set Vec2 -> Set Vec2)
-> [(Vec2, Set Vec2)] -> Map Vec2 (Set Vec2)
forall k a. Ord k => (a -> a -> a) -> [(k, a)] -> Map k a
M.fromListWith Set Vec2 -> Set Vec2 -> Set Vec2
forall a. Ord a => Set a -> Set a -> Set a
S.union [(Vec2, Set Vec2)]
families)

data ChildParent = ChildParent Vec2 (Maybe Vec2)

instance Eq ChildParent where
    ChildParent Vec2
child1 Maybe Vec2
_parent1 == :: ChildParent -> ChildParent -> Bool
== ChildParent Vec2
child2 Maybe Vec2
_parent2 = Vec2
child1 Vec2 -> Vec2 -> Bool
forall a. Eq a => a -> a -> Bool
== Vec2
child2

instance Ord ChildParent where
    ChildParent Vec2
child1 Maybe Vec2
_parent1 compare :: ChildParent -> ChildParent -> Ordering
`compare` ChildParent Vec2
child2 Maybe Vec2
_parent2 = Vec2
child1 Vec2 -> Vec2 -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` Vec2
child2

randomSetElement :: PrimMonad m => Gen (PrimState m) -> Set a -> m a
randomSetElement :: forall (m :: * -> *) a.
PrimMonad m =>
Gen (PrimState m) -> Set a -> m a
randomSetElement Gen (PrimState m)
gen Set a
set = do
    Int
i <- (Int, Int) -> Gen (PrimState m) -> m Int
forall a g (m :: * -> *).
(UniformRange a, StatefulGen g m) =>
(a, a) -> g -> m a
forall g (m :: * -> *). StatefulGen g m => (Int, Int) -> g -> m Int
uniformRM (Int
0, Set a -> Int
forall a. Set a -> Int
S.size Set a
setInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Gen (PrimState m)
gen
    a -> m a
forall a. a -> m a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Int -> Set a -> a
forall a. Int -> Set a -> a
S.elemAt Int
i Set a
set)

-- | http://extremelearning.com.au/an-improved-version-of-bridsons-algorithm-n-for-poisson-disc-sampling/
--
-- Enhanced by a random trial order so we don’t privilege clockwise selection
candidatesAroundSample
    :: PrimMonad m
    => Gen (PrimState m)
    -> Int -- ^ Number of attempts
    -> BoundingBox -- ^ Sampling region
    -> Double -- ^ Poisson radius
    -> Vec2 -- ^ Point to generate candidates around
    -> m (Vector Vec2)
candidatesAroundSample :: forall (m :: * -> *).
PrimMonad m =>
Gen (PrimState m)
-> Int -> BoundingBox -> Double -> Vec2 -> m (Vector Vec2)
candidatesAroundSample Gen (PrimState m)
gen Int
k BoundingBox
bb Double
r Vec2
v = do
    Angle
phi0 <- Double -> Angle
rad (Double -> Angle) -> m Double -> m Angle
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (Double, Double) -> Gen (PrimState m) -> m Double
forall a g (m :: * -> *).
(UniformRange a, StatefulGen g m) =>
(a, a) -> g -> m a
forall g (m :: * -> *).
StatefulGen g m =>
(Double, Double) -> g -> m Double
uniformRM (Double
0, Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi) Gen (PrimState m)
gen
    let deltaPhi :: Angle
deltaPhi = Double -> Angle
rad (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k)
        circle :: Vector Vec2
circle = Int -> (Int -> Vec2) -> Vector Vec2
forall a. Int -> (Int -> a) -> Vector a
V.generate Int
k (\Int
i -> Vec2
v Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
+. Angle -> Double -> Vec2
polar (Angle
phi0 Angle -> Angle -> Angle
forall v. VectorSpace v => v -> v -> v
+. Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Angle -> Angle
forall v. VectorSpace v => Double -> v -> v
*. Angle
deltaPhi) (Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1e-6))
        inside :: Vector Vec2
inside = (Vec2 -> Bool) -> Vector Vec2 -> Vector Vec2
forall a. (a -> Bool) -> Vector a -> Vector a
V.filter (Vec2 -> BoundingBox -> Bool
forall thing bigObject.
(HasBoundingBox thing, HasBoundingBox bigObject) =>
thing -> bigObject -> Bool
`insideBoundingBox` BoundingBox
bb) Vector Vec2
circle
    MVector (PrimState m) Vec2
inside' <- Vector Vec2 -> m (MVector (PrimState m) Vec2)
forall (m :: * -> *) a.
PrimMonad m =>
Vector a -> m (MVector (PrimState m) a)
V.thaw Vector Vec2
inside
    Gen (PrimState m) -> MVector (PrimState m) Vec2 -> m ()
forall (f :: * -> *) a.
PrimMonad f =>
Gen (PrimState f) -> MVector (PrimState f) a -> f ()
V.fisherYatesShuffle Gen (PrimState m)
gen MVector (PrimState m) Vec2
inside'
    MVector (PrimState m) Vec2 -> m (Vector Vec2)
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> m (Vector a)
V.freeze MVector (PrimState m) Vec2
inside'

-- A cell in the grid has a side length of r/sqrt(2). If we’re somewhere in the X
-- square and can only move at most a square diagonal, we only need to consider the
-- 21 cells marked with `?`, and X itself.
--
-- +---+---+---+---+---+
-- |   | ? | ? | ? |   |
-- +---+---+---+---+---+
-- | ? | ? | ? | ? | ? |
-- +---+---+---+---+---+
-- | ? | ? | X | ? | ? |
-- +---+---+---+---+---+
-- | ? | ? | ? | ? | ? |
-- +---+---+---+---+---+
-- |   | ? | ? | ? |   |
-- +---+---+---+---+---+
neighbouringPoints :: CellSize -> Map (Int, Int) Vec2 -> Vec2 -> [Vec2]
neighbouringPoints :: CellSize -> Map (Int, Int) Vec2 -> Vec2 -> [Vec2]
neighbouringPoints CellSize
cellSize Map (Int, Int) Vec2
grid Vec2
v =
    let (Int
x, Int
y) = CellSize -> Vec2 -> (Int, Int)
gridCell CellSize
cellSize Vec2
v
    in ((Int, Int) -> Maybe Vec2) -> [(Int, Int)] -> [Vec2]
forall a b. (a -> Maybe b) -> [a] -> [b]
mapMaybe (\(Int, Int)
cell -> (Int, Int) -> Map (Int, Int) Vec2 -> Maybe Vec2
forall k a. Ord k => k -> Map k a -> Maybe a
M.lookup (Int, Int)
cell Map (Int, Int) Vec2
grid)
        [             (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1, Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2), (Int
x  , Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2), (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1, Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2)
        , (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2, Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1), (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1, Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1), (Int
x  , Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1), (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1, Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1), (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2, Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
        , (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2, Int
y  ), (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1, Int
y  ), (Int
x  , Int
y  ), (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1, Int
y  ), (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2, Int
y  )
        , (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2, Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1), (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1, Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1), (Int
x  , Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1), (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1, Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1), (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2, Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
        ,             (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1, Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2), (Int
x  , Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2), (Int
xInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1, Int
yInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2)
        ]

gridCell :: CellSize -> Vec2 -> (Int, Int)
gridCell :: CellSize -> Vec2 -> (Int, Int)
gridCell (CellSize Double
cellSize) (Vec2 Double
x Double
y) = (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double
xDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
cellSize), Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double
yDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
cellSize))
{-# INLINE gridCell #-}