Add non-convex (star) regions with a general hit-and-run sampler
A 'non-convex region' checkbox switches make_region to a star polygon that is
star-shaped about the origin. Convex regions keep the fast JIT hit_and_run
(single chord); non-convex regions use the new hit_and_run_general, which finds
every crossing along the line, keeps the in-region segments, and samples
uniformly across their union — so concavities are handled correctly. It runs in
the interpreter, so N is capped lower for that mode.
The convex flag is threaded through the figure<->script protocol and the movie
overlay now draws the multiple in-region segments a line makes through a
non-convex region. README and CLAUDE.md updated to document the feature.
7 changed files+360−137
CLAUDE.mdmodified+29−17View file
@@ -8,9 +8,12 @@ figure is a 2D canvas instead of a three.js surface.
88
99 - `app/` — React single-file widget (`npm run build` → `app/dist/index.html`)
1010 - `*.m` at root — scripts users open and run in [numbl](https://numbl.org)
11-- `helpers/` — the region/sampling algorithm (`make_region`, `hit_and_run`),
12- put on the path by the driver via `addpath('helpers')`
13-- `numbl-project.json` — project metadata; `entry` is the landing page
11+- `helpers/` — the region/sampling algorithm (`make_region`, `hit_and_run`,
12+ `hit_and_run_general`), put on the path by the driver via `addpath('helpers')`
13+- `numbl-project.json` — project metadata; `entry` is the landing page. Its
14+ `figures` array declares the editor-less figure view the deploy exposes at
15+ `#figure/sampler` (runs `hitandrun_demo.m`, shows just the interactive figure);
16+ the org profile's "live" link points there.
1417
1518 ### Multi-file layout / addpath
1619
@@ -45,20 +48,28 @@ receive via `HTMLEventReceivedFcn`
4548 2. **(done)** Add controls to the figure panel: a **Samples** slider, a
4649 **Resample** button, and a **New region** button.
4750 3. **(done)** Wire those controls to the sampler (two-way bridge).
51+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.
4854
49-The figure → script protocol:
55+The figure → script protocol (each request carries `convex`):
5056
5157 - **Samples** slider (on release) / **Resample** →
52- `sendToMATLAB('resample', { n, x, y })`. `x, y` are the *current region's*
53- vertices — the script is **stateless**, so it samples the region it's given
54- and replies with `sendEventToHTMLSource(src, 'samples', { x, y, n })`.
55-- **New region** → `sendToMATLAB('newRegion', { n })`. The script builds a new
56- region and replies with a full `sendEventToHTMLSource(src, 'data', <payload>)`
57- (same shape as the initial `Data`).
58+ `sendToMATLAB('resample', { n, x, y, convex })`. `x, y` are the *current
59+ region's* vertices — the script is **stateless**, so it samples the region it
60+ is given and replies with `sendEventToHTMLSource(src, 'samples', { x, y, n })`.
61+- **New region** / **non-convex** toggle → `sendToMATLAB('newRegion', { n,
62+ convex })`. The script builds a new region of that type and replies with a full
63+ `sendEventToHTMLSource(src, 'data', <payload>)` (same shape as the initial
64+ `Data`, plus `convex`).
65+
66+`hitandrun_sampler.m`'s `sample_region(vx, vy, n, convex)` dispatches to
67+`hit_and_run` (convex, JIT) or `hit_and_run_general` (non-convex, interpreter).
68+The non-convex path can't JIT (sort + point-in-polygon), so `App.tsx` caps its
69+`N` lower via `SAMPLE_CHOICES_NONCONVEX`.
5870
5971 Ideas for the next iteration: show the chain path / burn-in, animate the walk,
60-a "uniform vs. non-uniform target" toggle, or a non-convex region (hit-and-run
61-chords become multiple segments).
72+or a "uniform vs. non-uniform target" toggle.
6273
6374 ## Client-side vs. server-side interactivity
6475
@@ -71,12 +82,13 @@ chords become multiple segments).
7182
7283 | File | Purpose |
7384 |------|---------|
74-| `app/src/App.tsx` | Main React component; reads data, hosts the panel |
75-| `app/src/render/RegionView.tsx` | 2D canvas renderer (region outline + sample dots) |
85+| `app/src/App.tsx` | Main React component; reads data, hosts the panel, non-convex toggle + movie `regionSegments` |
86+| `app/src/render/RegionView.tsx` | 2D canvas renderer (region outline + sample dots + movie `segments`) |
7687 | `app/src/bridge.ts` | `onData` / `onHostEvent` / `sendToMATLAB` helpers (generic; identical across these projects) |
77-| `hitandrun_sampler.m` | Opens the figure, sends data, handles `resample`/`newRegion` (figure plumbing: `on_event`, `pack_data`, `pack_samples`) |
78-| `helpers/make_region.m` | Random convex region (convex hull of disk points) |
79-| `helpers/hit_and_run.m` | The hit-and-run sampler (half-plane chord intersection) |
88+| `hitandrun_sampler.m` | Opens the figure, sends data, handles `resample`/`newRegion` (figure plumbing: `on_event`, `sample_region`, `pack_data`, `pack_samples`) |
89+| `helpers/make_region.m` | Random region — convex hull of disk points, or `make_region(false)` for a star polygon |
90+| `helpers/hit_and_run.m` | Convex hit-and-run sampler (half-plane chord intersection; JIT) |
91+| `helpers/hit_and_run_general.m` | Non-convex hit-and-run sampler (all in-region segments; interpreter) |
8092 | `hitandrun_demo.m` | User-facing driver: `addpath('helpers')` + seed + call the sampler |
8193
8294 ## Local iteration
README.mdmodified+35−17View file
@@ -10,39 +10,57 @@ developer view — file tree, editable code, console — is below.
1010
1111 ## ▶ [Open `hitandrun_demo.m`](hitandrun_demo.m) and click **Run**
1212
13-A random 2D convex region is generated and `N` points are drawn uniformly from
14-it by **hit-and-run**: from the current point, pick a random direction, take the
15-chord where that line crosses the region, and jump to a uniform point on it.
16-Repeat. The figure shows the region and the samples.
13+A random 2D region is generated and `N` points are drawn uniformly from it by
14+**hit-and-run**: from the current point, pick a random direction, take the chord
15+where that line crosses the region, and jump to a uniform point on it. Repeat.
16+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.
1719
1820 Controls:
1921
2022 - **Samples** — set `N` (re-runs the sampler).
2123 - **Resample** — new samples, same region.
22-- **New region** — a fresh region.
23-- **Play movie** — step through the algorithm: each step draws the chord and the
24+- **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.
28+- **Play movie** — step through the algorithm: each step draws the chord — one
29+ segment for a convex region, possibly several for a non-convex one — and the
2430 point that landed on it.
2531
2632 ## How it works
2733
2834 - [`hitandrun_demo.m`](hitandrun_demo.m) — driver: `addpath('helpers')`, seed, call the sampler.
29-- [`hitandrun_sampler.m`](hitandrun_sampler.m) — opens the figure, sends data, handles resample requests.
30-- `helpers/` — [`make_region.m`](helpers/make_region.m) and [`hit_and_run.m`](helpers/hit_and_run.m).
35+- [`hitandrun_sampler.m`](hitandrun_sampler.m) — opens the figure, sends data, handles resample requests, and dispatches to the convex or non-convex sampler.
36+- `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).
3137 - `app/` — a single-file React app that draws the region and samples on a canvas.
3238
3339 The script and figure talk both ways: the script sends the region + samples via
3440 `uihtml(..., 'Data', ...)`, and the controls call back with
35-`sendToMATLAB('resample' | 'newRegion', ...)`, which re-runs the sampler and
36-returns new points via `sendEventToHTMLSource`. The script is stateless — the
37-figure owns the region and passes it back with each request.
41+`sendToMATLAB('resample' | 'newRegion', ...)` (each carrying whether the region
42+is `convex`), which re-runs the sampler and returns new points via
43+`sendEventToHTMLSource`. The script is stateless — the figure owns the region
44+and passes it back with each request.
3845
39-## JIT-compiled kernel
46+## Convex vs. non-convex
4047
41-The loop in [`hit_and_run.m`](helpers/hit_and_run.m) runs once per sample, so
42-numbl JS-JIT-compiles it to JavaScript — about 30× faster than its interpreter,
43-which is what keeps large `N` instant. The `%!numbl:assert_jit` directive
44-asserts this happens (it errors rather than silently falling back). It relies on
45-numbl's scalar-`rand()` JIT support; `rng(seed)` still controls the shared PRNG.
48+`hitandrun_sampler.m` picks the sampler by region type:
49+
50+- **Convex** — [`hit_and_run.m`](helpers/hit_and_run.m) takes the single chord
51+ where the line crosses the region (an inward-half-plane intersection). Its
52+ loop runs once per sample, so numbl JS-JIT-compiles it to JavaScript — about
53+ 30× faster than its interpreter, which is what keeps large `N` instant. The
54+ `%!numbl:assert_jit` directive asserts this happens (it errors rather than
55+ silently falling back). It relies on numbl's scalar-`rand()` JIT support;
56+ `rng(seed)` still controls the shared PRNG.
57+- **Non-convex** — [`hit_and_run_general.m`](helpers/hit_and_run_general.m)
58+ finds *every* crossing along the line, keeps the segments whose midpoint is
59+ inside (a point-in-polygon test), and samples uniformly across their union, so
60+ concavities are handled correctly. The sort + polygon tests don't JIT, so this
61+ runs in the interpreter and `N` is capped lower. `make_region(false)` builds a
62+ star polygon that is star-shaped about the origin, so the origin is a valid
63+ interior start point.
4664
4765 ## Deploy
4866
app/src/App.tsxmodified+97−40View file
@@ -15,6 +15,7 @@ interface HitAndRunData {
1515 region: Points;
1616 samples: Points;
1717 n: number;
18+ convex?: boolean; // false → non-convex (star) region
1819 }
1920
2021 function isHitAndRunData(d: unknown): d is HitAndRunData {
@@ -35,8 +36,10 @@ interface SamplesEvent {
3536 n: number;
3637 }
3738
38-// Discrete sample-count choices; the slider indexes into this array.
39+// Discrete sample-count choices; the slider indexes into the active array.
3940 const SAMPLE_CHOICES = [10, 100, 1000, 10000, 100000, 1000000];
41+// Non-convex sampling runs in the interpreter (no JIT), so cap it lower.
42+const SAMPLE_CHOICES_NONCONVEX = [10, 100, 1000, 10000];
4043 const DEFAULT_N = 10000;
4144
4245 // Sampling movie: reveal points one at a time. Each frame shows the chord the
@@ -48,39 +51,60 @@ interface MovieState {
4851 step: number; // index of the point currently being sampled
4952 }
5053
51-/** The chord that hit-and-run samples along: the line through (px,py) with
52- * direction (dx,dy), clipped to the convex region. Same inward-half-plane
53- * intersection the sampler uses, so it reproduces the step exactly. */
54-function regionChord(
54+/** The in-region segment(s) hit-and-run samples along: the line through
55+ * (px,py) with direction (dx,dy), intersected with the polygon. One segment
56+ * for a convex region, possibly several for a non-convex one. Mirrors the
57+ * sampler's geometry, so it reproduces each step exactly. */
58+function regionSegments(
5559 region: Points,
5660 px: number,
5761 py: number,
5862 dx: number,
5963 dy: number
60-): Segment | null {
61- if (Math.hypot(dx, dy) < 1e-12) return null;
62- let tmin = -Infinity;
63- let tmax = Infinity;
64+): Segment[] {
65+ if (Math.hypot(dx, dy) < 1e-12) return [];
6466 const m = region.x.length;
67+ const ts: number[] = [];
6568 for (let i = 0; i < m; i++) {
6669 const j = (i + 1) % m;
6770 const ex = region.x[j] - region.x[i];
6871 const ey = region.y[j] - region.y[i];
69- // Inward normal of a CCW edge. inside: n·p >= n·v_i -> t·a >= rhs.
70- const nx = -ey;
71- const ny = ex;
72- const a = nx * dx + ny * dy;
73- const rhs = nx * region.x[i] + ny * region.y[i] - (nx * px + ny * py);
74- if (a > 1e-12) tmin = Math.max(tmin, rhs / a);
75- else if (a < -1e-12) tmax = Math.min(tmax, rhs / a);
72+ const denom = dy * ex - dx * ey;
73+ if (Math.abs(denom) < 1e-12) continue;
74+ const wx = region.x[i] - px;
75+ const wy = region.y[i] - py;
76+ const s = (dx * wy - dy * wx) / denom; // position along the edge
77+ if (s >= 0 && s < 1) ts.push((wy * ex - wx * ey) / denom);
7678 }
77- if (!(tmax > tmin) || !isFinite(tmin) || !isFinite(tmax)) return null;
78- return {
79- x0: px + tmin * dx,
80- y0: py + tmin * dy,
81- x1: px + tmax * dx,
82- y1: py + tmax * dy,
83- };
79+ ts.sort((a, b) => a - b);
80+ const segs: Segment[] = [];
81+ for (let k = 0; k < ts.length - 1; k++) {
82+ const tm = (ts[k] + ts[k + 1]) / 2;
83+ if (pointInPolygon(region, px + tm * dx, py + tm * dy)) {
84+ segs.push({
85+ x0: px + ts[k] * dx,
86+ y0: py + ts[k] * dy,
87+ x1: px + ts[k + 1] * dx,
88+ y1: py + ts[k + 1] * dy,
89+ });
90+ }
91+ }
92+ return segs;
93+}
94+
95+function pointInPolygon(region: Points, x: number, y: number): boolean {
96+ const n = region.x.length;
97+ let inside = false;
98+ for (let i = 0, j = n - 1; i < n; j = i++) {
99+ const xi = region.x[i];
100+ const yi = region.y[i];
101+ const xj = region.x[j];
102+ const yj = region.y[j];
103+ if (yi > y !== yj > y && x < ((xj - xi) * (y - yi)) / (yj - yi) + xi) {
104+ inside = !inside;
105+ }
106+ }
107+ return inside;
84108 }
85109
86110 const prefixPoints = (p: Points, k: number): Points => ({
@@ -144,17 +168,33 @@ export function App() {
144168 return () => clearTimeout(id);
145169 }, [movie, data]);
146170
171+ const nonConvex = !!data && data.convex === false;
172+ const choices = nonConvex ? SAMPLE_CHOICES_NONCONVEX : SAMPLE_CHOICES;
173+
147174 // Re-draw `count` samples in the current region (script round-trip).
148175 const resample = (count: number) => {
149176 if (!data) return;
150177 setBusy(true);
151- sendToMATLAB("resample", { n: count, x: data.region.x, y: data.region.y });
178+ sendToMATLAB("resample", {
179+ n: count,
180+ x: data.region.x,
181+ y: data.region.y,
182+ convex: data.convex !== false,
183+ });
152184 };
153185
154- // Generate a brand new region and sample it (script round-trip).
155- const newRegion = (count: number) => {
186+ // Generate a brand new region (convex or not) and sample it.
187+ const newRegion = (count: number, convex: boolean) => {
156188 setBusy(true);
157- sendToMATLAB("newRegion", { n: count });
189+ sendToMATLAB("newRegion", { n: count, convex });
190+ };
191+
192+ // Checkbox: switch region type. Clamp N to the active set's max first.
193+ const setNonConvex = (makeNonConvex: boolean) => {
194+ const c = makeNonConvex ? SAMPLE_CHOICES_NONCONVEX : SAMPLE_CHOICES;
195+ const clamped = Math.min(n, c[c.length - 1]);
196+ setN(clamped);
197+ newRegion(clamped, !makeNonConvex);
158198 };
159199
160200 const toggleMovie = () => {
@@ -165,9 +205,9 @@ export function App() {
165205 }
166206 };
167207
168- // Overlay for the current movie frame (settled points + chord + highlights).
208+ // Overlay for the current movie frame (settled points + segments + marks).
169209 let cloud: Points = data ? data.samples : { x: [], y: [] };
170- let chord: Segment | null = null;
210+ let segments: Segment[] | null = null;
171211 let from: Pt | null = null;
172212 let newPoint: Pt | null = null;
173213 if (movie && data) {
@@ -176,7 +216,7 @@ export function App() {
176216 const { x, y } = data.samples;
177217 cloud = prefixPoints(data.samples, i); // settled points 0..i-1
178218 from = { x: x[f], y: y[f] };
179- chord = regionChord(data.region, x[f], y[f], x[i] - x[f], y[i] - y[f]);
219+ segments = regionSegments(data.region, x[f], y[f], x[i] - x[f], y[i] - y[f]);
180220 newPoint = { x: x[i], y: y[i] };
181221 }
182222
@@ -188,7 +228,7 @@ export function App() {
188228 <RegionView
189229 region={data.region}
190230 samples={cloud}
191- chord={chord}
231+ segments={segments}
192232 from={from}
193233 newPoint={newPoint}
194234 />
@@ -202,19 +242,17 @@ export function App() {
202242 <input
203243 type="range"
204244 min={0}
205- max={SAMPLE_CHOICES.length - 1}
245+ max={choices.length - 1}
206246 step={1}
207- value={Math.max(0, SAMPLE_CHOICES.indexOf(n))}
247+ value={Math.max(0, choices.indexOf(n))}
208248 disabled={!data || busy || !!movie}
209249 // Drag updates the label live; the script round-trip fires on
210250 // release (and on arrow-key release) to avoid flooding it.
211- onChange={e => setN(SAMPLE_CHOICES[Number(e.target.value)])}
212- onPointerUp={e =>
213- resample(SAMPLE_CHOICES[Number(e.currentTarget.value)])
214- }
251+ onChange={e => setN(choices[Number(e.target.value)])}
252+ onPointerUp={e => resample(choices[Number(e.currentTarget.value)])}
215253 onKeyUp={e => {
216254 if (e.key.startsWith("Arrow")) {
217- resample(SAMPLE_CHOICES[Number(e.currentTarget.value)]);
255+ resample(choices[Number(e.currentTarget.value)]);
218256 }
219257 }}
220258 style={sliderStyle}
@@ -233,13 +271,23 @@ export function App() {
233271 <button
234272 style={btnStyle}
235273 disabled={busy || !!movie}
236- onClick={() => newRegion(n)}
237- title="Generate a new convex region and sample it"
274+ onClick={() => newRegion(n, !nonConvex)}
275+ title="Generate a new region and sample it"
238276 >
239277 New region
240278 </button>
241279 </div>
242280
281+ <label style={checkLabelStyle}>
282+ <input
283+ type="checkbox"
284+ checked={nonConvex}
285+ disabled={!data || busy || !!movie}
286+ onChange={e => setNonConvex(e.target.checked)}
287+ />
288+ non-convex region
289+ </label>
290+
243291 <button
244292 style={playBtnStyle}
245293 disabled={!canPlay || busy}
@@ -296,6 +344,15 @@ const labelStyle: CSSProperties = {
296344 fontSize: 11,
297345 };
298346
347+const checkLabelStyle: CSSProperties = {
348+ display: "flex",
349+ alignItems: "center",
350+ gap: 5,
351+ fontSize: 11,
352+ marginTop: 8,
353+ cursor: "pointer",
354+};
355+
299356 const sliderStyle: CSSProperties = {
300357 width: "100%",
301358 marginTop: 2,
app/src/render/RegionView.tsxmodified+28−28View file
@@ -22,9 +22,9 @@ interface RegionViewProps {
2222 region: Points;
2323 /** Sample points to scatter inside the region. */
2424 samples: Points;
25- /** Animation overlay: the chord the current step samples along (a line
26- * through `from` clipped to the region). */
27- chord?: Segment | null;
25+ /** Animation overlay: the in-region segments the current step samples along
26+ * (one for a convex region; possibly several for a non-convex one). */
27+ segments?: Segment[] | null;
2828 /** Animation overlay: the current point the step starts from (ringed). */
2929 from?: Pt | null;
3030 /** Animation overlay: the freshly sampled point (highlighted). */
@@ -49,7 +49,7 @@ const NEW = "#0f172a"; // near-black: the freshly sampled point (circled)
4949 export function RegionView({
5050 region,
5151 samples,
52- chord,
52+ segments,
5353 from,
5454 newPoint,
5555 }: RegionViewProps) {
@@ -128,31 +128,31 @@ export function RegionView({
128128 ctx.fill();
129129 }
130130
131- // Animation overlay (movie mode).
132- if (chord) {
133- const x0 = toPx(chord.x0);
134- const y0 = toPy(chord.y0);
135- const x1 = toPx(chord.x1);
136- const y1 = toPy(chord.y1);
137- // The chord itself.
138- ctx.beginPath();
139- ctx.moveTo(x0, y0);
140- ctx.lineTo(x1, y1);
141- ctx.strokeStyle = CHORD;
142- ctx.lineWidth = 1.5;
143- ctx.setLineDash([4, 3]);
144- ctx.stroke();
145- ctx.setLineDash([]);
146- // Open circles where the line crosses the region boundary.
147- ctx.strokeStyle = CHORD;
148- ctx.lineWidth = 1.25;
149- for (const [ex, ey] of [
150- [x0, y0],
151- [x1, y1],
152- ]) {
131+ // Animation overlay (movie mode): the in-region segment(s) of the line.
132+ if (segments) {
133+ for (const seg of segments) {
134+ const x0 = toPx(seg.x0);
135+ const y0 = toPy(seg.y0);
136+ const x1 = toPx(seg.x1);
137+ const y1 = toPy(seg.y1);
153138 ctx.beginPath();
154- ctx.arc(ex, ey, 3, 0, 2 * Math.PI);
139+ ctx.moveTo(x0, y0);
140+ ctx.lineTo(x1, y1);
141+ ctx.strokeStyle = CHORD;
142+ ctx.lineWidth = 1.5;
143+ ctx.setLineDash([4, 3]);
155144 ctx.stroke();
145+ ctx.setLineDash([]);
146+ // Open circles where the line crosses the region boundary.
147+ ctx.lineWidth = 1.25;
148+ for (const [ex, ey] of [
149+ [x0, y0],
150+ [x1, y1],
151+ ]) {
152+ ctx.beginPath();
153+ ctx.arc(ex, ey, 3, 0, 2 * Math.PI);
154+ ctx.stroke();
155+ }
156156 }
157157 }
158158 // Previous sample: the point the chord starts from — a filled amber dot.
@@ -182,7 +182,7 @@ export function RegionView({
182182 const ro = new ResizeObserver(draw);
183183 ro.observe(container);
184184 return () => ro.disconnect();
185- }, [region, samples, chord, from, newPoint]);
185+ }, [region, samples, segments, from, newPoint]);
186186
187187 return (
188188 <div ref={containerRef} style={containerStyle}>
helpers/hit_and_run_general.madded+94−0View file
@@ -0,0 +1,94 @@
1+function [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.
10+nv = numel(vx);
11+px = 0;
12+py = 0;
13+total = nBurn + N;
14+sx = zeros(N, 1);
15+sy = zeros(N, 1);
16+
17+for step = 1:total
18+ th = 2 * pi * rand;
19+ dx = cos(th);
20+ dy = sin(th);
21+
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));
45+
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
57+
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;
74+
75+ if step > nBurn
76+ sx(step - nBurn) = px;
77+ sy(step - nBurn) = py;
78+ end
79+end
80+end
81+
82+function inside = point_in_poly(x, y, vx, vy)
83+%POINT_IN_POLY Ray-casting test for a point against polygon (VX, VY).
84+n = numel(vx);
85+inside = false;
86+j = n;
87+for 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;
93+end
94+end
helpers/make_region.mmodified+41−13View file
@@ -1,23 +1,51 @@
1-function [vx, vy] = make_region()
2-%MAKE_REGION Random convex polygon = convex hull of random points in a disk.
3-% [VX, VY] = MAKE_REGION() returns the counterclockwise vertices of a random
4-% convex region. Anisotropic scaling makes it a bit more interesting than a
5-% circle. The vertices are ordered CCW so the interior is to the left of each
6-% edge (which hit_and_run relies on for its inward half-plane normals).
1+function [vx, vy] = make_region(convex)
2+%MAKE_REGION Random 2D sampling region, returned as CCW vertices (VX, VY).
3+% MAKE_REGION(true) — a convex polygon (convex hull of random disk points).
4+% MAKE_REGION(false) — a non-convex star polygon (star-shaped about the
5+% origin, so the origin is a valid interior start point).
6+if nargin < 1 || isempty(convex)
7+ convex = true;
8+end
9+
10+if convex
11+ [vx, vy] = convex_region();
12+else
13+ [vx, vy] = star_region();
14+end
15+
16+% CCW so the interior is to the left of each edge.
17+if signed_area(vx, vy) < 0
18+ vx = vx(end:-1:1);
19+ vy = vy(end:-1:1);
20+end
21+end
22+
23+function [vx, vy] = convex_region()
24+% Convex hull of random points in a disk; anisotropic scaling for variety.
725 m = 12;
826 ang = 2 * pi * rand(m, 1);
9-r = sqrt(rand(m, 1)); % sqrt -> uniform over the disk
27+r = sqrt(rand(m, 1));
1028 px = 1.4 * r .* cos(ang);
1129 py = 1.0 * r .* sin(ang);
12-k = convhull(px, py); % boundary indices, closed (last == first)
13-k = k(1:end-1); % drop the repeated closing vertex
30+k = convhull(px, py);
31+k = k(1:end-1);
1432 vx = px(k);
1533 vy = py(k);
16-% Ensure counterclockwise.
17-if signed_area(vx, vy) < 0
18- vx = vx(end:-1:1);
19- vy = vy(end:-1:1);
2034 end
35+
36+function [vx, vy] = star_region()
37+% Spikes at sorted angles, alternating outer/inner radius. Small angle jitter
38+% keeps the angles ordered, so the polygon stays simple and star-shaped about
39+% the origin (lines through interior points can still leave and re-enter it).
40+spikes = 4 + randi(4); % 5..8 spikes
41+k = 2 * spikes;
42+step = 2 * pi / k;
43+ang = (0:k - 1) * step + (rand(1, k) - 0.5) * step * 0.6;
44+r = zeros(1, k);
45+r(1:2:k) = 0.9 + 0.4 * rand(1, numel(1:2:k)); % outer
46+r(2:2:k) = 0.35 + 0.2 * rand(1, numel(2:2:k)); % inner
47+vx = (1.3 * r .* cos(ang)).';
48+vy = (1.0 * r .* sin(ang)).';
2149 end
2250
2351 function A = signed_area(vx, vy)
hitandrun_sampler.mmodified+36−22View file
@@ -19,14 +19,13 @@ function hitandrun_sampler(N)
1919 % server-side state to keep in sync across callbacks.
2020
2121 if nargin < 1 || isempty(N)
22- N = 1500;
22+ N = 10000;
2323 end
2424
25-% Build a random convex region and draw N uniform samples from it.
26-[vx, vy] = make_region();
27-tic;
28-[sx, sy] = hit_and_run(vx, vy, N, 50);
29-fprintf('hit_and_run: N=%d sampled in %.3f s\n', N, toc);
25+% Build a region (convex by default) and draw N uniform samples from it.
26+convex = true;
27+[vx, vy] = make_region(convex);
28+[sx, sy] = sample_region(vx, vy, N, convex);
3029
3130 % The figure app is the prebuilt single-file page, relative to the project root
3231 % (the current working directory when a top-level script is run).
@@ -35,45 +34,60 @@ html = fileread(fullfile('app', 'dist', 'index.html'));
3534 fig = figure;
3635 gl = uigridlayout(fig, [1 1], 'Padding', [0 0 0 0], ...
3736 'RowHeight', {'1x'}, 'ColumnWidth', {'1x'});
38-uihtml(gl, 'HTMLSource', html, 'Data', pack_data(vx, vy, sx, sy, N), ...
37+uihtml(gl, 'HTMLSource', html, 'Data', pack_data(vx, vy, sx, sy, N, convex), ...
3938 'HTMLEventReceivedFcn', @(src, ev) on_event(src, ev));
4039 end
4140
4241 function on_event(src, ev)
4342 % Figure -> script. Two requests the controls send:
44-% 'resample' {n, x, y} -> draw n fresh samples in the given region; reply
45-% with a 'samples' event (region is unchanged).
46-% 'newRegion' {n} -> build a new region, draw n samples in it; reply
47-% with a full 'data' event (region + samples).
43+% 'resample' {n, x, y, convex} -> draw n fresh samples in the given region;
44+% reply with a 'samples' event.
45+% 'newRegion' {n, convex} -> build a new region (convex or non-convex),
46+% draw n samples; reply with a 'data' event.
4847 d = ev.HTMLEventData;
4948 n = 10000;
50-if isstruct(d) && isfield(d, 'n')
51- n = max(1, round(d.n));
49+convex = true;
50+if isstruct(d)
51+ if isfield(d, 'n')
52+ n = max(1, round(d.n));
53+ end
54+ if isfield(d, 'convex')
55+ convex = logical(d.convex);
56+ end
5257 end
5358 switch ev.HTMLEventName
5459 case 'resample'
5560 vx = d.x(:);
5661 vy = d.y(:);
57- tic;
58- [sx, sy] = hit_and_run(vx, vy, n, 50);
59- fprintf('resample: N=%d sampled in %.3f s\n', n, toc);
62+ [sx, sy] = sample_region(vx, vy, n, convex);
6063 sendEventToHTMLSource(src, 'samples', pack_samples(sx, sy, n));
6164 case 'newRegion'
62- [vx, vy] = make_region();
63- tic;
64- [sx, sy] = hit_and_run(vx, vy, n, 50);
65- fprintf('newRegion: N=%d sampled in %.3f s\n', n, toc);
66- sendEventToHTMLSource(src, 'data', pack_data(vx, vy, sx, sy, n));
65+ [vx, vy] = make_region(convex);
66+ [sx, sy] = sample_region(vx, vy, n, convex);
67+ sendEventToHTMLSource(src, 'data', pack_data(vx, vy, sx, sy, n, convex));
68+end
69+end
70+
71+function [sx, sy] = sample_region(vx, vy, n, convex)
72+% Convex regions use the fast JIT-compiled chord sampler; non-convex regions
73+% use the general multi-segment sampler.
74+tic;
75+if convex
76+ [sx, sy] = hit_and_run(vx, vy, n, 50);
77+else
78+ [sx, sy] = hit_and_run_general(vx, vy, n, 50);
6779 end
80+fprintf('sample_region: N=%d convex=%d in %.3f s\n', n, convex, toc);
6881 end
6982
70-function data = pack_data(vx, vy, sx, sy, N)
83+function data = pack_data(vx, vy, sx, sy, N, convex)
7184 %PACK_DATA Full payload (region + samples) sent once when the figure opens.
7285 data = struct();
7386 data.type = 'hitandrun';
7487 data.region = struct('x', vx(:).', 'y', vy(:).');
7588 data.samples = struct('x', sx(:).', 'y', sy(:).');
7689 data.n = N;
90+data.convex = convex;
7791 end
7892
7993 function s = pack_samples(sx, sy, N)