diff --git a/.github/workflows/cibuildwheel.yml b/.github/workflows/cibuildwheel.yml
index 8c71da3f..2ce2a956 100644
--- a/.github/workflows/cibuildwheel.yml
+++ b/.github/workflows/cibuildwheel.yml
@@ -75,7 +75,7 @@ jobs:
           # regression tests unsuited to wheel smoke.
           CIBW_TEST_COMMAND: >
             python {project}/scripts/ci_test_setup.py {project} &&
-            python -c "import openptv2.algorithms.track_kernels_tracking, openptv2.algorithms.track3d; print('openmp import OK')" &&
+            python -c "import openptv2.algorithms.track_kernels_corr, openptv2.algorithms.track3d; print('openmp import OK')" &&
             pytest
             tests/unit/test_vec_utils.py
             tests/unit/test_correspondences.py
diff --git a/docs/index.md b/docs/index.md
index d8d55295..ba1aa091 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -13,6 +13,8 @@ Welcome to the openptv2 documentation.
 ### User Documentation
 
 - [Tracking Pipeline & Results Guide](tracking_guide.md) - Pipeline workflow, parameter guide, multi-pass tracking, and ptv_is.# output format
+- [Particle Trackers](trackers.md) - Which tracker to use and why: basics, all eight engines, caveats, tips, upstream credits
+- [Two-Phase Tracking](two-phase-tracking.md) - Two-Phase usage and parameters in depth
 - [Lid-Driven Cavity Flow Tutorial](tutorials/cavity_flow_tutorial.md) - End-to-end 3D-PTV case study: Autocalibration, Tracer Shaking, Warmup, and 3D Trajectories
 - [Aortic Pulsatile Flow Tutorial](aorta_tutorial.md) - Cloud-native 3D-PTV on complex aortic flow
 - [Auto-Calibration with `openptv warmup`](tutorials/warmup_tutorial.md) - Standalone parameter/engine auto-tuning before tracking
diff --git a/docs/plans/2026-09-20-tracking-kernels-dedup.md b/docs/plans/2026-09-20-tracking-kernels-dedup.md
new file mode 100644
index 00000000..3e4cd6a5
--- /dev/null
+++ b/docs/plans/2026-09-20-tracking-kernels-dedup.md
@@ -0,0 +1,89 @@
+# Tracking-kernels dedup (ponytail-audit items 1-9)
+
+Branch: `refactor/dedup-tracking-kernels` (off `main`). Committed and pushed (0e1200a0, 1b83cfa6); no PR yet.
+Goal: remove duplicated / orphaned code in `src/openptv2/algorithms/track_kernels_*.py`
+without changing behavior.
+
+## What was found
+
+The same Cython `cdef`/`nogil` kernels were copy-pasted across sibling modules
+(cross-module C calls need a `.pxd`, so people copied instead). Bodies were
+identical or differed only in decorators/comments. `track_kernels_pixel.pxd` and
+`track_kernels_position.pxd` already existed as the sharing mechanism, so I
+extended them rather than inventing anything.
+
+## Done
+
+Items 1-5, dedup. Owners:
+- `track_kernels_pixel.py`: `_multimed_r_nlay_1layer` (added `exceptval(check=False)` + pxd entry),
+  `_point_to_pixel_out`, `_candsearch_in_pix_rest_nogil`, `_pixel_to_metric_out`,
+  `_dist_to_flat_out`, `_sorted_candidates_fast_out_nogil`, `candsearch_in_pix_fast_nogil`
+  (the pixel copy is the newer one with `max_cands`/`out_dists`; search's was a stale, dead copy).
+- `track_kernels_position.py`: `_ray_tracing_out` (added exceptval + pxd), `_angle_acc_out`
+  (moved here from geom/corr, added exceptval + pxd), `_point_position_out`,
+  `assess_new_position_fast_nogil`.
+- Copies deleted from geom, search, transform, corr. Importers use the existing
+  `if cython.compiled: cimport ... else: from .x import ...` pattern.
+  Import direction is acyclic: pixel <- position <- geom/transform/corr/batch.
+
+Item 6: deleted the `track_kernels_tracking.py` shim. Repointed `track_kernels.py`,
+`setup.py` (`ALGORITHMS_MODULES`), `.github/workflows/cibuildwheel.yml` (import smoke test),
+and `tests/unit/test_track_kernels_tracking_coverage.py` (now imports from owners; shim
+constants defined locally; `_mod` = the pixel module, which keeps the old inert-patch behavior).
+
+Item 9: deleted 16 wrappers with no non-test callers, plus their tests and the two
+re-exports in `track_kernels.py`:
+`angle_acc_fast`, `_ray_tracing_fast`, `pixel_to_metric_fast`, `dist_to_flat_fast`,
+`metric_to_pixel_fast`, `_metric_to_pixel_out`, `_flat_image_coord_fast`, `_img_coord_fast`,
+`img_coord_batch_fast`, `flat_image_coord_batch_fast`, `point_position_fast`,
+`ray_tracing_batch_fast`, `point_position_batch_fast`, `pixel_to_metric_batch_fast`,
+`metric_to_pixel_batch_fast`, `sort_candidates_by_freq_fast`.
+Also dropped two search test classes that only tested the deleted dead copies.
+
+Net so far: roughly -2,200 lines from dedup, about -840 from item 9, plus tests.
+
+## Verification status (updated 2026-09-21)
+
+Clean Cython rebuild of the final tree: OK.
+- Hot-path tests (`test_track`, `test_track3d`, `test_correspondences`, `test_track4be`): 49 passed = baseline.
+- Full suite `uv run --no-sync pytest tests`: 2017 passed, 86 skipped, 39 deselected (12 min).
+- Pure-Python fallback, kernel coverage files only
+  (`test_track_kernels_*_coverage.py`, 6 files): 266 passed.
+- Pure-Python fallback over the whole `tests/unit/test_*_coverage.py` glob: 1227 passed, 24 failed
+  (38 min; the glob now matches 29 files, not the 16 CLAUDE.md mentions, so it is slow).
+  The failures I inspected (`test_epi_coverage`, `test_correspondences_coverage`, 14 of the 24)
+  are all `Coord2d.__init__() got an unexpected keyword argument 'pnr'` /
+  `Candidate.__init__() ... 'pnr'`: interpreted-mode constructor mismatch in modules this
+  branch does not touch. Not confirmed on `main`; the other 10 were not inspected.
+- Running the suite rewrites tracked `test_data/test_cavity/img/*_targets`; `git checkout -- test_data` before committing.
+
+## Remaining
+
+1. Optional: confirm the 24 fallback failures also occur on `main` (build `main`, run
+   `tests/unit/test_epi_coverage.py tests/unit/test_correspondences_coverage.py` interpreted).
+2. Optional perf sanity: `_angle_acc_out` was `ccall inline` inside corr and is now a cross-module
+   C call (lost inlining). Time a tracking run before/after; if it regressed, keep a private copy in corr.
+3. Cosmetic: stray banner comments in `test_track_kernels_batch_coverage.py` (about lines 84-102) and
+   `test_track_kernels_transform_coverage.py` (lines 22, 352).
+4. Open a PR.
+
+## Deliberately skipped
+
+- **Item 7** (merge forward/backward tracking loops in `track_kernels_corr.py`): not a clean
+  dedup. Normalised diff of `trackcorr_loop_fast` vs `trackback_loop_fast` shows about 670 of about 900
+  lines differ; forward uses the `_trackcorr_particle_fast` worker, backward is inline. Merging
+  changes core tracking logic and there is no golden-output regression data. Only attempt
+  with a recorded before/after trajectory comparison on a real dataset.
+- **Item 8** (delete 4BE tracker: `track4be_loop_fast`, `track4be.py`, `plugins/four_be_tracking.py`):
+  it is registered in `tracking_registry.py` and benchmarked in about 10 scripts/notebooks
+  (`bench_*`, `benchmark_*`, `tracker_tutorial_dashboard.py`). Needs an explicit product decision.
+  Do not delete without Alex saying so.
+
+## Gotchas
+
+- Cimported names are not importable from Python: tests and non-cimporting modules must import
+  the `cpdef` names from the owner module (pixel/position), not from geom/transform.
+- Cython pure mode + `.pxd`: `noexcept nogil` in the pxd requires `@cython.exceptval(check=False)`
+  on the `def`, or the signatures mismatch.
+- Shell cwd drift: `cd` to the repo root at the start of every command.
+- Use `git grep`, not `grep -r`, at the repo root (huge build dirs; a plain grep timed out).
diff --git a/docs/trackers.md b/docs/trackers.md
new file mode 100644
index 00000000..2ce7f169
--- /dev/null
+++ b/docs/trackers.md
@@ -0,0 +1,156 @@
+# Particle Trackers in openptv2
+
+Start here if you want to know which tracker to use and why. For measured
+numbers on a reference dataset, see [Tracker Tutorials](tracker-tutorials.md);
+for pipeline and file formats, see the [Tracking Guide](tracking_guide.md);
+for Two-Phase usage in depth, see [Two-Phase Tracking](two-phase-tracking.md).
+
+## 1. The basics: what tracking is
+
+A camera takes a picture every frame. In each picture, every particle is a
+dot. In the next picture the dots have moved a little. **Tracking means
+deciding which dot in picture 2 is the same particle as each dot in
+picture 1.**
+
+When particles are far apart and slow, this is easy: take the nearest dot.
+It gets hard when particles cross each other, when one is hidden for a
+frame, or when positions are noisy. Every tracker below is a different set
+of rules for that decision. They differ in only four things:
+
+1. **Guess** — where do you expect the particle next? (stay put, constant
+   velocity, smoothed curve, look at future frames, bigger time step)
+2. **Score** — how do you rank candidates? (distance, acceleration, angle,
+   camera-image agreement)
+3. **Conflicts** — two particles want the same dot: who wins? (first come,
+   cheapest first, best total pairing, nobody)
+4. **Evidence** — what do you look at? Only 3D points, or also the original
+   camera images? Trackers that check the images survive bad 3D points;
+   trackers that ignore them are faster.
+
+Two facts shape everything on this page. First, measured on synthetic flow
+(see `docs/tracking-benchmark-results.md`), most wrong links come from
+**missing detections across gaps (~60%)** and **cold starts (~20%)** — not
+from the scoring rule. So gap handling matters more than clever costs.
+Second, particles only ever **enter and exit** the volume; a track that ends
+mid-volume is almost always an occlusion or a dropout, not a real exit.
+
+## 2. All trackers at a glance
+
+Select with `plugins.selected_tracking` in your parameters (GUI: Plugins
+page). The `name` column is the exact preset string.
+
+| preset (`name`) | idea in one line | looks at |
+|---|---|---|
+| `default`, `standard_forward`, `full_multipass`, `two_directional` (trackcorr) | Guess forward, confirm in every camera, accept smooth links | 2D targets + 3D |
+| `priority_segment_3d` (Fast 3D / 3MA) | Cheapest (smoothest) links first, globally | 3D only |
+| `4be` | Peek at frame n+2 before accepting; conflicts link to nobody | 3D only |
+| `nearest_hungarian_3d` (MyPTV 3D) | Best total pairing per frame pair (Hungarian), survives gaps | 3D only |
+| `myptv_2d_tracking` (MyPTV 2D) | Each camera tracks its own images, then cameras vote | 2D per camera |
+| `predictive_gmm_3d` (proPTV) | Fit a smooth curve through history, predict from it | 3D only |
+| `two_phase` (Two-Phase) | 3D search for candidates, per-camera images for ranking | 3D + 2D |
+| `hybrid_deltat_3d` (Hybrid) | Match every N-th frame where motion beats noise, fill between | 3D only |
+
+Details per tracker live in `src/openptv2/tracking_registry.py`
+(`TRACKER_REGISTRY`) — the machine-readable version of this page.
+
+## 3. Each tracker: caveats, tips, tricks
+
+### trackcorr (`default`, `standard_forward`, …) — the original engine
+The safest default. Predicts forward, checks candidates in every camera
+image, accepts links with small acceleration and turning angle, and can run
+backward plus gap relinking. **Caveat:** slowest of the bunch, and three
+interacting parameters (`dvxmax`, `dacc`, `angle`). **Tip:** on noisy slow
+flow, switch the angle limit off — the measured turning angle is mostly
+noise then. If a particle is hidden in one camera, trackcorr usually still
+gets it through the others.
+
+### Fast 3D / 3MA (`priority_segment_3d`) — the fast one
+Hundreds of thousands of particles per second, simplest mental model
+(smoothest join wins). **Caveat:** blind to the images, so a ghost particle
+sitting where the guess expects is accepted without question. **Tip:** use
+it for clean, dense data and high throughput; distrust it where ghosts are
+likely (poor calibration, few cameras).
+
+### 4BE (`4be`) — the careful one
+Looks one frame into the future before committing, and refuses contested
+dots outright. **Caveat:** built for sparse clean data; on noisy/dense data
+it produces the worst accelerations of all engines and leaves gaps
+deliberately (gap bridging is off for its preset on purpose — it would
+rebuild exactly the links 4BE declined). **Tip:** sparse lab data with
+reliable detection; not turbulence.
+
+### MyPTV 3D (`nearest_hungarian_3d`) — the fair one
+Nobody grabs the nearest dot first: it finds the pairing with the lowest
+*total* distance, so nobody is paired badly. Tracks survive short gaps
+(`max_gap`), code is plain readable Python. **Caveat:** one frame pair at a
+time — it cannot use what happens next. **Tip:** good first alternative to
+the default; easiest engine to modify (`src/openptv2/plugins/myptv_3d_tracking.py`).
+
+### MyPTV 2D (`myptv_2d_tracking`) — the per-camera voter
+Each camera tracks its own movie; links with the most camera votes win.
+**Caveat:** a link needs only one vote, so a single confused camera can
+still create a bad link. **Tip:** reaches for it when 3D triangulation is
+unreliable but the raw images are clean.
+
+### proPTV (`predictive_gmm_3d`) — the smoother
+Fits smooth curves through each path, so speeds and accelerations stay
+sensible under noise. **Caveat (read first):** this port predicts from the
+*smoothed current position*, not from an extrapolated one, and its search
+radius is the same with or without history — so it currently behaves closer
+to smoothed nearest-neighbour than to the predictive scheme of the paper.
+Also dense data gets expensive. **Tip:** check the plugin README before
+trusting its "predictive" label; fix the extrapolation first if you build
+on it.
+
+### Two-Phase (`two_phase`) — the hybrid
+3D search lists candidates, per-camera image distances rank them, Hungarian
+per connected group decides. No motion model — immune to bad guesses, but
+fails once motion outruns particle spacing. Reported +74% multi-frame
+trajectories over the default on a poorly-conditioned aorta dataset.
+**Caveat:** needs Cet 2D targets per camera in the store; falls back to
+pure 3D (`leaf_weight=0`) without them. **Tips:** start with
+`leaf_weight=1`, `v_max` at ~3× your typical step; see
+[Two-Phase Tracking](two-phase-tracking.md) for the full parameter guide
+including the shared-observation prototype flags.
+
+### Hybrid multi-Δt (`hybrid_deltat_3d`) — the slow-flow specialist
+When particles crawl, frame-to-frame steps drown in noise — so it matches
+every N-th frame (where displacement beats noise) and fills the middle with
+a smooth curve. **Caveat:** the only engine that changes the signal-to-noise
+ratio instead of fighting it, but the smooth fill is wrong for fast or
+curved motion. **Tip:** high frame rate + slow flow; set `stride` so the
+coarse step clearly exceeds your 3D noise floor.
+
+## 4. Upstream credit: MyPTV and proPTV are plugins, not forks
+
+Two engines borrow ideas from outside projects. **We use them as plugins —
+adapted concepts on openptv2's own data structures, not modified copies of
+their code.** The full frameworks (triangulation pipelines, calibration,
+backtracking/repair, smoothing toolboxes) live only in their own
+repositories — use those projects directly if you need them. Both are
+permissively MIT-licensed.
+
+- **MyPTV** by Ron Shnapp — open-source Python 3D-PTV library.
+  Repository: <https://github.com/ronshnapp/MyPTV> ·
+  Paper: Shnapp, R. (2022). *MyPTV: A Python Package for 3D Particle
+  Tracking.* Journal of Open Source Software, 7(75), 4398.
+  <https://doi.org/10.21105/joss.04398> ·
+  What we adapted: per-camera 2D image-space tracking with multi-camera
+  consensus (`myptv_2d_tracking`), and kinematic prediction + assignment
+  matching in 3D (`nearest_hungarian_3d`).
+- **proPTV** by Robin Barta and colleagues (DLR) — probabilistic PTV
+  framework, Python.
+  Repository: <https://github.com/RobinBarta/proPTV> ·
+  Paper: Barta, R. et al. (2024). *proPTV: A probabilistic particle
+  tracking velocimetry framework.* Journal of Computational Physics, 514,
+  113212. <https://doi.org/10.1016/j.jcp.2024.113212> ·
+  What we adapted: the small pure-numpy core (Gaussian-mixture / basis
+  approximation and Savitzky–Golay smoothing, vendored under
+  `src/openptv2/plugins/proptv/`), wired into the `predictive_gmm_3d`
+  plugin. The original's triangulation, probability model, backtracking
+  and repair are *not* ported.
+- The classic engines descend from the OpenPTV/liboptv lineage
+  (<http://www.openptv.net>).
+
+If you publish with these methods, please cite the original authors above
+in addition to openptv2.
diff --git a/docs/two-phase-tracking.md b/docs/two-phase-tracking.md
new file mode 100644
index 00000000..97a9f933
--- /dev/null
+++ b/docs/two-phase-tracking.md
@@ -0,0 +1,88 @@
+# Two-Phase Tracking: usage and parameters
+
+Two-Phase (`selected_tracking: two_phase`,
+`src/openptv2/plugins/two_phase_tracking.py`) links particles in two steps:
+**Phase 1** finds candidates with a 3D KD-tree around each track's predicted
+position; **Phase 2** ranks them by per-camera 2D image distance and solves a
+Hungarian assignment per connected group. 3D proposes, the images dispose.
+
+## 1. How to use it
+
+**GUI:** Plugins page → tracking plugin → `two_phase`. Parameters live in
+the `track` section (same names as below).
+
+**Batch/YAML:** minimal setup —
+```yaml
+plugins:
+  selected_tracking: two_phase
+track:
+  v_max: 5.0        # mm/frame; or dvxmax, same meaning here
+  leaf_weight: 1.0
+  max_gap: 2
+```
+
+**Plain Python** (no experiment needed):
+```python
+import numpy as np
+from openptv2.plugins.two_phase_tracking import (
+    TwoPhaseTracker, TwoPhaseTrackerConfig)
+
+cfg = TwoPhaseTrackerConfig(v_max=5.0, leaf_weight=1.0, max_gap=2)
+links = TwoPhaseTracker(cfg).track_frames(
+    frame_particles,   # list of (N_i, 3) arrays, mm
+    frame_leaves,      # list of (N_i, 2*C) arrays, px (optional)
+    project_fn,        # (N,3) -> (N,2*C) re-projection (optional)
+)
+# links: list of (t0, row0, t1, row1); use return_chains=True to also get
+# per-track histories: links, chains = ...track_frames(..., return_chains=True)
+```
+
+Without leaves/`project_fn` it falls back to pure 3D distance costs
+(`cost_mode="3d"`); pass `leaf_weight=0` to force that explicitly.
+
+## 2. Parameters
+
+All live in the `track` YAML section (batch/GUI) or on
+`TwoPhaseTrackerConfig` (Python). Units in brackets.
+
+| parameter | default | what it does | how to set it |
+|---|---|---|---|
+| `v_max` (or `dvxmax`) [mm/frame] | 15.5 | 3D search radius around each prediction | ~3× your typical per-frame step. Too small: true links never become candidates. Too big: everything connects into giant groups (slow, sloppy) |
+| `leaf_weight` [–] | 1.0 | weight of 2D image distance in the cost | 1.0 normally; 0 = pure 3D (no leaves needed). Lower it when calibration is poor — bad projection poisons the ranking |
+| `use_velocity` [bool] | true | match predictions (`pos + vel·dt`), not positions | Keep on. Off = every crossing resolves as a bounce |
+| `cost_mode` [str] | `projected` | `projected`: rank by re-projected 2D distance (needs `project_fn`/cals); `3d`: rank by 3D distance | `projected` with good calibration, `3d` otherwise |
+| `max_gap` [frames] | 2 | a track survives this many unmatched frames | 2 covers single-frame dropouts (the biggest measured failure source). Higher = longer bridges, more impostors |
+| `dt` [–] | 1.0 | time step for velocity | 1.0 for consecutive frames |
+| `max_group_size` [nodes] | 128 | groups bigger than this skip the cubic Hungarian, greedy inside | Raise only if you can afford it; at production density frames percolate and this cap is what keeps a run from stalling |
+| `allow_shared` [bool] | false | **prototype:** in groups with more tracks than detections (occlusion), losers share the winner's detection instead of dying | Enable where occlusions matter; validated on synthetic crossings (0 switches). Shared points move position but never velocity |
+| `max_shared` [frames] | 2 | max consecutive shared frames per track | 2 covers brief overlaps; higher risks twin tracks that never separate |
+| `share_tol` [cost] | 1.0 | a shared claim needs an edge cost below this (mutual-prediction gate) | Without it, a stranded track hijacks strangers' detections (observed live). ~5–10× your position noise; `null` disables the gate (not advised) |
+
+## 3. How it behaves (caveats)
+
+- **No motion model to be wrong** — but also none to help: if particles move
+  farther per frame than the typical spacing, every tracker fails, this one
+  first. Check step-vs-spacing before blaming parameters.
+- **Gap survival is prediction-based:** a track coasts on `pos + vel`
+  through gaps up to `max_gap`; a maneuver inside the gap is lost. That is
+  by design — see `max_gap` above.
+- **Occlusions:** with `allow_shared`, one detection may serve two tracks
+  for up to `max_shared` frames (marked in chains when
+  `return_chains=True`). Owners of the other engines: this is the reference
+  implementation of the shared-observation rule — same idea ports to
+  `track3d_loop_fast` (recorded, not claimed) and to linkage postprocess
+  (mark, then bridge).
+- **Speed:** ~linear in particles (KD-tree search, small per-group
+  Hungarians); the `max_group_size` cap bounds the worst case.
+
+## 4. Tuning recipe
+
+1. Measure your typical step `s` (median linked displacement) and 3D noise
+   `n` (second-difference statistics — see `docs/algorithms/tracking.md`).
+2. `v_max` ≈ 3·s. `leaf_weight` = 1 with decent calibration, 0 without.
+3. Run; count short tracks (cold starts) vs gaps. More gaps than tracks →
+   raise `max_gap` to 3. More fragments at crossings → enable
+   `allow_shared`.
+4. Validate on synthetic ground truth with the same spacing/noise before
+   trusting a production run (`tests/helpers/synthetic_scene.py` generates
+   scenes; `scripts/proto_shared_validate.py` shows the scoring pattern).
diff --git a/scripts/proto_long_validate.py b/scripts/proto_long_validate.py
new file mode 100644
index 00000000..b5c83a96
--- /dev/null
+++ b/scripts/proto_long_validate.py
@@ -0,0 +1,239 @@
+"""Long-window validation (50 real frames) for the three shared prototypes.
+
+Data: HiDImaging wp1 test, frames 100001-100050, ~1000 pts/frame.
+Input: res_orig/rt_is.* 3D points. Reference: res_orig/ptv_is.* linkages
+(0-based prev/next chains) = the res_orig tracking output itself.
+Prototypes (worktree code): TwoPhase +/-share, Fast3D +/-share_tol,
+greedy baseline vs +mark/assemble (trackcorr-linkage path).
+Metrics: tracks, mean length, coverage, fragmentation of reference tracks,
+impurity of prototype tracks (switch/merge proxy), shared events, runtime.
+"""
+import sys
+import time
+
+import numpy as np
+
+TEST = r"C:\Users\alex\Downloads\HiDImaging\CompleteTest\wp1\test"
+F0, NF = 100001, 50
+TOL = 1.0  # mm
+
+
+def load_points():
+    frames = []
+    for f in range(F0, F0 + NF):
+        rows = []
+        with open(f"{TEST}/res_orig/rt_is.{f}") as fh:
+            n = int(fh.readline().split()[0])
+            for _ in range(n):
+                t = fh.readline().split()
+                rows.append([float(t[1]), float(t[2]), float(t[3])])
+        frames.append(np.array(rows))
+    return frames
+
+
+def load_linkages():
+    prev, nxt, pos = {}, {}, {}
+    for f in range(F0, F0 + NF):
+        with open(f"{TEST}/res_orig/ptv_is.{f}") as fh:
+            n = int(fh.readline().split()[0])
+            pv = np.full(n, -1, np.int32)
+            nx = np.full(n, -2, np.int32)
+            xy = np.zeros((n, 3))
+            for i in range(n):
+                t = fh.readline().split()
+                pv[i], nx[i] = int(t[0]), int(t[1])
+                xy[i] = [float(t[2]), float(t[3]), float(t[4])]
+        prev[f], nxt[f], pos[f] = pv, nx, xy
+    return prev, nxt, pos
+
+
+def assemble_ref(prev, nxt, pos):
+    chains = []
+    visited = set()
+    for f in range(F0, F0 + NF):
+        for i in range(len(pos[f])):
+            if (f, i) in visited or int(prev[f][i]) >= 0:
+                continue
+            chain = []
+            cf, ci = f, i
+            while True:
+                if (cf, ci) in visited:
+                    break
+                visited.add((cf, ci))
+                chain.append((cf, ci))
+                if cf not in nxt or ci >= len(nxt[cf]):
+                    break
+                ni = int(nxt[cf][ci])
+                if ni < 0 or cf + 1 not in pos or ni >= len(pos[cf + 1]):
+                    break
+                cf += 1
+                ci = ni
+            if len(chain) >= 1:
+                chains.append(chain)
+    # orphans (all-prev/no-next singletons not yet visited)
+    for f in range(F0, F0 + NF):
+        for i in range(len(pos[f])):
+            if (f, i) not in visited:
+                chains.append([(f, i)])
+                visited.add((f, i))
+    return chains, pos
+
+
+def calibrate(chains, pos):
+    steps = []
+    for c in chains:
+        for a in range(1, len(c)):
+            f0, i0 = c[a - 1]
+            f1, i1 = c[a]
+            if f1 == f0 + 1:
+                steps.append(float(np.linalg.norm(pos[f1][i1] - pos[f0][i0])))
+    steps = np.array(steps)
+    vm = float(np.percentile(steps, 99))
+    print(f"calibrate: {len(chains)} ref chains, step p50="
+          f"{np.median(steps):.3f} p99={vm:.3f}", flush=True)
+    return vm
+
+
+def score(proto_tracks, ref_chains, pos):
+    from scipy.spatial import cKDTree
+    # frame -> (points, proto_tid) / (points, ref_chain_idx) indexes
+    ptrees, rtrees = {}, {}
+    pmap, rmap = {}, {}
+    for ti, t in enumerate(proto_tracks):
+        for f, p in zip(t["frames"], t["pos"]):
+            pmap.setdefault(int(f), []).append((np.asarray(p), ti))
+    for ci, c in enumerate(ref_chains):
+        for (f, i) in c:
+            rmap.setdefault(f, []).append((pos[f][i], ci))
+    for f in set(list(pmap) + list(rmap)):
+        if f in pmap:
+            pts = np.array([p for p, _ in pmap[f]])
+            ptrees[f] = (cKDTree(pts), [ti for _, ti in pmap[f]])
+        if f in rmap:
+            pts = np.array([p for p, _ in rmap[f]])
+            rtrees[f] = (cKDTree(pts), [ci for _, ci in rmap[f]])
+    # fragmentation: ref chains split across >1 proto tracks
+    frag, frag_den, cover_pts, cover_den = 0, 0, 0, 0
+    for c in [c for c in ref_chains if len(c) >= 3]:
+        owners = set()
+        for (f, i) in c:
+            cover_den += 1
+            if f not in ptrees:
+                continue
+            tree, tids = ptrees[f]
+            d, k = tree.query(pos[f][i], k=1, distance_upper_bound=TOL)
+            if d <= TOL:
+                owners.add(tids[int(k)])
+                cover_pts += 1
+        frag_den += 1
+        if len(owners) > 1:
+            frag += 1
+    # impurity: proto tracks spanning >1 ref chain
+    impure = 0
+    for t in proto_tracks:
+        if len(t["frames"]) < 2:
+            continue
+        owners = set()
+        for f, p in zip(t["frames"], t["pos"]):
+            if int(f) not in rtrees:
+                continue
+            tree, cis = rtrees[int(f)]
+            d, k = tree.query(np.asarray(p), k=1, distance_upper_bound=TOL)
+            if d <= TOL:
+                owners.add(cis[int(k)])
+        if len(owners) > 1:
+            impure += 1
+    ntrk = sum(1 for t in proto_tracks if len(t["frames"]) >= 2)
+    ml = np.mean([len(t["frames"]) for t in proto_tracks
+                  if len(t["frames"]) >= 2]) if ntrk else 0.0
+    return {"ntrk": ntrk, "meanlen": round(float(ml), 2),
+            "cover": round(cover_pts / max(cover_den, 1), 3),
+            "frag%": round(100 * frag / max(frag_den, 1), 1),
+            "impure": impure}
+
+
+def main():
+    t0 = time.perf_counter()
+    frames = load_points()
+    prev, nxt, pos = load_linkages()
+    ref_chains, _ = assemble_ref(prev, nxt, pos)
+    v_max = calibrate(ref_chains, pos)
+    print(f"load done ({time.perf_counter() - t0:.1f}s)", flush=True)
+    ref_tracks = [{"frames": [f for f, _ in c],
+                   "pos": np.array([pos[f][i] for f, i in c])}
+                  for c in ref_chains]
+    print(f"res_orig: {score(ref_tracks, ref_chains, pos)}", flush=True)
+
+    from openptv2.plugins.two_phase_tracking import (
+        TwoPhaseTracker, TwoPhaseTrackerConfig)
+
+    for share in [False, True]:
+        t1 = time.perf_counter()
+        cfg = TwoPhaseTrackerConfig(v_max=v_max, max_gap=2, dt=1.0,
+                                    leaf_weight=0.0, cost_mode="3d",
+                                    use_velocity=True, allow_shared=share)
+        _, chains = TwoPhaseTracker(cfg).track_frames(
+            [np.asarray(p) for p in frames], return_chains=True)
+        tracks = [{"frames": [F0 + f for f in c["frames"]], "pos": c["pos"]}
+                  for c in chains if len(c["frames"]) >= 1]
+        ns = sum(sum(c["shared"]) for c in chains)
+        print(f"twophase share={share}: shared_pts={ns} "
+              f"{score(tracks, ref_chains, pos)} "
+              f"({time.perf_counter() - t1:.1f}s)", flush=True)
+
+    sys.path.insert(0, "scripts")
+    from scipy.spatial import cKDTree
+    from proto_shared_validate import fast3d_tracks, trackcorr_linkage_tracks
+
+    # Fast3D pure-Python kernel: 15 frames (50 would take too long in
+    # interpreted mode; the compiled path is unaffected).
+    NF_FAST = 15
+    pos15 = {f: pos[f] for f in range(F0, F0 + NF_FAST)}
+    chains15 = []
+    for c in ref_chains:
+        cc = [(f, i) for (f, i) in c if F0 <= f < F0 + NF_FAST]
+        if cc:
+            chains15.append(cc)
+    for tol in [0.0, 0.5]:
+        t1 = time.perf_counter()
+        tr, ns = fast3d_tracks([np.asarray(p) for p in frames[:NF_FAST]],
+                               v_max=v_max, share_tol=tol)
+        for t in tr:
+            t["frames"] = [F0 + f for f in t["frames"]]
+        print(f"fast3d[{NF_FAST}f] tol={tol}: shared={ns} "
+              f"{score(tr, chains15, pos15)} "
+              f"({time.perf_counter() - t1:.1f}s)", flush=True)
+
+    def greedy_links_kd(frame_particles, gate):
+        links = []
+        for t in range(len(frame_particles) - 1):
+            p0 = np.asarray(frame_particles[t])
+            p1 = np.asarray(frame_particles[t + 1])
+            if len(p0) == 0 or len(p1) == 0:
+                continue
+            tree = cKDTree(p1)
+            dists, idxs = tree.query(p0, k=1, distance_upper_bound=gate)
+            used = set()
+            for i in range(len(p0)):
+                j = int(idxs[i])
+                if j < len(p1) and j not in used and dists[i] <= gate:
+                    used.add(j)
+                    links.append((t, i, t + 1, j))
+        return links
+
+    t1 = time.perf_counter()
+    base = greedy_links_kd([np.asarray(p) for p in frames], gate=v_max)
+    plain, healed, nm = trackcorr_linkage_tracks(
+        [np.asarray(p) for p in frames], base, tol=1.0)
+    for name, tr in [("plain", plain), ("healed", healed)]:
+        for t in tr:
+            t["frames"] = [F0 + f for f in t["frames"]]
+        print(f"tclink {name}: marks={nm if name == 'healed' else 0} "
+              f"{score(tr, ref_chains, pos)} "
+              f"({time.perf_counter() - t1:.1f}s)", flush=True)
+
+    print("done.")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/proto_shared_validate.py b/scripts/proto_shared_validate.py
new file mode 100644
index 00000000..d4b69a24
--- /dev/null
+++ b/scripts/proto_shared_validate.py
@@ -0,0 +1,279 @@
+"""Shared-observation prototype validation: one 3D scene, three trackers.
+
+Scene (10 frames, mm units, constant velocity + jitter):
+  P0/P1: crossing pair, frames 4-5 occluded -> single midpoint point.
+  P2:    dropout gap at frames 3-4 (no point at all).
+  P3:    exits the volume after frame 6 (legal end).
+  P4/P5: clean crucirng traffic.
+Ground truth = 6 trajectories. Metrics: tracks, switches, missing, purity.
+"""
+import sys
+
+import numpy as np
+
+N_FRAMES = 10
+
+
+def make_scene(seed=0, jitter=0.05):
+    rng = np.random.default_rng(seed)
+    vel = {
+        0: np.array([2.0, 0.4, 0.0]),
+        1: np.array([2.0, -0.4, 0.0]),
+        2: np.array([1.0, 1.5, 0.2]),
+        3: np.array([3.0, 0.0, 0.0]),
+        4: np.array([-1.5, 1.0, 0.1]),
+        5: np.array([0.5, -1.2, -0.1]),
+    }
+    p0 = {
+        0: np.array([0.0, -2.0, 0.0]),
+        1: np.array([0.0, 2.0, 0.0]),
+        2: np.array([5.0, 5.0, 1.0]),
+        3: np.array([0.0, 10.0, -1.0]),
+        4: np.array([20.0, 0.0, 2.0]),
+        5: np.array([8.0, 12.0, 0.5]),
+    }
+    # P0/P1 cross near x=10 around frame 5: shift so midpoint lands ~frame 4-5
+    frames_pts, frames_ids = [], []  # per frame: (M,3) points, (M,) truth ids
+    for f in range(N_FRAMES):
+        pts, ids = [], []
+        for i in range(6):
+            if i == 3 and f > 6:
+                continue  # legal exit
+            if i == 2 and f in (3, 4):
+                continue  # dropout gap
+            pos = p0[i] + vel[i] * f + rng.normal(0, jitter, 3)
+            pts.append(pos)
+            ids.append(i)
+        pts = np.array(pts)
+        ids = np.array(ids)
+        if f in (4, 5):
+            # occlusion: P0 and P1 merge into their midpoint (single point)
+            m0 = np.where(ids == 0)[0]
+            m1 = np.where(ids == 1)[0]
+            if len(m0) and len(m1):
+                mid = 0.5 * (pts[m0[0]] + pts[m1[0]])
+                keep = [k for k in range(len(ids)) if k not in (m0[0], m1[0])]
+                pts = np.vstack([pts[keep], mid[None, :]])
+                ids = np.append(ids[keep], -1)  # -1 = merged, both truths
+        frames_pts.append(pts)
+        frames_ids.append(ids)
+    return frames_pts, frames_ids
+
+
+def score_tracks(tracks, frames_ids, frames_pts, tol=1.0):
+    """tracks: list of dicts with 'frames' (time idx) and 'pos' (xyz).
+
+    Coverage counts UNIQUE truth points (two tracks may honestly cover one
+    shared detection). Returns (n_tracks, switches, missing, pure)."""
+    truth_total = sum(int((ids >= 0).sum()) for ids in frames_ids)
+    # +2 merged points per occlusion frame counted separately
+    truth_total += sum(int((ids == -1).sum()) * 2 for ids in frames_ids)
+    switches, pure = 0, 0
+    covered: dict[tuple[int, int], set[int]] = {}  # (frame,row) -> track idx
+    for ti, t in enumerate(tracks):
+        seen: set[int] = set()
+        prev = None
+        for f, p in zip(t["frames"], t["pos"]):
+            ids = frames_ids[f]
+            pts = frames_pts[f]
+            if len(ids) == 0:
+                continue
+            d = np.linalg.norm(pts - p, axis=1)
+            b = int(np.argmin(d))
+            if d[b] > tol:
+                continue
+            covered.setdefault((int(f), int(b)), set()).add(ti)
+            if ids[b] == -1:
+                seen |= {0, 1}
+                truth = {0, 1}
+            else:
+                seen.add(int(ids[b]))
+                truth = {int(ids[b])}
+            if prev is not None and not truth & prev:
+                switches += 1
+            prev = truth if prev is None else (prev | truth)
+        if seen <= {0, 1} or len(seen) == 1:
+            pure += 1
+    # unique truth points covered: a merged row counts double only when two
+    # tracks honestly hold it (shared observation), else single.
+    uniq = 0
+    for (f, r), holders in covered.items():
+        if r == -1 and len(holders) >= 2:
+            uniq += 2
+        else:
+            uniq += 1
+    return {
+        "n_tracks": len(tracks),
+        "switches": switches,
+        "missing": max(0, truth_total - uniq),
+        "pure": pure,
+    }
+
+
+def to_links_format(tracks, key_frames="frames", key_pos="pos"):
+    return tracks
+
+
+def fast3d_tracks(frame_particles, v_max=3.0, share_tol=0.0):
+    """Drive track3d_loop_fast directly (mirrors Cython3DTracker.track_frames)
+    with extended per-frame arrays. Shared pairs materialize as virtual
+    carrier particles: shared coordinates, own history (prev -> frame-t
+    detection). Returns (tracks, n_shared)."""
+    from openptv2.algorithms.constants import NEXT_NONE, PREV_NONE
+    from openptv2.algorithms.track_kernels_track3d import track3d_loop_fast
+
+    nf = len(frame_particles)
+    ext_pos = [np.ascontiguousarray(p, dtype=np.float64)
+               for p in frame_particles]
+    ext_prev = [np.full(len(p), PREV_NONE, dtype=np.int32)
+                for p in frame_particles]
+    ext_next = [np.full(len(p), NEXT_NONE, dtype=np.int32)
+                for p in frame_particles]
+    n_shared_total = 0
+    for t in range(nf - 1):
+        if t == 0:
+            n0 = 0
+            pos_0 = np.empty((0, 3), dtype=np.float64)
+            prev_0 = np.empty(0, dtype=np.int32)
+        else:
+            n0 = len(ext_pos[t - 1])
+            pos_0, prev_0 = ext_pos[t - 1], ext_prev[t - 1]
+        n1, n2 = len(ext_pos[t]), len(ext_pos[t + 1])
+        if n1 == 0 or n2 == 0:
+            continue
+        pos_1, prev_1 = ext_pos[t], ext_prev[t]
+        nxt_1 = ext_next[t]
+        pos_2 = ext_pos[t + 1]
+        prev_2 = ext_prev[t + 1]
+        nxt_2 = ext_next[t + 1]
+        cap = max(n2, 1)
+        sc = np.zeros(1, dtype=np.int32)
+        si = np.full(cap, -1, dtype=np.int32)
+        sk = np.full(cap, -1, dtype=np.int32)
+        track3d_loop_fast(
+            n1, pos_0, prev_0, n0, pos_1, prev_1, nxt_1, n1,
+            pos_2, prev_2, nxt_2, n2,
+            v_max, v_max, v_max, 32, 0.0,
+            share_tol, sc, si, sk,
+        )
+        ext_next[t] = nxt_1
+        ext_prev[t + 1] = prev_2
+        # materialize virtual carriers for the NEXT steps (curr/prev roles)
+        for e in range(int(sc[0])):
+            i, k = int(si[e]), int(sk[e])
+            ext_pos[t + 1] = np.vstack(
+                [ext_pos[t + 1], pos_2[k][None, :]])
+            ext_prev[t + 1] = np.append(ext_prev[t + 1], np.int32(i))
+            ext_next[t + 1] = np.append(ext_next[t + 1], np.int32(NEXT_NONE))
+            n_shared_total += 1
+    # assemble (shared carriers emit their coords into their track)
+    counts = [len(p) for p in ext_pos]
+    visited = [np.zeros(n, dtype=bool) for n in counts]
+    tracks = []
+    for t in range(nf):
+        for i in range(counts[t]):
+            if visited[t][i]:
+                continue
+            tr_pos, tr_time, ct, ci = [], [], t, i
+            while ct < nf and ci != NEXT_NONE and not visited[ct][ci]:
+                visited[ct][ci] = True
+                tr_pos.append(ext_pos[ct][ci])
+                tr_time.append(ct)
+                nx = int(ext_next[ct][ci])
+                if nx >= 0 and ct + 1 < nf:
+                    ct, ci = ct + 1, nx
+                else:
+                    break
+            if len(tr_pos) >= 2:
+                tracks.append({"frames": tr_time,
+                               "pos": np.array(tr_pos)})
+    return tracks, n_shared_total
+
+
+def trackcorr_linkage_tracks(frame_particles, base_links, tol=1.0):
+    """Prototype trackcorr path: base_links (greedy forward output) -> plain
+    linkage arrays -> mark_shared_observations -> assemble_with_shared.
+    Returns (tracks_plain, tracks_shared, n_marks)."""
+    from openptv2.tracking_shared import (
+        mark_shared_observations,
+        assemble_with_shared,
+    )
+    from openptv2.algorithms.constants import PREV_NONE, NEXT_NONE
+
+    nf = len(frame_particles)
+    frames = {}
+    for f in range(nf):
+        n = len(frame_particles[f])
+        frames[f] = (np.full(n, PREV_NONE, dtype=np.int32),
+                     np.full(n, NEXT_NONE, dtype=np.int32),
+                     np.asarray(frame_particles[f], dtype=np.float64))
+    # row bookkeeping: links reference (frame, row); rows == indices here
+    for t0, r0, t1, r1 in base_links:
+        _, nxt, _ = frames[t0]
+        prv, _, _ = frames[t1]
+        if r0 < len(nxt) and r1 < len(prv):
+            nxt[r0] = r1
+            prv[r1] = r0
+    plain = assemble_with_shared(frames, 0, nf - 1, {})
+    shared = mark_shared_observations(frames, 0, nf - 1, tol=tol)
+    healed = assemble_with_shared(frames, 0, nf - 1, shared)
+    return plain, healed, len(shared)
+
+
+def greedy_links(frame_particles, gate=4.0):
+    """Baseline forward pass: nearest-neighbour chains (what trackcorr's
+    kernel would produce before postprocess). Returns link tuples."""
+    links = []
+    for t in range(len(frame_particles) - 1):
+        p0 = np.asarray(frame_particles[t])
+        p1 = np.asarray(frame_particles[t + 1])
+        used = set()
+        for i in range(len(p0)):
+            best, bj = gate, -1
+            for j in range(len(p1)):
+                if j in used:
+                    continue
+                d = float(np.linalg.norm(p1[j] - p0[i]))
+                if d < best:
+                    best, bj = d, j
+            if bj >= 0:
+                used.add(bj)
+                links.append((t, i, t + 1, bj))
+    return links
+
+
+if __name__ == "__main__":
+    from openptv2.plugins.two_phase_tracking import (
+        TwoPhaseTracker,
+        TwoPhaseTrackerConfig,
+    )
+
+    fp, fi = make_scene()
+    print("dets/frame:", [len(p) for p in fp], "truth(server): 6 trajs")
+    for share in [False, True]:
+        cfg = TwoPhaseTrackerConfig(v_max=4.0, max_gap=2, dt=1.0,
+                                    leaf_weight=0.0, cost_mode="3d",
+                                    use_velocity=True, allow_shared=share,
+                                    max_shared=3)
+        tr = TwoPhaseTracker(cfg)
+        links, chains = tr.track_frames([np.asarray(p) for p in fp],
+                                        return_chains=True)
+        tracks = [{"frames": c["frames"], "pos": c["pos"]} for c in chains
+                  if len(c["frames"]) >= 1]
+        n_shared = sum(sum(c["shared"]) for c in chains)
+        print(f"allow_shared={share}: links={len(links)} shared_pts={n_shared}",
+              score_tracks(tracks, fi, fp))
+
+    print("--- Fast3D (track3d_loop_fast, virtual carriers) ---")
+    for tol in [0.0, 0.3, 0.5, 1.0]:
+        tr, ns = fast3d_tracks([np.asarray(p) for p in fp], v_max=3.0,
+                               share_tol=tol)
+        print(f"share_tol={tol}: shared={ns}",
+              score_tracks(tr, fi, fp))
+
+    print("--- trackcorr-linkage (greedy links -> mark+assemble) ---")
+    base = greedy_links([np.asarray(p) for p in fp])
+    plain, healed, nm = trackcorr_linkage_tracks(
+        [np.asarray(p) for p in fp], base, tol=1.0)
+    print(f"plain: {score_tracks(plain, fi, fp)}")
+    print(f"shared marks={nm}: {score_tracks(healed, fi, fp)}")
diff --git a/setup.py b/setup.py
index d2e460e4..729f0f7a 100644
--- a/setup.py
+++ b/setup.py
@@ -58,7 +58,6 @@
     "track_kernels_position",  # 3D ray tracing + position reconstruction
     "track_kernels_track3d",  # stereo-3D tracking loop
     "track_kernels_corr",  # forward/backward tracking loops + CAS
-    "track_kernels_tracking",  # compatibility shim (re-exports above)
     "track_kernels_batch",  # batch processing + detection
 ]
 
diff --git a/src/openptv2/algorithms/track_kernels.py b/src/openptv2/algorithms/track_kernels.py
index 7862e433..ee29c7ec 100644
--- a/src/openptv2/algorithms/track_kernels.py
+++ b/src/openptv2/algorithms/track_kernels.py
@@ -105,6 +105,10 @@ def pack_mmlut(cal):
     init_mmlut_data_nlay_fast,
     targ_rec_fast,
 )
