Non-convex region: replace star polygon with a two-bulb dumbbell - concept-collection/hitandrun-interactive

concept-collection / hitandrun-interactive

Non-convex region: replace star polygon with a two-bulb dumbbell

Revamp the non-convex region to a dumbbell — two convex disks joined by a narrow tube. Unlike the star polygon it is genuinely non-star-shaped (the tube occludes each bulb from the other), which is the point of the local vs. union modes: the local walk mostly stays in one bulb while the union walk hops between them. Built as one CCW loop from each disk's major arc, with the tube edges closing the gaps (verified simple + non-star over 20k random instances). Since a dumbbell isn't star-shaped and need not contain the origin, hit_and_run_general now seeds itself with an interior point found by rejection instead of starting at the origin.

Jeremy Magland <jmagland@flatironinstitute.org> committed 1 month ago commit 76fb400a4633 parent a50c671 Browse files

4 changed files+77−36

CLAUDE.mdmodified+6−3View file

@@ -49,8 +49,11 @@ receive via `HTMLEventReceivedFcn`

4949 **Resample** button, and a **New region** button.

5050 3. **(done)** Wire those controls to the sampler (two-way bridge).

5151 4. **(done)** Non-convex regions. A **non-convex region** checkbox switches to a

52- star polygon (`make_region(false)`) sampled by `hit_and_run_general`; the

53- movie draws the multiple in-region segments a line makes through it.

52+ dumbbell — two convex bulbs joined by a narrow tube (`make_region(false)`),

53+ sampled by `hit_and_run_general`; the movie draws the multiple in-region

54+ segments a line makes through it. The dumbbell is non-star-shaped and need not

55+ contain the origin, so the general sampler seeds itself with an interior point

56+ found by rejection.

5457 5. **(done)** Local vs. union sampling. A **local segment only** checkbox (shown

5558 only when non-convex) switches `hit_and_run_general` between sampling the

5659 union of all in-region segments (default; uniform) and only the segment

@@ -95,7 +98,7 @@ walk.

9598 | `app/src/render/RegionView.tsx` | 2D canvas renderer (region outline + sample dots + movie `segments`) |

9699 | `app/src/bridge.ts` | `onData` / `onHostEvent` / `sendToMATLAB` helpers (generic; identical across these projects) |

97100 | `hitandrun_sampler.m` | Opens the figure, sends data, handles `resample`/`newRegion` (figure plumbing: `on_event`, `sample_region`, `pack_data`, `pack_samples`) |

