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). The region need not be star-shaped or contain the origin,
6% so it starts from an interior point found by rejection.
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);
21[px, py] = interior_seed(vx, vy);
22total = nBurn + N;
23sx = zeros(N, 1);
24sy = zeros(N, 1);
26for step = 1:total
27 th = 2 * pi * rand;
28 dx = cos(th);
29 dy = sin(th);
31 % All parameters t where the line p + t*d crosses the polygon boundary.
32 ts = zeros(1, nv);
33 m = 0;
34 for i = 1:nv
35 j = mod(i, nv) + 1;
36 ex = vx(j) - vx(i);
37 ey = vy(j) - vy(i);
38 denom = dy * ex - dx * ey;
39 if abs(denom) < 1e-12
40 continue
41 end
42 wx = vx(i) - px;
43 wy = vy(i) - py;
44 sParam = (dx * wy - dy * wx) / denom; % position along the edge
45 if sParam >= 0 && sParam < 1
46 m = m + 1;
47 ts(m) = (wy * ex - wx * ey) / denom; % position along the line
48 end
49 end
50 if m < 2
51 continue
52 end
53 ts = sort(ts(1:m));
55 if local
56 % Local segment: the in-region interval straddling t = 0, i.e. bounded
57 % by the nearest crossing on each side of the current (interior) point.
58 tlo = -inf;
59 thi = inf;
60 for k = 1:m
61 if ts(k) <= 0 && ts(k) > tlo
62 tlo = ts(k);
63 elseif ts(k) >= 0 && ts(k) < thi
64 thi = ts(k);
65 end
66 end
67 if ~isfinite(tlo) || ~isfinite(thi) || thi <= tlo
68 continue
69 end
70 tpick = tlo + (thi - tlo) * rand;
71 else
72 % In-region intervals are consecutive crossings whose midpoint is
73 % inside; sample uniformly across their union.
74 totalLen = 0;
75 for k = 1:m - 1
76 tm = (ts(k) + ts(k + 1)) / 2;
77 if point_in_poly(px + tm * dx, py + tm * dy, vx, vy)
78 totalLen = totalLen + (ts(k + 1) - ts(k));
79 end
80 end
81 if totalLen <= 0
82 continue
83 end
84 u = totalLen * rand;
85 tpick = 0;
86 for k = 1:m - 1
87 tm = (ts(k) + ts(k + 1)) / 2;
88 if point_in_poly(px + tm * dx, py + tm * dy, vx, vy)
89 len = ts(k + 1) - ts(k);
90 if u <= len
91 tpick = ts(k) + u;
92 break
93 end
94 u = u - len;
95 end
96 end
97 end
98 px = px + tpick * dx;
99 py = py + tpick * dy;
101 if step > nBurn
102 sx(step - nBurn) = px;
103 sy(step - nBurn) = py;
104 end
105end
106end
108function [px, py] = interior_seed(vx, vy)
109%INTERIOR_SEED A point strictly inside polygon (VX, VY), by rejection sampling
110% its bounding box. The polygon fills a large fraction of the box, so this
111% lands quickly; the vertex mean is a fallback if it somehow doesn't.
112minx = min(vx); maxx = max(vx);
113miny = min(vy); maxy = max(vy);
114px = mean(vx); py = mean(vy);
115for t = 1:5000
116 qx = minx + (maxx - minx) * rand;
117 qy = miny + (maxy - miny) * rand;
118 if point_in_poly(qx, qy, vx, vy)
119 px = qx;
120 py = qy;
121 return
122 end
123end
124end
126function inside = point_in_poly(x, y, vx, vy)
127%POINT_IN_POLY Ray-casting test for a point against polygon (VX, VY).
128n = numel(vx);
129inside = false;
130j = n;
131for i = 1:n
132 if ((vy(i) > y) ~= (vy(j) > y)) && ...
133 (x < (vx(j) - vx(i)) * (y - vy(i)) / (vy(j) - vy(i)) + vx(i))
134 inside = ~inside;
135 end
136 j = i;
137end
138end