From 5baf2a2cacf7005cdab0a753e69565649f422f04 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 13 Jul 2026 14:53:23 +0200 Subject: [PATCH 01/90] docker: static libdbus + disable Qt glib to shed the viewer's .so tail The jfjoch_viewer runtime pulled a large dynamic .so tail that traces to two Qt-side integrations, neither of which the viewer actually needs: - Qt6::DBus -> system libdbus-1.so -> libsystemd (libzstd/liblz4/libcap/ libgcrypt/libgpg-error) + libselinux (libpcre2). The viewer is a pure session-bus client (registers ch.psi.jfjoch_viewer + exports an adaptor for single-instance/remote-control), so the daemon-side systemd/selinux features are irrelevant. Build a static libdbus from source with those integrations disabled and link it into the viewer, the same from-source-static pattern already used for OpenSSL. Swap the dbus dev package for expat (dbus's configure-time dep; the client libdbus-1 links neither expat nor systemd) so nothing pulls the system shared dbus back in. Qt discovers dbus through find_package(DBus1) (its FindWrapDBus1.cmake), NOT pkg-config -- so point it at our build with -DDBus1_DIR=/lib/cmake/ DBus1 (the CMake package config dbus installs, which imports libdbus-1.a) and force QT_FEATURE_dbus_linked=ON so QtDBus links the archive instead of dlopen'ing a (now nonexistent) libdbus-1.so.3 at runtime. - Qt's glib event-dispatcher -> libglib-2.0/libgthread/libpcre. Nothing in the viewer (or any lib it links) drives a GLib main context, so switch Qt to its own QEventDispatcherUNIX via QT_FEATURE_glib=OFF and drop the glib dev package. Applied identically across rocky8, rocky9, ubuntu2204, ubuntu2404. Note: system libcrypto still appears at runtime, but it is dragged in by the dynamic system krb5/GSSAPI (built against OpenSSL on RHEL/Rocky), which we keep dynamic for the planned httpd-reverse-proxy Kerberos auth -- not from Qt or curl, which both statically link the from-source OpenSSL. Co-Authored-By: Claude Opus 4.8 (1M context) --- docker/rocky8/Dockerfile | 34 ++++++++++++++++++++++++++++++++-- docker/rocky9/Dockerfile | 34 ++++++++++++++++++++++++++++++++-- docker/ubuntu2204/Dockerfile | 34 ++++++++++++++++++++++++++++++++-- docker/ubuntu2404/Dockerfile | 34 ++++++++++++++++++++++++++++++++-- 4 files changed, 128 insertions(+), 8 deletions(-) diff --git a/docker/rocky8/Dockerfile b/docker/rocky8/Dockerfile index 5817d3e0..4d395230 100644 --- a/docker/rocky8/Dockerfile +++ b/docker/rocky8/Dockerfile @@ -3,6 +3,7 @@ LABEL authors="leonarski_f" ARG OPENSSL_VERSION=3.5.4 ARG QT_VERSION=6.9.1 +ARG DBUS_VERSION=1.14.10 ARG NODE_MAJOR=22 ARG EIGEN_VERSION=3.4.0 # HDF5, libtiff and libjpeg-turbo are built by the jungfraujoch CMake itself (FetchContent / @@ -47,11 +48,10 @@ RUN dnf -y install epel-release && \ libXinerama-devel \ mesa-libGL-devel \ mesa-libEGL-devel \ - dbus-devel \ + expat-devel \ krb5-devel \ zlib-devel \ zlib-static \ - glib2-devel \ fontconfig-devel \ zlib-static \ java-21-openjdk-headless \ @@ -96,6 +96,33 @@ ENV CXX=${GCC_TOOLSET_ROOT}/usr/bin/g++ ENV PKG_CONFIG_PATH=${OPENSSL_ROOT_DIR}/lib/pkgconfig:${OPENSSL_ROOT_DIR}/lib64/pkgconfig +# Static libdbus: link D-Bus into the viewer (Qt6::DBus) instead of pulling libdbus-1.so.3 at runtime. +# A distro libdbus-1.so drags in libsystemd (-> libzstd/liblz4/libcap/libgcrypt/libgpg-error) and +# libselinux (-> libpcre2). Disabling those integrations makes libdbus-1.a depend on ~libc alone, so +# that whole runtime .so tail disappears while the viewer keeps its single-instance / remote-control +# D-Bus feature -- it is a pure session-bus client (registers a name + exports an adaptor), so the +# daemon-side systemd/selinux features are irrelevant. Static-only install (no .so); the Qt build +# finds it through find_package(DBus1) via the -DDBus1_DIR hint (the CMake package config dbus +# installs, which imports libdbus-1.a). expat is only a configure-time dep of +# the dbus daemon; the client libdbus-1 links neither expat nor systemd. +RUN set -eux; \ + cd /tmp; \ + curl -LO https://dbus.freedesktop.org/releases/dbus/dbus-${DBUS_VERSION}.tar.xz; \ + tar -xf dbus-${DBUS_VERSION}.tar.xz; \ + cd dbus-${DBUS_VERSION}; \ + ./configure --prefix=/opt/dbus-${DBUS_VERSION}-static \ + --enable-static --disable-shared \ + --disable-systemd --without-systemdsystemunitdir \ + --disable-selinux --disable-apparmor --disable-libaudit \ + --disable-tests --disable-asserts \ + --disable-doxygen-docs --disable-xml-docs --disable-ducktype-docs; \ + make -j"$(nproc)"; \ + make install; \ + cd /; rm -rf /tmp/dbus-${DBUS_VERSION} /tmp/dbus-${DBUS_VERSION}.tar.xz + +# Put the static libdbus pkgconfig ahead of the system one so Qt6::DBus resolves to the .a +ENV PKG_CONFIG_PATH=/opt/dbus-${DBUS_VERSION}-static/lib/pkgconfig:${PKG_CONFIG_PATH} + ARG QT_PREFIX=/opt/qt-${QT_VERSION}-static RUN set -eux; \ cd /tmp; \ @@ -108,6 +135,8 @@ RUN set -eux; \ -DQT_BUILD_TESTS=OFF \ -DQT_BUILD_EXAMPLES=OFF \ -DQT_FEATURE_dbus=ON \ + -DQT_FEATURE_dbus_linked=ON \ + -DDBus1_DIR=/opt/dbus-${DBUS_VERSION}-static/lib/cmake/DBus1 \ -DQT_FEATURE_xcb=ON \ -DQT_FEATURE_xcb_xlib=OFF \ -DQT_FEATURE_xkbcommon_x11=ON \ @@ -115,6 +144,7 @@ RUN set -eux; \ -DQT_FEATURE_opengl_desktop=ON \ -DQT_FEATURE_opengl_dynamic=OFF \ -DQT_FEATURE_vulkan=OFF \ + -DQT_FEATURE_glib=OFF \ -DCMAKE_BUILD_TYPE=Release \ -DCMAKE_INSTALL_PREFIX=${QT_PREFIX} \ -DCMAKE_C_COMPILER=${CC} \ diff --git a/docker/rocky9/Dockerfile b/docker/rocky9/Dockerfile index 32c0e94d..d0e27deb 100644 --- a/docker/rocky9/Dockerfile +++ b/docker/rocky9/Dockerfile @@ -3,6 +3,7 @@ LABEL authors="leonarski_f" ARG OPENSSL_VERSION=3.5.4 ARG QT_VERSION=6.9.1 +ARG DBUS_VERSION=1.14.10 ARG NODE_MAJOR=22 ARG EIGEN_VERSION=3.4.0 # HDF5, libtiff and libjpeg-turbo are built by the jungfraujoch CMake itself (FetchContent / @@ -57,11 +58,10 @@ RUN dnf -y update && \ libXinerama-devel \ mesa-libGL-devel \ mesa-libEGL-devel \ - dbus-devel \ + expat-devel \ krb5-devel \ zlib-devel \ zlib-static \ - glib2-devel \ zlib-static \ fontconfig-devel \ mesa-dri-drivers \ @@ -99,6 +99,33 @@ ENV PATH=${GCC_TOOLSET_ROOT}/usr/bin:${PATH} ENV CC=${GCC_TOOLSET_ROOT}/usr/bin/gcc ENV CXX=${GCC_TOOLSET_ROOT}/usr/bin/g++ +# Static libdbus: link D-Bus into the viewer (Qt6::DBus) instead of pulling libdbus-1.so.3 at runtime. +# A distro libdbus-1.so drags in libsystemd (-> libzstd/liblz4/libcap/libgcrypt/libgpg-error) and +# libselinux (-> libpcre2). Disabling those integrations makes libdbus-1.a depend on ~libc alone, so +# that whole runtime .so tail disappears while the viewer keeps its single-instance / remote-control +# D-Bus feature -- it is a pure session-bus client (registers a name + exports an adaptor), so the +# daemon-side systemd/selinux features are irrelevant. Static-only install (no .so); the Qt build +# finds it through find_package(DBus1) via the -DDBus1_DIR hint (the CMake package config dbus +# installs, which imports libdbus-1.a). expat is only a configure-time dep of +# the dbus daemon; the client libdbus-1 links neither expat nor systemd. +RUN set -eux; \ + cd /tmp; \ + curl -LO https://dbus.freedesktop.org/releases/dbus/dbus-${DBUS_VERSION}.tar.xz; \ + tar -xf dbus-${DBUS_VERSION}.tar.xz; \ + cd dbus-${DBUS_VERSION}; \ + ./configure --prefix=/opt/dbus-${DBUS_VERSION}-static \ + --enable-static --disable-shared \ + --disable-systemd --without-systemdsystemunitdir \ + --disable-selinux --disable-apparmor --disable-libaudit \ + --disable-tests --disable-asserts \ + --disable-doxygen-docs --disable-xml-docs --disable-ducktype-docs; \ + make -j"$(nproc)"; \ + make install; \ + cd /; rm -rf /tmp/dbus-${DBUS_VERSION} /tmp/dbus-${DBUS_VERSION}.tar.xz + +# Put the static libdbus pkgconfig ahead of the system one so Qt6::DBus resolves to the .a +ENV PKG_CONFIG_PATH=/opt/dbus-${DBUS_VERSION}-static/lib/pkgconfig:${PKG_CONFIG_PATH} + # Build and install static Qt 6.9 with Core, Gui, Widgets, Charts, DBus ARG QT_PREFIX=/opt/qt-${QT_VERSION}-static RUN set -eux; \ @@ -112,6 +139,8 @@ RUN set -eux; \ -DQT_BUILD_TESTS=OFF \ -DQT_BUILD_EXAMPLES=OFF \ -DQT_FEATURE_dbus=ON \ + -DQT_FEATURE_dbus_linked=ON \ + -DDBus1_DIR=/opt/dbus-${DBUS_VERSION}-static/lib/cmake/DBus1 \ -DQT_FEATURE_xcb=ON \ -DQT_FEATURE_xcb_xlib=OFF \ -DQT_FEATURE_xkbcommon_x11=ON \ @@ -119,6 +148,7 @@ RUN set -eux; \ -DQT_FEATURE_opengl_desktop=ON \ -DQT_FEATURE_opengl_dynamic=OFF \ -DQT_FEATURE_vulkan=OFF \ + -DQT_FEATURE_glib=OFF \ -DCMAKE_BUILD_TYPE=Release \ -DCMAKE_INSTALL_PREFIX=${QT_PREFIX} \ -DCMAKE_C_COMPILER=${CC} \ diff --git a/docker/ubuntu2204/Dockerfile b/docker/ubuntu2204/Dockerfile index 7dbdb35a..b0c790b8 100644 --- a/docker/ubuntu2204/Dockerfile +++ b/docker/ubuntu2204/Dockerfile @@ -6,6 +6,7 @@ ENV DEBIAN_FRONTEND=noninteractive ARG OPENSSL_VERSION=3.5.4 ARG QT_VERSION=6.9.1 +ARG DBUS_VERSION=1.14.10 ARG NODE_MAJOR=22 ARG EIGEN_VERSION=3.4.0 # HDF5, libtiff and libjpeg-turbo are built by the jungfraujoch CMake itself (FetchContent / @@ -71,9 +72,8 @@ RUN set -eux; \ libgl1-mesa-dev \ libglx-dev \ libegl1-mesa-dev \ - libdbus-1-dev \ + libexpat1-dev \ zlib1g-dev \ - libglib2.0-dev \ libfontconfig1-dev \ libdrm-dev \ mesa-utils \ @@ -123,6 +123,33 @@ RUN set -eux; \ rm -rf /var/lib/apt/lists/*; \ node --version; npm --version; (corepack enable || true) +# Static libdbus: link D-Bus into the viewer (Qt6::DBus) instead of pulling libdbus-1.so.3 at runtime. +# A distro libdbus-1.so drags in libsystemd (-> libzstd/liblz4/libcap/libgcrypt/libgpg-error) and +# libselinux (-> libpcre2). Disabling those integrations makes libdbus-1.a depend on ~libc alone, so +# that whole runtime .so tail disappears while the viewer keeps its single-instance / remote-control +# D-Bus feature -- it is a pure session-bus client (registers a name + exports an adaptor), so the +# daemon-side systemd/selinux features are irrelevant. Static-only install (no .so); the Qt build +# finds it through find_package(DBus1) via the -DDBus1_DIR hint (the CMake package config dbus +# installs, which imports libdbus-1.a). expat is only a configure-time dep of +# the dbus daemon; the client libdbus-1 links neither expat nor systemd. +RUN set -eux; \ + cd /tmp; \ + curl -LO https://dbus.freedesktop.org/releases/dbus/dbus-${DBUS_VERSION}.tar.xz; \ + tar -xf dbus-${DBUS_VERSION}.tar.xz; \ + cd dbus-${DBUS_VERSION}; \ + ./configure --prefix=/opt/dbus-${DBUS_VERSION}-static \ + --enable-static --disable-shared \ + --disable-systemd --without-systemdsystemunitdir \ + --disable-selinux --disable-apparmor --disable-libaudit \ + --disable-tests --disable-asserts \ + --disable-doxygen-docs --disable-xml-docs --disable-ducktype-docs; \ + make -j"$(nproc)"; \ + make install; \ + cd /; rm -rf /tmp/dbus-${DBUS_VERSION} /tmp/dbus-${DBUS_VERSION}.tar.xz + +# Put the static libdbus pkgconfig ahead of the system one so Qt6::DBus resolves to the .a +ENV PKG_CONFIG_PATH=/opt/dbus-${DBUS_VERSION}-static/lib/pkgconfig:${PKG_CONFIG_PATH} + # Build and install static Qt with Core, Gui, Widgets, Charts, DBus ARG QT_PREFIX=/opt/qt-${QT_VERSION}-static RUN set -eux; \ @@ -136,6 +163,8 @@ RUN set -eux; \ -DQT_BUILD_TESTS=OFF \ -DQT_BUILD_EXAMPLES=OFF \ -DQT_FEATURE_dbus=ON \ + -DQT_FEATURE_dbus_linked=ON \ + -DDBus1_DIR=/opt/dbus-${DBUS_VERSION}-static/lib/cmake/DBus1 \ -DQT_FEATURE_xcb=ON \ -DQT_FEATURE_xcb_xlib=OFF \ -DQT_FEATURE_xkbcommon_x11=ON \ @@ -143,6 +172,7 @@ RUN set -eux; \ -DQT_FEATURE_opengl_desktop=ON \ -DQT_FEATURE_opengl_dynamic=OFF \ -DQT_FEATURE_vulkan=OFF \ + -DQT_FEATURE_glib=OFF \ -DCMAKE_BUILD_TYPE=Release \ -DCMAKE_INSTALL_PREFIX=${QT_PREFIX} \ -DCMAKE_C_COMPILER=${CC} \ diff --git a/docker/ubuntu2404/Dockerfile b/docker/ubuntu2404/Dockerfile index 796ee67b..bd7be0b1 100644 --- a/docker/ubuntu2404/Dockerfile +++ b/docker/ubuntu2404/Dockerfile @@ -6,6 +6,7 @@ ENV DEBIAN_FRONTEND=noninteractive ARG OPENSSL_VERSION=3.5.4 ARG QT_VERSION=6.9.1 +ARG DBUS_VERSION=1.14.10 ARG NODE_MAJOR=22 ARG EIGEN_VERSION=3.4.0 # HDF5, libtiff and libjpeg-turbo are built by the jungfraujoch CMake itself (FetchContent / @@ -71,9 +72,8 @@ RUN set -eux; \ libgl1-mesa-dev \ libglx-dev \ libegl1-mesa-dev \ - libdbus-1-dev \ + libexpat1-dev \ zlib1g-dev \ - libglib2.0-dev \ libfontconfig1-dev \ libdrm-dev \ libglvnd-dev \ @@ -110,6 +110,33 @@ RUN set -eux; \ rm -rf /var/lib/apt/lists/*; \ node --version; npm --version; (corepack enable || true) +# Static libdbus: link D-Bus into the viewer (Qt6::DBus) instead of pulling libdbus-1.so.3 at runtime. +# A distro libdbus-1.so drags in libsystemd (-> libzstd/liblz4/libcap/libgcrypt/libgpg-error) and +# libselinux (-> libpcre2). Disabling those integrations makes libdbus-1.a depend on ~libc alone, so +# that whole runtime .so tail disappears while the viewer keeps its single-instance / remote-control +# D-Bus feature -- it is a pure session-bus client (registers a name + exports an adaptor), so the +# daemon-side systemd/selinux features are irrelevant. Static-only install (no .so); the Qt build +# finds it through find_package(DBus1) via the -DDBus1_DIR hint (the CMake package config dbus +# installs, which imports libdbus-1.a). expat is only a configure-time dep of +# the dbus daemon; the client libdbus-1 links neither expat nor systemd. +RUN set -eux; \ + cd /tmp; \ + curl -LO https://dbus.freedesktop.org/releases/dbus/dbus-${DBUS_VERSION}.tar.xz; \ + tar -xf dbus-${DBUS_VERSION}.tar.xz; \ + cd dbus-${DBUS_VERSION}; \ + ./configure --prefix=/opt/dbus-${DBUS_VERSION}-static \ + --enable-static --disable-shared \ + --disable-systemd --without-systemdsystemunitdir \ + --disable-selinux --disable-apparmor --disable-libaudit \ + --disable-tests --disable-asserts \ + --disable-doxygen-docs --disable-xml-docs --disable-ducktype-docs; \ + make -j"$(nproc)"; \ + make install; \ + cd /; rm -rf /tmp/dbus-${DBUS_VERSION} /tmp/dbus-${DBUS_VERSION}.tar.xz + +# Put the static libdbus pkgconfig ahead of the system one so Qt6::DBus resolves to the .a +ENV PKG_CONFIG_PATH=/opt/dbus-${DBUS_VERSION}-static/lib/pkgconfig:${PKG_CONFIG_PATH} + # Build and install static Qt with Core, Gui, Widgets, Charts, DBus ARG QT_PREFIX=/opt/qt-${QT_VERSION}-static RUN set -eux; \ @@ -123,6 +150,8 @@ RUN set -eux; \ -DQT_BUILD_TESTS=OFF \ -DQT_BUILD_EXAMPLES=OFF \ -DQT_FEATURE_dbus=ON \ + -DQT_FEATURE_dbus_linked=ON \ + -DDBus1_DIR=/opt/dbus-${DBUS_VERSION}-static/lib/cmake/DBus1 \ -DQT_FEATURE_xcb=ON \ -DQT_FEATURE_xcb_xlib=OFF \ -DQT_FEATURE_xkbcommon_x11=ON \ @@ -130,6 +159,7 @@ RUN set -eux; \ -DQT_FEATURE_opengl_desktop=ON \ -DQT_FEATURE_opengl_dynamic=OFF \ -DQT_FEATURE_vulkan=OFF \ + -DQT_FEATURE_glib=OFF \ -DCMAKE_BUILD_TYPE=Release \ -DCMAKE_INSTALL_PREFIX=${QT_PREFIX} \ -DCMAKE_C_COMPILER=${CC} \ -- 2.54.0 From 254bd82d52df6add80a83e30a94f732643e2a5b3 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 13 Jul 2026 15:43:18 +0200 Subject: [PATCH 02/90] docker: add parallel image-build + in-container viewer-verify scripts build_images.sh builds (and optionally pushes) the four build-environment images as gitea.psi.ch/leonarski_f/jfjoch_:. Runs a JOBS-capped parallel job pool (default 2; each build internally runs make -j$(nproc), so a low cap avoids CPU/RAM thrash while still overlapping the network-bound base pull / install / download phases), streams each build to docker/build-logs/-.log, prints an OK/FAIL summary, pushes only the ones that succeeded, and exits non-zero on any failure. Context is the tiny per-variant dir (no COPY in the Dockerfiles), so the repo is never sent to the daemon. build-logs/ is already covered by the root .gitignore build*/ rule. build_in_rocky9.sh builds the viewer (JFJOCH_VIEWER_ONLY) inside a chosen variant's image and ldd-checks that the dbus/systemd/glib/selinux tail is gone -- the quick verification for the static-libdbus + glib-off changes. The repo is mounted but the build tree lives in the container's /tmp, so nothing root-owned lands in the working copy. Co-Authored-By: Claude Opus 4.8 (1M context) --- docker/build_images.sh | 98 +++++++++++++++++++++++++++++++++++++++ docker/build_in_rocky9.sh | 48 +++++++++++++++++++ 2 files changed, 146 insertions(+) create mode 100755 docker/build_images.sh create mode 100755 docker/build_in_rocky9.sh diff --git a/docker/build_images.sh b/docker/build_images.sh new file mode 100755 index 00000000..a464ef50 --- /dev/null +++ b/docker/build_images.sh @@ -0,0 +1,98 @@ +#!/usr/bin/env bash +# Build (and optionally push) the jungfraujoch build-environment images, in parallel. +# +# Each docker//Dockerfile is a *self-contained build environment* (toolchain, static Qt, +# static libdbus, static OpenSSL, Eigen, DIALS, XDS, Node, ...). It does NOT copy the project in -- +# the source is mounted and built at run time (see build_in_rocky9.sh). So the build context is just +# the per-variant directory; the huge repo is never sent to the Docker daemon. +# +# Images are tagged gitea.psi.ch/leonarski_f/jfjoch_: (matches the CI runner images). +# +# Concurrency (JOBS): +# JOBS controls how many image builds run at once (default 2). Each build internally runs +# `make -j$(nproc)` for Qt/dbus, so J parallel builds ~= J*nproc compile threads and several GB of +# RAM each at link time. JOBS=2 is a safe default that overlaps the (network-bound) base-image +# pull / package-install / source-download phases across builds without thrashing the CPU. Bump to +# 4 for full parallelism on a big machine (enough RAM + disk for 4 CUDA bases), or JOBS=1 for serial. +# +# Usage: +# docker/build_images.sh # all four, JOBS=2 +# docker/build_images.sh rocky9 rocky8 # subset +# JOBS=4 docker/build_images.sh # all four at once +# JOBS=1 docker/build_images.sh # serial +# TAG=2607b docker/build_images.sh # override the tag (default 2607b) +# PUSH=1 docker/build_images.sh # build (parallel) then push (serial; needs docker login) +# +# Parallel stdout would be an unreadable interleave, so each build streams to its own log file: +# docker/build-logs/-.log (follow live with: tail -f docker/build-logs/*.log) +set -euo pipefail + +REGISTRY="${REGISTRY:-gitea.psi.ch/leonarski_f}" +TAG="${TAG:-2607b}" +PUSH="${PUSH:-0}" +JOBS="${JOBS:-2}" + +here="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +logdir="$here/build-logs" +mkdir -p "$logdir" + +all=(rocky8 rocky9 ubuntu2204 ubuntu2404) +variants=("$@"); [ "${#variants[@]}" -eq 0 ] && variants=("${all[@]}") + +for v in "${variants[@]}"; do + [ -f "$here/$v/Dockerfile" ] || { echo "!! no Dockerfile for variant '$v' ($here/$v/Dockerfile)"; exit 1; } +done + +# One image build. Records OK/FAIL to a status file so the parent can summarise after wait; never +# returns non-zero itself, so the job-pool's `wait -n` under `set -e` stays happy. +build_one() { + local v="$1" + local img="$REGISTRY/jfjoch_$v:$TAG" + local log="$logdir/${v}-${TAG}.log" + if docker build --pull -t "$img" "$here/$v" >"$log" 2>&1; then + echo OK > "$logdir/${v}-${TAG}.status" + else + echo FAIL > "$logdir/${v}-${TAG}.status" + fi +} + +echo "==> building ${#variants[@]} image(s): ${variants[*]}" +echo "==> JOBS=$JOBS TAG=$TAG REGISTRY=$REGISTRY" +echo "==> logs in $logdir/ (follow: tail -f $logdir/*.log)" +echo + +# Launch with a concurrency cap of JOBS. `wait -n` reaps exactly one finished build before the next +# is started once we are at the cap; the trailing `wait` drains the rest. +running=0 +for v in "${variants[@]}"; do + if [ "$running" -ge "$JOBS" ]; then wait -n || true; running=$((running - 1)); fi + : > "$logdir/${v}-${TAG}.status" # clear any stale status from a previous run + echo "==> [$v] started -> $REGISTRY/jfjoch_$v:$TAG" + build_one "$v" & + running=$((running + 1)) +done +wait + +echo +echo "===================== build summary =====================" +rc=0; ok_variants=() +for v in "${variants[@]}"; do + s="$(cat "$logdir/${v}-${TAG}.status" 2>/dev/null || echo '??')" + printf " %-12s %s\n" "$v" "$s" + if [ "$s" = OK ]; then ok_variants+=("$v"); else rc=1; fi +done +[ "$rc" -eq 0 ] || echo "!! one or more builds FAILED -- see the per-variant logs above" + +echo +if [ "$PUSH" = "1" ]; then + for v in "${ok_variants[@]}"; do + img="$REGISTRY/jfjoch_$v:$TAG" + echo "==> push $img" + docker push "$img" + done +else + echo "Built locally. To push the successful ones (after 'docker login gitea.psi.ch'):" + for v in "${ok_variants[@]}"; do echo " docker push $REGISTRY/jfjoch_$v:$TAG"; done +fi + +exit "$rc" diff --git a/docker/build_in_rocky9.sh b/docker/build_in_rocky9.sh new file mode 100755 index 00000000..921f285d --- /dev/null +++ b/docker/build_in_rocky9.sh @@ -0,0 +1,48 @@ +#!/usr/bin/env bash +# Manually build the viewer inside one of the build-environment images and inspect its runtime .so +# deps -- the quick way to confirm the static-libdbus + glib-off changes actually dropped the +# dbus/systemd/glib/selinux tail. Defaults to the rocky9 image built by build_images.sh. +# +# The repo is mounted read-write at /workspace, but the CMake build tree lives in the container's +# /tmp (thrown away with --rm), so nothing root-owned is left in your working copy. Building needs +# network (FetchContent pulls curl/zstd/hdf5/...) but NOT a GPU. +# +# Usage: +# docker/build_in_rocky9.sh # build viewer in jfjoch_rocky9:$TAG, then ldd it +# docker/build_in_rocky9.sh rocky8 # use a different variant's image +# TAG=2607b docker/build_in_rocky9.sh +set -euo pipefail + +REGISTRY="${REGISTRY:-gitea.psi.ch/leonarski_f}" +TAG="${TAG:-2607b}" +VARIANT="${1:-rocky9}" +IMG="$REGISTRY/jfjoch_$VARIANT:$TAG" + +repo="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" # repo root = docker/.. + +echo "==> image $IMG" +echo "==> repo $repo (mounted at /workspace; build tree in container /tmp)" + +docker run --rm \ + -v "$repo":/workspace -w /workspace \ + "$IMG" bash -lc ' + set -eux + git config --global --add safe.directory /workspace || true + cmake -G Ninja -S /workspace -B /tmp/build \ + -DJFJOCH_VIEWER_ONLY=ON -DJFJOCH_USE_CUDA=ON \ + -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_CXX_FLAGS="-march=x86-64-v3" -DCMAKE_C_FLAGS="-march=x86-64-v3" + ninja -C /tmp/build -j"$(nproc)" jfjoch_viewer + bin=/tmp/build/viewer/jfjoch_viewer + + echo "===================== ldd jfjoch_viewer =====================" + ldd "$bin" | sort + echo "===================== shed-lib check ========================" + # These should ALL be gone now (libpcre2-16 is Qt s own and legitimately stays). + shed="dbus|systemd|zstd|lz4|selinux|glib|gthread|libpcre\.|libpcre2-8|gcrypt|gpg-error" + if ldd "$bin" | grep -Ei "$shed"; then + echo ">>> UNEXPECTED: the above libs are still linked" + else + echo ">>> CLEAN: dbus / systemd / glib / selinux tail all gone" + fi + ' -- 2.54.0 From b6c2d5848da9d36681001c25c01184a415c15521 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 13 Jul 2026 15:52:50 +0200 Subject: [PATCH 03/90] docker: build the static libdbus with -fPIC so Qt's PIE tools link it The static libdbus linked fine into libQt6DBus.a, but Qt links its D-Bus code generators (qdbusxml2cpp/qdbuscpp2xml) as position-independent executables (-fPIE -pie), and the default autotools dbus build produces non-PIC objects: ld: libdbus-1.a(...dbus-address.o): relocation R_X86_64_32 against `.rodata' can not be used when making a PIE object; recompile with -fPIE Build libdbus with --with-pic and CFLAGS=-fPIC so its archive objects are position-independent and link into both the static Qt libs and the PIE tools. Reproduced (non-PIC -> R_X86_64_32 fail) and fixed (PIC -> PIE whole-archive link passes) in a rocky8 container before applying. Co-Authored-By: Claude Opus 4.8 (1M context) --- docker/rocky8/Dockerfile | 10 ++++++---- docker/rocky9/Dockerfile | 10 ++++++---- docker/ubuntu2204/Dockerfile | 10 ++++++---- docker/ubuntu2404/Dockerfile | 10 ++++++---- 4 files changed, 24 insertions(+), 16 deletions(-) diff --git a/docker/rocky8/Dockerfile b/docker/rocky8/Dockerfile index 4d395230..357d031f 100644 --- a/docker/rocky8/Dockerfile +++ b/docker/rocky8/Dockerfile @@ -103,19 +103,21 @@ ENV PKG_CONFIG_PATH=${OPENSSL_ROOT_DIR}/lib/pkgconfig:${OPENSSL_ROOT_DIR}/lib64/ # D-Bus feature -- it is a pure session-bus client (registers a name + exports an adaptor), so the # daemon-side systemd/selinux features are irrelevant. Static-only install (no .so); the Qt build # finds it through find_package(DBus1) via the -DDBus1_DIR hint (the CMake package config dbus -# installs, which imports libdbus-1.a). expat is only a configure-time dep of -# the dbus daemon; the client libdbus-1 links neither expat nor systemd. +# installs, which imports libdbus-1.a). Built -fPIC (--with-pic) so the archive links into Qt's PIE +# tools (qdbusxml2cpp). expat is only a configure-time dep of the dbus daemon; the client libdbus-1 +# links neither expat nor systemd. RUN set -eux; \ cd /tmp; \ curl -LO https://dbus.freedesktop.org/releases/dbus/dbus-${DBUS_VERSION}.tar.xz; \ tar -xf dbus-${DBUS_VERSION}.tar.xz; \ cd dbus-${DBUS_VERSION}; \ ./configure --prefix=/opt/dbus-${DBUS_VERSION}-static \ - --enable-static --disable-shared \ + --enable-static --disable-shared --with-pic \ --disable-systemd --without-systemdsystemunitdir \ --disable-selinux --disable-apparmor --disable-libaudit \ --disable-tests --disable-asserts \ - --disable-doxygen-docs --disable-xml-docs --disable-ducktype-docs; \ + --disable-doxygen-docs --disable-xml-docs --disable-ducktype-docs \ + CFLAGS="-fPIC"; \ make -j"$(nproc)"; \ make install; \ cd /; rm -rf /tmp/dbus-${DBUS_VERSION} /tmp/dbus-${DBUS_VERSION}.tar.xz diff --git a/docker/rocky9/Dockerfile b/docker/rocky9/Dockerfile index d0e27deb..f5328bff 100644 --- a/docker/rocky9/Dockerfile +++ b/docker/rocky9/Dockerfile @@ -106,19 +106,21 @@ ENV CXX=${GCC_TOOLSET_ROOT}/usr/bin/g++ # D-Bus feature -- it is a pure session-bus client (registers a name + exports an adaptor), so the # daemon-side systemd/selinux features are irrelevant. Static-only install (no .so); the Qt build # finds it through find_package(DBus1) via the -DDBus1_DIR hint (the CMake package config dbus -# installs, which imports libdbus-1.a). expat is only a configure-time dep of -# the dbus daemon; the client libdbus-1 links neither expat nor systemd. +# installs, which imports libdbus-1.a). Built -fPIC (--with-pic) so the archive links into Qt's PIE +# tools (qdbusxml2cpp). expat is only a configure-time dep of the dbus daemon; the client libdbus-1 +# links neither expat nor systemd. RUN set -eux; \ cd /tmp; \ curl -LO https://dbus.freedesktop.org/releases/dbus/dbus-${DBUS_VERSION}.tar.xz; \ tar -xf dbus-${DBUS_VERSION}.tar.xz; \ cd dbus-${DBUS_VERSION}; \ ./configure --prefix=/opt/dbus-${DBUS_VERSION}-static \ - --enable-static --disable-shared \ + --enable-static --disable-shared --with-pic \ --disable-systemd --without-systemdsystemunitdir \ --disable-selinux --disable-apparmor --disable-libaudit \ --disable-tests --disable-asserts \ - --disable-doxygen-docs --disable-xml-docs --disable-ducktype-docs; \ + --disable-doxygen-docs --disable-xml-docs --disable-ducktype-docs \ + CFLAGS="-fPIC"; \ make -j"$(nproc)"; \ make install; \ cd /; rm -rf /tmp/dbus-${DBUS_VERSION} /tmp/dbus-${DBUS_VERSION}.tar.xz diff --git a/docker/ubuntu2204/Dockerfile b/docker/ubuntu2204/Dockerfile index b0c790b8..7fcadefd 100644 --- a/docker/ubuntu2204/Dockerfile +++ b/docker/ubuntu2204/Dockerfile @@ -130,19 +130,21 @@ RUN set -eux; \ # D-Bus feature -- it is a pure session-bus client (registers a name + exports an adaptor), so the # daemon-side systemd/selinux features are irrelevant. Static-only install (no .so); the Qt build # finds it through find_package(DBus1) via the -DDBus1_DIR hint (the CMake package config dbus -# installs, which imports libdbus-1.a). expat is only a configure-time dep of -# the dbus daemon; the client libdbus-1 links neither expat nor systemd. +# installs, which imports libdbus-1.a). Built -fPIC (--with-pic) so the archive links into Qt's PIE +# tools (qdbusxml2cpp). expat is only a configure-time dep of the dbus daemon; the client libdbus-1 +# links neither expat nor systemd. RUN set -eux; \ cd /tmp; \ curl -LO https://dbus.freedesktop.org/releases/dbus/dbus-${DBUS_VERSION}.tar.xz; \ tar -xf dbus-${DBUS_VERSION}.tar.xz; \ cd dbus-${DBUS_VERSION}; \ ./configure --prefix=/opt/dbus-${DBUS_VERSION}-static \ - --enable-static --disable-shared \ + --enable-static --disable-shared --with-pic \ --disable-systemd --without-systemdsystemunitdir \ --disable-selinux --disable-apparmor --disable-libaudit \ --disable-tests --disable-asserts \ - --disable-doxygen-docs --disable-xml-docs --disable-ducktype-docs; \ + --disable-doxygen-docs --disable-xml-docs --disable-ducktype-docs \ + CFLAGS="-fPIC"; \ make -j"$(nproc)"; \ make install; \ cd /; rm -rf /tmp/dbus-${DBUS_VERSION} /tmp/dbus-${DBUS_VERSION}.tar.xz diff --git a/docker/ubuntu2404/Dockerfile b/docker/ubuntu2404/Dockerfile index bd7be0b1..f0b649e7 100644 --- a/docker/ubuntu2404/Dockerfile +++ b/docker/ubuntu2404/Dockerfile @@ -117,19 +117,21 @@ RUN set -eux; \ # D-Bus feature -- it is a pure session-bus client (registers a name + exports an adaptor), so the # daemon-side systemd/selinux features are irrelevant. Static-only install (no .so); the Qt build # finds it through find_package(DBus1) via the -DDBus1_DIR hint (the CMake package config dbus -# installs, which imports libdbus-1.a). expat is only a configure-time dep of -# the dbus daemon; the client libdbus-1 links neither expat nor systemd. +# installs, which imports libdbus-1.a). Built -fPIC (--with-pic) so the archive links into Qt's PIE +# tools (qdbusxml2cpp). expat is only a configure-time dep of the dbus daemon; the client libdbus-1 +# links neither expat nor systemd. RUN set -eux; \ cd /tmp; \ curl -LO https://dbus.freedesktop.org/releases/dbus/dbus-${DBUS_VERSION}.tar.xz; \ tar -xf dbus-${DBUS_VERSION}.tar.xz; \ cd dbus-${DBUS_VERSION}; \ ./configure --prefix=/opt/dbus-${DBUS_VERSION}-static \ - --enable-static --disable-shared \ + --enable-static --disable-shared --with-pic \ --disable-systemd --without-systemdsystemunitdir \ --disable-selinux --disable-apparmor --disable-libaudit \ --disable-tests --disable-asserts \ - --disable-doxygen-docs --disable-xml-docs --disable-ducktype-docs; \ + --disable-doxygen-docs --disable-xml-docs --disable-ducktype-docs \ + CFLAGS="-fPIC"; \ make -j"$(nproc)"; \ make install; \ cd /; rm -rf /tmp/dbus-${DBUS_VERSION} /tmp/dbus-${DBUS_VERSION}.tar.xz -- 2.54.0 From d90d498cdbce8f3bc71053a8a9efb2d0ce324c2a Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 13 Jul 2026 19:13:29 +0200 Subject: [PATCH 04/90] docker: put the static libdbus on CMAKE_PREFIX_PATH for Qt6::DBus consumers The Qt build found dbus via -DDBus1_DIR, but a static Qt records the dependency and every consumer re-runs find_package(DBus1) through Qt's exported config (Qt6Gui -> Qt6DBus -> find_dependency(WrapDBus1)). Without the dbus prefix on the consumer's search path that fails: Could NOT find WrapDBus1 (missing: DBus1_LIBRARY DBus1_INCLUDE_DIR ...) Qt6DBus could not be found because dependency WrapDBus1 could not be found. Add /opt/dbus--static to the image's CMAKE_PREFIX_PATH env so the viewer build (and the CI .deb/.rpm and tarball builds) resolve it. Verified against the rocky9 image: find_package(Qt6 ... DBus) fails with the old path and passes with dbus on the env CMAKE_PREFIX_PATH. Also make build_in_rocky9.sh self-sufficient on images built before this fix: it discovers /opt/dbus-*-static/lib/cmake/DBus1 and passes -DDBus1_DIR, so the viewer verification works without rebuilding the image first. Co-Authored-By: Claude Opus 4.8 (1M context) --- docker/build_in_rocky9.sh | 5 +++++ docker/rocky8/Dockerfile | 6 ++++-- docker/rocky9/Dockerfile | 6 ++++-- docker/ubuntu2204/Dockerfile | 6 ++++-- docker/ubuntu2404/Dockerfile | 6 ++++-- 5 files changed, 21 insertions(+), 8 deletions(-) diff --git a/docker/build_in_rocky9.sh b/docker/build_in_rocky9.sh index 921f285d..71aa71b0 100755 --- a/docker/build_in_rocky9.sh +++ b/docker/build_in_rocky9.sh @@ -28,8 +28,13 @@ docker run --rm \ "$IMG" bash -lc ' set -eux git config --global --add safe.directory /workspace || true + # If this image predates the CMAKE_PREFIX_PATH fix, point cmake at the static libdbus config + # dir so the viewer (Qt6::DBus -> find_package(DBus1)) resolves WrapDBus1. Harmless once the + # image ENV already carries it. + dbus_dir=$(ls -d /opt/dbus-*-static/lib/cmake/DBus1 2>/dev/null | head -1) cmake -G Ninja -S /workspace -B /tmp/build \ -DJFJOCH_VIEWER_ONLY=ON -DJFJOCH_USE_CUDA=ON \ + ${dbus_dir:+-DDBus1_DIR=$dbus_dir} \ -DCMAKE_BUILD_TYPE=Release \ -DCMAKE_CXX_FLAGS="-march=x86-64-v3" -DCMAKE_C_FLAGS="-march=x86-64-v3" ninja -C /tmp/build -j"$(nproc)" jfjoch_viewer diff --git a/docker/rocky8/Dockerfile b/docker/rocky8/Dockerfile index 357d031f..f3db6b7f 100644 --- a/docker/rocky8/Dockerfile +++ b/docker/rocky8/Dockerfile @@ -168,8 +168,10 @@ RUN set -eux; \ cmake --install eigen-build; \ cd /; rm -rf /tmp/eigen /tmp/eigen-build -# Make Qt and Eigen discoverable by CMake -ENV CMAKE_PREFIX_PATH=/opt/qt-${QT_VERSION}-static:/opt/eigen-3.4 +# Make Qt, the static libdbus, and Eigen discoverable by CMake. The static libdbus prefix is +# required here (not only at Qt build time): every Qt6::DBus consumer re-runs find_package(DBus1) +# through Qt's exported config, so DBus1Config.cmake must be on CMAKE_PREFIX_PATH for the viewer build. +ENV CMAKE_PREFIX_PATH=/opt/qt-${QT_VERSION}-static:/opt/dbus-${DBUS_VERSION}-static:/opt/eigen-3.4 # Set workdir for your project WORKDIR /workspace diff --git a/docker/rocky9/Dockerfile b/docker/rocky9/Dockerfile index f5328bff..57a513f8 100644 --- a/docker/rocky9/Dockerfile +++ b/docker/rocky9/Dockerfile @@ -207,8 +207,10 @@ RUN set -eux; \ cmake --install eigen-build; \ cd /; rm -rf /tmp/eigen /tmp/eigen-build -# Make Qt and Eigen discoverable by CMake -ENV CMAKE_PREFIX_PATH=/opt/qt-${QT_VERSION}-static:/opt/eigen-3.4 +# Make Qt, the static libdbus, and Eigen discoverable by CMake. The static libdbus prefix is +# required here (not only at Qt build time): every Qt6::DBus consumer re-runs find_package(DBus1) +# through Qt's exported config, so DBus1Config.cmake must be on CMAKE_PREFIX_PATH for the viewer build. +ENV CMAKE_PREFIX_PATH=/opt/qt-${QT_VERSION}-static:/opt/dbus-${DBUS_VERSION}-static:/opt/eigen-3.4 ENV LANG=en_US.UTF-8 ENV LANGUAGE=en_US:en diff --git a/docker/ubuntu2204/Dockerfile b/docker/ubuntu2204/Dockerfile index 7fcadefd..56588301 100644 --- a/docker/ubuntu2204/Dockerfile +++ b/docker/ubuntu2204/Dockerfile @@ -196,8 +196,10 @@ RUN set -eux; \ cmake --install eigen-build; \ cd /; rm -rf /tmp/eigen /tmp/eigen-build -# Make Qt and Eigen discoverable by CMake -ENV CMAKE_PREFIX_PATH=/opt/qt-${QT_VERSION}-static:/opt/eigen-3.4 +# Make Qt, the static libdbus, and Eigen discoverable by CMake. The static libdbus prefix is +# required here (not only at Qt build time): every Qt6::DBus consumer re-runs find_package(DBus1) +# through Qt's exported config, so DBus1Config.cmake must be on CMAKE_PREFIX_PATH for the viewer build. +ENV CMAKE_PREFIX_PATH=/opt/qt-${QT_VERSION}-static:/opt/dbus-${DBUS_VERSION}-static:/opt/eigen-3.4 # Set workdir for your project WORKDIR /workspace diff --git a/docker/ubuntu2404/Dockerfile b/docker/ubuntu2404/Dockerfile index f0b649e7..ca9678c7 100644 --- a/docker/ubuntu2404/Dockerfile +++ b/docker/ubuntu2404/Dockerfile @@ -183,8 +183,10 @@ RUN set -eux; \ cmake --install eigen-build; \ cd /; rm -rf /tmp/eigen /tmp/eigen-build -# Make Qt and Eigen discoverable by CMake -ENV CMAKE_PREFIX_PATH=/opt/qt-${QT_VERSION}-static:/opt/eigen-3.4 +# Make Qt, the static libdbus, and Eigen discoverable by CMake. The static libdbus prefix is +# required here (not only at Qt build time): every Qt6::DBus consumer re-runs find_package(DBus1) +# through Qt's exported config, so DBus1Config.cmake must be on CMAKE_PREFIX_PATH for the viewer build. +ENV CMAKE_PREFIX_PATH=/opt/qt-${QT_VERSION}-static:/opt/dbus-${DBUS_VERSION}-static:/opt/eigen-3.4 # Set workdir for your project WORKDIR /workspace -- 2.54.0 From 3df9ad631bed8741cc0f0293406d3e4e3daec875 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 13 Jul 2026 19:26:45 +0200 Subject: [PATCH 05/90] docker: fix build_in_rocky9.sh shed-check for RHEL9's glib/selinux sources The old check flagged glib/selinux/pcre as failures, inheriting a rocky8-host assumption (glib is Qt-only, selinux is dbus-only) that does not hold on RHEL9: - libglib is pulled by the system harfbuzz/fontconfig/freetype text stack (glib-linked on RHEL9), not by Qt. QT_FEATURE_glib=OFF is honoured (Qt's enabled_features has no glib) but cannot drop the font stack's libglib. - libselinux/libpcre2-8 come from the krb5/GSSAPI stack (libgssapi_krb5 -> libselinux) that static curl links for auth -- kept for the reverse-proxy Kerberos plan, same bucket as libcrypto. So hard-fail only on the dbus/systemd tail we actually remove (libdbus/ libsystemd/zstd/lz4/cap/gcrypt/gpg-error), and just report glib/selinux as expected residuals with the reason. Verified in the rocky9 image: that tail is gone; the survivors trace to the font stack and krb5. Co-Authored-By: Claude Opus 4.8 (1M context) --- docker/build_in_rocky9.sh | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/docker/build_in_rocky9.sh b/docker/build_in_rocky9.sh index 71aa71b0..6f02a7a9 100755 --- a/docker/build_in_rocky9.sh +++ b/docker/build_in_rocky9.sh @@ -42,12 +42,20 @@ docker run --rm \ echo "===================== ldd jfjoch_viewer =====================" ldd "$bin" | sort - echo "===================== shed-lib check ========================" - # These should ALL be gone now (libpcre2-16 is Qt s own and legitimately stays). - shed="dbus|systemd|zstd|lz4|selinux|glib|gthread|libpcre\.|libpcre2-8|gcrypt|gpg-error" - if ldd "$bin" | grep -Ei "$shed"; then - echo ">>> UNEXPECTED: the above libs are still linked" + echo "========== dbus/systemd tail (what the static libdbus removes) ==========" + # These only ever entered via the system libdbus-1.so -> libsystemd chain, so they MUST be gone. + tail_re="libdbus-1|libsystemd|libzstd|liblz4|libgcrypt|libgpg-error|libcap\." + if ldd "$bin" | grep -Ei "$tail_re"; then + echo ">>> FAIL: the dbus/systemd tail is still linked" else - echo ">>> CLEAN: dbus / systemd / glib / selinux tail all gone" + echo ">>> OK: libdbus/libsystemd/zstd/lz4/cap/gcrypt/gpg-error all gone" fi + echo "========== residual glib/selinux (expected on RHEL9, NOT from us) ==========" + # libglib here is pulled by the system harfbuzz/fontconfig/freetype text stack (glib-linked on + # RHEL9), not by Qt -- QT_FEATURE_glib=OFF is honoured (Qt is off the glib event loop) but cannot + # drop the font-stack glib. libselinux/libpcre2-8 come from the krb5/GSSAPI stack that static + # curl links for auth (kept for the reverse-proxy Kerberos plan) -- same bucket as libcrypto. + # Shedding these would mean static-linking the whole font stack / dropping krb5. (libpcre2-16 is + # Qt-owned and legitimately stays.) + ldd "$bin" | grep -Ei "libglib|gthread|selinux|libpcre\.|libpcre2-8" || echo "(none)" ' -- 2.54.0 From 87778c6500e0a919ed9c0a9dbd38a738901060a5 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 13 Jul 2026 19:46:48 +0200 Subject: [PATCH 06/90] reader: treat negative total_flux as unknown/absent Serial-stills master files can store /entry/instrument/beam/total_flux as a -1.0 "unknown flux" sentinel (seen in the OCP dark dataset). DatasetSettings rejects values below 0, so opening such a file hard-failed with "Input parameter below min (Total flux)". Reset a negative flux to nullopt so it is treated as absent instead of aborting the read. Co-Authored-By: Claude Opus 4.8 (1M context) --- reader/HDF5MetadataSource.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 61076f08..95bfc19a 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -611,7 +611,10 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen if (master_file->Exists("/entry/instrument/attenuator")) dataset->experiment.AttenuatorTransmission( master_file->GetOptFloat("/entry/instrument/attenuator/attenuator_transmission")); - dataset->experiment.TotalFlux(master_file->GetOptFloat("/entry/instrument/beam/total_flux")); + auto total_flux = master_file->GetOptFloat("/entry/instrument/beam/total_flux"); + if (total_flux.has_value() && total_flux.value() < 0) + total_flux.reset(); // negative value is an "unknown flux" sentinel; treat as absent + dataset->experiment.TotalFlux(total_flux); if (master_file->Exists("/entry/azint") && master_file->Exists("/entry/azint/bin_to_q")) { HDF5DataSet bin_to_q_dataset(*master_file, "/entry/azint/bin_to_q"); -- 2.54.0 From 2d7ddb1b678b578b344de12a25e77a31d398d19c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 13 Jul 2026 19:46:48 +0200 Subject: [PATCH 07/90] ScaleOnTheFly: resolve the stills indexing ambiguity per image vs a reference Serial stills index each crystal in one of the merohedrally-equivalent hands at random, so the ambiguity must be broken per image, not globally as the rotation path (RotationScaleMerge/ChooseReindex) does. When an external reference MTZ is supplied, IndexAndRefine::ReferenceIntensities now builds the per-image scaler with resolve_ambiguity=true: for each image it picks the reindexing (identity or a twin law from ReindexAmbiguityOperators) whose partiality/Lorentz-corrected intensities correlate best with the reference, then reindexes that image's reflections before scaling. Off for the self-merge scaling pass (no trusted reference to break the tie). Holohedral crystals (no twin law, e.g. lysozyme P43212) are a no-op. Validated on OCP P3221 (twin law -x,-y,z) with a 7ZSJ reference: CCref 24->48% and CC1/2 24->64% at 30k images; CC1/2 96.9% / CCref 76% at 150k. Co-Authored-By: Claude Opus 4.8 (1M context) --- image_analysis/IndexAndRefine.cpp | 4 +- image_analysis/scale_merge/ScaleOnTheFly.cpp | 59 +++++++++++++++++++- image_analysis/scale_merge/ScaleOnTheFly.h | 13 ++++- 3 files changed, 73 insertions(+), 3 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 40bfaf4f..b0de29fa 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -431,7 +431,9 @@ std::optional IndexAndRefine::FinalizeRotationIndexing() } IndexAndRefine &IndexAndRefine::ReferenceIntensities(std::vector &reference) { - scaling_engine = std::make_unique(experiment, reference); + // An external reference is trusted to be in the correct hand, so let the per-image scaler also break + // the merohedral indexing ambiguity against it (serial stills index each crystal independently). + scaling_engine = std::make_unique(experiment, reference, /*resolve_ambiguity=*/true); return *this; } diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index 175b2aac..a2b3c8f0 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -96,13 +96,67 @@ namespace { }; } -ScaleOnTheFly::ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref) +ScaleOnTheFly::ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref, + bool resolve_ambiguity) : s(x.GetScalingSettings()), hkl_key_generator(s.GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1)) { for (const auto &r: ref) { const auto key = hkl_key_generator(r); reference_data[key] = r.I; } + if (resolve_ambiguity && x.GetUnitCell().has_value() && x.GetSpaceGroupNumber().has_value()) + ambiguity_ops = ReindexAmbiguityOperators(*x.GetUnitCell(), + static_cast(*x.GetSpaceGroupNumber())); +} + +// Serial stills index each crystal in one of the merohedrally-equivalent hands at random; pick, for this +// image alone, the reindexing whose intensities correlate best with the external reference and apply it. +// The correlation uses the partiality/Lorentz-corrected intensity (scale-invariant, so the unknown +// per-image scale G cancels) matched to the reference by ASU index under the candidate hand. +void ScaleOnTheFly::ResolveIndexingAmbiguity(std::vector &reflections) const { + auto reference_cc = [&](const gemmi::Op &op) -> std::pair { + double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0; + size_t n = 0; + for (const auto &r: reflections) { + if (!Accept(r) || r.partiality < s.GetMinPartiality()) + continue; + if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f || r.partiality <= 0.0f) + continue; + const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}}); + const auto it = reference_data.find(hkl_key_generator(h[0], h[1], h[2])); + if (it == reference_data.end()) + continue; + const double x = static_cast(r.I) * SafeInv(r.rlp, 1.0) / r.partiality; + const double y = it->second; + if (!std::isfinite(x) || !std::isfinite(y)) + continue; + sx += x; sy += y; sxx += x * x; syy += y * y; sxy += x * y; + ++n; + } + if (n < MIN_REFLECTIONS) + return {NAN, n}; + const double nd = static_cast(n); + const double cov = sxy - sx * sy / nd; + const double vx = sxx - sx * sx / nd; + const double vy = syy - sy * sy / nd; + return {(vx > 0 && vy > 0) ? cov / std::sqrt(vx * vy) : NAN, n}; + }; + + const double identity_cc = reference_cc(gemmi::Op::identity()).first; + const gemmi::Op *best = nullptr; + double best_cc = identity_cc; + for (const auto &op: ambiguity_ops) { + const double cc = reference_cc(op).first; + if (std::isfinite(cc) && (!std::isfinite(best_cc) || cc > best_cc)) { + best_cc = cc; + best = &op; + } + } + if (best != nullptr) + for (auto &r: reflections) { + const gemmi::Op::Miller h = best->apply_to_hkl({{r.h, r.k, r.l}}); + r.h = h[0]; r.k = h[1]; r.l = h[2]; + } } bool ScaleOnTheFly::Accept(const Reflection &r) const { @@ -169,6 +223,9 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { if (integration_outcome.reflections.empty()) return; + if (!ambiguity_ops.empty()) + ResolveIndexingAmbiguity(integration_outcome.reflections); + auto start = std::chrono::steady_clock::now(); ScaleOnTheFlyResult result{ .B = 0.0, .G = 1.0 }; diff --git a/image_analysis/scale_merge/ScaleOnTheFly.h b/image_analysis/scale_merge/ScaleOnTheFly.h index 5e5f4454..6166b4f1 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.h +++ b/image_analysis/scale_merge/ScaleOnTheFly.h @@ -5,10 +5,12 @@ #include "HKLKey.h" #include "Merge.h" +#include "ReindexAmbiguity.h" #include "../../common/DiffractionExperiment.h" #include "../IntegrationOutcome.h" #include +#include struct ScaleOnTheFlyResult { @@ -31,11 +33,20 @@ class ScaleOnTheFly { const ScalingSettings s; const HKLKeyGenerator hkl_key_generator; std::map reference_data; + // Twin-law reindexing operators for the merohedral indexing ambiguity (empty = no ambiguity, or + // resolution disabled). Populated only when the ref is a trusted external reference (see constructor). + std::vector ambiguity_ops; bool Accept(const Reflection &r) const; [[nodiscard]] std::pair CalculateGlobalCC(const std::vector &reflections) const; + // Reindex this image's reflections into the ambiguity hand that best correlates with the reference. + void ResolveIndexingAmbiguity(std::vector &reflections) const; public: - ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref); + // resolve_ambiguity: when the reference is an external, correctly-handed dataset (not the data's own + // running merge), pick per image the reindexing that agrees best with it - serial stills index each + // crystal independently, so the merohedral ambiguity must be broken per image, not globally. + ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref, + bool resolve_ambiguity = false); void Scale(IntegrationOutcome &integration_outcome) const; void Scale(std::vector &integration_outcome, size_t nthreads = 0) const; -- 2.54.0 From 2f3041c44e416b939347dcafd6c1267c74f96b19 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 13 Jul 2026 22:02:49 +0200 Subject: [PATCH 08/90] rugnux: add --min-pix-per-spot CLI option (spot finding) The minimum connected strong-pixel count per spot was hardcoded to 2 (SpotFindingSettings::min_pix_per_spot). Expose it on the CLI so serial-stills data can be run with min-pix 1 (plus a higher --spot-threshold), matching the CrystFEL peakfinder8 practice that trades noisier spots for a higher indexing rate. Default unchanged (2). Co-Authored-By: Claude Opus 4.8 (1M context) --- rugnux/rugnux_cli.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 30402b28..9fe6ff3e 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -58,6 +58,7 @@ void print_usage() { std::cout << " Spot finding" << std::endl; std::cout << " --spot-sigma Noise sigma level for spot finding (default: 3.0)" << std::endl; std::cout << " --spot-threshold Photon count threshold for spot finding (default: 10)" << std::endl; + std::cout << " --min-pix-per-spot Minimum connected strong pixels per spot (default: 2; serial data can index better with 1 + a higher --spot-threshold)" << std::endl; std::cout << " --spot-high-resolution High resolution limit for spot finding (default: 1.5)" << std::endl; std::cout << " --max-spots Max spot count (default: 250)" << std::endl; std::cout << " --detect-ice-rings[=on|off] Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling/merging; overrides the dataset/master-file setting (default: use dataset value)" << std::endl; @@ -128,6 +129,7 @@ void print_usage() { enum { OPT_SPOT_SIGMA = 1000, OPT_SPOT_THRESHOLD, + OPT_MIN_PIX_PER_SPOT, OPT_SPOT_RESOLUTION, OPT_MAX_SPOTS, OPT_MIN_PARTIALITY, @@ -224,6 +226,7 @@ static option long_options[] = { {"spot-sigma", required_argument, nullptr, OPT_SPOT_SIGMA}, {"spot-threshold", required_argument, nullptr, OPT_SPOT_THRESHOLD}, + {"min-pix-per-spot", required_argument, nullptr, OPT_MIN_PIX_PER_SPOT}, {"spot-high-resolution", required_argument, nullptr, OPT_SPOT_RESOLUTION}, {"max-spots", required_argument, nullptr, OPT_MAX_SPOTS}, {"min-partiality", required_argument, nullptr, OPT_MIN_PARTIALITY}, @@ -454,6 +457,7 @@ int main(int argc, char **argv) { std::optional max_spot_count_override; float sigma_spot_finding = 3.0; int64_t photon_count_threshold_spot_finding = 10; + int64_t min_pix_per_spot = 2; bool refine_bfactor = false; std::string ref_mtz; std::string ref_column; @@ -647,6 +651,10 @@ int main(int argc, char **argv) { logger.Info("Photon-count threshold level for spot finding set to {:d}", photon_count_threshold_spot_finding); break; + case OPT_MIN_PIX_PER_SPOT: + min_pix_per_spot = atoi(optarg); + logger.Info("Minimum pixels per spot set to {:d}", min_pix_per_spot); + break; case OPT_SPOT_RESOLUTION: d_min_spot_finding = atof(optarg); logger.Info("High resolution limit for spot finding set to {:.2f} A", d_min_spot_finding); @@ -1299,6 +1307,7 @@ int main(int argc, char **argv) { spot_settings.high_resolution_limit = d_min_spot_finding; spot_settings.signal_to_noise_threshold = sigma_spot_finding; spot_settings.photon_count_threshold = photon_count_threshold_spot_finding; + spot_settings.min_pix_per_spot = min_pix_per_spot; if (d_min_spot_finding > 0.0f) spot_settings.high_resolution_limit = d_min_spot_finding; -- 2.54.0 From 097b9a29d16072aaf26d67d303c0d565d0746329 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 13 Jul 2026 23:14:27 +0200 Subject: [PATCH 09/90] rugnux: add --spot-low-resolution CLI option (exclude beam halo) The spot-finding low-resolution limit was hardcoded at 50 A (r ~36 px), which admits the direct-beam halo/overload zone (24-50 A) on weakly-diffracting serial data. On KR2 (I222, 235 A axis, ~3.9 A diffraction) that halo is ~387 strong pixels/frame vs ~25 real Bragg pixels, swamping the signal ~15:1 (CrystFEL cuts it with min-res=75 px = 24 A). Expose the limit so it can be tightened, e.g. --spot-low-resolution 24. Default unchanged (50 A). Co-Authored-By: Claude Opus 4.8 (1M context) --- rugnux/rugnux_cli.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 9fe6ff3e..33d55ebc 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -60,6 +60,7 @@ void print_usage() { std::cout << " --spot-threshold Photon count threshold for spot finding (default: 10)" << std::endl; std::cout << " --min-pix-per-spot Minimum connected strong pixels per spot (default: 2; serial data can index better with 1 + a higher --spot-threshold)" << std::endl; std::cout << " --spot-high-resolution High resolution limit for spot finding (default: 1.5)" << std::endl; + std::cout << " --spot-low-resolution Low resolution limit for spot finding, in A (default: 50; lower it, e.g. 24, to exclude the direct-beam halo on weakly-diffracting serial data)" << std::endl; std::cout << " --max-spots Max spot count (default: 250)" << std::endl; std::cout << " --detect-ice-rings[=on|off] Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling/merging; overrides the dataset/master-file setting (default: use dataset value)" << std::endl; std::cout << std::endl; @@ -131,6 +132,7 @@ enum { OPT_SPOT_THRESHOLD, OPT_MIN_PIX_PER_SPOT, OPT_SPOT_RESOLUTION, + OPT_SPOT_LOW_RESOLUTION, OPT_MAX_SPOTS, OPT_MIN_PARTIALITY, OPT_MIN_IMAGE_CC, @@ -228,6 +230,7 @@ static option long_options[] = { {"spot-threshold", required_argument, nullptr, OPT_SPOT_THRESHOLD}, {"min-pix-per-spot", required_argument, nullptr, OPT_MIN_PIX_PER_SPOT}, {"spot-high-resolution", required_argument, nullptr, OPT_SPOT_RESOLUTION}, + {"spot-low-resolution", required_argument, nullptr, OPT_SPOT_LOW_RESOLUTION}, {"max-spots", required_argument, nullptr, OPT_MAX_SPOTS}, {"min-partiality", required_argument, nullptr, OPT_MIN_PARTIALITY}, {"capture-uncertainty", required_argument, nullptr, OPT_CAPTURE_UNCERTAINTY}, @@ -479,6 +482,7 @@ int main(int argc, char **argv) { std::optional intensity_format; // --scaling-output override; default lives in ScalingSettings (mmCIF) float d_min_spot_finding = 1.5; + float d_max_spot_finding = 0; // 0 = keep the SpotFindingSettings default (50 A) std::optional d_min_scale_merge; std::optional resolution_cutoff_method; // --resolution-cutoff cc-logistic|off std::optional resolution_cc_target; // --resolution-cc-target @@ -655,6 +659,10 @@ int main(int argc, char **argv) { min_pix_per_spot = atoi(optarg); logger.Info("Minimum pixels per spot set to {:d}", min_pix_per_spot); break; + case OPT_SPOT_LOW_RESOLUTION: + d_max_spot_finding = atof(optarg); + logger.Info("Low resolution limit for spot finding set to {:.1f} A", d_max_spot_finding); + break; case OPT_SPOT_RESOLUTION: d_min_spot_finding = atof(optarg); logger.Info("High resolution limit for spot finding set to {:.2f} A", d_min_spot_finding); @@ -1310,6 +1318,8 @@ int main(int argc, char **argv) { spot_settings.min_pix_per_spot = min_pix_per_spot; if (d_min_spot_finding > 0.0f) spot_settings.high_resolution_limit = d_min_spot_finding; + if (d_max_spot_finding > 0.0f) + spot_settings.low_resolution_limit = d_max_spot_finding; // Run the shared full-analysis workflow (rotation indexing + scaling/merging live in // Rugnux; the experiment above carries all algorithm settings). -- 2.54.0 From db8ae0452e39609a842c59e8093ed6bc1ee439ca Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 10:06:07 +0200 Subject: [PATCH 10/90] rugnux: fix --scale on a self-contained _process.h5 (reflections + error model) Offline --scale re-scales the reflections stored in a _process.h5 without re-integration, but was broken for a self-contained integrated snapshot: 1. HDF5MetadataSource::ReadReflections fell back per-image to the linked source pixel files for every non-indexed frame (a snapshot's master holds /entry/reflections only for the sparse indexed images). With the raw _data_NNNNNN.h5 absent or not co-located this threw an HDF5 error; when present it needlessly reopened multi-GB files thousands of times. Decide once whether the master is the authoritative reflection store and, if so, never fall back - a missing per-image group just means the frame has none. Legacy/VDS acquisitions (no /entry/reflections in the master) still resolve reflections lazily to the source data files. 2. The scale-only stills branch never fit the (a,b) error model, so it merged with the identity model - far worse intensities than the full pipeline (lyso8 CC1/2 76%->21%). Fit RefineErrorModel and honour --reject-outliers, mirroring Rugnux.cpp. Round-trip validated: full merge == integrate(--no-merge) + --scale (identical error model a=0.793 b=2.287, matching unique reflection counts), and --scale now reads reflections in ~2.5s with no raw-file access. Co-Authored-By: Claude Opus 4.8 (1M context) --- reader/HDF5MetadataSource.cpp | 9 ++++++++- rugnux/rugnux_cli.cpp | 12 ++++++++++++ 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 95bfc19a..40469725 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -1099,6 +1099,13 @@ std::vector HDF5MetadataSource::ReadReflections(size_t start std::vector ret; ret.reserve(end_image_val - start_image + 1); + // A self-contained integrated _process.h5 keeps all reflections in this master (one group per + // indexed image), so a missing per-image group means that image simply has none - never fall + // back to the linked source pixel files (which may be absent, and never hold a snapshot's + // reflections). A legacy/VDS acquisition has no /entry/reflections in the master and resolves + // reflections lazily from the source data files instead. + const bool master_reflections_authoritative = master_file->Exists("/entry/reflections"); + for (size_t img = start_image; img <= end_image_val; img++) { IntegrationOutcome outcome; @@ -1112,7 +1119,7 @@ std::vector HDF5MetadataSource::ReadReflections(size_t start HDF5ReadOnlyFile *meta_file = master_file.get(); size_t meta_image_id = img; std::string refl_group = fmt::format("/entry/reflections/image_{:06d}", img); - if (!master_file->Exists(refl_group)) { + if (!master_reflections_authoritative && !master_file->Exists(refl_group)) { const auto loc = ResolveMeta(static_cast(img)); meta_file = loc.file.get(); meta_image_id = loc.local_index; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 33d55ebc..e9f86dc1 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -986,6 +986,9 @@ int main(int argc, char **argv) { scaling_settings.MinCapturedFraction(min_captured_fraction_arg.value_or( (experiment.GetGoniometer().has_value() && !force_still) ? 0.7 : 0.0)); scaling_settings.MinCCForImage(min_image_cc / 100.0); // --min-image-cc is percent; the setting is a fraction + scaling_settings.OutlierRejectNsigma( + outlier_reject_nsigma.value_or( + (experiment.GetGoniometer().has_value() && !force_still) ? REJECT_OUTLIERS_DEFAULT_NSIGMA : 0.0)); if (intensity_format) scaling_settings.FileFormat(intensity_format.value()); experiment.ImportScalingSettings(scaling_settings); @@ -1031,6 +1034,15 @@ int main(int argc, char **argv) { } MergeOnTheFly merge_engine(experiment); merge_engine.ReferenceCell(experiment.GetUnitCell()); + // Fit the (a, b) error model from symmetry-mate scatter before merging, exactly as the full + // pipeline does (Rugnux.cpp). Without this the offline --scale merge would use the identity + // model and produce much worse stills intensities (no (b*I)^2 systematic term, no sigma floor). + merge_engine.RefineErrorModel(reflections); + if (merge_engine.ErrorModelActive()) + logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", merge_engine.ErrorModelA(), + merge_engine.ErrorModelB(), + merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0, + merge_engine.ErrorModelChi2()); for (size_t i = 0; i < reflections.size(); ++i) merge_engine.AddImage(reflections[i], static_cast(i)); merged_reflections = merge_engine.ExportReflections(); -- 2.54.0 From b0202a3aa2846c2c2587e44a5c3d4155e590e0c1 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 12:00:26 +0200 Subject: [PATCH 11/90] rugnux: default the stills integration box to r=6 (6,8,12) Serial-stills spots span a range of crystal orientations within a single shot, so they land wider on the detector than the r1=4 monochromatic-rotation default assumes, and the r1=4 signal box truncates the (bandwidth/mosaicity-broadened) spot wings. The profile-fit integrator can use a larger box for free - its profile weighting drives far-out background pixels to ~zero weight, so the estimator variance saturates rather than growing (a plain box-sum, which sums every pixel at weight 1, degrades). Gated on stills (!rotation_indexing); rotation stays at 4,6,10 (a larger box there only adds background - measured wash-to-negative on the XDS rotation battery, with an SG regression); an explicit --integration-radius wins. Validated: OCP R-free 0.406->0.400, LOV 0.354->0.336; lyso8 1.60->1.57 A; lyso5 (1% DMM) 1.80->1.72 A with +13% reflections at equal CC1/2 - the broadened case benefits most. Co-Authored-By: Claude Opus 4.8 (1M context) --- rugnux/rugnux_cli.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index e9f86dc1..40621736 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1310,6 +1310,16 @@ int main(int argc, char **argv) { bis.R1(r1).R2(r2).R3(r3); experiment.ImportBraggIntegrationSettings(bis); logger.Info("Integration radii set to r1={:.1f} r2={:.1f} r3={:.1f}", r1, r2, r3); + } else if (!rotation_indexing) { + // Stills spots span a range of crystal orientations captured in a single shot, so they land + // wider on the detector than the r1=4 monochromatic-rotation default assumes. A larger signal + // box lets the profile-fit integrator capture the whole spot while its profile weighting keeps + // the extra background from adding noise (a plain box-sum degrades with it). Measured R-free + // gains on serial stills (OCP, LOV). An explicit --integration-radius always wins. + BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); + bis.R1(6.0f).R2(8.0f).R3(12.0f); + experiment.ImportBraggIntegrationSettings(bis); + logger.Info("Stills integration radii default to r1=6.0 r2=8.0 r3=12.0 (override with --integration-radius)"); } if (integrator_mode) { -- 2.54.0 From fb399480e5c93808aa93cdb3ff9ee6fab888d735 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 12:03:21 +0200 Subject: [PATCH 12/90] rugnux: add opt-in stills partiality (--still-partiality) Serial stills currently treat every reflection as a full (partiality hardcoded to 1). Add an opt-in Gaussian excitation-error partiality set at prediction time (CPU + CUDA): p = exp(-dist_ewald^2 / (2*sigma^2)), sigma^2 = profile_radius^2 + (bandwidth_sigma*|recip_z|)^2, with sigma = the per-image profile radius (ewald_dist_cutoff/2), so an edge-of-acceptance reflection keeps p ~ exp(-2). Off by default; the merge weight (~p^2) then down-weights far-from-Ewald partials instead of trusting them as fulls. Validated: helps medium/strong stills (LOV R-free 0.336->0.329, lyso8 0.433->0.410, lowers the systematic error-model b in both) but HURTS weak OCP (dividing by a small, uncertain p amplifies orientation error -> high-res noise, resolution collapse), so it is left opt-in. A static forward p explains only ~6-10% of the partiality scatter; the full win needs per-image post-refinement (future work), for which this is the prediction-side groundwork. Co-Authored-By: Claude Opus 4.8 (1M context) --- common/BraggIntegrationSettings.cpp | 9 +++++++++ common/BraggIntegrationSettings.h | 3 +++ image_analysis/IndexAndRefine.cpp | 6 +++++- .../bragg_prediction/BraggPrediction.cpp | 9 ++++++++- .../bragg_prediction/BraggPrediction.h | 10 ++++++++-- .../bragg_prediction/BraggPredictionGPU.cu | 17 ++++++++++++++--- .../bragg_prediction/BraggPredictionGPU.h | 2 ++ rugnux/rugnux_cli.cpp | 14 ++++++++++++++ 8 files changed, 63 insertions(+), 7 deletions(-) diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index d69a0028..3436852f 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -96,3 +96,12 @@ float BraggIntegrationSettings::GetDMinLimit_A() const { float BraggIntegrationSettings::GetMinimumSigmaInRegardsToI() const { return minimum_sigma_in_regards_to_i; } + +BraggIntegrationSettings &BraggIntegrationSettings::StillPartiality(bool input) { + still_partiality = input; + return *this; +} + +bool BraggIntegrationSettings::GetStillPartiality() const { + return still_partiality; +} diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index bd53b7db..e0598756 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -20,6 +20,7 @@ class BraggIntegrationSettings { float d_min_limit_A = 1.0; std::optional fixed_profile_radius; float minimum_sigma_in_regards_to_i = 0.02; + bool still_partiality = false; // experimental stills excitation-error partiality (rugnux --still-partiality) public: BraggIntegrationSettings& R1(float input); @@ -28,6 +29,7 @@ public: BraggIntegrationSettings& DMinLimit_A(float input); BraggIntegrationSettings& FixedProfileRadius_recipA(std::optional input); BraggIntegrationSettings& Integrator(IntegratorMode input); + BraggIntegrationSettings& StillPartiality(bool input); [[nodiscard]] IntegratorMode GetIntegrator() const; @@ -38,4 +40,5 @@ public: [[nodiscard]] float GetDMinLimit_A() const; [[nodiscard]] float GetMinimumSigmaInRegardsToI() const; + [[nodiscard]] bool GetStillPartiality() const; }; diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index b0de29fa..8d542eac 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -324,7 +324,11 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, .wedge_deg = std::fabs(wedge_deg), .mosaicity_deg = std::fabs(mos_deg), // FWHM -> sigma; 0 when monochromatic, leaving the prediction unchanged. - .bandwidth_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f + .bandwidth_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f, + // Experimental stills partiality (off by default): sigma = ewald_dist_cutoff/2 = the per-image + // profile radius, so a reflection at the acceptance edge (dist_ewald ~ 2*sigma) keeps p ~ exp(-2). + .still_partiality = experiment.GetBraggIntegrationSettings().GetStillPartiality(), + .profile_radius_recipA = ewald_dist_cutoff * 0.5f }; // Predict, then integrate with the selected integrator (box-sum or profile-fit). diff --git a/image_analysis/bragg_prediction/BraggPrediction.cpp b/image_analysis/bragg_prediction/BraggPrediction.cpp index 701da2b9..e19d0249 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.cpp +++ b/image_analysis/bragg_prediction/BraggPrediction.cpp @@ -156,6 +156,13 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal continue; float d = 1.0f / sqrtf(recip_sq); + float partiality = 1.0f; + if (settings.still_partiality && settings.profile_radius_recipA > 0.0f) { + const float sig_bw = settings.bandwidth_sigma * std::fabs(recip_z); + const float sigma2 = settings.profile_radius_recipA * settings.profile_radius_recipA + + sig_bw * sig_bw; + partiality = std::exp(-0.5f * dist_ewald_sphere * dist_ewald_sphere / sigma2); + } reflections[i] = Reflection{ .h = h, .k = k, @@ -168,7 +175,7 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal .d = d, .dist_ewald = dist_ewald_sphere, .rlp = 1.0, - .partiality = 1.0, + .partiality = partiality, .zeta = 1.0, .image_scale_corr = 1.0 }; diff --git a/image_analysis/bragg_prediction/BraggPrediction.h b/image_analysis/bragg_prediction/BraggPrediction.h index a2a5b52a..5406d23e 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.h +++ b/image_analysis/bragg_prediction/BraggPrediction.h @@ -21,9 +21,15 @@ struct BraggPredictionSettings { // Relative X-ray bandwidth Δλ/λ expressed as a Gaussian sigma (0 = monochromatic). // When > 0 the Ewald-shell acceptance is thickened radially per reflection by // σ_bw = |recip_z|·bandwidth_sigma (= bλ/2d²), so the 1/d² pink-beam smear no - // longer clips high-resolution reflections. Must stay the last member so existing - // designated initializers remain valid. + // longer clips high-resolution reflections. float bandwidth_sigma = 0.0f; + // Experimental stills partiality (rugnux --still-partiality). When still_partiality is set and + // profile_radius_recipA > 0, each reflection gets a Gaussian excitation-error partiality + // p = exp(-dist_ewald^2 / (2*sigma^2)), sigma^2 = profile_radius_recipA^2 + (bandwidth_sigma*|recip_z|)^2, + // instead of the fixed 1.0 (full). Off by default. Keep these two as the trailing members so the + // existing designated initializers (which stop at bandwidth_sigma) remain valid. + bool still_partiality = false; + float profile_radius_recipA = 0.0f; }; class BraggPrediction { diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.cu b/image_analysis/bragg_prediction/BraggPredictionGPU.cu index 615eb794..8e172867 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.cu @@ -127,7 +127,13 @@ namespace { out.d = 1.0f / sqrtf(recip_sq); out.dist_ewald = dist_ewald; out.rlp = 1.0f; - out.partiality = 1.0f; + float partiality = 1.0f; + if (C.still_partiality && C.profile_radius_recipA > 0.0f) { + const float sig_bw = C.bandwidth_sigma * fabsf(recip_z); + const float sigma2 = C.profile_radius_recipA * C.profile_radius_recipA + sig_bw * sig_bw; + partiality = expf(-0.5f * dist_ewald * dist_ewald / sigma2); + } + out.partiality = partiality; out.zeta = 1.0f; out.image_scale_corr = 1.0f; return true; @@ -158,7 +164,9 @@ namespace { float high_res_A, float ewald_dist_cutoff, char centering, - float bandwidth_sigma) { + float bandwidth_sigma, + bool still_partiality, + float profile_radius_recipA) { KernelConsts kc{}; auto geom = experiment.GetDiffractionGeometry(); kc.det_width_pxl = static_cast(experiment.GetXPixelsNum()); @@ -171,6 +179,8 @@ namespace { kc.one_over_wavelength = 1.0f / geom.GetWavelength_A(); kc.ewald_cutoff = ewald_dist_cutoff; kc.bandwidth_sigma = bandwidth_sigma; + kc.still_partiality = still_partiality; + kc.profile_radius_recipA = profile_radius_recipA; kc.Astar = lattice.Astar(); kc.Bstar = lattice.Bstar(); kc.Cstar = lattice.Cstar(); @@ -193,7 +203,8 @@ int BraggPredictionGPU::Calc(const DiffractionExperiment &experiment, const BraggPredictionSettings &settings) { // Build constants on host KernelConsts hK = BuildKernelConsts(experiment, lattice, settings.high_res_A, settings.ewald_dist_cutoff, - settings.centering, settings.bandwidth_sigma); + settings.centering, settings.bandwidth_sigma, + settings.still_partiality, settings.profile_radius_recipA); cudaMemcpyAsync(dK, &hK, sizeof(KernelConsts), cudaMemcpyHostToDevice, stream); cudaMemsetAsync(d_count, 0, sizeof(int), stream); diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.h b/image_analysis/bragg_prediction/BraggPredictionGPU.h index ba0fc1f5..a7d2032d 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.h +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.h @@ -19,6 +19,8 @@ struct KernelConsts { float one_over_dmax_sq; float ewald_cutoff; float bandwidth_sigma; // relative Δλ/λ (sigma); 0 = monochromatic + bool still_partiality; // experimental stills excitation-error partiality (off => p = 1, fulls) + float profile_radius_recipA; // Gaussian sigma [1/A] for the stills partiality Coord Astar, Bstar, Cstar, S0; float rot[9]; char centering; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 40621736..42446fec 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -109,6 +109,7 @@ void print_usage() { std::cout << " --bandwidth Relative X-ray bandwidth FWHM (e.g. 0.01 for 1% DMM); default from file or 0" << std::endl; std::cout << " --integration-radius Signal-box radius r1, or r1,r2,r3 (px). One value => r2=r1+2, r3=r1+4" << std::endl; std::cout << " --integrator Spot integrator boxsum|gaussian|empirical (default: gaussian profile-fit; boxsum is the classical fallback)" << std::endl; + std::cout << " --still-partiality Experimental: weight stills reflections by a Gaussian excitation-error partiality exp(-dist_ewald^2/2sigma^2) instead of treating each as a full" << std::endl; std::cout << " -q, --azim-q-spacing Azimuthal-integration Q bin spacing (1/A) (default: 0.01)" << std::endl; std::cout << " --azim-min-q Azimuthal-integration minimum Q (1/A)" << std::endl; std::cout << " --azim-max-q Azimuthal-integration maximum Q (1/A)" << std::endl; @@ -153,6 +154,7 @@ enum { OPT_MODEL, OPT_DUMP_OBSERVATIONS, OPT_INTEGRATOR, + OPT_STILL_PARTIALITY, OPT_SCALE_FULLS, OPT_CAPTURE_UNCERTAINTY, OPT_MIN_CAPTURED_FRACTION, @@ -246,6 +248,7 @@ static option long_options[] = { {"bandwidth", required_argument, nullptr, OPT_BANDWIDTH}, {"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS}, {"integrator", required_argument, nullptr, OPT_INTEGRATOR}, + {"still-partiality", no_argument, nullptr, OPT_STILL_PARTIALITY}, {"detect-ice-rings", optional_argument, nullptr, OPT_DETECT_ICE_RINGS}, {"reject-outliers", required_argument, nullptr, OPT_REJECT_OUTLIERS}, {"reject-delta-cchalf", required_argument, nullptr, OPT_REJECT_DELTA_CCHALF}, @@ -489,6 +492,7 @@ int main(int argc, char **argv) { std::optional report_shell_count; // --resolution-shells std::optional integration_radius_arg; // "r1" or "r1,r2,r3" std::optional integrator_mode; // --integrator boxsum|gaussian|empirical + bool still_partiality_flag = false; // --still-partiality (experimental stills partiality) std::optional outlier_reject_nsigma; // merge per-observation outlier rejection std::optional delta_cchalf_nsigma; // per-crystal CC1/2-delta rejection @@ -726,6 +730,9 @@ int main(int argc, char **argv) { else if (strcmp(optarg, "empirical") == 0) integrator_mode = IntegratorMode::ProfileEmpirical; else { logger.Error("--integrator expects boxsum|gaussian|empirical"); return 1; } break; + case OPT_STILL_PARTIALITY: + still_partiality_flag = true; + break; case OPT_REJECT_OUTLIERS: outlier_reject_nsigma = parse_double_arg(optarg, "--reject-outliers", logger); break; @@ -1331,6 +1338,13 @@ int main(int argc, char **argv) { : "profile (empirical)"); } + if (still_partiality_flag) { + BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); + bis.StillPartiality(true); + experiment.ImportBraggIntegrationSettings(bis); + logger.Info("Stills partiality enabled (experimental Gaussian excitation-error weighting)"); + } + SpotFindingSettings spot_settings; spot_settings.enable = true; spot_settings.indexing = true; -- 2.54.0 From 4a5028d98d5bd94bd3d0b88be5987b604551b7c0 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 12:04:32 +0200 Subject: [PATCH 13/90] rugnux: add -r multi geometry refinement (per-image best of three) Per-image geometry refinement is a tradeoff that flips per dataset: beam+cell refinement extends OCP merge quality but DIVERGES on sparse-spot stills (e.g. KR2's ~20 spots around a 232 A axis), pushing good lattices out of tolerance so they fail the acceptance floor. `-r multi` runs all three (none / orientation / beam_and_lattice) on a copy per image and keeps whichever indexes the most spots (scorer = fractional-Miller-within-tolerance count, mirroring AnalyzeIndexing); ties prefer less refinement, to avoid overfitting the sparse list. Validated: KR2 index 7.08% (-r beam_and_lattice default) -> 10.05% (matches the -r none best), while OCP R-free stays ~equal to beam_and_lattice. New GeomRefinementAlgorithmEnum::Multi handled in the CLI/HDF5/command-line echoes; the API convert maps it to beam_and_lattice (offline only). Co-Authored-By: Claude Opus 4.8 (1M context) --- broker/OpenAPIConvert.cpp | 4 +++ common/IndexingSettings.h | 2 +- image_analysis/IndexAndRefine.cpp | 50 +++++++++++++++++++++++++++++++ rugnux/RugnuxCommandLine.cpp | 1 + rugnux/rugnux_cli.cpp | 4 ++- writer/HDF5NXmx.cpp | 3 ++ 6 files changed, 62 insertions(+), 2 deletions(-) diff --git a/broker/OpenAPIConvert.cpp b/broker/OpenAPIConvert.cpp index be884867..61b02aba 100644 --- a/broker/OpenAPIConvert.cpp +++ b/broker/OpenAPIConvert.cpp @@ -1023,6 +1023,10 @@ org::openapitools::server::model::Indexing_settings Convert(const IndexingSettin case GeomRefinementAlgorithmEnum::BeamCenter: refinement.setValue(org::openapitools::server::model::Geom_refinement_algorithm::eGeom_refinement_algorithm::BEAMCENTER); break; + case GeomRefinementAlgorithmEnum::Multi: + // No dedicated API value (offline rugnux only); report it as beam+lattice. + refinement.setValue(org::openapitools::server::model::Geom_refinement_algorithm::eGeom_refinement_algorithm::BEAMCENTER); + break; } ret.setGeomRefinementAlgorithm(refinement); diff --git a/common/IndexingSettings.h b/common/IndexingSettings.h index 6e900440..7b92297d 100644 --- a/common/IndexingSettings.h +++ b/common/IndexingSettings.h @@ -6,7 +6,7 @@ #include enum class IndexingAlgorithmEnum {FFBIDX, FFT, FFTW, Auto, None}; -enum class GeomRefinementAlgorithmEnum {None, OrientationOnly, BeamCenter}; +enum class GeomRefinementAlgorithmEnum {None, OrientationOnly, BeamCenter, Multi}; class IndexingSettings { IndexingAlgorithmEnum algorithm; diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 8d542eac..b73a97ab 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -176,6 +176,23 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data return outcome; } +namespace { +// Count spots whose fractional Miller index falls within the indexing tolerance of an integer for a +// given lattice + geometry - the "how well does this model explain the spots" score used by -r multi. +int CountIndexedSpots(const DiffractionGeometry &geom, const CrystalLattice &latt, + const std::vector &spots, float tol_sq) { + const Coord a = latt.Vec0(), b = latt.Vec1(), c = latt.Vec2(); + int n = 0; + for (const auto &s : spots) { + const Coord recip = s.ReciprocalCoord(geom); + const float hf = recip * a, kf = recip * b, lf = recip * c; + const float dh = hf - std::round(hf), dk = kf - std::round(kf), dl = lf - std::round(lf); + if (dh * dh + dk * dk + dl * dl < tol_sq) ++n; + } + return n; +} +} // namespace + void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::IndexingOutcome &outcome) { if (!outcome.lattice_candidate) return; @@ -213,6 +230,39 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In outcome.beam_center_updated = true; } break; + case GeomRefinementAlgorithmEnum::Multi: { + // Try all three refinements per image and keep whichever indexes the most spots. Beam+cell + // refinement helps some stills but diverges on sparse spot lists (few spots, long axes), + // where it pushes a good lattice out of tolerance; scoring by indexed-spot count lets each + // image fall back to orientation-only or no refinement when refinement would hurt. Ties + // prefer less refinement (strict >, not >=) to avoid overfitting. + const float tol = experiment.GetIndexingSettings().GetTolerance(); + const float tol_sq = tol * tol; + + XtalOptimizerData d_none = data; + XtalOptimizerData d_orient = data; + XtalOptimizerRotationOnly(d_orient, msg.spots, 0.2); + XtalOptimizerRotationOnly(d_orient, msg.spots, 0.1); + XtalOptimizerRotationOnly(d_orient, msg.spots, 0.05); + XtalOptimizerData d_beam = data; + const bool beam_ok = XtalOptimizer(d_beam, {msg.spots}); + + const int s_none = CountIndexedSpots(d_none.geom, d_none.latt, msg.spots, tol_sq); + const int s_orient = CountIndexedSpots(d_orient.geom, d_orient.latt, msg.spots, tol_sq); + const int s_beam = beam_ok ? CountIndexedSpots(d_beam.geom, d_beam.latt, msg.spots, tol_sq) : -1; + + if (s_beam > s_none && s_beam > s_orient) { + data = d_beam; + outcome.experiment.BeamX_pxl(data.geom.GetBeamX_pxl()) + .BeamY_pxl(data.geom.GetBeamY_pxl()); + outcome.beam_center_updated = true; + } else if (s_orient > s_none) { + data = d_orient; + } else { + data = d_none; + } + break; + } } outcome.lattice_candidate = data.latt; diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 3ea8f241..b8390bde 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -36,6 +36,7 @@ namespace { switch (r) { case GeomRefinementAlgorithmEnum::None: return "none"; case GeomRefinementAlgorithmEnum::OrientationOnly: return "orientation"; + case GeomRefinementAlgorithmEnum::Multi: return "multi"; case GeomRefinementAlgorithmEnum::BeamCenter: default: return "beam_and_lattice"; } diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 42446fec..fa417766 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -75,7 +75,7 @@ void print_usage() { std::cout << " -X, --indexing-algorithm Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)" << std::endl; std::cout << " -S, --space-group Space group number - used for both indexing and scaling" << std::endl; std::cout << " -C, --unit-cell Fix reference unit cell: \"a,b,c,alpha,beta,gamma\"" << std::endl; - std::cout << " -r, --refine Geometry refinement algorithm (none|orientation|beam_and_lattice)" << std::endl; + std::cout << " -r, --refine Geometry refinement algorithm (none|orientation|beam_and_lattice|multi); multi tries all three per image and keeps whichever indexes the most spots" << std::endl; std::cout << std::endl; std::cout << " Scaling and merging (on by default)" << std::endl; @@ -604,6 +604,8 @@ int main(int argc, char **argv) { refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter; else if (alg == "orientation") refinement_algorithm = GeomRefinementAlgorithmEnum::OrientationOnly; + else if (alg == "multi") + refinement_algorithm = GeomRefinementAlgorithmEnum::Multi; else { logger.Error("Invalid geom refinement algorithm: {}", alg); print_usage(); diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index 0497e96e..44ecfeea 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -455,6 +455,9 @@ void NXmx::MX(const StartMessage &start) { case GeomRefinementAlgorithmEnum::BeamCenter: hdf5_file->SaveScalar("/entry/MX/geom_refinement_algorithm", "beam_center"); break; + case GeomRefinementAlgorithmEnum::Multi: + hdf5_file->SaveScalar("/entry/MX/geom_refinement_algorithm", "multi"); + break; default: break; } -- 2.54.0 From 9527fa6e2c9698ec48e13769bc6ff3f6ca05bb94 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 12:05:08 +0200 Subject: [PATCH 14/90] rugnux: scale self-referenced stills in a single pass (fix weak-data collapse) The stills self-scaling loop rebuilt the merge reference FROM the just-scaled reflections every iteration, so on weak data each pass re-fit its own noise and the merged CC1/2 collapsed (PfluDING combine: 0.3% at the default 3 iters). It is provably pointless too: ScaleOnTheFly::Scale re-solves the per-image G from raw I and never reads the prior correction, so against a fixed reference passes 2..N are bit-identical - the single G is the exact one-pass solution and iteration only ever does the harmful reference rebuild (N=1 is best on every dataset, including strong lyso 82% vs 76% at N=3; leave-one-out and external anchors were tested and do not beat it). Do one self-reference pass; rotation (RotationScaleMerge) and the external-reference branch are untouched. Validated: PfluDING-combined CC1/2 0.3% -> 42.0% by default (pinned 2.5 A); strong data unchanged. Co-Authored-By: Claude Opus 4.8 (1M context) --- rugnux/Rugnux.cpp | 12 ++++++------ rugnux/rugnux_cli.cpp | 7 ++++++- 2 files changed, 12 insertions(+), 7 deletions(-) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 5ceb6a80..242b60ee 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -594,12 +594,12 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { // running merge with ScaleOnTheFly (fixed partiality), then merge directly. With an external // reference each image is already scaled against it during the per-image pass, so skip. if (config_.reference_data.empty()) { - for (int i = 0; i < config_.scaling_iter; i++) { - phase("Scaling images (" + label + ", iteration " + std::to_string(i + 1) + "/" - + std::to_string(config_.scaling_iter) + ")"); - auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome(), false); - indexer->ScaleAllImages(merge_result); - } + // One pass: the per-image scale G is the exact one-pass solution, so iterating only + // rebuilds the reference from the freshly-scaled (noisy) data - degrading weak stills + // instead of converging (measured CC1/2 collapse at the default 3 iters on weak data). + phase("Scaling images (" + label + ")"); + auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome(), false); + indexer->ScaleAllImages(merge_result); } const std::vector &merge_input = indexer->GetIntegrationOutcome(); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index fa417766..78cde433 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1035,7 +1035,12 @@ int main(int argc, char **argv) { merged_statistics = std::move(r.statistics); error_model_isa = r.isa; } else { - for (int i = 0; i < scaling_iter; i++) { + // The per-image scale G is the exact one-pass solution (Scale re-solves from raw I and never + // reads the prior correction), so with a self-rebuilt reference iterating only re-fits the + // freshly-scaled noise and degrades weak stills (measured CC1/2 collapse at the default 3 + // iters). Scale the self-reference merge once; a fixed external reference keeps scaling_iter. + const int n_self = reference_data.empty() ? 1 : scaling_iter; + for (int i = 0; i < n_self; i++) { if (reference_data.empty()) ScaleOnTheFly(experiment, MergeAll(experiment, reflections)).Scale(reflections, nthreads); else -- 2.54.0 From 95af137146cdba8dc58cc99f9ce462b6c6f8d07c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 12:18:19 +0200 Subject: [PATCH 15/90] rugnux: the external reference MTZ is never a scale anchor (stills + rotation) An external reference MTZ is a poor per-image scale reference: it is a different crystal/dataset, so scaling against it injects cross-dataset systematics rather than removing per-image scale error (measured: reference-scaled R-free 0.410 vs self-scaled 0.388 on OCP; 0.378 vs 0.350 on LOV). Its real value is fixing the cell/space group, breaking the merohedral indexing ambiguity, reporting CCref and providing the R-free test set - none of which need it to be a scale anchor. Make that consistent across both workflows: - Stills: the per-image pass (ScaleImage) now uses the reference ONLY to break the indexing ambiguity (once, for good), not to fit G. Scaling self-references at the post-measurement merge for every run, whether or not a reference is given (Rugnux + rugnux_cli --scale). - Rotation: ScaleImage no longer scales against the reference either (it was dead code - RotationScaleMerge recomputes the per-frame scale, so the reference never actually moved the result). Rotation resolves the ambiguity globally (ChooseReindex / ReferenceIntensityCC) and self-scales in RotationScaleMerge. - A reference-as-scale mode, if ever wanted, would be a dedicated later step, not this per-image pass. Validated: OCP CCref 56.5%->64.6%, R-free 0.400->0.393 (honest 2.36 A cutoff instead of a reference- propped 1.5 A); LOV CCref 88.3%->90.1%, R-free 0.336->0.331 with more reflections; rotation lyso_ref with the reference gives ISa 17.8 == de-novo 15.6 (same 15975 reflections, CC1/2 99.8%) - the reference demonstrably does not touch rotation scaling, only cell/SG + ambiguity + CCref. Co-Authored-By: Claude Opus 4.8 (1M context) --- image_analysis/IndexAndRefine.cpp | 20 ++++++++++++-------- image_analysis/scale_merge/ScaleOnTheFly.h | 6 ++++-- rugnux/Rugnux.cpp | 14 ++++++++------ rugnux/rugnux_cli.cpp | 17 ++++++----------- 4 files changed, 30 insertions(+), 27 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index b73a97ab..f2178a4b 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -495,15 +495,19 @@ void IndexAndRefine::ScaleImage(DataMessage &msg, IntegrationOutcome& outcome) { if (!scaling_engine) return; - auto scaling_start_time = std::chrono::steady_clock::now(); - scaling_engine->Scale(outcome); - auto scaling_end_time = std::chrono::steady_clock::now(); + // The external reference fixes the cell/space group, breaks the indexing ambiguity and reports CCref, + // but is NEVER a scale anchor: scaling an image against a foreign dataset injects cross-dataset + // systematics and is a worse reference than the data's own merge, so scaling self-references at the + // post-measurement merge for both workflows. Rotation resolves the ambiguity globally and self-scales + // in RotationScaleMerge (ChooseReindex / ReferenceIntensityCC), so there is nothing to do per image. + // Stills resolve the merohedral ambiguity per image here (each crystal indexes in a random hand; pick + // the hand best-correlated with the reference, once and for good). + if (experiment.IsRotationIndexing()) + return; - scale_cc[msg.number] = outcome.image_scale_cc.value_or(NAN); - msg.image_scale_b_factor = outcome.image_scale_b_factor_Ang2; - msg.image_scale_factor = outcome.image_scale_g; - msg.image_scale_mosaicity = outcome.mosaicity_deg; - msg.image_scale_cc = outcome.image_scale_cc; + auto scaling_start_time = std::chrono::steady_clock::now(); + scaling_engine->ResolveIndexingAmbiguity(outcome.reflections); + auto scaling_end_time = std::chrono::steady_clock::now(); msg.image_scale_time_s = std::chrono::duration(scaling_end_time - scaling_start_time).count(); } diff --git a/image_analysis/scale_merge/ScaleOnTheFly.h b/image_analysis/scale_merge/ScaleOnTheFly.h index 6166b4f1..63d9e587 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.h +++ b/image_analysis/scale_merge/ScaleOnTheFly.h @@ -39,8 +39,6 @@ class ScaleOnTheFly { bool Accept(const Reflection &r) const; [[nodiscard]] std::pair CalculateGlobalCC(const std::vector &reflections) const; - // Reindex this image's reflections into the ambiguity hand that best correlates with the reference. - void ResolveIndexingAmbiguity(std::vector &reflections) const; public: // resolve_ambiguity: when the reference is an external, correctly-handed dataset (not the data's own // running merge), pick per image the reindexing that agrees best with it - serial stills index each @@ -48,6 +46,10 @@ public: ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref, bool resolve_ambiguity = false); + // Reindex this image's reflections into the merohedral-ambiguity hand that best correlates with the + // external reference (stills; done once per image at integration time, decoupled from scaling). + void ResolveIndexingAmbiguity(std::vector &reflections) const; + void Scale(IntegrationOutcome &integration_outcome) const; void Scale(std::vector &integration_outcome, size_t nthreads = 0) const; }; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 242b60ee..210ff371 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -591,12 +591,14 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { return ScaleMergeResult{std::move(r.merged), std::move(r.statistics)}; } // Stills (rotation goes through RotationScaleMerge above): self-scale each image against the - // running merge with ScaleOnTheFly (fixed partiality), then merge directly. With an external - // reference each image is already scaled against it during the per-image pass, so skip. - if (config_.reference_data.empty()) { - // One pass: the per-image scale G is the exact one-pass solution, so iterating only - // rebuilds the reference from the freshly-scaled (noisy) data - degrading weak stills - // instead of converging (measured CC1/2 collapse at the default 3 iters on weak data). + // running merge with ScaleOnTheFly (fixed partiality), then merge directly. This runs even + // with an external reference: the reference is used only to break the per-image indexing + // ambiguity and to report CCref / inherit the R-free set, NOT as a scale anchor (scaling each + // image against a foreign dataset injects cross-dataset systematics and is a worse reference + // than the data's own merge). One pass: the per-image scale G is the exact one-pass solution, + // so iterating would only rebuild the reference from the freshly-scaled (noisy) data and + // degrade weak stills (measured CC1/2 collapse at the default 3 iters). + { phase("Scaling images (" + label + ")"); auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome(), false); indexer->ScaleAllImages(merge_result); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 78cde433..371cc39c 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1035,17 +1035,12 @@ int main(int argc, char **argv) { merged_statistics = std::move(r.statistics); error_model_isa = r.isa; } else { - // The per-image scale G is the exact one-pass solution (Scale re-solves from raw I and never - // reads the prior correction), so with a self-rebuilt reference iterating only re-fits the - // freshly-scaled noise and degrades weak stills (measured CC1/2 collapse at the default 3 - // iters). Scale the self-reference merge once; a fixed external reference keeps scaling_iter. - const int n_self = reference_data.empty() ? 1 : scaling_iter; - for (int i = 0; i < n_self; i++) { - if (reference_data.empty()) - ScaleOnTheFly(experiment, MergeAll(experiment, reflections)).Scale(reflections, nthreads); - else - ScaleOnTheFly(experiment, reference_data).Scale(reflections, nthreads); - } + // Scaling self-references: the reference MTZ (if any) fixes the cell/space group, reports + // CCref and provides the R-free test set, but is NOT a scale anchor - scaling each image + // against a foreign dataset injects cross-dataset systematics and is a worse reference than + // the data's own merge. The per-image scale G is the exact one-pass solution, so one pass + // (iterating a self-rebuilt reference only re-fits the freshly-scaled noise on weak stills). + ScaleOnTheFly(experiment, MergeAll(experiment, reflections)).Scale(reflections, nthreads); MergeOnTheFly merge_engine(experiment); merge_engine.ReferenceCell(experiment.GetUnitCell()); // Fit the (a, b) error model from symmetry-mate scatter before merging, exactly as the full -- 2.54.0 From 3f6462b30f96223b478694071cfb6c76c1ac1f91 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 15:18:58 +0200 Subject: [PATCH 16/90] rugnux: rotation indexing always uses FFT/FFTW, never FFBIDX The rotation two-pass accumulates a dense multi-frame reciprocal-space cloud (all first-pass frames de-rotated into one common frame) and needs a GLOBAL lattice finder - a 3D FFT whose peaks are the basis vectors. ffbidx is a single-still, known-cell ORIENTATION finder: it assumes the spots lie on one Ewald shell (one crystal, one orientation) and maximizes a near-integer-hkl count for a fixed cell. On the volume-filling rotation cloud that objective has no dominant global maximum, so ffbidx's coarse orientation scan locks onto an orientation valid for only a local wedge (~18%) -> garbage merge. GetIndexingAlgorithm resolved Auto->FFBIDX whenever a cell was known, even for rotation - and a reference MTZ (-z) or -C silently supplies the cell. So a routine rotation run with a reference for CCref/R-free flipped the indexer to the one solver that cannot handle its input (lyso_ref: ISa 17.8 with FFT -> 0.5 with FFBIDX, no warning). Force rotation to FFT (GPU) / FFTW (CPU) regardless of the requested algorithm or a known cell. Validated: lyso_ref -R with -z reference MTZ now indexes via FFT -> ISa 17.8 (was 0.5), 15975 refl. Co-Authored-By: Claude Opus 4.8 (1M context) --- common/DiffractionExperiment.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/common/DiffractionExperiment.cpp b/common/DiffractionExperiment.cpp index d7b0f661..81593874 100644 --- a/common/DiffractionExperiment.cpp +++ b/common/DiffractionExperiment.cpp @@ -1510,6 +1510,15 @@ DiffractionExperiment &DiffractionExperiment::IndexingAlgorithm(IndexingAlgorith } IndexingAlgorithmEnum DiffractionExperiment::GetIndexingAlgorithm() const { + // Rotation indexing accumulates a dense multi-frame reciprocal-space cloud that needs a GLOBAL + // lattice finder (3D FFT). ffbidx is a single-still, known-cell ORIENTATION finder: it cannot + // recover the global orientation from that cloud (it locks onto a local wedge -> garbage). So + // rotation always uses FFT (GPU) or FFTW (CPU), regardless of the requested algorithm or whether a + // cell is known - a reference MTZ (-z) or -C supplies a cell and used to silently flip Auto->FFBIDX, + // switching the rotation indexer to the one solver that cannot handle its input. + if (IsRotationIndexing()) + return get_gpu_count() == 0 ? IndexingAlgorithmEnum::FFTW : IndexingAlgorithmEnum::FFT; + auto cell = GetUnitCell().has_value(); switch (indexing.GetAlgorithm()) { -- 2.54.0 From 49963ab855e1e9bc02247fd6db9ec96716353019 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 15:35:40 +0200 Subject: [PATCH 17/90] rugnux: move indexing-ambiguity resolution out of ScaleOnTheFly Indexing-ambiguity detection/resolution is no longer part of scaling: it is a dedicated per-image step for stills (or a post-full-merge step for rotation). Pull the stills per-image resolver out of ScaleOnTheFly into a new ReindexAmbiguityResolver in ReindexAmbiguity.{h,cpp}, alongside the existing rotation free functions (ChooseReindex / ReferenceIntensityCC). Both workflows now share the operator generation and the best-op selection (new file-local PickBestReindex helper); ChooseReindex is rewired onto it. ScaleOnTheFly is now purely a scaling engine (no ambiguity_ops, no resolve_ambiguity flag, no ResolveIndexingAmbiguity). Pure refactor: behaviour is unchanged. IndexAndRefine now holds a ReindexAmbiguityResolver and calls Resolve() at the same call site. Co-Authored-By: Claude Opus 4.8 (1M context) --- image_analysis/IndexAndRefine.cpp | 11 +- image_analysis/IndexAndRefine.h | 3 +- .../scale_merge/ReindexAmbiguity.cpp | 124 ++++++++++++++++-- image_analysis/scale_merge/ReindexAmbiguity.h | 32 +++++ image_analysis/scale_merge/ScaleOnTheFly.cpp | 59 +-------- image_analysis/scale_merge/ScaleOnTheFly.h | 14 +- rugnux/Rugnux.cpp | 2 +- 7 files changed, 153 insertions(+), 92 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index f2178a4b..036e00ae 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -10,6 +10,7 @@ #include "indexing/FFTIndexer.h" #include "indexing/MultiLatticeSearch.h" #include "lattice_search/LatticeSearch.h" +#include "scale_merge/ReindexAmbiguity.h" #include "scale_merge/ScaleOnTheFly.h" IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer, bool retain_outcomes) @@ -485,14 +486,14 @@ std::optional IndexAndRefine::FinalizeRotationIndexing() } IndexAndRefine &IndexAndRefine::ReferenceIntensities(std::vector &reference) { - // An external reference is trusted to be in the correct hand, so let the per-image scaler also break - // the merohedral indexing ambiguity against it (serial stills index each crystal independently). - scaling_engine = std::make_unique(experiment, reference, /*resolve_ambiguity=*/true); + // An external reference is trusted to be in the correct hand, so use it to break the merohedral + // indexing ambiguity per image (serial stills index each crystal independently). + reindex_resolver = std::make_unique(experiment, reference); return *this; } void IndexAndRefine::ScaleImage(DataMessage &msg, IntegrationOutcome& outcome) { - if (!scaling_engine) + if (!reindex_resolver) return; // The external reference fixes the cell/space group, breaks the indexing ambiguity and reports CCref, @@ -506,7 +507,7 @@ void IndexAndRefine::ScaleImage(DataMessage &msg, IntegrationOutcome& outcome) { return; auto scaling_start_time = std::chrono::steady_clock::now(); - scaling_engine->ResolveIndexingAmbiguity(outcome.reflections); + reindex_resolver->Resolve(outcome.reflections); auto scaling_end_time = std::chrono::steady_clock::now(); msg.image_scale_time_s = std::chrono::duration(scaling_end_time - scaling_start_time).count(); } diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 3f912a8b..ed6b85d3 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -17,6 +17,7 @@ #include "lattice_search/LatticeSearch.h" #include "rotation_indexer/RotationIndexer.h" #include "rotation_indexer/RotationIndexerCounter.h" +#include "scale_merge/ReindexAmbiguity.h" #include "scale_merge/ScaleOnTheFly.h" #include "scale_merge/ScalingResult.h" #include "IntegrationOutcome.h" @@ -81,7 +82,7 @@ class IndexAndRefine { const BraggIntegrateFn &integrate, const IndexingOutcome &outcome); - std::unique_ptr scaling_engine; + std::unique_ptr reindex_resolver; void ScaleImage(DataMessage &msg, IntegrationOutcome& outcome); std::optional RotationAngle(int64_t image) const; // mid-exposure angle for the indexer diff --git a/image_analysis/scale_merge/ReindexAmbiguity.cpp b/image_analysis/scale_merge/ReindexAmbiguity.cpp index 929f641e..59e53558 100644 --- a/image_analysis/scale_merge/ReindexAmbiguity.cpp +++ b/image_analysis/scale_merge/ReindexAmbiguity.cpp @@ -10,6 +10,42 @@ #include "HKLKey.h" +namespace { + constexpr size_t MIN_REFLECTIONS = 20; + + double SafeInv(double x, double fallback) { + if (!std::isfinite(x) || x == 0.0) + return fallback; + return 1.0 / x; + } + + struct BestReindex { + gemmi::Op op = gemmi::Op::identity(); + bool is_identity = true; + double score = 0.0; + double identity_score = 0.0; + }; + + // Among identity (the baseline, scored `identity_score`) and `ops`, return the operator with the + // highest `score_op(op)`. A non-finite score never wins (too little overlap to decide). Shared by the + // stills per-image resolver and the rotation post-merge ChooseReindex. + BestReindex PickBestReindex(double identity_score, const std::vector &ops, + const std::function &score_op) { + BestReindex best; + best.identity_score = identity_score; + best.score = identity_score; + for (const auto &op : ops) { + const double sc = score_op(op); + if (std::isfinite(sc) && (!std::isfinite(best.score) || sc > best.score)) { + best.score = sc; + best.op = op; + best.is_identity = false; + } + } + return best; + } +} + std::vector ReindexAmbiguityOperators(const UnitCell &cell, int space_group_number, double max_obliquity_deg) { const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number); @@ -37,21 +73,18 @@ ReindexChoice ChooseReindex(const std::vector &merged, const UnitCell &cell, int space_group_number, const std::function &)> &score, double max_obliquity_deg) { - ReindexChoice best; - best.identity_score = score(merged); - best.score = best.identity_score; - const auto ops = ReindexAmbiguityOperators(cell, space_group_number, max_obliquity_deg); - best.n_candidates = 1 + static_cast(ops.size()); - for (const auto &op : ops) { - const double s = score(ReindexReflections(merged, op)); - if (s > best.score) { - best.score = s; - best.op = op; - best.is_identity = false; - } - } - return best; + const BestReindex best = PickBestReindex( + score(merged), ops, + [&](const gemmi::Op &op) { return score(ReindexReflections(merged, op)); }); + + ReindexChoice choice; + choice.op = best.op; + choice.is_identity = best.is_identity; + choice.score = best.score; + choice.identity_score = best.identity_score; + choice.n_candidates = 1 + static_cast(ops.size()); + return choice; } double ReferenceIntensityCC(const std::vector &merged, @@ -87,3 +120,66 @@ double ReferenceIntensityCC(const std::vector &merged, const double vy = n * syy - sy * sy; return (vx > 0 && vy > 0) ? cov / std::sqrt(vx * vy) : 0.0; } + +ReindexAmbiguityResolver::ReindexAmbiguityResolver(const DiffractionExperiment &x, + const std::vector &reference) + : s(x.GetScalingSettings()), + hkl_key_generator(s.GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1)) { + for (const auto &r : reference) + reference_data[hkl_key_generator(r)] = r.I; + if (x.GetUnitCell().has_value() && x.GetSpaceGroupNumber().has_value()) + ops = ReindexAmbiguityOperators(*x.GetUnitCell(), static_cast(*x.GetSpaceGroupNumber())); +} + +bool ReindexAmbiguityResolver::Accept(const Reflection &r) const { + if (r.on_ice_ring) // ice-contaminated intensity would bias the correlation; keep it out + return false; + return AcceptReflection(r, s.GetHighResolutionLimit_A()); +} + +double ReindexAmbiguityResolver::ReferenceCC(const std::vector &reflections, + const gemmi::Op &op) const { + double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0; + size_t n = 0; + for (const auto &r : reflections) { + if (!Accept(r) || r.partiality < s.GetMinPartiality()) + continue; + if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f || r.partiality <= 0.0f) + continue; + const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}}); + const auto it = reference_data.find(hkl_key_generator(h[0], h[1], h[2])); + if (it == reference_data.end()) + continue; + const double x = static_cast(r.I) * SafeInv(r.rlp, 1.0) / r.partiality; + const double y = it->second; + if (!std::isfinite(x) || !std::isfinite(y)) + continue; + sx += x; sy += y; sxx += x * x; syy += y * y; sxy += x * y; + ++n; + } + if (n < MIN_REFLECTIONS) + return NAN; + const double nd = static_cast(n); + const double cov = sxy - sx * sy / nd; + const double vx = sxx - sx * sx / nd; + const double vy = syy - sy * sy / nd; + return (vx > 0 && vy > 0) ? cov / std::sqrt(vx * vy) : NAN; +} + +// Serial stills index each crystal in one of the merohedrally-equivalent hands at random; pick, for this +// image alone, the reindexing whose intensities correlate best with the external reference and apply it. +void ReindexAmbiguityResolver::Resolve(std::vector &reflections) const { + if (ops.empty()) + return; + + const BestReindex best = PickBestReindex( + ReferenceCC(reflections, gemmi::Op::identity()), ops, + [&](const gemmi::Op &op) { return ReferenceCC(reflections, op); }); + if (best.is_identity) + return; + + for (auto &r : reflections) { + const gemmi::Op::Miller h = best.op.apply_to_hkl({{r.h, r.k, r.l}}); + r.h = h[0]; r.k = h[1]; r.l = h[2]; + } +} diff --git a/image_analysis/scale_merge/ReindexAmbiguity.h b/image_analysis/scale_merge/ReindexAmbiguity.h index 300276d3..93f44ce5 100644 --- a/image_analysis/scale_merge/ReindexAmbiguity.h +++ b/image_analysis/scale_merge/ReindexAmbiguity.h @@ -4,10 +4,13 @@ #pragma once #include +#include #include #include "gemmi/symmetry.hpp" +#include "HKLKey.h" +#include "../../common/DiffractionExperiment.h" #include "../../common/Reflection.h" #include "../../common/UnitCell.h" @@ -50,3 +53,32 @@ ReindexChoice ChooseReindex(const std::vector &merged, double ReferenceIntensityCC(const std::vector &merged, const std::vector &reference, int space_group_number); + +// Per-image resolver for serial stills. Each crystal is indexed independently in one of the +// merohedrally-equivalent hands at random, so the ambiguity has to be broken per image (unlike rotation +// data, which resolve it once for the whole merged set with ChooseReindex above). Holds the external +// reference intensities and the crystal's twin-law operators, and reindexes one image's reflections into +// the hand best-correlated with the reference. This runs in the per-image pipeline at integration time, +// decoupled from scaling. +class ReindexAmbiguityResolver { + const ScalingSettings s; + const HKLKeyGenerator hkl_key_generator; + std::map reference_data; + // Twin-law reindexing operators (identity excluded). Empty = no ambiguity (holohedral crystal) or the + // cell/space group is unknown, in which case Resolve() is a no-op. + std::vector ops; + + bool Accept(const Reflection &r) const; + // Pearson correlation of this image's partiality/Lorentz-corrected intensities against the reference, + // matched by ASU index under `op`. Scale-invariant (the unknown per-image scale cancels). NaN when + // too few reflections match. + double ReferenceCC(const std::vector &reflections, const gemmi::Op &op) const; + +public: + // The reference must be an external, correctly-handed dataset (not the data's own running merge). + ReindexAmbiguityResolver(const DiffractionExperiment &x, const std::vector &reference); + + // Reindex this image's reflections in place into the hand best-correlated with the reference (done + // once per image, for good). No-op when there is no ambiguity. + void Resolve(std::vector &reflections) const; +}; diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index a2b3c8f0..175b2aac 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -96,67 +96,13 @@ namespace { }; } -ScaleOnTheFly::ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref, - bool resolve_ambiguity) +ScaleOnTheFly::ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref) : s(x.GetScalingSettings()), hkl_key_generator(s.GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1)) { for (const auto &r: ref) { const auto key = hkl_key_generator(r); reference_data[key] = r.I; } - if (resolve_ambiguity && x.GetUnitCell().has_value() && x.GetSpaceGroupNumber().has_value()) - ambiguity_ops = ReindexAmbiguityOperators(*x.GetUnitCell(), - static_cast(*x.GetSpaceGroupNumber())); -} - -// Serial stills index each crystal in one of the merohedrally-equivalent hands at random; pick, for this -// image alone, the reindexing whose intensities correlate best with the external reference and apply it. -// The correlation uses the partiality/Lorentz-corrected intensity (scale-invariant, so the unknown -// per-image scale G cancels) matched to the reference by ASU index under the candidate hand. -void ScaleOnTheFly::ResolveIndexingAmbiguity(std::vector &reflections) const { - auto reference_cc = [&](const gemmi::Op &op) -> std::pair { - double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0; - size_t n = 0; - for (const auto &r: reflections) { - if (!Accept(r) || r.partiality < s.GetMinPartiality()) - continue; - if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f || r.partiality <= 0.0f) - continue; - const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}}); - const auto it = reference_data.find(hkl_key_generator(h[0], h[1], h[2])); - if (it == reference_data.end()) - continue; - const double x = static_cast(r.I) * SafeInv(r.rlp, 1.0) / r.partiality; - const double y = it->second; - if (!std::isfinite(x) || !std::isfinite(y)) - continue; - sx += x; sy += y; sxx += x * x; syy += y * y; sxy += x * y; - ++n; - } - if (n < MIN_REFLECTIONS) - return {NAN, n}; - const double nd = static_cast(n); - const double cov = sxy - sx * sy / nd; - const double vx = sxx - sx * sx / nd; - const double vy = syy - sy * sy / nd; - return {(vx > 0 && vy > 0) ? cov / std::sqrt(vx * vy) : NAN, n}; - }; - - const double identity_cc = reference_cc(gemmi::Op::identity()).first; - const gemmi::Op *best = nullptr; - double best_cc = identity_cc; - for (const auto &op: ambiguity_ops) { - const double cc = reference_cc(op).first; - if (std::isfinite(cc) && (!std::isfinite(best_cc) || cc > best_cc)) { - best_cc = cc; - best = &op; - } - } - if (best != nullptr) - for (auto &r: reflections) { - const gemmi::Op::Miller h = best->apply_to_hkl({{r.h, r.k, r.l}}); - r.h = h[0]; r.k = h[1]; r.l = h[2]; - } } bool ScaleOnTheFly::Accept(const Reflection &r) const { @@ -223,9 +169,6 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { if (integration_outcome.reflections.empty()) return; - if (!ambiguity_ops.empty()) - ResolveIndexingAmbiguity(integration_outcome.reflections); - auto start = std::chrono::steady_clock::now(); ScaleOnTheFlyResult result{ .B = 0.0, .G = 1.0 }; diff --git a/image_analysis/scale_merge/ScaleOnTheFly.h b/image_analysis/scale_merge/ScaleOnTheFly.h index 63d9e587..67bd5623 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.h +++ b/image_analysis/scale_merge/ScaleOnTheFly.h @@ -5,7 +5,6 @@ #include "HKLKey.h" #include "Merge.h" -#include "ReindexAmbiguity.h" #include "../../common/DiffractionExperiment.h" #include "../IntegrationOutcome.h" @@ -33,22 +32,11 @@ class ScaleOnTheFly { const ScalingSettings s; const HKLKeyGenerator hkl_key_generator; std::map reference_data; - // Twin-law reindexing operators for the merohedral indexing ambiguity (empty = no ambiguity, or - // resolution disabled). Populated only when the ref is a trusted external reference (see constructor). - std::vector ambiguity_ops; bool Accept(const Reflection &r) const; [[nodiscard]] std::pair CalculateGlobalCC(const std::vector &reflections) const; public: - // resolve_ambiguity: when the reference is an external, correctly-handed dataset (not the data's own - // running merge), pick per image the reindexing that agrees best with it - serial stills index each - // crystal independently, so the merohedral ambiguity must be broken per image, not globally. - ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref, - bool resolve_ambiguity = false); - - // Reindex this image's reflections into the merohedral-ambiguity hand that best correlates with the - // external reference (stills; done once per image at integration time, decoupled from scaling). - void ResolveIndexingAmbiguity(std::vector &reflections) const; + ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref); void Scale(IntegrationOutcome &integration_outcome) const; void Scale(std::vector &integration_outcome, size_t nthreads = 0) const; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 210ff371..72406ddb 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -825,7 +825,7 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { // the crystal has an indexing ambiguity (merohedral: lattice symmetry higher than the Laue group), // pick the reindexing that best correlates with the reference intensities and re-merge in it. The // twin-law reindex is metric-preserving, so only the hkl labels change (the cell is unchanged). - // Stills resolve the ambiguity per image in ScaleOnTheFly, not here. + // Stills resolve the ambiguity per image with ReindexAmbiguityResolver, not here. if (rsm && !config_.reference_data.empty() && result.consensus_cell && experiment_.GetSpaceGroupNumber().has_value()) { const int sg_num = static_cast(*experiment_.GetSpaceGroupNumber()); -- 2.54.0 From 2f5ed411e242983e280938c5934f4b0bc32244aa Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 18:17:58 +0200 Subject: [PATCH 18/90] rugnux: --refine-geometry stills two-pass (global bundle adjustment + re-index) Weak/jet serial stills are often geometry-limited: a few-px beam error or a mm-scale detector distance error (unreachable per-image) fails many frames, but is well-determined jointly from the strong frames. Mirror the rotation two-pass: an index-only first pass over a spread sample collects each indexed frame's spots + assigned HKL + orientation; the strongest ~N (default 200) feed one Ceres bundle adjustment; the refined geometry is applied and the main pass re-indexes + integrates + merges every frame from scratch. GeometryRefiner (reusing the extracted XtalResidual - the RecipToDetector geometry residual pulled out of XtalOptimizer, behaviour-preserving): one problem with SHARED beam(2)/distance(1)/cell-length(3) blocks + a PER-FRAME orientation(3) block, robust Cauchy loss, a cell-length regularizer anchoring the known cell to break the low-resolution distance<->cell-scale degeneracy, DENSE_SCHUR eliminating the per-frame orientations, and a 3-round HKL-reassignment / tolerance-tightening loop. Tilt is not refined (gauge-coupled, zero gain). Opt-in via --refine-geometry[=N]; stills only (rotation untouched). Validated: KR2 7.58% -> 21.85% (matches CrystFEL's 21.5%; a real ~1.4mm distance error + ~3px beam), OCP 2.92% -> 4.19% (~3px beam). OFF runs are bit-identical to baseline (XtalResidual extraction non-regressing; rotation lyso_ref de-novo ISa 17.3 unchanged). Co-Authored-By: Claude Opus 4.8 (1M context) --- image_analysis/geom_refinement/CMakeLists.txt | 3 + .../geom_refinement/GeometryRefiner.cpp | 311 ++++++++++++++++++ .../geom_refinement/GeometryRefiner.h | 59 ++++ .../geom_refinement/XtalOptimizer.cpp | 194 +---------- image_analysis/geom_refinement/XtalResidual.h | 214 ++++++++++++ rugnux/Rugnux.cpp | 147 +++++++++ rugnux/Rugnux.h | 13 + rugnux/rugnux_cli.cpp | 14 + 8 files changed, 762 insertions(+), 193 deletions(-) create mode 100644 image_analysis/geom_refinement/GeometryRefiner.cpp create mode 100644 image_analysis/geom_refinement/GeometryRefiner.h create mode 100644 image_analysis/geom_refinement/XtalResidual.h diff --git a/image_analysis/geom_refinement/CMakeLists.txt b/image_analysis/geom_refinement/CMakeLists.txt index 51df98fa..dcf9936f 100644 --- a/image_analysis/geom_refinement/CMakeLists.txt +++ b/image_analysis/geom_refinement/CMakeLists.txt @@ -6,6 +6,9 @@ ADD_LIBRARY(JFJochGeomRefinement STATIC AssignSpotsToRings.h XtalOptimizer.cpp XtalOptimizer.h + XtalResidual.h + GeometryRefiner.cpp + GeometryRefiner.h LatticeReduction.cpp LatticeReduction.h ) diff --git a/image_analysis/geom_refinement/GeometryRefiner.cpp b/image_analysis/geom_refinement/GeometryRefiner.cpp new file mode 100644 index 00000000..b10ca118 --- /dev/null +++ b/image_analysis/geom_refinement/GeometryRefiner.cpp @@ -0,0 +1,311 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "GeometryRefiner.h" + +#include +#include +#include +#include + +#include "../../common/JFJochMath.h" // PI +#include "LatticeReduction.h" +#include "XtalResidual.h" +#include "ceres/ceres.h" + +namespace { + +// Trigonal is described with the hexagonal cell + B, same as the per-image XtalOptimizer. +gemmi::CrystalSystem MapSystem(gemmi::CrystalSystem s) { + return (s == gemmi::CrystalSystem::Trigonal) ? gemmi::CrystalSystem::Hexagonal : s; +} + +// The cell angles (radians) that XtalResidual bakes into its B matrix for a given system. Only +// monoclinic (beta) and triclinic (all three) read p2; the rest use fixed symmetry angles. We keep the +// angles constant in v1, so these are set once from the input cell. +void InitCellAngles(gemmi::CrystalSystem system, const UnitCell &cell, double ang[3]) { + const double d2r = PI / 180.0; + switch (system) { + case gemmi::CrystalSystem::Hexagonal: + ang[0] = PI / 2.0; ang[1] = PI / 2.0; ang[2] = 2.0 * PI / 3.0; + break; + case gemmi::CrystalSystem::Monoclinic: + // XtalResidual reads p2[0] as beta. + ang[0] = cell.beta * d2r; ang[1] = PI / 2.0; ang[2] = PI / 2.0; + break; + case gemmi::CrystalSystem::Triclinic: + ang[0] = cell.alpha * d2r; ang[1] = cell.beta * d2r; ang[2] = cell.gamma * d2r; + break; + default: // orthorhombic / tetragonal / cubic + ang[0] = PI / 2.0; ang[1] = PI / 2.0; ang[2] = PI / 2.0; + break; + } +} + +// Effective (lengths, alpha,beta,gamma) that reproduce the lattice XtalResidual builds from (len, ang), +// so a frame's lattice can be reconstructed with AngleAxisAndCellToLattice for HKL re-assignment. +void EffectiveCell(gemmi::CrystalSystem system, const double len[3], const double ang[3], + double lengths[3], double &alpha, double &beta, double &gamma) { + switch (system) { + case gemmi::CrystalSystem::Hexagonal: + lengths[0] = len[0]; lengths[1] = len[0]; lengths[2] = len[2]; + alpha = PI / 2.0; beta = PI / 2.0; gamma = 2.0 * PI / 3.0; + break; + case gemmi::CrystalSystem::Tetragonal: + lengths[0] = len[0]; lengths[1] = len[0]; lengths[2] = len[2]; + alpha = PI / 2.0; beta = PI / 2.0; gamma = PI / 2.0; + break; + case gemmi::CrystalSystem::Cubic: + lengths[0] = len[0]; lengths[1] = len[0]; lengths[2] = len[0]; + alpha = PI / 2.0; beta = PI / 2.0; gamma = PI / 2.0; + break; + case gemmi::CrystalSystem::Monoclinic: + lengths[0] = len[0]; lengths[1] = len[1]; lengths[2] = len[2]; + alpha = PI / 2.0; beta = ang[0]; gamma = PI / 2.0; + break; + case gemmi::CrystalSystem::Triclinic: + lengths[0] = len[0]; lengths[1] = len[1]; lengths[2] = len[2]; + alpha = ang[0]; beta = ang[1]; gamma = ang[2]; + break; + default: // orthorhombic + lengths[0] = len[0]; lengths[1] = len[1]; lengths[2] = len[2]; + alpha = PI / 2.0; beta = PI / 2.0; gamma = PI / 2.0; + break; + } +} + +// Orientation seed (angle-axis) for a frame's lattice, using the per-system decomposition the per-image +// optimizer uses so it matches XtalResidual's B convention. +void OrientationSeed(gemmi::CrystalSystem system, const CrystalLattice &latt, double rod[3]) { + double scratch[3]; + if (system == gemmi::CrystalSystem::Hexagonal) { + LatticeToRodriguesAndLengths_Hex(latt, rod, scratch); + } else if (system == gemmi::CrystalSystem::Monoclinic) { + double beta; + LatticeToRodriguesLengthsBeta_Mono(latt, rod, scratch, beta); + } else { + LatticeToRodriguesAndLengths_GS(latt, rod, scratch); + } +} + +UnitCell BuildUnitCell(gemmi::CrystalSystem system, const double len[3], const double ang[3]) { + const double r2d = 180.0 / PI; + UnitCell uc{}; + switch (system) { + case gemmi::CrystalSystem::Hexagonal: + uc = {(float)len[0], (float)len[0], (float)len[2], 90.0f, 90.0f, 120.0f}; + break; + case gemmi::CrystalSystem::Tetragonal: + uc = {(float)len[0], (float)len[0], (float)len[2], 90.0f, 90.0f, 90.0f}; + break; + case gemmi::CrystalSystem::Cubic: + uc = {(float)len[0], (float)len[0], (float)len[0], 90.0f, 90.0f, 90.0f}; + break; + case gemmi::CrystalSystem::Monoclinic: + uc = {(float)len[0], (float)len[1], (float)len[2], 90.0f, (float)(ang[0] * r2d), 90.0f}; + break; + case gemmi::CrystalSystem::Triclinic: + uc = {(float)len[0], (float)len[1], (float)len[2], + (float)(ang[0] * r2d), (float)(ang[1] * r2d), (float)(ang[2] * r2d)}; + break; + default: // orthorhombic + uc = {(float)len[0], (float)len[1], (float)len[2], 90.0f, 90.0f, 90.0f}; + break; + } + return uc; +} + +// Anchors the cell lengths to the input cell. Penalises each length's deviation from its input value +// with a large weight, breaking the low-resolution distance<->cell-scale degeneracy (position ~ +// distance / cell_scale) so distance is determined by the data with the cell pinned. +struct CellLengthRegularizer { + double w0, w1, w2; + double l0, l1, l2; + template + bool operator()(const T *const len, T *r) const { + r[0] = T(w0) * (len[0] - T(l0)); + r[1] = T(w1) * (len[1] - T(l1)); + r[2] = T(w2) * (len[2] - T(l2)); + return true; + } +}; + +} // namespace + +GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_geom, + const UnitCell &input_cell, + const std::vector &frames, + const GeometryRefinerSettings &settings) { + GeometryRefinerResult result; + try { + const gemmi::CrystalSystem system = MapSystem(settings.crystal_system); + + // Shared geometry parameter blocks (refined). Detector tilt and the rotation axis are constant. + double beam[2] = {nominal_geom.GetBeamX_pxl(), nominal_geom.GetBeamY_pxl()}; + double distance_mm = nominal_geom.GetDetectorDistance_mm(); + double detector_rot[2] = {nominal_geom.GetPoniRot1_rad(), nominal_geom.GetPoniRot2_rad()}; + double rot_vec[3] = {1.0, 0.0, 0.0}; + const double rot3 = nominal_geom.GetPoniRot3_rad(); + const double lambda = nominal_geom.GetWavelength_A(); + const double pixel = nominal_geom.GetPixelSize_mm(); + + // Shared cell blocks. Lengths refined (regularized); angles constant. + double cell_len[3] = {input_cell.a, input_cell.b, input_cell.c}; + const double cell_len0[3] = {input_cell.a, input_cell.b, input_cell.c}; + double cell_ang[3]; + InitCellAngles(system, input_cell, cell_ang); + + double eff_len[3], eff_alpha, eff_beta, eff_gamma; // recomputed each round from cell_len/ang + + // Per-frame orientation seeds. + std::vector> orient(frames.size()); + for (size_t i = 0; i < frames.size(); ++i) + OrientationSeed(system, frames[i].lattice, orient[i].data()); + + // Reciprocal length of one detector pixel at the nominal geometry, used to scale the robust + // loss (~3 px) into the reciprocal-space residual units. + const double recip_per_px = pixel / (distance_mm * lambda); + const double loss_scale = 3.0 * recip_per_px; + + // Tolerance schedule (fractional-hkl norm), tightening each round like the per-image optimizer. + auto round_tol = [&](int r) { + if (settings.rounds <= 1) return 0.2; + const double t = static_cast(r) / (settings.rounds - 1); + return 0.3 - (0.3 - 0.1) * t; + }; + + for (int round = 0; round < settings.rounds; ++round) { + EffectiveCell(system, cell_len, cell_ang, eff_len, eff_alpha, eff_beta, eff_gamma); + + DiffractionGeometry cur = nominal_geom; + cur.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]).DetectorDistance_mm(distance_mm); + + const double tol = round_tol(round); + const double tol_sq = tol * tol; + + ceres::Problem problem; + auto ordering = std::make_shared(); + + int spots_used = 0, frames_used = 0; + for (size_t i = 0; i < frames.size(); ++i) { + const CrystalLattice latt = AngleAxisAndCellToLattice( + orient[i].data(), eff_len, eff_alpha, eff_beta, eff_gamma); + const Coord v0 = latt.Vec0(), v1 = latt.Vec1(), v2 = latt.Vec2(); + + int used_in_frame = 0; + for (const auto &s : frames[i].spots) { + const Coord recip = cur.DetectorToRecip(s.x, s.y); + const double hf = recip * v0, kf = recip * v1, lf = recip * v2; + const double h = std::round(hf), k = std::round(kf), l = std::round(lf); + const double dn = (hf - h) * (hf - h) + (kf - k) * (kf - k) + (lf - l) * (lf - l); + if (dn > tol_sq) + continue; + + problem.AddResidualBlock( + new ceres::AutoDiffCostFunction( + new XtalResidual(s.x, s.y, lambda, pixel, rot3, 0.0, h, k, l, system)), + new ceres::CauchyLoss(loss_scale), + beam, &distance_mm, detector_rot, rot_vec, orient[i].data(), cell_len, cell_ang); + ++used_in_frame; + ++spots_used; + } + if (used_in_frame > 0) { + ordering->AddElementToGroup(orient[i].data(), 0); // eliminate orientations first + ++frames_used; + } + } + + if (spots_used < settings.min_spots_total) { + if (round == 0) + return result; // ok stays false + break; // keep the previous round's solution + } + + // Cell-length regularizer (weight ~ coeff * sqrt(N_obs) / L0^2 -> effectively pins the cell). + const double sqrtN = std::sqrt(static_cast(spots_used)); + CellLengthRegularizer *reg = new CellLengthRegularizer{ + settings.cell_reg_coeff * sqrtN / (cell_len0[0] * cell_len0[0]), + settings.cell_reg_coeff * sqrtN / (cell_len0[1] * cell_len0[1]), + settings.cell_reg_coeff * sqrtN / (cell_len0[2] * cell_len0[2]), + cell_len0[0], cell_len0[1], cell_len0[2]}; + problem.AddResidualBlock( + new ceres::AutoDiffCostFunction(reg), nullptr, cell_len); + + // Shared blocks go in the second elimination group; constant ones are stripped by Ceres. + for (double *b : {beam, &distance_mm, detector_rot, rot_vec, cell_len, cell_ang}) + ordering->AddElementToGroup(b, 1); + + // Beam + distance bounds around the nominal value; cell lengths loosely bounded (the + // regularizer does the anchoring). + problem.SetParameterLowerBound(beam, 0, nominal_geom.GetBeamX_pxl() - settings.beam_range_px); + problem.SetParameterUpperBound(beam, 0, nominal_geom.GetBeamX_pxl() + settings.beam_range_px); + problem.SetParameterLowerBound(beam, 1, nominal_geom.GetBeamY_pxl() - settings.beam_range_px); + problem.SetParameterUpperBound(beam, 1, nominal_geom.GetBeamY_pxl() + settings.beam_range_px); + problem.SetParameterLowerBound(&distance_mm, 0, + std::max(1.0, nominal_geom.GetDetectorDistance_mm() - settings.distance_range_mm)); + problem.SetParameterUpperBound(&distance_mm, 0, + nominal_geom.GetDetectorDistance_mm() + settings.distance_range_mm); + for (int j = 0; j < 3; ++j) { + problem.SetParameterLowerBound(cell_len, j, 0.7 * cell_len0[j]); + problem.SetParameterUpperBound(cell_len, j, 1.3 * cell_len0[j]); + } + + problem.SetParameterBlockConstant(detector_rot); // tilt not refined + problem.SetParameterBlockConstant(rot_vec); + problem.SetParameterBlockConstant(cell_ang); // angles not refined in v1 + + ceres::Solver::Options options; + options.linear_solver_type = ceres::DENSE_SCHUR; + options.linear_solver_ordering = ordering; + options.num_threads = std::max(1, settings.num_threads); + options.max_num_iterations = 100; + options.max_solver_time_in_seconds = 60.0; + options.logging_type = ceres::LoggingType::SILENT; + options.minimizer_progress_to_stdout = false; + + ceres::Solver::Summary summary; + ceres::Solve(options, &problem, &summary); + + result.frames_used = frames_used; + result.spots_used = spots_used; + } + + // Final fit quality: median detector residual of the retained spots at the refined geometry. + EffectiveCell(system, cell_len, cell_ang, eff_len, eff_alpha, eff_beta, eff_gamma); + DiffractionGeometry cur = nominal_geom; + cur.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]).DetectorDistance_mm(distance_mm); + std::vector residuals; + for (size_t i = 0; i < frames.size(); ++i) { + const CrystalLattice latt = AngleAxisAndCellToLattice( + orient[i].data(), eff_len, eff_alpha, eff_beta, eff_gamma); + const Coord a = latt.Astar(), b = latt.Bstar(), c = latt.Cstar(); + const Coord v0 = latt.Vec0(), v1 = latt.Vec1(), v2 = latt.Vec2(); + for (const auto &s : frames[i].spots) { + const Coord recip = cur.DetectorToRecip(s.x, s.y); + const double hf = recip * v0, kf = recip * v1, lf = recip * v2; + const double h = std::round(hf), k = std::round(kf), l = std::round(lf); + const double dn = (hf - h) * (hf - h) + (kf - k) * (kf - k) + (lf - l) * (lf - l); + if (dn > 0.1 * 0.1) + continue; + const Coord recip_pred = a * (float)h + b * (float)k + c * (float)l; + const auto [px, py] = cur.RecipToDetector(recip_pred); + if (std::isfinite(px) && std::isfinite(py)) + residuals.push_back(std::hypot(px - s.x, py - s.y)); + } + } + if (!residuals.empty()) { + std::nth_element(residuals.begin(), residuals.begin() + residuals.size() / 2, residuals.end()); + result.median_residual_px = residuals[residuals.size() / 2]; + } + + result.ok = result.spots_used >= settings.min_spots_total && result.frames_used > 0; + result.beam_x_px = beam[0]; + result.beam_y_px = beam[1]; + result.distance_mm = distance_mm; + result.cell = BuildUnitCell(system, cell_len, cell_ang); + return result; + } catch (...) { + result.ok = false; + return result; + } +} diff --git a/image_analysis/geom_refinement/GeometryRefiner.h b/image_analysis/geom_refinement/GeometryRefiner.h new file mode 100644 index 00000000..8bb3b787 --- /dev/null +++ b/image_analysis/geom_refinement/GeometryRefiner.h @@ -0,0 +1,59 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "../../common/CrystalLattice.h" +#include "../../common/DiffractionGeometry.h" +#include "../../common/UnitCell.h" +#include "gemmi/symmetry.hpp" + +// Offline global (multi-frame) geometry refinement for serial stills. Pass 1 indexes many frames +// independently; this joint bundle adjustment then determines the ONE shared detector geometry (beam +// centre, distance, cell scale) that per-image refinement cannot reach from a single sparse still, so +// re-indexing with it converts more frames. Mirrors the validated Python prototype: one Ceres problem +// with shared beam/distance/cell blocks and a per-frame orientation block, robust loss, and a cell +// regularizer that anchors the (known) cell to break the distance<->cell-scale degeneracy at low +// resolution. Detector tilt (rot1/rot2) is NOT refined (gauge-coupled with the beam; zero indexing gain). + +struct GeomRefineSpot { + float x = 0, y = 0; // observed (raw) spot position, pixels + int32_t h = 0, k = 0, l = 0; // integer Miller index assigned in pass 1 +}; + +struct GeomRefineFrame { + CrystalLattice lattice; // per-frame indexed lattice (real space), for the orientation seed + std::vector spots; // indexed spots on this frame +}; + +struct GeometryRefinerSettings { + gemmi::CrystalSystem crystal_system = gemmi::CrystalSystem::Triclinic; + // Cell-scale regularizer strength (anchors the cell lengths to the input cell so the otherwise + // degenerate distance is determined). Large => cell effectively fixed (the validated safe regime). + double cell_reg_coeff = 30.0; + int rounds = 3; // hkl-reassignment / tolerance-tightening rounds + double beam_range_px = 8.0; // beam bound around the nominal value + double distance_range_mm = 8.0; // distance bound around the nominal value + int min_spots_total = 100; // refuse to refine below this many indexed spots + int num_threads = 1; +}; + +struct GeometryRefinerResult { + bool ok = false; + double beam_x_px = 0, beam_y_px = 0; // refined beam centre (== direct beam, tilt not refined) + double distance_mm = 0; + UnitCell cell{}; // refined cell + int frames_used = 0; + int spots_used = 0; + double median_residual_px = 0; // final per-spot detector residual (fit quality) +}; + +// nominal_geom supplies the starting beam/distance/tilt/wavelength/pixel size; input_cell is anchored by +// the regularizer. Returns ok=false (and leaves the caller's geometry untouched) on any failure. +GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_geom, + const UnitCell &input_cell, + const std::vector &frames, + const GeometryRefinerSettings &settings); diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp index aa7f4b95..acf4071f 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.cpp +++ b/image_analysis/geom_refinement/XtalOptimizer.cpp @@ -5,203 +5,11 @@ #include #include "XtalOptimizer.h" +#include "XtalResidual.h" #include "ceres/ceres.h" #include "ceres/rotation.h" #include "LatticeReduction.h" -struct XtalResidual { - XtalResidual(double x, double y, - double lambda, - double pixel_size, - double rot3, - double angle_rad, - double exp_h, double exp_k, - double exp_l, - gemmi::CrystalSystem symmetry) - : obs_x(x), obs_y(y), - inv_lambda(1.0/lambda), - pixel_size(pixel_size), - rot3(rot3), - exp_h(exp_h), - exp_k(exp_k), - exp_l(exp_l), - angle_rad(angle_rad), - symmetry(symmetry) { - if (std::fabs(lambda) < 1e-6) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "Lambda cannot be close to zero"); - } - - template - bool operator()(const T *const beam, - const T *const distance_mm, - const T *const detector_rot, - const T *const rotation_axis, - const T *const p0, - const T *const p1, - const T *const p2, - T *residual) const { - // PyFAI convention: poni_rot = Rz(-rot3) * Rx(-rot2) * Ry(+rot1). - // detector_rot[0] = rot1, detector_rot[1] = rot2 are refined; rot3 is fixed - // (e.g. from a PONI import) and baked in here as a constant so that a non-zero - // rot3 is not silently dropped during refinement. - const T rot1 = detector_rot[0]; - const T rot2 = detector_rot[1]; - - // Ry(+rot1): rotation around Y-axis - const T c1 = ceres::cos(rot1); - const T s1 = ceres::sin(rot1); - - // Rx(-rot2): rotation around X-axis with inverted sign (PyFAI left-handed) - const T c2 = ceres::cos(rot2); - const T s2 = ceres::sin(rot2); - - // Rz(-rot3): rotation around Z (beam); constant, identity when rot3 == 0 - const T c3 = T(cos(rot3)); - const T s3 = T(sin(rot3)); - - // Detector coordinates in mm - const T det_x = (T(obs_x) - beam[0]) * T(pixel_size); - const T det_y = (T(obs_y) - beam[1]) * T(pixel_size); - const T det_z = T(distance_mm[0]); - - // Apply Ry(rot1) first: rotate around Y - const T t1_x = c1 * det_x + s1 * det_z; - const T t1_y = det_y; - const T t1_z = -s1 * det_x + c1 * det_z; - - // Then apply Rx(-rot2): rotate around X - const T t2_x = t1_x; - const T t2_y = c2 * t1_y + s2 * t1_z; - const T t2_z = -s2 * t1_y + c2 * t1_z; - - // Then apply Rz(-rot3): rotate around Z (beam) - const T x = c3 * t2_x + s3 * t2_y; - const T y = -s3 * t2_x + c3 * t2_y; - const T z = t2_z; - - // convert to recip space - const T lab_norm = ceres::sqrt(x * x + y * y + z * z); - const T inv_norm = T(1) / lab_norm; - - T recip_raw[3]; - recip_raw[0] = x * inv_norm * T(inv_lambda); - recip_raw[1] = y * inv_norm * T(inv_lambda); - recip_raw[2] = (z * inv_norm - T(1.0)) * T(inv_lambda); - - // Apply goniometer "back-to-start" rotation: - // brings observed reciprocal from image orientation into reference crystal frame - const T aa_back[3] = { - T(angle_rad) * rotation_axis[0], - T(angle_rad) * rotation_axis[1], - T(angle_rad) * rotation_axis[2] - }; - - T recip_obs[3]; - ceres::AngleAxisRotatePoint(aa_back, recip_raw, recip_obs); - - const Eigen::Matrix e_obs_recip(recip_obs[0], recip_obs[1], recip_obs[2]); - - // Build unit cell lengths and B (convention: columns are a, b, c prior to global rotation) - Eigen::Matrix e_uc_len = Eigen::Matrix::Zero(); - Eigen::Matrix B = Eigen::Matrix::Identity(); - - if (symmetry == gemmi::CrystalSystem::Hexagonal) { - e_uc_len << p1[0], p1[0], p1[2]; - B(0, 1) = T(-0.5); // cos(120) - B(1, 1) = T(sqrt(3.0) / 2.0); // sin(120) - } else if (symmetry == gemmi::CrystalSystem::Orthorhombic) { - e_uc_len << p1[0], p1[1], p1[2]; - } else if (symmetry == gemmi::CrystalSystem::Tetragonal) { - e_uc_len << p1[0], p1[0], p1[2]; - } else if (symmetry == gemmi::CrystalSystem::Cubic) { - e_uc_len << p1[0], p1[0], p1[0]; - } else if (symmetry == gemmi::CrystalSystem::Monoclinic) { - // Unique axis b: alpha = gamma = 90°, beta free (angle between a and c) - e_uc_len << p1[0], p1[1], p1[2]; - B(0, 2) = ceres::cos(p2[0]); - B(2, 2) = ceres::sin(p2[0]); - } else { - // Triclinic: p1 = (a,b,c), p2 = (alpha, beta, gamma) in radians - const T ca = ceres::cos(p2[0]); - const T cb = ceres::cos(p2[1]); - const T cg = ceres::cos(p2[2]); - const T sg = ceres::sin(p2[2]); - - e_uc_len << p1[0], p1[1], p1[2]; - - B(0, 0) = T(1); - B(1, 0) = T(0); - B(2, 0) = T(0); - B(0, 1) = cg; - B(1, 1) = sg; - B(2, 1) = T(0); - - // c vector components: - const T cx = cb; - const T cy = (ca - cb * cg) / sg; - const T v = T(1) - cx * cx - cy * cy; - const T cz = (v >= T(0)) ? ceres::sqrt(v) : T(0); - - B(0, 2) = cx; - B(1, 2) = cy; - B(2, 2) = cz; - } - - // Build unrotated direct lattice columns: (B * D), then rotate them by p0. - // This avoids AngleAxisToRotationMatrix + matrix multiplications. - const T L0 = e_uc_len[0]; - const T L1 = e_uc_len[1]; - const T L2 = e_uc_len[2]; - - T col0_unrot[3] = {B(0, 0) * L0, B(1, 0) * L0, B(2, 0) * L0}; - T col1_unrot[3] = {B(0, 1) * L1, B(1, 1) * L1, B(2, 1) * L1}; - T col2_unrot[3] = {B(0, 2) * L2, B(1, 2) * L2, B(2, 2) * L2}; - - T col0_rot[3], col1_rot[3], col2_rot[3]; - ceres::AngleAxisRotatePoint(p0, col0_unrot, col0_rot); - ceres::AngleAxisRotatePoint(p0, col1_unrot, col1_rot); - ceres::AngleAxisRotatePoint(p0, col2_unrot, col2_rot); - - const Eigen::Matrix A(col0_rot[0], col0_rot[1], col0_rot[2]); - const Eigen::Matrix Bv(col1_rot[0], col1_rot[1], col1_rot[2]); - const Eigen::Matrix C(col2_rot[0], col2_rot[1], col2_rot[2]); - - const Eigen::Matrix BxC = Bv.cross(C); - const Eigen::Matrix CxA = C.cross(A); - const Eigen::Matrix AxB = A.cross(Bv); - - const T V = A.dot(BxC); - const T invV = T(1) / V; - - const Eigen::Matrix Astar = BxC * invV; - const Eigen::Matrix Bstar = CxA * invV; - const Eigen::Matrix Cstar = AxB * invV; - - const T h = T(exp_h); - const T k = T(exp_k); - const T l = T(exp_l); - - const Eigen::Matrix e_pred_recip = Astar * h + Bstar * k + Cstar * l; - - residual[0] = e_obs_recip[0] - e_pred_recip[0]; - residual[1] = e_obs_recip[1] - e_pred_recip[1]; - residual[2] = e_obs_recip[2] - e_pred_recip[2]; - - return true; - } - - const double obs_x, obs_y; - const double inv_lambda; - const double pixel_size; - const double rot3; - const double exp_h; - const double exp_k; - const double exp_l; - const double angle_rad; - gemmi::CrystalSystem symmetry; -}; - struct XtalResidualRotationOnlyPrecomp { XtalResidualRotationOnlyPrecomp(const Coord &recip_obs, const CrystalLattice &latt, diff --git a/image_analysis/geom_refinement/XtalResidual.h b/image_analysis/geom_refinement/XtalResidual.h new file mode 100644 index 00000000..6bdcbc3c --- /dev/null +++ b/image_analysis/geom_refinement/XtalResidual.h @@ -0,0 +1,214 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include + +#include "ceres/ceres.h" +#include "ceres/rotation.h" +#include "gemmi/symmetry.hpp" + +#include "../../common/JFJochException.h" + +// Detector -> reciprocal geometry residual, shared by the per-image XtalOptimizer (one lattice, one +// frame) and the offline GeometryRefiner (shared beam/distance/cell blocks, one orientation block per +// frame). Parameter blocks: beam(2), distance_mm(1), detector_rot(2 = rot1,rot2), rotation_axis(3), +// p0 = lattice orientation angle-axis(3), p1 = unit-cell lengths(3), p2 = unit-cell angles(3, rad). +// The residual is the difference between the observed reciprocal vector (from the detector position via +// the current beam/distance/tilt) and the predicted one (h*a* + k*b* + l*c* for the current cell + +// orientation), in Angstrom^-1. rot3 and the goniometer angle_rad are baked in as constants. +struct XtalResidual { + XtalResidual(double x, double y, + double lambda, + double pixel_size, + double rot3, + double angle_rad, + double exp_h, double exp_k, + double exp_l, + gemmi::CrystalSystem symmetry) + : obs_x(x), obs_y(y), + inv_lambda(1.0/lambda), + pixel_size(pixel_size), + rot3(rot3), + exp_h(exp_h), + exp_k(exp_k), + exp_l(exp_l), + angle_rad(angle_rad), + symmetry(symmetry) { + if (std::fabs(lambda) < 1e-6) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Lambda cannot be close to zero"); + } + + template + bool operator()(const T *const beam, + const T *const distance_mm, + const T *const detector_rot, + const T *const rotation_axis, + const T *const p0, + const T *const p1, + const T *const p2, + T *residual) const { + // PyFAI convention: poni_rot = Rz(-rot3) * Rx(-rot2) * Ry(+rot1). + // detector_rot[0] = rot1, detector_rot[1] = rot2 are refined; rot3 is fixed + // (e.g. from a PONI import) and baked in here as a constant so that a non-zero + // rot3 is not silently dropped during refinement. + const T rot1 = detector_rot[0]; + const T rot2 = detector_rot[1]; + + // Ry(+rot1): rotation around Y-axis + const T c1 = ceres::cos(rot1); + const T s1 = ceres::sin(rot1); + + // Rx(-rot2): rotation around X-axis with inverted sign (PyFAI left-handed) + const T c2 = ceres::cos(rot2); + const T s2 = ceres::sin(rot2); + + // Rz(-rot3): rotation around Z (beam); constant, identity when rot3 == 0 + const T c3 = T(cos(rot3)); + const T s3 = T(sin(rot3)); + + // Detector coordinates in mm + const T det_x = (T(obs_x) - beam[0]) * T(pixel_size); + const T det_y = (T(obs_y) - beam[1]) * T(pixel_size); + const T det_z = T(distance_mm[0]); + + // Apply Ry(rot1) first: rotate around Y + const T t1_x = c1 * det_x + s1 * det_z; + const T t1_y = det_y; + const T t1_z = -s1 * det_x + c1 * det_z; + + // Then apply Rx(-rot2): rotate around X + const T t2_x = t1_x; + const T t2_y = c2 * t1_y + s2 * t1_z; + const T t2_z = -s2 * t1_y + c2 * t1_z; + + // Then apply Rz(-rot3): rotate around Z (beam) + const T x = c3 * t2_x + s3 * t2_y; + const T y = -s3 * t2_x + c3 * t2_y; + const T z = t2_z; + + // convert to recip space + const T lab_norm = ceres::sqrt(x * x + y * y + z * z); + const T inv_norm = T(1) / lab_norm; + + T recip_raw[3]; + recip_raw[0] = x * inv_norm * T(inv_lambda); + recip_raw[1] = y * inv_norm * T(inv_lambda); + recip_raw[2] = (z * inv_norm - T(1.0)) * T(inv_lambda); + + // Apply goniometer "back-to-start" rotation: + // brings observed reciprocal from image orientation into reference crystal frame + const T aa_back[3] = { + T(angle_rad) * rotation_axis[0], + T(angle_rad) * rotation_axis[1], + T(angle_rad) * rotation_axis[2] + }; + + T recip_obs[3]; + ceres::AngleAxisRotatePoint(aa_back, recip_raw, recip_obs); + + const Eigen::Matrix e_obs_recip(recip_obs[0], recip_obs[1], recip_obs[2]); + + // Build unit cell lengths and B (convention: columns are a, b, c prior to global rotation) + Eigen::Matrix e_uc_len = Eigen::Matrix::Zero(); + Eigen::Matrix B = Eigen::Matrix::Identity(); + + if (symmetry == gemmi::CrystalSystem::Hexagonal) { + e_uc_len << p1[0], p1[0], p1[2]; + B(0, 1) = T(-0.5); // cos(120) + B(1, 1) = T(sqrt(3.0) / 2.0); // sin(120) + } else if (symmetry == gemmi::CrystalSystem::Orthorhombic) { + e_uc_len << p1[0], p1[1], p1[2]; + } else if (symmetry == gemmi::CrystalSystem::Tetragonal) { + e_uc_len << p1[0], p1[0], p1[2]; + } else if (symmetry == gemmi::CrystalSystem::Cubic) { + e_uc_len << p1[0], p1[0], p1[0]; + } else if (symmetry == gemmi::CrystalSystem::Monoclinic) { + // Unique axis b: alpha = gamma = 90°, beta free (angle between a and c) + e_uc_len << p1[0], p1[1], p1[2]; + B(0, 2) = ceres::cos(p2[0]); + B(2, 2) = ceres::sin(p2[0]); + } else { + // Triclinic: p1 = (a,b,c), p2 = (alpha, beta, gamma) in radians + const T ca = ceres::cos(p2[0]); + const T cb = ceres::cos(p2[1]); + const T cg = ceres::cos(p2[2]); + const T sg = ceres::sin(p2[2]); + + e_uc_len << p1[0], p1[1], p1[2]; + + B(0, 0) = T(1); + B(1, 0) = T(0); + B(2, 0) = T(0); + B(0, 1) = cg; + B(1, 1) = sg; + B(2, 1) = T(0); + + // c vector components: + const T cx = cb; + const T cy = (ca - cb * cg) / sg; + const T v = T(1) - cx * cx - cy * cy; + const T cz = (v >= T(0)) ? ceres::sqrt(v) : T(0); + + B(0, 2) = cx; + B(1, 2) = cy; + B(2, 2) = cz; + } + + // Build unrotated direct lattice columns: (B * D), then rotate them by p0. + // This avoids AngleAxisToRotationMatrix + matrix multiplications. + const T L0 = e_uc_len[0]; + const T L1 = e_uc_len[1]; + const T L2 = e_uc_len[2]; + + T col0_unrot[3] = {B(0, 0) * L0, B(1, 0) * L0, B(2, 0) * L0}; + T col1_unrot[3] = {B(0, 1) * L1, B(1, 1) * L1, B(2, 1) * L1}; + T col2_unrot[3] = {B(0, 2) * L2, B(1, 2) * L2, B(2, 2) * L2}; + + T col0_rot[3], col1_rot[3], col2_rot[3]; + ceres::AngleAxisRotatePoint(p0, col0_unrot, col0_rot); + ceres::AngleAxisRotatePoint(p0, col1_unrot, col1_rot); + ceres::AngleAxisRotatePoint(p0, col2_unrot, col2_rot); + + const Eigen::Matrix A(col0_rot[0], col0_rot[1], col0_rot[2]); + const Eigen::Matrix Bv(col1_rot[0], col1_rot[1], col1_rot[2]); + const Eigen::Matrix C(col2_rot[0], col2_rot[1], col2_rot[2]); + + const Eigen::Matrix BxC = Bv.cross(C); + const Eigen::Matrix CxA = C.cross(A); + const Eigen::Matrix AxB = A.cross(Bv); + + const T V = A.dot(BxC); + const T invV = T(1) / V; + + const Eigen::Matrix Astar = BxC * invV; + const Eigen::Matrix Bstar = CxA * invV; + const Eigen::Matrix Cstar = AxB * invV; + + const T h = T(exp_h); + const T k = T(exp_k); + const T l = T(exp_l); + + const Eigen::Matrix e_pred_recip = Astar * h + Bstar * k + Cstar * l; + + residual[0] = e_obs_recip[0] - e_pred_recip[0]; + residual[1] = e_obs_recip[1] - e_pred_recip[1]; + residual[2] = e_obs_recip[2] - e_pred_recip[2]; + + return true; + } + + const double obs_x, obs_y; + const double inv_lambda; + const double pixel_size; + const double rot3; + const double exp_h; + const double exp_k; + const double exp_l; + const double angle_rad; + gemmi::CrystalSystem symmetry; +}; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 72406ddb..656e0d38 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -23,6 +24,7 @@ #include "../writer/FileWriter.h" #include "../image_analysis/MXAnalysisWithoutFPGA.h" #include "../image_analysis/IndexAndRefine.h" +#include "../image_analysis/geom_refinement/GeometryRefiner.h" #include "../image_analysis/indexing/IndexerThreadPool.h" #include "../image_analysis/azint/AzIntEngineCPU.h" #include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" @@ -74,6 +76,144 @@ Rugnux::Rugnux(JFJochHDF5Reader &reader, DiffractionExperiment experiment, : reader_(reader), experiment_(std::move(experiment)), pixel_mask_(std::move(pixel_mask)), config_(std::move(config)) {} +void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_process, + RugnuxObserver *observer) { + Logger logger("Rugnux"); + + // Rotation has its own two-pass; this is stills only. + if (experiment_.IsRotationIndexing()) { + logger.Warning("--refine-geometry is stills-only; rotation data uses two-pass indexing. Ignoring."); + return; + } + // The bundle anchors a known cell; without one there is nothing to refine against. + const auto cell = experiment_.GetUnitCell(); + if (!cell.has_value()) { + logger.Warning("--refine-geometry needs a reference unit cell (-C / -S / reference MTZ). Skipping."); + return; + } + + const int refine_frames = std::max(1, config_.refine_geometry.value()); + if (observer) + observer->OnPhase("Geometry refinement (first pass)"); + + const auto dataset = reader_.GetDataset(); + + // Index a spread sample large enough to yield ~refine_frames strong frames even at a low hit rate. + const int sample_budget = std::min(images_to_process, std::max(refine_frames * 50, 8000)); + const std::vector sample = select_equally_spaced_image_ordinals(images_to_process, sample_budget); + + // First-pass engines, built from the current (nominal) geometry: an index-only pass (no + // integration) to obtain each frame's spots + assigned HKL + orientation. + AzimuthalIntegrationMapping mapping(experiment_, pixel_mask_); + IndexerThreadPool pool(experiment_.GetIndexingSettings()); + IndexAndRefine indexer(experiment_, &pool, /*retain_outcomes=*/false); + + auto pass_settings = config_.spot_finding; + pass_settings.enable = true; + pass_settings.indexing = true; + pass_settings.quick_integration = false; + + std::vector frames; + std::mutex frames_mutex; + std::atomic next_i = 0; + + auto worker = [&]() { + pin_gpu(); // round-robin per worker thread; must precede engine construction + MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, indexer); + AzimuthalIntegrationProfile profile(mapping); + + while (!cancelled_) { + const int idx = next_i.fetch_add(1); + if (idx >= static_cast(sample.size())) break; + const int ordinal = sample[idx]; + const int image_idx = start_image + ordinal * config_.stride; + + std::shared_ptr img; + try { + img = reader_.GetRawImage(image_idx); + } catch (const std::exception &e) { + logger.Warning("Geometry refinement: failed to load image {}: {}", image_idx, e.what()); + continue; + } + if (!img) continue; + + DataMessage msg{}; + msg.image = img->image; + msg.number = ordinal; + msg.original_number = image_idx; + if (dataset->efficiency.size() > image_idx) + msg.image_collection_efficiency = dataset->efficiency[image_idx]; + + try { + analysis.Analyze(msg, profile, pass_settings); + } catch (const std::exception &e) { + continue; + } + if (!msg.indexing_result.value_or(false) || !msg.indexing_lattice.has_value()) + continue; + + GeomRefineFrame f; + f.lattice = *msg.indexing_lattice; + for (const auto &s : msg.spots) + if (s.indexed && s.lattice == 0) + f.spots.push_back(GeomRefineSpot{s.x, s.y, + static_cast(s.h), static_cast(s.k), static_cast(s.l)}); + if (f.spots.size() >= 6) { + const std::unique_lock ul(frames_mutex); + frames.push_back(std::move(f)); + } + } + }; + + std::vector> futures; + futures.reserve(config_.nthreads); + for (int i = 0; i < config_.nthreads; ++i) + futures.push_back(std::async(std::launch::async, worker)); + for (auto &fut : futures) + fut.get(); + + if (cancelled_) + return; + + // Keep the strongest frames (most indexed spots) for the bundle - the true cell indexes many + // spots per frame, and orientation diversity comes for free from independent serial stills. + constexpr int MIN_STRONG_SPOTS = 20; + std::vector strong; + for (auto &f : frames) + if (static_cast(f.spots.size()) >= MIN_STRONG_SPOTS) + strong.push_back(std::move(f)); + std::sort(strong.begin(), strong.end(), + [](const GeomRefineFrame &a, const GeomRefineFrame &b) { return a.spots.size() > b.spots.size(); }); + if (static_cast(strong.size()) > refine_frames) + strong.resize(refine_frames); + + logger.Info("Geometry refinement: indexed {} of {} sampled frames, bundling {} strong frames (>= {} spots)", + frames.size(), sample.size(), strong.size(), MIN_STRONG_SPOTS); + + GeometryRefinerSettings gr_settings; + gr_settings.crystal_system = experiment_.GetCrystalSystem(); + gr_settings.num_threads = config_.nthreads; + + const GeometryRefinerResult r = RefineGlobalGeometry( + experiment_.GetDiffractionGeometry(), *cell, strong, gr_settings); + if (!r.ok) { + logger.Warning("Geometry refinement did not run (too few strong frames / did not converge); " + "keeping the input geometry"); + return; + } + + logger.Info("Geometry refinement: beam ({:.2f}, {:.2f}) -> ({:.2f}, {:.2f}) px, distance {:.4f} -> {:.4f} mm", + experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(), r.beam_x_px, r.beam_y_px, + experiment_.GetDetectorDistance_mm(), r.distance_mm); + logger.Info("Geometry refinement: cell a/b/c {:.3f}/{:.3f}/{:.3f} -> {:.3f}/{:.3f}/{:.3f} A " + "(median residual {:.3f} px, {} frames, {} spots)", + cell->a, cell->b, cell->c, r.cell.a, r.cell.b, r.cell.c, + r.median_residual_px, r.frames_used, r.spots_used); + + experiment_.BeamX_pxl(r.beam_x_px).BeamY_pxl(r.beam_y_px).DetectorDistance_mm(r.distance_mm); + experiment_.SetUnitCell(r.cell); +} + ProcessResult Rugnux::Run(RugnuxObserver *observer) { Logger logger("Rugnux"); ProcessResult result; @@ -130,6 +270,13 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { experiment_.Goniometer(shifted); } + // Stills global geometry-refinement first pass: bundle-adjust the shared beam/distance/cell from a + // sample of strongly-indexed frames and apply it to experiment_, so everything below (azimuthal + // mapping, the main index+integrate pass, scaling) uses the refined geometry. Rotation is skipped + // inside (it has its own two-pass). + if (full && config_.refine_geometry.has_value()) + RefineStillsGeometry(start_image, end_image, images_to_process, observer); + AzimuthalIntegrationMapping mapping(experiment_, pixel_mask_); JFJochReceiverPlots plots; diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index da918b96..f8aa5dd5 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -61,6 +61,12 @@ struct ProcessConfig { int rotation_indexing_image_count = 100; std::optional forced_rotation_lattice; + // Stills global geometry-refinement two-pass (FullAnalysis, stills only). When set, an extra first + // pass indexes a spread sample of frames, a single joint bundle adjustment over the strongest ones + // determines the shared beam/distance/cell that per-image refinement cannot reach, and the main pass + // re-indexes with it. The value is the number of strong frames fed to the bundle (--refine-geometry). + std::optional refine_geometry; + // Scaling / merging (FullAnalysis). reference_data (from a reference MTZ) also drives the // per-image live scaling path when present. bool run_scaling = false; @@ -129,6 +135,13 @@ class Rugnux { ProcessConfig config_; std::atomic cancelled_{false}; + // Stills global geometry-refinement first pass (config_.refine_geometry): index a spread sample of + // frames, bundle-adjust the shared beam/distance/cell from the strongest ones, and apply the result + // to experiment_ so the main pass re-indexes with it. No-op (leaves experiment_ unchanged) if it + // cannot run. Stills only; rotation has its own two-pass. + void RefineStillsGeometry(int start_image, int end_image, int images_to_process, + RugnuxObserver *observer); + public: Rugnux(JFJochHDF5Reader &reader, DiffractionExperiment experiment, PixelMask pixel_mask, ProcessConfig config); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 371cc39c..a31bdaa8 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -76,6 +76,7 @@ void print_usage() { std::cout << " -S, --space-group Space group number - used for both indexing and scaling" << std::endl; std::cout << " -C, --unit-cell Fix reference unit cell: \"a,b,c,alpha,beta,gamma\"" << std::endl; std::cout << " -r, --refine Geometry refinement algorithm (none|orientation|beam_and_lattice|multi); multi tries all three per image and keeps whichever indexes the most spots" << std::endl; + std::cout << " --refine-geometry[=N] Stills: extra first pass that bundle-adjusts the shared beam/distance/cell from N strongly-indexed frames (default: 200), then re-indexes with it (lifts weak-stills indexing)" << std::endl; std::cout << std::endl; std::cout << " Scaling and merging (on by default)" << std::endl; @@ -146,6 +147,7 @@ enum { OPT_SINGLE_PASS_ROTATION, OPT_REDO_ROTATION_SPOTS, OPT_FORCE_ROTATION_LATTICE, + OPT_REFINE_GEOMETRY, OPT_BANDWIDTH, OPT_INTEGRATION_RADIUS, OPT_REJECT_OUTLIERS, @@ -226,6 +228,7 @@ static option long_options[] = { {"polarization", required_argument, nullptr, OPT_POLARIZATION}, {"redo-rotation-spots", no_argument, nullptr, OPT_REDO_ROTATION_SPOTS}, {"force-rotation-lattice", required_argument, nullptr, OPT_FORCE_ROTATION_LATTICE}, + {"refine-geometry", optional_argument, nullptr, OPT_REFINE_GEOMETRY}, {"spot-sigma", required_argument, nullptr, OPT_SPOT_SIGMA}, @@ -476,6 +479,7 @@ int main(int argc, char **argv) { double min_image_cc = 0.0; int64_t scaling_iter = 3; std::optional forced_rotation_lattice; + std::optional refine_geometry; // --refine-geometry[=N]: stills global geometry-refinement pass std::optional bandwidth_fwhm; // relative FWHM of dlambda/lambda @@ -550,6 +554,15 @@ int main(int argc, char **argv) { case OPT_REDO_ROTATION_SPOTS: reuse_rotation_spots = false; break; + case OPT_REFINE_GEOMETRY: { + int n = optarg ? atoi(optarg) : 200; + if (n <= 0) { + logger.Error("Invalid --refine-geometry frame count: {}", optarg ? optarg : ""); + exit(EXIT_FAILURE); + } + refine_geometry = n; + break; + } case OPT_FORCE_ROTATION_LATTICE: { if (rotation_indexing) { logger.Error("Rotation indexing already enabled"); @@ -1374,6 +1387,7 @@ int main(int argc, char **argv) { config.reuse_rotation_spots = reuse_rotation_spots; config.rotation_indexing_image_count = rotation_indexing_image_count; config.forced_rotation_lattice = forced_rotation_lattice; + config.refine_geometry = refine_geometry; config.run_scaling = run_scaling; config.scaling_iter = scaling_iter; config.reference_data = reference_data; -- 2.54.0 From 4c5da20eda87662deebd5ed2bd587d5427204948 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 19:52:46 +0200 Subject: [PATCH 19/90] rugnux: add --partiality-uncertainty stills merge term (partiality-model error) A stills reflection recorded at partiality p carries a systematic intensity error ~(dp/p) that is proportional to and grows as p falls; plain counting sigma misses it, so strong low-p partials are over-trusted in the merge. Add sigma^2 += (c**(1-p))^2 in MergeOnTheFly::CorrectedSigma - the stills-partiality analog of the existing rotation --capture-uncertainty term ((1-captured_fraction)*I in RotationScaleMerge). Inert when partiality==1 (no --still-partiality), and gated on a real systematic (error_model_b > 1, i.e. ISa < 1) so it fires on strong/medium stills but auto-skips weak counting-limited data where it would only over-concentrate the merge and hurt. Opt-in via --partiality-uncertainty (default 0; ~2.5 recommended with --still-partiality). Prototyped (Python replica of the merge): lyso8 CC1/2 +1.7, CCref +4.8, R-free -0.023; LOV CCref +2; harmful on weak OCP (hence the b>1 gate). Full in-binary validation across all serial targets pending. Co-Authored-By: Claude Opus 4.8 (1M context) --- common/ScalingSettings.cpp | 9 +++++++++ common/ScalingSettings.h | 3 +++ image_analysis/scale_merge/Merge.cpp | 21 +++++++++++++++++---- image_analysis/scale_merge/Merge.h | 2 +- rugnux/rugnux_cli.cpp | 8 ++++++++ 5 files changed, 38 insertions(+), 5 deletions(-) diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index f3341e1d..c57ff626 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -191,6 +191,15 @@ double ScalingSettings::GetCaptureUncertaintyCoeff() const { return capture_uncertainty_coeff; } +ScalingSettings &ScalingSettings::PartialityUncertaintyCoeff(double input) { + partiality_uncertainty_coeff = input; + return *this; +} + +double ScalingSettings::GetPartialityUncertaintyCoeff() const { + return partiality_uncertainty_coeff; +} + ScalingSettings &ScalingSettings::MinCapturedFraction(double input) { if (input < 0.0) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index 1056bb1b..9995ecf2 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -29,6 +29,7 @@ class ScalingSettings { // a fraction f<1 of its rocking curve is extrapolated, and the unobserved (1-f) carries a // systematic error ~coeff*(1-f)*I that plain counting sigma misses. 0 = off (baseline). double capture_uncertainty_coeff = 0.0; + double partiality_uncertainty_coeff = 0.0; // Full-level captured-fraction floor for the rot3d combine: drop a reconstructed full whose rocking // curve was only fractionally captured (sum of its partials' partiality < this). Distinct from // min_partiality, which gates individual partials; this gates the assembled full. 0 = off (baseline). @@ -75,6 +76,7 @@ public: ScalingSettings& MinPartiality(double min_partiality); ScalingSettings& ForcedMosaicity(std::optional input); ScalingSettings& CaptureUncertaintyCoeff(double input); + ScalingSettings& PartialityUncertaintyCoeff(double input); ScalingSettings& MinCapturedFraction(double input); ScalingSettings& MinCCForImage(double min_cc_for_image); ScalingSettings& OutlierRejectNsigma(double input); @@ -112,6 +114,7 @@ public: [[nodiscard]] double GetMinPartiality() const; [[nodiscard]] std::optional GetForcedMosaicity() const; [[nodiscard]] double GetCaptureUncertaintyCoeff() const; + [[nodiscard]] double GetPartialityUncertaintyCoeff() const; [[nodiscard]] double GetMinCapturedFraction() const; [[nodiscard]] double GetMinCCForImage() const; [[nodiscard]] double GetOutlierRejectNsigma() const; diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 405a9b26..916e72e9 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -93,7 +93,7 @@ void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id continue; auto hkl = generator(r); auto hkl_key = hkl.pack(); - sigma_corr = CorrectedSigma(I_corr, sigma_corr, hkl_key); + sigma_corr = CorrectedSigma(I_corr, sigma_corr, hkl_key, r.partiality); // Robust outlier rejection: drop this observation if it sits more than // reject_nsigma error-model sigmas from the reflection's median. Needs the active @@ -131,7 +131,7 @@ void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id } } -float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, uint64_t hkl_key) const { +float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, uint64_t hkl_key, float partiality) const { if (!error_model_active) return sigma_corr; @@ -140,8 +140,21 @@ float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, uint64_t hkl const auto it = error_model_mean_I.find(hkl_key); const double I_for_b = (it != error_model_mean_I.end()) ? it->second : I_corr; - const double v = error_model_a * static_cast(sigma_corr) * sigma_corr - + (error_model_b * I_for_b) * (error_model_b * I_for_b); + double v = error_model_a * static_cast(sigma_corr) * sigma_corr + + (error_model_b * I_for_b) * (error_model_b * I_for_b); + + // Partiality-model uncertainty: a reflection recorded at fraction p carries a systematic intensity + // error ~ (dp/p) that is proportional to and grows as p falls - plain counting sigma misses it, + // so strong low-p partials would otherwise be over-trusted. This is the stills-partiality analog of + // the rotation --capture-uncertainty term ((1-captured_fraction)*I in RotationScaleMerge). Inert when + // partiality == 1 (no stills partiality model). Gated on a real systematic (error_model_b > 1, i.e. + // ISa < 1): on weak counting-limited data (small b) it would only over-concentrate the merge and hurt. + const double c = scaling_settings.GetPartialityUncertaintyCoeff(); + if (c > 0.0 && error_model_b > 1.0) { + const double one_minus_p = std::max(0.0, std::min(1.0, 1.0 - static_cast(partiality))); + const double t = c * I_for_b * one_minus_p; + v += t * t; + } return (v > 0.0) ? static_cast(std::sqrt(v)) : sigma_corr; } diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 9c0ee572..32f5c58a 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -109,7 +109,7 @@ class MergeOnTheFly { // observations), so it inflates sigma without biasing the inverse-variance weights - // using the per-observation I_i instead would over-weight down-fluctuated points. std::unordered_map error_model_mean_I; - [[nodiscard]] float CorrectedSigma(float I_corr, float sigma_corr, uint64_t hkl_key) const; + [[nodiscard]] float CorrectedSigma(float I_corr, float sigma_corr, uint64_t hkl_key, float partiality) const; // Optional per-observation outlier rejection: drop observations whose corrected // intensity lies more than reject_nsigma error-model sigmas from the reflection's diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index a31bdaa8..2ce2dab7 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -94,6 +94,7 @@ void print_usage() { std::cout << " --resolution-shells Number of resolution shells in the reported statistics table (default: 10)" << std::endl; std::cout << " --min-partiality Minimum partiality to accept reflection (default: 0.02)" << std::endl; std::cout << " --capture-uncertainty rot3d: systematic sigma ~num*(1-captured_fraction)*I on under-captured fulls (default: 1.0 for rot3d, 0 otherwise)" << std::endl; + std::cout << " --partiality-uncertainty stills: extra merge sigma ~num*(1-partiality)* on partials (use with --still-partiality; auto-gated to error-model b>1 / ISa<1; default 0, ~2.5 recommended)" << std::endl; std::cout << " --min-captured-fraction rot3d: drop a combined full whose rocking curve was captured below this fraction (edge-of-sweep truncated fulls) (default: 0.7 for rotation, 0 otherwise; 0 = off)" << std::endl; std::cout << " --mosaicity Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl; std::cout << " --reject-outliers Per-observation merge outlier rejection, N sigma from the per-reflection median (default: 6 for rot3d, XDS/DIALS-style; 0 = off)" << std::endl; @@ -159,6 +160,7 @@ enum { OPT_STILL_PARTIALITY, OPT_SCALE_FULLS, OPT_CAPTURE_UNCERTAINTY, + OPT_PARTIALITY_UNCERTAINTY, OPT_MIN_CAPTURED_FRACTION, OPT_MOSAICITY, OPT_SMOOTH_G, @@ -239,6 +241,7 @@ static option long_options[] = { {"max-spots", required_argument, nullptr, OPT_MAX_SPOTS}, {"min-partiality", required_argument, nullptr, OPT_MIN_PARTIALITY}, {"capture-uncertainty", required_argument, nullptr, OPT_CAPTURE_UNCERTAINTY}, + {"partiality-uncertainty", required_argument, nullptr, OPT_PARTIALITY_UNCERTAINTY}, {"min-captured-fraction", required_argument, nullptr, OPT_MIN_CAPTURED_FRACTION}, {"mosaicity", required_argument, nullptr, OPT_MOSAICITY}, {"min-image-cc", required_argument, nullptr, OPT_MIN_IMAGE_CC}, @@ -475,6 +478,7 @@ int main(int argc, char **argv) { double min_partiality = 0.02; std::optional min_captured_fraction_arg; // explicit --min-captured-fraction; default depends on rotation std::optional capture_uncertainty_arg; // explicit --capture-uncertainty; default depends on rot3d + std::optional partiality_uncertainty_arg; // --partiality-uncertainty (stills partiality merge sigma) std::optional forced_mosaicity_arg; // diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed double min_image_cc = 0.0; int64_t scaling_iter = 3; @@ -730,6 +734,9 @@ int main(int argc, char **argv) { case OPT_CAPTURE_UNCERTAINTY: capture_uncertainty_arg = parse_double_arg(optarg, "--capture-uncertainty", logger); break; + case OPT_PARTIALITY_UNCERTAINTY: + partiality_uncertainty_arg = parse_double_arg(optarg, "--partiality-uncertainty", logger); + break; case OPT_MIN_CAPTURED_FRACTION: min_captured_fraction_arg = parse_double_arg(optarg, "--min-captured-fraction", logger); break; @@ -1306,6 +1313,7 @@ int main(int argc, char **argv) { // over-extrapolated under-captured fulls and, with the mosaicity fix, lifts rotation ISa/anomalous // substantially. Off for non-rot3d (no combine). An explicit --capture-uncertainty always wins. scaling_settings.CaptureUncertaintyCoeff(capture_uncertainty_arg.value_or(rotation_indexing ? 1.0 : 0.0)); + scaling_settings.PartialityUncertaintyCoeff(partiality_uncertainty_arg.value_or(0.0)); scaling_settings.ForcedMosaicity(forced_mosaicity_arg); scaling_settings.MinCCForImage(min_image_cc / 100.0); // --min-image-cc is in percent; the setting is a fraction scaling_settings.OutlierRejectNsigma( -- 2.54.0 From 4c050b576ceb5cba912e9d5727ecd39d566a676b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 21:53:08 +0200 Subject: [PATCH 20/90] rugnux: de-double spurious rotation supercell in candidate selection The two-pass rotation indexer could keep a doubled-axis supercell over its own primitive sub-cell: when the supercell won candidate (ci) order and the primitive cleared the indexed-fraction hysteresis by less than the 0.05 "prefer-earlier" margin, the supercell survived. On pding4_003 full data this gave P222 on a 65.6x131.3x173 cell instead of P422 on 65.6x65.6x173 - the orthorhombic metric forecloses the 4-fold before the symmetry search ever runs. The FFT peak-finding is correct (the primitive cell is among the candidates, just out-ranked). Add a sub-cell override to the selection loop: adopt a later candidate that is a genuinely smaller cell (> ROT_SUBCELL_VOLUME_RATIO=1.5x smaller volume; a doubling is 2x) indexing at least as many spots (within ROT_SUBCELL_FRAC_SLACK=0.02). That is the signature of a spurious doubling - the primitive always indexes >= its integer multiple, whereas a real superstructure's larger cell indexes MORE (kept by the existing clearly-more branch) and twins share the cell volume (untouched). Full 24-crystal rotation battery (pre-fix vs post-fix): pding4_003 P222->P422 (cell halves to 65.6x65.6x173, R_meas 6.6% CC1/2 99.9% ISa 11.9), pding4_001 holds P422, all 22 other crystals bit-identical in space group. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../rotation_indexer/RotationIndexer.cpp | 24 ++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index 540fe022..355e5b00 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/image_analysis/rotation_indexer/RotationIndexer.cpp @@ -12,6 +12,14 @@ #include namespace { + // Sub-cell override thresholds used in candidate selection to undo a spurious axis doubling: + // a later candidate replaces the chosen cell when it is smaller by more than this volume ratio + // (a doubling is 2x, well past 1.5) and indexes within this fraction slack of it. The slack is + // far below the indexed-fraction gap a real superstructure opens between its true cell and its + // sub-cell, so genuine large cells are kept. + constexpr float ROT_SUBCELL_VOLUME_RATIO = 1.5f; + constexpr float ROT_SUBCELL_FRAC_SLACK = 0.02f; + // Re-express a primitive hexagonal/trigonal lattice in the conventional hexagonal setting // (a = b, gamma = 120). The Niggli-reduced primitive cell carries the two equal-length axes // at gamma = 60; replacing b with b - a opens that angle to 120 without changing the lattice. @@ -224,6 +232,7 @@ void RotationIndexer::RunIndexing() { // Assemble and select serially, in candidate order - identical to refining them one by one. float best_frac = -1.0f; + float best_vol = 0.0f; bool have_best = false; size_t best_ci = 0; XtalOptimizerData best_data; @@ -264,8 +273,21 @@ void RotationIndexer::RunIndexing() { // marginally-higher alternative from displacing the primary when both index poorly (e.g. a // twin, where the accumulated-spot fraction is a noisy proxy) - only a decisively better // cell (a superstructure's true cell vs its sublattice) takes over. - if (!have_best || (frac > best_frac + 0.05f && frac > 0.15f)) { + // Displace the current best when the candidate indexes clearly more, OR when it is a + // genuine sub-cell: a meaningfully smaller cell that still indexes at least as many spots. + // The sub-cell branch unmasks a spurious supercell (axis doubling): the primitive cell + // always indexes >= its integer multiple, so a doubled cell that wins ci-order by the + // hysteresis margin is overridden by its own primitive. A real superstructure's true + // (larger) cell indexes MORE than its sub-cell and is kept by the clearly-more branch; + // twin lattices share the cell volume, so this never disturbs twin selection. + const float cand_vol = std::abs(data.latt.CalcVolume()); + const bool clearly_more = frac > best_frac + 0.05f && frac > 0.15f; + const bool smaller_subcell = have_best && frac > 0.15f + && frac >= best_frac - ROT_SUBCELL_FRAC_SLACK + && cand_vol < best_vol / ROT_SUBCELL_VOLUME_RATIO; + if (!have_best || clearly_more || smaller_subcell) { best_frac = frac; + best_vol = cand_vol; have_best = true; best_data = std::move(data); best_sr = sr; -- 2.54.0 From 768a926e0c9c210e240f1d2b55dcb36d59ce9492 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 23:22:38 +0200 Subject: [PATCH 21/90] rugnux: reject 2nd-lattice overlaps in the space-group search (E^2 cap) On a two-lattice crystal a second lattice deposits intensity on one reciprocal position but not its symmetry mate, so a contaminated reflection is a one-sided resolution-normalised-E outlier that poisons the operator I(h)/I(Rh) correlation in SearchSpaceGroup Stage A. On EP_cs_02-424 under -A this dropped the monoclinic 2-fold CC to 0.326 (below the 0.5 gate) -> the crystal was under-called P1 instead of P2_1. The existing --reject-outliers cannot see it: it rejects WITHIN a P1 orbit, but the poison is a BETWEEN-orbit effect (inflates orbit h, not its mate Rh). Add SearchSpaceGroupOptions::max_e_squared_for_cc (default 9.0, i.e. E>3): after the resolution-normalised Esq[] is computed, drop the extreme-E tail from the correlation pairs only (the absence stage keeps the full range - that is where the screw signal lives). Clean Wilson-distributed data almost never reaches E^2=9 (P ~ 0.01-0.3%), so it is self-targeting - it trims the overlap tail without touching genuine reflections. EP_cs_02-424 -A: 2-fold CC 0.326 -> 0.64 -> P2_1 (= XDS). Full 24-crystal rotation battery (pre/post): only EP_cs_02-424 changed (P1 -> P2_1, SG match 20 -> 21/24); all 23 others bit-identical, including the marginal cubic Ins_I_3 (still I23, carried by the systematic-b rescue). A Friedel-merge-the-search alternative was tried and rejected (net-negative: did not fix EP424 and regressed Ins_I_3 -> I222). Co-Authored-By: Claude Opus 4.8 (1M context) --- image_analysis/scale_merge/SearchSpaceGroup.cpp | 8 ++++++++ image_analysis/scale_merge/SearchSpaceGroup.h | 9 +++++++++ 2 files changed, 17 insertions(+) diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index b5da6777..8f581bc8 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -222,6 +222,14 @@ SearchSpaceGroupResult SearchSpaceGroup( } } + // Overlap guard (Stage A / correlation only): drop the extreme resolution-normalised-E tail, which + // on a two-lattice crystal is the one-sided overlap contamination that poisons the operator CC. + // See SearchSpaceGroupOptions::max_e_squared_for_cc. Absences (pass_absence) keep the full range. + if (opt.max_e_squared_for_cc > 0.0) + for (size_t i = 0; i < n; ++i) + if (pass_cc[i] && Esq[i] > opt.max_e_squared_for_cc) + pass_cc[i] = false; + std::unordered_map key_to_index; key_to_index.reserve(n * 2); for (size_t i = 0; i < n; ++i) diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index c8e9ca18..ce0a506e 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -65,6 +65,15 @@ struct SearchSpaceGroupOptions { // keep weak reflections - that is where the screw-axis signal lives). double min_i_over_sigma = 0.0; + // Drop reflections STRONGER than this resolution-normalised E^2 = I/(shell) from the correlation + // stage only. A second lattice deposits intensity on one reciprocal position but not its symmetry + // mate, so an overlap-contaminated reflection is a one-sided E^2 outlier that poisons an operator's + // I(h)/I(Rh) correlation (it flipped EP_cs_02-424 P2_1->P1 under -A: excluding the E^2>9 tail, ~0.4% + // of reflections, lifted the 2-fold CC 0.33->0.53 back over the gate). Clean Wilson-distributed data + // almost never reaches E^2=9 (P(E^2>9) ~ 0.01-0.3%), so this removes essentially nothing there and + // only trims the overlap tail. 0 disables. Absences are unaffected (they need the weak tail). + double max_e_squared_for_cc = 9.0; + // --- Stage A: point group --- // A rotation is accepted as a real symmetry when its I(h)/I(Rh) correlation reaches this over // at least min_pairs_per_operator independent pairs. -- 2.54.0 From 135765e5b58708b67adfad548f68fb1602908cc5 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 15 Jul 2026 07:27:25 +0200 Subject: [PATCH 22/90] rugnux: default --refine-geometry on for stills-with-cell; accept F/FP reference columns Two stills usability fixes. 1. --refine-geometry now defaults ON for stills whenever a reference cell is available (-C or a reference MTZ) - exactly the case where the geometry bundle-adjust can act (it anchors on a known cell) and where it lifts weak/ sparse-stills indexing (OCP +42%, KR2 +108% indexed in the target study). It stays a no-op for rotation (own two-pass) and de-novo stills (no cell yet), so auto-enabling it only where it does something avoids spurious "skipping" warnings. --refine-geometry=off opts out; explicit --refine-geometry[=N] forces it. 2. The reference-MTZ auto column selection now falls back to a plain amplitude column (FP / FOBS / F / FC, squared to an intensity) after F-model and the intensity (J) columns. A deposition that carries only structure-factor amplitudes (e.g. KR2's 8cl8, whose only usable column is FP) now seeds CCref without needing an explicit --reference-column FP. Co-Authored-By: Claude Opus 4.8 (1M context) --- image_analysis/LoadFCalcFromMtz.cpp | 17 ++++++++++++++++- rugnux/rugnux_cli.cpp | 24 +++++++++++++++++++++--- 2 files changed, 37 insertions(+), 4 deletions(-) diff --git a/image_analysis/LoadFCalcFromMtz.cpp b/image_analysis/LoadFCalcFromMtz.cpp index 00d36b5a..49dc4549 100644 --- a/image_analysis/LoadFCalcFromMtz.cpp +++ b/image_analysis/LoadFCalcFromMtz.cpp @@ -36,6 +36,21 @@ const gemmi::Mtz::Column* SelectDefaultColumn(const gemmi::Mtz& mtz, bool& squar square = false; return &c; } + // Fall back to a plain structure-factor amplitude (type F): a deposited/observed FP or FOBS, or a + // lone FC, squared to an intensity like F-model. Lets a reference MTZ that carries only amplitudes + // (no F-model or IMEAN/I - e.g. the KR2 8cl8 deposition, which has only FP) still seed CCref + // without an explicit --reference-column. + for (const char* label : {"FP", "FOBS", "F", "FC", "Fobs", "F-obs"}) { + if (const auto* col = mtz.column_with_label(label, nullptr, 'F')) { + square = true; + return col; + } + } + for (const auto& c : mtz.columns) + if (c.type == 'F') { + square = true; + return &c; + } return nullptr; } @@ -63,7 +78,7 @@ ReferenceMtzData LoadReferenceMtz(const std::string& path, } else { col = SelectDefaultColumn(mtz, out.squared); if (col == nullptr) - throw std::runtime_error("MTZ has no F-model or intensity (J) column to use as reference"); + throw std::runtime_error("MTZ has no F-model, intensity (J) or amplitude (F) column to use as reference"); out.default_column = true; } out.used_column = col->label; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 2ce2dab7..6ce8e0de 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -76,7 +76,7 @@ void print_usage() { std::cout << " -S, --space-group Space group number - used for both indexing and scaling" << std::endl; std::cout << " -C, --unit-cell Fix reference unit cell: \"a,b,c,alpha,beta,gamma\"" << std::endl; std::cout << " -r, --refine Geometry refinement algorithm (none|orientation|beam_and_lattice|multi); multi tries all three per image and keeps whichever indexes the most spots" << std::endl; - std::cout << " --refine-geometry[=N] Stills: extra first pass that bundle-adjusts the shared beam/distance/cell from N strongly-indexed frames (default: 200), then re-indexes with it (lifts weak-stills indexing)" << std::endl; + std::cout << " --refine-geometry[=N|off] Stills: extra first pass that bundle-adjusts the shared beam/distance/cell from N strongly-indexed frames (default: 200) then re-indexes (lifts weak-stills indexing). Default ON for stills when a reference cell is given (-C / reference MTZ); =off disables" << std::endl; std::cout << std::endl; std::cout << " Scaling and merging (on by default)" << std::endl; @@ -103,7 +103,7 @@ void print_usage() { std::cout << " --scaling-iterations Number of scaling iterations with no reference data (default: 3)" << std::endl; std::cout << " --scaling-output Output format for scaling results mtz|cif|txt (default: cif)" << std::endl; std::cout << " -z, --reference-mtz Reference MTZ file" << std::endl; - std::cout << " --reference-column