+from .track_kernels_corr import (  # noqa: E402, F401
+    trackback_loop_fast,
+    trackcorr_loop_fast,
+)
 from .track_kernels_geom import (  # noqa: E402, F401
     point_to_pixel_fast,
     searchquader_fast,
@@ -112,15 +116,9 @@ def pack_mmlut(cal):
 from .track_kernels_search import (  # noqa: E402, F401
     candsearch_in_pix_fast,
     candsearch_in_pix_rest_fast,
-    sort_candidates_by_freq_fast,
     sorted_candidates_fast,
 )
-from .track_kernels_tracking import (  # noqa: E402, F401
+from .track_kernels_track3d import (  # noqa: E402, F401
     track3d_loop_fast,
     track4be_loop_fast,
-    trackback_loop_fast,
-    trackcorr_loop_fast,
-)
-from .track_kernels_transform import (  # noqa: E402, F401
-    point_position_fast,
 )
diff --git a/src/openptv2/algorithms/track_kernels_batch.py b/src/openptv2/algorithms/track_kernels_batch.py
index 899910f8..7ebec371 100644
--- a/src/openptv2/algorithms/track_kernels_batch.py
+++ b/src/openptv2/algorithms/track_kernels_batch.py
@@ -37,145 +37,14 @@
 _M_PI: cython.double = 3.141592653589793
 
 
-from .track_kernels_geom import (
-    _multimed_r_nlay_1layer,
-    _ray_tracing_out,
-)
-from .track_kernels_transform import (
-    _metric_to_pixel_out,
-    _pixel_to_metric_out,
-    point_position_fast,
-)
-
-
-def ray_tracing_batch_fast(xy: cython.double[:, :], cal: cython.double[:]):
-    """Trace N rays through multi-media interface.
-
-    Args:
-        xy: (N, 2) float64 — metric image coordinates.
-        cal: (31,) float64 — packed calibration.
-
-    Returns:
-        (positions, directions) each (N, 3) float64.
-    """
-    n: cython.Py_ssize_t
-    i: cython.Py_ssize_t
-    _ray_out = np.empty(6, dtype=np.float64)
-    _ray_out_mv: cython.double[:] = _ray_out
-    n = xy.shape[0]
-    positions = np.empty((n, 3), dtype=np.float64)
-    directions = np.empty((n, 3), dtype=np.float64)
-    for i in range(n):
-        _ray_tracing_out(xy[i, 0], xy[i, 1], cal, _ray_out_mv)
-        positions[i, 0] = _ray_out_mv[0]
-        positions[i, 1] = _ray_out_mv[1]
-        positions[i, 2] = _ray_out_mv[2]
-        directions[i, 0] = _ray_out_mv[3]
-        directions[i, 1] = _ray_out_mv[4]
-        directions[i, 2] = _ray_out_mv[5]
-    return positions, directions
-
-
-@cython.ccall
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def point_position_batch_fast(
-    all_targets: cython.double[:, :, :],
-    num_pts: cython.int,
-    num_cams: cython.int,
-    cal_arrays,
-):
-    """Triangulate M targets from N cameras.
-
-    Args:
-        all_targets: (M, num_cams, 2) float64.
-        num_pts: M.
-        num_cams: N.
-        cal_arrays: tuple of (31,) float64 arrays.
-
-    Returns:
-        (positions, distances) — (M, 3) and (M,) float64.
-    """
-    i: cython.Py_ssize_t
-    dist: cython.double
-    positions = np.empty((num_pts, 3), dtype=np.float64)
-    distances = np.empty(num_pts, dtype=np.float64)
-    for i in range(num_pts):
-        _cal_arr = np.asarray(list(cal_arrays), dtype=np.float64)
-        pos, dist = point_position_fast(all_targets[i], num_cams, _cal_arr)
-        positions[i, 0] = pos[0]
-        positions[i, 1] = pos[1]
-        positions[i, 2] = pos[2]
-        distances[i] = dist
-    return positions, distances
-
-
-@cython.ccall
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def pixel_to_metric_batch_fast(
-    xy: cython.double[:, :],
-    imx: cython.int,
-    imy: cython.int,
-    pix_x: cython.double,
-    pix_y: cython.double,
-    chfield: cython.int,
-):
-    """Convert N pixel coordinates to metric."""
-    n: cython.Py_ssize_t
-    i: cython.Py_ssize_t
-    _pp = np.empty(2, dtype=np.float64)
-    _pp_mv: cython.double[:] = _pp
-    n = xy.shape[0]
-    result = np.empty((n, 2), dtype=np.float64)
-    for i in range(n):
-        _pixel_to_metric_out(
-            xy[i, 0],
-            xy[i, 1],
-            imx,
-            imy,
-            pix_x,
-            pix_y,
-            chfield,
-            _pp_mv,
-        )
-        result[i, 0] = _pp_mv[0]
-        result[i, 1] = _pp_mv[1]
-    return result
-
-
-@cython.ccall
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def metric_to_pixel_batch_fast(
-    xy: cython.double[:, :],
-    imx: cython.int,
-    imy: cython.int,
-    pix_x: cython.double,
-    pix_y: cython.double,
-    chfield: cython.int,
-):
-    """Convert N metric coordinates to pixel."""
-    n: cython.Py_ssize_t
-    i: cython.Py_ssize_t
-    _pp = np.empty(2, dtype=np.float64)
-    _pp_mv: cython.double[:] = _pp
-    n = xy.shape[0]
-    result = np.empty((n, 2), dtype=np.float64)
-    for i in range(n):
-        _metric_to_pixel_out(
-            xy[i, 0],
-            xy[i, 1],
-            imx,
-            imy,
-            pix_x,
-            pix_y,
-            chfield,
-            _pp_mv,
-        )
-        result[i, 0] = _pp_mv[0]
-        result[i, 1] = _pp_mv[1]
-    return result
+if cython.compiled:
+    from cython.cimports.openptv2.algorithms.track_kernels_pixel import (
+        _multimed_r_nlay_1layer,
+    )
+else:
+    from .track_kernels_pixel import (
+        _multimed_r_nlay_1layer,
+    )
 
 
 @cython.ccall
diff --git a/src/openptv2/algorithms/track_kernels_corr.py b/src/openptv2/algorithms/track_kernels_corr.py
index 94ea23d7..668ef679 100644
--- a/src/openptv2/algorithms/track_kernels_corr.py
+++ b/src/openptv2/algorithms/track_kernels_corr.py
@@ -4,16 +4,10 @@
 import numpy as np
 
 if cython.compiled:
-    from cython.cimports.libc.math import (
-        acos as c_acos,
-    )
     from cython.cimports.libc.math import (
         sqrt as c_sqrt,
     )
 else:
-    from math import (
-        acos as c_acos,
-    )
     from math import (
         sqrt as c_sqrt,
     )
@@ -24,6 +18,7 @@
         _sorted_candidates_fast_out_nogil,
     )
     from cython.cimports.openptv2.algorithms.track_kernels_position import (
+        _angle_acc_out,
         _point_position_out,
         assess_new_position_fast_nogil,
     )
@@ -33,6 +28,7 @@
         _sorted_candidates_fast_out_nogil,
     )
     from .track_kernels_position import (
+        _angle_acc_out,
         _point_position_out,
         assess_new_position_fast_nogil,
     )
@@ -61,73 +57,6 @@
 ADD_PART_K = 3.0
 
 
-@cython.ccall
-@cython.inline
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.cdivision(True)
-@cython.profile(False)
-@cython.nogil
-def _angle_acc_out(
-    start_x: cython.double,
-    start_y: cython.double,
-    start_z: cython.double,
-    pred_x: cython.double,
-    pred_y: cython.double,
-    pred_z: cython.double,
-    cand_x: cython.double,
-    cand_y: cython.double,
-    cand_z: cython.double,
-    out: cython.double[:],
-) -> cython.int:
-    """Write angle and acc to out[0], out[1] — no tuple creation."""
-    v0x: cython.double
-    v0y: cython.double
-    v0z: cython.double
-    v1x: cython.double
-    v1y: cython.double
-    v1z: cython.double
-    angle: cython.double
-    norm0: cython.double
-    norm1: cython.double
-    dot: cython.double
-    dx: cython.double
-    dy: cython.double
-    dz: cython.double
-    acc: cython.double
-    v0x = pred_x - start_x
-    v0y = pred_y - start_y
-    v0z = pred_z - start_z
-    v1x = cand_x - start_x
-    v1y = cand_y - start_y
-    v1z = cand_z - start_z
-
-    if v0x == -v1x and v0y == -v1y and v0z == -v1z:
-        angle = 200.0
-    elif v0x == v1x and v0y == v1y and v0z == v1z:
-        angle = 0.0
-    else:
-        norm0 = c_sqrt(v0x * v0x + v0y * v0y + v0z * v0z)
-        norm1 = c_sqrt(v1x * v1x + v1y * v1y + v1z * v1z)
-        if norm0 == 0.0 or norm1 == 0.0:
-            angle = 0.0
-        else:
-            dot = (v0x * v1x + v0y * v1y + v0z * v1z) / (norm0 * norm1)
-            if dot > 1.0:
-                dot = 1.0
-            elif dot < -1.0:
-                dot = -1.0
-            angle = c_acos(dot) * 200.0 / 3.141592653589793
-
-    dx = v1x - v0x
-    dy = v1y - v0y
-    dz = v1z - v0z
-    acc = c_sqrt(dx * dx + dy * dy + dz * dz)
-    out[0] = angle
-    out[1] = acc
-    return 0
-
-
 @cython.cfunc
 @cython.inline
 @cython.boundscheck(False)
diff --git a/src/openptv2/algorithms/track_kernels_geom.py b/src/openptv2/algorithms/track_kernels_geom.py
index e1f8a82e..487047db 100644
--- a/src/openptv2/algorithms/track_kernels_geom.py
+++ b/src/openptv2/algorithms/track_kernels_geom.py
@@ -8,15 +8,6 @@
 import numpy as np
 
 if cython.compiled:
-    from cython.cimports.libc.math import (
-        acos as c_acos,
-    )
-    from cython.cimports.libc.math import (
-        asin as c_asin,
-    )
-    from cython.cimports.libc.math import (
-        atan as c_atan,
-    )
     from cython.cimports.libc.math import (
         cos as c_cos,
     )
@@ -26,19 +17,7 @@
     from cython.cimports.libc.math import (
         sqrt as c_sqrt,
     )
-    from cython.cimports.libc.math import (
-        tan as c_tan,
-    )
 else:
-    from math import (
-        acos as c_acos,
-    )
-    from math import (
-        asin as c_asin,
-    )
-    from math import (
-        atan as c_atan,
-    )
     from math import (
         cos as c_cos,
     )
@@ -48,12 +27,18 @@
     from math import (
         sqrt as c_sqrt,
     )
-    from math import (
-        tan as c_tan,
-    )
 
 _M_PI: cython.double = 3.141592653589793
 
+if cython.compiled:
+    from cython.cimports.openptv2.algorithms.track_kernels_pixel import (
+        _multimed_r_nlay_1layer,
+    )
+else:
+    from .track_kernels_pixel import (
+        _multimed_r_nlay_1layer,
+    )
+
 
 # Cal array layout (31 float64):
 #  0-2:   ext_x0, ext_y0, ext_z0
@@ -69,85 +54,6 @@
 CAL_ARRAY_SIZE = 31
 
 
