sampling polygons

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

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)

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

hgeometry/src/HGeometry/Polygon/Simple/Sample.hs

@@ -0,0 +1,112 @@
+module HGeometry.Polygon.Simple.Sample
+  ( samplePolygon
+  , samplePolygons
+  ) where
+
+import           Control.Lens
+import qualified Data.Foldable as F
+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 dat a 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
+
+
+-- | Helper data type
+data Weighted w v = Weighted !w v deriving (Show)
+
+
+-- | 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) ((w,x):acc)
+                  in Sampler total (Map.fromAscList xs')
+
+-- | 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)
+
+-- | 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.