From 2ba28aea0e8cbda7229ee8378249705a4f47cf43 Mon Sep 17 00:00:00 2001
From: leonarski_f <filip.leonarski@psi.ch>
Date: Tue, 16 Jun 2026 14:19:28 +0200
Subject: [PATCH] LoadFCalcFromMtz: fall back to an intensity column when
 F-model is absent

PixelRefine's reference can now be a merged/observed-intensity MTZ, not only a
calculated F-model. If F-model is missing, read a mean-intensity column (IMEAN,
which a jfjoch merge writes; also I/IOBS/Iobs/I-obs, or any 'J'-type column) and
use it directly as I_ref (F-model is squared, intensities are not).

This enables a SELF-SEEDED / EM-style workflow: run a first pass (traditional
integration, or PixelRefine against an external reference) to produce a merge,
then re-run PixelRefine with `-z <that merge>.mtz` so it scales against the data's
own intensities - free of the reference structure's non-isomorphism bias - and
iterate. Verified: a jfjoch IMEAN merge loads as the reference and PixelRefine
runs against it (crystal 2, 55279 reflections, CC1/2 94.3).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 image_analysis/LoadFCalcFromMtz.cpp | 30 +++++++++++++++++++++++------
 1 file changed, 24 insertions(+), 6 deletions(-)

diff --git a/image_analysis/LoadFCalcFromMtz.cpp b/image_analysis/LoadFCalcFromMtz.cpp
index b7153f32..503ed9fd 100644
--- a/image_analysis/LoadFCalcFromMtz.cpp
+++ b/image_analysis/LoadFCalcFromMtz.cpp
@@ -16,9 +16,27 @@ std::vector<MergedReflection> LoadFCalcFromMtz(const std::string& path) {
     gemmi::Mtz mtz;
     mtz.read_file_gz(path, true);
 
-    const gemmi::Mtz::Column* fc = mtz.column_with_label("F-model", nullptr, 'F');
-    if (fc == nullptr)
-        throw std::runtime_error("MTZ does not contain F-model column");
+    // Prefer a calculated structure factor F-model (e.g. a reference structure): I_ref = F^2.
+    const gemmi::Mtz::Column* col = mtz.column_with_label("F-model", nullptr, 'F');
+    const bool square = (col != nullptr);
+
+    // Otherwise fall back to a merged/observed intensity column, used directly as I_ref. This
+    // lets PixelRefine be SELF-SEEDED from the data's own previous merge (a jfjoch merge writes
+    // IMEAN) instead of an external reference structure - the EM-style outer loop, free of the
+    // reference's non-isomorphism bias.
+    if (col == nullptr) {
+        for (const char* label : {"IMEAN", "I", "IOBS", "Iobs", "I-obs"}) {
+            col = mtz.column_with_label(label, nullptr, 'J');
+            if (col != nullptr)
+                break;
+        }
+    }
+    if (col == nullptr) {
+        for (const auto& c : mtz.columns)
+            if (c.type == 'J') { col = &c; break; } // any mean-intensity column
+    }
+    if (col == nullptr)
+        throw std::runtime_error("MTZ has no F-model or intensity (J) column to use as reference");
 
     std::vector<MergedReflection> result;
     result.reserve(static_cast<std::size_t>(mtz.nreflections));
@@ -27,16 +45,16 @@ std::vector<MergedReflection> LoadFCalcFromMtz(const std::string& path) {
 
     for (int i = 0; i < mtz.nreflections; ++i) {
         const std::size_t row = static_cast<std::size_t>(i) * stride;
-        const float f = (*fc)[static_cast<std::size_t>(i)];
+        const float v = (*col)[static_cast<std::size_t>(i)];
 
-        if (std::isnan(f))
+        if (std::isnan(v))
             continue;
 
         MergedReflection r;
         r.h = static_cast<int32_t>(mtz.data[row + 0]);
         r.k = static_cast<int32_t>(mtz.data[row + 1]);
         r.l = static_cast<int32_t>(mtz.data[row + 2]);
-        r.I = f * f;
+        r.I = square ? v * v : v;
         r.sigma = NAN;
         r.d = 0.0f;