added random points with minimum distance between them

image-triangles.cabal

@@ -50,7 +50,7 @@ executable image-triangles
   main-is:             Main.hs
 
   -- Modules included in this executable, other than Main.
-  other-modules: Triangles
+  other-modules: Triangles, MinDistanceSample
 
   -- LANGUAGE extensions used by modules in this package.
   default-extensions:  ScopedTypeVariables,
@@ -76,6 +76,7 @@ executable image-triangles
   -- Other library packages from which modules are imported.
   build-depends:       base
                      , random
+                     , array
                      , hip
                      , vector-th-unbox
                      , colour

src/MinDistanceSample.hs

@@ -0,0 +1,81 @@
+module MinDistanceSample where
+
+import qualified Data.Array as A
+import qualified Data.Bifunctor as B
+import qualified Data.List as L
+import qualified Data.List.NonEmpty as NE
+import qualified Data.Maybe as My
+import Diagrams.Prelude
+import System.Random.Stateful
+import qualified Control.Monad as M
+
+k :: Int
+k = 30
+
+randomPolarCoord :: StdGen -> (V2 Double, StdGen)
+randomPolarCoord gen = (view (from r2PolarIso) (distance, angle @@ turn), gen'')
+  where
+    angle :: Double
+    (angle, gen') = uniformR (0, 1) gen
+
+    distance :: Double
+    (distance, gen'') = uniformR (1, 2) gen
+
+randomValueFrom :: NE.NonEmpty a -> StdGen -> (a, StdGen)
+randomValueFrom xs gen = B.first (xs NE.!!) (uniformR (0, pred . length $ xs) gen)
+
+randomPoints :: P2 Double -> Double -> StdGen -> [P2 Double]
+randomPoints dims minDistance gen = map (fmap (* bucketSize)) $ randomPointsRec [] (A.listArray (mkP2 0 0, fmap (ceiling . (/ bucketSize)) dims) (repeat Nothing)) gen
+  where
+
+    bucketSize :: Double
+    bucketSize = minDistance / (sqrt 2)
+
+    randomPointsRec :: [Point V2 Double] -> A.Array (P2 Word) (Maybe (P2 Double)) -> StdGen -> [P2 Double]
+    randomPointsRec [] _ _ = []
+    randomPointsRec (x : xs') grid gen =
+        maybe [] randomPoints'
+            . L.find isValidPoint
+            $ candidates
+      where
+
+        xs :: NE.NonEmpty (Point V2 Double)
+        xs = x NE.:| xs'
+
+        randomPoints' :: P2 Double -> [P2 Double]
+        randomPoints' p = (p :) $ randomPointsRec xs'' (grid A.// [(floor <$> p, Just p)]) gen''
+
+        (candidate, gen') = randomValueFrom xs gen
+
+        xs'' = NE.filter (/= candidate) xs
+
+        unitVectorsAround =
+            V2
+                <$> [-1, 0, 1]
+                <*> [-1, 0, 1]
+
+        candidates :: [P2 Double]
+        (candidates, (_, gen'') : _) =
+            B.first (map ((candidate .-^) . fst))
+                . splitAt k
+                . iterate (randomPolarCoord . snd)
+                $ (unitX, gen')
+
+        newPoint :: Maybe (P2 Double)
+        newPoint = L.find isValidPoint candidates
+
+        isValidPoint :: Point V2 Double -> Bool
+        isValidPoint p =
+            all ((>= 1) . abs . norm . (p .-.))
+                . My.mapMaybe
+                    ( M.join
+                        . (grid !?)
+                        . fmap floor
+                        . (p .-^)
+                    )
+                $ unitVectorsAround
+
+arr !? index | index >= minBound && (index <= maxBound) = Just (arr A.! index)
+  where
+    (minBound, maxBound) = A.bounds arr
+_ !? _ = Nothing