d9d003fAdd non-convex (star) regions with a general hit-and-run samplerJeremy Magland 1function [sx, sy] = hit_and_run_general(vx, vy, N, nBurn)
2%HIT_AND_RUN_GENERAL Hit-and-run for an arbitrary simple polygon (convex or
3% non-convex). Along each random line it finds *every* in-region segment and
4% samples uniformly across their union, so concavities are handled correctly
5% (unlike the convex-only chord in hit_and_run.m). Starts at the origin —
6% make_region's non-convex regions are star-shaped about it.
7%
8% This runs in the numbl interpreter (sort + point-in-polygon tests don't
9% JIT), so it's used only for the non-convex demo at modest N.
10nv = numel(vx);
11px = 0;
12py = 0;
13total = nBurn + N;
14sx = zeros(N, 1);
15sy = zeros(N, 1);
17for step = 1:total
18 th = 2 * pi * rand;
19 dx = cos(th);
20 dy = sin(th);
22 % All parameters t where the line p + t*d crosses the polygon boundary.
23 ts = zeros(1, nv);
24 m = 0;
25 for i = 1:nv
26 j = mod(i, nv) + 1;
27 ex = vx(j) - vx(i);
28 ey = vy(j) - vy(i);
29 denom = dy * ex - dx * ey;
30 if abs(denom) < 1e-12
31 continue
32 end
33 wx = vx(i) - px;
34 wy = vy(i) - py;
35 sParam = (dx * wy - dy * wx) / denom; % position along the edge
36 if sParam >= 0 && sParam < 1
37 m = m + 1;
38 ts(m) = (wy * ex - wx * ey) / denom; % position along the line
39 end
40 end
41 if m < 2
42 continue
43 end
44 ts = sort(ts(1:m));
46 % In-region intervals are consecutive crossings whose midpoint is inside.
47 totalLen = 0;
48 for k = 1:m - 1
49 tm = (ts(k) + ts(k + 1)) / 2;
50 if point_in_poly(px + tm * dx, py + tm * dy, vx, vy)
51 totalLen = totalLen + (ts(k + 1) - ts(k));
52 end
53 end
54 if totalLen <= 0
55 continue
56 end
58 % Pick a point uniformly across the union of in-region intervals.
59 u = totalLen * rand;
60 tpick = 0;
61 for k = 1:m - 1
62 tm = (ts(k) + ts(k + 1)) / 2;
63 if point_in_poly(px + tm * dx, py + tm * dy, vx, vy)
64 len = ts(k + 1) - ts(k);
65 if u <= len
66 tpick = ts(k) + u;
67 break
68 end
69 u = u - len;
70 end
71 end
72 px = px + tpick * dx;
73 py = py + tpick * dy;
75 if step > nBurn
76 sx(step - nBurn) = px;
77 sy(step - nBurn) = py;
78 end
79end
80end
82function inside = point_in_poly(x, y, vx, vy)
83%POINT_IN_POLY Ray-casting test for a point against polygon (VX, VY).
84n = numel(vx);
85inside = false;
86j = n;
87for i = 1:n
88 if ((vy(i) > y) ~= (vy(j) > y)) && ...
89 (x < (vx(j) - vx(i)) * (y - vy(i)) / (vy(j) - vy(i)) + vx(i))
90 inside = ~inside;
91 end
92 j = i;
93end
94end