9fbdc3eInteractive hit-and-run sampling of a 2D convex regionJeremy Magland 1function [sx, sy] = hit_and_run(vx, vy, N, nBurn)
2%HIT_AND_RUN Draw N uniform samples from a convex polygon by hit-and-run.
3% [SX, SY] = HIT_AND_RUN(VX, VY, N, NBURN) draws N samples (after NBURN
4% burn-in steps) from the uniform distribution on the convex polygon with
5% counterclockwise vertices (VX, VY).
6%
7% Represents the polygon as a set of half-planes (inside iff n_i . p >= c_i
8% for every edge i). From the current interior point, pick a random
9% direction, intersect the line with every half-plane to get the chord
10% [tmin, tmax], then jump to a uniform-random point on it.
11%
13% Guard: this hot loop must JS-JIT-compile (it's ~30x slower in the
14% interpreter).
15%!numbl:assert_jit
17nv = numel(vx);
19% Inward half-plane form for each edge: the left normal of a CCW edge points
20% into the interior.
21nx = zeros(nv, 1);
22ny = zeros(nv, 1);
23c = zeros(nv, 1);
24for i = 1:nv
25 j = mod(i, nv) + 1;
26 ex = vx(j) - vx(i);
27 ey = vy(j) - vy(i);
28 len = hypot(ex, ey);
29 n1 = -ey / len;
30 n2 = ex / len;
31 nx(i) = n1;
32 ny(i) = n2;
33 c(i) = n1 * vx(i) + n2 * vy(i);
34end
36% Start at the centroid (always interior for a convex polygon).
37px = mean(vx);
38py = mean(vy);
40total = nBurn + N;
41sx = zeros(N, 1);
42sy = zeros(N, 1);
43for s = 1:total
44 th = 2 * pi * rand;
45 dx = cos(th);
46 dy = sin(th);
47 % Chord [tmin, tmax] of the line p + t*d that stays inside the region.
48 tmin = -inf;
49 tmax = inf;
50 for i = 1:nv
51 a = nx(i) * dx + ny(i) * dy;
52 rhs = c(i) - (nx(i) * px + ny(i) * py); % <= 0 since p is interior
53 if a > 1e-12
54 tmin = max(tmin, rhs / a);
55 elseif a < -1e-12
56 tmax = min(tmax, rhs / a);
57 end
58 end
59 t = tmin + (tmax - tmin) * rand;
60 px = px + t * dx;
61 py = py + t * dy;
62 if s > nBurn
63 sx(s - nBurn) = px;
64 sy(s - nBurn) = py;
65 end
66end
67end