Implement uniform random point in polygon, triangle and circle

examples/svg/uniform_distribution.clj

@@ -0,0 +1,90 @@
+(ns thi.ng.geom.examples.svg.uniform-distribution
+  (:require [thi.ng.geom.circle :as c]
+            [thi.ng.geom.core :as g]
+            [thi.ng.geom.polygon :as p]
+            [thi.ng.geom.rect :as r]
+            [thi.ng.geom.svg.adapter :as adapt]
+            [thi.ng.geom.svg.core :as svg]
+            [thi.ng.geom.triangle :as t]
+            [thi.ng.geom.vector :as v]))
+
+(defn sample-points
+  ([shape sample-method]
+   (sample-points shape sample-method 400))
+  ([shape sample-method n]
+   (let [centered (g/center shape)]
+     (repeatedly n #(sample-method centered)))))
+
+(defn example [pos shape points description]
+  (let [shape-centered (g/center shape)]
+    (svg/group {}
+               (svg/text (g/translate pos (v/vec2 0 -70))
+                         description
+                         {:text-anchor "middle"})
+               (svg/group {:fill "none" :stroke "red"} (g/translate shape-centered pos))
+               (svg/group {:opacity 0.8}
+                          (for [[x y] points]
+                            (g/translate (c/circle x y 0.5) pos))))))
+
+(defn scene []
+  (let [circle (c/circle 0 0 50)
+        ;; ellipse is not implemented?
+        rectangle (r/rect 0 0 100 100)
+        rotated-rectangle (g/rotate rectangle 0.25)
+        triangle (t/triangle2 (v/vec2 0 0) (v/vec2 100 50) (v/vec2 25 100))
+        polygon (p/polygon2 [0 0] [50 75] [100 100] [100 50] [75 25])]
+    (svg/svg {:width 800 :height 600 :stroke "black"}
+             (example (v/vec2 100 100) circle
+                      (sample-points circle g/random-point-inside)
+                      "g/random-point-inside")
+             (example (v/vec2 300 100) circle
+                      (sample-points circle c/random-point-in-circle)
+                      "c/random-point-in-circle")
+             (example (v/vec2 500 100) circle
+                      (sample-points circle g/random-point 200)
+                      "g/random-point")
+             (example (v/vec2 700 100) circle
+                      (g/sample-uniform (g/center circle) 10 true)
+                      "g/sample-uniform")
+
+             (example (v/vec2 100 250) rectangle
+                      (sample-points rectangle g/random-point-inside)
+                      "g/random-point-inside")
+             (example (v/vec2 300 250) rotated-rectangle
+                      (sample-points rotated-rectangle g/random-point-inside)
+                      "g/random-point-inside (polygon)")
+             (example (v/vec2 500 250) rectangle
+                      (sample-points rectangle g/random-point 200)
+                      "g/random-point")
+             (example (v/vec2 700 250) rectangle
+                      (g/sample-uniform (g/center rectangle) 10 true)
+                      "g/sample-uniform")
+
+             (example (v/vec2 100 400) triangle
+                      (sample-points triangle g/random-point-inside)
+                      "g/random-point-inside")
+             (example (v/vec2 300 400) triangle
+                      (sample-points triangle t/random-point-in-triangle2)
+                      "t/random-point-in-triangle2")
+             (example (v/vec2 500 400) triangle
+                      (sample-points triangle g/random-point 200)
+                      "g/random-point")
+             (example (v/vec2 700 400) triangle
+                      (g/sample-uniform (g/center triangle) 10 true)
+                      "g/sample-uniform")
+
+             (example (v/vec2 300 550) polygon
+                      (sample-points polygon g/random-point-inside)
+                      "g/random-point-inside")
+             (example (v/vec2 500 550) polygon
+                      (sample-points polygon g/random-point 200)
+                      "g/random-point")
+             (example (v/vec2 700 550) polygon
+                      (g/sample-uniform (g/center polygon) 10 true)
+                      "g/sample-uniform")
+             )))
+
+(->> (scene)
+     adapt/all-as-svg
+     svg/serialize
+     (spit "out/uniform-distribution.svg"))

src/thi/ng/geom/circle.cljc

@@ -24,6 +24,15 @@
   (let [m (m/mix p q)]
     (isec/intersect-circle-circle? _ (circle m (g/dist m p)))))
 
+;; https://stats.stackexchange.com/questions/481543/generating-random-points-uniformly-on-a-disk
+(defn random-point-in-circle
+  "Randomly generate a point inside of circle with uniform distribution."
+  [{p :p r :r}]
+  (-> (vec2 (* r (Math/sqrt (m/random)))
+            (* m/TWO_PI (m/random)))
+      g/as-cartesian
+      (m/+ p)))
+
 ;; Even though a circle is a specialization of an ellipse, we define
 ;; an extra type for performance reasons.
 

src/thi/ng/geom/polygon.cljc

@@ -3,6 +3,7 @@
    [thi.ng.geom.core :as g]
    [thi.ng.geom.utils :as gu]
    [thi.ng.geom.utils.intersect :as isec]
+   [thi.ng.geom.utils.probability :as prob]
    [thi.ng.geom.vector :as v :refer [vec2 vec3]]
    [thi.ng.geom.line :as l]
    [thi.ng.geom.triangle :as t]
@@ -364,7 +365,19 @@
     [{points :points} t] (gu/point-at t (conj points (first points))))
   (random-point
     [_] (g/point-at _ (m/random)))
-  (random-point-inside [_] nil) ; TODO
+  ;; Uniformly sample point from tessellated triangles of polygon
+  ;; https://blogs.sas.com/content/iml/2020/10/21/random-points-in-polygon.html
+  ;; https://observablehq.com/@scarysize/finding-random-points-in-a-polygon
+  (random-point-inside
+    [_]
+    "Uniformly sample point from inside polygon
+
+    Tessellates into triangles, randomly selects a triangle weighted by area,
+    and then samples a point inside of that triangle uniformly."
+    (->> (g/tessellate _)
+         (map t/triangle2)
+         (prob/rand-weighted-by g/area)
+         t/random-point-in-triangle2))
   (sample-uniform
     [{points :points} udist include-last?]
     (gu/sample-uniform udist include-last? (conj points (first points))))

src/thi/ng/geom/triangle.cljc

@@ -164,6 +164,17 @@
    (let [[o r] (circumcircle-raw a b c)]
      (Circle2. o r))))
 
