diff --git a/common/ParallelFor.h b/common/ParallelFor.h
new file mode 100644
index 00000000..767a070f
--- /dev/null
+++ b/common/ParallelFor.h
@@ -0,0 +1,55 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <algorithm>
+#include <atomic>
+#include <future>
+#include <vector>
+
+// Two shapes of "run this over a range on several threads", used by the analysis code. Both take the
+// worker count from the caller rather than asking the hardware, so a run that was told how many
+// threads to use keeps to it.
+
+// Chunked: each worker gets one contiguous [lo, hi) range and there is no per-item synchronisation.
+// Right for millions of cheap uniform items - the CPU stand-in for a flat CUDA grid-stride kernel.
+// The split is fixed and deterministic, so a pass whose per-element work is independent gives the
+// same answer as the serial loop, bit for bit.
+template <typename Fn>
+void ParallelChunks(int n, size_t nthreads, Fn fn) {
+    if (n <= 0) return;
+    const int nt = static_cast<int>(std::max<size_t>(1, std::min(nthreads, static_cast<size_t>(n))));
+    if (nt == 1) { fn(0, n); return; }
+    const int chunk = (n + nt - 1) / nt;
+    std::vector<std::future<void>> futures;
+    futures.reserve(nt);
+    for (int t = 0; t < nt; t++) {
+        const int lo = t * chunk, hi = std::min(n, lo + chunk);
+        if (lo >= hi) break;
+        futures.emplace_back(std::async(std::launch::async, [&fn, lo, hi] { fn(lo, hi); }));
+    }
+    for (auto &f : futures) f.get();
+}
+
+// Work-stealing per item, off a shared atomic counter: one atomic per item, so use it only where the
+// per-item work is heavy and uneven (per-frame fits, per-ring selections) and the atomic amortises.
+// For millions of tiny uniform items a per-item atomic is pure contention - use ParallelChunks.
+template <typename Fn>
+void ParallelFor(int n, size_t nthreads, Fn fn) {
+    if (n <= 0) return;
+    if (nthreads <= 1 || n == 1) {
+        for (int i = 0; i < n; i++) fn(i);
+        return;
+    }
+    const size_t local = std::min(nthreads, static_cast<size_t>(n));
+    std::atomic<int> next = 0;
+    std::vector<std::future<void>> futures;
+    futures.reserve(local);
+    for (size_t t = 0; t < local; t++)
+        futures.emplace_back(std::async(std::launch::async, [&] {
+            for (int i = next.fetch_add(1); i < n; i = next.fetch_add(1))
+                fn(i);
+        }));
+    for (auto &f : futures) f.get();
+}
diff --git a/image_analysis/beam_stop/ShadowFinder.cpp b/image_analysis/beam_stop/ShadowFinder.cpp
index 233a8d74..29eb36bc 100644
--- a/image_analysis/beam_stop/ShadowFinder.cpp
+++ b/image_analysis/beam_stop/ShadowFinder.cpp
@@ -2,6 +2,7 @@
 // SPDX-License-Identifier: GPL-3.0-only
 
 #include "ShadowFinder.h"
+#include "../../common/ParallelFor.h"
 
 #include <algorithm>
 #include <cmath>
