diff --git a/image_analysis/beam_stop/ShadowFinder.cpp b/image_analysis/beam_stop/ShadowFinder.cpp
index 5ed5799d..233a8d74 100644
--- a/image_analysis/beam_stop/ShadowFinder.cpp
+++ b/image_analysis/beam_stop/ShadowFinder.cpp
@@ -6,6 +6,9 @@
 #include <algorithm>
 #include <cmath>
 #include <limits>
+#include <atomic>
+#include <future>
+#include <thread>
 #include <queue>
 #include <type_traits>
 
@@ -44,43 +47,113 @@ constexpr int MIN_RING_PIXELS = 32;
 // at GetMask() time; the BFS forms keep them O(pixels) rather than O(pixels * radius).
 namespace {
 
-// 8-connected dilation by `r` pixels (Chebyshev), via a multi-source BFS.
-std::vector<char> dilate(const std::vector<char> &in, int W, int H, int r) {
+// 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
+// with no obstacles the 8-connected graph distance IS the Chebyshev distance (a path stepping
+// towards the target never has to leave the frame), so the result is a dilation by the (2r+1) square
+// clipped to the frame. A square is the product of a horizontal and a vertical segment and max is
+// associative, so it separates into one pass along x and one along y - O(1) per pixel whatever r is,
+// no queue, and no 4-bytes-per-pixel distance array. The values are 0/1, so the running "is any set"
+// is just a count of set pixels in the window.
+std::vector<char> dilate(const std::vector<char> &in, int W, int H, int r, size_t nthreads) {
     if (r <= 0)
         return in;
-    std::vector<int> dist(in.size(), -1);
-    std::queue<int> q;
-    for (size_t i = 0; i < in.size(); i++)
-        if (in[i]) { dist[i] = 0; q.push(static_cast<int>(i)); }
-    while (!q.empty()) {
-        const int i = q.front(); q.pop();
-        if (dist[i] >= r)
-            continue;
-        const int y = i / W, x = i % W;
-        for (int dy = -1; dy <= 1; dy++)
-            for (int dx = -1; dx <= 1; dx++) {
-                const int yy = y + dy, xx = x + dx;
-                if (yy < 0 || yy >= H || xx < 0 || xx >= W)
-                    continue;
-                const int j = yy * W + xx;
-                if (dist[j] < 0) { dist[j] = dist[i] + 1; q.push(j); }
+    std::vector<char> tmp(in.size()), out(in.size());
+
+    ParallelChunks(H, nthreads, [&](int ylo, int yhi) {
+        for (int y = ylo; y < yhi; y++) {
+            const char *src = in.data() + static_cast<size_t>(y) * W;
+            char *dst = tmp.data() + static_cast<size_t>(y) * W;
+            int count = 0;
+            for (int x = 0; x <= std::min(r, W - 1); x++)
+                count += src[x];
+            for (int x = 0; x < W; x++) {
+                dst[x] = count > 0;
+                if (x + r + 1 < W) count += src[x + r + 1];
+                if (x - r >= 0)    count -= src[x - r];
             }
-    }
-    std::vector<char> out(in.size());
-    for (size_t i = 0; i < out.size(); i++)
-        out[i] = (dist[i] >= 0) ? 1 : 0;
+        }
+    });
+
+    // The vertical pass walks a column, which strides by a whole row. Take a strip of columns at a
+    // time so each row it touches is read contiguously instead of one byte per cache line.
+    constexpr int STRIP = 64;
+    ParallelChunks((W + STRIP - 1) / STRIP, nthreads, [&](int slo, int shi) {
+        std::vector<int> count(STRIP);
+        for (int st = slo; st < shi; st++) {
+            const int x0 = st * STRIP, xn = std::min(STRIP, W - x0);
+            std::fill(count.begin(), count.begin() + xn, 0);
+            for (int y = 0; y <= std::min(r, H - 1); y++)
+                for (int i = 0; i < xn; i++)
+                    count[i] += tmp[static_cast<size_t>(y) * W + x0 + i];
+            for (int y = 0; y < H; y++) {
+                for (int i = 0; i < xn; i++)
+                    out[static_cast<size_t>(y) * W + x0 + i] = count[i] > 0;
+                if (y + r + 1 < H)
+                    for (int i = 0; i < xn; i++)
+                        count[i] += tmp[static_cast<size_t>(y + r + 1) * W + x0 + i];
+                if (y - r >= 0)
+                    for (int i = 0; i < xn; i++)
+                        count[i] -= tmp[static_cast<size_t>(y - r) * W + x0 + i];
+            }
+        }
+    });
     return out;
 }
 
-// Erosion by `r` = dilation of the complement; outside the frame counts as complement.
-std::vector<char> erode(const std::vector<char> &in, int W, int H, int r) {
+// Erosion by `r` = dilation of the complement. The dilation is clipped to the frame and cannot seed
+// outside it, so outside the frame contributes nothing - a pixel within r of the edge is eroded only
+// by what the frame actually holds.
+std::vector<char> erode(const std::vector<char> &in, int W, int H, int r, size_t nthreads) {
     std::vector<char> comp(in.size());
-    for (size_t i = 0; i < in.size(); i++)
-        comp[i] = !in[i];
-    const auto grown = dilate(comp, W, H, r);
+    ParallelChunks(static_cast<int>(in.size()), nthreads, [&](int lo, int hi) {
+        for (int i = lo; i < hi; i++) comp[i] = !in[i];
+    });
+    const auto grown = dilate(comp, W, H, r, nthreads);
     std::vector<char> out(in.size());
-    for (size_t i = 0; i < out.size(); i++)
-        out[i] = !grown[i];
+    ParallelChunks(static_cast<int>(out.size()), nthreads, [&](int lo, int hi) {
+        for (int i = lo; i < hi; i++) out[i] = !grown[i];
+    });
     return out;
 }
 
@@ -106,75 +179,140 @@ std::vector<char> flood(const std::vector<char> &passable, int W, int H, const s
 }
 
 // Fill holes: background not reachable from the image border becomes region.
+//
+// The flood is run over the bounding box of `region` grown by one, not the whole detector. Outside
+// that box every pixel is background and the box's surrounding ring is background too, so the whole
+// outside is one border-connected component: a background pixel inside the box is border-connected
+// exactly when it reaches the ring. The beam stop occupies a small part of a detector, so this is
+// the same answer over a fraction of the pixels.
 std::vector<char> fill_holes(const std::vector<char> &region, int W, int H) {
-    std::vector<char> bg_visited(region.size(), 0);
-    std::queue<int> q;
-    auto push = [&](int i) { if (!region[i] && !bg_visited[i]) { bg_visited[i] = 1; q.push(i); } };
-    for (int x = 0; x < W; x++) { push(x); push((H - 1) * W + x); }
-    for (int y = 0; y < H; y++) { push(y * W); push(y * W + W - 1); }
+    int x0 = W, x1 = -1, y0 = H, y1 = -1;
+    for (int y = 0; y < H; y++)
+        for (int x = 0; x < W; x++)
+            if (region[static_cast<size_t>(y) * W + x]) {
+                x0 = std::min(x0, x); x1 = std::max(x1, x);
+                y0 = std::min(y0, y); y1 = std::max(y1, y);
+            }
+    if (x1 < 0)
+        return region;   // nothing to enclose
+    x0 = std::max(0, x0 - 1); x1 = std::min(W - 1, x1 + 1);
+    y0 = std::max(0, y0 - 1); y1 = std::min(H - 1, y1 + 1);
+
+    const int BW = x1 - x0 + 1, BH = y1 - y0 + 1;
+    std::vector<char> bg_visited(static_cast<size_t>(BW) * BH, 0);
+    std::queue<int> q;   // indices into the box
+    auto push = [&](int bx, int by) {
+        const int j = by * BW + bx;
+        if (!region[static_cast<size_t>(by + y0) * W + bx + x0] && !bg_visited[j]) {
+            bg_visited[j] = 1; q.push(j);
+        }
+    };
+    for (int bx = 0; bx < BW; bx++) { push(bx, 0); push(bx, BH - 1); }
+    for (int by = 0; by < BH; by++) { push(0, by); push(BW - 1, by); }
     while (!q.empty()) {
         const int i = q.front(); q.pop();
-        const int y = i / W, x = i % W;
+        const int by = i / BW, bx = i % BW;
         for (int dy = -1; dy <= 1; dy++)
             for (int dx = -1; dx <= 1; dx++) {
-                const int yy = y + dy, xx = x + dx;
-                if (yy < 0 || yy >= H || xx < 0 || xx >= W)
+                const int yy = by + dy, xx = bx + dx;
+                if (yy < 0 || yy >= BH || xx < 0 || xx >= BW)
                     continue;
-                const int j = yy * W + xx;
-                if (!region[j] && !bg_visited[j]) { bg_visited[j] = 1; q.push(j); }
+                const int j = yy * BW + xx;
+                if (!region[static_cast<size_t>(yy + y0) * W + xx + x0] && !bg_visited[j]) {
+                    bg_visited[j] = 1; q.push(j);
+                }
             }
     }
+
     std::vector<char> out = region;
-    for (size_t i = 0; i < out.size(); i++)
-        if (!region[i] && !bg_visited[i])
-            out[i] = 1;
+    for (int by = 0; by < BH; by++)
+        for (int bx = 0; bx < BW; bx++) {
+            const size_t i = static_cast<size_t>(by + y0) * W + bx + x0;
+            if (!region[i] && !bg_visited[by * BW + bx])
+                out[i] = 1;
+        }
     return out;
 }
 
 // Sum of `in` over the k x k box centred on each pixel, zero outside the frame.
-std::vector<double> box_sum(const std::vector<double> &in, int W, int H, int k) {
+//
+// Each row's and each column's running sum keeps exactly the terms it had, in the order it had them,
+// so the floating-point form rounds identically - only the traversal changes. Walking one column at
+// a time strides a whole row per step and misses on every access, so the vertical pass takes a strip
+// of columns together and reads each row it touches contiguously. Rows, and strips, are independent.
+template <typename T>
+std::vector<T> box_sum(const std::vector<T> &in, int W, int H, int k, size_t nthreads) {
     const int half = k / 2;
-    std::vector<double> row(in.size(), 0.0), out(in.size(), 0.0);
-    for (int y = 0; y < H; y++) {
-        double s = 0;
-        for (int x = 0; x <= std::min(half, W - 1); x++)
-            s += in[y * W + x];
-        for (int x = 0; x < W; x++) {
-            row[y * W + x] = s;
-            if (x + half + 1 < W) s += in[y * W + x + half + 1];
-            if (x - half >= 0)    s -= in[y * W + x - half];
+    std::vector<T> row(in.size()), out(in.size());
+
+    ParallelChunks(H, nthreads, [&](int ylo, int yhi) {
+        for (int y = ylo; y < yhi; y++) {
+            const T *src = in.data() + static_cast<size_t>(y) * W;
+            T *dst = row.data() + static_cast<size_t>(y) * W;
+            T s = 0;
+            for (int x = 0; x <= std::min(half, W - 1); x++)
+                s += src[x];
+            for (int x = 0; x < W; x++) {
+                dst[x] = s;
+                if (x + half + 1 < W) s += src[x + half + 1];
+                if (x - half >= 0)    s -= src[x - half];
+            }
         }
-    }
-    for (int x = 0; x < W; x++) {
-        double s = 0;
-        for (int y = 0; y <= std::min(half, H - 1); y++)
-            s += row[y * W + x];
-        for (int y = 0; y < H; y++) {
-            out[y * W + x] = s;
-            if (y + half + 1 < H) s += row[(y + half + 1) * W + x];
-            if (y - half >= 0)    s -= row[(y - half) * W + x];
+    });
+
+    constexpr int STRIP = 64;
+    ParallelChunks((W + STRIP - 1) / STRIP, nthreads, [&](int slo, int shi) {
+        std::vector<T> s(STRIP);
+        for (int st = slo; st < shi; st++) {
+            const int x0 = st * STRIP, xn = std::min(STRIP, W - x0);
+            std::fill(s.begin(), s.begin() + xn, T{0});
+            for (int y = 0; y <= std::min(half, H - 1); y++)
+                for (int i = 0; i < xn; i++)
+                    s[i] += row[static_cast<size_t>(y) * W + x0 + i];
+            for (int y = 0; y < H; y++) {
+                for (int i = 0; i < xn; i++)
+                    out[static_cast<size_t>(y) * W + x0 + i] = s[i];
+                if (y + half + 1 < H)
+                    for (int i = 0; i < xn; i++)
+                        s[i] += row[static_cast<size_t>(y + half + 1) * W + x0 + i];
+                if (y - half >= 0)
+                    for (int i = 0; i < xn; i++)
+                        s[i] -= row[static_cast<size_t>(y - half) * W + x0 + i];
+            }
         }
-    }
+    });
     return out;
 }
 
-// Median of `values` per integer radius, over the pixels flagged in `use`.
-std::vector<float> ring_median(const std::vector<float> &values, const std::vector<char> &use,
-                               const std::vector<int> &radius, int max_radius) {
-    std::vector<std::vector<float>> bins(max_radius + 1);
+// The values of each ring, laid end to end, with the ring's slice given by offset[r]..offset[r+1].
+// radius and pooled do not change over the three baseline iterations, so this is built once and
+// every iteration is an order statistic of the same, already-sorted, ring.
+struct RingValues {
+    std::vector<float> values;
+    std::vector<int>   offset;
+};
+
+RingValues bin_by_ring(const std::vector<float> &values, const std::vector<char> &valid,
+                       const std::vector<int> &radius, int max_radius, size_t nthreads) {
+    RingValues rv;
+    rv.offset.assign(max_radius + 2, 0);
     for (size_t i = 0; i < values.size(); i++)
-        if (use[i])
-            bins[radius[i]].push_back(values[i]);
-    std::vector<float> out(max_radius + 1, 0.0f);
-    for (int r = 0; r <= max_radius; r++) {
-        auto &b = bins[r];
-        if (!b.empty()) {
-            const size_t k = b.size() / 2;
-            std::nth_element(b.begin(), b.begin() + k, b.end());
-            out[r] = b[k];
-        }
-    }
-    return out;
+        if (valid[i])
+            rv.offset[radius[i] + 1]++;
+    for (int r = 0; r <= max_radius; r++)
+        rv.offset[r + 1] += rv.offset[r];
+
+    rv.values.resize(rv.offset[max_radius + 1]);
+    std::vector<int> cursor(rv.offset.begin(), rv.offset.end() - 1);
+    for (size_t i = 0; i < values.size(); i++)
+        if (valid[i])
+            rv.values[cursor[radius[i]]++] = values[i];
+
+    // Sorted once; the three iterations then only pick a rank and count a prefix.
+    ParallelFor(max_radius + 1, nthreads, [&](int r) {
+        std::sort(rv.values.begin() + rv.offset[r], rv.values.begin() + rv.offset[r + 1]);
+    });
+    return rv;
 }
 
 }  // namespace
@@ -212,17 +350,34 @@ ShadowFinder::Projection ShadowFinder::Reduce() const {
     out.max_value.assign(npixels, 0);
     out.sum_value.assign(npixels, 0);
     out.valid_count.assign(npixels, 0);
-    for (const auto &p : shards) {
+    for (const auto &p : shards)
         out.frames += p.frames;
-        for (size_t i = 0; i < npixels; i++) {
-            if (p.valid_count[i] == 0)
-                continue;
-            if (out.valid_count[i] == 0 || p.max_value[i] > out.max_value[i])
-                out.max_value[i] = p.max_value[i];
-            out.sum_value[i] += p.sum_value[i];
-            out.valid_count[i] += p.valid_count[i];
-        }
+
+    // Each worker owns a slice of the pixels and folds every shard into it. The sums and counts are
+    // integers and a pixel is touched by one worker only, so the result is the same as folding them
+    // one shard at a time on one thread - this is several hundred megabytes per shard and is limited
+    // by memory rather than by arithmetic.
+    const size_t nthreads = std::max<size_t>(1, std::min<size_t>(std::thread::hardware_concurrency(),
+                                                                 shards.size() * 2));
+    const size_t chunk = (npixels + nthreads - 1) / nthreads;
+    std::vector<std::future<void>> futures;
+    futures.reserve(nthreads);
+    for (size_t t = 0; t < nthreads; t++) {
+        const size_t lo = t * chunk, hi = std::min(npixels, lo + chunk);
+        if (lo >= hi) break;
+        futures.emplace_back(std::async(std::launch::async, [&, lo, hi] {
+            for (const auto &p : shards)
+                for (size_t i = lo; i < hi; i++) {
+                    if (p.valid_count[i] == 0)
+                        continue;
+                    if (out.valid_count[i] == 0 || p.max_value[i] > out.max_value[i])
+                        out.max_value[i] = p.max_value[i];
+                    out.sum_value[i] += p.sum_value[i];
+                    out.valid_count[i] += p.valid_count[i];
+                }
+        }));
     }
+    for (auto &f : futures) f.get();
     return out;
 }
 
@@ -295,8 +450,10 @@ std::vector<float> ShadowFinder::GetMeanProjection() const {
     return mean;
 }
 
-std::vector<uint32_t> ShadowFinder::GetMask() const {
+std::vector<uint32_t> ShadowFinder::GetMask(size_t nthreads) const {
     std::unique_lock ul(m);
+    if (nthreads == 0)
+        nthreads = std::max(1u, std::thread::hardware_concurrency());
     const Projection p = Reduce();
     const auto &max_value = p.max_value;
     const auto &sum_value = p.sum_value;
@@ -314,51 +471,90 @@ std::vector<uint32_t> ShadowFinder::GetMask() const {
     std::vector<float> mean(n_pixels, 0.0f);
     std::vector<char>  valid(n_pixels, 0);
     std::vector<int>   radius(n_pixels, 0);
-    int max_radius = 0;
-    for (int y = 0; y < H; y++)
-        for (int x = 0; x < W; x++) {
-            const int i = y * W + x;
-            if (valid_count[i] > 0 && pixel_mask[i] == 0) {
-                mean[i] = static_cast<float>(static_cast<double>(sum_value[i]) / valid_count[i]);
-                valid[i] = 1;
+    std::atomic<int> max_radius_atomic{0};
+    ParallelChunks(H, nthreads, [&](int ylo, int yhi) {
+        int local_max = 0;
+        for (int y = ylo; y < yhi; y++)
+            for (int x = 0; x < W; x++) {
+                const int i = y * W + x;
+                if (valid_count[i] > 0 && pixel_mask[i] == 0) {
+                    mean[i] = static_cast<float>(static_cast<double>(sum_value[i]) / valid_count[i]);
+                    valid[i] = 1;
+                }
+                const float dx = x - beam_x, dy = y - beam_y;
+                radius[i] = static_cast<int>(std::lround(std::sqrt(dx * dx + dy * dy)));
+                local_max = std::max(local_max, radius[i]);
             }
-            const float dx = x - beam_x, dy = y - beam_y;
-            radius[i] = static_cast<int>(std::lround(std::sqrt(dx * dx + dy * dy)));
-            max_radius = std::max(max_radius, radius[i]);
-        }
+        // max is associative, so folding the per-worker maxima gives the same answer whatever
+        // order they finish in.
+        int prev = max_radius_atomic.load();
+        while (prev < local_max && !max_radius_atomic.compare_exchange_weak(prev, local_max)) {}
+    });
+    const int max_radius = max_radius_atomic.load();
+
 
     // Pool the background over a small box before testing it. A background of a fraction of
     // a count per pixel per frame gives no single pixel enough counts to tell a shadow from
     // a Poisson hole; the stop and its arm are wider than the box, so pooling costs no
     // resolution that matters and multiplies the statistics by the pixels in the box.
-    std::vector<double> num(n_pixels), den(n_pixels);
-    for (int i = 0; i < n_pixels; i++) {
-        num[i] = valid[i] ? mean[i] : 0.0;
-        den[i] = valid[i] ? 1.0 : 0.0;
-    }
-    const auto pooled_sum = box_sum(num, W, H, POOL_PX);
-    const auto pooled_count = box_sum(den, W, H, POOL_PX);
+    // The count is a count: at most 25 pixels, so an integer box sum is exact and costs half the
+    // memory of the floating-point one it replaces. The background itself stays in double - its
+    // running sum adds and subtracts across a whole row, and in float the rounding of the two would
+    // not cancel, which moves pixels across the shadow threshold below.
+    std::vector<double>  num(n_pixels);
+    std::vector<int32_t> den(n_pixels);
+    ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
+        for (int i = lo; i < hi; i++) {
+            num[i] = valid[i] ? mean[i] : 0.0;
+            den[i] = valid[i] ? 1 : 0;
+        }
+    });
+    const auto pooled_sum = box_sum(num, W, H, POOL_PX, nthreads);
+    const auto pooled_count = box_sum(den, W, H, POOL_PX, nthreads);
     std::vector<float> pooled(n_pixels, 0.0f);
-    for (int i = 0; i < n_pixels; i++)
-        if (pooled_count[i] > 0)
-            pooled[i] = static_cast<float>(pooled_sum[i] / pooled_count[i]);
+    ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
+        for (int i = lo; i < hi; i++)
+            if (pooled_count[i] > 0)
+                pooled[i] = static_cast<float>(pooled_sum[i] / pooled_count[i]);
+    });
+
 
     // Azimuthal comparison: the median of the ring, iterated so the shadow stays out of the
     // baseline it is measured against.
+    //
+    // The iteration only ever excludes pixels whose pooled background is below a cut, and dividing by
+    // a positive baseline is monotone - so the excluded pixels of a ring are exactly the lowest ones,
+    // and the pixels the next median is taken over are exactly the rest. Each iteration's median is
+    // therefore an order statistic of the ring's values, which do not change: bin and sort the rings
+    // once, then each iteration picks a rank and counts a prefix. That replaces nine full-image
+    // passes (use / ratio / excluded, three times) with one, and three re-binnings with none.
+    const RingValues rings = bin_by_ring(pooled, valid, radius, max_radius, nthreads);
+    std::vector<float> baseline(max_radius + 1, 0.0f);
+    {
+        std::vector<int> excluded_in_ring(max_radius + 1, 0);
+        for (int iter = 0; iter < 3; iter++)
+            ParallelFor(max_radius + 1, nthreads, [&](int r) {
+                const int lo = rings.offset[r], hi = rings.offset[r + 1];
+                const int n = hi - lo;
+                const int m = excluded_in_ring[r];
+                const int avail = n - m;
+                baseline[r] = (avail <= 0) ? 0.0f : rings.values[lo + m + avail / 2];
+                // Counted the same way the per-pixel test below is written, so the two agree bit for
+                // bit; the ring is sorted, so this is the length of a prefix.
+                const float d = std::max(baseline[r], 1e-6f);
+                int excl = 0;
+                while (excl < n && rings.values[lo + excl] / d < SHADOW_RATIO)
+                    excl++;
+                excluded_in_ring[r] = excl;
+            });
+    }
     std::vector<float> ratio(n_pixels, 1.0f);
-    std::vector<char>  excluded(n_pixels, 0);
-    std::vector<float> baseline;
-    for (int iter = 0; iter < 3; iter++) {
-        std::vector<char> use(n_pixels);
-        for (int i = 0; i < n_pixels; i++)
-            use[i] = valid[i] && !excluded[i];
-        baseline = ring_median(pooled, use, radius, max_radius);
-        for (int i = 0; i < n_pixels; i++)
+    ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
+        for (int i = lo; i < hi; i++)
             if (valid[i])
                 ratio[i] = pooled[i] / std::max(baseline[radius[i]], 1e-6f);
-        for (int i = 0; i < n_pixels; i++)
-            excluded[i] = valid[i] && ratio[i] < SHADOW_RATIO;
-    }
+    });
+
 
     // A ring whose background was never counted carries no information to test a pixel against.
     // Walking outward, every ring before the first countable one lies wholly inside the stop - a
@@ -394,19 +590,22 @@ std::vector<uint32_t> ShadowFinder::GetMask() const {
     }
 
     std::vector<char> low(n_pixels, 0);
-    for (int i = 0; i < n_pixels; i++) {
-        if (!valid[i])
-            continue;
-        if (radius[i] <= blocked_out_to) {
-            low[i] = 1;
-            continue;
+    ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
+        for (int i = lo; i < hi; i++) {
+            if (!valid[i])
+                continue;
+            if (radius[i] <= blocked_out_to) {
+                low[i] = 1;
+                continue;
+            }
+            const double counted = static_cast<double>(frames) * pooled_count[i];
+            low[i] = ratio[i] < SHADOW_RATIO && baseline[radius[i]] * counted >= MIN_EXPECTED_COUNTS;
         }
-        const double counted = frames * pooled_count[i];
-        low[i] = ratio[i] < SHADOW_RATIO && baseline[radius[i]] * counted >= MIN_EXPECTED_COUNTS;
-    }
+    });
+
 
     // The shadow is the low region connected to the beam centre, bridging the gaps it crosses.
-    const std::vector<char> bridged = dilate(low, W, H, BRIDGE_PX);
+    const std::vector<char> bridged = dilate(low, W, H, BRIDGE_PX, nthreads);
     std::vector<int> seeds;
     for (int i = 0; i < n_pixels; i++)
         if (radius[i] < 4)
@@ -416,37 +615,46 @@ std::vector<uint32_t> ShadowFinder::GetMask() const {
     for (int i = 0; i < n_pixels; i++)
         region[i] = low[i] && connected[i];
 
+
     // Recorded reflections. A small cluster is required so a single-frame zinger does not count.
     std::vector<char> lit(n_pixels, 0);
-    for (int i = 0; i < n_pixels; i++)
-        lit[i] = (valid_count[i] > 0) && (max_value[i] >= MIN_REFLECTION);
+    ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
+        for (int i = lo; i < hi; i++)
+            lit[i] = (valid_count[i] > 0) && (max_value[i] >= MIN_REFLECTION);
+    });
     std::vector<char> reflection(n_pixels, 0);
-    for (int y = 0; y < H; y++)
-        for (int x = 0; x < W; x++) {
-            const int i = y * W + x;
-            if (!lit[i]) continue;
-            int neighbours = 0;
-            for (int dy = -1; dy <= 1; dy++)
-                for (int dx = -1; dx <= 1; dx++) {
-                    const int yy = y + dy, xx = x + dx;
-                    if ((dx || dy) && yy >= 0 && yy < H && xx >= 0 && xx < W && lit[yy * W + xx])
-                        neighbours++;
-                }
-            reflection[i] = (neighbours >= 2);
-        }
+    ParallelChunks(H, nthreads, [&](int ylo, int yhi) {
+        for (int y = ylo; y < yhi; y++)
+            for (int x = 0; x < W; x++) {
+                const int i = y * W + x;
+                if (!lit[i]) continue;
+                int neighbours = 0;
+                for (int dy = -1; dy <= 1; dy++)
+                    for (int dx = -1; dx <= 1; dx++) {
+                        const int yy = y + dy, xx = x + dx;
+                        if ((dx || dy) && yy >= 0 && yy < H && xx >= 0 && xx < W && lit[yy * W + xx])
+                            neighbours++;
+                    }
+                reflection[i] = (neighbours >= 2);
+            }
+    });
+
 
     // Grow the soft boundary, round it and fill the disk interior.
-    const std::vector<char> penumbra = dilate(region, W, H, PENUMBRA_MAX_PX);
+    const std::vector<char> penumbra = dilate(region, W, H, PENUMBRA_MAX_PX, nthreads);
     for (int i = 0; i < n_pixels; i++)
         if (penumbra[i] && valid[i] && ratio[i] < PENUMBRA_RATIO)
             region[i] = 1;
 
-    region = erode(dilate(region, W, H, 2), W, H, 2);
+
+    region = erode(dilate(region, W, H, 2, nthreads), W, H, 2, nthreads);
+
     region = fill_holes(region, W, H);
 
+
     // Expose recorded reflections - done last, with no fill afterwards, so a spot the shadow
     // still covered is given back rather than re-enclosed.
-    const std::vector<char> reflection_grown = dilate(reflection, W, H, 1);
+    const std::vector<char> reflection_grown = dilate(reflection, W, H, 1, nthreads);
     for (int i = 0; i < n_pixels; i++)
         if (reflection_grown[i])
             region[i] = 0;
diff --git a/image_analysis/beam_stop/ShadowFinder.h b/image_analysis/beam_stop/ShadowFinder.h
index 66740237..4899efb3 100644
--- a/image_analysis/beam_stop/ShadowFinder.h
+++ b/image_analysis/beam_stop/ShadowFinder.h
@@ -74,7 +74,9 @@ public:
 
     // Compute the shadow mask (1 = shadow, 0 = keep), of the converted pixel count.
     // Recomputed from the accumulators on each call - meant to be called once at the end.
-    [[nodiscard]] std::vector<uint32_t> GetMask() const;
+    // nthreads = 0 asks for all hardware threads. The per-pixel passes over a 16M-pixel detector
+    // dominate this, and they are all exactly parallel.
+    [[nodiscard]] std::vector<uint32_t> GetMask(size_t nthreads = 0) const;
 
     // Mean counts per pixel over the frames added, NAN where nothing was counted. This is the
     // projection GetMask() tests, so anything else that wants the background before indexing
diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp
index 5771053c..38969704 100644
--- a/rugnux/Rugnux.cpp
+++ b/rugnux/Rugnux.cpp
@@ -412,7 +412,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
     }
 
     if (want_shadow) {
-        const auto shadow = finder.GetMask();
+        const auto shadow = finder.GetMask(config_.nthreads);
         const auto shadowed = std::count(shadow.begin(), shadow.end(), 1u);
         pixel_mask_.LoadBeamStopMask(experiment_, shadow);
         logger.Info("Beam stop shadow: {} pixels ({:.2f}% of the detector) found in {} images",