+;; Kraemer Method
+;; http://extremelearning.com.au/evenly-distributing-points-in-a-triangle/
+;; https://stackoverflow.com/questions/47410054/generate-random-locations-within-a-triangular-domain/47418580#47418580
+(defn random-point-in-triangle2
+  "Randomly generate a point inside a 2d triangle with uniform distribution."
+  [_]
+  (let [[a b] [(m/random) (m/random)]
+        [s t] (if (< a b) [a b] [b a])
+        weighting [s (- t s) (- 1 t)]]
+    (gu/from-barycentric (get _ :points) weighting)))
+
 (extend-type Triangle2
 
   g/IArea

src/thi/ng/geom/utils/probability.cljc

@@ -0,0 +1,35 @@
+(ns thi.ng.geom.utils.probability
+  (:require [thi.ng.math.core :as m]))
+
+;; TODO: move namespace to thi.ng.math.core?
+
+(defn- mapping
+  "Create a mapping of every value in collection to f(value)."
+  [f coll]
+  (reduce (fn [m x] (assoc m x (f x))) {} coll))
+
+(defn rand-weighted
+  "Given a mapping of values to weights, randomly choose a value biased by weight"
+  [weights]
+  (let [sample (m/random (apply + (vals weights)))]
+    (loop [cumulative 0.0
+           [[choice weight] & remaining] weights]
+      (when weight
+        (let [sum (+ cumulative weight)]
+          (if (< sample sum)
+            choice
+            (recur sum remaining)))))))
+
+(defn rand-weighted-by
+  "Given a sequence of values `xs`, weight each value by a function `f` and return
+  a weighted random selection."
+  [f xs]
+  (rand-weighted (mapping f xs)))
+
+(comment
+  (rand-weighted {})
+  (frequencies (repeatedly 1000 #(rand-weighted {:a 0.2 :b 0.8})))
+  (frequencies (repeatedly 1000 #(rand-weighted {:a 0.1 :b 0.7 :c 0.2})))
+  (frequencies (repeatedly 1000 #(rand-weighted {:a 2 :b 8 :c 32})))
+  (rand-weighted-by inc [])
+  (frequencies (repeatedly 1000 #(rand-weighted-by inc [0 2 5]))))