diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp
index cfb829e9..08f875ba 100644
--- a/image_analysis/scale_merge/Merge.cpp
+++ b/image_analysis/scale_merge/Merge.cpp
@@ -65,6 +65,19 @@ void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, bool cc_mask) {
         auto hkl_key = hkl.pack();
         sigma_corr = CorrectedSigma(I_corr, sigma_corr, hkl_key);
 
+        // Robust outlier rejection: drop this observation if it sits more than
+        // reject_nsigma error-model sigmas from the reflection's median. Needs the active
+        // error model so sigma_corr reflects the real scatter (else the threshold is the
+        // bare counting sigma and would cull good partials).
+        if (reject_outliers && error_model_active) {
+            const auto mit = reject_median_I.find(hkl_key);
+            if (mit != reject_median_I.end() &&
+                std::fabs(I_corr - mit->second) > reject_nsigma * sigma_corr) {
+                ++reject_count;
+                continue;
+            }
+        }
+
         auto it = accumulator.find(hkl_key);
         if (it == accumulator.end())
             it = accumulator.emplace(hkl_key, MergeAccum{
@@ -110,6 +123,8 @@ void MergeOnTheFly::RefineErrorModel(const std::vector<IntegrationOutcome> &outc
     error_model_a = 1.0;
     error_model_b = 0.0;
     error_model_mean_I.clear();
+    reject_median_I.clear();
+    reject_count = 0;
 
     // --- 1. Collect accepted, scaled observations grouped by symmetry-equivalent hkl,
     //        applying exactly the filters AddImage uses. ---
@@ -156,6 +171,16 @@ void MergeOnTheFly::RefineErrorModel(const std::vector<IntegrationOutcome> &outc
             continue;
         const double mean = sum_wI / sum_w;
         error_model_mean_I[key] = static_cast<float>(mean);
+        // Robust centre for outlier rejection: the median intensity (resists the very
+        // outliers the inverse-variance mean is being protected from).
+        {
+            std::vector<float> iv;
+            iv.reserve(obs.size());
+            for (const auto &o: obs)
+                iv.push_back(o.I);
+            std::nth_element(iv.begin(), iv.begin() + iv.size() / 2, iv.end());
+            reject_median_I[key] = iv[iv.size() / 2];
+        }
         const double I2 = mean * mean;
         for (const auto &o: obs) {
             const double w = 1.0 / (static_cast<double>(o.sigma) * o.sigma);
diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h
index d69ee817..a1ea2776 100644
--- a/image_analysis/scale_merge/Merge.h
+++ b/image_analysis/scale_merge/Merge.h
@@ -87,6 +87,16 @@ class MergeOnTheFly {
     std::unordered_map<uint64_t, float> error_model_mean_I;
     [[nodiscard]] float CorrectedSigma(float I_corr, float sigma_corr, uint64_t hkl_key) const;
 
+    // Optional per-observation outlier rejection: drop observations whose corrected
+    // intensity lies more than reject_nsigma error-model sigmas from the reflection's
+    // *median* (a robust centre). The error-model sigma already captures the genuine
+    // (e.g. partiality) scatter, so this removes only the tail beyond it - zingers,
+    // overlaps, mis-indexed frames - not good partials. Populated by RefineErrorModel.
+    bool   reject_outliers = false;
+    double reject_nsigma = 6.0;
+    std::unordered_map<uint64_t, float> reject_median_I;
+    size_t reject_count = 0;
+
     bool Mask(const IntegrationOutcome &outcome, bool cc_mask);
 public:
     MergeOnTheFly(const DiffractionExperiment &x);
@@ -99,6 +109,10 @@ public:
     [[nodiscard]] double ErrorModelA() const { return error_model_a; }
     [[nodiscard]] double ErrorModelB() const { return error_model_b; }
 
+    // Enable per-observation outlier rejection (call after RefineErrorModel, before AddImage).
+    void SetRejectOutliers(double nsigma) { reject_outliers = true; reject_nsigma = nsigma; }
+    [[nodiscard]] size_t RejectedCount() const { return reject_count; }
+
     void AddImage(const IntegrationOutcome& outcome, bool cc_mask = false);
 
     MergeStatistics MergeStats(const std::vector<MergedReflection> &merged,
diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp
index 01ed1bf6..ef847d86 100644
--- a/tools/jfjoch_process.cpp
+++ b/tools/jfjoch_process.cpp
@@ -1034,8 +1034,15 @@ int main(int argc, char **argv) {
             const double isa = (b > 0.0) ? 1.0 / b : std::numeric_limits<double>::infinity();
             logger.Info("Error model: sigma'^2 = {:.3f} sigma^2 + ({:.4f} I)^2   ISa = {:.1f}", a, b, isa);
         }
+        if (const char *rj = std::getenv("PR_REJECT")) { // TEMP A/B: per-observation outlier rejection
+            const double nsig = std::atof(rj);
+            merge_engine.SetRejectOutliers(nsig);
+            logger.Info("Outlier rejection enabled: dropping observations > {:.1f} sigma from the per-reflection median", nsig);
+        }
         for (auto &i : indexer.GetIntegrationOutcome())
             merge_engine.AddImage(i);
+        if (merge_engine.RejectedCount() > 0)
+            logger.Info("Outlier rejection removed {} observations", merge_engine.RejectedCount());
 
         auto merged_reflections = merge_engine.ExportReflections();
         auto merged_statistics = merge_engine.MergeStats(merged_reflections, indexer.GetIntegrationOutcome(), reference_data);