diff --git a/hgeometry-showcase/hgeometry-showcase.cabal b/hgeometry-showcase/hgeometry-showcase.cabal
index cb77dba92..6ee03fbec 100644
--- a/hgeometry-showcase/hgeometry-showcase.cabal
+++ b/hgeometry-showcase/hgeometry-showcase.cabal
@@ -45,7 +45,7 @@ executable showcase
                        RandomMonotone,
                        ZHashing
                        SvgPolygons
-                       -- UniformSampling
+                       UniformSampling
   -- other-extensions:
   build-depends:       base,
                        reanimate >= 1.1.3.0, reanimate-svg,
diff --git a/hgeometry-showcase/lib/Common.hs b/hgeometry-showcase/lib/Common.hs
index a5a3f3804..6c8434fae 100644
--- a/hgeometry-showcase/lib/Common.hs
+++ b/hgeometry-showcase/lib/Common.hs
@@ -2,6 +2,7 @@
 {-# LANGUAGE RankNTypes          #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeOperators       #-}
+{-# LANGUAGE PatternSynonyms #-}
 module Common where
 
 import           Codec.Picture             (PixelRGBA8 (..))
@@ -27,7 +28,7 @@ import Data.Ext
 import Data.Geometry.LineSegment
 import Data.Geometry.Point
 import Data.Geometry.Polygon.Bezier
-import Data.Geometry.BezierSpline
+import Data.Geometry.BezierSpline (pattern Bezier3, pattern Bezier2, quadToCubic)
 import Data.Geometry.Polygon
 import Data.Geometry.Box (Rectangle, box, minPoint, maxPoint)
 import Data.Geometry.Transformation
@@ -263,6 +264,12 @@ type SomeBezierPolygon r = SomePolygon (PathJoin r) r
 
 svgToPolygons :: SVG -> [SomeBezierPolygon Double]
 svgToPolygons = plGroupShapes . cmdsToBezier . toLineCommands . extractPath
+-- svgToPolygons svg =
+--   [ case cmdsToBezier $ toLineCommands $ extractPath (ctx glyph) of
+--       [] -> error "No polygon in glyph."
+--       [x] -> Left x
+--       (x:hs) -> Right (MultiPolygon x hs)
+--   | (ctx, _attr, glyph) <- svgGlyphs svg ]
 
 cmdsToBezier :: [LineCommand] -> [SimpleBezierPolygon Double]
 cmdsToBezier [] = []
@@ -320,7 +327,9 @@ plGroupShapes = worker
     -- This is super fragile. I wish there was a way to check for bezier-line intersections without using sqrt.
     isInsideOf :: SimpleBezierPolygon Double -> SimpleBezierPolygon Double -> Bool
     isInsideOf a b =
-      _core (CV.head (a^.outerBoundaryVector)) `insidePolygon` b
+      CV.all (\elt -> _core elt `insidePolygon` approximate eps b) ((approximate eps a)^.outerBoundaryVector)
+    eps = 0.001
+      -- _core (CV.head (a^.outerBoundaryVector)) `insidePolygon` b
 
 ------------------------------------------------------------------
 -- Random data
diff --git a/hgeometry-showcase/src/Showcase.hs b/hgeometry-showcase/src/Showcase.hs
index 9126ddbb2..c7c219ecb 100644
--- a/hgeometry-showcase/src/Showcase.hs
+++ b/hgeometry-showcase/src/Showcase.hs
@@ -12,7 +12,7 @@ import FastVisibility
 import SSSP
 import RandomMonotone
 import ZHashing
--- import UniformSampling
+import UniformSampling
 import SvgPolygons
 
 import Reanimate
@@ -28,7 +28,7 @@ main = do
     "fast_visibility":rest -> withArgs rest $ runReanimate fastVisibilityShowcase
     "random_monotone":rest -> withArgs rest $ runReanimate randomMonotoneShowcase
     "zhashing":rest        -> withArgs rest $ runReanimate zHashingShowcase
-    -- "sampling":rest        -> withArgs rest $ runReanimate uniformSamplingShowcase
+    "sampling":rest        -> withArgs rest $ runReanimate uniformSamplingShowcase
     "svg":rest             -> withArgs rest $ runReanimate svgPolygonsShowcase
     _                      -> printUsage
 
diff --git a/hgeometry-showcase/src/UniformSampling.hs b/hgeometry-showcase/src/UniformSampling.hs
new file mode 100644
index 000000000..4da31bb49
--- /dev/null
+++ b/hgeometry-showcase/src/UniformSampling.hs
@@ -0,0 +1,52 @@
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE OverloadedStrings #-}
+module UniformSampling (uniformSamplingShowcase) where
+
+import           Data.Geometry.Polygon.Sample
+import           Common
+import           Control.Monad.Random                       (evalRand, forM_, mkStdGen, replicateM)
+import           Data.Geometry.Point                        (Point (Point2))
+import           Data.Geometry.Polygon
+import           Data.Geometry.Polygon.Bezier               (approximateSome)
+import           Reanimate
+
+eps :: Double
+eps = 0.001
+
+svgPolygons :: [SomePolygon () Double]
+svgPolygons = map (approximateSome eps) $
+  svgToPolygons $ lowerTransformations $ scale 2 $ center $ latex "Uniform\n\nSampling"
+
+-- targetPolygon :: SomePolygon () Double
+-- targetPolygon = scaleUniformlyBy 2 $ pAtCenter $ simpleFromPoints $ map ext
+--   [ 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 ]
+
+uniformSamplingShowcase :: Animation
+uniformSamplingShowcase = scene $ do
+  forM_ svgPolygons $ \svgPolygon -> fork $ do 
+    newSpriteSVG_ $ overAny (ppPolygonBody' grey) svgPolygon
+    -- newSpriteSVG_ $ overAny (ppPolygonOutline' black) targetPolygon
+    wait 1
+    forM_ (overAny points svgPolygon) $ \point -> do
+      fork $ play $ playThenReverseA $ mkAnimation 0.5 $ \t ->
+        aroundCenter (scale t) $ ppPoint point
+      wait 0.025
+
+nPoints :: Int
+nPoints = 100
+
+points :: (Show p) => Polygon t p Double -> [Point 2 Double]
+points p = flip evalRand (mkStdGen seed) $ replicateM nPoints $ samplePolygon p
+
+seed :: Int
+seed = 0xDEADBEEF
+
+ppPoint :: Real r => Point 2 r -> SVG
+ppPoint (Point2 x y) =
+  withFillColorPixel red $
+  translate (realToFrac x) (realToFrac y) $
+  mkCircle (nodeRadius/5)
diff --git a/hgeometry/hgeometry.cabal b/hgeometry/hgeometry.cabal
index c443a33ba..3f1113598 100644
--- a/hgeometry/hgeometry.cabal
+++ b/hgeometry/hgeometry.cabal
@@ -85,6 +85,7 @@ library
                     Data.Geometry.Ellipse
 
                     Data.Geometry.Polygon
+                    Data.Geometry.Polygon.Sample
                     Data.Geometry.Polygon.Bezier
                     Data.Geometry.Polygon.Inflate
                     Data.Geometry.Polygon.Convex
diff --git a/hgeometry/src/Data/Geometry/Polygon/Sample.hs b/hgeometry/src/Data/Geometry/Polygon/Sample.hs
new file mode 100644
index 000000000..e873046f3
--- /dev/null
+++ b/hgeometry/src/Data/Geometry/Polygon/Sample.hs
@@ -0,0 +1,50 @@
+module Data.Geometry.Polygon.Sample where
+
+import           Algorithms.Geometry.PolygonTriangulation.Triangulate
+import           Control.Lens
+import           Control.Monad.Random
+import           Data.Ext
+import           Data.Geometry.PlanarSubdivision                      (PolygonFaceData (..))
+import           Data.Geometry.Point
+import           Data.Geometry.Polygon.Core                           as Polygon
+import           Data.Geometry.Triangle                               as Triangle
+import qualified Data.List.NonEmpty                                   as NonEmpty
+import           Data.PlaneGraph
+import           Data.Proxy
+import qualified Data.Vector                                          as V
+import           Linear.Affine                                        hiding (Point)
+import           Linear.Vector
+
+-- | Uniformly samples from a set of polygon in \(O(n \log n)\) triangles.
+samplePolygons :: (RandomGen g, Random r, Fractional r, Ord r, Real r) => NonEmpty.NonEmpty (Polygon t p r) -> Rand g (Point 2 r)
+samplePolygons ps = do
+   randTri <- fromList $ map (\tri -> (tri, toRational $ areaRatio tri)) $ concatMap toTriangles ps
+   sampleTriangle randTri
+   where
+      areaRatio tri = Triangle.area tri / sum (fmap Polygon.area ps)
+
+-- | Uniformly samples a polygon in \(O(n \log n)\).
+samplePolygon :: (RandomGen g, Random r, Fractional r, Ord r, Real r) => Polygon t p r -> Rand g (Point 2 r)
+samplePolygon p = samplePolygons $ p NonEmpty.:| []
+
+-- | Uniformly samples a triangle in \(O(1)\).
+sampleTriangle :: (RandomGen g, Random r, Fractional r, Ord r) => Triangle 2 p r -> Rand g (Point 2 r)
+sampleTriangle (Triangle v1 v2 v3) = do
+  a' <- getRandomR (0, 1)
+  b' <- getRandomR (0, 1)
+  let (a, b) = if a' + b' > 1 then (1 - a', 1 - b') else (a', b')
+  return $ v1^.core .+^ a*^u .+^ b*^v
+  where
+    u = v2^.core .-. v1^.core
+    v = v3^.core .-. v1^.core
+
+-- | Triangulates a polygon \(O(n \log n)\).
+toTriangles :: (Fractional r, Ord r) => Polygon t p r -> [Triangle 2 p r]
+toTriangles p =
+    [ (polygonToTriangle . view core . (`rawFacePolygon` g)) f
+    | (f, Inside) <- V.toList (faces g) ]
+  where
+    g = triangulate' Proxy p
+    polygonToTriangle poly = case NonEmpty.toList $ polygonVertices poly of
+      [v1, v2, v3] -> Triangle v1 v2 v3
+      _            -> error "Invalid Triangulation"
diff --git a/hie.yaml b/hie.yaml
index f00d434de..2e130e144 100644
--- a/hie.yaml
+++ b/hie.yaml
@@ -135,6 +135,9 @@ cradle:
     - path: "hgeometry-showcase/src/RandomMonotone.hs"
       component: "hgeometry-showcase:exe:showcase"
 
+    - path: "hgeometry-showcase/src/UniformSampling.hs"
+      component: "hgeometry-showcase:exe:showcase"
+
     - path: "hgeometry-showcase/lib/Common.hs"
       component: "hgeometry-showcase:exe:showcase"