minimum distance sampling

https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf

That paper provided a wonderful way to do minimum distance sampling.

.gitignore

@@ -1,3 +1,4 @@
 /dist-newstyle/
 /result
 /.direnv/
+*/.DS_Store

README.md

@@ -1,14 +1,14 @@
 # image-triangles
 
-Makes a [voroni diagram](https://en.wikipedia.org/wiki/Voronoi_diagram) and fills each cell with the average color the image below in. Please note, there's currently a bug. Some cells will be blank.
+Makes a [voroni diagram](https://en.wikipedia.org/wiki/Voronoi_diagram) and fills each cell with the average color the image below in.
 
 ### examples
 ![Sierra mountains original](examples/sierra.jpg)
 ![Sierra mountains post-filter](examples/sierra-result.svg)
-![Hawaii original](examples/birds-eye-view.webp)
+![Hawaii original](examples/birds-eye-view.png)
 ![Hawaii post-filter](examples/birds-eye-view.svg)
-![Dog original](examples/luna-but-square.jpeg)
-![Dog post-filter](examples/luna-but-square.svg)
+![Dog original](examples/luna.jpeg)
+![Dog post-filter](examples/luna.svg)
 
 
 ### to run:
@@ -16,7 +16,7 @@ Makes a [voroni diagram](https://en.wikipedia.org/wiki/Voronoi_diagram) and fill
 #### with nix
 
 ```
-nix run  --experimental-features 'nix-command flakes' -- --cornerCount 800 --input examples/birds-eye-view.webp --output output.svg
+nix run  --experimental-features 'nix-command flakes' -- --minDistance 0.05 --input examples/birds-eye-view.png --output output.svg
 ```
 
 #### with cabal
@@ -24,6 +24,6 @@ Install [cabal & ghc](https://www.haskell.org/ghcup/) if you don't have them.
 
 ```
 cabal update
-cabal run image-triangles -- --cornerCount 800 --input examples/birds-eye-view.webp --output output.svg
+cabal run image-triangles -- --minDistance 0.05 --input examples/birds-eye-view.png --output output.svg
 ```
 

src/Main.hs

@@ -33,18 +33,10 @@ import System.Random.Internal
 import System.Random.SplitMix
 import Triangles (getTriangleAverageRGB)
 import qualified Triangles as Tri
-
-data Options = Options
-  { numPoints :: Int
-  , gen :: Maybe StdGen
-  }
-
--- -- modify this to your liking
-defaultOpts =
-  Options
-    { numPoints = 10
-    , gen = Nothing
-    }
+import qualified MinDistanceSample as MDS
+import Control.Monad.Zip (MonadZip(mzipWith))
+import qualified Data.Bifunctor as Bi
+import qualified Debug.Trace as D
 
 -- CL.rgb might be the wrong fn...
 tosRGB' :: (Ord b, Floating b) => Pixel G.RGB b -> CL.Colour b
@@ -58,12 +50,11 @@ corners = (,) <$> [0, 1] <*> [0, 1]
 shapeCircumference :: [Point V2 Double] -> Double
 shapeCircumference = Data.List.sum . map D.norm . loopOffsets . fromVertices
 
-genImage' :: Image VU G.RGB Double -> StdGen -> Int -> QDiagram SVG V2 Double Any
-genImage' image gen cornerCount =
-  scaleY widthHeightRatio
+genImage :: Image VU G.RGB Double -> Double -> StdGen  -> QDiagram SVG V2 Double Any
+genImage image minDistance gen  =
+  scaleX widthHeightRatio
     . reflectY
     . rectEnvelope (mkP2 0 0) (1 ^& 1)
-    -- . atop visualizeGraph
     . mconcat
     . map drawVoroniRegion
     . sortOn shapeCircumference
@@ -79,10 +70,11 @@ genImage' image gen cornerCount =
       $ shape
 
   widthHeightRatio :: Double
-  widthHeightRatio = (fromIntegral . snd $ dimensions) / (fromIntegral . fst $ dimensions)
+  widthHeightRatio = (fromIntegral . fst $ dimensions) / (fromIntegral . snd $ dimensions)
 
   img' = convImage image
   dimensions = uncurry (flip (,)) . Img.dims $ image
+  dimensionsVec = fromIntegral <$> uncurry V2 dimensions
 
   singleVoroni = last voroni
 
@@ -100,7 +92,13 @@ genImage' image gen cornerCount =
       . Tri.toPlanarGraph
       $ corners'
 
-  corners' = take cornerCount . map p2 $ corners ++ Tri.randomPoints gen
+
+  averageSideSize = (fromIntegral (uncurry (+) dimensions)) / 2
+
+  padding = (/) 10 . (*) widthHeightRatio <$> V2 1 1
+
+  corners' :: [P2 Double]
+  corners' = map (p2 . Bi.first (/ widthHeightRatio) . unp2 . (.-^ ((/ 2) <$> padding))) . MDS.randomPoints ((mkP2 widthHeightRatio 1) .+^ padding) ( minDistance * widthHeightRatio) $ gen
 
 deriving instance Generic (CL.RGB a)
 deriving instance NFData a => NFData (CL.RGB a)
@@ -113,7 +111,7 @@ toDimensionVector image =
 data CLIOptions = CLIOptions
   { input :: FilePath
   , output :: FilePath
-  , cornerCount :: Int
+  , minDistance :: Double
   }
   deriving (Generic)
 
@@ -126,4 +124,4 @@ main = do
   print gen'
   image <- Img.readImageRGB VU input
   let dimVector = toDimensionVector image
-  renderSVG output dimVector (genImage' image gen' cornerCount)
+  renderSVG output dimVector (genImage image minDistance gen')

src/MinDistanceSample.hs

@@ -1,81 +1,73 @@
 module MinDistanceSample where
 
+import qualified Control.Monad as M
 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 qualified Debug.Trace as D
 import Diagrams.Prelude
 import System.Random.Stateful
-import qualified Control.Monad as M
+import qualified Data.Map.Strict as M
+import qualified Data.Ix as Ix
 
 k :: Int
-k = 30
+k = 10
 
-randomPolarCoord :: StdGen -> (V2 Double, StdGen)
-randomPolarCoord gen = (view (from r2PolarIso) (distance, angle @@ turn), gen'')
+randomPolarCoord :: (StatefulGen g m, UniformRange n, RealFloat n) => g -> m (V2 n)
+randomPolarCoord gen = do
+    angle <- uniformRM (0, 1) gen
+    distance <- uniformRM (1, 2) gen
+    pure $ view (from r2PolarIso) (distance, angle @@ turn)
+
+randomPoints :: RandomGen r => Point V2 Double -> Double -> r -> [Point V2 Double]
+randomPoints dims minDistance gen' = runStateGen_ gen' randomPointsM
   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
+    randomPointsM gen = do
+        startingPoint <- randomPoint
+        map (fmap (* bucketSize)) <$> randomPointsRec [startingPoint] (initialGrid startingPoint)
       where
+        initialGrid = addPointToGrid M.empty
+
+        gridBounds = (mkP2 0 0, fmap (floor . (/ bucketSize)) dims)
+
+        randomPoint = p2 <$> randomRM ((0, 0), unp2 . fmap fromIntegral . snd $ gridBounds) gen
 
         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
+        addPointToGrid grid p = M.insert (floor <$> p) (p) grid
+
+        randomValueFrom xs = (xs NE.!!) <$> uniformRM (0, pred . length $ xs) gen
+
+        randomPointsRec [] grid = pure . M.elems $ grid
+        randomPointsRec (x : xs') grid = do
+            startingPoint <- randomValueFrom xs
+            newPoint <- L.find isValidPoint <$> candidates startingPoint
+            randomPointsRec (prependIfJust newPoint $ (if My.isJust newPoint then x : xs' else xsWithout startingPoint)) . gridPlus $ newPoint
           where
+            prependIfJust (Just a) xs = a : xs
+            prependIfJust Nothing xs = xs
+
+            gridPlus = maybe grid (addPointToGrid grid)
 
             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
+            xsWithout x = NE.filter (/= x) 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
+            candidates startingPoint = do
+                candidatePoints <- M.replicateM k . randomPolarCoord $ gen
+                pure $ map (startingPoint .-^) candidatePoints
 
             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
+                (Ix.inRange gridBounds (floor <$> p)) && (all ((>= 1) . abs . norm . (p .-.))
+                    . My.mapMaybe ((grid M.!?) . fmap floor . (p .-^))
+                    $ unitVectorsAround)

src/Triangles.hs

@@ -65,7 +65,7 @@ combinations xs =
     . filter (uncurry (/=))
     . concat
     . withStrategy (parListChunk 50 rdeepseq)
-    . map (\(x : xs) -> take 15 . sortOn (abs . uncurry distanceA) . map (x,) $ xs)
+    . map (\(x : xs) -> take 10 . sortOn (abs . uncurry distanceA) . map (x,) $ xs)
     . init -- last output of tails is empty list
     . tails
     $ xs