-@cython.ccall
-@cython.nogil
-def _multimed_r_nlay_1layer(
-    pos_x: cython.double,
-    pos_y: cython.double,
-    pos_z: cython.double,
-    ext_x0: cython.double,
-    ext_y0: cython.double,
-    ext_z0: cython.double,
-    mm_n1: cython.double,
-    mm_n2_0: cython.double,
-    mm_n3: cython.double,
-    mm_d0: cython.double,
-) -> cython.double:
-    """Single-layer iterative radial shift."""
-    zout: cython.double
-    dx: cython.double
-    dy: cython.double
-    r: cython.double
-    rq: cython.double
-    it: cython.int
-    denom: cython.double
-    beta1: cython.double
-    sin_beta1: cython.double
-    arg: cython.double
-    beta2_0: cython.double
-    arg3: cython.double
-    beta3: cython.double
-    rbeta: cython.double
-    rdiff: cython.double
-    if mm_n1 == 1.0 and mm_n2_0 == 1.0 and mm_n3 == 1.0:
-        return 1.0
-
-    zout = pos_z
-    dx = pos_x - ext_x0
-    dy = pos_y - ext_y0
-    r = c_sqrt(dx * dx + dy * dy)
-    rq = r
-
-    for it in range(40):
-        denom = ext_z0 - pos_z
-        if denom == 0.0:
-            return 1.0
-        beta1 = c_atan(rq / denom)
-        sin_beta1 = c_sin(beta1)
-
-        arg = sin_beta1 * mm_n1 / mm_n2_0
-        if arg > 1.0:
-            arg = 1.0
-        elif arg < -1.0:
-            arg = -1.0
-        beta2_0 = c_asin(arg)
-
-        arg3 = sin_beta1 * mm_n1 / mm_n3
-        if arg3 > 1.0:
-            arg3 = 1.0
-        elif arg3 < -1.0:
-            arg3 = -1.0
-        beta3 = c_asin(arg3)
-
-        rbeta = (
-            (ext_z0 - mm_d0) * c_tan(beta1)
-            + mm_d0 * c_tan(beta2_0)
-            - zout * c_tan(beta3)
-        )
-
-        rdiff = r - rbeta
-        rq += rdiff
-
-        if abs(rdiff) < 0.001:
-            break
-    else:
-        return 1.0
-
-    if r != 0.0:
-        return rq / r
-    return 1.0
-
-
 @cython.ccall
 def point_to_pixel_fast(
     pos: cython.double[:],
@@ -427,277 +333,6 @@ def point_to_pixel_fast(
     return x_pixel, y_pixel
 
 
-@cython.ccall
-@cython.inline
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.cdivision(True)
-@cython.profile(False)
-@cython.nogil
-def _point_to_pixel_out(
-    pos: cython.double[:],
-    cal: cython.double[:],
-    mmlut_data: cython.double[:],
-    mmlut_origin: cython.double[:],
-    mmlut_nr: cython.int,
-    mmlut_nz: cython.int,
-    mmlut_rw: cython.double,
-    has_mmlut: cython.int,
-    imx_half: cython.double,
-    imy_half: cython.double,
-    inv_pix_x: cython.double,
-    inv_pix_y: cython.double,
-    chfield: cython.int,
-    out: cython.double[:],
-) -> cython.int:
-    """Write pixel coordinates to out[0], out[1] — no tuple creation."""
-    pos0: cython.double
-    pos1: cython.double
-    pos2: cython.double
-    ext_x0: cython.double
-    ext_y0: cython.double
-    ext_z0: cython.double
-    dm00: cython.double
-    dm10: cython.double
-    dm20: cython.double
-    dm01: cython.double
-    dm11: cython.double
-    dm21: cython.double
-    dm02: cython.double
-    dm12: cython.double
-    dm22: cython.double
-    int_cc: cython.double
-    xh: cython.double
-    yh: cython.double
-    gx: cython.double
-    gy: cython.double
-    gz: cython.double
-    inv_dog: cython.double
-    mm_n1: cython.double
-    mm_n2_0: cython.double
-    mm_n3: cython.double
-    mm_d0: cython.double
-    k1: cython.double
-    k2: cython.double
-    k3: cython.double
-    p1: cython.double
-    p2: cython.double
-    scx: cython.double
-    she: cython.double
-    dot_cam: cython.double
-    dist_o_glas: cython.double
-    dist_cam_glas: cython.double
-    dot_pos: cython.double
-    dist_point_glas: cython.double
-    s_cam: cython.double
-    cc_x: cython.double
-    cc_y: cython.double
-    cc_z: cython.double
-    s_pt: cython.double
-    cp_x: cython.double
-    cp_y: cython.double
-    cp_z: cython.double
-    ext_t_z0: cython.double
-    s_d: cython.double
-    ag_x: cython.double
-    ag_y: cython.double
-    ag_z: cython.double
-    tmp_x: cython.double
-    tmp_y: cython.double
-    tmp_z: cython.double
-    pos_t_0: cython.double
-    pos_t_2: cython.double
-    radial_shift: cython.double
-    tx: cython.double
-    ty: cython.double
-    tz: cython.double
-    sz: cython.double
-    iz: cython.int
-    R: cython.double
-    sr: cython.double
-    ir: cython.int
-    v0: cython.int
-    v3: cython.int
-    mmf: cython.double
-    X_t: cython.double
-    s_z: cython.double
-    bx: cython.double
-    by: cython.double
-    bz: cython.double
-    s_x: cython.double
-    dx: cython.double
-    dy: cython.double
-    dz: cython.double
-    deno: cython.double
-    x: cython.double
-    y: cython.double
-    r: cython.double
-    r2: cython.double
-    r4: cython.double
-    radial_factor: cython.double
-    xd: cython.double
-    yd: cython.double
-    sin_she: cython.double
-    cos_she: cython.double
-    x_dist: cython.double
-    y_dist: cython.double
-    x_pixel: cython.double
-    y_pixel: cython.double
-    pos0 = pos[0]
-    pos1 = pos[1]
-    pos2 = pos[2]
-
-    ext_x0 = cal[0]
-    ext_y0 = cal[1]
-    ext_z0 = cal[2]
-    dm00 = cal[3]
-    dm10 = cal[4]
-    dm20 = cal[5]
-    dm01 = cal[6]
-    dm11 = cal[7]
-    dm21 = cal[8]
-    dm02 = cal[9]
-    dm12 = cal[10]
-    dm22 = cal[11]
-    int_cc = cal[12]
-    xh = cal[13]
-    yh = cal[14]
-    gx = cal[15]
-    gy = cal[16]
-    gz = cal[17]
-    inv_dog = cal[19]
-    mm_n1 = cal[20]
-    mm_n2_0 = cal[21]
-    mm_n3 = cal[22]
-    mm_d0 = cal[23]
-    k1 = cal[24]
-    k2 = cal[25]
-    k3 = cal[26]
-    p1 = cal[27]
-    p2 = cal[28]
-    scx = cal[29]
-    she = cal[30]
-
-    # trans_cam_point
-    dot_cam = ext_x0 * gx + ext_y0 * gy + ext_z0 * gz
-    dist_o_glas = cal[18]
-    dist_cam_glas = dot_cam * inv_dog - dist_o_glas - mm_d0
-
-    dot_pos = pos0 * gx + pos1 * gy + pos2 * gz
-    dist_point_glas = dot_pos * inv_dog - dist_o_glas
-
-    s_cam = dist_cam_glas * inv_dog
-    cc_x = ext_x0 - gx * s_cam
-    cc_y = ext_y0 - gy * s_cam
-    cc_z = ext_z0 - gz * s_cam
-
-    s_pt = dist_point_glas * inv_dog
-    cp_x = pos0 - gx * s_pt
-    cp_y = pos1 - gy * s_pt
-    cp_z = pos2 - gz * s_pt
-
-    ext_t_z0 = dist_cam_glas + mm_d0
-
-    s_d = mm_d0 * inv_dog
-    ag_x = cc_x - gx * s_d
-    ag_y = cc_y - gy * s_d
-    ag_z = cc_z - gz * s_d
-    tmp_x = cp_x - ag_x
-    tmp_y = cp_y - ag_y
-    tmp_z = cp_z - ag_z
-
-    pos_t_0 = c_sqrt(tmp_x * tmp_x + tmp_y * tmp_y + tmp_z * tmp_z)
-    pos_t_2 = dist_point_glas
-
-    # mmlut lookup + multimed_nlay
-    radial_shift = 1.0
-    if has_mmlut:
-        tx = pos_t_0 - mmlut_origin[0]
-        ty = -mmlut_origin[1]
-        tz = pos_t_2 - mmlut_origin[2]
-        sz = tz / mmlut_rw
-        iz = int(sz)
-        sz -= iz
-        R = c_sqrt(tx * tx + ty * ty)
-        sr = R / mmlut_rw
-        ir = int(sr)
-        sr -= ir
-        if ir <= mmlut_nr and iz >= 0 and iz <= mmlut_nz:
-            v0 = ir * mmlut_nz + iz
-            v3 = v0 + mmlut_nz + 1
-            if v0 >= 0 and v3 <= mmlut_nr * mmlut_nz:
-                mmf = (
-                    mmlut_data[v0] * (1.0 - sr) * (1.0 - sz)
-                    + mmlut_data[v0 + 1] * (1.0 - sr) * sz
-                    + mmlut_data[v0 + mmlut_nz] * sr * (1.0 - sz)
-                    + mmlut_data[v3] * sr * sz
-                )
-                if mmf > 0.0:
-                    radial_shift = mmf
-    if radial_shift == 1.0:
-        radial_shift = _multimed_r_nlay_1layer(
-            pos_t_0,
-            0.0,
-            pos_t_2,
-            0.0,
-            0.0,
-            ext_t_z0,
-            mm_n1,
-            mm_n2_0,
-            mm_n3,
-            mm_d0,
-        )
-    X_t = pos_t_0 * radial_shift
-
-    # back_trans_point
-    s_z = -pos_t_2 * inv_dog
-    bx = ag_x - gx * s_z
-    by = ag_y - gy * s_z
-    bz = ag_z - gz * s_z
-    if pos_t_0 > 0.0:
-        s_x = -X_t / pos_t_0
-        bx -= tmp_x * s_x
-        by -= tmp_y * s_x
-        bz -= tmp_z * s_x
-
-    # perspective projection
-    dx = bx - ext_x0
-    dy = by - ext_y0
-    dz = bz - ext_z0
-    deno = dm02 * dx + dm12 * dy + dm22 * dz
-    x = -int_cc * (dm00 * dx + dm10 * dy + dm20 * dz) / deno
-    y = -int_cc * (dm01 * dx + dm11 * dy + dm21 * dz) / deno
-
-    # flat_to_dist + distort_brown_affin
-    x += xh
-    y += yh
-    r = c_sqrt(x * x + y * y)
-    if r < 1e-10:
-        x_dist = 0.0
-        y_dist = 0.0
-    else:
-        r2 = r * r
-        r4 = r2 * r2
-        radial_factor = 1.0 + k1 * r2 + k2 * r4 + k3 * r4 * r2
-        xd = x * radial_factor + p1 * (r2 + 2.0 * x * x) + 2.0 * p2 * x * y
-        yd = y * radial_factor + p2 * (r2 + 2.0 * y * y) + 2.0 * p1 * x * y
-        sin_she = c_sin(she)
-        cos_she = c_cos(she)
-        x_dist = scx * (xd - sin_she * yd)
-        y_dist = scx * cos_she * yd
-
-    # metric_to_pixel
-    x_pixel = x_dist * inv_pix_x + imx_half
-    y_pixel = imy_half - y_dist * inv_pix_y
-    if chfield == 1:
-        y_pixel = (y_pixel - 1.0) * 0.5
-    elif chfield == 2:
-        y_pixel = y_pixel * 0.5
-    out[0] = x_pixel
-    out[1] = y_pixel
-    return 0
-
-
 PT_UNUSED = -999
 
 
@@ -836,492 +471,3 @@ def searchquader_fast(
 
     return xr, xl, yd, yu
 
-
-@cython.ccall
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def angle_acc_fast(
-    start_x: cython.double,
-    start_y: cython.double,
-    start_z: cython.double,
-    pred_x: cython.double,
-    pred_y: cython.double,
-    pred_z: cython.double,
-    cand_x: cython.double,
-    cand_y: cython.double,
-    cand_z: cython.double,
-):
-    """Compute angle and acceleration between predicted and candidate."""
-    v0x: cython.double
-    v0y: cython.double
-    v0z: cython.double
-    v1x: cython.double
-    v1y: cython.double
-    v1z: cython.double
-    angle: cython.double
-    norm0: cython.double
-    norm1: cython.double
-    dot: cython.double
-    dx: cython.double
-    dy: cython.double
-    dz: cython.double
-    acc: cython.double
-    v0x = pred_x - start_x
-    v0y = pred_y - start_y
-    v0z = pred_z - start_z
-    v1x = cand_x - start_x
-    v1y = cand_y - start_y
-    v1z = cand_z - start_z
-
-    if v0x == -v1x and v0y == -v1y and v0z == -v1z:
-        angle = 200.0
-    elif v0x == v1x and v0y == v1y and v0z == v1z:
-        angle = 0.0
-    else:
-        norm0 = c_sqrt(v0x * v0x + v0y * v0y + v0z * v0z)
-        norm1 = c_sqrt(v1x * v1x + v1y * v1y + v1z * v1z)
-        if norm0 == 0.0 or norm1 == 0.0:
-            angle = 0.0
-        else:
-            dot = (v0x * v1x + v0y * v1y + v0z * v1z) / (norm0 * norm1)
-            if dot > 1.0:
-                dot = 1.0
-            elif dot < -1.0:
-                dot = -1.0
-            angle = c_acos(dot) * 200.0 / _M_PI
-
-    dx = v1x - v0x
-    dy = v1y - v0y
-    dz = v1z - v0z
-    acc = c_sqrt(dx * dx + dy * dy + dz * dz)
-    return angle, acc
-
-
-@cython.ccall
-@cython.inline
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.cdivision(True)
-@cython.profile(False)
-@cython.nogil
-def _angle_acc_out(
-    start_x: cython.double,
-    start_y: cython.double,
-    start_z: cython.double,
-    pred_x: cython.double,
-    pred_y: cython.double,
-    pred_z: cython.double,
-    cand_x: cython.double,
-    cand_y: cython.double,
-    cand_z: cython.double,
-    out: cython.double[:],
-) -> cython.int:
-    """Write angle and acc to out[0], out[1] — no tuple creation."""
-    v0x: cython.double
-    v0y: cython.double
-    v0z: cython.double
-    v1x: cython.double
-    v1y: cython.double
-    v1z: cython.double
-    angle: cython.double
-    norm0: cython.double
-    norm1: cython.double
-    dot: cython.double
-    dx: cython.double
-    dy: cython.double
-    dz: cython.double
-    acc: cython.double
-    v0x = pred_x - start_x
-    v0y = pred_y - start_y
-    v0z = pred_z - start_z
-    v1x = cand_x - start_x
-    v1y = cand_y - start_y
-    v1z = cand_z - start_z
-
-    if v0x == -v1x and v0y == -v1y and v0z == -v1z:
-        angle = 200.0
-    elif v0x == v1x and v0y == v1y and v0z == v1z:
-        angle = 0.0
-    else:
-        norm0 = c_sqrt(v0x * v0x + v0y * v0y + v0z * v0z)
-        norm1 = c_sqrt(v1x * v1x + v1y * v1y + v1z * v1z)
-        if norm0 == 0.0 or norm1 == 0.0:
-            angle = 0.0
-        else:
-            dot = (v0x * v1x + v0y * v1y + v0z * v1z) / (norm0 * norm1)
-            if dot > 1.0:
-                dot = 1.0
-            elif dot < -1.0:
-                dot = -1.0
-            angle = c_acos(dot) * 200.0 / 3.141592653589793
-
-    dx = v1x - v0x
-    dy = v1y - v0y
-    dz = v1z - v0z
-    acc = c_sqrt(dx * dx + dy * dy + dz * dz)
-    out[0] = angle
-    out[1] = acc
-    return 0
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def _ray_tracing_fast(x: cython.double, y: cython.double, cal: cython.double[:]):
-    """Trace ray through multi-media interface.
-
-    Returns (Xx, Xy, Xz, ox, oy, oz) — crossing point and direction.
-    """
-    ext_x0: cython.double
-    ext_y0: cython.double
-    ext_z0: cython.double
-    dm00: cython.double
-    dm10: cython.double
-    dm20: cython.double
-    dm01: cython.double
-    dm11: cython.double
-    dm21: cython.double
-    dm02: cython.double
-    dm12: cython.double
-    dm22: cython.double
-    int_cc: cython.double
-    gx: cython.double
-    gy: cython.double
-    gz: cython.double
-    mm_n1: cython.double
-    mm_n2_0: cython.double
-    mm_n3: cython.double
-    mm_d0: cython.double
-    t0: cython.double
-    t1: cython.double
-    t2: cython.double
-    tn: cython.double
-    sd0: cython.double
-    sd1: cython.double
-    sd2: cython.double
-    gn: cython.double
-    gd0: cython.double
-    gd1: cython.double
-    gd2: cython.double
-    c: cython.double
-    dcg: cython.double
-    denom: cython.double
-    d1: cython.double
-    Xb0: cython.double
-    Xb1: cython.double
-    Xb2: cython.double
-    n: cython.double
-    bp0: cython.double
-    bp1: cython.double
-    bp2: cython.double
-    bpn: cython.double
-    p: cython.double
-    n_glass: cython.double
-    a2_0: cython.double
-    a2_1: cython.double
-    a2_2: cython.double
-    d2_denom: cython.double
-    d2: cython.double
-    Xx: cython.double
-    Xy: cython.double
-    Xz: cython.double
-    n_a2: cython.double
-    p2: cython.double
-    n_final: cython.double
-    ox: cython.double
-    oy: cython.double
-    oz: cython.double
-    ext_x0 = cal[0]
-    ext_y0 = cal[1]
-    ext_z0 = cal[2]
-    dm00 = cal[3]
-    dm10 = cal[4]
-    dm20 = cal[5]
-    dm01 = cal[6]
-    dm11 = cal[7]
-    dm21 = cal[8]
-    dm02 = cal[9]
-    dm12 = cal[10]
-    dm22 = cal[11]
-    int_cc = cal[12]
-    gx = cal[15]
-    gy = cal[16]
-    gz = cal[17]
-    mm_n1 = cal[20]
-    mm_n2_0 = cal[21]
-    mm_n3 = cal[22]
-    mm_d0 = cal[23]
-
-    # tmp1 = unit_vector([x, y, -int_cc])
-    t0 = x
-    t1 = y
-    t2 = -int_cc
-    tn = c_sqrt(t0 * t0 + t1 * t1 + t2 * t2)
-    if tn > 0.0:
-        t0 /= tn
-        t1 /= tn
-        t2 /= tn
-
-    # start_dir = dm @ tmp1
-    sd0 = dm00 * t0 + dm01 * t1 + dm02 * t2
-    sd1 = dm10 * t0 + dm11 * t1 + dm12 * t2
-    sd2 = dm20 * t0 + dm21 * t1 + dm22 * t2
-
-    # glass_dir = unit_vector(glass_vec)
-    gn = c_sqrt(gx * gx + gy * gy + gz * gz)
-    if gn > 0.0:
-        gd0 = gx / gn
-        gd1 = gy / gn
-        gd2 = gz / gn
-    else:
-        gd0 = 0.0
-        gd1 = 0.0
-        gd2 = 0.0
-    c = gn + mm_d0
-
-    # dist_cam_glass, d1
-    dcg = gd0 * ext_x0 + gd1 * ext_y0 + gd2 * ext_z0 - c
-    denom = gd0 * sd0 + gd1 * sd1 + gd2 * sd2
-    d1 = -dcg / denom
-
-    # Xb = primary_point + start_dir * d1
-    Xb0 = ext_x0 + sd0 * d1
-    Xb1 = ext_y0 + sd1 * d1
-    Xb2 = ext_z0 + sd2 * d1
-
-    # Decompose ray: n = dot(start_dir, glass_dir)
-    n = sd0 * gd0 + sd1 * gd1 + sd2 * gd2
-    # bp = unit_vector(start_dir - glass_dir * n)
-    bp0 = sd0 - gd0 * n
-    bp1 = sd1 - gd1 * n
-    bp2 = sd2 - gd2 * n
-    bpn = c_sqrt(bp0 * bp0 + bp1 * bp1 + bp2 * bp2)
-    if bpn > 0.0:
-        bp0 /= bpn
-        bp1 /= bpn
-        bp2 /= bpn
-
-    # Snell's law: air -> glass
-    p = c_sqrt(1.0 - n * n) * mm_n1 / mm_n2_0
-    n_glass = c_sqrt(1.0 - p * p) if n >= 0 else -c_sqrt(1.0 - p * p)
-
-    # a2 = bp * p + glass_dir * n_glass
-    a2_0 = bp0 * p + gd0 * n_glass
-    a2_1 = bp1 * p + gd1 * n_glass
-    a2_2 = bp2 * p + gd2 * n_glass
-
-    d2_denom = gd0 * a2_0 + gd1 * a2_1 + gd2 * a2_2
-    d2 = mm_d0 / abs(d2_denom)
-
-    # X = Xb + a2 * d2
-    Xx = Xb0 + a2_0 * d2
-    Xy = Xb1 + a2_1 * d2
-    Xz = Xb2 + a2_2 * d2
-
-    # Direction in next medium: Snell glass -> water
-    n_a2 = a2_0 * gd0 + a2_1 * gd1 + a2_2 * gd2
-    # bp = unit_vector(a2 - glass_dir * n_glass)
-    bp0 = a2_0 - gd0 * n_glass
-    bp1 = a2_1 - gd1 * n_glass
-    bp2 = a2_2 - gd2 * n_glass
-    bpn = c_sqrt(bp0 * bp0 + bp1 * bp1 + bp2 * bp2)
-    if bpn > 0.0:
-        bp0 /= bpn
-        bp1 /= bpn
-        bp2 /= bpn
-
-    p2 = c_sqrt(1.0 - n_a2 * n_a2) * mm_n2_0 / mm_n3
-    n_final = c_sqrt(1.0 - p2 * p2) if n_a2 >= 0 else -c_sqrt(1.0 - p2 * p2)
-
-    ox = bp0 * p2 + gd0 * n_final
-    oy = bp1 * p2 + gd1 * n_final
-    oz = bp2 * p2 + gd2 * n_final
-
-    return Xx, Xy, Xz, ox, oy, oz
-
-
-@cython.ccall
-@cython.inline
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.cdivision(True)
-@cython.profile(False)
-@cython.nogil
-def _ray_tracing_out(
-    x: cython.double,
-    y: cython.double,
-    cal: cython.double[:],
-    out: cython.double[:],
-) -> cython.int:
-    """Write ray tracing results into out[0:6] — no tuple creation."""
-    ext_x0: cython.double
-    ext_y0: cython.double
-    ext_z0: cython.double
-    dm00: cython.double
-    dm10: cython.double
-    dm20: cython.double
-    dm01: cython.double
-    dm11: cython.double
-    dm21: cython.double
-    dm02: cython.double
-    dm12: cython.double
-    dm22: cython.double
-    int_cc: cython.double
-    gx: cython.double
-    gy: cython.double
-    gz: cython.double
-    mm_n1: cython.double
-    mm_n2_0: cython.double
-    mm_n3: cython.double
-    mm_d0: cython.double
-    t0: cython.double
-    t1: cython.double
-    t2: cython.double
-    tn: cython.double
-    sd0: cython.double
-    sd1: cython.double
-    sd2: cython.double
-    gn: cython.double
-    gd0: cython.double
-    gd1: cython.double
-    gd2: cython.double
-    c: cython.double
-    dcg: cython.double
-    denom: cython.double
-    d1: cython.double
-    Xb0: cython.double
-    Xb1: cython.double
-    Xb2: cython.double
-    n: cython.double
-    bp0: cython.double
-    bp1: cython.double
-    bp2: cython.double
-    bpn: cython.double
-    p: cython.double
-    n_glass: cython.double
-    a2_0: cython.double
-    a2_1: cython.double
-    a2_2: cython.double
-    d2_denom: cython.double
-    d2: cython.double
-    Xx: cython.double
-    Xy: cython.double
-    Xz: cython.double
-    n_a2: cython.double
-    p2: cython.double
-    n_final: cython.double
-    ox: cython.double
-    oy: cython.double
-    oz: cython.double
-    ext_x0 = cal[0]
-    ext_y0 = cal[1]
-    ext_z0 = cal[2]
-    dm00 = cal[3]
-    dm10 = cal[4]
-    dm20 = cal[5]
-    dm01 = cal[6]
-    dm11 = cal[7]
-    dm21 = cal[8]
-    dm02 = cal[9]
-    dm12 = cal[10]
-    dm22 = cal[11]
-    int_cc = cal[12]
-    gx = cal[15]
-    gy = cal[16]
-    gz = cal[17]
-    mm_n1 = cal[20]
-    mm_n2_0 = cal[21]
-    mm_n3 = cal[22]
-    mm_d0 = cal[23]
-
-    # tmp1 = unit_vector([x, y, -int_cc])
-    t0 = x
-    t1 = y
-    t2 = -int_cc
-    tn = c_sqrt(t0 * t0 + t1 * t1 + t2 * t2)
-    if tn > 0.0:
-        t0 /= tn
-        t1 /= tn
-        t2 /= tn
-
-    # start_dir = dm @ tmp1
-    sd0 = dm00 * t0 + dm01 * t1 + dm02 * t2
-    sd1 = dm10 * t0 + dm11 * t1 + dm12 * t2
-    sd2 = dm20 * t0 + dm21 * t1 + dm22 * t2
-
-    # glass_dir = unit_vector(glass_vec)
-    gn = c_sqrt(gx * gx + gy * gy + gz * gz)
-    if gn > 0.0:
-        gd0 = gx / gn
-        gd1 = gy / gn
-        gd2 = gz / gn
-    else:
-        gd0 = 0.0
-        gd1 = 0.0
-        gd2 = 0.0
-    c = gn + mm_d0
-
-    # dist_cam_glass, d1
-    dcg = gd0 * ext_x0 + gd1 * ext_y0 + gd2 * ext_z0 - c
-    denom = gd0 * sd0 + gd1 * sd1 + gd2 * sd2
-    d1 = -dcg / denom
-
-    # Xb = primary_point + start_dir * d1
-    Xb0 = ext_x0 + sd0 * d1
-    Xb1 = ext_y0 + sd1 * d1
-    Xb2 = ext_z0 + sd2 * d1
-
-    # Decompose ray: n = dot(start_dir, glass_dir)
-    n = sd0 * gd0 + sd1 * gd1 + sd2 * gd2
-    # bp = unit_vector(start_dir - glass_dir * n)
-    bp0 = sd0 - gd0 * n
-    bp1 = sd1 - gd1 * n
-    bp2 = sd2 - gd2 * n
-    bpn = c_sqrt(bp0 * bp0 + bp1 * bp1 + bp2 * bp2)
-    if bpn > 0.0:
-        bp0 /= bpn
-        bp1 /= bpn
-        bp2 /= bpn
-
-    # Snell's law: air -> glass
-    p = c_sqrt(1.0 - n * n) * mm_n1 / mm_n2_0
-    n_glass = c_sqrt(1.0 - p * p) if n >= 0 else -c_sqrt(1.0 - p * p)
-
-    # a2 = bp * p + glass_dir * n_glass
-    a2_0 = bp0 * p + gd0 * n_glass
-    a2_1 = bp1 * p + gd1 * n_glass
-    a2_2 = bp2 * p + gd2 * n_glass
-
-    d2_denom = gd0 * a2_0 + gd1 * a2_1 + gd2 * a2_2
-    d2 = mm_d0 / abs(d2_denom)
-
-    # X = Xb + a2 * d2
-    Xx = Xb0 + a2_0 * d2
-    Xy = Xb1 + a2_1 * d2
-    Xz = Xb2 + a2_2 * d2
-
-    # Direction in next medium: Snell glass -> water
-    n_a2 = a2_0 * gd0 + a2_1 * gd1 + a2_2 * gd2
-    # bp = unit_vector(a2 - glass_dir * n_glass)
-    bp0 = a2_0 - gd0 * n_glass
-    bp1 = a2_1 - gd1 * n_glass
-    bp2 = a2_2 - gd2 * n_glass
-    bpn = c_sqrt(bp0 * bp0 + bp1 * bp1 + bp2 * bp2)
-    if bpn > 0.0:
-        bp0 /= bpn
-        bp1 /= bpn
-        bp2 /= bpn
-
-    p2 = c_sqrt(1.0 - n_a2 * n_a2) * mm_n2_0 / mm_n3
-    n_final = c_sqrt(1.0 - p2 * p2) if n_a2 >= 0 else -c_sqrt(1.0 - p2 * p2)
-
-    ox = bp0 * p2 + gd0 * n_final
-    oy = bp1 * p2 + gd1 * n_final
-    oz = bp2 * p2 + gd2 * n_final
-
-    out[0] = Xx
-    out[1] = Xy
-    out[2] = Xz
-    out[3] = ox
-    out[4] = oy
-    out[5] = oz
-    return 0
diff --git a/src/openptv2/algorithms/track_kernels_pixel.pxd b/src/openptv2/algorithms/track_kernels_pixel.pxd
index ebc0ae23..b8171f54 100644
--- a/src/openptv2/algorithms/track_kernels_pixel.pxd
+++ b/src/openptv2/algorithms/track_kernels_pixel.pxd
@@ -105,4 +105,15 @@ cpdef int _sorted_candidates_fast_out_nogil(
     double[:] _pp,
 ) noexcept nogil
 
-
+cpdef double _multimed_r_nlay_1layer(
+    double pos_x,
+    double pos_y,
+    double pos_z,
+    double ext_x0,
+    double ext_y0,
+    double ext_z0,
+    double mm_n1,
+    double mm_n2_0,
+    double mm_n3,
+    double mm_d0,
+) noexcept nogil
diff --git a/src/openptv2/algorithms/track_kernels_pixel.py b/src/openptv2/algorithms/track_kernels_pixel.py
index a120e3ef..72bb695c 100644
--- a/src/openptv2/algorithms/track_kernels_pixel.py
+++ b/src/openptv2/algorithms/track_kernels_pixel.py
@@ -45,6 +45,7 @@
 
 @cython.ccall
 @cython.nogil
+@cython.exceptval(check=False)
 def _multimed_r_nlay_1layer(
     pos_x: cython.double,
     pos_y: cython.double,
diff --git a/src/openptv2/algorithms/track_kernels_position.pxd b/src/openptv2/algorithms/track_kernels_position.pxd
index 2de4ac27..65638675 100644
--- a/src/openptv2/algorithms/track_kernels_position.pxd
+++ b/src/openptv2/algorithms/track_kernels_position.pxd
@@ -41,4 +41,22 @@ cpdef int assess_new_position_fast_nogil(
     double[:] scratch,
 ) noexcept nogil
 
+cpdef int _ray_tracing_out(
+    double x,
+    double y,
+    double[:] cal,
+    double[:] out,
+) noexcept nogil
 
+cpdef int _angle_acc_out(
+    double start_x,
+    double start_y,
+    double start_z,
+    double pred_x,
+    double pred_y,
+    double pred_z,
+    double cand_x,
+    double cand_y,
+    double cand_z,
+    double[:] out,
+) noexcept nogil
diff --git a/src/openptv2/algorithms/track_kernels_position.py b/src/openptv2/algorithms/track_kernels_position.py
index 5c3d7d81..51958775 100644
--- a/src/openptv2/algorithms/track_kernels_position.py
+++ b/src/openptv2/algorithms/track_kernels_position.py
@@ -5,10 +5,16 @@
 import numpy as np
 
 if cython.compiled:
+    from cython.cimports.libc.math import (
+        acos as c_acos,
+    )
     from cython.cimports.libc.math import (
         sqrt as c_sqrt,
     )
 else:
+    from math import (
+        acos as c_acos,
+    )
     from math import (
         sqrt as c_sqrt,
     )
@@ -55,6 +61,75 @@
 @cython.cdivision(True)
 @cython.profile(False)
 @cython.nogil
+@cython.exceptval(check=False)
+def _angle_acc_out(
+    start_x: cython.double,
+    start_y: cython.double,
+    start_z: cython.double,
+    pred_x: cython.double,
+    pred_y: cython.double,
+    pred_z: cython.double,
+    cand_x: cython.double,
+    cand_y: cython.double,
+    cand_z: cython.double,
+    out: cython.double[:],
+) -> cython.int:
+    """Write angle and acc to out[0], out[1] — no tuple creation."""
+    v0x: cython.double
+    v0y: cython.double
+    v0z: cython.double
+    v1x: cython.double
+    v1y: cython.double
+    v1z: cython.double
+    angle: cython.double
+    norm0: cython.double
+    norm1: cython.double
+    dot: cython.double
+    dx: cython.double
+    dy: cython.double
+    dz: cython.double
+    acc: cython.double
+    v0x = pred_x - start_x
+    v0y = pred_y - start_y
+    v0z = pred_z - start_z
+    v1x = cand_x - start_x
+    v1y = cand_y - start_y
+    v1z = cand_z - start_z
+
+    if v0x == -v1x and v0y == -v1y and v0z == -v1z:
+        angle = 200.0
+    elif v0x == v1x and v0y == v1y and v0z == v1z:
+        angle = 0.0
+    else:
+        norm0 = c_sqrt(v0x * v0x + v0y * v0y + v0z * v0z)
+        norm1 = c_sqrt(v1x * v1x + v1y * v1y + v1z * v1z)
+        if norm0 == 0.0 or norm1 == 0.0:
+            angle = 0.0
+        else:
+            dot = (v0x * v1x + v0y * v1y + v0z * v1z) / (norm0 * norm1)
+            if dot > 1.0:
+                dot = 1.0
+            elif dot < -1.0:
+                dot = -1.0
+            angle = c_acos(dot) * 200.0 / 3.141592653589793
+
+    dx = v1x - v0x
+    dy = v1y - v0y
+    dz = v1z - v0z
+    acc = c_sqrt(dx * dx + dy * dy + dz * dz)
+    out[0] = angle
+    out[1] = acc
+    return 0
+
+
+@cython.ccall
+@cython.inline
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.cdivision(True)
+@cython.profile(False)
+@cython.nogil
+@cython.exceptval(check=False)
 def _ray_tracing_out(
     x: cython.double,
     y: cython.double,
diff --git a/src/openptv2/algorithms/track_kernels_search.py b/src/openptv2/algorithms/track_kernels_search.py
index eee820bd..63153d2a 100644
--- a/src/openptv2/algorithms/track_kernels_search.py
+++ b/src/openptv2/algorithms/track_kernels_search.py
@@ -7,392 +7,24 @@
 import numpy as np
 
 if cython.compiled:
-    from cython.cimports.libc.math import (
-        asin as c_asin,
-    )
-    from cython.cimports.libc.math import (
-        atan as c_atan,
-    )
-    from cython.cimports.libc.math import (
-        cos as c_cos,
-    )
-    from cython.cimports.libc.math import (
-        sin as c_sin,
-    )
     from cython.cimports.libc.math import (
         sqrt as c_sqrt,
     )
-    from cython.cimports.libc.math import (
-        tan as c_tan,
-    )
 else:
-    from math import (
-        asin as c_asin,
-    )
-    from math import (
-        atan as c_atan,
-    )
-    from math import (
-        cos as c_cos,
-    )
-    from math import (
-        sin as c_sin,
-    )
     from math import (
         sqrt as c_sqrt,
     )
-    from math import (
-        tan as c_tan,
-    )
 
 _M_PI: cython.double = 3.141592653589793
 
-
-@cython.cfunc
-@cython.nogil
-def _multimed_r_nlay_1layer(
-    pos_x: cython.double,
-    pos_y: cython.double,
-    pos_z: cython.double,
-    ext_x0: cython.double,
-    ext_y0: cython.double,
-    ext_z0: cython.double,
-    mm_n1: cython.double,
-    mm_n2_0: cython.double,
-    mm_n3: cython.double,
-    mm_d0: cython.double,
-) -> cython.double:
-    """Single-layer iterative radial shift."""
-    zout: cython.double
-    dx: cython.double
-    dy: cython.double
-    r: cython.double
-    rq: cython.double
-    it: cython.int
-    denom: cython.double
-    beta1: cython.double
-    sin_beta1: cython.double
-    arg: cython.double
-    beta2_0: cython.double
-    arg3: cython.double
-    beta3: cython.double
-    rbeta: cython.double
-    rdiff: cython.double
-    if mm_n1 == 1.0 and mm_n2_0 == 1.0 and mm_n3 == 1.0:
-        return 1.0
-
-    zout = pos_z
-    dx = pos_x - ext_x0
-    dy = pos_y - ext_y0
-    r = c_sqrt(dx * dx + dy * dy)
-    rq = r
-
-    for it in range(40):
-        denom = ext_z0 - pos_z
-        if denom == 0.0:
-            return 1.0
-        beta1 = c_atan(rq / denom)
-        sin_beta1 = c_sin(beta1)
-
-        arg = sin_beta1 * mm_n1 / mm_n2_0
-        if arg > 1.0:
-            arg = 1.0
-        elif arg < -1.0:
-            arg = -1.0
-        beta2_0 = c_asin(arg)
-
-        arg3 = sin_beta1 * mm_n1 / mm_n3
-        if arg3 > 1.0:
-            arg3 = 1.0
-        elif arg3 < -1.0:
-            arg3 = -1.0
-        beta3 = c_asin(arg3)
-
-        rbeta = (
-            (ext_z0 - mm_d0) * c_tan(beta1)
-            + mm_d0 * c_tan(beta2_0)
-            - zout * c_tan(beta3)
-        )
-
-        rdiff = r - rbeta
-        rq += rdiff
-
-        if abs(rdiff) < 0.001:
-            break
-    else:
-        return 1.0
-
-    if r != 0.0:
-        return rq / r
-    return 1.0
-
-
-@cython.cfunc
-@cython.profile(False)
-@cython.nogil
-def _point_to_pixel_out(
-    pos: cython.double[:],
-    cal: cython.double[:],
-    mmlut_data: cython.double[:],
-    mmlut_origin: cython.double[:],
-    mmlut_nr: cython.int,
-    mmlut_nz: cython.int,
-    mmlut_rw: cython.double,
-    has_mmlut: cython.int,
-    imx_half: cython.double,
-    imy_half: cython.double,
-    inv_pix_x: cython.double,
-    inv_pix_y: cython.double,
-    chfield: cython.int,
-    out: cython.double[:],
-) -> cython.int:
-    """Write pixel coordinates to out[0], out[1] — no tuple creation."""
-    pos0: cython.double
-    pos1: cython.double
-    pos2: cython.double
-    ext_x0: cython.double
-    ext_y0: cython.double
-    ext_z0: cython.double
-    dm00: cython.double
-    dm10: cython.double
-    dm20: cython.double
-    dm01: cython.double
-    dm11: cython.double
-    dm21: cython.double
-    dm02: cython.double
-    dm12: cython.double
-    dm22: cython.double
-    int_cc: cython.double
-    xh: cython.double
-    yh: cython.double
-    gx: cython.double
-    gy: cython.double
-    gz: cython.double
-    inv_dog: cython.double
-    mm_n1: cython.double
-    mm_n2_0: cython.double
-    mm_n3: cython.double
-    mm_d0: cython.double
-    k1: cython.double
-    k2: cython.double
-    k3: cython.double
-    p1: cython.double
-    p2: cython.double
-    scx: cython.double
-    she: cython.double
-    dot_cam: cython.double
-    dist_o_glas: cython.double
-    dist_cam_glas: cython.double
-    dot_pos: cython.double
-    dist_point_glas: cython.double
-    s_cam: cython.double
-    cc_x: cython.double
-    cc_y: cython.double
-    cc_z: cython.double
-    s_pt: cython.double
-    cp_x: cython.double
-    cp_y: cython.double
-    cp_z: cython.double
-    ext_t_z0: cython.double
-    s_d: cython.double
-    ag_x: cython.double
-    ag_y: cython.double
-    ag_z: cython.double
-    tmp_x: cython.double
-    tmp_y: cython.double
-    tmp_z: cython.double
-    pos_t_0: cython.double
-    pos_t_2: cython.double
-    radial_shift: cython.double
-    tx: cython.double
-    ty: cython.double
-    tz: cython.double
-    sz: cython.double
-    iz: cython.int
-    R: cython.double
-    sr: cython.double
-    ir: cython.int
-    v0: cython.int
-    v3: cython.int
-    mmf: cython.double
-    X_t: cython.double
-    s_z: cython.double
-    bx: cython.double
-    by: cython.double
-    bz: cython.double
-    s_x: cython.double
-    dx: cython.double
-    dy: cython.double
-    dz: cython.double
-    deno: cython.double
-    x: cython.double
-    y: cython.double
-    r: cython.double
-    r2: cython.double
-    r4: cython.double
-    radial_factor: cython.double
-    xd: cython.double
-    yd: cython.double
-    sin_she: cython.double
-    cos_she: cython.double
-    x_dist: cython.double
-    y_dist: cython.double
-    x_pixel: cython.double
-    y_pixel: cython.double
-
-    pos0 = pos[0]
-    pos1 = pos[1]
-    pos2 = pos[2]
-
-    ext_x0 = cal[0]
-    ext_y0 = cal[1]
-    ext_z0 = cal[2]
-    dm00 = cal[3]
-    dm10 = cal[4]
-    dm20 = cal[5]
-    dm01 = cal[6]
-    dm11 = cal[7]
-    dm21 = cal[8]
-    dm02 = cal[9]
-    dm12 = cal[10]
-    dm22 = cal[11]
-    int_cc = cal[12]
-    xh = cal[13]
-    yh = cal[14]
-    gx = cal[15]
-    gy = cal[16]
-    gz = cal[17]
-    inv_dog = cal[19]
-    mm_n1 = cal[20]
-    mm_n2_0 = cal[21]
-    mm_n3 = cal[22]
-    mm_d0 = cal[23]
-    k1 = cal[24]
-    k2 = cal[25]
-    k3 = cal[26]
-    p1 = cal[27]
-    p2 = cal[28]
-    scx = cal[29]
-    she = cal[30]
-
-    # trans_cam_point
-    dot_cam = ext_x0 * gx + ext_y0 * gy + ext_z0 * gz
-    dist_o_glas = cal[18]
-    dist_cam_glas = dot_cam * inv_dog - dist_o_glas - mm_d0
-
-    dot_pos = pos0 * gx + pos1 * gy + pos2 * gz
-    dist_point_glas = dot_pos * inv_dog - dist_o_glas
-
-    s_cam = dist_cam_glas * inv_dog
-    cc_x = ext_x0 - gx * s_cam
-    cc_y = ext_y0 - gy * s_cam
-    cc_z = ext_z0 - gz * s_cam
-
-    s_pt = dist_point_glas * inv_dog
-    cp_x = pos0 - gx * s_pt
-    cp_y = pos1 - gy * s_pt
-    cp_z = pos2 - gz * s_pt
-
-    ext_t_z0 = dist_cam_glas + mm_d0
-
-    s_d = mm_d0 * inv_dog
-    ag_x = cc_x - gx * s_d
-    ag_y = cc_y - gy * s_d
-    ag_z = cc_z - gz * s_d
-    tmp_x = cp_x - ag_x
-    tmp_y = cp_y - ag_y
-    tmp_z = cp_z - ag_z
-
-    pos_t_0 = c_sqrt(tmp_x * tmp_x + tmp_y * tmp_y + tmp_z * tmp_z)
-    pos_t_2 = dist_point_glas
-
-    # mmlut lookup + multimed_nlay
-    radial_shift = 1.0
-    if has_mmlut:
-        tx = pos_t_0 - mmlut_origin[0]
-        ty = -mmlut_origin[1]
-        tz = pos_t_2 - mmlut_origin[2]
-        sz = tz / mmlut_rw
-        iz = int(sz)
-        sz -= iz
-        R = c_sqrt(tx * tx + ty * ty)
-        sr = R / mmlut_rw
-        ir = int(sr)
-        sr -= ir
-        if ir <= mmlut_nr and iz >= 0 and iz <= mmlut_nz:
-            v0 = ir * mmlut_nz + iz
-            v3 = v0 + mmlut_nz + 1
-            if v0 >= 0 and v3 <= mmlut_nr * mmlut_nz:
-                mmf = (
-                    mmlut_data[v0] * (1.0 - sr) * (1.0 - sz)
-                    + mmlut_data[v0 + 1] * (1.0 - sr) * sz
-                    + mmlut_data[v0 + mmlut_nz] * sr * (1.0 - sz)
-                    + mmlut_data[v3] * sr * sz
-                )
-                if mmf > 0.0:
-                    radial_shift = mmf
-    if radial_shift == 1.0:
-        radial_shift = _multimed_r_nlay_1layer(
-            pos_t_0,
-            0.0,
-            pos_t_2,
-            0.0,
-            0.0,
-            ext_t_z0,
-            mm_n1,
-            mm_n2_0,
-            mm_n3,
-            mm_d0,
-        )
-    X_t = pos_t_0 * radial_shift
-
-    # back_trans_point
-    s_z = -pos_t_2 * inv_dog
-    bx = ag_x - gx * s_z
-    by = ag_y - gy * s_z
-    bz = ag_z - gz * s_z
-    if pos_t_0 > 0.0:
-        s_x = -X_t / pos_t_0
-        bx -= tmp_x * s_x
-        by -= tmp_y * s_x
-        bz -= tmp_z * s_x
-
-    # perspective projection
-    dx = bx - ext_x0
-    dy = by - ext_y0
-    dz = bz - ext_z0
-    deno = dm02 * dx + dm12 * dy + dm22 * dz
-    x = -int_cc * (dm00 * dx + dm10 * dy + dm20 * dz) / deno
-    y = -int_cc * (dm01 * dx + dm11 * dy + dm21 * dz) / deno
-
-    # flat_to_dist + distort_brown_affin
-    x += xh
-    y += yh
-    r = c_sqrt(x * x + y * y)
-    if r < 1e-10:
-        x_dist = 0.0
-        y_dist = 0.0
-    else:
-        r2 = r * r
-        r4 = r2 * r2
-        radial_factor = 1.0 + k1 * r2 + k2 * r4 + k3 * r4 * r2
-        xd = x * radial_factor + p1 * (r2 + 2.0 * x * x) + 2.0 * p2 * x * y
-        yd = y * radial_factor + p2 * (r2 + 2.0 * y * y) + 2.0 * p1 * x * y
-        sin_she = c_sin(she)
-        cos_she = c_cos(she)
-        x_dist = scx * (xd - sin_she * yd)
-        y_dist = scx * cos_she * yd
-
-    # metric_to_pixel
-    x_pixel = x_dist * inv_pix_x + imx_half
-    y_pixel = imy_half - y_dist * inv_pix_y
-    if chfield == 1:
-        y_pixel = (y_pixel - 1.0) * 0.5
-    elif chfield == 2:
-        y_pixel = y_pixel * 0.5
-    out[0] = x_pixel
-    out[1] = y_pixel
-    return 0
+if cython.compiled:
+    from cython.cimports.openptv2.algorithms.track_kernels_pixel import (
+        _point_to_pixel_out,
+    )
+else:
+    from .track_kernels_pixel import (
+        _point_to_pixel_out,
+    )
 
 
 # Sentinel values for unused particle/candidate indices — typed C int
@@ -622,89 +254,6 @@ def candsearch_in_pix_rest_fast(
     return best, counter
 
 
-@cython.ccall
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def sort_candidates_by_freq_fast(
-    ftnr: cython.int[:],
-    freq: cython.int[:],
-    whichcam: cython.int[:, ::1],
-    n: cython.int,
-    num_cams: cython.int,
-    max_cands: cython.int,
-):
-    """Sort candidates by frequency, matches C algorithm.
-
-    Args:
-        ftnr: (n,) int32 — candidate target numbers (TR_UNUSED = -1).
-        freq: (n,) int32 — frequency counts (zeroed on entry).
-        whichcam: (n, num_cams) int32 — camera flags.
-        n: total number of entries (num_cams * max_cands).
-        num_cams: number of cameras.
-        max_cands: candidates per camera (4).
-
-    Returns:
-        num_valid: number of valid candidates after sort.
-    """
-    i: cython.int
-    j: cython.int
-    m: cython.int
-    k: cython.int
-    ftnr_i: cython.int
-    num_valid: cython.int
-    tr_unused = -1
-
-    for i in range(n):
-        ftnr_i = ftnr[i]
-        if ftnr_i == tr_unused:
-            continue
-        for j in range(num_cams):
-            for m in range(max_cands):
-                if ftnr_i == ftnr[max_cands * j + m]:
-                    whichcam[i, j] = 1
-
-    for i in range(n):
-        if ftnr[i] != tr_unused:
-            for j in range(num_cams):
-                if whichcam[i, j] == 1:
-                    freq[i] += 1
-
-    for i in range(1, n):
-        for j in range(n - 1, i - 1, -1):
-            if freq[j - 1] < freq[j]:
-                ftnr[j - 1], ftnr[j] = ftnr[j], ftnr[j - 1]
-                freq[j - 1], freq[j] = freq[j], freq[j - 1]
-                for k in range(num_cams):
-                    whichcam[j - 1, k], whichcam[j, k] = (
-                        whichcam[j, k],
-                        whichcam[j - 1, k],
-                    )
-
-    for i in range(n):
-        ftnr_i = ftnr[i]
-        for j in range(i + 1, n):
-            if ftnr[j] == ftnr_i or freq[j] < 2:
-                freq[j] = 0
-                ftnr[j] = tr_unused
-
-    for i in range(1, n):
-        for j in range(n - 1, i - 1, -1):
-            if freq[j - 1] < freq[j]:
-                ftnr[j - 1], ftnr[j] = ftnr[j], ftnr[j - 1]
-                freq[j - 1], freq[j] = freq[j], freq[j - 1]
-                for k in range(num_cams):
-                    whichcam[j - 1, k], whichcam[j, k] = (
-                        whichcam[j, k],
-                        whichcam[j - 1, k],
-                    )
-
-    num_valid = 0
-    for i in range(n):
-        if freq[i] != 0:
-            num_valid += 1
-    return num_valid
-
-
 @cython.ccall
 @cython.boundscheck(False)
 @cython.wraparound(False)
@@ -1028,406 +577,3 @@ def _sorted_candidates_fast_out(
             num_valid += 1
     return num_valid
 
-
-@cython.ccall
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.nogil
-def candsearch_in_pix_fast_nogil(
-    targ_x: cython.double[:],
-    targ_y: cython.double[:],
-    targ_tnr: cython.int[:],
-    num_targets: cython.int,
-    cent_x: cython.double,
-    cent_y: cython.double,
-    dl: cython.double,
-    dr: cython.double,
-    du: cython.double,
-    dd: cython.double,
-    imx: cython.double,
-    imy: cython.double,
-    tr_unused: cython.int,
-    out_indices: cython.int[:],
-) -> cython.int:
-    xmin: cython.double
-    xmax: cython.double
-    ymin: cython.double
-    ymax: cython.double
-    p1: cython.int
-    p2: cython.int
-    p3: cython.int
-    p4: cython.int
-    d1: cython.double
-    d2: cython.double
-    d3: cython.double
-    d4: cython.double
-    j0: cython.int
-    dj: cython.int
-    j: cython.int
-    ty: cython.double
-    tx: cython.double
-    dx: cython.double
-    dy: cython.double
-    d: cython.double
-
-    xmin = cent_x - dl
-    xmax = cent_x + dr
-    ymin = cent_y - du
-    ymax = cent_y + dd
-
-    if xmin < 0.0:
-        xmin = 0.0
-    if xmax > imx:
-        xmax = imx
-    if ymin < 0.0:
-        ymin = 0.0
-    if ymax > imy:
-        ymax = imy
-
-    p1 = -999  # PT_UNUSED is -999
-    p2 = -999
-    p3 = -999
-    p4 = -999
-    d1 = 1e20
-    d2 = 1e20
-    d3 = 1e20
-    d4 = 1e20
-
-    if not (0.0 <= cent_x <= imx and 0.0 <= cent_y <= imy):
-        out_indices[0] = p1
-        out_indices[1] = p2
-        out_indices[2] = p3
-        out_indices[3] = p4
-        return 0
-
-    j0 = num_targets // 2
-    dj = num_targets // 4
-    while dj > 1:
-        if targ_y[j0] < ymin:
-            j0 += dj
-        else:
-            j0 -= dj
-        dj //= 2
-
-    j0 -= 12
-    if j0 < 0:
-        j0 = 0
-
-    for j in range(j0, num_targets):
-        ty = targ_y[j]
-        if targ_tnr[j] != tr_unused:
-            if ty > ymax:
-                break
-            tx = targ_x[j]
-            if tx > xmin and tx < xmax and ty > ymin and ty < ymax:
-                dx = cent_x - tx
-                dy = cent_y - ty
-                d = c_sqrt(dx * dx + dy * dy)
-
-                if d < d1:
-                    p4 = p3
-                    p3 = p2
-                    p2 = p1
-                    p1 = j
-                    d4 = d3
-                    d3 = d2
-                    d2 = d1
-                    d1 = d
-                elif d < d2:
-                    p4 = p3
-                    p3 = p2
-                    p2 = j
-                    d4 = d3
-                    d3 = d2
-                    d2 = d
-                elif d < d3:
-                    p4 = p3
-                    p3 = j
-                    d4 = d3
-                    d3 = d
-                elif d < d4:
-                    p4 = j
-                    d4 = d
-
-    out_indices[0] = p1
-    out_indices[1] = p2
-    out_indices[2] = p3
-    out_indices[3] = p4
-    return 0
-
-
-@cython.ccall
-@cython.inline
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.nogil
-def _sorted_candidates_fast_out_nogil(
-    center: cython.double[:],
-    center_proj_x: cython.double[:],
-    center_proj_y: cython.double[:],
-    num_cams: cython.int,
-    max_cands: cython.int,
-    cal_arr: cython.double[:, ::1],
-    md0: cython.double[:],
-    md1: cython.double[:],
-    md2: cython.double[:],
-    md3: cython.double[:],
-    md4: cython.double[:],
-    md5: cython.double[:],
-    md6: cython.double[:],
-    md7: cython.double[:],
-    mo_arr: cython.double[:, ::1],
-    mnr_arr: cython.int[:],
-    mnz_arr: cython.int[:],
-    mrw_arr: cython.double[:],
-    targ_x: cython.double[:, ::1],
-    targ_y: cython.double[:, ::1],
-    targ_tnr: cython.int[:, ::1],
-    num_targets: cython.int[:],
-    dvxmin: cython.double,
-    dvxmax: cython.double,
-    dvymin: cython.double,
-    dvymax: cython.double,
-    dvzmin: cython.double,
-    dvzmax: cython.double,
-    imx_half: cython.double,
-    imy_half: cython.double,
-    inv_pix_x: cython.double,
-    inv_pix_y: cython.double,
-    chfield: cython.int,
-    imx: cython.double,
-    imy: cython.double,
-    tr_unused: cython.int,
-    ftnr_out: cython.int[:],
-    freq_out: cython.int[:],
-    whichcam_out: cython.int[:, :],
-) -> cython.int:
-    n: cython.int
-    px: cython.double
-    py: cython.double
-    pz: cython.double
-    i: cython.int
-    pt: cython.int
-    xr_i: cython.double
-    xl_i: cython.double
-    yd_i: cython.double
-    yu_i: cython.double
-    cx: cython.double
-    cy: cython.double
-    corner_x: cython.double
-    corner_y: cython.double
-    mrw: cython.double
-    mnr: cython.int
-    mnz: cython.int
-    has_mmlut: cython.int
-    cam: cython.int
-    base: cython.int
-    ci: cython.int
-    idx: cython.int
-    ftnr_i: cython.int
-    num_valid: cython.int
-    j: cython.int
-    m: cython.int
-    k: cython.int
-    _pp: cython.double[:]
-    quader_buf: cython.double[:]
-    pt_buf: cython.double[:]
-    with cython.gil:
-        _pp_buf = np.zeros(2, dtype=np.float64)
-        _pp = _pp_buf
-        _quader_buf = np.zeros(24, dtype=np.float64)
-        quader_buf = _quader_buf
-        _pt_buf = np.zeros(3, dtype=np.float64)
-        pt_buf = _pt_buf
-
-    n = num_cams * max_cands
-
-    # --- searchquader inlined ---
-    px = center[0]
-    py = center[1]
-    pz = center[2]
-    for pt in range(8):
-        quader_buf[pt * 3 + 0] = px + (dvxmax if pt & 1 else dvxmin)
-        quader_buf[pt * 3 + 1] = py + (dvymax if pt & 2 else dvymin)
-        quader_buf[pt * 3 + 2] = pz + (dvzmax if pt & 4 else dvzmin)
-
-    xr: cython.double[:]
-    xl: cython.double[:]
-    yd: cython.double[:]
-    yu: cython.double[:]
-    with cython.gil:
-        _xr_buf = np.zeros(8, dtype=np.float64)
-        xr = _xr_buf
-        _xl_buf = np.zeros(8, dtype=np.float64)
-        xl = _xl_buf
-        _yd_buf = np.zeros(8, dtype=np.float64)
-        yd = _yd_buf
-        _yu_buf = np.zeros(8, dtype=np.float64)
-        yu = _yu_buf
-
-    for i in range(num_cams):
-        cal = cal_arr[i]
-        mo = mo_arr[i]
-        mnr = mnr_arr[i]
-        mnz = mnz_arr[i]
-        mrw = mrw_arr[i]
-        has_mmlut = mnr > 0
-
-        # Select pre-unpacked md memoryview without GIL
-        md: cython.double[:]
-        if i == 0:
-            md = md0
-        elif i == 1:
-            md = md1
-        elif i == 2:
-            md = md2
-        elif i == 3:
-            md = md3
-        elif i == 4:
-            md = md4
-        elif i == 5:
-            md = md5
-        elif i == 6:
-            md = md6
-        else:
-            md = md7
-
-        xr_i = 0.0
-        xl_i = float(imx)
-        yd_i = 0.0
-        yu_i = float(imy)
-        # Use pre-computed center projection
-        cx = center_proj_x[i]
-        cy = center_proj_y[i]
-        for pt in range(8):
-            pt_buf[0] = quader_buf[pt * 3 + 0]
-            pt_buf[1] = quader_buf[pt * 3 + 1]
-            pt_buf[2] = quader_buf[pt * 3 + 2]
-            _point_to_pixel_out(
-                pt_buf,
-                cal,
-                md,
-                mo,
-                mnr,
-                mnz,
-                mrw,
-                has_mmlut,
-                imx_half,
-                imy_half,
-                inv_pix_x,
-                inv_pix_y,
-                chfield,
-                _pp,
-            )
-            corner_x = _pp[0]
-            corner_y = _pp[1]
-            if corner_x < xl_i:
-                xl_i = corner_x
-            if corner_y < yu_i:
-                yu_i = corner_y
-            if corner_x > xr_i:
-                xr_i = corner_x
-            if corner_y > yd_i:
-                yd_i = corner_y
-        if xl_i < 0.0:
-            xl_i = 0.0
-        if yu_i < 0.0:
-            yu_i = 0.0
-        if xr_i > imx:
-            xr_i = imx
-        if yd_i > imy:
-            yd_i = imy
-        xr[i] = xr_i - cx
-        xl[i] = cx - xl_i
-        yd[i] = yd_i - cy
-        yu[i] = cy - yu_i
-
-    # --- initialize output buffers ---
-    for i in range(n):
-        ftnr_out[i] = tr_unused
-        freq_out[i] = 0
-        for j in range(num_cams):
-            whichcam_out[i, j] = 0
-
-    # Local buffer for candsearch_in_pix_fast_nogil
-    cands_buf: cython.int[:]
-    with cython.gil:
-        _cands_buf = np.zeros(4, dtype=np.int32)
-        cands_buf = _cands_buf
-
-    # --- candsearch per camera, write directly into ftnr_out/whichcam_out ---
-    for cam in range(num_cams):
-        candsearch_in_pix_fast_nogil(
-            targ_x[cam],
-            targ_y[cam],
-            targ_tnr[cam],
-            num_targets[cam],
-            center_proj_x[cam],
-            center_proj_y[cam],
-            xl[cam],
-            xr[cam],
-            yu[cam],
-            yd[cam],
-            imx,
-            imy,
-            tr_unused,
-            cands_buf,
-        )
-
-        base = cam * max_cands
-        for ci in range(4):
-            idx = cands_buf[ci]
-            if idx != -999:  # PT_UNUSED is -999
-                whichcam_out[base + ci, cam] = 1
-                ftnr_out[base + ci] = targ_tnr[cam, idx]
-
-    # --- sort_candidates_by_freq inlined ---
-    for i in range(n):
-        ftnr_i = ftnr_out[i]
-        if ftnr_i == tr_unused:
-            continue
-        for j in range(num_cams):
-            for m in range(max_cands):
-                if ftnr_i == ftnr_out[max_cands * j + m]:
-                    whichcam_out[i, j] = 1
-
-    for i in range(n):
-        if ftnr_out[i] != tr_unused:
-            for j in range(num_cams):
-                if whichcam_out[i, j] == 1:
-                    freq_out[i] += 1
-
-    for i in range(1, n):
-        for j in range(n - 1, i - 1, -1):
-            if freq_out[j - 1] < freq_out[j]:
-                ftnr_out[j - 1], ftnr_out[j] = ftnr_out[j], ftnr_out[j - 1]
-                freq_out[j - 1], freq_out[j] = freq_out[j], freq_out[j - 1]
-                for k in range(num_cams):
-                    whichcam_out[j - 1, k], whichcam_out[j, k] = (
-                        whichcam_out[j, k],
-                        whichcam_out[j - 1, k],
-                    )
-
-    for i in range(n):
-        ftnr_i = ftnr_out[i]
-        for j in range(i + 1, n):
-            if ftnr_out[j] == ftnr_i or freq_out[j] < 2:
-                freq_out[j] = 0
-                ftnr_out[j] = tr_unused
-
-    for i in range(1, n):
-        for j in range(n - 1, i - 1, -1):
-            if freq_out[j - 1] < freq_out[j]:
-                ftnr_out[j - 1], ftnr_out[j] = ftnr_out[j], ftnr_out[j - 1]
-                freq_out[j - 1], freq_out[j] = freq_out[j], freq_out[j - 1]
-                for k in range(num_cams):
-                    whichcam_out[j - 1, k], whichcam_out[j, k] = (
-                        whichcam_out[j, k],
-                        whichcam_out[j - 1, k],
-                    )
-
-    num_valid = 0
-    for i in range(n):
-        if freq_out[i] != 0:
-            num_valid += 1
-    return num_valid
diff --git a/src/openptv2/algorithms/track_kernels_track3d.py b/src/openptv2/algorithms/track_kernels_track3d.py
index 7918b580..b0b167fe 100644
--- a/src/openptv2/algorithms/track_kernels_track3d.py
+++ b/src/openptv2/algorithms/track_kernels_track3d.py
@@ -8,9 +8,13 @@
 UNSUPPORTED_PENALTY = 1e6
 
 if cython.compiled:
-    from cython.cimports.libc.math import floor as c_floor, sqrt as c_sqrt, isnan as c_isnan
+    from cython.cimports.libc.math import floor as c_floor
+    from cython.cimports.libc.math import isnan as c_isnan
+    from cython.cimports.libc.math import sqrt as c_sqrt
 else:
-    from math import floor as c_floor, sqrt as c_sqrt, isnan as c_isnan
+    from math import floor as c_floor
+    from math import isnan as c_isnan
+    from math import sqrt as c_sqrt
 
 
 @cython.cfunc
@@ -137,9 +141,9 @@ def _find_closest_in_3d(
 ) -> cython.int:
     """Find up to max_cands closest candidates by distance within a 3D box.
 
-    @cython.ccall rather than @cython.cfunc: track_kernels_tracking re-exports
-    this one, so it has to stay importable from Python while still being
-    C-callable from _find_closest_in_3d_grid's small-frame fallback.
+    @cython.ccall rather than @cython.cfunc: this one has to stay importable
+    from Python (tests import it) while still being C-callable from
+    _find_closest_in_3d_grid's small-frame fallback.
     """
     s: cython.int
     k: cython.int
@@ -202,6 +206,19 @@ def track3d_loop_fast(
     dz: cython.double,
     max_cands: cython.int,
     cold_start_gate: cython.double = 1.0,
+    # Prototype shared-observation (Level 1 only):
+    # share_tol > 0 enables recording (not claiming) of contested edges:
+    # when the cheapest remaining edge (i -> k) finds k already taken AND
+    # both this edge and the winner's edge cost less than share_tol, the
+    # pair is recorded into shared_i/shared_k (capacity len(shared_i))
+    # instead of stealing. The driver materializes each pair as a virtual
+    # carrier particle (shared coordinates, own history) so the track
+    # continues with an honest velocity. Cold levels never share: a track
+    # with no history of its own must not twin.
+    share_tol: cython.double = 0.0,
+    shared_count: cython.int[:] = None,
+    shared_i: cython.int[:] = None,
+    shared_k: cython.int[:] = None,
 ):
     """Full track3d loop (3 levels) — single compiled entry.
 
@@ -254,9 +271,22 @@ def track3d_loop_fast(
     oi: cython.int
     e: cython.int
     order: cython.int[:]
+    use_share: cython.bint
+    share_cap: cython.int
+    w: cython.int
 
     count1 = 0
     np2 = num_parts_2
+    use_share = (
+        share_tol > 0.0
+        and shared_count is not None
+        and shared_i is not None
+        and shared_k is not None
+    )
+    share_cap = shared_i.shape[0] if use_share else 0
+
+    _claim_cost_2 = np.full(np2 if np2 > 0 else 1, np.inf, dtype=np.float64)
+    claim_cost_2: cython.double[:] = _claim_cost_2
 
     _cand_inds = np.empty(max_cands, dtype=np.int32)
     _cand_dists = np.empty(max_cands, dtype=np.float64)
@@ -370,7 +400,33 @@ def track3d_loop_fast(
             if path_next_1[i] < 0 and path_prev_2[k] < 0:
                 path_next_1[i] = k
                 path_prev_2[k] = i
+                claim_cost_2[k] = edge_cost[e]
                 count1 += 1
+            elif (
+                use_share
+                and path_next_1[i] < 0
+                and path_prev_2[k] >= 0
+                and edge_cost[e] < share_tol
+                and claim_cost_2[k] < share_tol
+                and shared_count[0] < share_cap
+            ):
+                # Prototype shared-observation (Level 1 only): contested and
+                # mutually well-predicted -> share, don't steal. Cold levels
+                # never share: a track with no history of its own must not
+                # twin. A sharer that later claims normally keeps the claim;
+                # its share record dies below (no forks).
+                shared_i[shared_count[0]] = i
+                shared_k[shared_count[0]] = k
+                shared_count[0] += 1
+
+    if use_share:
+        w = 0
+        for e in range(shared_count[0]):
+            if path_next_1[shared_i[e]] < 0:
+                shared_i[w] = shared_i[e]
+                shared_k[w] = shared_k[e]
+                w += 1
+        shared_count[0] = w
 
     # ===== Level 2: No previous link, neighbor velocity =====
     n_edges = 0
diff --git a/src/openptv2/algorithms/track_kernels_tracking.py b/src/openptv2/algorithms/track_kernels_tracking.py
deleted file mode 100644
index 0dfc50e1..00000000
--- a/src/openptv2/algorithms/track_kernels_tracking.py
+++ /dev/null
@@ -1,38 +0,0 @@
-# ruff: noqa: E402
-"""Compatibility re-export — content split into focused sub-modules 2026-07-10."""
-
-# These mirror the cython.declare() C-level constants in track_kernels_corr,
-# which are not importable from Python when compiled.
-PT_UNUSED = -999
-POSI_K = 80
-MAX_CANDS_K = 32
-TR_UNUSED_K = -1  # noqa: E702
-CORRES_NONE_K = -1
-PREV_NONE_K = -1
-NEXT_NONE_K = -2  # noqa: E702
-COORD_UNUSED_K = -1e10
-ADD_PART_K = 3.0  # noqa: E702
-
-from .track_kernels_corr import (  # noqa: F401, E402
-    trackback_loop_fast,
-    trackcorr_loop_fast,
-)
-from .track_kernels_geom import _angle_acc_out, _ray_tracing_out  # noqa: F401
-from .track_kernels_pixel import (  # noqa: F401
-    _candsearch_in_pix_rest_nogil,
-    _dist_to_flat_out,
-    _multimed_r_nlay_1layer,
-    _pixel_to_metric_out,
-    _point_to_pixel_out,
-    _sorted_candidates_fast_out_nogil,
-    candsearch_in_pix_fast_nogil,
-)
-from .track_kernels_position import (  # noqa: F401
-    _point_position_out,
-    assess_new_position_fast_nogil,
-)
-from .track_kernels_track3d import (  # noqa: F401
-    _find_closest_in_3d,
-    track3d_loop_fast,
-    track4be_loop_fast,
-)
diff --git a/src/openptv2/algorithms/track_kernels_transform.py b/src/openptv2/algorithms/track_kernels_transform.py
index 2c81b607..bd897238 100644
--- a/src/openptv2/algorithms/track_kernels_transform.py
+++ b/src/openptv2/algorithms/track_kernels_transform.py
@@ -8,50 +8,25 @@
 import numpy as np
 
 if cython.compiled:
-    from cython.cimports.libc.math import (
-        asin as c_asin,
-    )
-    from cython.cimports.libc.math import (
-        atan as c_atan,
-    )
-    from cython.cimports.libc.math import (
-        cos as c_cos,
-    )
-    from cython.cimports.libc.math import (
-        sin as c_sin,
-    )
-    from cython.cimports.libc.math import (
-        sqrt as c_sqrt,
-    )
-    from cython.cimports.libc.math import (
-        tan as c_tan,
-    )
+    pass
 else:
-    from math import (
-        asin as c_asin,
-    )
-    from math import (
-        atan as c_atan,
-    )
-    from math import (
-        cos as c_cos,
-    )
-    from math import (
-        sin as c_sin,
-    )
-    from math import (
-        sqrt as c_sqrt,
-    )
-    from math import (
-        tan as c_tan,
-    )
+    pass
 
 _M_PI: cython.double = 3.141592653589793
 
 
-from .track_kernels_geom import (
-    _point_to_pixel_out,
-)
+if cython.compiled:
+    from cython.cimports.openptv2.algorithms.track_kernels_pixel import (
+        _dist_to_flat_out,
+        _pixel_to_metric_out,
+        _point_to_pixel_out,
+    )
+else:
+    from .track_kernels_pixel import (
+        _dist_to_flat_out,
+        _pixel_to_metric_out,
+        _point_to_pixel_out,
+    )
 from .track_kernels_search import (
     candsearch_in_pix_rest_fast,
 )
@@ -65,655 +40,6 @@
 COORD_UNUSED = -1e10
 
 
-@cython.ccall
-@cython.nogil
-def _multimed_r_nlay_1layer(
-    pos_x: cython.double,
-    pos_y: cython.double,
-    pos_z: cython.double,
-    ext_x0: cython.double,
-    ext_y0: cython.double,
-    ext_z0: cython.double,
-    mm_n1: cython.double,
-    mm_n2_0: cython.double,
-    mm_n3: cython.double,
-    mm_d0: cython.double,
-) -> cython.double:
-    """Single-layer iterative radial shift."""
-    zout: cython.double
-    dx: cython.double
-    dy: cython.double
-    r: cython.double
-    rq: cython.double
-    it: cython.int
-    denom: cython.double
-    beta1: cython.double
-    sin_beta1: cython.double
-    arg: cython.double
-    beta2_0: cython.double
-    arg3: cython.double
-    beta3: cython.double
-    rbeta: cython.double
-    rdiff: cython.double
-    if mm_n1 == 1.0 and mm_n2_0 == 1.0 and mm_n3 == 1.0:
-        return 1.0
-
-    zout = pos_z
-    dx = pos_x - ext_x0
-    dy = pos_y - ext_y0
-    r = c_sqrt(dx * dx + dy * dy)
-    rq = r
-
-    for it in range(40):
-        denom = ext_z0 - pos_z
-        if denom == 0.0:
-            return 1.0
-        beta1 = c_atan(rq / denom)
-        sin_beta1 = c_sin(beta1)
-
-        arg = sin_beta1 * mm_n1 / mm_n2_0
-        if arg > 1.0:
-            arg = 1.0
-        elif arg < -1.0:
-            arg = -1.0
-        beta2_0 = c_asin(arg)
-
-        arg3 = sin_beta1 * mm_n1 / mm_n3
-        if arg3 > 1.0:
-            arg3 = 1.0
-        elif arg3 < -1.0:
-            arg3 = -1.0
-        beta3 = c_asin(arg3)
-
-        rbeta = (
-            (ext_z0 - mm_d0) * c_tan(beta1)
-            + mm_d0 * c_tan(beta2_0)
-            - zout * c_tan(beta3)
-        )
-
-        rdiff = r - rbeta
-        rq += rdiff
-
-        if abs(rdiff) < 0.001:
-            break
-    else:
-        return 1.0
-
-    if r != 0.0:
-        return rq / r
-    return 1.0
-
-
-@cython.ccall
-@cython.inline
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.cdivision(True)
-@cython.profile(False)
-@cython.nogil
-def _ray_tracing_out(
-    x: cython.double,
-    y: cython.double,
-    cal: cython.double[:],
-    out: cython.double[:],
-) -> cython.int:
-    """Write ray tracing results into out[0:6] — no tuple creation."""
-    ext_x0: cython.double
-    ext_y0: cython.double
-    ext_z0: cython.double
-    dm00: cython.double
-    dm10: cython.double
-    dm20: cython.double
-    dm01: cython.double
-    dm11: cython.double
-    dm21: cython.double
-    dm02: cython.double
-    dm12: cython.double
-    dm22: cython.double
-    int_cc: cython.double
-    gx: cython.double
-    gy: cython.double
-    gz: cython.double
-    mm_n1: cython.double
-    mm_n2_0: cython.double
-    mm_n3: cython.double
-    mm_d0: cython.double
-    t0: cython.double
-    t1: cython.double
-    t2: cython.double
-    tn: cython.double
-    sd0: cython.double
-    sd1: cython.double
-    sd2: cython.double
-    gn: cython.double
-    gd0: cython.double
-    gd1: cython.double
-    gd2: cython.double
-    c: cython.double
-    dcg: cython.double
-    denom: cython.double
-    d1: cython.double
-    Xb0: cython.double
-    Xb1: cython.double
-    Xb2: cython.double
-    n: cython.double
-    bp0: cython.double
-    bp1: cython.double
-    bp2: cython.double
-    bpn: cython.double
-    p: cython.double
-    n_glass: cython.double
-    a2_0: cython.double
-    a2_1: cython.double
-    a2_2: cython.double
-    d2_denom: cython.double
-    d2: cython.double
-    Xx: cython.double
-    Xy: cython.double
-    Xz: cython.double
-    n_a2: cython.double
-    p2: cython.double
-    n_final: cython.double
-    ox: cython.double
-    oy: cython.double
-    oz: cython.double
-    ext_x0 = cal[0]
-    ext_y0 = cal[1]
-    ext_z0 = cal[2]
-    dm00 = cal[3]
-    dm10 = cal[4]
-    dm20 = cal[5]
-    dm01 = cal[6]
-    dm11 = cal[7]
-    dm21 = cal[8]
-    dm02 = cal[9]
-    dm12 = cal[10]
-    dm22 = cal[11]
-    int_cc = cal[12]
-    gx = cal[15]
-    gy = cal[16]
-    gz = cal[17]
-    mm_n1 = cal[20]
-    mm_n2_0 = cal[21]
-    mm_n3 = cal[22]
-    mm_d0 = cal[23]
-
-    t0 = x
-    t1 = y
-    t2 = -int_cc
-    tn = c_sqrt(t0 * t0 + t1 * t1 + t2 * t2)
-    if tn > 0.0:
-        t0 /= tn
-        t1 /= tn
-        t2 /= tn
-
-    sd0 = dm00 * t0 + dm01 * t1 + dm02 * t2
-    sd1 = dm10 * t0 + dm11 * t1 + dm12 * t2
-    sd2 = dm20 * t0 + dm21 * t1 + dm22 * t2
-
-    gn = c_sqrt(gx * gx + gy * gy + gz * gz)
-    if gn > 0.0:
-        gd0 = gx / gn
-        gd1 = gy / gn
-        gd2 = gz / gn
-    else:
-        gd0 = 0.0
-        gd1 = 0.0
-        gd2 = 0.0
-    c = gn + mm_d0
-
-    dcg = gd0 * ext_x0 + gd1 * ext_y0 + gd2 * ext_z0 - c
-    denom = gd0 * sd0 + gd1 * sd1 + gd2 * sd2
-    d1 = -dcg / denom
-
-    Xb0 = ext_x0 + sd0 * d1
-    Xb1 = ext_y0 + sd1 * d1
-    Xb2 = ext_z0 + sd2 * d1
-
-    n = sd0 * gd0 + sd1 * gd1 + sd2 * gd2
-    bp0 = sd0 - gd0 * n
-    bp1 = sd1 - gd1 * n
-    bp2 = sd2 - gd2 * n
-    bpn = c_sqrt(bp0 * bp0 + bp1 * bp1 + bp2 * bp2)
-    if bpn > 0.0:
-        bp0 /= bpn
-        bp1 /= bpn
-        bp2 /= bpn
-
-    p = c_sqrt(1.0 - n * n) * mm_n1 / mm_n2_0
-    n_glass = c_sqrt(1.0 - p * p) if n >= 0 else -c_sqrt(1.0 - p * p)
-
-    a2_0 = bp0 * p + gd0 * n_glass
-    a2_1 = bp1 * p + gd1 * n_glass
-    a2_2 = bp2 * p + gd2 * n_glass
-
-    d2_denom = gd0 * a2_0 + gd1 * a2_1 + gd2 * a2_2
-    d2 = mm_d0 / abs(d2_denom)
-
-    Xx = Xb0 + a2_0 * d2
-    Xy = Xb1 + a2_1 * d2
-    Xz = Xb2 + a2_2 * d2
-
-    n_a2 = a2_0 * gd0 + a2_1 * gd1 + a2_2 * gd2
-    bp0 = a2_0 - gd0 * n_glass
-    bp1 = a2_1 - gd1 * n_glass
-    bp2 = a2_2 - gd2 * n_glass
-    bpn = c_sqrt(bp0 * bp0 + bp1 * bp1 + bp2 * bp2)
-    if bpn > 0.0:
-        bp0 /= bpn
-        bp1 /= bpn
-        bp2 /= bpn
-
-    p2 = c_sqrt(1.0 - n_a2 * n_a2) * mm_n2_0 / mm_n3
-    n_final = c_sqrt(1.0 - p2 * p2) if n_a2 >= 0 else -c_sqrt(1.0 - p2 * p2)
-
-    ox = bp0 * p2 + gd0 * n_final
-    oy = bp1 * p2 + gd1 * n_final
-    oz = bp2 * p2 + gd2 * n_final
-
-    out[0] = Xx
-    out[1] = Xy
-    out[2] = Xz
-    out[3] = ox
-    out[4] = oy
-    out[5] = oz
-    return 0
-
-
-@cython.ccall
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.nogil
-def _point_position_out(
-    targets: cython.double[:, ::1],
-    num_cams: cython.int,
-    cal_arr: cython.double[:, ::1],
-    out: cython.double[:],
-    scratch_ray: cython.double[:],
-) -> cython.double:
-    """Internal — writes 3D position to out[0:3], returns avg_dist (scalar).
-
-    Pure C entry — zero Python object creation, zero tuple overhead.
-    """
-    cam: cython.int
-    pair: cython.int
-    tx: cython.double
-    ty: cython.double
-    Xx: cython.double
-    Xy: cython.double
-    Xz: cython.double
-    ox: cython.double
-    oy: cython.double
-    oz: cython.double
-    dtot: cython.double
-    num_used: cython.int
-    px: cython.double
-    py: cython.double
-    pz: cython.double
-    v1x: cython.double
-    v1y: cython.double
-    v1z: cython.double
-    d1x: cython.double
-    d1y: cython.double
-    d1z: cython.double
-    v2x: cython.double
-    v2y: cython.double
-    v2z: cython.double
-    d2x: cython.double
-    d2y: cython.double
-    d2z: cython.double
-    sp0: cython.double
-    sp1: cython.double
-    sp2: cython.double
-    pb0: cython.double
-    pb1: cython.double
-    pb2: cython.double
-    scale: cython.double
-    dist: cython.double
-    mx: cython.double
-    my: cython.double
-    mz: cython.double
-    t0: cython.double
-    t1: cython.double
-    t2: cython.double
-    s1: cython.double
-    on1x: cython.double
-    on1y: cython.double
-    on1z: cython.double
-    s2: cython.double
-    on2x: cython.double
-    on2y: cython.double
-    on2z: cython.double
-    ddx: cython.double
-    ddy: cython.double
-    ddz: cython.double
-    verts_x: cython.double[:]
-    verts_y: cython.double[:]
-    verts_z: cython.double[:]
-    dirs_x: cython.double[:]
-    dirs_y: cython.double[:]
-    dirs_z: cython.double[:]
-    valid: cython.int[:]
-    _vi: cython.int
-    with cython.gil:
-        _verts_x_buf = np.zeros(8, dtype=np.float64)
-        verts_x = _verts_x_buf
-        _verts_y_buf = np.zeros(8, dtype=np.float64)
-        verts_y = _verts_y_buf
-        _verts_z_buf = np.zeros(8, dtype=np.float64)
-        verts_z = _verts_z_buf
-        _dirs_x_buf = np.zeros(8, dtype=np.float64)
-        dirs_x = _dirs_x_buf
-        _dirs_y_buf = np.zeros(8, dtype=np.float64)
-        dirs_y = _dirs_y_buf
-        _dirs_z_buf = np.zeros(8, dtype=np.float64)
-        dirs_z = _dirs_z_buf
-        _valid_buf = np.zeros(8, dtype=np.int32)
-        valid = _valid_buf
-
-    for _vi in range(8):
-        valid[_vi] = 0
-
-    for cam in range(num_cams):
-        tx = targets[cam, 0]
-        ty = targets[cam, 1]
-        if tx == COORD_UNUSED:
-            continue
-        _ray_tracing_out(tx, ty, cal_arr[cam], scratch_ray)
-        verts_x[cam] = scratch_ray[0]
-        verts_y[cam] = scratch_ray[1]
-        verts_z[cam] = scratch_ray[2]
-        dirs_x[cam] = scratch_ray[3]
-        dirs_y[cam] = scratch_ray[4]
-        dirs_z[cam] = scratch_ray[5]
-        valid[cam] = 1
-
-    dtot = 0.0
-    num_used = 0
-    px = 0.0
-    py = 0.0
-    pz = 0.0
-
-    for cam in range(num_cams):
-        if valid[cam] == 0:
-            continue
-        for pair in range(cam + 1, num_cams):
-            if valid[pair] == 0:
-                continue
-
-            v1x = verts_x[cam]
-            v1y = verts_y[cam]
-            v1z = verts_z[cam]
-            d1x = dirs_x[cam]
-            d1y = dirs_y[cam]
-            d1z = dirs_z[cam]
-            v2x = verts_x[pair]
-            v2y = verts_y[pair]
-            v2z = verts_z[pair]
-            d2x = dirs_x[pair]
-            d2y = dirs_y[pair]
-            d2z = dirs_z[pair]
-
-            sp0 = v2x - v1x
-            sp1 = v2y - v1y
-            sp2 = v2z - v1z
-
-            pb0 = d1y * d2z - d1z * d2y
-            pb1 = d1z * d2x - d1x * d2z
-            pb2 = d1x * d2y - d1y * d2x
-            scale = pb0 * pb0 + pb1 * pb1 + pb2 * pb2
-
-            if scale < 1e-20:
-                dist = c_sqrt(sp0 * sp0 + sp1 * sp1 + sp2 * sp2)
-                mx = (v1x + v2x) * 0.5
-                my = (v1y + v2y) * 0.5
-                mz = (v1z + v2z) * 0.5
-            else:
-                t0 = sp1 * d2z - sp2 * d2y
-                t1 = sp2 * d2x - sp0 * d2z
-                t2 = sp0 * d2y - sp1 * d2x
-                s1 = (pb0 * t0 + pb1 * t1 + pb2 * t2) / scale
-                on1x = v1x + d1x * s1
-                on1y = v1y + d1y * s1
-                on1z = v1z + d1z * s1
-
-                t0 = sp1 * d1z - sp2 * d1y
-                t1 = sp2 * d1x - sp0 * d1z
-                t2 = sp0 * d1y - sp1 * d1x
-                s2 = (pb0 * t0 + pb1 * t1 + pb2 * t2) / scale
-                on2x = v2x + d2x * s2
-                on2y = v2y + d2y * s2
-                on2z = v2z + d2z * s2
-
-                ddx = on1x - on2x
-                ddy = on1y - on2y
-                ddz = on1z - on2z
-                dist = c_sqrt(ddx * ddx + ddy * ddy + ddz * ddz)
-                mx = (on1x + on2x) * 0.5
-                my = (on1y + on2y) * 0.5
-                mz = (on1z + on2z) * 0.5
-
-            num_used += 1
-            dtot += dist
-            px += mx
-            py += my
-            pz += mz
-
-    if num_used > 0:
-        inv = 1.0 / num_used
-        out[0] = px * inv
-        out[1] = py * inv
-        out[2] = pz * inv
-        return dtot * inv
-    else:
-        out[0] = 0.0
-        out[1] = 0.0
-        out[2] = 0.0
-        return 0.0
-
-
-@cython.ccall
-def point_position_fast(
-    targets: cython.double[:, ::1], num_cams: cython.int, cal_arr: cython.double[:, ::1]
-):
-    """Compute 3D position from multiple camera rays.
-
-    Returns:
-        (pos, avg_dist) — (3,) float64 position and average ray distance.
-    """
-    pos = np.zeros(3, dtype=np.float64)
-    pos_mv: cython.double[:] = pos
-    scratch_ray = np.zeros(6, dtype=np.float64)
-    dtot = _point_position_out(targets, num_cams, cal_arr, pos_mv, scratch_ray)
-    return pos, dtot
-
-
-@cython.ccall
-def pixel_to_metric_fast(
-    x_pixel: cython.double,
-    y_pixel: cython.double,
-    imx: cython.int,
-    imy: cython.int,
-    pix_x: cython.double,
-    pix_y: cython.double,
-    chfield: cython.int,
-):
-    """Convert pixel to metric coordinates."""
-    yp: cython.double
-    x_metric: cython.double
-    y_metric: cython.double
-    yp = y_pixel
-    if chfield == 1:
-        yp = 2.0 * yp + 1.0
-    elif chfield == 2:
-        yp = 2.0 * yp
-    x_metric = (x_pixel - imx * 0.5) * pix_x
-    y_metric = (imy * 0.5 - yp) * pix_y
-    return x_metric, y_metric
-
-
-@cython.ccall
-@cython.inline
-@cython.cdivision(True)
-@cython.profile(False)
-@cython.nogil
-def _pixel_to_metric_out(
-    x_pixel: cython.double,
-    y_pixel: cython.double,
-    imx: cython.int,
-    imy: cython.int,
-    pix_x: cython.double,
-    pix_y: cython.double,
-    chfield: cython.int,
-    out: cython.double[:],
-) -> cython.int:
-    """Write pixel-to-metric coords to out[0], out[1]."""
-    yp: cython.double
-    yp = y_pixel
-    if chfield == 1:
-        yp = 2.0 * yp + 1.0
-    elif chfield == 2:
-        yp = 2.0 * yp
-    out[0] = (x_pixel - imx * 0.5) * pix_x
-    out[1] = (imy * 0.5 - yp) * pix_y
-    return 0
-
-
-@cython.ccall
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def dist_to_flat_fast(
-    dist_x: cython.double,
-    dist_y: cython.double,
-    xh: cython.double,
-    yh: cython.double,
-    k1: cython.double,
-    k2: cython.double,
-    k3: cython.double,
-    p1: cython.double,
-    p2: cython.double,
-    scx: cython.double,
-    she: cython.double,
-    tol: cython.double,
-):
-    """Inverse Brown distortion."""
-    r_init: cython.double
-    sin_she: cython.double
-    cos_she: cython.double
-    inv_scx: cython.double
-    xq: cython.double
-    yq: cython.double
-    _: cython.int
-    r2: cython.double
-    r4: cython.double
-    r6: cython.double
-    radial_factor: cython.double
-    dx: cython.double
-    dy: cython.double
-    xq_new: cython.double
-    yq_new: cython.double
-    dx_change: cython.double
-    dy_change: cython.double
-    r_init = c_sqrt(dist_x * dist_x + dist_y * dist_y)
-    if r_init < 1e-10:
-        return -xh, -yh
-
-    sin_she = c_sin(she)
-    cos_she = c_cos(she)
-    inv_scx = 1.0 / scx
-
-    xq = (dist_x + dist_y * sin_she) * inv_scx
-    yq = dist_y / cos_she
-
-    for _ in range(50):
-        r2 = xq * xq + yq * yq
-        r4 = r2 * r2
-        r6 = r4 * r2
-
-        radial_factor = k1 * r2 + k2 * r4 + k3 * r6
-
-        dx = xq * radial_factor + p1 * (r2 + 2.0 * xq * xq) + 2.0 * p2 * xq * yq
-        dy = yq * radial_factor + p2 * (r2 + 2.0 * yq * yq) + 2.0 * p1 * xq * yq
-
-        xq_new = (dist_x + dist_y * sin_she) * inv_scx - dx
-        yq_new = dist_y / cos_she - dy
-
-        dx_change = xq_new - xq
-        dy_change = yq_new - yq
-
-        xq += 0.5 * dx_change
-        yq += 0.5 * dy_change
-
-        if c_sqrt(dx_change * dx_change + dy_change * dy_change) < tol:
-            break
-
-    return xq - xh, yq - yh
-
-
-@cython.cfunc
-@cython.inline
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.cdivision(True)
-@cython.profile(False)
-@cython.nogil
-def _dist_to_flat_out(
-    dist_x: cython.double,
-    dist_y: cython.double,
-    xh: cython.double,
-    yh: cython.double,
-    k1: cython.double,
-    k2: cython.double,
-    k3: cython.double,
-    p1: cython.double,
-    p2: cython.double,
-    scx: cython.double,
-    she: cython.double,
-    tol: cython.double,
-    out: cython.double[:],
-) -> cython.int:
-    """Write dist-to-flat coords to out[0], out[1]."""
-    r_init: cython.double = c_sqrt(dist_x * dist_x + dist_y * dist_y)
-    if r_init < 1e-10:
-        out[0] = -xh
-        out[1] = -yh
-        return 0
-    sin_she: cython.double = c_sin(she)
-    cos_she: cython.double = c_cos(she)
-    inv_scx: cython.double = 1.0 / scx
-    xq: cython.double = (dist_x + dist_y * sin_she) * inv_scx
-    yq: cython.double = dist_y / cos_she
-    _: cython.int
-    r2: cython.double
-    r4: cython.double
-    r6: cython.double
-    radial_factor: cython.double
-    dx: cython.double
-    dy: cython.double
-    xq_new: cython.double
-    yq_new: cython.double
-    dx_change: cython.double
-    dy_change: cython.double
-    for _ in range(50):
-        r2 = xq * xq + yq * yq
-        r4 = r2 * r2
-        r6 = r4 * r2
-        radial_factor = k1 * r2 + k2 * r4 + k3 * r6
-        dx = xq * radial_factor + p1 * (r2 + 2.0 * xq * xq) + 2.0 * p2 * xq * yq
-        dy = yq * radial_factor + p2 * (r2 + 2.0 * yq * yq) + 2.0 * p1 * xq * yq
-        xq_new = (dist_x + dist_y * sin_she) * inv_scx - dx
-        yq_new = dist_y / cos_she - dy
-        dx_change = xq_new - xq
-        dy_change = yq_new - yq
-        xq += 0.5 * dx_change
-        yq += 0.5 * dy_change
-        if c_sqrt(dx_change * dx_change + dy_change * dy_change) < tol:
-            break
-    out[0] = xq - xh
-    out[1] = yq - yh
-    return 0
-
-
 @cython.ccall
 @cython.boundscheck(False)
 @cython.wraparound(False)
@@ -875,204 +201,6 @@ def assess_new_position_fast(
     return targ_pos, cand_inds, valid_cams
 
 
-@cython.ccall
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.nogil
-def _candsearch_in_pix_rest_nogil(
-    targ_x: cython.double[:],
-    targ_y: cython.double[:],
-    targ_tnr: cython.int[:],
-    num_targets: cython.int,
-    cent_x: cython.double,
-    cent_y: cython.double,
-    dl: cython.double,
-    dr: cython.double,
-    du: cython.double,
-    dd: cython.double,
-    imx: cython.double,
-    imy: cython.double,
-    tr_unused: cython.int,
-) -> cython.int:
-    """Find closest unused candidate GIL-free."""
-    xmin: cython.double
-    xmax: cython.double
-    ymin: cython.double
-    ymax: cython.double
-    best: cython.int
-    dmin: cython.double
-    j0: cython.int
-    dj: cython.int
-    j: cython.int
-    ty: cython.double
-    tx: cython.double
-    dx: cython.double
-    dy: cython.double
-    d: cython.double
-    xmin = cent_x - dl
-    xmax = cent_x + dr
-    ymin = cent_y - du
-    ymax = cent_y + dd
-
-    if xmin < 0.0:
-        xmin = 0.0
-    if xmax > imx:
-        xmax = imx
-    if ymin < 0.0:
-        ymin = 0.0
-    if ymax > imy:
-        ymax = imy
-
-    best = tr_unused
-    dmin = 1e20
-
-    if not (0.0 <= cent_x <= imx and 0.0 <= cent_y <= imy):
-        return best
-
-    j0 = num_targets // 2
-    dj = num_targets // 4
-    while dj > 1:
-        if targ_y[j0] < ymin:
-            j0 += dj
-        else:
-            j0 -= dj
-        dj //= 2
-
-    j0 -= 12
-    if j0 < 0:
-        j0 = 0
-
-    for j in range(j0, num_targets):
-        ty = targ_y[j]
-        if targ_tnr[j] == tr_unused:
-            if ty > ymax:
-                break
-            tx = targ_x[j]
-            if tx > xmin and tx < xmax and ty > ymin and ty < ymax:
-                dx = cent_x - tx
-                dy = cent_y - ty
-                d = c_sqrt(dx * dx + dy * dy)
-                if d < dmin:
-                    dmin = d
-                    best = j
-
-    return best
-
-
-@cython.ccall
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.nogil
-def assess_new_position_fast_nogil(
-    pos: cython.double[:],
-    num_cams: cython.int,
-    add_part: cython.double,
-    cal_arr: cython.double[:, ::1],
-    mo_arr: cython.double[:, ::1],
-    mnr_arr: cython.int[:],
-    mnz_arr: cython.int[:],
-    mrw_arr: cython.double[:],
-    targ_x: cython.double[:, ::1],
-    targ_y: cython.double[:, ::1],
-    targ_tnr: cython.int[:, ::1],
-    num_targets: cython.int[:],
-    imx_half: cython.double,
-    imy_half: cython.double,
-    inv_pix_x: cython.double,
-    inv_pix_y: cython.double,
-    chfield: cython.int,
-    imx: cython.int,
-    imy: cython.int,
-    pix_x: cython.double,
-    pix_y: cython.double,
-    flatten_tol: cython.double,
-    tr_unused: cython.int,
-    coord_unused: cython.double,
-    proj_x: cython.double[:],
-    proj_y: cython.double[:],
-    targ_pos_out: cython.double[:, :],
-    cand_inds_out: cython.int[:],
-    scratch: cython.double[:],
-) -> cython.int:
-    """Assess new position GIL-free. Assumes use_proj=True."""
-    cam: cython.int
-    valid_cams: cython.int
-    best: cython.int
-    px: cython.double
-    py: cython.double
-    mx: cython.double
-    my: cython.double
-
-    for cam in range(num_cams):
-        cand_inds_out[cam] = tr_unused
-        targ_pos_out[cam, 0] = coord_unused
-        targ_pos_out[cam, 1] = coord_unused
-
-    for cam in range(num_cams):
-        px = proj_x[cam]
-        py = proj_y[cam]
-
-        best = _candsearch_in_pix_rest_nogil(
-            targ_x[cam],
-            targ_y[cam],
-            targ_tnr[cam],
-            num_targets[cam],
-            px,
-            py,
-            add_part,
-            add_part,
-            add_part,
-            add_part,
-            imx,
-            imy,
-            tr_unused,
-        )
-
-        if best != tr_unused:
-            cand_inds_out[cam] = best
-            targ_pos_out[cam, 0] = targ_x[cam, best]
-            targ_pos_out[cam, 1] = targ_y[cam, best]
-
-    valid_cams = 0
-    for cam in range(num_cams):
-        if targ_pos_out[cam, 0] != coord_unused:
-            _pixel_to_metric_out(
-                targ_pos_out[cam, 0],
-                targ_pos_out[cam, 1],
-                imx,
-                imy,
-                pix_x,
-                pix_y,
-                chfield,
-                scratch,
-            )
-            mx = scratch[0]
-            my = scratch[1]
-
-            cal = cal_arr[cam]
-            _dist_to_flat_out(
-                mx,
-                my,
-                cal[13],
-                cal[14],
-                cal[24],
-                cal[25],
-                cal[26],
-                cal[27],
-                cal[28],
-                cal[29],
-                cal[30],
-                flatten_tol,
-                scratch,
-            )
-
-            targ_pos_out[cam, 0] = scratch[0]
-            targ_pos_out[cam, 1] = scratch[1]
-            valid_cams += 1
-
-    return valid_cams
-
-
 POSI_K = 80
 MAX_CANDS_K = 32
 TR_UNUSED_K = -1
@@ -1082,368 +210,3 @@ def assess_new_position_fast_nogil(
 COORD_UNUSED_K = -1e10
 ADD_PART_K = 3.0
 
-
-@cython.ccall
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def metric_to_pixel_fast(
-    x_metric: cython.double,
-    y_metric: cython.double,
-    imx: cython.int,
-    imy: cython.int,
-    pix_x: cython.double,
-    pix_y: cython.double,
-    chfield: cython.int,
-):
-    """Convert metric to pixel coordinates."""
-    x_pixel: cython.double
-    y_pixel: cython.double
-    x_pixel = x_metric / pix_x + imx * 0.5
-    y_pixel = imy * 0.5 - y_metric / pix_y
-    if chfield == 1:
-        y_pixel = (y_pixel - 1.0) * 0.5
-    elif chfield == 2:
-        y_pixel = y_pixel * 0.5
-    return x_pixel, y_pixel
-
-
-@cython.ccall
-@cython.inline
-@cython.cdivision(True)
-@cython.profile(False)
-def _metric_to_pixel_out(
-    x_metric: cython.double,
-    y_metric: cython.double,
-    imx: cython.int,
-    imy: cython.int,
-    pix_x: cython.double,
-    pix_y: cython.double,
-    chfield: cython.int,
-    out: cython.double[:],
-):
-    """Write metric-to-pixel coords to out[0], out[1]."""
-    x_pixel: cython.double = x_metric / pix_x + imx * 0.5
-    y_pixel: cython.double = imy * 0.5 - y_metric / pix_y
-    if chfield == 1:
-        y_pixel = (y_pixel - 1.0) * 0.5
-    elif chfield == 2:
-        y_pixel = y_pixel * 0.5
-    out[0] = x_pixel
-    out[1] = y_pixel
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def _flat_image_coord_fast(
-    pos: cython.double[:],
-    cal: cython.double[:],
-    mmlut_data: cython.double[:],
-    mmlut_origin: cython.double[:],
-    mmlut_nr: cython.int,
-    mmlut_nz: cython.int,
-    mmlut_rw: cython.double,
-):
-    """Project 3D to flat metric image coordinates.
-
-    Returns (x, y) without distortion or pixel conversion.
-    """
-    pos0: cython.double
-    pos1: cython.double
-    pos2: cython.double
-    ext_x0: cython.double
-    ext_y0: cython.double
-    ext_z0: cython.double
-    dm00: cython.double
-    dm10: cython.double
-    dm20: cython.double
-    dm01: cython.double
-    dm11: cython.double
-    dm21: cython.double
-    dm02: cython.double
-    dm12: cython.double
-    dm22: cython.double
-    int_cc: cython.double
-    gx: cython.double
-    gy: cython.double
-    gz: cython.double
-    inv_dog: cython.double
-    mm_n1: cython.double
-    mm_n2_0: cython.double
-    mm_n3: cython.double
-    mm_d0: cython.double
-    dot_cam: cython.double
-    dist_o_glas: cython.double
-    dist_cam_glas: cython.double
-    dot_pos: cython.double
-    dist_point_glas: cython.double
-    s_cam: cython.double
-    cc_x: cython.double
-    cc_y: cython.double
-    cc_z: cython.double
-    s_pt: cython.double
-    cp_x: cython.double
-    cp_y: cython.double
-    cp_z: cython.double
-    ext_t_z0: cython.double
-    s_d: cython.double
-    ag_x: cython.double
-    ag_y: cython.double
-    ag_z: cython.double
-    tmp_x: cython.double
-    tmp_y: cython.double
-    tmp_z: cython.double
-    pos_t_0: cython.double
-    pos_t_2: cython.double
-    radial_shift: cython.double
-    has_mmlut: cython.bint
-    tx: cython.double
-    ty: cython.double
-    tz: cython.double
-    sz: cython.double
-    iz: cython.int
-    R: cython.double
-    sr: cython.double
-    ir: cython.int
-    v0: cython.int
-    v3: cython.int
-    mmf: cython.double
-    X_t: cython.double
-    s_z: cython.double
-    bx: cython.double
-    by: cython.double
-    bz: cython.double
-    s_x: cython.double
-    dx: cython.double
-    dy: cython.double
-    dz: cython.double
-    deno: cython.double
-    x: cython.double
-    y: cython.double
-    pos0 = pos[0]
-    pos1 = pos[1]
-    pos2 = pos[2]
-
-    ext_x0 = cal[0]
-    ext_y0 = cal[1]
-    ext_z0 = cal[2]
-    dm00 = cal[3]
-    dm10 = cal[4]
-    dm20 = cal[5]
-    dm01 = cal[6]
-    dm11 = cal[7]
-    dm21 = cal[8]
-    dm02 = cal[9]
-    dm12 = cal[10]
-    dm22 = cal[11]
-    int_cc = cal[12]
-    gx = cal[15]
-    gy = cal[16]
-    gz = cal[17]
-    inv_dog = cal[19]
-    mm_n1 = cal[20]
-    mm_n2_0 = cal[21]
-    mm_n3 = cal[22]
-    mm_d0 = cal[23]
-
-    dot_cam = ext_x0 * gx + ext_y0 * gy + ext_z0 * gz
-    dist_o_glas = cal[18]
-    dist_cam_glas = dot_cam * inv_dog - dist_o_glas - mm_d0
-
-    dot_pos = pos0 * gx + pos1 * gy + pos2 * gz
-    dist_point_glas = dot_pos * inv_dog - dist_o_glas
-
-    s_cam = dist_cam_glas * inv_dog
-    cc_x = ext_x0 - gx * s_cam
-    cc_y = ext_y0 - gy * s_cam
-    cc_z = ext_z0 - gz * s_cam
-
-    s_pt = dist_point_glas * inv_dog
-    cp_x = pos0 - gx * s_pt
-    cp_y = pos1 - gy * s_pt
-    cp_z = pos2 - gz * s_pt
-
-    ext_t_z0 = dist_cam_glas + mm_d0
-
-    s_d = mm_d0 * inv_dog
-    ag_x = cc_x - gx * s_d
-    ag_y = cc_y - gy * s_d
-    ag_z = cc_z - gz * s_d
-    tmp_x = cp_x - ag_x
-    tmp_y = cp_y - ag_y
-    tmp_z = cp_z - ag_z
-
-    pos_t_0 = c_sqrt(tmp_x * tmp_x + tmp_y * tmp_y + tmp_z * tmp_z)
-    pos_t_2 = dist_point_glas
-
-    radial_shift = 1.0
-    has_mmlut = len(mmlut_data) > 0
-    if has_mmlut:
-        tx = pos_t_0 - mmlut_origin[0]
-        ty = -mmlut_origin[1]
-        tz = pos_t_2 - mmlut_origin[2]
-        sz = tz / mmlut_rw
-        iz = int(sz)
-        sz -= iz
-        R = c_sqrt(tx * tx + ty * ty)
-        sr = R / mmlut_rw
-        ir = int(sr)
-        sr -= ir
-        if ir <= mmlut_nr and iz >= 0 and iz <= mmlut_nz:
-            v0 = ir * mmlut_nz + iz
-            v3 = v0 + mmlut_nz + 1
-            if v0 >= 0 and v3 <= mmlut_nr * mmlut_nz:
-                mmf = (
-                    mmlut_data[v0] * (1.0 - sr) * (1.0 - sz)
-                    + mmlut_data[v0 + 1] * (1.0 - sr) * sz
-                    + mmlut_data[v0 + mmlut_nz] * sr * (1.0 - sz)
-                    + mmlut_data[v3] * sr * sz
-                )
-                if mmf > 0.0:
-                    radial_shift = mmf
-    if radial_shift == 1.0:
-        radial_shift = _multimed_r_nlay_1layer(
-            pos_t_0,
-            0.0,
-            pos_t_2,
-            0.0,
-            0.0,
-            ext_t_z0,
-            mm_n1,
-            mm_n2_0,
-            mm_n3,
-            mm_d0,
-        )
-    X_t = pos_t_0 * radial_shift
-
-    s_z = -pos_t_2 * inv_dog
-    bx = ag_x - gx * s_z
-    by = ag_y - gy * s_z
-    bz = ag_z - gz * s_z
-    if pos_t_0 > 0.0:
-        s_x = -X_t / pos_t_0
-        bx -= tmp_x * s_x
-        by -= tmp_y * s_x
-        bz -= tmp_z * s_x
-
-    dx = bx - ext_x0
-    dy = by - ext_y0
-    dz = bz - ext_z0
-    deno = dm02 * dx + dm12 * dy + dm22 * dz
-    x = -int_cc * (dm00 * dx + dm10 * dy + dm20 * dz) / deno
-    y = -int_cc * (dm01 * dx + dm11 * dy + dm21 * dz) / deno
-
-    return x, y
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def _img_coord_fast(
-    pos: cython.double[:],
-    cal: cython.double[:],
-    mmlut_data: cython.double[:],
-    mmlut_origin: cython.double[:],
-    mmlut_nr: cython.int,
-    mmlut_nz: cython.int,
-    mmlut_rw: cython.double,
-):
-    """Project 3D to distorted metric image coordinates."""
-    xh: cython.double
-    yh: cython.double
-    k1: cython.double
-    k2: cython.double
-    k3: cython.double
-    p1: cython.double
-    p2: cython.double
-    scx: cython.double
-    she: cython.double
-    x: cython.double
-    y: cython.double
-    r: cython.double
-    r2: cython.double
-    r4: cython.double
-    radial_factor: cython.double
-    xd: cython.double
-    yd: cython.double
-    sin_she: cython.double
-    cos_she: cython.double
-    x_dist: cython.double
-    y_dist: cython.double
-    x, y = _flat_image_coord_fast(
-        pos, cal, mmlut_data, mmlut_origin, mmlut_nr, mmlut_nz, mmlut_rw
-    )
-
-    xh = cal[13]
-    yh = cal[14]
-    k1 = cal[24]
-    k2 = cal[25]
-    k3 = cal[26]
-    p1 = cal[27]
-    p2 = cal[28]
-    scx = cal[29]
-    she = cal[30]
-
-    x += xh
-    y += yh
-    r = c_sqrt(x * x + y * y)
-    if r < 1e-10:
-        return 0.0, 0.0
-
-    r2 = r * r
-    r4 = r2 * r2
-    radial_factor = 1.0 + k1 * r2 + k2 * r4 + k3 * r4 * r2
-    xd = x * radial_factor + p1 * (r2 + 2.0 * x * x) + 2.0 * p2 * x * y
-    yd = y * radial_factor + p2 * (r2 + 2.0 * y * y) + 2.0 * p1 * x * y
-    sin_she = c_sin(she)
-    cos_she = c_cos(she)
-    x_dist = scx * (xd - sin_she * yd)
-    y_dist = scx * cos_she * yd
-
-    return x_dist, y_dist
-
-
-@cython.ccall
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def img_coord_batch_fast(
-    positions: cython.double[:, ::1],
-    cal: cython.double[:],
-    mmlut_data: cython.double[:],
-    mmlut_origin: cython.double[:],
-    mmlut_nr: cython.int,
-    mmlut_nz: cython.int,
-    mmlut_rw: cython.double,
-):
-    """Project N 3D positions to distorted metric coords."""
-    n: cython.Py_ssize_t
-    i: cython.Py_ssize_t
-    n = positions.shape[0]
-    result = np.empty((n, 2), dtype=np.float64)
-    for i in range(n):
-        result[i, 0], result[i, 1] = _img_coord_fast(
-            positions[i], cal, mmlut_data, mmlut_origin, mmlut_nr, mmlut_nz, mmlut_rw
-        )
-    return result
-
-
-@cython.ccall
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def flat_image_coord_batch_fast(
-    positions: cython.double[:, ::1],
-    cal: cython.double[:],
-    mmlut_data: cython.double[:],
-    mmlut_origin: cython.double[:],
-    mmlut_nr: cython.int,
-    mmlut_nz: cython.int,
-    mmlut_rw: cython.double,
-):
-    """Project N 3D positions to flat metric coords."""
-    n: cython.Py_ssize_t
-    i: cython.Py_ssize_t
-    n = positions.shape[0]
-    result = np.empty((n, 2), dtype=np.float64)
-    for i in range(n):
-        result[i, 0], result[i, 1] = _flat_image_coord_fast(
-            positions[i], cal, mmlut_data, mmlut_origin, mmlut_nr, mmlut_nz, mmlut_rw
-        )
-    return result
diff --git a/src/openptv2/plugins/two_phase_tracking.py b/src/openptv2/plugins/two_phase_tracking.py
index b2c8c965..af4b0ea9 100644
--- a/src/openptv2/plugins/two_phase_tracking.py
+++ b/src/openptv2/plugins/two_phase_tracking.py
@@ -36,12 +36,49 @@ class TwoPhaseTrackerConfig:
     leaf_weight : float
         Weight for 2D leaf distances in the cost matrix. If 0, falls back
         to pure 3D matching.
+    use_velocity : bool
+        Predict each track forward with its constant-velocity estimate and
+        match predictions (not current positions) against new detections.
+        Required to cross steady trajectories correctly; without it every
+        X-crossing resolves as a bounce. Default True.
+    cost_mode : str
+        "projected": leaf costs are evaluated at re-projected predicted
+        positions (needs project_fn) -- the benchmarked fix for
+        maneuver-at-crossing scenes. "3d": cost is the 3D distance between
+        prediction and candidate. Falls back to "3d" when no project_fn
+        is available.
+    allow_shared : bool
+        Prototype (shared-observation): in a contested component with more
+        tracks than candidates (detector undercount = occlusion), let the
+        losing tracks SHARE the winner's detection instead of dying.
+        Shared links update position but never velocity (each track's speed
+        comes only from its own points). Default False (legacy behaviour).
+    max_shared : int
+        Maximum consecutive shared frames per track before it must match
+        alone again. Default 2.
+    share_tol : float | None
+        Maximum edge cost (same units as the cost matrix) for a shared
+        claim. Sharing without it hijacks strangers: any unassigned track
+        inside the gate would co-opt a foreign detection (seen live: a
+        gap-stranded track shared two unrelated detections). None disables
+        the gate. Default 1.0.
+    max_group_size : int
+        Groups bigger than this skip the cubic Hungarian and fall back to
+        greedy claiming inside the group (sharing still applies). At
+        production density the frame percolates into giant components;
+        without the cap one frame stalls the run. Default 128.
     """
 
     v_max: float = 5.0
     max_gap: int = 2
     dt: float = 1.0
     leaf_weight: float = 1.0
+    use_velocity: bool = True
+    cost_mode: str = "projected"
+    allow_shared: bool = False
+    max_shared: int = 2
+    share_tol: float | None = 1.0
+    max_group_size: int = 128
 
 
 def _match_two_phase_frame(
@@ -53,25 +90,42 @@ def _match_two_phase_frame(
     p1: np.ndarray,
     radius: float,
     leaf_weight: float = 1.0,
-) -> set[tuple[int, int]]:
+    cost_mode: str = "projected",
+    allow_shared: bool = False,
+    share_tol: float | None = None,
+    max_group_size: int = 128,
+) -> tuple[set[tuple[int, int]], set[tuple[int, int]]]:
     """Two-phase frame-to-frame matching: 3D search + 2D ranking.
 
     Parameters
     ----------
-    pts0, pts1 : (N, 3) and (M, 3) — 3D positions in mm
+    pts0, pts1 : (N, 3) and (M, 3) — 3D positions in mm. With velocity
+        prediction enabled these are PREDICTED positions; pass the
+        re-projected pixel positions as xy0 so costs are evaluated at the
+        prediction (stale appearance nulls the prediction -- the bounce
+        bias returns).
     xy0, xy1 : (N, D) and (M, D) — flattened 2D leaf features
     p0, p1 : particle IDs for frame 0 and 1
     radius : float — 3D search radius in mm
     leaf_weight : float — weight for 2D distances in cost matrix
+    cost_mode : str — "projected" (2D leaf costs) or "3d" (3D distance
+        between pts0 and candidates; needs no calibration).
+    allow_shared : bool — prototype: in components with more predictors
+        than candidates, unassigned predictors share their best candidate
+        instead of going unmatched.
+    share_tol : float | None — maximum edge cost for a shared claim
+        (mutual-good-prediction gate). None disables the gate.
 
     Returns
     -------
     links : set of (pid0, pid1) pairs
+    shared : set of (pid0, pid1) pairs, subset of links, observed jointly
+        (empty unless allow_shared)
     """
     n_pred = len(pts0)
     n_cand = len(pts1)
     if n_pred == 0 or n_cand == 0:
-        return set()
+        return set(), set()
 
     # Phase 1: 3D KD-tree candidate search
     tree3d = cKDTree(pts1)
@@ -79,11 +133,13 @@ def _match_two_phase_frame(
 
     # Build edge list with 2D costs
     rows, cols, costs = [], [], []
+    use_leaves = (cost_mode == "projected" and leaf_weight > 0
+                  and xy0.shape[1] > 0)
     for pi in range(n_pred):
         cands = neighbours[pi]
         if len(cands) == 0:
             continue
-        if leaf_weight > 0 and xy0.shape[1] > 0:
+        if use_leaves:
             # 2D cost: mean Euclidean distance per camera, weighted by overlap count
             C = xy0.shape[1] // 2
             xy0_cam = xy0[pi].reshape(C, 2)
@@ -110,7 +166,7 @@ def _match_two_phase_frame(
                 costs.append(np.linalg.norm(pts0[pi] - pts1[ci]))
 
     if len(rows) == 0:
-        return set()
+        return set(), set()
 
     rows = np.array(rows)
     cols = np.array(cols)
@@ -125,6 +181,7 @@ def _match_two_phase_frame(
     n_comp, labels = connected_components(graph, directed=False)
 
     links = set()
+    shared: set[tuple[int, int]] = set()
     edge_comp = labels[rows]
     comp_edges = np.bincount(edge_comp, minlength=n_comp)
 
@@ -135,6 +192,8 @@ def _match_two_phase_frame(
 
     # Non-trivial: small dense Hungarian per component
     rest = np.flatnonzero(~trivial)
+    assigned_rows: set[int] = set()
+    assigned_cols: set[int] = set()
     if len(rest):
         rest = rest[np.argsort(edge_comp[rest], kind="stable")]
         splits = np.flatnonzero(np.diff(edge_comp[rest])) + 1
@@ -144,6 +203,29 @@ def _match_two_phase_frame(
             c_costs = costs[group].tolist()
             uniq_r = sorted(set(c_rows))
             uniq_c = sorted(set(c_cols))
+            if len(uniq_r) + len(uniq_c) > max_group_size:
+                # Production-density percolation: cubic Hungarian would stall
+                # the frame. Greedy inside the group, sharing still applies.
+                order = np.argsort(np.array(c_costs), kind="stable")
+                for k in order.tolist():
+                    r, c = c_rows[k], c_cols[k]
+                    if r not in assigned_rows and c not in assigned_cols:
+                        links.add((int(p0[r]), int(p1[c])))
+                        assigned_rows.add(r)
+                        assigned_cols.add(c)
+                if allow_shared:
+                    for r in uniq_r:
+                        if r in assigned_rows:
+                            continue
+                        best_c, best_d = None, np.inf
+                        for k, (rr, cc) in enumerate(zip(c_rows, c_cols)):
+                            if rr == r and c_costs[k] < best_d:
+                                best_d, best_c = c_costs[k], cc
+                        if best_c is not None and \
+                                (share_tol is None or best_d < share_tol):
+                            links.add((int(p0[r]), int(p1[best_c])))
+                            shared.add((int(p0[r]), int(p1[best_c])))
+                continue
             r_local = {v: i for i, v in enumerate(uniq_r)}
             c_local = {v: i for i, v in enumerate(uniq_c)}
             max_cost = max(c_costs) if c_costs else 1.0
@@ -153,10 +235,38 @@ def _match_two_phase_frame(
                 sub[r_local[rr], c_local[cc]] = dd
             r_ind, c_ind = linear_sum_assignment(sub)
             real = sub[r_ind, c_ind] < sentinel
+            group_winners = []
             for r_i, c_i in zip(r_ind[real], c_ind[real]):
                 links.add((int(p0[uniq_r[r_i]]), int(p1[uniq_c[c_i]])))
-
-    return links
+                assigned_rows.add(int(uniq_r[r_i]))
+                assigned_cols.add(int(uniq_c[c_i]))
+                group_winners.append((int(p0[uniq_r[r_i]]),
+                                      int(p1[uniq_c[c_i]])))
+            # Prototype shared-observation: more predictors than candidates
+            # = detector undercount (occlusion). Unassigned predictors share
+            # their best candidate instead of dying. Tracked by the caller
+            # via the streak cap.
+            group_shared = []
+            if allow_shared and len(uniq_r) > len(uniq_c):
+                for r_i, r in enumerate(uniq_r):
+                    if r in assigned_rows:
+                        continue
+                    best_c, best_d = None, sentinel
+                    for c_i, c in enumerate(uniq_c):
+                        if sub[r_i, c_i] < best_d:
+                            best_d, best_c = sub[r_i, c_i], c
+                    if best_c is not None and best_d < sentinel and \
+                            (share_tol is None or best_d < share_tol):
+                        links.add((int(p0[r]), int(p1[best_c])))
+                        group_shared.append((int(p0[r]), int(p1[best_c])))
+            if group_shared:
+                # The winners' points in a sharing group are joint evidence
+                # too: nobody updates velocity from a merged point, or the
+                # winner predicts from poisoned history at separation.
+                shared.update(group_winners)
+                shared.update(group_shared)
+
+    return links, shared
 
 
 class TwoPhaseTracker:
@@ -175,9 +285,18 @@ def track_frames(
         self,
         frame_particles: list[np.ndarray],
         frame_leaves: list[np.ndarray] | None = None,
+        project_fn=None,
+        return_chains: bool = False,
     ) -> list[tuple[int, int, int, int]]:
         """Track particles across frames using two-phase matching.
 
+        Stateful: every live track carries a constant-velocity estimate.
+        Each step matches velocity PREDICTIONS (not current positions)
+        against new detections, so steady crossings resolve correctly.
+        New detections spawn zero-velocity tracks (cold start); tracks
+        unmatched for more than ``max_gap`` frames retire, so a particle
+        occluded for a frame is re-caught by gap-spanning prediction.
+
         Parameters
         ----------
         frame_particles : list of (N_i, 3) arrays
@@ -185,11 +304,22 @@ def track_frames(
         frame_leaves : list of (N_i, D) arrays, optional
             Flattened 2D leaf features per frame. If None, falls back to
             pure 3D matching.
+        project_fn : callable, optional
+            ``(N, 3) -> (N, D)`` mapping predicted 3D positions to leaf
+            features (re-projection through the camera models). Required
+            for ``cost_mode="projected"``; without it costs fall back to
+            3D distance (see ``cost_mode``).
 
         Returns
         -------
         links : list of (t0, pid0, t1, pid1) tuples
-            Frame-to-frame particle links (0-based time indices).
+            Frame-to-frame particle links (0-based time indices, row ids
+            within each frame's arrays). Gap-spanning links reference the
+            track's last seen frame/row.
+        (if return_chains) chains : list of dicts with keys tid, frames,
+            rows, pos, shared — per-track point histories; shared flags
+            mark jointly-observed points. Row-links alone cannot represent
+            sharing (one node, two owners), hence chains.
         """
         num_frames = len(frame_particles)
         if num_frames < 2:
@@ -198,28 +328,138 @@ def track_frames(
         if frame_leaves is None:
             frame_leaves = [np.zeros((len(p), 0)) for p in frame_particles]
 
+        cost_mode = self.cfg.cost_mode
+        if cost_mode == "projected" and project_fn is None:
+            cost_mode = "3d"
+
+        next_tid = 0
+        # tid -> dict(pos, vel, last_t, last_row, misses, shared_streak)
+        # shared_streak counts consecutive shared observations; velocity is
+        # NEVER updated from a shared point (each track's speed comes only
+        # from its own points).
+        tracks: dict[int, dict] = {}
+        # tid -> list of (frame, row, is_shared): full point history, kept
+        # for retired tracks too (row-links cannot represent sharing).
+        hist: dict[int, list[tuple[int, int, bool]]] = {}
+        # (frame_idx, row) -> tid, for emitting row-based links
+        loc2tid: dict[tuple[int, int], int] = {}
+        for i, p in enumerate(frame_particles[0]):
+            tracks[next_tid] = {
+                "pos": np.asarray(p, dtype=np.float64),
+                "vel": np.zeros(3),
+                "last_t": 0,
+                "last_row": i,
+                "misses": 0,
+                "shared_streak": 0,
+            }
+            loc2tid[(0, i)] = next_tid
+            hist[next_tid] = [(0, i, False)]
+            next_tid += 1
+
         all_links = []
-        for i in range(num_frames - 1):
-            t0, t1 = i, i + 1
-            pts0, pts1 = frame_particles[t0], frame_particles[t1]
-            lf0, lf1 = frame_leaves[t0], frame_leaves[t1]
-            p0 = np.arange(len(pts0), dtype=np.int32)
-            p1 = np.arange(len(pts1), dtype=np.int32)
-
-            links = _match_two_phase_frame(
-                pts0,
-                pts1,
-                lf0,
-                lf1,
-                p0,
-                p1,
+        for t in range(num_frames - 1):
+            pts1 = np.asarray(frame_particles[t + 1], dtype=np.float64)
+            lf1 = frame_leaves[t + 1]
+            n1 = len(pts1)
+
+            # Active tracks: seen within max_gap frames.
+            active = [tid for tid, tr in tracks.items()
+                      if t + 1 - tr["last_t"] <= self.cfg.max_gap]
+            if not active or n1 == 0:
+                pred_pts = np.zeros((0, 3))
+                pred_xy = np.zeros((0, lf1.shape[1] if n1 else 0))
+                tids: list[int] = []
+            else:
+                steps = np.array([t + 1 - tracks[tid]["last_t"]
+                                  for tid in active], dtype=np.float64)
+                if self.cfg.use_velocity:
+                    pred_pts = np.array(
+                        [tracks[tid]["pos"]
+                         + tracks[tid]["vel"] * steps[k] * self.cfg.dt
+                         for k, tid in enumerate(active)])
+                else:
+                    pred_pts = np.array([tracks[tid]["pos"] for tid in active])
+                if cost_mode == "projected":
+                    pred_xy = np.asarray(project_fn(pred_pts))
+                else:
+                    pred_xy = np.zeros((len(active), 0))
+                tids = active
+
+            got, got_shared = _match_two_phase_frame(
+                np.asarray(pred_pts, dtype=np.float64),
+                np.asarray(pts1, dtype=np.float64),
+                np.asarray(pred_xy, dtype=np.float64),
+                np.asarray(lf1, dtype=np.float64),
+                np.arange(len(tids), dtype=np.int32),
+                np.arange(n1, dtype=np.int32),
                 self.cfg.v_max,
                 self.cfg.leaf_weight,
+                cost_mode=cost_mode,
+                allow_shared=self.cfg.allow_shared,
+                share_tol=self.cfg.share_tol,
             )
-            for pid0, pid1 in links:
-                all_links.append((t0, pid0, t1, pid1))
-
-        return all_links
+            matched_det = set()
+            for ai, det in got:
+                tid = tids[int(ai)]
+                tr = tracks[tid]
+                gap = t + 1 - tr["last_t"]
+                old_pos = tr["pos"]
+                if (int(ai), int(det)) in got_shared and \
+                        tr.get("shared_streak", 0) < self.cfg.max_shared:
+                    # Shared observation: follow the point, keep own speed.
+                    tr["pos"] = pts1[det].copy()
+                    tr["shared_streak"] = tr.get("shared_streak", 0) + 1
+                    is_shared = True
+                else:
+                    tr["vel"] = (pts1[det] - old_pos) / (gap * self.cfg.dt)
+                    tr["pos"] = pts1[det].copy()
+                    tr["shared_streak"] = 0
+                    is_shared = False
+                all_links.append((tr["last_t"], tr["last_row"], t + 1, det))
+                tr["last_t"] = t + 1
+                tr["last_row"] = int(det)
+                tr["misses"] = 0
+                hist[tid].append((t + 1, int(det), is_shared))
+                loc2tid[(t + 1, int(det))] = tid
+                matched_det.add(int(det))
+
+            # Age every unseen track (including ones already outside the
+            # active window) and retire the exhausted.
+            for tid in list(tracks.keys()):
+                if tracks[tid]["last_t"] <= t:
+                    tracks[tid]["misses"] += 1
+                    if tracks[tid]["misses"] > self.cfg.max_gap:
+                        del tracks[tid]
+
+            # Cold start: unmatched detections become zero-velocity tracks.
+            for det in range(n1):
+                if det not in matched_det:
+                    tracks[next_tid] = {
+                        "pos": pts1[det].copy(),
+                        "vel": np.zeros(3),
+                        "last_t": t + 1,
+                        "last_row": det,
+                        "misses": 0,
+                        "shared_streak": 0,
+                    }
+                    loc2tid[(t + 1, det)] = next_tid
+                    hist[next_tid] = [(t + 1, det, False)]
+                    next_tid += 1
+
+        if not return_chains:
+            return all_links
+        chains = []
+        fp_arr = [np.asarray(p, dtype=np.float64) for p in frame_particles]
+        for tid, pts in hist.items():
+            if len(pts) == 0:
+                continue
+            chains.append({
+                "tid": tid,
+                "frames": [f for f, _, _ in pts],
+                "pos": np.array([fp_arr[f][r] for f, r, _ in pts]),
+                "shared": [s for _, _, s in pts],
+            })
+        return all_links, chains
 
 
 class Tracking:
@@ -234,6 +474,44 @@ def __init__(self, ptv=None, exp=None):
         self.ptv = ptv
         self.exp = exp
 
+    def _build_project_fn(self):
+        """Re-project predicted 3D positions to leaf pixels via exp cals.
+
+        Returns None when calibrations are unavailable; the tracker then
+        falls back to 3D-distance costs (see ``cost_mode``).
+        """
+        try:
+            cals = list(getattr(self.exp, "cals", None) or [])
+            cpar = getattr(self.exp, "cpar", None)
+            if not cals or cpar is None:
+                return None
+            mm = cpar.mm
+            imx, imy = float(cpar.imx), float(cpar.imy)
+            pix_x, pix_y = float(cpar.pix_x), float(cpar.pix_y)
+
+            from openptv2.algorithms.imgcoord import img_coord_batch
+
+            def project_fn(pred):
+                pred = np.asarray(pred, dtype=np.float64)
+                n = len(pred)
+                nc = len(cals)
+                xy = np.full((n, nc * 2), np.nan)
+                for i in range(n):
+                    for ci, cal in enumerate(cals):
+                        m = img_coord_batch(pred[i : i + 1], cal, mm)[0]
+                        xy[i, 2 * ci] = m[0] / pix_x + imx / 2
+                        xy[i, 2 * ci + 1] = imy / 2 - m[1] / pix_y
+                return np.nan_to_num(xy)
+
+            # Smoke-test on one point so a broken model fails here, not
+            # mid-run.
+            project_fn(np.zeros((1, 3)))
+            return project_fn
+        except Exception as exc:
+            print(f"TwoPhaseTracker: no projection ({exc}); "
+                  f"falling back to 3D costs.")
+            return None
+
     def do_tracking(self) -> None:
         if self.exp is None:
             raise ValueError("No experiment object provided")
@@ -244,7 +522,15 @@ def do_tracking(self) -> None:
 
         track_cfg = pm.parameters.get("track", {}) if pm else {}
         leaf_weight = float(track_cfg.get("leaf_weight", 1.0))
-        v_max = float(track_cfg.get("dvxmax", 15.5))
+        v_max = float(track_cfg.get("v_max", track_cfg.get("dvxmax", 15.5)))
+        use_velocity = bool(track_cfg.get("use_velocity", True))
+        cost_mode = str(track_cfg.get("cost_mode", "projected"))
+        max_gap = int(track_cfg.get("max_gap", 2))
+        allow_shared = bool(track_cfg.get("allow_shared", False))
+        max_shared = int(track_cfg.get("max_shared", 2))
+        share_tol_raw = track_cfg.get("share_tol", 1.0)
+        share_tol = None if share_tol_raw is None else float(share_tol_raw)
+        max_group_size = int(track_cfg.get("max_group_size", 128))
 
         store = getattr(self.exp, "_store", None)
         if store is None:
@@ -329,9 +615,17 @@ def do_tracking(self) -> None:
                     xy[valid, c] = t[cam_ids[valid, c], 1:3]
             frame_leaves.append(np.nan_to_num(xy.reshape(n, -1)))
 
-        cfg = TwoPhaseTrackerConfig(v_max=v_max, leaf_weight=leaf_weight)
+        cfg = TwoPhaseTrackerConfig(v_max=v_max, leaf_weight=leaf_weight,
+                                      use_velocity=use_velocity,
+                                      cost_mode=cost_mode, max_gap=max_gap,
+                                      allow_shared=allow_shared,
+                                      max_shared=max_shared,
+                                      share_tol=share_tol,
+                                      max_group_size=max_group_size)
         tracker = TwoPhaseTracker(cfg)
-        links = tracker.track_frames(frame_particles, frame_leaves)
+        project_fn = self._build_project_fn()
+        links = tracker.track_frames(frame_particles, frame_leaves,
+                                     project_fn=project_fn)
 
         # Per-step progress like trackcorr (track3d step: curr/next/links)
         from collections import Counter
@@ -375,5 +669,8 @@ def do_tracking(self) -> None:
 
         print(
             f"TwoPhaseTracker: {len(links)} links across {len(frames)} frames "
-            f"(leaf_weight={leaf_weight}, v_max={v_max})"
+            f"(leaf_weight={leaf_weight}, v_max={v_max}, "
+            f"use_velocity={use_velocity}, cost_mode={cost_mode}, "
+            f"max_gap={max_gap}, "
+            f"project_fn={'yes' if project_fn is not None else 'no'})"
         )
diff --git a/src/openptv2/tracking_shared.py b/src/openptv2/tracking_shared.py
new file mode 100644
index 00000000..b2d039a7
--- /dev/null
+++ b/src/openptv2/tracking_shared.py
@@ -0,0 +1,185 @@
+"""Shared-observation prototype for trackcorr, at linkage-file level.
+
+trackcorr's kernel (``track_kernels_corr._trackcorr_particle_fast``) is
+compiled per-particle greedy code -- the plan is to prototype the lesson in
+pure Python on its linkage files first and move it into the kernel (behind a
+flag) only if it wins:
+
+1. :func:`mark_shared_observations` -- after the forward pass, find frames
+   where one 3D point serves two tracks (detector undercount = occlusion)
+   and mark it shared instead of letting one track die.
+2. Downstream, shared points update position but never velocity (each
+   track's speed comes only from its own points); gap bridging
+   (``tracking_postprocess.relink_trajectory_gaps``) then stitches with
+   honest speeds.
+
+Linkage convention (mirrors ``tracking_postprocess``):
+    frames[k] = (prev_k, next_k, xyz_k) with prev/next int arrays
+    (PREV_NONE/NEXT_NONE = unlinked) and xyz_k (N,3) float positions.
+Shared marks: dict mapping (frame, idx) -> list of track-end ids sharing it.
+A "track end" here is identified by (frame, idx) of its last point.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+PREV_NONE = -1
+NEXT_NONE = -2
+
+
+def _velocity(frames, k, idx, shared, back_steps=1):
+    """Velocity at (k, idx) from own points only: walk back over prev links,
+    skipping shared marks. Returns None when fewer than two independent
+    points exist."""
+    pts = [(k, idx)]
+    ck, ci = k, idx
+    while True:
+        prev, _, _ = frames[ck]
+        pi = int(prev[ci])
+        if pi < 0:
+            break
+        ck -= 1
+        if ck not in frames:
+            break
+        ci = pi
+        pts.append((ck, ci))
+        if len(pts) >= back_steps + 2:
+            break
+    indep = [(fk, fi) for (fk, fi) in pts if (fk, fi) not in shared]
+    if len(indep) < 2:
+        return None
+    (f1, i1), (f0, i0) = indep[-1], indep[-2]
+    dt = max(f1 - f0, 1)
+    _, _, xyz1 = frames[f1]
+    _, _, xyz0 = frames[f0]
+    return (xyz1[i1] - xyz0[i0]) / dt
+
+
+def mark_shared_observations(frames, first, last, tol, max_share=2):
+    """Find occlusions: a frame-(k+1) detection predicted well by >= 2
+    distinct frame-k tracks that have a velocity (prev link), where at most
+    one of them actually claimed it.
+
+    Returns shared: dict (frame, idx) -> list of (frame_k, idx_k) sharers.
+    Does not modify linkages; assembly/relink consume the marks.
+    """
+    shared: dict[tuple[int, int], list[tuple[int, int]]] = {}
+    for k in range(first, last):
+        if k not in frames or (k + 1) not in frames:
+            continue
+        prev_k, next_k, xyz_k = frames[k]
+        _, _, xyz_n = frames[k + 1]
+        # predictors: frame-k particles with a prev link (have velocity)
+        for j in range(len(xyz_n)):
+            claimants = []
+            for i in range(len(xyz_k)):
+                pi = int(prev_k[i])
+                if pi < 0 or k - 1 not in frames:
+                    continue
+                _, _, xyz_p = frames[k - 1]
+                if pi >= len(xyz_p):
+                    continue
+                v = xyz_k[i] - xyz_p[pi]
+                if np.linalg.norm(xyz_n[j] - (xyz_k[i] + v)) < tol:
+                    claimants.append((i, float(
+                        np.linalg.norm(xyz_n[j] - (xyz_k[i] + v)))))
+            # undercount signature: >= 2 predictors, <= 1 actual claim.
+            # actual claims: frame-k particles with next_k[i] == j
+            actual = [i for i in range(len(xyz_k)) if int(next_k[i]) == j]
+            if len(claimants) >= 2 and len(actual) <= 1:
+                key = (k + 1, j)
+                shared.setdefault(key, [])
+                for i, _ in sorted(claimants, key=lambda t: t[1]):
+                    if (k, i) not in shared[key]:
+                        shared[key].append((k, i))
+                shared[key] = shared[key][: max_share + 1]
+    return shared
+
+
+def assemble_with_shared(frames, first, last, shared):
+    """Assemble tracks from linkage arrays, emitting a shared point into
+    every sharing track (position shared, histories stay distinct because
+    each track keeps its own prev chain).
+
+    Returns list of dicts: {frames: [...], pos: (L,3)}.
+    """
+    # forward chains from every unclaimed start; shared points entered
+    # once per sharing track.
+    tracks = []
+    visited = set()  # (frame, idx, owner-key) to allow shared re-entry
+    # 1. ordinary chains from particles with no prev link
+    for k in range(first, last + 1):
+        if k not in frames:
+            continue
+        prev_k, next_k, xyz_k = frames[k]
+        for i in range(len(xyz_k)):
+            if int(prev_k[i]) >= 0:
+                continue
+            chain = []
+            ck, ci = k, i
+            while True:
+                chain.append((ck, ci))
+                if ck not in frames:
+                    break
+                _, nxt, _ = frames[ck]
+                ni = int(nxt[ci]) if ci < len(nxt) else NEXT_NONE
+                if ni < 0:
+                    break
+                ck += 1
+                if ck not in frames:
+                    break
+                ci = ni
+            tracks.append(chain)
+            for node in chain:
+                visited.add((node, "main"))
+    # 2. sharing tracks: re-walk from each sharer end through the shared
+    # point, then STOP. The shared point is observed jointly, but what
+    # follows belongs to whoever links there next -- riding the winner's
+    # tail would graft the wrong identity onto the sharer (a switch).
+    # Continuation past separation is gap-relink's job, not assembly's.
+    for (fk, fj), sharers in shared.items():
+        for (sk, si) in sharers:
+            # walk the sharer's own history up to (sk, si)
+            hist = []
+            ck, ci = sk, si
+            while True:
+                hist.append((ck, ci))
+                if ck not in frames:
+                    break
+                pv, _, _ = frames[ck]
+                pi = int(pv[ci]) if ci < len(pv) else PREV_NONE
+                if pi < 0:
+                    break
+                ck -= 1
+                if ck not in frames:
+                    break
+                ci = pi
+            hist.reverse()
+            # then the shared point itself -- and stop. Whatever continues
+            # past separation belongs to whoever links there next; following
+            # the winner's tail would graft the wrong identity (a switch).
+            full = hist + [(fk, fj)] if (fk, fj) not in hist else list(hist)
+            key = ("shared", sk, si, fk, fj)
+            if key in visited:
+                continue
+            visited.add(key)
+            # skip if an identical main chain already covers it
+            if any(all(n in c for n in full) for c in tracks
+                   if len(c) >= len(full)):
+                continue
+            tracks.append(full)
+    out = []
+    for chain in tracks:
+        fr, ps = [], []
+        for (ck, ci) in chain:
+            if ck not in frames:
+                continue
+            _, _, xyz = frames[ck]
+            if ci >= len(xyz):
+                continue
+            fr.append(ck)
+            ps.append(xyz[ci])
+        if len(fr) >= 1:
+            out.append({"frames": fr, "pos": np.array(ps)})
+    return out
diff --git a/tests/unit/test_track_kernels_batch_coverage.py b/tests/unit/test_track_kernels_batch_coverage.py
index 7234ff45..1a7ecd76 100644
--- a/tests/unit/test_track_kernels_batch_coverage.py
+++ b/tests/unit/test_track_kernels_batch_coverage.py
@@ -25,10 +25,6 @@
 
 from openptv2.algorithms.track_kernels_batch import (  # noqa: E402
     init_mmlut_data_fast,
-    metric_to_pixel_batch_fast,
-    pixel_to_metric_batch_fast,
-    point_position_batch_fast,
-    ray_tracing_batch_fast,
     targ_rec_fast,
 )
 
@@ -89,46 +85,6 @@ def _make_cal(
 # ─────────────────────────────────────────────────────────────────────────────
 
 
-def test_ray_tracing_batch_empty():
-    """N=0 input → (0, 3) position and direction arrays."""
-    xy = np.empty((0, 2), dtype=np.float64)
-    pos, dirs = ray_tracing_batch_fast(xy, _CAL)
-    assert pos.shape == (0, 3)
-    assert dirs.shape == (0, 3)
-
-
-def test_ray_tracing_batch_single():
-    """N=1 → (1, 3) outputs; values are finite."""
-    xy = np.array([[0.0, 0.0]], dtype=np.float64)
-    pos, dirs = ray_tracing_batch_fast(xy, _CAL)
-    assert pos.shape == (1, 3)
-    assert dirs.shape == (1, 3)
-    assert np.all(np.isfinite(pos))
-    assert np.all(np.isfinite(dirs))
-
-
-def test_ray_tracing_batch_multiple():
-    """N=5 → (5, 3) outputs; all finite."""
-    xy = np.array(
-        [[-5.0, -5.0], [-2.0, 0.0], [0.0, 0.0], [2.0, 0.0], [5.0, 5.0]],
-        dtype=np.float64,
-    )
-    pos, dirs = ray_tracing_batch_fast(xy, _CAL)
-    assert pos.shape == (5, 3)
-    assert dirs.shape == (5, 3)
-    assert np.all(np.isfinite(pos))
-    assert np.all(np.isfinite(dirs))
-
-
-def test_ray_tracing_batch_off_axis():
-    """Non-zero x0,y0 camera position still produces finite rays."""
-    cal_off = _make_cal(x0=10.0, y0=5.0, z0=80.0)
-    xy = np.array([[1.0, -1.0], [0.5, 0.5]], dtype=np.float64)
-    pos, dirs = ray_tracing_batch_fast(xy, cal_off)
-    assert pos.shape == (2, 3)
-    assert np.all(np.isfinite(pos))
-
-
 # ─────────────────────────────────────────────────────────────────────────────
 # pixel_to_metric_batch_fast
 # ─────────────────────────────────────────────────────────────────────────────
@@ -137,89 +93,11 @@ def test_ray_tracing_batch_off_axis():
 _PIXX, _PIXY = 0.017, 0.017
 
 
-def test_pixel_to_metric_batch_empty():
-    """N=0 → (0, 2) result, no crash."""
-    xy = np.empty((0, 2), dtype=np.float64)
-    result = pixel_to_metric_batch_fast(xy, _IMX, _IMY, _PIXX, _PIXY, 0)
-    assert result.shape == (0, 2)
-
-
-def test_pixel_to_metric_batch_chfield0():
-    """chfield=0 → standard pixel-to-metric; shape (3, 2), finite."""
-    xy = np.array([[320.0, 240.0], [0.0, 0.0], [640.0, 480.0]], dtype=np.float64)
-    result = pixel_to_metric_batch_fast(xy, _IMX, _IMY, _PIXX, _PIXY, 0)
-    assert result.shape == (3, 2)
-    assert np.all(np.isfinite(result))
-
-
-def test_pixel_to_metric_batch_chfield1():
-    """chfield=1 → yp = 2*y + 1 branch executed."""
-    xy = np.array([[100.0, 100.0]], dtype=np.float64)
-    result = pixel_to_metric_batch_fast(xy, _IMX, _IMY, _PIXX, _PIXY, 1)
-    assert result.shape == (1, 2)
-    assert np.isfinite(result[0, 1])
-
-
-def test_pixel_to_metric_batch_chfield2():
-    """chfield=2 → yp = 2*y branch executed."""
-    xy = np.array([[100.0, 100.0]], dtype=np.float64)
-    result = pixel_to_metric_batch_fast(xy, _IMX, _IMY, _PIXX, _PIXY, 2)
-    assert result.shape == (1, 2)
-    assert np.isfinite(result[0, 1])
-
-
-def test_pixel_to_metric_batch_center():
-    """Image centre maps to metric origin (0, 0) for chfield=0."""
-    xy = np.array([[_IMX / 2.0, _IMY / 2.0]], dtype=np.float64)
-    result = pixel_to_metric_batch_fast(xy, _IMX, _IMY, _PIXX, _PIXY, 0)
-    assert abs(result[0, 0]) < 1e-10
-    assert abs(result[0, 1]) < 1e-10
-
-
 # ─────────────────────────────────────────────────────────────────────────────
 # metric_to_pixel_batch_fast
 # ─────────────────────────────────────────────────────────────────────────────
 
 
-def test_metric_to_pixel_batch_chfield0():
-    """chfield=0 → standard metric-to-pixel; shape (3, 2), finite."""
-    xy = np.array([[0.0, 0.0], [1.0, 1.0], [-1.0, -1.0]], dtype=np.float64)
-    result = metric_to_pixel_batch_fast(xy, _IMX, _IMY, _PIXX, _PIXY, 0)
-    assert result.shape == (3, 2)
-    assert np.all(np.isfinite(result))
-
-
-def test_metric_to_pixel_batch_chfield1():
-    """chfield=1 → y_pixel = (y_pixel - 1) * 0.5 branch executed."""
-    xy = np.array([[0.5, 0.5]], dtype=np.float64)
-    result = metric_to_pixel_batch_fast(xy, _IMX, _IMY, _PIXX, _PIXY, 1)
-    assert result.shape == (1, 2)
-    assert np.isfinite(result[0, 1])
-
-
-def test_metric_to_pixel_batch_chfield2():
-    """chfield=2 → y_pixel = y_pixel * 0.5 branch executed."""
-    xy = np.array([[0.5, 0.5]], dtype=np.float64)
-    result = metric_to_pixel_batch_fast(xy, _IMX, _IMY, _PIXX, _PIXY, 2)
-    assert result.shape == (1, 2)
-    assert np.isfinite(result[0, 1])
-
-
-def test_pixel_metric_roundtrip():
-    """pixel→metric→pixel recovers original coordinates (chfield=0)."""
-    pts_px = np.array([[320.0, 240.0], [100.0, 380.0]], dtype=np.float64)
-    metric = pixel_to_metric_batch_fast(pts_px, _IMX, _IMY, _PIXX, _PIXY, 0)
-    back = metric_to_pixel_batch_fast(metric, _IMX, _IMY, _PIXX, _PIXY, 0)
-    assert np.allclose(back, pts_px, atol=1e-8)
-
-
-def test_metric_pixel_empty():
-    """N=0 metric_to_pixel → (0, 2) result."""
-    xy = np.empty((0, 2), dtype=np.float64)
-    result = metric_to_pixel_batch_fast(xy, _IMX, _IMY, _PIXX, _PIXY, 0)
-    assert result.shape == (0, 2)
-
-
 # ─────────────────────────────────────────────────────────────────────────────
 # point_position_batch_fast
 # ─────────────────────────────────────────────────────────────────────────────
@@ -232,47 +110,6 @@ def _two_cams():
     return (cal1, cal2)
 
 
-def test_point_position_batch_empty():
-    """num_pts=0 → (0, 3) positions and (0,) distances."""
-    all_targets = np.empty((0, 2, 2), dtype=np.float64)
-    cal_arrays = _two_cams()
-    positions, distances = point_position_batch_fast(all_targets, 0, 2, cal_arrays)
-    assert positions.shape == (0, 3)
-    assert distances.shape == (0,)
-
-
-def test_point_position_batch_one_point():
-    """num_pts=1, num_cams=2 → (1, 3) and (1,); finite values."""
-    all_targets = np.zeros((1, 2, 2), dtype=np.float64)
-    cal_arrays = _two_cams()
-    positions, distances = point_position_batch_fast(all_targets, 1, 2, cal_arrays)
-    assert positions.shape == (1, 3)
-    assert distances.shape == (1,)
-    assert np.all(np.isfinite(positions))
-    assert np.isfinite(distances[0])
-
-
-def test_point_position_batch_multiple_points():
-    """num_pts=3 → (3, 3) positions and (3,) distances."""
-    all_targets = np.zeros((3, 2, 2), dtype=np.float64)
-    cal_arrays = _two_cams()
-    positions, distances = point_position_batch_fast(all_targets, 3, 2, cal_arrays)
-    assert positions.shape == (3, 3)
-    assert distances.shape == (3,)
-
-
-def test_point_position_batch_nonzero_targets():
-    """Finite target coords still yield finite positions."""
-    all_targets = np.array(
-        [[[1.0, 2.0], [-1.0, 2.0]], [[0.5, 0.5], [-0.5, 0.5]]],
-        dtype=np.float64,
-    )
-    cal_arrays = _two_cams()
-    positions, distances = point_position_batch_fast(all_targets, 2, 2, cal_arrays)
-    assert positions.shape == (2, 3)
-    assert np.all(np.isfinite(positions))
-
-
 # ─────────────────────────────────────────────────────────────────────────────
 # targ_rec_fast — helpers
 # ─────────────────────────────────────────────────────────────────────────────
diff --git a/tests/unit/test_track_kernels_coverage.py b/tests/unit/test_track_kernels_coverage.py
index e1aafce9..5adab7f4 100644
--- a/tests/unit/test_track_kernels_coverage.py
+++ b/tests/unit/test_track_kernels_coverage.py
@@ -31,11 +31,9 @@
     is_compiled,
     pack_cal_array,
     pack_mmlut,
-    point_position_fast,
     # re-exports from sub-modules
     point_to_pixel_fast,
     searchquader_fast,
-    sort_candidates_by_freq_fast,
     sorted_candidates_fast,
     targ_rec_fast,
     track3d_loop_fast,
@@ -301,9 +299,7 @@ def test_reexported_callables_exist():
         searchquader_fast,
         candsearch_in_pix_fast,
         candsearch_in_pix_rest_fast,
-        sort_candidates_by_freq_fast,
         sorted_candidates_fast,
-        point_position_fast,
         trackcorr_loop_fast,
         trackback_loop_fast,
         track3d_loop_fast,
diff --git a/tests/unit/test_track_kernels_geom_coverage.py b/tests/unit/test_track_kernels_geom_coverage.py
index 564862e4..3f9db554 100644
--- a/tests/unit/test_track_kernels_geom_coverage.py
+++ b/tests/unit/test_track_kernels_geom_coverage.py
@@ -19,14 +19,16 @@
 from openptv2.algorithms.track_kernels_geom import (
     CAL_ARRAY_SIZE,
     PT_UNUSED,
-    _angle_acc_out,
+    point_to_pixel_fast,
+    searchquader_fast,
+)
+from openptv2.algorithms.track_kernels_pixel import (
     _multimed_r_nlay_1layer,
     _point_to_pixel_out,
-    _ray_tracing_fast,
+)
+from openptv2.algorithms.track_kernels_position import (
+    _angle_acc_out,
     _ray_tracing_out,
-    angle_acc_fast,
-    point_to_pixel_fast,
-    searchquader_fast,
 )
 
 # ---------------------------------------------------------------------------
@@ -1006,160 +1008,6 @@ def test_searchquader_zero_quader():
     assert xr.shape == (1,)
 
 
-# ---------------------------------------------------------------------------
-# angle_acc_fast
-# ---------------------------------------------------------------------------
-
-
-def test_angle_acc_fast_same_vectors_zero():
-    """v0 == v1 → angle = 0.0, acc = 0.0."""
-    angle, acc = angle_acc_fast(
-        0.0,
-        0.0,
-        0.0,  # start
-        1.0,
-        0.0,
-        0.0,  # pred
-        1.0,
-        0.0,
-        0.0,  # cand  (same as pred)
-    )
-    assert angle == 0.0
-    assert acc == 0.0
-
-
-def test_angle_acc_fast_opposite_vectors_200():
-    """v0 == -v1 → angle = 200.0."""
-    angle, acc = angle_acc_fast(
-        0.0,
-        0.0,
-        0.0,
-        1.0,
-        0.0,
-        0.0,
-        -1.0,
-        0.0,
-        0.0,
-    )
-    assert angle == 200.0
-
-
-def test_angle_acc_fast_90_degrees():
-    """Perpendicular vectors → angle ≈ 100.0 (90° scaled to 200/π·rad)."""
-    angle, acc = angle_acc_fast(
-        0.0,
-        0.0,
-        0.0,
-        1.0,
-        0.0,
-        0.0,
-        0.0,
-        1.0,
-        0.0,
-    )
-    assert abs(angle - 100.0) < 0.1
-    assert np.isfinite(acc)
-
-
-def test_angle_acc_fast_norm0_zero():
-    """start == pred → v0 = (0,0,0) → norm0 = 0 → angle = 0.0."""
-    angle, acc = angle_acc_fast(
-        0.0,
-        0.0,
-        0.0,
-        0.0,
-        0.0,
-        0.0,
-        1.0,
-        0.0,
-        0.0,
-    )
-    assert angle == 0.0
-
-
-def test_angle_acc_fast_norm1_zero():
-    """start == cand → v1 = (0,0,0) → norm1 = 0 → angle = 0.0."""
-    angle, acc = angle_acc_fast(
-        0.0,
-        0.0,
-        0.0,
-        1.0,
-        0.0,
-        0.0,
-        0.0,
-        0.0,
-        0.0,
-    )
-    assert angle == 0.0
-
-
-def test_angle_acc_fast_nearly_parallel():
-    """Almost parallel vectors — dot may be > 1 in floating point → clamped."""
-    eps = 1e-14
-    angle, acc = angle_acc_fast(
-        0.0,
-        0.0,
-        0.0,
-        1.0,
-        0.0,
-        0.0,
-        1.0 + eps,
-        0.0,
-        0.0,
-    )
-    assert 0.0 <= angle <= 200.0
-
-
-def test_angle_acc_fast_acceleration_value():
-    """Acc is the distance between v1 and v0."""
-    angle, acc = angle_acc_fast(
-        0.0,
-        0.0,
-        0.0,
-        1.0,
-        0.0,
-        0.0,
-        2.0,
-        0.0,
-        0.0,
-    )
-    # v0=(1,0,0), v1=(2,0,0) → dx=1 → acc=1
-    assert abs(acc - 1.0) < 1e-10
-
-
-def test_angle_acc_fast_3d_vectors():
-    angle, acc = angle_acc_fast(
-        0.0,
-        0.0,
-        0.0,
-        1.0,
-        1.0,
-        1.0,
-        1.0,
-        -1.0,
-        0.0,
-    )
-    assert 0.0 <= angle <= 200.0
-    assert np.isfinite(acc)
-
-
-def test_angle_acc_fast_negative_dot_clamped():
-    """Antiparallel but not exact → dot < -1 gets clamped to -1."""
-    # Make two nearly-opposite unit vectors with floating-point excess
-    angle, acc = angle_acc_fast(
-        0.0,
-        0.0,
-        0.0,
-        1.0,
-        0.0,
-        0.0,
-        -1.0,
-        1e-15,
-        0.0,  # nearly opposite, not exactly
-    )
-    assert 0.0 <= angle <= 200.0
-
-
 # ---------------------------------------------------------------------------
 # _angle_acc_out
 # ---------------------------------------------------------------------------
@@ -1186,14 +1034,6 @@ def test_angle_acc_out_90_degrees():
     assert np.isfinite(out[1])
 
 
-def test_angle_acc_out_matches_fast():
-    out = np.zeros(2, dtype=np.float64)
-    _angle_acc_out(0.0, 0.0, 0.0, 2.0, 1.0, 0.0, 1.0, 2.0, 0.5, out)
-    angle, acc = angle_acc_fast(0.0, 0.0, 0.0, 2.0, 1.0, 0.0, 1.0, 2.0, 0.5)
-    assert abs(out[0] - angle) < 1e-10
-    assert abs(out[1] - acc) < 1e-10
-
-
 def test_angle_acc_out_norm0_zero():
     out = np.zeros(2, dtype=np.float64)
     _angle_acc_out(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, out)
@@ -1213,75 +1053,6 @@ def test_angle_acc_out_3d():
     assert np.isfinite(out[1])
 
 
-# ---------------------------------------------------------------------------
-# _ray_tracing_fast
-# ---------------------------------------------------------------------------
-
-
-def test_ray_tracing_fast_on_axis_tuple_len():
-    cal = _make_cal_array()
-    result = _ray_tracing_fast(0.0, 0.0, cal)
-    assert len(result) == 6
-    assert all(np.isfinite(v) for v in result)
-
-
-def test_ray_tracing_fast_off_axis():
-    cal = _make_cal_array()
-    Xx, Xy, Xz, ox, oy, oz = _ray_tracing_fast(1.0, 0.5, cal)
-    assert all(np.isfinite(v) for v in [Xx, Xy, Xz, ox, oy, oz])
-
-
-def test_ray_tracing_fast_negative_xy():
-    cal = _make_cal_array()
-    result = _ray_tracing_fast(-2.0, -1.0, cal)
-    assert all(np.isfinite(v) for v in result)
-
-
-def test_ray_tracing_fast_tilted_glass():
-    """Tilted glass → non-trivial Snell refraction."""
-    cal = _make_cal_array(gx=0.5, gy=0.0, gz=1.0)
-    Xx, Xy, Xz, ox, oy, oz = _ray_tracing_fast(1.0, 0.0, cal)
-    assert np.isfinite(Xx)
-    assert np.isfinite(ox)
-
-
-def test_ray_tracing_fast_on_axis_bpn_zero():
-    """x=0, y=0, dm=identity, glass=[0,0,1] → start_dir·glass parallel → bpn=0."""
-    cal = _make_cal_array(gx=0.0, gy=0.0, gz=1.0)
-    result = _ray_tracing_fast(0.0, 0.0, cal)
-    assert len(result) == 6
-
-
-def test_ray_tracing_fast_varied_indices():
-    cal = _make_cal_array(mm_n1=1.33, mm_n2_0=1.5, mm_n3=1.33, mm_d0=5.0)
-    result = _ray_tracing_fast(2.0, 1.0, cal)
-    assert all(np.isfinite(v) for v in result)
-
-
-def test_ray_tracing_fast_large_xy():
-    cal = _make_cal_array()
-    result = _ray_tracing_fast(8.0, 6.0, cal)
-    assert len(result) == 6
-
-
-def test_ray_tracing_fast_symmetry_x():
-    """_ray_tracing_fast(-x, y) mirrors _ray_tracing_fast(x, y) in X."""
-    cal = _make_cal_array()
-    Xx_p, Xy_p, Xz_p, _, _, _ = _ray_tracing_fast(2.0, 0.0, cal)
-    Xx_n, Xy_n, Xz_n, _, _, _ = _ray_tracing_fast(-2.0, 0.0, cal)
-    assert abs(Xx_p + Xx_n) < 1e-10
-    assert abs(Xy_p - Xy_n) < 1e-10
-
-
-def test_ray_tracing_fast_zero_glass_gn_zero_branch():
-    """gx=gy=gz=0 → gn=0 branch executed (raises ZeroDivision later — acceptable)."""
-    cal = _make_cal_array(gx=0.0, gy=0.0, gz=0.0)
-    try:
-        _ray_tracing_fast(1.0, 0.0, cal)
-    except (ZeroDivisionError, ValueError):
-        pass  # branch covered; exception is expected
-
-
 # ---------------------------------------------------------------------------
 # _ray_tracing_out
 # ---------------------------------------------------------------------------
@@ -1302,19 +1073,6 @@ def test_ray_tracing_out_off_axis():
     assert all(np.isfinite(out[i]) for i in range(6))
 
 
-def test_ray_tracing_out_matches_fast():
-    cal = _make_cal_array()
-    out = np.zeros(6, dtype=np.float64)
-    _ray_tracing_out(1.0, 0.5, cal, out)
-    Xx, Xy, Xz, ox, oy, oz = _ray_tracing_fast(1.0, 0.5, cal)
-    assert abs(out[0] - Xx) < 1e-10
-    assert abs(out[1] - Xy) < 1e-10
-    assert abs(out[2] - Xz) < 1e-10
-    assert abs(out[3] - ox) < 1e-10
-    assert abs(out[4] - oy) < 1e-10
-    assert abs(out[5] - oz) < 1e-10
-
-
 def test_ray_tracing_out_negative_xy():
     cal = _make_cal_array()
     out = np.zeros(6, dtype=np.float64)
diff --git a/tests/unit/test_track_kernels_search_coverage.py b/tests/unit/test_track_kernels_search_coverage.py
index 66a6f951..e27dbbfb 100644
--- a/tests/unit/test_track_kernels_search_coverage.py
+++ b/tests/unit/test_track_kernels_search_coverage.py
@@ -74,22 +74,17 @@
 # ---------------------------------------------------------------------------
 # Imports from the module under test
 # ---------------------------------------------------------------------------
+from openptv2.algorithms.track_kernels_pixel import (
+    _multimed_r_nlay_1layer,
+    _point_to_pixel_out,
+)
 from openptv2.algorithms.track_kernels_search import (
     _sorted_candidates_fast_out,
-    _sorted_candidates_fast_out_nogil,
     candsearch_in_pix_fast,
-    candsearch_in_pix_fast_nogil,
     candsearch_in_pix_rest_fast,
-    sort_candidates_by_freq_fast,
     sorted_candidates_fast,
 )
 
-if not _is_compiled():
-    from openptv2.algorithms.track_kernels_search import (
-        _multimed_r_nlay_1layer,
-        _point_to_pixel_out,
-    )
-
 EPS = 1e-8
 
 # ---------------------------------------------------------------------------
@@ -887,504 +882,6 @@ def test_large_num_targets(self):
         assert cnt == 1
 
 
-# ---------------------------------------------------------------------------
-# sort_candidates_by_freq_fast
-# ---------------------------------------------------------------------------
-
-
-class TestSortCandidatesByFreqFast:
-    NUM_CAMS = 4
-    MAX_CANDS = 4
-
-    def _make_arrays(self, ftnr_vals):
-        n = len(ftnr_vals)
-        ftnr = np.asarray(ftnr_vals, dtype=np.int32)
-        freq = np.zeros(n, dtype=np.int32)
-        whichcam = np.zeros((n, self.NUM_CAMS), dtype=np.int32)
-        return ftnr, freq, whichcam
-
-    def test_all_unused(self):
-        ftnr, freq, whichcam = self._make_arrays([-1] * 16)
-        nv = sort_candidates_by_freq_fast(
-            ftnr, freq, whichcam, 16, self.NUM_CAMS, self.MAX_CANDS
-        )
-        assert nv == 0
-
-    def test_single_candidate_one_cam(self):
-        # One candidate in camera 0 slot, rest unused
-        vals = [-1] * 16
-        vals[0] = 5  # cam0, slot0
-        ftnr, freq, whichcam = self._make_arrays(vals)
-        nv = sort_candidates_by_freq_fast(
-            ftnr, freq, whichcam, 16, self.NUM_CAMS, self.MAX_CANDS
-        )
-        # Position 0 entry keeps freq=1 (dedup zeroes only j > i slots)
-        assert nv == 1
-
-    def test_same_candidate_two_cams(self):
-        # Target 10 appears in cam0 slot0 and cam1 slot0
-        vals = [-1] * 16
-        vals[0] = 10  # cam0 slot0
-        vals[4] = 10  # cam1 slot0
-        ftnr, freq, whichcam = self._make_arrays(vals)
-        nv = sort_candidates_by_freq_fast(
-            ftnr, freq, whichcam, 16, self.NUM_CAMS, self.MAX_CANDS
-        )
-        assert nv >= 1
-
-    def test_same_candidate_all_cams(self):
-        # Target 7 seen in all 4 cameras
-        vals = [-1] * 16
-        vals[0] = 7  # cam0
-        vals[4] = 7  # cam1
-        vals[8] = 7  # cam2
-        vals[12] = 7  # cam3
-        ftnr, freq, whichcam = self._make_arrays(vals)
-        nv = sort_candidates_by_freq_fast(
-            ftnr, freq, whichcam, 16, self.NUM_CAMS, self.MAX_CANDS
-        )
-        assert nv >= 1
-        # First entry should be target 7 with freq=4
-        assert ftnr[0] == 7
-
-    def test_duplicate_elimination(self):
-        # Same target twice in same camera → should be deduped
-        vals = [-1] * 16
-        vals[0] = 3
-        vals[1] = 3  # duplicate
-        vals[4] = 3  # also in cam1
-        ftnr, freq, whichcam = self._make_arrays(vals)
-        nv = sort_candidates_by_freq_fast(
-            ftnr, freq, whichcam, 16, self.NUM_CAMS, self.MAX_CANDS
-        )
-        assert nv >= 1
-        # Only one unique entry should remain
-        active = [(ftnr[i], freq[i]) for i in range(16) if freq[i] > 0]
-        tnrs = [t for t, f in active]
-        assert tnrs.count(3) <= 1
-
-    def test_sorting_higher_freq_first(self):
-        # Target A in 3 cams, target B in 2 cams — A should sort first
-        vals = [-1] * 16
-        vals[0] = 20  # cam0: target A
-        vals[4] = 20  # cam1: target A
-        vals[8] = 20  # cam2: target A
-        vals[1] = 30  # cam0 slot1: target B
-        vals[5] = 30  # cam1 slot1: target B
-        ftnr, freq, whichcam = self._make_arrays(vals)
-        nv = sort_candidates_by_freq_fast(
-            ftnr, freq, whichcam, 16, self.NUM_CAMS, self.MAX_CANDS
-        )
-        assert nv >= 1
-        if nv >= 2:
-            assert freq[0] >= freq[1]
-
-    def test_whichcam_swap(self):
-        # Test that whichcam is correctly swapped during sort
-        vals = [-1] * 16
-        vals[0] = 5
-        vals[4] = 5
-        vals[1] = 6
-        vals[5] = 6
-        vals[9] = 6
-        vals[13] = 6  # target 6 in all 4 cams → higher freq
-        ftnr, freq, whichcam = self._make_arrays(vals)
-        nv = sort_candidates_by_freq_fast(
-            ftnr, freq, whichcam, 16, self.NUM_CAMS, self.MAX_CANDS
-        )
-        # Target 6 (freq=4) should sort before target 5 (freq=2)
-        if nv >= 1:
-            assert ftnr[0] == 6
-
-    def test_returns_nonzero_for_freq_one_at_position_zero(self):
-        # Dedup loop only eliminates j > i, so index 0 always kept if freq > 0
-        vals = [-1] * 16
-        vals[0] = 99  # only in 1 camera → freq=1
-        ftnr, freq, whichcam = self._make_arrays(vals)
-        nv = sort_candidates_by_freq_fast(
-            ftnr, freq, whichcam, 16, self.NUM_CAMS, self.MAX_CANDS
-        )
-        # freq[0]=1 ≠ 0 → counted as valid
-        assert nv == 1
-
-    def test_freq_one_entry_after_another_gets_zeroed(self):
-        # Second unique candidate at freq=1 is eliminated by dedup loop at j>i
-        vals = [-1] * 16
-        vals[0] = 10  # cam0: target 10 in 2 cams → freq=2
-        vals[4] = 10  # cam1
-        vals[1] = 99  # cam0 slot1: target 99 in 1 cam → freq=1
-        ftnr, freq, whichcam = self._make_arrays(vals)
-        nv = sort_candidates_by_freq_fast(
-            ftnr, freq, whichcam, 16, self.NUM_CAMS, self.MAX_CANDS
-        )
-        # Target 10 (freq=2) survives; target 99 (freq=1) is zeroed
-        assert nv >= 1
-
-
-# ---------------------------------------------------------------------------
-# candsearch_in_pix_fast_nogil (cfunc — callable in pure Python)
-# ---------------------------------------------------------------------------
-
-
-@_needs_pure_python_loose_types
-class TestCandsearchInPixFastNogil:
-    """Tests for the nogil cfunc variant. Pass a list for out_indices."""
-
-    IMX, IMY = 1024.0, 1024.0
-
-    def test_center_out_of_bounds(self):
-        tx = np.array([100.0], dtype=np.float64)
-        ty = np.array([100.0], dtype=np.float64)
-        tnr = np.array([0], dtype=np.int32)
-        out = [0, 0, 0, 0]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            1,
-            2000.0,
-            512.0,
-            10.0,
-            10.0,
-            10.0,
-            10.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        assert all(v == -999 for v in out)
-
-    def test_empty_targets(self):
-        tx = np.zeros(0, dtype=np.float64)
-        ty = np.zeros(0, dtype=np.float64)
-        tnr = np.zeros(0, dtype=np.int32)
-        out = [0, 0, 0, 0]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            0,
-            512.0,
-            512.0,
-            10.0,
-            10.0,
-            10.0,
-            10.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        assert out[0] == -999
-
-    def test_single_target_found(self):
-        tx, ty, tnr = _make_targets([500.0], [500.0], [5])
-        out = [-999, -999, -999, -999]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            1,
-            500.0,
-            500.0,
-            20.0,
-            20.0,
-            20.0,
-            20.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        assert out[0] == 0  # index 0 in sorted arrays
-        assert out[1] == -999
-
-    def test_unused_tnr_skipped(self):
-        tx, ty, tnr = _make_targets([500.0], [500.0], [TR])
-        out = [-999, -999, -999, -999]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            1,
-            500.0,
-            500.0,
-            20.0,
-            20.0,
-            20.0,
-            20.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        assert out[0] == -999
-
-    def test_four_targets_fills_all_slots(self):
-        n = 5
-        xs = [500.0, 501.0, 502.0, 503.0, 504.0]
-        ys = [500.0, 500.1, 500.2, 500.3, 500.4]
-        tnrs = [10, 11, 12, 13, 14]
-        tx, ty, tnr = _make_targets(xs, ys, tnrs)
-        out = [-999, -999, -999, -999]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            n,
-            500.0,
-            500.0,
-            10.0,
-            10.0,
-            10.0,
-            10.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        found = [v for v in out if v != -999]
-        assert len(found) == 4
-
-    def test_d2_replacement(self):
-        # Target 1 is closest, target 2 is second closest
-        tx, ty, tnr = _make_targets([500.0, 500.5], [500.0, 500.5], [1, 2])
-        out = [-999, -999, -999, -999]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            2,
-            500.0,
-            500.0,
-            10.0,
-            10.0,
-            10.0,
-            10.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        found = [v for v in out if v != -999]
-        assert len(found) == 2
-
-    def test_d3_replacement(self):
-        tx, ty, tnr = _make_targets(
-            [500.0, 500.5, 501.0], [500.0, 500.5, 501.0], [1, 2, 3]
-        )
-        out = [-999, -999, -999, -999]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            3,
-            500.0,
-            500.0,
-            10.0,
-            10.0,
-            10.0,
-            10.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        found = [v for v in out if v != -999]
-        assert len(found) == 3
-
-    def test_ymax_early_break(self):
-        tx, ty, tnr = _make_targets([500.0, 500.0], [502.0, 600.0], [1, 2])
-        out = [-999, -999, -999, -999]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            2,
-            500.0,
-            500.0,
-            5.0,
-            5.0,
-            5.0,
-            5.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        found = [v for v in out if v != -999]
-        assert len(found) <= 1
-
-    def test_large_num_targets(self):
-        n = 60
-        xs = np.linspace(490.0, 510.0, n)
-        ys = np.linspace(490.0, 510.0, n)
-        tnrs = np.arange(n, dtype=np.int32)
-        tx, ty, tnr = _make_targets(xs, ys, tnrs)
-        out = [-999, -999, -999, -999]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            n,
-            500.0,
-            500.0,
-            12.0,
-            12.0,
-            12.0,
-            12.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        found = [v for v in out if v != -999]
-        assert len(found) == 4
-
-    def test_ymin_clamp(self):
-        """cent_y < du → ymin = 0.0 (covers line 1057)."""
-        tx, ty, tnr = _make_targets([512.0], [5.0], [1])
-        out = [-999, -999, -999, -999]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            1,
-            512.0,
-            5.0,
-            10.0,
-            10.0,
-            20.0,
-            10.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        assert out[0] != -999
-
-    def test_ymax_clamp(self):
-        """cent_y + dd > imy → ymax = imy (covers line 1059)."""
-        tx, ty, tnr = _make_targets([512.0], [1020.0], [1])
-        out = [-999, -999, -999, -999]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            1,
-            512.0,
-            1020.0,
-            10.0,
-            10.0,
-            10.0,
-            20.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        assert out[0] != -999
-
-    def test_xmin_clamp(self):
-        """cent_x < dl → xmin = 0.0 (covers line 1053)."""
-        tx, ty, tnr = _make_targets([5.0], [512.0], [1])
-        out = [-999, -999, -999, -999]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            1,
-            5.0,
-            512.0,
-            100.0,
-            10.0,
-            10.0,
-            10.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        assert out[0] != -999
-
-    def test_binary_search_j0_increment(self):
-        """Binary search j0 += dj (covers line 1081) and j0 no-clamp (1087->1090)."""
-        n = 40
-        xs = np.full(n, 512.0)
-        ys = np.linspace(100.0, 200.0, n)
-        tnrs = np.arange(n, dtype=np.int32)
-        tx, ty, tnr = _make_targets(xs, ys, tnrs)
-        out = [-999, -999, -999, -999]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            n,
-            512.0,
-            190.0,
-            5.0,
-            5.0,
-            5.0,
-            5.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        found = [v for v in out if v != -999]
-        assert len(found) >= 1
-
-    def test_j0_no_clamp(self):
-        """num_targets >= 24 → j0 - 12 >= 0 (covers branch 1087->1090)."""
-        n = 50
-        xs = np.full(n, 512.0)
-        ys = np.linspace(500.0, 510.0, n)
-        tnrs = np.arange(n, dtype=np.int32)
-        tx, ty, tnr = _make_targets(xs, ys, tnrs)
-        out = [-999, -999, -999, -999]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            n,
-            512.0,
-            505.0,
-            5.0,
-            5.0,
-            5.0,
-            5.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        found = [v for v in out if v != -999]
-        assert len(found) >= 1
-
-    def test_target_in_y_out_of_x(self):
-        """Target in y range but not x → False branch at tx check (covers 1096->1090)."""
-        tx, ty, tnr = _make_targets([900.0], [500.0], [1])
-        out = [-999, -999, -999, -999]
-        candsearch_in_pix_fast_nogil(
-            tx,
-            ty,
-            tnr,
-            1,
-            500.0,
-            500.0,
-            10.0,
-            10.0,
-            10.0,
-            10.0,
-            self.IMX,
-            self.IMY,
-            TR,
-            out,
-        )
-        assert all(v == -999 for v in out)
-
-
 # ---------------------------------------------------------------------------
 # sorted_candidates_fast — calls _sorted_candidates_fast_out which crashes
 # ---------------------------------------------------------------------------
@@ -1529,82 +1026,6 @@ def test_runs_with_two_cams(self):
         _sorted_candidates_fast_out(*args)
 
 
-@_needs_pure_python_loose_types
-class TestSortedCandidatesFastOutNogil:
-    """_sorted_candidates_fast_out_nogil runs in pure Python (C-array bug fixed 2026-07-10)."""
-
-    def _make_nogil_args(self, num_cams=1, max_cands=4, n_targ=5):
-        n = num_cams * max_cands
-        center = np.array([0.0, 0.0, 100.0], dtype=np.float64)
-        center_proj_x = np.zeros(num_cams, dtype=np.float64)
-        center_proj_y = np.zeros(num_cams, dtype=np.float64)
-
-        cal_arr = np.zeros((num_cams, 31), dtype=np.float64, order="C")
-        for i in range(num_cams):
-            cal_arr[i] = _make_cal()
-
-        # nogil variant takes up to 8 separate md arrays
-        empty_md = np.zeros(4, dtype=np.float64)
-        md_list = [empty_md] * 8
-
-        mo_arr = np.zeros((num_cams, 3), dtype=np.float64, order="C")
-        mnr_arr = np.zeros(num_cams, dtype=np.int32)
-        mnz_arr = np.zeros(num_cams, dtype=np.int32)
-        mrw_arr = np.ones(num_cams, dtype=np.float64)
-
-        targ_x = np.zeros((num_cams, n_targ), dtype=np.float64, order="C")
-        targ_y = np.zeros((num_cams, n_targ), dtype=np.float64, order="C")
-        targ_tnr = np.full((num_cams, n_targ), -1, dtype=np.int32, order="C")
-        num_targets = np.zeros(num_cams, dtype=np.int32)
-
-        ftnr_out = np.full(n, -1, dtype=np.int32)
-        freq_out = np.zeros(n, dtype=np.int32)
-        whichcam_out = np.zeros((n, num_cams), dtype=np.int32)
-
-        return (
-            center,
-            center_proj_x,
-            center_proj_y,
-            num_cams,
-            max_cands,
-            cal_arr,
-            *md_list,  # md0..md7
-            mo_arr,
-            mnr_arr,
-            mnz_arr,
-            mrw_arr,
-            targ_x,
-            targ_y,
-            targ_tnr,
-            num_targets,
-            -1.0,
-            1.0,
-            -1.0,
-            1.0,
-            -1.0,
-            1.0,
-            512.0,
-            512.0,
-            1.0,
-            1.0,
-            0,
-            1024.0,
-            1024.0,
-            -1,
-            ftnr_out,
-            freq_out,
-            whichcam_out,
-        )
-
-    def test_runs_without_error(self):
-        args = self._make_nogil_args()
-        _sorted_candidates_fast_out_nogil(*args)
-
-    def test_runs_with_two_cams(self):
-        args = self._make_nogil_args(num_cams=2)
-        _sorted_candidates_fast_out_nogil(*args)
-
-
 # ---------------------------------------------------------------------------
 # Compiled-mode sanity check (this module should be skipped when compiled)
 # ---------------------------------------------------------------------------
diff --git a/tests/unit/test_track_kernels_tracking_coverage.py b/tests/unit/test_track_kernels_tracking_coverage.py
index 9fa217ed..c80221c4 100644
--- a/tests/unit/test_track_kernels_tracking_coverage.py
+++ b/tests/unit/test_track_kernels_tracking_coverage.py
@@ -1,18 +1,8 @@
-"""Pure-Python coverage tests for track_kernels_tracking.py.
+"""Pure-Python coverage tests for the tracking kernels (corr/pixel/position/track3d).
 
 Skip when the compiled .so is active (coverage measures the .py source).
 
-Verification command (from repo root):
-    cp src/openptv2/algorithms/track_kernels_tracking.py \
-        /tmp/ppsrc/openptv2/algorithms/track_kernels_tracking.py
-    COVERAGE_FILE=/tmp/.cov_track_kernels_tracking \
-    uv run pytest tests/unit/test_track_kernels_tracking_coverage.py \
-      -o pythonpath=/tmp/ppsrc \
-      -p no:cacheprovider \
-      --cov=/tmp/ppsrc/openptv2 \
-      --cov-config=/tmp/covrc \
-      --cov-report=term-missing \
-      -q 2>&1 | grep -E '(algorithms/track_kernels_tracking\\.|TOTAL|passed|failed|error)'
+Run the interpreted-source variant per openptv2/CLAUDE.md (pure-Python fallback tests).
 """
 
 import numpy as np
@@ -34,34 +24,44 @@
 )
 
 import openptv2.algorithms.track_kernels_corr as _corr_mod
-import openptv2.algorithms.track_kernels_tracking as _mod
-from openptv2.algorithms.track_kernels_tracking import (
-    ADD_PART_K,
-    COORD_UNUSED_K,
-    CORRES_NONE_K,
-    MAX_CANDS_K,
-    NEXT_NONE_K,
-    POSI_K,
-    PREV_NONE_K,
-    PT_UNUSED,
-    TR_UNUSED_K,
-    _angle_acc_out,
+from openptv2.algorithms.track_kernels_corr import (
+    trackback_loop_fast,
+    trackcorr_loop_fast,
+)
+from openptv2.algorithms.track_kernels_pixel import (
     _candsearch_in_pix_rest_nogil,
     _dist_to_flat_out,
-    _find_closest_in_3d,
     _multimed_r_nlay_1layer,
     _pixel_to_metric_out,
-    _point_position_out,
     _point_to_pixel_out,
-    _ray_tracing_out,
     _sorted_candidates_fast_out_nogil,
-    assess_new_position_fast_nogil,
     candsearch_in_pix_fast_nogil,
+)
+from openptv2.algorithms.track_kernels_position import (
+    _angle_acc_out,
+    _point_position_out,
+    _ray_tracing_out,
+    assess_new_position_fast_nogil,
+)
+from openptv2.algorithms.track_kernels_track3d import (
+    _find_closest_in_3d,
     track3d_loop_fast,
-    trackback_loop_fast,
-    trackcorr_loop_fast,
 )
 
+# These mirror the cython.declare() C-level constants in track_kernels_corr,
+# which are not importable from Python when compiled.
+PT_UNUSED = -999
+POSI_K = 80
+MAX_CANDS_K = 32
+TR_UNUSED_K = -1
+CORRES_NONE_K = -1
+PREV_NONE_K = -1
+NEXT_NONE_K = -2
+COORD_UNUSED_K = -1e10
+ADD_PART_K = 3.0
+
+import openptv2.algorithms.track_kernels_pixel as _mod  # noqa: E402
+
 # ---------------------------------------------------------------------------
 # Helpers
 # ---------------------------------------------------------------------------
diff --git a/tests/unit/test_track_kernels_transform_coverage.py b/tests/unit/test_track_kernels_transform_coverage.py
index 196d8f7f..91e6eea2 100644
--- a/tests/unit/test_track_kernels_transform_coverage.py
+++ b/tests/unit/test_track_kernels_transform_coverage.py
@@ -52,27 +52,20 @@
 # ---------------------------------------------------------------------------
 # Imports under test
 # ---------------------------------------------------------------------------
-from openptv2.algorithms.track_kernels_transform import (
+from openptv2.algorithms.track_kernels_pixel import (
     _candsearch_in_pix_rest_nogil,
-    _flat_image_coord_fast,
-    _img_coord_fast,
-    _metric_to_pixel_out,
+    _dist_to_flat_out,
     _multimed_r_nlay_1layer,
     _pixel_to_metric_out,
+)
+from openptv2.algorithms.track_kernels_position import (
     _point_position_out,
     _ray_tracing_out,
-    assess_new_position_fast,
     assess_new_position_fast_nogil,
-    dist_to_flat_fast,
-    flat_image_coord_batch_fast,
-    img_coord_batch_fast,
-    metric_to_pixel_fast,
-    pixel_to_metric_fast,
-    point_position_fast,
 )
-
-if not _is_compiled():
-    from openptv2.algorithms.track_kernels_transform import _dist_to_flat_out
+from openptv2.algorithms.track_kernels_transform import (
+    assess_new_position_fast,
+)
 
 # ---------------------------------------------------------------------------
 # Helpers: build a minimal 31-element calibration flat array
@@ -443,64 +436,6 @@ def test_num_used_positive_branch(self):
         assert out is not None
 
 
-class TestPointPositionFast:
-    """Cover point_position_fast (line 495) — wraps _point_position_out."""
-
-    def test_returns_pos_and_dist(self):
-        """Smoke test: returns (pos_array, dist_float)."""
-        cal_arr = _make_cal_arr_batch(2)
-        targets = np.full((2, 2), COORD_UNUSED, dtype=np.float64)
-        pos, dist = point_position_fast(targets, 2, cal_arr)
-        assert pos.shape == (3,)
-        assert dist == 0.0
-
-    def test_two_cams_valid_targets(self):
-        """Two cameras with valid targets → non-trivial result."""
-        cal_arr = _make_cal_arr_batch(2)
-        targets = np.zeros((2, 2), dtype=np.float64)
-        pos, dist = point_position_fast(targets, 2, cal_arr)
-        assert pos.shape == (3,)
-        assert isinstance(dist, float)
-
-
-# ---------------------------------------------------------------------------
-# 4. pixel_to_metric_fast
-# ---------------------------------------------------------------------------
-
-
-class TestPixelToMetricFast:
-    def test_chfield_zero(self):
-        x_m, y_m = pixel_to_metric_fast(512.0, 384.0, 1024, 768, 0.01, 0.01, 0)
-        assert math.isclose(x_m, 0.0, abs_tol=1e-10)
-        assert math.isclose(y_m, 0.0, abs_tol=1e-10)
-
-    def test_chfield_one(self):
-        """chfield==1: yp = 2*y_pixel + 1."""
-        x_m, y_m = pixel_to_metric_fast(512.0, 100.0, 1024, 768, 0.01, 0.01, 1)
-        # yp = 2*100 + 1 = 201
-        expected_y = (768 * 0.5 - 201) * 0.01
-        assert math.isclose(y_m, expected_y, rel_tol=1e-9)
-
-    def test_chfield_two(self):
-        """chfield==2: yp = 2*y_pixel."""
-        x_m, y_m = pixel_to_metric_fast(512.0, 100.0, 1024, 768, 0.01, 0.01, 2)
-        yp = 2.0 * 100.0
-        expected_y = (768 * 0.5 - yp) * 0.01
-        assert math.isclose(y_m, expected_y, rel_tol=1e-9)
-
-    def test_origin_pixel(self):
-        """Pixel at image centre → metric (0,0)."""
-        x_m, y_m = pixel_to_metric_fast(512.0, 384.0, 1024, 768, 0.01, 0.01, 0)
-        assert abs(x_m) < 1e-10
-        assert abs(y_m) < 1e-10
-
-    def test_corner_pixel(self):
-        """Pixel at (0,0) → negative metric coords."""
-        x_m, y_m = pixel_to_metric_fast(0.0, 0.0, 1024, 768, 0.01, 0.01, 0)
-        assert x_m < 0.0
-        assert y_m > 0.0
-
-
 # ---------------------------------------------------------------------------
 # 5. _pixel_to_metric_out
 # ---------------------------------------------------------------------------
@@ -533,64 +468,6 @@ def test_returns_zero(self):
         assert ret == 0
 
 
-# ---------------------------------------------------------------------------
-# 6. dist_to_flat_fast
-# ---------------------------------------------------------------------------
-
-
-class TestDistToFlatFast:
-    def test_r_near_zero_returns_minus_xh_yh(self):
-        """Very small dist_x/dist_y → returns (-xh, -yh)."""
-        xh, yh = 0.5, -0.3
-        x, y = dist_to_flat_fast(
-            1e-15, 0.0, xh, yh, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1e-6
-        )
-        assert math.isclose(x, -xh, rel_tol=1e-9)
-        assert math.isclose(y, -yh, rel_tol=1e-9)
-
-    def test_zero_distortion_identity(self):
-        """k=p=0, scx=1, she=0 → output ≈ input (modulo xh/yh)."""
-        x, y = dist_to_flat_fast(
-            1.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1e-8
-        )
-        assert math.isclose(x, 1.0, rel_tol=1e-5)
-        assert math.isclose(y, 2.0, rel_tol=1e-5)
-
-    def test_with_principal_point_offset(self):
-        """Non-zero xh/yh shifts the result."""
-        xh, yh = 0.1, 0.2
-        x, y = dist_to_flat_fast(
-            1.0, 2.0, xh, yh, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1e-8
-        )
-        # xq starts at 1.0 and converges; result should be near 1.0 - xh
-        assert abs(x - (1.0 - xh)) < 0.01
-
-    def test_with_k1_distortion(self):
-        """Non-zero k1 changes the result noticeably."""
-        x0, y0 = dist_to_flat_fast(
-            2.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1e-8
-        )
-        x1, y1 = dist_to_flat_fast(
-            2.0, 1.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1e-8
-        )
-        assert x0 != x1
-
-    def test_she_nonzero(self):
-        """Non-zero shear angle exercises sin/cos branches."""
-        x, y = dist_to_flat_fast(
-            1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.01, 1e-8
-        )
-        assert isinstance(x, float)
-        assert isinstance(y, float)
-
-    def test_convergence_tol(self):
-        """Tight tolerance converges; result shifts slightly due to k1 correction."""
-        x, y = dist_to_flat_fast(
-            0.5, 0.5, 0.0, 0.0, 0.001, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1e-12
-        )
-        assert math.isclose(x, 0.5, abs_tol=0.01)
-
-
 # ---------------------------------------------------------------------------
 # 7. _dist_to_flat_out
 # ---------------------------------------------------------------------------
@@ -626,269 +503,6 @@ def test_with_distortion(self):
         assert math.isclose(out[0], 1.5, rel_tol=0.05)
 
 
-# ---------------------------------------------------------------------------
-# 8. metric_to_pixel_fast
-# ---------------------------------------------------------------------------
-
-
-class TestMetricToPixelFast:
-    def test_chfield_zero(self):
-        """metric_to_pixel_fast is inverse of pixel_to_metric_fast."""
-        xp0, yp0 = 200.0, 300.0
-        xm, ym = pixel_to_metric_fast(xp0, yp0, 1024, 768, 0.01, 0.01, 0)
-        xp1, yp1 = metric_to_pixel_fast(xm, ym, 1024, 768, 0.01, 0.01, 0)
-        assert math.isclose(xp1, xp0, rel_tol=1e-9)
-        assert math.isclose(yp1, yp0, rel_tol=1e-9)
-
-    def test_chfield_one(self):
-        xp, yp = metric_to_pixel_fast(0.0, 0.0, 1024, 768, 0.01, 0.01, 1)
-        # x_pixel = 0/0.01 + 512 = 512
-        assert math.isclose(xp, 512.0, rel_tol=1e-9)
-        # y_pixel raw = 768/2 - 0/0.01 = 384 → (384 - 1) / 2 = 191.5
-        assert math.isclose(yp, (384.0 - 1.0) * 0.5, rel_tol=1e-9)
-
-    def test_chfield_two(self):
-        xp, yp = metric_to_pixel_fast(0.0, 0.0, 1024, 768, 0.01, 0.01, 2)
-        # y_pixel raw = 384 → 384 * 0.5 = 192
-        assert math.isclose(yp, 192.0, rel_tol=1e-9)
-
-    def test_off_centre(self):
-        xp, yp = metric_to_pixel_fast(1.0, 0.0, 1024, 768, 0.01, 0.01, 0)
-        assert math.isclose(xp, 512.0 + 100.0, rel_tol=1e-9)
-
-
-# ---------------------------------------------------------------------------
-# 9. _metric_to_pixel_out
-# ---------------------------------------------------------------------------
-
-
-class TestMetricToPixelOut:
-    def test_chfield_zero(self):
-        out = np.zeros(2, dtype=np.float64)
-        _metric_to_pixel_out(0.0, 0.0, 1024, 768, 0.01, 0.01, 0, out)
-        assert math.isclose(out[0], 512.0, rel_tol=1e-9)
-        assert math.isclose(out[1], 384.0, rel_tol=1e-9)
-
-    def test_chfield_one(self):
-        out = np.zeros(2, dtype=np.float64)
-        _metric_to_pixel_out(0.0, 0.0, 1024, 768, 0.01, 0.01, 1, out)
-        assert math.isclose(out[1], (384.0 - 1.0) * 0.5, rel_tol=1e-9)
-
-    def test_chfield_two(self):
-        out = np.zeros(2, dtype=np.float64)
-        _metric_to_pixel_out(0.0, 0.0, 1024, 768, 0.01, 0.01, 2, out)
-        assert math.isclose(out[1], 192.0, rel_tol=1e-9)
-
-
-# ---------------------------------------------------------------------------
-# 10. _flat_image_coord_fast
-# ---------------------------------------------------------------------------
-
-
-class TestFlatImageCoordFast:
-    def _pos(self, x=0.0, y=0.0, z=0.0):
-        return np.array([x, y, z], dtype=np.float64)
-
-    def _empty_mmlut(self):
-        return np.array([], dtype=np.float64), np.zeros(3, dtype=np.float64), 0, 0, 1.0
-
-    def _filled_mmlut(self, factor=1.0):
-        """2x2 LUT with constant factor."""
-        data = np.full(4, factor, dtype=np.float64)
-        origin = np.zeros(3, dtype=np.float64)
-        return data, origin, 2, 2, 1000.0
-
-    def test_basic_no_mmlut(self):
-        cal = _make_cal_arr(x0=0.0, y0=0.0, z0=100.0, gz=50.0)
-        pos = self._pos(0.0, 0.0, 0.0)
-        mmlut_data, mmlut_origin, nr, nz, rw = self._empty_mmlut()
-        x, y = _flat_image_coord_fast(pos, cal, mmlut_data, mmlut_origin, nr, nz, rw)
-        assert isinstance(x, float)
-        assert isinstance(y, float)
-
-    def test_with_mmlut_in_bounds(self):
-        """LUT in-bounds path (has_mmlut=True, mmf > 0)."""
-        cal = _make_cal_arr(x0=0.0, y0=0.0, z0=100.0, gz=50.0)
-        pos = self._pos(1.0, 0.0, -10.0)
-        mmlut_data, mmlut_origin, nr, nz, rw = self._filled_mmlut(1.2)
-        x, y = _flat_image_coord_fast(pos, cal, mmlut_data, mmlut_origin, nr, nz, rw)
-        assert isinstance(x, float)
-
-    def test_with_mmlut_zero_factor(self):
-        """LUT path where mmf == 0 → falls back to _multimed_r_nlay_1layer."""
-        cal = _make_cal_arr(x0=0.0, y0=0.0, z0=100.0, gz=50.0)
-        pos = self._pos(1.0, 0.0, -10.0)
-        mmlut_data = np.zeros(4, dtype=np.float64)  # mmf == 0
-        mmlut_origin = np.zeros(3, dtype=np.float64)
-        x, y = _flat_image_coord_fast(pos, cal, mmlut_data, mmlut_origin, 2, 2, 1000.0)
-        assert isinstance(x, float)
-
-    def test_pos_t_0_zero_branch(self):
-        """pos_t_0 == 0 → the s_x branch is skipped."""
-        # Place the point directly along the glass normal from the camera projection
-        cal = _make_cal_arr(x0=0.0, y0=0.0, z0=100.0, gz=50.0)
-        pos = self._pos(0.0, 0.0, 0.0)
-        mmlut_data, mmlut_origin, nr, nz, rw = self._empty_mmlut()
-        x, y = _flat_image_coord_fast(pos, cal, mmlut_data, mmlut_origin, nr, nz, rw)
-        assert isinstance(x, float)
-
-    def test_mmlut_out_of_bounds(self):
-        """LUT v3 > nr*nz → skip LUT, fall back to iterative solver."""
-        cal = _make_cal_arr(x0=0.0, y0=0.0, z0=100.0, gz=50.0)
-        pos = self._pos(500.0, 0.0, -10.0)  # large R
-        # Small LUT so ir > mmlut_nr
-        data = np.ones(4, dtype=np.float64)
-        origin = np.zeros(3, dtype=np.float64)
-        x, y = _flat_image_coord_fast(pos, cal, data, origin, 1, 2, 0.001)
-        assert isinstance(x, float)
-
-    def test_radial_shift_one_fallback(self):
-        """When mmlut lookup gives radial_shift still == 1.0, falls through to
-        _multimed_r_nlay_1layer."""
-        cal = _make_cal_arr(
-            x0=0.0, y0=0.0, z0=100.0, gz=50.0, n1=1.0, n2_0=1.5, n3=1.33, d0=2.0
-        )
-        pos = self._pos(2.0, 1.0, -5.0)
-        mmlut_data, mmlut_origin, nr, nz, rw = self._empty_mmlut()
-        x, y = _flat_image_coord_fast(pos, cal, mmlut_data, mmlut_origin, nr, nz, rw)
-        assert isinstance(x, float)
-
-
-# ---------------------------------------------------------------------------
-# 11. _img_coord_fast
-# ---------------------------------------------------------------------------
-
-
-class TestImgCoordFast:
-    def _empty_mmlut(self):
-        return np.array([], dtype=np.float64), np.zeros(3, dtype=np.float64), 0, 0, 1.0
-
-    def test_r_near_zero_returns_zero(self):
-        """_flat_image_coord_fast returns x≈0, y≈0 → r < 1e-10 → (0,0)."""
-        cal = _make_cal_arr(x0=0.0, y0=0.0, z0=100.0, gz=50.0, xh=0.0, yh=0.0)
-        pos = np.array([0.0, 0.0, 0.0], dtype=np.float64)
-        mmlut_data, mmlut_origin, nr, nz, rw = self._empty_mmlut()
-        xd, yd = _img_coord_fast(pos, cal, mmlut_data, mmlut_origin, nr, nz, rw)
-        assert xd == 0.0 and yd == 0.0
-
-    def test_normal_case(self):
-        cal = _make_cal_arr(
-            x0=0.0,
-            y0=0.0,
-            z0=100.0,
-            gz=50.0,
-            xh=0.0,
-            yh=0.0,
-            k1=0.001,
-            scx=1.0,
-            she=0.0,
-        )
-        pos = np.array([1.0, 2.0, 0.0], dtype=np.float64)
-        mmlut_data, mmlut_origin, nr, nz, rw = self._empty_mmlut()
-        xd, yd = _img_coord_fast(pos, cal, mmlut_data, mmlut_origin, nr, nz, rw)
-        assert isinstance(xd, float)
-        assert isinstance(yd, float)
-
-    def test_with_she_nonzero(self):
-        cal = _make_cal_arr(x0=0.0, y0=0.0, z0=100.0, gz=50.0, she=0.05)
-        pos = np.array([2.0, 1.0, 0.0], dtype=np.float64)
-        mmlut_data, mmlut_origin, nr, nz, rw = self._empty_mmlut()
-        xd, yd = _img_coord_fast(pos, cal, mmlut_data, mmlut_origin, nr, nz, rw)
-        assert isinstance(xd, float)
-
-
-# ---------------------------------------------------------------------------
-# 12. img_coord_batch_fast
-# ---------------------------------------------------------------------------
-
-
-class TestImgCoordBatchFast:
-    def _empty_mmlut(self):
-        return np.array([], dtype=np.float64), np.zeros(3, dtype=np.float64), 0, 0, 1.0
-
-    def test_empty_batch(self):
-        cal = _make_cal_arr()
-        positions = np.empty((0, 3), dtype=np.float64, order="C")
-        mmlut_data, mmlut_origin, nr, nz, rw = self._empty_mmlut()
-        result = img_coord_batch_fast(
-            positions, cal, mmlut_data, mmlut_origin, nr, nz, rw
-        )
-        assert result.shape == (0, 2)
-
-    def test_single_point(self):
-        cal = _make_cal_arr(gz=50.0)
-        positions = np.array([[1.0, 2.0, 0.0]], dtype=np.float64, order="C")
-        mmlut_data, mmlut_origin, nr, nz, rw = self._empty_mmlut()
-        result = img_coord_batch_fast(
-            positions, cal, mmlut_data, mmlut_origin, nr, nz, rw
-        )
-        assert result.shape == (1, 2)
-        assert result.dtype == np.float64
-
-    def test_multiple_points(self):
-        cal = _make_cal_arr(gz=50.0)
-        positions = np.array(
-            [
-                [1.0, 0.0, 0.0],
-                [2.0, 1.0, -5.0],
-                [0.5, -0.5, 3.0],
-            ],
-            dtype=np.float64,
-            order="C",
-        )
-        mmlut_data, mmlut_origin, nr, nz, rw = self._empty_mmlut()
-        result = img_coord_batch_fast(
-            positions, cal, mmlut_data, mmlut_origin, nr, nz, rw
-        )
-        assert result.shape == (3, 2)
-
-
-# ---------------------------------------------------------------------------
-# 13. flat_image_coord_batch_fast
-# ---------------------------------------------------------------------------
-
-
-class TestFlatImageCoordBatchFast:
-    def _empty_mmlut(self):
-        return np.array([], dtype=np.float64), np.zeros(3, dtype=np.float64), 0, 0, 1.0
-
-    def test_empty_batch(self):
-        cal = _make_cal_arr()
-        positions = np.empty((0, 3), dtype=np.float64, order="C")
-        mmlut_data, mmlut_origin, nr, nz, rw = self._empty_mmlut()
-        result = flat_image_coord_batch_fast(
-            positions, cal, mmlut_data, mmlut_origin, nr, nz, rw
-        )
-        assert result.shape == (0, 2)
-
-    def test_single_point(self):
-        cal = _make_cal_arr(gz=50.0)
-        positions = np.array([[0.0, 0.0, 0.0]], dtype=np.float64, order="C")
-        mmlut_data, mmlut_origin, nr, nz, rw = self._empty_mmlut()
-        result = flat_image_coord_batch_fast(
-            positions, cal, mmlut_data, mmlut_origin, nr, nz, rw
-        )
-        assert result.shape == (1, 2)
-
-    def test_multiple_points(self):
-        cal = _make_cal_arr(gz=50.0)
-        positions = np.array(
-            [
-                [1.0, 0.0, 0.0],
-                [0.0, 1.0, 0.0],
-                [2.0, 2.0, -10.0],
-                [0.0, 0.0, 5.0],
-            ],
-            dtype=np.float64,
-            order="C",
-        )
-        mmlut_data, mmlut_origin, nr, nz, rw = self._empty_mmlut()
-        result = flat_image_coord_batch_fast(
-            positions, cal, mmlut_data, mmlut_origin, nr, nz, rw
-        )
-        assert result.shape == (4, 2)
-
-
 # ---------------------------------------------------------------------------
 # 14. _candsearch_in_pix_rest_nogil
 # ---------------------------------------------------------------------------
@@ -1568,130 +1182,11 @@ def test_constants():
 # ---------------------------------------------------------------------------
 
 
-def test_pixel_metric_pixel_roundtrip():
-    """metric_to_pixel(pixel_to_metric(px, py)) == (px, py) for chfield 0."""
-    for chfield in [0]:
-        px0, py0 = 300.0, 200.0
-        xm, ym = pixel_to_metric_fast(px0, py0, 1024, 768, 0.01, 0.01, chfield)
-        px1, py1 = metric_to_pixel_fast(xm, ym, 1024, 768, 0.01, 0.01, chfield)
-        assert math.isclose(px1, px0, rel_tol=1e-9)
-        assert math.isclose(py1, py0, rel_tol=1e-9)
-
-
 # ---------------------------------------------------------------------------
 # 19. dist_to_flat / _dist_to_flat_out consistency
 # ---------------------------------------------------------------------------
 
 
-@_needs_pure_python
-def test_dist_to_flat_fast_and_out_agree():
-    """fast and _out variants give same result."""
-    args = (2.5, -1.0, 0.1, -0.2, 0.001, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1e-8)
-    x1, y1 = dist_to_flat_fast(*args)
-    out = np.zeros(2, dtype=np.float64)
-    _dist_to_flat_out(*args, out)
-    assert math.isclose(x1, out[0], rel_tol=1e-12)
-    assert math.isclose(y1, out[1], rel_tol=1e-12)
-
-
-# ---------------------------------------------------------------------------
-# 20. _pixel_to_metric_out / pixel_to_metric_fast consistency
-# ---------------------------------------------------------------------------
-
-
-def test_pixel_to_metric_out_and_fast_agree():
-    """_out and fast variants give same result for all chfields."""
-    for chfield in [0, 1, 2]:
-        x_f, y_f = pixel_to_metric_fast(400.0, 300.0, 1024, 768, 0.01, 0.01, chfield)
-        out = np.zeros(2, dtype=np.float64)
-        _pixel_to_metric_out(400.0, 300.0, 1024, 768, 0.01, 0.01, chfield, out)
-        assert math.isclose(x_f, out[0], rel_tol=1e-12)
-        assert math.isclose(y_f, out[1], rel_tol=1e-12)
-
-
-# ---------------------------------------------------------------------------
-# 21. _metric_to_pixel_out / metric_to_pixel_fast consistency
-# ---------------------------------------------------------------------------
-
-
-def test_metric_to_pixel_out_and_fast_agree():
-    """_out and fast variants give same result for all chfields."""
-    for chfield in [0, 1, 2]:
-        xp_f, yp_f = metric_to_pixel_fast(0.5, -0.3, 1024, 768, 0.01, 0.01, chfield)
-        out = np.zeros(2, dtype=np.float64)
-        _metric_to_pixel_out(0.5, -0.3, 1024, 768, 0.01, 0.01, chfield, out)
-        assert math.isclose(xp_f, out[0], rel_tol=1e-12)
-        assert math.isclose(yp_f, out[1], rel_tol=1e-12)
-
-
-# ---------------------------------------------------------------------------
-# 22. img_coord_batch_fast vs _img_coord_fast element-wise
-# ---------------------------------------------------------------------------
-
-
-def test_img_coord_batch_matches_elementwise():
-    """batch result == repeated scalar calls."""
-    cal = _make_cal_arr(gz=50.0, k1=0.001)
-    mmlut_data = np.array([], dtype=np.float64)
-    mmlut_origin = np.zeros(3, dtype=np.float64)
-    positions = np.array(
-        [
-            [1.0, 0.0, 0.0],
-            [0.0, 2.0, -5.0],
-        ],
-        dtype=np.float64,
-        order="C",
-    )
-    result = img_coord_batch_fast(positions, cal, mmlut_data, mmlut_origin, 0, 0, 1.0)
-    for i in range(len(positions)):
-        xi, yi = _img_coord_fast(positions[i], cal, mmlut_data, mmlut_origin, 0, 0, 1.0)
-        assert math.isclose(result[i, 0], xi, rel_tol=1e-12)
-        assert math.isclose(result[i, 1], yi, rel_tol=1e-12)
-
-
-# ---------------------------------------------------------------------------
-# 23. flat_image_coord_batch_fast vs _flat_image_coord_fast element-wise
-# ---------------------------------------------------------------------------
-
-
-def test_flat_image_coord_batch_matches_elementwise():
-    """batch result == repeated scalar calls."""
-    cal = _make_cal_arr(gz=50.0)
-    mmlut_data = np.array([], dtype=np.float64)
-    mmlut_origin = np.zeros(3, dtype=np.float64)
-    positions = np.array(
-        [
-            [0.5, -0.5, 0.0],
-            [2.0, 1.0, -3.0],
-        ],
-        dtype=np.float64,
-        order="C",
-    )
-    result = flat_image_coord_batch_fast(
-        positions, cal, mmlut_data, mmlut_origin, 0, 0, 1.0
-    )
-    for i in range(len(positions)):
-        xi, yi = _flat_image_coord_fast(
-            positions[i], cal, mmlut_data, mmlut_origin, 0, 0, 1.0
-        )
-        assert math.isclose(result[i, 0], xi, rel_tol=1e-12)
-        assert math.isclose(result[i, 1], yi, rel_tol=1e-12)
-
-
-# ---------------------------------------------------------------------------
-# 24. dist_to_flat_fast / _dist_to_flat_out loop-exhaustion branches (596->618, 665->680)
-# ---------------------------------------------------------------------------
-
-
-def test_dist_to_flat_fast_loop_exhaustion():
-    """tol=0.0 prevents the break from firing → all 50 iterations run (596->618 branch)."""
-    x, y = dist_to_flat_fast(
-        1.0, 1.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0
-    )
-    # Result exists even without convergence
-    assert isinstance(x, float)
-    assert isinstance(y, float)
-
 
 @_needs_pure_python
 def test_dist_to_flat_out_loop_exhaustion():
@@ -1733,30 +1228,6 @@ def test_candsearch_xmax_ymax_clamp():
     assert result == 0
 
 
-# ---------------------------------------------------------------------------
-# 26. _flat_image_coord_fast: branch 1262->1271 (v3 > mmlut_nr*mmlut_nz)
-# ---------------------------------------------------------------------------
-
-
-def test_flat_image_coord_fast_lut_boundary_branch():
-    """ir == mmlut_nr AND iz in [0,nz] → v3 > mmlut_nr*mmlut_nz → inner if False → 1262->1271.
-
-    Geometry (x0=0, y0=0, z0=100, gz=50):
-      With pos=[2.5, 0, 51] and rw=1.0, nr=2, nz=2:
-        pos_t_0 = 2.5 → R=2.5 → ir=2 = nr=2 (outer condition True: 2<=2)
-        dist_point_glas = 51-50 = 1 → iz=1 (in [0,2])
-        v3 = 2*2+1+2+1 = 8 > 4 = nr*nz → inner condition False → 1262->1271
-    """
-    cal = _make_cal_arr(x0=0.0, y0=0.0, z0=100.0, gz=50.0, d0=0.0)
-    pos = np.array([2.5, 0.0, 51.0], dtype=np.float64)
-    nr, nz = 2, 2
-    mmlut_data = np.ones(nr * nz, dtype=np.float64)  # non-empty LUT
-    mmlut_origin = np.zeros(3, dtype=np.float64)
-    rw = 1.0
-    x, y = _flat_image_coord_fast(pos, cal, mmlut_data, mmlut_origin, nr, nz, rw)
-    assert isinstance(x, float)
-
-
 # ---------------------------------------------------------------------------
 # 27. assess_new_position_fast: use_proj=False path (lines 765-783)
 # ---------------------------------------------------------------------------