98-| `helpers/make_region.m` | Random region — convex polygon on an ellipse (no convhull, so the auto-running figure view doesn't race the qhull backend load), or `make_region(false)` for a star polygon |

101+| `helpers/make_region.m` | Random region — convex polygon on an ellipse (no convhull, so the auto-running figure view doesn't race the qhull backend load), or `make_region(false)` for a dumbbell (two disks + narrow tube) |

99102 | `helpers/hit_and_run.m` | Convex hit-and-run sampler (half-plane chord intersection; JIT) |

100103 | `helpers/hit_and_run_general.m` | Non-convex hit-and-run sampler (union of all in-region segments, or `local` = only the one through the current point; interpreter) |

101104 | `hitandrun_demo.m` | User-facing driver: `addpath('helpers')` + seed + call the sampler |

README.mdmodified+14−13View file

@@ -14,22 +14,22 @@ A random 2D region is generated and `N` points are drawn uniformly from it by

1414 **hit-and-run**: from the current point, pick a random direction, take the chord

1515 where that line crosses the region, and jump to a uniform point on it. Repeat.

1616 The figure shows the region and the samples. The region is convex by default,

17-but you can switch to a non-convex (star-shaped) one, where a single line can

18-enter and leave the region several times.

17+but you can switch to a non-convex **dumbbell** — two convex bulbs joined by a

18+narrow tube — where a single line can enter and leave the region more than once.

1919

2020 Controls:

2121

2222 - **Samples** — set `N` (re-runs the sampler).

2323 - **Resample** — new samples, same region.

2424 - **New region** — a fresh region (of the current type).

25-- **non-convex region** — toggle between a convex region and a non-convex

26- star-shaped one. The non-convex region uses the slower general sampler, so

27- `N` is capped lower.

25+- **non-convex region** — toggle between a convex region and the non-convex

26+ dumbbell. The non-convex region uses the slower general sampler, so `N` is

27+ capped lower.

2828 - **local segment only** *(non-convex only)* — by default the step samples

2929 across *every* segment where the line crosses the region (the standard walk,

3030 uniform on the whole region); check this to instead sample only within the

31- segment that contains the current point, a local walk that can't hop across a

32- concavity and so does not sample uniformly.

31+ segment that contains the current point — a local walk that can't jump across

32+ the gap between two bulbs, so it does not sample uniformly.

3333 - **Play movie** — step through the algorithm: each step draws the chord — one

3434 segment for a convex region, possibly several for a non-convex one (or just

3535 the local segment) — and the point that landed on it.

@@ -38,7 +38,7 @@ Controls:

3838

3939 - [`hitandrun_demo.m`](hitandrun_demo.m) — driver: `addpath('helpers')`, seed, call the sampler.

4040 - [`hitandrun_sampler.m`](hitandrun_sampler.m) — opens the figure, sends data, handles resample requests, and dispatches to the convex or non-convex sampler.

41-- `helpers/` — [`make_region.m`](helpers/make_region.m) (convex or star region), [`hit_and_run.m`](helpers/hit_and_run.m) (convex chord), and [`hit_and_run_general.m`](helpers/hit_and_run_general.m) (arbitrary simple polygon).

41+- `helpers/` — [`make_region.m`](helpers/make_region.m) (convex ellipse or dumbbell region), [`hit_and_run.m`](helpers/hit_and_run.m) (convex chord), and [`hit_and_run_general.m`](helpers/hit_and_run_general.m) (arbitrary simple polygon).

4242 - `app/` — a single-file React app that draws the region and samples on a canvas.

4343

4444 The script and figure talk both ways: the script sends the region + samples via

@@ -65,11 +65,12 @@ stateless — the figure owns the region and passes it back with each request.

6565 correctly. It then samples in one of two modes: **union** (default) picks a

6666 point uniformly across all those segments, which samples the whole region

6767 uniformly; **local** (the *local segment only* checkbox) restricts to the

68- single segment straddling the current point — a local walk that can't cross a

69- concavity, so it does *not* sample uniformly. The sort + polygon tests don't

70- JIT, so this runs in the interpreter and `N` is capped lower.

71- `make_region(false)` builds a star polygon that is star-shaped about the

72- origin, so the origin is a valid interior start point.

68+ single segment straddling the current point — a local walk that can't jump the

69+ gap between the two bulbs, so it does *not* sample uniformly. The sort +

70+ polygon tests don't JIT, so this runs in the interpreter and `N` is capped

71+ lower. `make_region(false)` builds the dumbbell, which is not star-shaped and

72+ need not contain the origin, so the sampler finds an interior start by

73+ rejection rather than assuming one.

7374

7475 ## Deploy

7576

helpers/hit_and_run_general.mmodified+21−4View file

@@ -2,8 +2,8 @@ function [sx, sy] = hit_and_run_general(vx, vy, N, nBurn, local)

22 %HIT_AND_RUN_GENERAL Hit-and-run for an arbitrary simple polygon (convex or

33 % non-convex). Along each random line it finds *every* crossing of the polygon

44 % 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.

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.

77 %

88 % LOCAL (default false) selects how the step samples along that line:

99 % false — sample uniformly across the *union* of all in-region segments the

@@ -18,8 +18,7 @@ if nargin < 5 || isempty(local)

1818 local = false;

1919 end

2020 nv = numel(vx);

21-px = 0;

22-py = 0;

21+[px, py] = interior_seed(vx, vy);

2322 total = nBurn + N;

2423 sx = zeros(N, 1);

2524 sy = zeros(N, 1);

@@ -106,6 +105,24 @@ for step = 1:total

106105 end

107106 end

108107

108+function [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.

112+minx = min(vx); maxx = max(vx);

113+miny = min(vy); maxy = max(vy);

114+px = mean(vx); py = mean(vy);

115+for 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

123+end

124+end

125+

109126 function inside = point_in_poly(x, y, vx, vy)

110127 %POINT_IN_POLY Ray-casting test for a point against polygon (VX, VY).

111128 n = numel(vx);

helpers/make_region.mmodified+36−16View file

@@ -1,8 +1,9 @@

11 function [vx, vy] = make_region(convex)

22 %MAKE_REGION Random 2D sampling region, returned as CCW vertices (VX, VY).

33 % MAKE_REGION(true) — a convex polygon (vertices on a random ellipse).

4-% MAKE_REGION(false) — a non-convex star polygon (star-shaped about the

5-% origin, so the origin is a valid interior start point).

4+% MAKE_REGION(false) — a non-convex "dumbbell": two convex disks joined by a

5+% narrow tube. It is not star-shaped and need not contain

6+% the origin — the sampler finds its own interior start.

67 if nargin < 1 || isempty(convex)

78 convex = true;

89 end

@@ -10,7 +11,7 @@ end

1011 if convex

1112 [vx, vy] = convex_region();

1213 else

13- [vx, vy] = star_region();

14+ [vx, vy] = nonconvex_region();

1415 end

1516

1617 % CCW so the interior is to the left of each edge.

@@ -37,19 +38,38 @@ vx = cos(phi) * ex - sin(phi) * ey;

3738 vy = sin(phi) * ex + cos(phi) * ey;

3839 end

3940

40-function [vx, vy] = star_region()

41-% Spikes at sorted angles, alternating outer/inner radius. Small angle jitter

42-% keeps the angles ordered, so the polygon stays simple and star-shaped about

43-% the origin (lines through interior points can still leave and re-enter it).

44-spikes = 4 + randi(4); % 5..8 spikes

45-k = 2 * spikes;

46-step = 2 * pi / k;

47-ang = (0:k - 1) * step + (rand(1, k) - 0.5) * step * 0.6;

48-r = zeros(1, k);

49-r(1:2:k) = 0.9 + 0.4 * rand(1, numel(1:2:k)); % outer

50-r(2:2:k) = 0.35 + 0.2 * rand(1, numel(2:2:k)); % inner

51-vx = (1.3 * r .* cos(ang)).';

52-vy = (1.0 * r .* sin(ang)).';

41+function [vx, vy] = nonconvex_region()

42+% A "dumbbell": two convex disks (radius r, centers at x = +/-cx) joined by a

43+% narrow tube of half-width w < r. This is non-convex and non-star-shaped — the

44+% tube occludes each bulb from the other, so no single point sees the whole

45+% region — which is what separates the two sampling modes (the local walk mostly

46+% stays in one bulb, escaping only along a line down the tube; the union walk

47+% also hops between bulbs along a line that clips both without the tube). Built

48+% as one CCW loop: the outer (major) arc of each disk, with the straight tube

49+% edges closing the gaps between the arc ends.

50+r = 0.52 + 0.12 * rand; % bulb radius

51+cx = 0.82 + 0.24 * rand; % half-distance between bulb centers (> r: disjoint)

52+w = 0.07 + 0.06 * rand; % tube half-width (a narrow neck)

53+w = min(w, 0.5 * r); % keep the tube clearly narrower than a bulb

54+gam = asin(w / r); % half-angle each bulb's tube opening subtends

55+K = 16; % points per bulb arc

56+

57+% Right bulb: major arc from the bottom opening CCW round to the top opening.

58+tR = linspace(pi + gam, 3 * pi - gam, K);

59+rxx = cx + r * cos(tR);

60+ryy = r * sin(tR);

61+% Left bulb: major arc from the top opening CCW round to the bottom opening.

62+tL = linspace(gam, 2 * pi - gam, K);

63+lxx = -cx + r * cos(tL);

64+lyy = r * sin(tL);

65+

66+% Concatenate; the jumps arc-end -> next-arc-start are the straight tube edges.

67+px = [rxx, lxx].';

68+py = [ryy, lyy].';

69+

70+phi = 2 * pi * rand; % random overall orientation

71+vx = cos(phi) * px - sin(phi) * py;

72+vy = sin(phi) * px + cos(phi) * py;

5373 end

5474

5575 function A = signed_area(vx, vy)