From f0756e2b6b4ef8c898506ee55f66cad2b018e4ec Mon Sep 17 00:00:00 2001
From: Frank Staals <991345+noinia@users.noreply.github.com>
Date: Fri, 16 Aug 2024 21:00:55 +0200
Subject: [PATCH] Sampling points in polygons (adapted/ported from #129) (#238)

* sampling polygons

* sampling polygons

* some additional documentation

* some cleaning up of the imports

* no longer load semigroupids?
---
 hgeometry-examples/hgeometry-examples.cabal   |  13 ++
 hgeometry-examples/sampler/Main.hs            |  48 +++++++
 hgeometry/doctests.hs                         |   1 -
 hgeometry/hgeometry.cabal                     |   2 +
 hgeometry/kernel/src/HGeometry/Triangle.hs    |  11 ++
 hgeometry/src/HGeometry/Polygon.hs            |  26 ++++
 .../src/HGeometry/Polygon/Simple/Sample.hs    | 125 ++++++++++++++++++
 7 files changed, 225 insertions(+), 1 deletion(-)
 create mode 100644 hgeometry-examples/sampler/Main.hs
 create mode 100644 hgeometry/src/HGeometry/Polygon.hs
 create mode 100644 hgeometry/src/HGeometry/Polygon/Simple/Sample.hs

diff --git a/hgeometry-examples/hgeometry-examples.cabal b/hgeometry-examples/hgeometry-examples.cabal
index 07d67e006..aa525b481 100644
--- a/hgeometry-examples/hgeometry-examples.cabal
+++ b/hgeometry-examples/hgeometry-examples.cabal
@@ -204,6 +204,19 @@ executable hgeometry-geojson
           Paths_hgeometry_examples
   --      Miso.Event.Extra
 
+--------------------------------------------------------------------------------
+-- * Polygon Sampler example
+
+executable hgeometry-sampler
+  import:         setup, miso-setup
+  hs-source-dirs: sampler
+  main-is:        Main.hs
+  other-modules:
+          -- Paths_hgeometry_examples
+  --      Miso.Event.Extra
+
+
+
 --------------------------------------------------------------------------------
 -- * Polyline Drawing
 
diff --git a/hgeometry-examples/sampler/Main.hs b/hgeometry-examples/sampler/Main.hs
new file mode 100644
index 000000000..979eb8744
--- /dev/null
+++ b/hgeometry-examples/sampler/Main.hs
@@ -0,0 +1,48 @@
+{-# LANGUAGE QuasiQuotes #-}
+module Main(main) where
+
+import Control.Monad (replicateM)
+import Data.List.NonEmpty (NonEmpty(..))
+import Data.Maybe
+import HGeometry.Point
+import HGeometry.Polygon.Simple
+import HGeometry.Polygon.Simple.Sample
+import HGeometry.Transformation
+import HGeometry.Triangle
+import Ipe
+import System.OsPath
+import System.Random.Stateful
+--------------------------------------------------------------------------------
+
+type R = Double
+
+targetPolygon :: SimplePolygon (Point 2 R)
+targetPolygon = scaleUniformlyBy 20 $ fromJust $ fromPoints
+  [ Point2 0 0
+  , Point2 1 0
+  , Point2 1 1, Point2 2 1, Point2 2 (-1)
+  , Point2 0 (-1), Point2 0 (-2)
+  , Point2 3 (-2), Point2 3 2, Point2 0 2
+  ]
+
+sampler :: Sampler R (Triangle (Point 2 R))
+sampler = triangleSampler $ targetPolygon :| []
+
+numPoints :: Int
+numPoints = 1000
+
+samples   :: StatefulGen g IO => g -> IO [Point 2 Double]
+samples g = replicateM numPoints (samplePoint sampler g)
+
+-- | Example that samples 1000 points inside some polygon, and writes the ouput to an ipe
+-- file.
+main :: IO ()
+main = do
+  pts <- samples globalStdGen
+  let outFp = [osp|foo.ipe|]
+      out   = [ iO $ defIO targetPolygon ]
+               <>
+              [ iO $ defIO p
+              | p <- pts
+              ]
+  writeIpeFile outFp . singlePageFromContent $ out
diff --git a/hgeometry/doctests.hs b/hgeometry/doctests.hs
index f244e1f81..77fd2d263 100644
--- a/hgeometry/doctests.hs
+++ b/hgeometry/doctests.hs
@@ -29,7 +29,6 @@
                 -optF -XDerivingStrategies
                 -optF -XDerivingVia
 
-                -optF -package=semigroupoids
                 -optF -package=hgeometry-combinatorial
                 -optF -package=hiraffe
                 -optF -package=hgeometry
diff --git a/hgeometry/hgeometry.cabal b/hgeometry/hgeometry.cabal
index cfb850b6a..d3035b020 100644
--- a/hgeometry/hgeometry.cabal
+++ b/hgeometry/hgeometry.cabal
@@ -345,10 +345,12 @@ library
         HGeometry.PlaneGraph
         HGeometry.PlaneGraph.Class
 
+        HGeometry.Polygon
         HGeometry.Polygon.Class
 
         HGeometry.Polygon.Simple
         HGeometry.Polygon.Simple.Class
+        HGeometry.Polygon.Simple.Sample
 --
         HGeometry.Polygon.Convex
         HGeometry.Polygon.Convex.Class
diff --git a/hgeometry/kernel/src/HGeometry/Triangle.hs b/hgeometry/kernel/src/HGeometry/Triangle.hs
index 4c541ab79..03f671d32 100644
--- a/hgeometry/kernel/src/HGeometry/Triangle.hs
+++ b/hgeometry/kernel/src/HGeometry/Triangle.hs
@@ -15,6 +15,7 @@ module HGeometry.Triangle
   ) where
 
 import Control.Lens
+import Data.Foldable1
 import GHC.Generics (Generic)
 import HGeometry.Box.Boxable
 import HGeometry.Intersection
@@ -31,12 +32,22 @@ import Text.Read
 -- | Triangles in d-dimensional space
 newtype Triangle point = MkTriangle (Vector 3 point)
   deriving (Generic)
+  deriving newtype (Functor,Foldable)
 
 -- | Construct a triangle from its three points
 pattern Triangle       :: point -> point -> point -> Triangle point
 pattern Triangle a b c = MkTriangle (Vector3 a b c)
 {-# COMPLETE Triangle #-}
 
+instance Traversable Triangle where
+  traverse f (MkTriangle v) = MkTriangle <$> traverse f v
+
+instance Foldable1 Triangle where
+  foldMap1 f (MkTriangle v) = foldMap1 f v
+
+instance Traversable1 Triangle where
+  traverse1 f (MkTriangle v) = MkTriangle <$> traverse1 f v
+
 
 deriving instance Eq  (Vector 3 point) => Eq (Triangle point)
 deriving instance Ord (Vector 3 point) => Ord (Triangle point)
diff --git a/hgeometry/src/HGeometry/Polygon.hs b/hgeometry/src/HGeometry/Polygon.hs
new file mode 100644
index 000000000..e729337b1
--- /dev/null
+++ b/hgeometry/src/HGeometry/Polygon.hs
@@ -0,0 +1,26 @@
+--------------------------------------------------------------------------------
+-- |
+-- Module      :  HGeometry.Polygon
+-- Copyright   :  (C) Frank Staals
+-- License     :  see the LICENSE file
+-- Maintainer  :  Frank Staals
+--
+-- A polygon and some basic functions to interact with them.
+--
+--------------------------------------------------------------------------------
+module HGeometry.Polygon
+  ( module HGeometry.Polygon.Class
+  , asTriangle
+  ) where
+
+import Control.Lens
+import HGeometry.Polygon.Class
+import HGeometry.Triangle
+
+--------------------------------------------------------------------------------
+
+-- | Try to convert the polygon into a triangle (whose vertices are given in CCW order).
+asTriangle    :: Polygon_ polygon point r => polygon -> Maybe (Triangle point)
+asTriangle pg = case pg^..vertices of
+                  [u,v,w] -> Just $ Triangle u v w
+                  _       -> Nothing
diff --git a/hgeometry/src/HGeometry/Polygon/Simple/Sample.hs b/hgeometry/src/HGeometry/Polygon/Simple/Sample.hs
new file mode 100644
index 000000000..e0daa1e01
--- /dev/null
+++ b/hgeometry/src/HGeometry/Polygon/Simple/Sample.hs
@@ -0,0 +1,125 @@
+--------------------------------------------------------------------------------
+-- |
+-- Module      :  HGeometry.Polygon.Simple.Sample
+-- Copyright   :  (C) Frank Staals, Owen Graves,  David Himmelstrup
+-- License     :  see the LICENSE file
+-- Maintainer  :  Frank Staals
+--
+-- Functionality to sample points uniformly at random from within a simple polygon.
+--
+--------------------------------------------------------------------------------
+module HGeometry.Polygon.Simple.Sample
+  ( samplePolygon
+  , samplePolygons
+
+  , Sampler
+  , samplePoint
+  , triangleSampler
+  ) where
+
+import           Control.Lens
+import           Data.Foldable1
+import qualified Data.List.NonEmpty as NonEmpty
+import qualified Data.Map as Map
+import           Data.Maybe (mapMaybe)
+import           HGeometry.Ext
+import           HGeometry.PlaneGraph
+import           HGeometry.Point
+import           HGeometry.Polygon
+import           HGeometry.Polygon.Simple.Class
+import           HGeometry.Polygon.Triangulation
+import           HGeometry.Triangle
+import           HGeometry.Vector
+import           System.Random.Stateful
+
+--------------------------------------------------------------------------------
+
+-- | A data structure that can be used to efficiently sample values of type v.
+data Sampler w v = Sampler !w -- ^ the total weight
+                           (Map.Map w v)
+  deriving (Show,Read,Eq,Functor,Foldable)
+
+instance Traversable (Sampler w) where
+  traverse f (Sampler w m) = Sampler w <$> traverse f m
+
+-- | Build a sampler
+--
+-- O(n)
+buildSampler    :: (Foldable1 nonEmpty, Num w, Ord w) => nonEmpty (w, v) -> Sampler w v
+buildSampler xs = let Weighted total xs'       = foldr f (Weighted 0 []) xs
+                      f (w,x) (Weighted t acc) = Weighted (w+t) ((t,x):acc)
+                  in Sampler total (Map.fromDescList xs')
+
+-- | Helper data type
+data Weighted w v = Weighted !w v deriving (Show)
+
+-- | Sample a value from the sampler
+--
+-- \(O(\log n)\)
+sample                     :: (StatefulGen g m, Ord w, Num w, UniformRange w)
+                           => Sampler w v -> g -> m v
+sample (Sampler total m) g = (maybe err snd . flip Map.lookupLE m) <$> uniformRM (0, total) g
+  where err = error "sample: absurd."
+
+
+--------------------------------------------------------------------------------
+
+
+-- | Build a triangle sampler; i.e. samples a triagnle from the polygons
+-- with probability proportional to its area.
+--
+-- \(O(N\log N\), where \(N\) is the total size of all polygons.
+triangleSampler :: (SimplePolygon_ polygon point r, Num r, Ord r, Foldable1 nonEmpty)
+                => nonEmpty polygon -> Sampler r (Triangle point)
+triangleSampler = buildSampler . fmap (\tri -> (triangleSignedArea2X tri, tri))
+                . foldMap1 toTriangles
+
+-- | Sample a point
+samplePoint           :: (Point_ point 2 r, StatefulGen g m, Real r, Ord r, UniformRange r)
+                      => Sampler r (Triangle point) -> g -> m (Point 2 Double)
+samplePoint sampler g = sample sampler g >>= flip sampleTriangle g
+
+-- | Unfiormly samples a point from a set of polygons. You may want to build a
+-- pointSampler if the goal is to sample multiple points.
+--
+-- \(O(N\log N\), where \(N\) is the total size of all polygons.
+samplePolygons       :: ( SimplePolygon_ polygon point r, StatefulGen g m
+                        , Foldable1 nonEmpty, Real r, Ord r, UniformRange r
+                        )
+                     => nonEmpty polygon -> g -> m (Point 2 Double)
+samplePolygons pgs g = flip samplePoint g $ triangleSampler pgs
+
+-- | Uniformly samples a point in a polygon.
+--
+-- \(O(n\log n )\)
+samplePolygon   :: ( SimplePolygon_ polygon point r, Ord r, Real r, UniformRange r
+                   , StatefulGen g m)
+                => polygon -> g -> m (Point 2 Double)
+samplePolygon p = samplePolygons $ p NonEmpty.:| []
+
+-- | Uniformly samples a triangle in \(O(1)\) time.
+sampleTriangle                                       :: ( Point_ point 2 r, Real r
+                                                       , StatefulGen g m)
+                                                     => Triangle point -> g -> m (Point 2 Double)
+sampleTriangle (fmap doubleP -> Triangle v1 v2 v3) g =
+    f <$> uniformRM (0, 1) g <*> uniformRM (0, 1) g
+  where
+    f       :: Double -> Double -> Point 2 Double
+    f a' b' = let (a, b) = if a' + b' > 1 then (1 - a', 1 - b') else (a', b')
+                  u = v2 .-. v1
+                  v = v3 .-. v1
+              in v1 .+^ (a*^u) .+^ (b*^v)
+-- the idea is based on
+-- https://blogs.sas.com/content/iml/2020/10/19/random-points-in-triangle.html
+
+-- | Convert to a point with double coordiantes
+doubleP :: (Point_ point 2 r, Real r) => point -> Point 2 Double
+doubleP = over coordinates realToFrac . view asPoint
+
+-- | Triangulates a polygon \(O(n \log n)\).
+toTriangles :: (SimplePolygon_ polygon point r, Num r, Ord r)
+            => polygon -> NonEmpty.NonEmpty (Triangle point)
+toTriangles = NonEmpty.fromList . fmap (fmap (view core))
+            . mapMaybe asTriangle . toListOf interiorFacePolygons . triangulate @()
+  -- any valid simple polygon produces at least one triangle, so the
+  -- NonEmpty.fromList is safe.