diff --git a/common/mathlib.cc b/common/mathlib.cc
index 0a7a3de9..9072b193 100644
--- a/common/mathlib.cc
+++ b/common/mathlib.cc
@@ -465,6 +465,42 @@ glm::vec3 GLM_PolyCentroid(const std::vector<glm::vec3> &points)
     return poly_centroid;
 }
 
+glm::vec3 GLM_PolyRandomPoint(const std::vector<glm::vec3> &points)
+{
+    Q_assert(points.size() >= 3);
+    
+    // FIXME: Precompute this
+    float poly_area = 0;
+    std::vector<float> triareas;
+    
+    const vec3 v0 = points.at(0);
+    for (int i = 2; i < points.size(); i++) {
+        const vec3 v1 = points.at(i-1);
+        const vec3 v2 = points.at(i);
+        
+        const float triarea = GLM_TriangleArea(v0, v1, v2);
+        Q_assert(triarea >= 0.0f);
+        
+        triareas.push_back(triarea);
+        poly_area += triarea;
+    }
+    
+    // Pick a random triangle, with probability proportional to triangle area
+    const float uniformRandom = Random();
+    const std::vector<float> cdf = MakeCDF(triareas);
+    const int whichTri = SampleCDF(cdf, uniformRandom);
+    
+    Q_assert(whichTri >= 0 && whichTri < triareas.size());
+
+    const tri_t tri { points.at(0), points.at(1 + whichTri), points.at(2 + whichTri) };
+    
+    // Pick random barycentric coords.
+    const glm::vec3 bary = Barycentric_Random(Random(), Random());
+    const glm::vec3 point = Barycentric_ToPoint(bary, tri);
+    
+    return point;
+}
+
 std::pair<int, glm::vec3> GLM_ClosestPointOnPolyBoundary(const std::vector<glm::vec3> &poly, const vec3 &point)
 {
     const int N = static_cast<int>(poly.size());
diff --git a/include/common/mathlib.hh b/include/common/mathlib.hh
index e5a83dca..0aedc0a4 100644
--- a/include/common/mathlib.hh
+++ b/include/common/mathlib.hh
@@ -343,6 +343,7 @@ std::pair<bool, glm::vec3> GLM_InterpolateNormal(const std::vector<glm::vec3> &p
 std::vector<glm::vec3> GLM_ShrinkPoly(const std::vector<glm::vec3> &poly, const float amount);
 /// Returns (front part, back part)
 std::pair<std::vector<glm::vec3>,std::vector<glm::vec3>> GLM_ClipPoly(const std::vector<glm::vec3> &poly, const glm::vec4 &plane);
+glm::vec3 GLM_PolyRandomPoint(const std::vector<glm::vec3> &points);
 
 // Returns weights for f(0,0), f(1,0), f(0,1), f(1,1)
 // from: https://en.wikipedia.org/wiki/Bilinear_interpolation#Unit_Square
diff --git a/light/test_light.cc b/light/test_light.cc
index 1b74beaa..cc7411ea 100644
--- a/light/test_light.cc
+++ b/light/test_light.cc
@@ -459,3 +459,42 @@ TEST(mathlib, pointsAlongLine) {
     ASSERT_TRUE(pointsEqualEpsilon(vec3(1,0,0), res[0], POINT_EQUAL_EPSILON));
     ASSERT_TRUE(pointsEqualEpsilon(vec3(2.5,0,0), res[1], POINT_EQUAL_EPSILON));
 }
+
+TEST(mathlib, RandomPointInPoly) {
+    const vector<vec3> poly {
+        { 0,0,0 },
+        { 0,32,0 }, // colinear point
+        { 0,64,0 },
+        { 64,64,0 },
+        { 64,0,0 }
+    };
+    
+    const auto edgeplanes = GLM_MakeInwardFacingEdgePlanes(poly);
+    
+    glm::vec3 min(FLT_MAX);
+    glm::vec3 max(-FLT_MAX);
+    glm::vec3 avg(0);
+    
+    const int N=100;
+    for (int i=0; i<N; i++) {
+        const glm::vec3 point = GLM_PolyRandomPoint(poly);
+        ASSERT_TRUE(GLM_EdgePlanes_PointInside(edgeplanes, point));
+        
+        //std::cout << "point: " << glm::to_string(point) << std::endl;
+        
+        min = glm::min(min, point);
+        max = glm::max(max, point);
+        avg += point;
+    }
+    avg /= N;
+    
+    ASSERT_LT(min.x, 4);
+    ASSERT_LT(min.y, 4);
+    ASSERT_EQ(min.z, 0);
+    
+    ASSERT_GT(max.x, 60);
+    ASSERT_GT(max.y, 60);
+    ASSERT_EQ(max.z, 0);
+    
+    ASSERT_LT(glm::length(avg - glm::vec3(32, 32, 0)), 4);
+}