1function [sx, sy] = hit_and_run_general(vx, vy, N, nBurn, local)
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* crossing of the polygon
4% boundary, so concavities are handled correctly (unlike the convex-only chord
5% in hit_and_run.m). Starts at the origin — make_region's non-convex regions
6% are star-shaped about it.
7%
8% LOCAL (default false) selects how the step samples along that line:
9% false — sample uniformly across the *union* of all in-region segments the
10% line makes (standard hit-and-run; uniform on the whole region).
11% true — sample only within the single in-region segment that contains the
12% current point (a local walk that can't jump across a concavity;
13% it does not sample the region uniformly).
14%
15% This runs in the numbl interpreter (sort + point-in-polygon tests don't
16% JIT), so it's used only for the non-convex demo at modest N.
17if nargin < 5 || isempty(local)
18 local = false;
19end
20nv = numel(vx);
21px = 0;
22py = 0;
23total = nBurn + N;
24sx = zeros(N, 1);
25sy = zeros(N, 1);
27for step = 1:total
28 th = 2 * pi * rand;
29 dx = cos(th);
30 dy = sin(th);
32 % All parameters t where the line p + t*d crosses the polygon boundary.
33 ts = zeros(1, nv);
34 m = 0;
35 for i = 1:nv
36 j = mod(i, nv) + 1;
37 ex = vx(j) - vx(i);
38 ey = vy(j) - vy(i);
39 denom = dy * ex - dx * ey;
40 if abs(denom) < 1e-12
41 continue
42 end
43 wx = vx(i) - px;
44 wy = vy(i) - py;
45 sParam = (dx * wy - dy * wx) / denom; % position along the edge
46 if sParam >= 0 && sParam < 1
47 m = m + 1;
48 ts(m) = (wy * ex - wx * ey) / denom; % position along the line
49 end
50 end
51 if m < 2
52 continue
53 end
54 ts = sort(ts(1:m));
56 if local
57 % Local segment: the in-region interval straddling t = 0, i.e. bounded
58 % by the nearest crossing on each side of the current (interior) point.
59 tlo = -inf;
60 thi = inf;
61 for k = 1:m
62 if ts(k) <= 0 && ts(k) > tlo
63 tlo = ts(k);
64 elseif ts(k) >= 0 && ts(k) < thi
65 thi = ts(k);
66 end
67 end
68 if ~isfinite(tlo) || ~isfinite(thi) || thi <= tlo
69 continue
70 end
71 tpick = tlo + (thi - tlo) * rand;
72 else
73 % In-region intervals are consecutive crossings whose midpoint is
74 % inside; sample uniformly across their union.
75 totalLen = 0;
76 for k = 1:m - 1
77 tm = (ts(k) + ts(k + 1)) / 2;
78 if point_in_poly(px + tm * dx, py + tm * dy, vx, vy)
79 totalLen = totalLen + (ts(k + 1) - ts(k));
80 end
81 end
82 if totalLen <= 0
83 continue
84 end
85 u = totalLen * rand;
86 tpick = 0;
87 for k = 1:m - 1
88 tm = (ts(k) + ts(k + 1)) / 2;
89 if point_in_poly(px + tm * dx, py + tm * dy, vx, vy)
90 len = ts(k + 1) - ts(k);
91 if u <= len
92 tpick = ts(k) + u;
93 break
94 end
95 u = u - len;
96 end
97 end
98 end
99 px = px + tpick * dx;
100 py = py + tpick * dy;
102 if step > nBurn
103 sx(step - nBurn) = px;
104 sy(step - nBurn) = py;
105 end
106end
107end
109function inside = point_in_poly(x, y, vx, vy)
110%POINT_IN_POLY Ray-casting test for a point against polygon (VX, VY).
111n = numel(vx);
112inside = false;
113j = n;
114for i = 1:n
115 if ((vy(i) > y) ~= (vy(j) > y)) && ...
116 (x < (vx(j) - vx(i)) * (y - vy(i)) / (vy(j) - vy(i)) + vx(i))
117 inside = ~inside;
118 end
119 j = i;
120end
121end