@@ -47,45 +48,6 @@ constexpr int MIN_RING_PIXELS = 32;
 // at GetMask() time; the BFS forms keep them O(pixels) rather than O(pixels * radius).
 namespace {
 
-// Run fn(lo, hi) over contiguous slices of [0, n), one per worker. For the per-pixel passes: the
-// same self-load-balancing shape the rest of the codebase uses, kept local because this is the only
-// other file that wants it.
-template <typename Fn>
-void ParallelChunks(int n, size_t nthreads, Fn fn) {
-    if (n <= 0) return;
-    const int nt = static_cast<int>(std::max<size_t>(1, std::min(nthreads, static_cast<size_t>(n))));
-    if (nt == 1) { fn(0, n); return; }
-    const int chunk = (n + nt - 1) / nt;
-    std::vector<std::future<void>> futures;
-    futures.reserve(nt);
-    for (int t = 0; t < nt; t++) {
-        const int lo = t * chunk, hi = std::min(n, lo + chunk);
-        if (lo >= hi) break;
-        futures.emplace_back(std::async(std::launch::async, [&fn, lo, hi] { fn(lo, hi); }));
-    }
-    for (auto &f : futures) f.get();
-}
-
-// Work-stealing per-item parallel, for items whose cost is very uneven - the rings, whose sizes go
-// as the circumference and then fall away at the detector corners.
-template <typename Fn>
-void ParallelFor(int n, size_t nthreads, Fn fn) {
-    if (n <= 0) return;
-    if (nthreads <= 1 || n == 1) {
-        for (int i = 0; i < n; i++) fn(i);
-        return;
-    }
-    const size_t local = std::min(nthreads, static_cast<size_t>(n));
-    std::atomic<int> next = 0;
-    std::vector<std::future<void>> futures;
-    futures.reserve(local);
-    for (size_t t = 0; t < local; t++)
-        futures.emplace_back(std::async(std::launch::async, [&] {
-            for (int i = next.fetch_add(1); i < n; i = next.fetch_add(1)) fn(i);
-        }));
-    for (auto &f : futures) f.get();
-}
-
 // 8-connected dilation by `r` pixels, i.e. every pixel within Chebyshev distance r of a set one.
 //
 // This was a multi-source BFS, which is what the distance is defined by - but on a full rectangle
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 790296a7..9892c0b5 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -51,6 +51,7 @@ ADD_EXECUTABLE(jfjoch_test
         ZMQMetadataSocketTest.cpp
         JFJochReaderTest.cpp
         ShadowFinderTest.cpp
+        ParallelForTest.cpp
         RugnuxTest.cpp
         ResultReportTest.cpp
         RugnuxLargeTest.cpp
diff --git a/tests/ParallelForTest.cpp b/tests/ParallelForTest.cpp
new file mode 100644
index 00000000..cec2988a
--- /dev/null
+++ b/tests/ParallelForTest.cpp
@@ -0,0 +1,127 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+
+#include <algorithm>
+#include <mutex>
+#include <numeric>
+#include <stdexcept>
+#include <vector>
+
+#include "../common/ParallelFor.h"
+
+namespace {
+    const std::vector<size_t> THREAD_COUNTS = {0, 1, 2, 3, 8, 64};
+    const std::vector<int> SIZES = {0, 1, 2, 7, 64, 1000};
+}
+
+// The split is fixed and contiguous, which is what lets a pass whose per-element work is independent
+// give the serial answer bit for bit. Every element has to land in exactly one slice, the slices have
+// to tile [0, n) in order, and there must never be an empty one.
+TEST_CASE("ParallelChunks_SlicesTileTheRange", "[ParallelFor]") {
+    for (int n: SIZES) {
+        for (size_t nthreads: THREAD_COUNTS) {
+            CAPTURE(n, nthreads);
+            std::mutex m;
+            std::vector<std::pair<int, int>> slices;
+            std::vector<int> visits(static_cast<size_t>(std::max(n, 1)), 0);
+
+            ParallelChunks(n, nthreads, [&](int lo, int hi) {
+                {
+                    std::scoped_lock lock(m);
+                    slices.emplace_back(lo, hi);
+                }
+                for (int i = lo; i < hi; i++)
+                    visits[static_cast<size_t>(i)]++;   // slices are disjoint, so no lock needed
+            });
+
+            std::sort(slices.begin(), slices.end());
+            int expected_lo = 0;
+            for (const auto &[lo, hi]: slices) {
+                CHECK(lo == expected_lo);
+                CHECK(hi > lo);                          // never an empty slice
+                expected_lo = hi;
+            }
+            CHECK(expected_lo == n);                     // and they reach the end
+            CHECK(slices.size() <= static_cast<size_t>(std::max(n, 0)));
+            for (int i = 0; i < n; i++)
+                CHECK(visits[static_cast<size_t>(i)] == 1);
+        }
+    }
+}
+
+// Work stealing, so a worker takes whatever is next rather than a fixed slice - but every item still
+// has to be done exactly once, whatever order they come out in.
+TEST_CASE("ParallelFor_VisitsEveryItemOnce", "[ParallelFor]") {
+    for (int n: SIZES) {
+        for (size_t nthreads: THREAD_COUNTS) {
+            CAPTURE(n, nthreads);
+            std::vector<std::atomic<int>> visits(static_cast<size_t>(std::max(n, 1)));
+            for (auto &v: visits)
+                v.store(0);
+
+            ParallelFor(n, nthreads, [&](int i) { visits[static_cast<size_t>(i)].fetch_add(1); });
+
+            for (int i = 0; i < n; i++)
+                CHECK(visits[static_cast<size_t>(i)].load() == 1);
+        }
+    }
+}
+
+// One thread means the caller's loop, in order - the property the serial fallbacks rely on.
+TEST_CASE("ParallelFor_IsSerialAndInOrderForOneThread", "[ParallelFor]") {
+    for (size_t nthreads: {size_t{0}, size_t{1}}) {
+        CAPTURE(nthreads);
+        std::vector<int> order;
+        ParallelFor(16, nthreads, [&](int i) { order.push_back(i); });
+        std::vector<int> expected(16);
+        std::iota(expected.begin(), expected.end(), 0);
+        CHECK(order == expected);
+    }
+}
+
+// Splitting the work must not change the answer. Each element is written by exactly one worker, so
+// the result has to match the serial loop element for element, at every thread count.
+TEST_CASE("ParallelFor_SplitDoesNotChangeTheResult", "[ParallelFor]") {
+    constexpr int N = 1000;
+    std::vector<double> serial(N);
+    for (int i = 0; i < N; i++)
+        serial[static_cast<size_t>(i)] = std::sin(i * 0.001) * 1e6 + i;
+
+    for (size_t nthreads: THREAD_COUNTS) {
+        CAPTURE(nthreads);
+        std::vector<double> chunked(N, 0.0), stolen(N, 0.0);
+        ParallelChunks(N, nthreads, [&](int lo, int hi) {
+            for (int i = lo; i < hi; i++)
+                chunked[static_cast<size_t>(i)] = std::sin(i * 0.001) * 1e6 + i;
+        });
+        ParallelFor(N, nthreads, [&](int i) {
+            stolen[static_cast<size_t>(i)] = std::sin(i * 0.001) * 1e6 + i;
+        });
+        CHECK(chunked == serial);       // bit for bit, not approximately
+        CHECK(stolen == serial);
+    }
+}
+
+// A negative or zero count is a no-op rather than an error - callers pass a computed size.
+TEST_CASE("ParallelFor_DoesNothingForAnEmptyRange", "[ParallelFor]") {
+    int calls = 0;
+    for (int n: {0, -1, -1000}) {
+        ParallelChunks(n, 8, [&](int, int) { calls++; });
+        ParallelFor(n, 8, [&](int) { calls++; });
+    }
+    CHECK(calls == 0);
+}
+
+// An exception thrown in a worker reaches the caller rather than terminating: the futures are waited
+// on, so the other workers finish first and only then does it propagate.
+TEST_CASE("ParallelFor_PropagatesAnException", "[ParallelFor]") {
+    CHECK_THROWS_AS(ParallelChunks(64, 4, [](int lo, int) {
+        if (lo == 0) throw std::runtime_error("from a chunk");
+    }), std::runtime_error);
+
+    CHECK_THROWS_AS(ParallelFor(64, 4, [](int i) {
+        if (i == 0) throw std::runtime_error("from an item");
+    }), std::runtime_error);
+}