diff --git a/examples/svg/uniform_distribution.clj b/examples/svg/uniform_distribution.clj new file mode 100644 index 00000000..3fafda1f --- /dev/null +++ b/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")) diff --git a/src/thi/ng/geom/circle.cljc b/src/thi/ng/geom/circle.cljc index 45a21995..c5f1465c 100644 --- a/src/thi/ng/geom/circle.cljc +++ b/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. diff --git a/src/thi/ng/geom/polygon.cljc b/src/thi/ng/geom/polygon.cljc index c6758a5d..9c829d16 100644 --- a/src/thi/ng/geom/polygon.cljc +++ b/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)))) diff --git a/src/thi/ng/geom/triangle.cljc b/src/thi/ng/geom/triangle.cljc index 5a0aff1e..c383d1b5 100644 --- a/src/thi/ng/geom/triangle.cljc +++ b/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 diff --git a/src/thi/ng/geom/utils/probability.cljc b/src/thi/ng/geom/utils/probability.cljc new file mode 100644 index 00000000..f5e3ca98 --- /dev/null +++ b/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]))))