diff --git a/common/Reflection.h b/common/Reflection.h
index 4f7b4251..bf78ee49 100644
--- a/common/Reflection.h
+++ b/common/Reflection.h
@@ -43,6 +43,8 @@ struct MergedReflection {
     float sigma_half[2] = {NAN, NAN};
     float d = 0.0;
     bool rfree_flag = false;
+    float F = NAN;       // French-Wilson amplitude |F| (filled by ApplyFrenchWilson at end of merge)
+    float sigmaF = NAN;  // its sigma
 };
 
 
diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp
index c2e10d7d..b01bc0a7 100644
--- a/image_analysis/WriteReflections.cpp
+++ b/image_analysis/WriteReflections.cpp
@@ -181,6 +181,8 @@ void WriteMmcifReflections(const std::vector<MergedReflection> &reflections,
     out << "_refln.index_l\n";
     out << "_refln.intensity_meas\n";
     out << "_refln.intensity_sigma\n";
+    out << "_refln.F_meas_au\n";
+    out << "_refln.F_meas_sigma_au\n";
     out << "_refln.status_free\n";
     out << "_refln.status\n";
 
@@ -190,6 +192,8 @@ void WriteMmcifReflections(const std::vector<MergedReflection> &reflections,
             << std::setw(5) << r.l << " "
             << std::setw(14) << Fmt(r.I, 4) << " "
             << std::setw(14) << Fmt(r.sigma, 4) << " "
+            << std::setw(14) << Fmt(r.F, 4) << " "
+            << std::setw(14) << Fmt(r.sigmaF, 4) << " "
             << (r.rfree_flag ? 1 : 0) << " "
             << "o"  // 'o' = observed
             << "\n";
@@ -224,13 +228,15 @@ void WriteMtzReflections(const std::vector<MergedReflection> &reflections,
     mtz.add_column("L", 'H', dataset_id, -1, false);
     mtz.add_column("IMEAN", 'J', dataset_id, -1, false);
     mtz.add_column("SIGIMEAN", 'Q', dataset_id, -1, false);
+    mtz.add_column("F", 'F', dataset_id, -1, false);      // French-Wilson amplitude
+    mtz.add_column("SIGF", 'Q', dataset_id, -1, false);
     mtz.add_column("FreeR_flag", 'I', dataset_id, -1, false);
 
     // Number of rows
     mtz.nreflections = reflections.size();
 
     // Allocate data table
-    mtz.data.reserve(reflections.size() * 6);
+    mtz.data.reserve(reflections.size() * 8);
 
     for (const auto& r : reflections) {
         mtz.data.push_back(static_cast<float>(r.h));
@@ -238,6 +244,8 @@ void WriteMtzReflections(const std::vector<MergedReflection> &reflections,
         mtz.data.push_back(static_cast<float>(r.l));
         mtz.data.push_back(r.I);
         mtz.data.push_back(r.sigma);
+        mtz.data.push_back(r.F);
+        mtz.data.push_back(r.sigmaF);
         mtz.data.push_back(r.rfree_flag ? 1 : 0);
     }
 
@@ -254,7 +262,8 @@ void WriteHKLReflections(const std::vector<MergedReflection> &reflections,
     for (const auto &r: reflections)
         hkl_file << r.h << " " << r.k << " " << r.l << " "
                 << r.I << " " << r.sigma << " "
-                << (r.rfree_flag ? 1 : 0) << std::endl;
+                << (r.rfree_flag ? 1 : 0) << " "
+                << r.F << " " << r.sigmaF << std::endl;
 
     hkl_file.close();
 }
diff --git a/image_analysis/scale_merge/CMakeLists.txt b/image_analysis/scale_merge/CMakeLists.txt
index ea3e655a..9eb12b27 100644
--- a/image_analysis/scale_merge/CMakeLists.txt
+++ b/image_analysis/scale_merge/CMakeLists.txt
@@ -15,6 +15,8 @@ ADD_LIBRARY(JFJochScaleMerge
         HKLKey.h
         RfreeFlags.cpp
         RfreeFlags.h
+        FrenchWilson.cpp
+        FrenchWilson.h
         ScalingResult.h
         ScalingResult.cpp)
 TARGET_LINK_LIBRARIES(JFJochScaleMerge Ceres::ceres Eigen3::Eigen JFJochCommon)
diff --git a/image_analysis/scale_merge/FrenchWilson.cpp b/image_analysis/scale_merge/FrenchWilson.cpp
new file mode 100644
index 00000000..0f3a5d8c
--- /dev/null
+++ b/image_analysis/scale_merge/FrenchWilson.cpp
@@ -0,0 +1,131 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "FrenchWilson.h"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include <vector>
+
+#include "../../common/ResolutionShells.h"
+#include "gemmi/symmetry.hpp"
+
+namespace {
+
+struct Posterior {
+    double mean_I;   // <J>   (posterior mean true intensity)
+    double mean_F;   // <|F|> (posterior mean amplitude)
+};
+
+// Posterior moments of the true intensity J >= 0 given a measurement I +/- sigma and the Wilson
+// prior with mean sigma_wilson. Integrated numerically over J in [0, I + 8 sigma] with a log-shift
+// so the exponentials never overflow/underflow. acentric: p(J) ~ exp(-J/S); centric:
+// p(J) ~ exp(-J/2S)/sqrt(J).
+Posterior integrate_posterior(double I, double sigma, double sigma_wilson, bool centric, int npts) {
+    const double inv_2s2 = 1.0 / (2.0 * sigma * sigma);
+    const double j_max = std::max(I, 0.0) + 8.0 * sigma;
+    const double dj = j_max / npts;
+
+    std::vector<double> logw(npts);
+    double max_logw = -std::numeric_limits<double>::infinity();
+    for (int i = 0; i < npts; ++i) {
+        const double j = (i + 0.5) * dj;
+        const double diff = I - j;
+        const double log_prior = centric ? (-j / (2.0 * sigma_wilson) - 0.5 * std::log(j))
+                                          : (-j / sigma_wilson);
+        logw[i] = log_prior - diff * diff * inv_2s2;
+        max_logw = std::max(max_logw, logw[i]);
+    }
+
+    double sum_w = 0, sum_wI = 0, sum_wF = 0;
+    for (int i = 0; i < npts; ++i) {
+        const double j = (i + 0.5) * dj;
+        const double w = std::exp(logw[i] - max_logw);
+        if (!std::isfinite(w))
+            continue;
+        sum_w += w;
+        sum_wI += w * j;
+        sum_wF += w * std::sqrt(j);
+    }
+    if (sum_w <= 0.0) {
+        const double j = std::max(I, 0.0);
+        return {j, std::sqrt(j)};
+    }
+    return {sum_wI / sum_w, sum_wF / sum_w};
+}
+
+}  // namespace
+
+void ApplyFrenchWilson(std::vector<MergedReflection> &merged, int32_t space_group_number,
+                       const FrenchWilsonOptions &opts) {
+    auto naive = [](MergedReflection &r) {
+        const double ip = std::max(r.I, 0.0f);
+        r.F = static_cast<float>(std::sqrt(ip));
+        r.sigmaF = (ip > 0.0 && std::isfinite(r.sigma)) ? static_cast<float>(r.sigma / (2.0 * std::sqrt(ip)))
+                                                        : NAN;
+    };
+
+    const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number);
+    if (sg == nullptr || merged.empty()) {
+        for (auto &r : merged) naive(r);
+        return;
+    }
+    const gemmi::GroupOps gops = sg->operations();
+
+    float d_min = std::numeric_limits<float>::max(), d_max = 0.0f;
+    for (const auto &r : merged)
+        if (std::isfinite(r.d) && r.d > 0.0f) {
+            d_min = std::min(d_min, r.d);
+            d_max = std::max(d_max, r.d);
+        }
+    if (!(d_min < d_max && d_min > 0.0f)) {
+        for (auto &r : merged) naive(r);
+        return;
+    }
+
+    // Wilson mean intensity <I/epsilon> per resolution shell.
+    ResolutionShells shells(d_min * 0.999f, d_max * 1.001f, opts.num_shells);
+    std::vector<double> shell_sum(opts.num_shells, 0.0);
+    std::vector<int> shell_count(opts.num_shells, 0);
+    double global_sum = 0.0;
+    int global_count = 0;
+    auto epsilon = [&](const MergedReflection &r) {
+        return std::max(1, gops.epsilon_factor_without_centering({{r.h, r.k, r.l}}));
+    };
+    for (const auto &r : merged) {
+        if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f)
+            continue;
+        const double i_over_eps = r.I / epsilon(r);
+        global_sum += i_over_eps;
+        ++global_count;
+        if (const auto s = shells.GetShell(r.d)) {
+            shell_sum[*s] += i_over_eps;
+            ++shell_count[*s];
+        }
+    }
+    const double global_mean = global_count > 0 ? std::max(global_sum / global_count, 1e-10) : 1.0;
+    std::vector<double> shell_mean(opts.num_shells, global_mean);
+    for (int s = 0; s < opts.num_shells; ++s)
+        if (shell_count[s] >= opts.min_reflections_per_shell)
+            shell_mean[s] = std::max(shell_sum[s] / shell_count[s], 1e-10);
+
+    for (auto &r : merged) {
+        if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f) {
+            naive(r);
+            continue;
+        }
+        // Strong reflections: the FW correction is negligible, <|F|> = sqrt(I).
+        if (r.I > opts.strong_cutoff * r.sigma) {
+            naive(r);
+            continue;
+        }
+        const auto s = shells.GetShell(r.d);
+        const double sigma_wilson = epsilon(r) * (s ? shell_mean[*s] : global_mean);
+        const bool centric = gops.is_reflection_centric({{r.h, r.k, r.l}});
+        const Posterior post = integrate_posterior(r.I, r.sigma, sigma_wilson, !centric,
+                                                   opts.integration_points);
+        r.F = static_cast<float>(post.mean_F);
+        r.sigmaF = static_cast<float>(std::sqrt(std::max(0.0, post.mean_I - post.mean_F * post.mean_F)));
+    }
+}
diff --git a/image_analysis/scale_merge/FrenchWilson.h b/image_analysis/scale_merge/FrenchWilson.h
new file mode 100644
index 00000000..9ea38350
--- /dev/null
+++ b/image_analysis/scale_merge/FrenchWilson.h
@@ -0,0 +1,24 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <cstdint>
+#include <vector>
+
+#include "../../common/Reflection.h"
+
+struct FrenchWilsonOptions {
+    int num_shells = 20;                 // resolution shells for the Wilson mean <I>
+    int min_reflections_per_shell = 20;  // below this, fall back to the overall mean for the shell
+    int integration_points = 400;        // posterior integration steps (weak reflections only)
+    double strong_cutoff = 4.0;          // I/sigma above which <|F|> = sqrt(I) (FW bias negligible)
+};
+
+// Fill F and sigmaF on every merged reflection with the French-Wilson estimate of the amplitude:
+// the posterior mean |F| given the measured intensity I and its sigma under the Wilson prior. The
+// prior uses the resolution-shell mean intensity, the correct centric/acentric form, and the
+// reflection's epsilon (symmetry-enhancement) multiplicity. Strong reflections reduce to sqrt(I);
+// reflections with an unusable I/sigma fall back to sqrt(max(I,0)) with propagated sigma.
+void ApplyFrenchWilson(std::vector<MergedReflection> &merged, int32_t space_group_number,
+                       const FrenchWilsonOptions &opts = {});
diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp
index 9bfa9b62..405a9b26 100644
--- a/image_analysis/scale_merge/Merge.cpp
+++ b/image_analysis/scale_merge/Merge.cpp
@@ -18,6 +18,7 @@
 #include "../../common/Definitions.h"
 #include "HKLKey.h"
 #include "RfreeFlags.h"
+#include "FrenchWilson.h"
 
 namespace {
     // Deterministic CC1/2 half-set assignment: a splitmix64 bit-mix of the image's stable index.
@@ -472,6 +473,7 @@ std::vector<MergedReflection> MergeOnTheFly::ExportReflections() {
     }
 
     AssignRfreeFlags(out, space_group_number, scaling_settings.GetRfreeFraction());
+    ApplyFrenchWilson(out, space_group_number);
     return out;
 }
 
diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp
index 71ec34d9..af6db35d 100644
--- a/image_analysis/scale_merge/RotationScaleMerge.cpp
+++ b/image_analysis/scale_merge/RotationScaleMerge.cpp
@@ -17,6 +17,7 @@
 
 #include "HKLKey.h"
 #include "RfreeFlags.h"
+#include "FrenchWilson.h"
 #include "ResolutionCutoff.h"
 #include "../../common/CorrelationCoefficient.h"
 #include "../../common/CrystalLattice.h"
@@ -1300,6 +1301,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
         result.merged, d_min_limit, resolution_cutoff_method, resolution_cc_target, for_search, logger);
 
     AssignRfreeFlags(result.merged, x.GetSpaceGroupNumber().value_or(1), rfree_fraction);
+    ApplyFrenchWilson(result.merged, x.GetSpaceGroupNumber().value_or(1));
 
     if (reject_count > 0)
         logger.Info("Merge outlier rejection: dropped {} observations", reject_count);