JFJochServices: Protect detector with mutex to avoid problems while detector is being configured

Updates to CI pipeline

* New docker image for Ubuntu 22.04 with CMake 3.26
* New docker image for Rocky 9 with DIALS 3.27
* New automated test to check for DIALS processing with xia2.ssx

Reviewed-on: #42

.gitea/workflows/build_and_test.yml

@@ -76,27 +76,16 @@ jobs:
             upload_url: https://gitea.psi.ch/api/packages/mx/debian/pool/noble/nocuda/upload
     steps:
       - uses: actions/checkout@v4
-      - name: Install CMake 3.26 on Ubuntu 22.04
-        if: matrix.distro == 'ubuntu2204' || matrix.distro == 'ubuntu2204_nocuda'
-        shell: bash
-        run: |
-          set -euo pipefail
-          apt-get update
-          apt-get install -y wget gpg ca-certificates
-          wget -qO- https://apt.kitware.com/keys/kitware-archive-latest.asc \
-            | gpg --dearmor \
-            | tee /usr/share/keyrings/kitware-archive-keyring.gpg > /dev/null
-          echo "deb [signed-by=/usr/share/keyrings/kitware-archive-keyring.gpg] https://apt.kitware.com/ubuntu/ jammy main" \
-            | tee /etc/apt/sources.list.d/kitware.list > /dev/null
-          apt-get update
-          apt-get install -y cmake=3.26.* cmake-data=3.26.* kitware-archive-keyring
-          cmake --version
-      - name: Build packages
+      - name: Setup build (cmake)
         shell: bash
         run: |
           mkdir -p build
           cd build
           cmake -G Ninja -DJFJOCH_INSTALL_DRIVER_SOURCE=ON -DJFJOCH_VIEWER_BUILD=ON -DCMAKE_BUILD_TYPE=Release ${{ matrix.cmake_flags }} ..
+      - name: Build packages
+        shell: bash
+        run: |
+          cd build
           ninja frontend
           ninja -j16 package
       - name: Upload packages
@@ -133,6 +122,32 @@ jobs:
                  python3 gitea_upload_file.py "$file"
             done
           fi
+  dials-test:
+    name: DIALS processing test
+    runs-on: jfjoch_rocky9
+    steps:
+      - uses: actions/checkout@v4
+      - name: Build processing test
+        shell: bash
+        run: |
+          mkdir -p build
+          cd build
+          cmake -G Ninja -DCMAKE_BUILD_TYPE=Release ..
+          ninja -j16 jfjoch_hdf5_test
+      - name: Generate test data (with virtual data set and 4 linked image files)
+        shell: bash
+        run: |
+          set -euo pipefail
+          mkdir -p dials_test
+          cd dials_test
+          ../build/tools/jfjoch_hdf5_test ../tests/test_data/compression_benchmark.h5 -n100 -f25 -V
+      - name: Run DIALS processing
+        shell: bash
+        run: |
+          source /opt/dials-v3-27-0/dials_env.sh
+          set -euo pipefail
+          cd dials_test
+          xia2.ssx image=writing_test_master.h5 space_group=P43212 unit_cell=78.551,78.551,36.914,90.000,90.000,90.000
   python-client:
     name: Generate python client
     runs-on: jfjoch_rocky8

CMakeLists.txt

@@ -175,7 +175,7 @@ ENDIF()
 
 IF (NOT JFJOCH_WRITER_ONLY)
     ADD_CUSTOM_COMMAND(OUTPUT frontend/dist/index.html
-            COMMAND npm install
+            COMMAND npm ci
             COMMAND npm run build
             COMMAND npm run redocly
             WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/frontend)

broker/JFJochServices.cpp

@@ -18,9 +18,12 @@ void JFJochServices::Start(const DiffractionExperiment& experiment,
     if (receiver != nullptr) {
         logger.Info("   ... receiver start");
         receiver->Start(experiment, pixel_mask, calibration);
-        if (detector && !experiment.IsUsingInternalPacketGen()) {
-            logger.Info("   ... detector start");
-            detector->Start(experiment);
+        {
+            std::shared_lock ul(detector_mutex);
+            if (detector && !experiment.IsUsingInternalPacketGen()) {
+                logger.Info("   ... detector start");
+                detector->Start(experiment);
+            }
         }
     }
 
@@ -28,27 +31,38 @@ void JFJochServices::Start(const DiffractionExperiment& experiment,
 }
 
 void JFJochServices::Off() {
-    if (detector) {
-        detector->Deactivate();
-        detector.reset();
+    std::unique_ptr<DetectorWrapper> old_detector;
+    {
+        std::unique_lock ul(detector_mutex);
+        old_detector = std::move(detector);
+    }
+    if (old_detector) {
+        old_detector->Deactivate();
+        // destroyed here, outside lock
     }
 }
 
 void JFJochServices::On(DiffractionExperiment &x) {
     if (x.IsUsingInternalPacketGen() || (receiver == nullptr))  {
+        std::unique_lock ul(detector_mutex);
         detector.reset();
     } else {
         logger.Info("Detector on");
+        std::unique_ptr<DetectorWrapper> new_detector;
         switch (x.GetDetectorType()) {
             case DetectorType::EIGER:
             case DetectorType::JUNGFRAU:
-                detector = std::make_unique<SLSDetectorWrapper>();
+                new_detector = std::make_unique<SLSDetectorWrapper>();
                 break;
             case DetectorType::DECTRIS:
-                detector = std::make_unique<DectrisDetectorWrapper>();
+                new_detector = std::make_unique<DectrisDetectorWrapper>();
                 break;
         }
-        detector->Initialize(x, receiver->GetNetworkConfig());
+        new_detector->Initialize(x, receiver->GetNetworkConfig());
+        {
+            std::unique_lock ul(detector_mutex);
+            detector = std::move(new_detector);
+        }
         logger.Info("   ... done");
     }
 }
@@ -62,24 +76,28 @@ JFJochServicesOutput JFJochServices::Stop() {
 
     if (receiver != nullptr) {
         try {
-            if (detector) {
-                logger.Info("Wait for detector idle");
-                DetectorState state = detector->GetState();
-                while ((!cannot_stop_detector)
-                       && ((state == DetectorState::WAITING) || (state == DetectorState::BUSY))) {
-                    // check detector state every 5 ms
-                    std::this_thread::sleep_for(std::chrono::milliseconds(5));
-                    state = detector->GetState();
-                }
-                if (state == DetectorState::IDLE) {
-                    logger.Info("   ... detector idle");
-                    receiver->Cancel(true); // cancel silently
-                } else {
-                    logger.Error("   ... detector in error state");
-                    receiver->Cancel(false);
-                    detector_error = true;
+            {
+                std::shared_lock ul(detector_mutex);
+                if (detector) {
+                    logger.Info("Wait for detector idle");
+                    DetectorState state = detector->GetState();
+                    while ((!cannot_stop_detector)
+                           && ((state == DetectorState::WAITING) || (state == DetectorState::BUSY))) {
+                        // check detector state every 5 ms
+                        std::this_thread::sleep_for(std::chrono::milliseconds(5));
+                        state = detector->GetState();
+                    }
+                    if (state == DetectorState::IDLE) {
+                        logger.Info("   ... detector idle");
+                        receiver->Cancel(true); // cancel silently
+                    } else {
+                        logger.Error("   ... detector in error state");
+                        receiver->Cancel(false);
+                        detector_error = true;
+                    }
                 }
             }
+
             logger.Info("Wait for receiver done");
             ret.receiver_output = receiver->Stop();
 
@@ -117,6 +135,7 @@ JFJochServicesOutput JFJochServices::Stop() {
 }
 
 void JFJochServices::Cancel() {
+    std::shared_lock ul(detector_mutex);
     if (detector) {
         // Best effort - if detector cannot be stopped, this is OK, important to still stop receiver
         try {
@@ -163,15 +182,16 @@ void JFJochServices::SetSpotFindingSettings(const SpotFindingSettings &settings)
 }
 
 void JFJochServices::Trigger() {
+    std::shared_lock ul(detector_mutex);
     if (detector && (receiver != nullptr))
         detector->Trigger();
 }
 
 std::optional<DetectorStatus> JFJochServices::GetDetectorStatus() const {
-    if (detector)
+    std::shared_lock ul(detector_mutex, std::defer_lock);
+    if (ul.try_lock_for(std::chrono::milliseconds(500)) && detector)
         return detector->GetStatus();
-    else
-        return {};
+    return {};
 }
 
 std::string JFJochServices::GetPreviewJPEG(const PreviewImageSettings &settings, int64_t image_number) const {
@@ -189,6 +209,7 @@ std::string JFJochServices::GetPreviewTIFF(int64_t image_number) const {
 }
 
 void JFJochServices::ConfigureDetector(const DiffractionExperiment &experiment) {
+    std::unique_lock ul(detector_mutex); // While configuring detector ensure exclusive access (even though pointer is not modified here)
     if (detector)
         detector->Configure(experiment);
 }
@@ -262,6 +283,7 @@ void JFJochServices::ClearImageBuffer() const {
 }
 
 void JFJochServices::LoadDetectorPixelMask(PixelMask &mask) {
+    std::shared_lock ul(detector_mutex);
     if (detector)
         detector->LoadPixelMask(mask);
 }

broker/JFJochServices.h

@@ -4,6 +4,7 @@
 #ifndef JUNGFRAUJOCH_JFJOCHSERVICES_H
 #define JUNGFRAUJOCH_JFJOCHSERVICES_H
 
+#include <shared_mutex>
 #include "../common/DiffractionExperiment.h"
 #include "../jungfrau/JFCalibration.h"
 #include "../common/Logger.h"
@@ -15,6 +16,8 @@ struct JFJochServicesOutput {
 };
 
 class JFJochServices {
+    mutable std::shared_timed_mutex detector_mutex;
+
     JFJochReceiverService *receiver = nullptr;
     std::unique_ptr<DetectorWrapper> detector;
 

docker/rocky9/Dockerfile

@@ -8,6 +8,9 @@ ARG EIGEN_TAG=3.4.0
 ARG LIBJPEG_TURBO_VERSION=3.1.2
 ARG LIBTIFF_VERSION=v4.7.1
 ARG HDF5_TAG="hdf5_1.14.6"
+ARG DIALS_VERSION=3.27.0
+ARG DIALS_TARBALL_URL=https://github.com/dials/dials/releases/download/v3.27.0/dials-v3-27-0-linux-x86_64.tar.xz
+ARG DIALS_PREFIX=/opt
 
 # Update base, enable CRB (RHEL/Rocky 9), and install toolchain + Qt deps
 RUN dnf -y update && \
@@ -65,6 +68,8 @@ RUN dnf -y update && \
       libdrm-devel \
       libglvnd-core-devel \
       libglvnd-devel \
+      glibc-langpack-en \
+      glibc-locale-source \
       freetype-devel && \
     dnf clean all && rm -rf /var/cache/dnf
 
@@ -185,9 +190,23 @@ RUN set -eux; \
     cmake --install .; \
     cd /; rm -rf /tmp/qt-everywhere-src-${QT_VERSION} /tmp/qt-everywhere-src-${QT_VERSION}.tar.xz
 
+# Install DIALS using the official binary tarball workflow
+RUN set -eux; \
+    cd /tmp; \
+    curl -fL -o dials.tar.xz "${DIALS_TARBALL_URL}"; \
+    tar -xJf dials.tar.xz; \
+    cd dials-installer; \
+    ./install --prefix="${DIALS_PREFIX}"; \
+    cd /; \
+    rm -rf /tmp/dials.tar.xz /tmp/dials-*-linux-x86_64
+
 # Make Qt and Eigen discoverable by CMake
 ENV CMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH}:/opt/hdf5-${HDF5_TAG}-static:/opt/qt-${QT_VERSION}-static
 
+ENV LANG=en_US.UTF-8
+ENV LANGUAGE=en_US:en
+ENV LC_ALL=en_US.UTF-8
+
 # Set workdir for your project
 WORKDIR /workspace
 

docker/ubuntu2204/Dockerfile

@@ -22,11 +22,11 @@ RUN set -eux; \
       ca-certificates \
       curl \
       wget \
+      gpg \
       git \
       tar \
       xz-utils \
       build-essential \
-      cmake \
       ninja-build \
       python3 \
       python3-requests \
@@ -78,10 +78,25 @@ RUN set -eux; \
       libassimp-dev \
       libglvnd-dev \
       libfreetype6-dev; \
+    wget -qO- https://apt.kitware.com/keys/kitware-archive-latest.asc \
+      | gpg --dearmor \
+      | tee /usr/share/keyrings/kitware-archive-keyring.gpg > /dev/null; \
+    echo "deb [signed-by=/usr/share/keyrings/kitware-archive-keyring.gpg] https://apt.kitware.com/ubuntu/ jammy main" \
+      | tee /etc/apt/sources.list.d/kitware.list > /dev/null; \
+    apt-get update; \
+    apt-get install -y --no-install-recommends \
+      cmake=3.26.* \
+      cmake-data=3.26.* \
+      kitware-archive-keyring; \
     apt-get -y install gcc-12 g++-12; \
     apt-get clean; \
     rm -rf /var/lib/apt/lists/*
 
+# Use GCC/G++ 12 for builds
+ENV CC=/usr/bin/gcc-12
+ENV CXX=/usr/bin/g++-12
+ENV PATH=/usr/bin:${PATH}
+
 # Build a static OpenSSL
 RUN set -eux; \
     cd /tmp; \
@@ -166,11 +181,6 @@ RUN set -eux; \
 ENV CMAKE_PREFIX_PATH=/opt/libtiff-${LIBTIFF_VERSION}-static:/opt/libjpeg-turbo-${LIBJPEG_TURBO_VERSION}-static
 ENV PKG_CONFIG_PATH=/opt/hdf5-${HDF5_TAG}-static/lib/pkgconfig:/opt/libjpeg-turbo-${LIBJPEG_TURBO_VERSION}-static/lib/pkgconfig:/opt/libtiff-${LIBTIFF_VERSION}-static/lib/pkgconfig:${OPENSSL_ROOT_DIR}/lib/pkgconfig:${OPENSSL_ROOT_DIR}/lib64/pkgconfig
 
-# Use GCC/G++ 12 for builds
-ENV CC=/usr/bin/gcc-12
-ENV CXX=/usr/bin/g++-12
-ENV PATH=/usr/bin:${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; \

image_analysis/lattice_search/LatticeSearch.cpp

@@ -19,10 +19,9 @@ struct NiggliClass {
     gemmi::Mat33 reindex;
     gemmi::CrystalSystem system;
     char centering;
-    int lowest_spacegroup_number;
 };
 
-LatticeSearchResult LatticeSearch(const CrystalLattice &L, double dist_tolerance, double angle_tolerance, int ngc_id) {
+LatticeSearchResult LatticeSearch(const CrystalLattice &L, double dist_tolerance, double angle_tolerance) {
     UnitCell uc = L.GetUnitCell();
     gemmi::UnitCell g_uc(uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
     gemmi::GruberVector g_vec(g_uc, 'P', true);
@@ -44,274 +43,268 @@ LatticeSearchResult LatticeSearch(const CrystalLattice &L, double dist_tolerance
             1, 1,
             true, true, A / 2, A / 2, A / 2, false, false,
             gemmi::Mat33{1, -1, 1, 1, 1, -1, -1, 1, 1},
-            gemmi::CrystalSystem::Cubic, 'F', 196
+            gemmi::CrystalSystem::Cubic, 'F'
         },
         {
             2, 1,
             true, true, D, D, D, false, false,
             {1, -1, 0, -1, 0, 1, -1, -1, -1},
-            gemmi::CrystalSystem::Trigonal, 'R', 146
+            gemmi::CrystalSystem::Trigonal, 'R'
         },
         {
             3, 2,
             true, true, 0, 0, 0, false, false,
             gemmi::Mat33{1, 0, 0, 0, 1, 0, 0, 0, 1},
-            gemmi::CrystalSystem::Cubic, 'P', 195
+            gemmi::CrystalSystem::Cubic, 'P'
         },
         {
             5, 2,
             true, true, -A / 3, -A / 3, -A / 3, false, false,
             gemmi::Mat33{1, 0, 1, 1, 1, 0, 0, 1, 1},
-            gemmi::CrystalSystem::Cubic, 'I', 197
+            gemmi::CrystalSystem::Cubic, 'I'
         },
         {
             4, 2,
             true, true, D, D, D, false, false,
             {1, -1, 0, -1, 0, 1, -1, -1, -1},
-            gemmi::CrystalSystem::Trigonal, 'R', 146
+            gemmi::CrystalSystem::Trigonal, 'R'
         },
         {
             6, 2,
             true, true, D, D, F, true, false,
             {0, 1, 1, 1, 0, 1, 1, 1, 0},
-            gemmi::CrystalSystem::Tetragonal, 'I', 79
+            gemmi::CrystalSystem::Tetragonal, 'I'
         },
         {
             7, 2,
             true, true, D, E, E, true, false,
             {1, 0, 1, 1, 1, 0, 0, 1, 1},
-            gemmi::CrystalSystem::Tetragonal, 'I', 79
+            gemmi::CrystalSystem::Tetragonal, 'I'
         },
         {
             8, 2,
             true, true, D, E, F, true, false,
             {-1, -1, 0, -1, 0, -1, 0, -1, -1},
-            gemmi::CrystalSystem::Orthorhombic, 'I', 23
+            gemmi::CrystalSystem::Orthorhombic, 'I'
         },
         {
             9, 1,
             true, false, A / 2, A / 2, A / 2, false, false,
             {1, 0, 0, -1, 1, 0, -1, -1, -3},
-            gemmi::CrystalSystem::Trigonal, 'R', 146
+            gemmi::CrystalSystem::Trigonal, 'R'
         },
 
         {
             10, 1,
             true, false, D, D, F, false, false,
             {1, 1, 0, 1, -1, 0, 0, 0, -1},
-            gemmi::CrystalSystem::Monoclinic, 'C', 5
+            gemmi::CrystalSystem::Monoclinic, 'C'
         },
         {
             11, 2,
             true, false, 0, 0, 0, false, false,
             {1, 0, 0, 0, 1, 0, 0, 0, 1},
-            gemmi::CrystalSystem::Tetragonal, 'P', 75
+            gemmi::CrystalSystem::Tetragonal, 'P'
         },
         {
             12, 2,
             true, false, 0, 0, -A / 2, false, false,
             {1, 0, 0, 0, 1, 0, 0, 0, 1},
-            gemmi::CrystalSystem::Hexagonal, 'P', 143
+            gemmi::CrystalSystem::Hexagonal, 'P'
         },
         {
             13, 2,
             true, false, 0, 0, F, false, false,
             {1, 1, 0, -1, 1, 0, 0, 0, 1},
-            gemmi::CrystalSystem::Orthorhombic, 'C', 20
+            gemmi::CrystalSystem::Orthorhombic, 'C'
         },
         {
             15, 2,
             true, false, -A / 2, -A / 2, 0, false, false,
             {1, 0, 0, 0, 1, 0, 1, 1, 2},
-            gemmi::CrystalSystem::Tetragonal, 'I', 79
+            gemmi::CrystalSystem::Tetragonal, 'I'
         },
         {
             16, 2,
             true, false, D, D, F, true, false,
             {-1, -1, 0, 1, -1, 0, 1, 1, 2},
-            gemmi::CrystalSystem::Orthorhombic, 'F', 22
+            gemmi::CrystalSystem::Orthorhombic, 'F'
         },
         {
             14, 2,
             true, false, D, D, F, false, false,
             {1, 1, 0, -1, 1, 0, 0, 0, 1},
-            gemmi::CrystalSystem::Monoclinic, 'C', 5
+            gemmi::CrystalSystem::Monoclinic, 'C'
         },
         {
             17, 2,
             true, false, D, E, F, true, false,
             {1, -1, 0, 1, 1, 0, -1, 0, -1},
-            gemmi::CrystalSystem::Monoclinic, 'C', 5
+            gemmi::CrystalSystem::Monoclinic, 'C'
         },
         {
             18, 1,
             false, true, A / 4, A / 2, A / 2, false, false,
             {0, -1, 1, 1, -1, -1, 1, 0, 0},
-            gemmi::CrystalSystem::Tetragonal, 'I', 79
+            gemmi::CrystalSystem::Tetragonal, 'I'
         },
         {
             19, 1,
             false, true, D, A / 2, A / 2, false, false,
             {-1, 0, 0, 0, -1, 1, -1, 1, 1},
-            gemmi::CrystalSystem::Orthorhombic, 'I', 23
+            gemmi::CrystalSystem::Orthorhombic, 'I'
         },
         {
             20, 1,
             false, true, D, E, E, false, false,
             {0, 1, 1, 0, 1, -1, -1, 0, 0},
-            gemmi::CrystalSystem::Monoclinic, 'C', 5
+            gemmi::CrystalSystem::Monoclinic, 'C'
         },
         {
             21, 2,
             false, true, 0, 0, 0, false, false,
             {0, 1, 0, 0, 0, 1, 1, 0, 0},
-            gemmi::CrystalSystem::Tetragonal, 'P', 75
+            gemmi::CrystalSystem::Tetragonal, 'P'
         },
         {
             22, 2,
             false, true, -B / 2, 0, 0, false, false,
             {0, 1, 0, 0, 0, 1, 1, 0, 0},
-            gemmi::CrystalSystem::Hexagonal, 'P', 143
+            gemmi::CrystalSystem::Hexagonal, 'P'
         },
         {
             23, 2,
             false, true, D, 0, 0, false, false,
             {0, 1, 1, 0, -1, 1, 1, 0, 0},
-            gemmi::CrystalSystem::Orthorhombic, 'C', 20
+            gemmi::CrystalSystem::Orthorhombic, 'C'
         },
         {
             24, 2,
             false, true, D, -A / 3, -A / 3, true, false,
             {1, 2, 1, 0, -1, 1, 1, 0, 0},
-            gemmi::CrystalSystem::Trigonal, 'R', 146
+            gemmi::CrystalSystem::Trigonal, 'R'
         },
         {
             25, 2,
             false, true, D, E, E, false, false,
             {0, 1, 1, 0, -1, 1, 1, 0, 0},
-            gemmi::CrystalSystem::Monoclinic, 'C', 5
+            gemmi::CrystalSystem::Monoclinic, 'C'
         },
         {
             26, 1,
             false, false, A / 4, A / 2, A / 2, false, false,
             {1, 0, 0, -1, 2, 0, -1, 0, 2},
-            gemmi::CrystalSystem::Orthorhombic, 'F', 22
+            gemmi::CrystalSystem::Orthorhombic, 'F'
         },
         {
             27, 1,
             false, false, D, A / 2, A / 2, false, false,
             {-1, 2, 0, -1, 0, 0, 0, -1, 1},
-            gemmi::CrystalSystem::Monoclinic, 'C', 5
+            gemmi::CrystalSystem::Monoclinic, 'C'
         },
         {
             28, 1,
             false, false, D, A / 2, 2 * D, false, false,
             {-1, 0, 0, -1, 0, 2, 0, 1, 0},
-            gemmi::CrystalSystem::Monoclinic, 'C', 5
+            gemmi::CrystalSystem::Monoclinic, 'C'
         },
         {
             29, 1,
             false, false, D, 2 * D, A / 2, false, false,
             {1, 0, 0, 1, -2, 0, 0, 0, -1},
-            gemmi::CrystalSystem::Monoclinic, 'C', 5
+            gemmi::CrystalSystem::Monoclinic, 'C'
         },
         {
             30, 1,
             false, false, B / 2, E, 2 * E, false, false,
             {0, 1, 0, 0, 1, -2, -1, 0, 0},
-            gemmi::CrystalSystem::Monoclinic, 'C', 5
+            gemmi::CrystalSystem::Monoclinic, 'C'
         },
 
         {
             31, 1,
             false, false, D, E, F, false, false,
             {1, 0, 0, 0, 1, 0, 0, 0, 1},
-            gemmi::CrystalSystem::Triclinic, 'P', 75
+            gemmi::CrystalSystem::Triclinic, 'P'
         },
 
         {
             32, 2,
             false, false, 0, 0, 0, false, false,
             {1, 0, 0, 0, 1, 0, 0, 0, 1},
-            gemmi::CrystalSystem::Orthorhombic, 'P', 16
+            gemmi::CrystalSystem::Orthorhombic, 'P'
         },
         {
             40, 2,
             false, false, -B / 2, 0, 0, false, false,
             {0, -1, 0, -1, 0, 0, 0, 0, -1},
-            gemmi::CrystalSystem::Orthorhombic, 'C', 20
+            gemmi::CrystalSystem::Orthorhombic, 'C'
         },
         {
             35, 2,
             false, false, D, 0, 0, false, false,
             {0, -1, 0, -1, 0, 0, 0, 0, -1},
-            gemmi::CrystalSystem::Monoclinic, 'P', 3
+            gemmi::CrystalSystem::Monoclinic, 'P'
         },
         {
             36, 2,
             false, false, 0, -A / 2, 0, false, false,
             {1, 0, 0, -1, 0, -2, 0, 1, 0},
-            gemmi::CrystalSystem::Orthorhombic, 'C', 20
+            gemmi::CrystalSystem::Orthorhombic, 'C'
         },
         {
             33, 2,
             false, false, 0, E, 0, false, false,
             {1, 0, 0, 0, 1, 0, 0, 0, 1},
-            gemmi::CrystalSystem::Monoclinic, 'P', 3
+            gemmi::CrystalSystem::Monoclinic, 'P'
         },
         {
             38, 2,
             false, false, 0, 0, -A / 2, false, false,
             {-1, 0, 0, 1, 2, 0, 0, 0, -1},
-            gemmi::CrystalSystem::Orthorhombic, 'C', 20
+            gemmi::CrystalSystem::Orthorhombic, 'C'
         },
         {
             34, 2,
             false, false, 0, 0, F, false, false,
             {-1, 0, 0, 0, 0, -1, 0, -1, 0},
-            gemmi::CrystalSystem::Monoclinic, 'P', 3
+            gemmi::CrystalSystem::Monoclinic, 'P'
         },
         {
             42, 2,
             false, false, -B / 2, -A / 2, 0, false, false,
             {-1, 0, 0, 0, -1, 0, 1, 1, 2},
-            gemmi::CrystalSystem::Orthorhombic, 'I', 23
+            gemmi::CrystalSystem::Orthorhombic, 'I'
         },
         {
             41, 2,
            false, false, -B / 2, E, 0, false, false,
             {0, -1, -2, 0, -1, 0, -1, 0, 0},
-            gemmi::CrystalSystem::Monoclinic, 'C', 5
+            gemmi::CrystalSystem::Monoclinic, 'C'
         },
         {
             37, 2,
             false, false, D, -A / 2, 0, false, false,
             {1, 0, 2, 1, 0, 0, 0, 1, 0},
-            gemmi::CrystalSystem::Monoclinic, 'C', 5
+            gemmi::CrystalSystem::Monoclinic, 'C'
         },
         {
             39, 2,
             false, false, D, 0, -A / 2, false, false,
             {-1, -2, 0, -1, 0, 0, 0, 0, -1},
-            gemmi::CrystalSystem::Monoclinic, 'C', 5
+            gemmi::CrystalSystem::Monoclinic, 'C'
         },
         {
             44, 2,
             false, false, D, E, F, false, false,
             {1, 0, 0, 0, 1, 0, 0, 0, 1},
-            gemmi::CrystalSystem::Triclinic, 'P', 1
+            gemmi::CrystalSystem::Triclinic, 'P'
         }
     };
 
     const auto uc_reduced = L_niggli.GetUnitCell();
 
-    std::vector<NiggliClass> niggli_class_selected;
-    if (ngc_id != -1)
-        niggli_class_selected = {niggli_classes[ngc_id]};
-    else
-        niggli_class_selected = niggli_classes;
-
-    for (const auto &c: niggli_class_selected) {
+    for (const auto &c: niggli_classes) {
         if (c.type == 1 && uc_reduced.beta >= 90 - angle_tolerance )
             continue;
 
@@ -346,7 +339,6 @@ LatticeSearchResult LatticeSearch(const CrystalLattice &L, double dist_tolerance
                 .conventional = L_niggli.Multiply(c.reindex),
                 .system = c.system,
                 .centering = c.centering,
-                .lowest_spacegroup_number = c.lowest_spacegroup_number,
             };
         }
     }
@@ -357,19 +349,5 @@ LatticeSearchResult LatticeSearch(const CrystalLattice &L, double dist_tolerance
         .conventional = L,
         .system = gemmi::CrystalSystem::Triclinic,
         .centering = 'P',
-        .lowest_spacegroup_number = 1,
     };
 }
-
-std::vector<LatticeSearchResult> LatticeCandidates(const CrystalLattice &L, double dist_tolerance, double angle_tolerance) {
-    std::vector<LatticeSearchResult> result;
-    int n_possible_niggli_classes = 44;
-
-    for (int id = 1; id <= n_possible_niggli_classes; id++) {
-        LatticeSearchResult lattice_found = LatticeSearch(L, dist_tolerance, angle_tolerance, id);
-        if (lattice_found.niggli_class != 44)
-            result.push_back(lattice_found);
-    }
-
-    return result;
-}

image_analysis/lattice_search/LatticeSearch.h

@@ -16,11 +16,8 @@ struct LatticeSearchResult {
     CrystalLattice conventional;
     gemmi::CrystalSystem system = gemmi::CrystalSystem::Triclinic;
     char centering = 'P'; // 'P','A','B','C','I','F','R'
-    int lowest_spacegroup_number = 1;
 };
 
-LatticeSearchResult LatticeSearch(const CrystalLattice& L, double dist_tolerance = 0.03, double angle_tolerance_deg = 3, int ngc_id = -1);
-
-std::vector<LatticeSearchResult> LatticeCandidates(const CrystalLattice& L, double dist_tolerance = 0.03, double angle_tolerance_deg = 3);
+LatticeSearchResult LatticeSearch(const CrystalLattice& L, double dist_tolerance = 0.03, double angle_tolerance_deg = 3);
 
 #endif //JFJOCH_LATTICESEARCH_H

image_analysis/scale_merge/ScaleAndMerge.cpp

@@ -7,7 +7,6 @@
 
 #include <thread>
 #include <iostream>
-#include <iomanip>
 #include <algorithm>
 #include <cmath>
 #include <limits>
@@ -18,7 +17,6 @@
 #include <vector>
 
 #include "../../common/ResolutionShells.h"
-#include "../../common/Logger.h"
 
 namespace {
     struct HKLKey {
@@ -237,113 +235,14 @@ namespace {
         int img_id = 0;
         int hkl_slot = -1;
         double sigma = 0.0;
-        mutable bool selected = true;
     };
 
     struct CorrectedObs {
         int hkl_slot;
         double I_corr;
         double sigma_corr;
-        double weight;
     };
 
-    void select_reflections_by_quasi_random(const std::vector<ObsRef> &obs,
-        const ScaleMergeOptions &opt,
-        std::vector<bool> &hkl_selected,
-        Logger &logger) {
-        float stat_d_min = std::numeric_limits<float>::max();
-        float stat_d_max = 0.0f;
-
-        const gemmi::SpaceGroup &sg = *opt.space_group;
-        const gemmi::GroupOps gops = sg.operations();
-        int n_operator = gops.order();
-
-        struct HKLStats {
-            int n = 0;
-            float d = std::numeric_limits<float>::max();
-            int shell_id = 0;
-        };
-        const int nhkl = static_cast<int>(hkl_selected.size());
-        std::vector<HKLStats> per_hkl(nhkl);
-        int reflection_above_cutoff = 0;
-
-        for (const auto &o: obs) {
-            if (o.r->I / o.r->sigma < opt.filter_sigma_cutoff) {
-                o.selected = false;
-                continue;
-            }
-
-            const auto d = o.r->d;
-            reflection_above_cutoff += 1;
-            if (std::isfinite(d) && d > 0.0f) {
-                if (opt.d_min_limit_A > 0.0 && d < static_cast<float>(opt.d_min_limit_A))
-                    continue;
-                stat_d_min = std::min(stat_d_min, d);
-                stat_d_max = std::max(stat_d_max, d);
-                auto &hs = per_hkl[o.hkl_slot];
-                hs.n += 1;
-                hs.d = d;
-            }
-        }
-
-        logger.Info("Reflections of I/sigma > {} : {}", opt.filter_sigma_cutoff, reflection_above_cutoff);
-        if (reflection_above_cutoff < opt.filter_min_per_bin * opt.filter_n_resolution_bins) {
-            logger.Info("No additional selection applied before scaling");
-            return;
-        }
-
-        if (stat_d_min < stat_d_max && stat_d_min > 0.0f) {
-            const float d_min_pad = stat_d_min * 0.999f;
-            const float d_max_pad = stat_d_max * 1.001f;
-            ResolutionShells scaling_shells(d_min_pad, d_max_pad, opt.filter_n_resolution_bins);
-
-            for (int h = 0; h < nhkl; ++h) {
-                const auto d = per_hkl[h].d;
-                if (std::isfinite(d) && d > 0.0f) {
-                    if (opt.d_min_limit_A > 0.0 && d < static_cast<float>(opt.d_min_limit_A))
-                        continue;
-                    auto s = scaling_shells.GetShell(d);
-                    if (s.has_value())
-                        per_hkl[h].shell_id = s.value();
-                }
-            }
-
-            // Accumulators per shell
-            struct ShellAccum {
-                int obs_unique = 0;
-                int obs_total = 0;
-                bool selected = true;
-            };
-            std::vector<ShellAccum> shell_acc(opt.filter_n_resolution_bins);
-
-            for (int h = 0; h < nhkl; ++h) {
-                auto &sa = shell_acc[per_hkl[h].shell_id];
-                if (sa.obs_unique  * n_operator  > opt.filter_min_per_bin * 1.2 || sa.obs_total > opt.filter_max_per_bin)
-                    hkl_selected[h] = false;
-                else
-                    sa.obs_unique += 1;
-                    sa.obs_total += per_hkl[h].n;
-            }
-
-            const auto shell_min_res = scaling_shells.GetShellMinRes();
-
-            logger.Info("| d-mean | n_refl_uni | n_refl_tot | selection |");
-            for (int n=0; n < opt.filter_n_resolution_bins; ++n) {
-                if (shell_acc[n].obs_unique * n_operator < opt.filter_min_per_bin)
-                    shell_acc[n].selected = false;
-                if (shell_min_res[n] <= 0.0f) continue;
-                logger.Info("| {:6.3f} | {:10d} | {:10d} | {:9d} |",
-                    shell_min_res[n], shell_acc[n].obs_unique, shell_acc[n].obs_total, shell_acc[n].selected);
-            }
-
-            for (int h = 0; h < nhkl; ++h) {
-                auto &sa = shell_acc[per_hkl[h].shell_id];
-                if (!sa.selected)
-                    hkl_selected[h] = false;
-            }
-       }
-    }
-
     void scale(const ScaleMergeOptions &opt,
                std::vector<double> &g,
                std::vector<double> &mosaicity,
@@ -351,13 +250,10 @@ namespace {
                const std::vector<uint8_t> &image_slot_used,
                bool rotation_crystallography,
                size_t nhkl,
-               const std::vector<ObsRef> &obs,
-               bool selection,
-               Logger &logger) {
+               const std::vector<ObsRef> &obs) {
         ceres::Problem problem;
 
         std::vector<double> Itrue(nhkl, 0.0);
-        std::vector<bool> hkl_selected(nhkl, true);
 
         // Initialize Itrue from per-HKL median of observed intensities
         {
@@ -381,12 +277,8 @@ namespace {
 
         double wedge = opt.wedge_deg.value_or(0.0);
 
-        if (selection) select_reflections_by_quasi_random(obs, opt, hkl_selected, logger);
-
         std::vector<bool> is_valid_hkl_slot(nhkl, false);
         for (const auto &o: obs) {
-            if (!o.selected) continue;
-            if (!hkl_selected[o.hkl_slot]) continue;
             switch (opt.partiality_model) {
                 case ScaleMergeOptions::PartialityModel::Rotation: {
                     auto *cost = new ceres::AutoDiffCostFunction<IntensityRotResidual, 1, 1, 1, 1, 1>(
@@ -506,9 +398,8 @@ namespace {
         options.function_tolerance = 1e-4;
 
         ceres::Solver::Summary summary;
-        logger.Info("Now start the ceres-solver with residual blocks: {}", problem.NumResidualBlocks());
         ceres::Solve(options, &problem, &summary);
-        logger.Info(summary.FullReport());
+        std::cout << summary.FullReport() << std::endl;
     }
 
     void merge(size_t nhkl, ScaleMergeResult &out, const std::vector<CorrectedObs> &corr_obs) {
@@ -517,19 +408,18 @@ namespace {
         struct HKLAccum {
             double sum_wI = 0.0;
             double sum_w = 0.0;
-            // double sum_wI2 = 0.0;
-            double sum_wsigma2 = 0.0;
-            // int count = 0;
+            double sum_wI2 = 0.0;
+            int count = 0;
         };
         std::vector<HKLAccum> accum(nhkl);
 
         for (const auto &co: corr_obs) {
-            // const double w = 1.0 / (co.sigma_corr * co.sigma_corr);
-            const double w = co.weight / (co.sigma_corr * co.sigma_corr);
+            const double w = 1.0 / (co.sigma_corr * co.sigma_corr);
             auto &a = accum[co.hkl_slot];
             a.sum_wI += w * co.I_corr;
             a.sum_w += w;
-            a.sum_wsigma2 += w * w * co.sigma_corr * co.sigma_corr;
+            a.sum_wI2 += w * co.I_corr * co.I_corr;
+            a.count++;
         }
 
         for (int h = 0; h < nhkl; ++h) {
@@ -538,7 +428,15 @@ namespace {
                 continue;
 
             out.merged[h].I = a.sum_wI / a.sum_w;
-            out.merged[h].sigma = std::sqrt(a.sum_wsigma2) / a.sum_w;
+            const double sigma_stat = std::sqrt(1.0 / a.sum_w);
+
+            if (a.count >= 2) {
+                const double var_internal = a.sum_wI2 - (a.sum_wI * a.sum_wI) / a.sum_w;
+                const double sigma_int = std::sqrt(var_internal / (a.sum_w * (a.count - 1)));
+                out.merged[h].sigma = std::max(sigma_stat, sigma_int);
+            } else {
+                out.merged[h].sigma = sigma_stat;
+            }
         }
     }
 
@@ -745,8 +643,7 @@ namespace {
             corr_obs.push_back({
                 o.hkl_slot,
                 static_cast<double>(r.I) / correction,
-                o.sigma / correction,
-                1 / correction,
+                o.sigma / correction
             });
         }
     }
@@ -797,7 +694,6 @@ namespace {
                 o.img_id = img_id;
                 o.hkl_slot = hkl_slot;
                 o.sigma = sigma;
-                o.selected = true;
                 obs.push_back(o);
             }
         }
@@ -807,14 +703,10 @@ namespace {
 
 ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Reflection> > &observations,
                                                const ScaleMergeOptions &opt) {
-
-    Logger logger("ScaleAndMergeReflections");
-
     if (opt.image_cluster <= 0)
         throw std::invalid_argument("image_cluster must be positive");
 
     const bool rotation_crystallography = opt.wedge_deg.has_value();
-    const bool selection = opt.selection;
 
     size_t nrefl = 0;
     for (const auto &i: observations)
@@ -852,8 +744,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Ref
         }
     }
 
-    logger.Info("Now scale the reflections: {}", nrefl);
-    scale(opt, g, mosaicity, R_sq, image_slot_used, rotation_crystallography, nhkl, obs, selection, logger);
+    scale(opt, g, mosaicity, R_sq, image_slot_used, rotation_crystallography, nhkl, obs);
 
     ScaleMergeResult out;
 

image_analysis/scale_merge/ScaleAndMerge.h

@@ -50,12 +50,6 @@ struct ScaleMergeOptions {
 
     bool refine_wedge = false;
 
-    bool selection = true;
-    double filter_sigma_cutoff = 1.0;
-    int filter_min_per_bin = 10000; // needs optimization
-    int filter_max_per_bin = 80000;
-    int filter_n_resolution_bins = 20;
-
     enum class PartialityModel {Fixed, Rotation, Unity, Still} partiality_model = PartialityModel::Fixed;
 };
 

image_analysis/scale_merge/SearchSpaceGroup.cpp

@@ -8,7 +8,6 @@
 #include <cmath>
 #include <cctype>
 #include <iomanip>
-#include <iostream>
 #include <limits>
 #include <sstream>
 #include <tuple>
@@ -395,16 +394,8 @@ SearchSpaceGroupResult SearchSpaceGroup(
     if (merged.empty())
         return result;
 
-    std::vector<gemmi::SpaceGroup> candidates;
-
     const auto lookup = BuildMergedLookup(merged, opt);
-
-    if (opt.candidate_space_groups.has_value())
-        candidates = opt.candidate_space_groups.value();
-    else
-        candidates = EnumerateCandidateSpaceGroups(opt.crystal_system, opt.centering);
-
-    std::cout << "No. of SG-search candidates: " << candidates.size() << std::endl;
+    const auto candidates = EnumerateCandidateSpaceGroups(opt.crystal_system, opt.centering);
 
     for (const auto& sg : candidates) {
         SpaceGroupCandidateScore score{.space_group = sg};

image_analysis/scale_merge/SearchSpaceGroup.h

@@ -91,8 +91,6 @@ struct SearchSpaceGroupOptions {
     // Acceptance thresholds for absent/allowed <I/sig> ratios.
     double max_absent_to_allowed_i_over_sigma_ratio = 0.20;
     double max_absent_to_allowed_i_over_sigma_ratio_in_any_bin = 0.50;
-
-    std::optional<std::vector<gemmi::SpaceGroup>> candidate_space_groups;
 };
 
 struct SearchSpaceGroupResult {

make_doc.sh

@@ -1,6 +1,8 @@
 #!/bin/bash
 
-python3.12 -m venv tmp_venv/
+set -euo pipefail
+
+python3.11 -m venv tmp_venv/
 source tmp_venv/bin/activate
 
 pip install -r docs/requirements.txt

tools/jfjoch_hdf5_test.cpp

@@ -18,26 +18,48 @@ int main(int argc, char **argv) {
 
     RegisterHDF5Filter();
 
-    if ((argc < 2) || (argc > 4)) {
-        std::cout << "Usage: ./jfjoch_hdf5_test <JF4M hdf5 file> {{<#images>} <rate in Hz>}" << std::endl;
-        std::cout << std::endl;
-        std::cout << "Env. variables:" << std::endl;
-        std::cout << "HDF5DATASET_WRITE_TEST_PREFIX" << std::endl;
-        std::cout << "HDF5MASTER_NEW_FORMAT" << std::endl;
+    int64_t nimages_out = 100;
+
+    std::string prefix = "writing_test";
+    FileWriterFormat format = FileWriterFormat::NXmxLegacy;
+    std::optional<int64_t> images_per_file;
+    std::optional<float> rotation;
+
+    int opt;
+    while ((opt = getopt(argc, argv, "o:n:Vf:R:")) != -1) {
+        switch (opt) {
+            case 'o':
+                prefix = optarg;
+                break;
+            case 'n':
+                nimages_out = atoll(optarg);
+                break;
+            case 'V':
+                format = FileWriterFormat::NXmxVDS;
+                break;
+            case 'R':
+                rotation = atof(optarg);
+                break;
+            case 'f':
+                images_per_file = atoll(optarg);
+                if (images_per_file.value() <= 0) {
+                    std::cerr << "Invalid number of images per file: " << optarg << std::endl;
+                    exit(EXIT_FAILURE);
+                }
+                images_per_file = atoll(optarg);
+                break;
+            default:
+                std::cout << "Usage: ./jfjoch_hdf5_test <JF4M hdf5 file> [-o <prefix>] [-n <num images>] [-V] [-f <num images per file>] [-R<rotation angle>]" << std::endl;
+                exit(EXIT_FAILURE);
+        }
+    }
+
+    if (optind >= argc) {
+        std::cout << "Usage: ./jfjoch_hdf5_test <JF4M hdf5 file> [-o <prefix>] [-n <num images>] [-V] [-f <num images per file>] [-R<rotation angle>]" << std::endl;
         exit(EXIT_FAILURE);
     }
 
-    int64_t nimages_out = 100;
-    double rate = 2200;
-
-    if (argc >= 3)
-        nimages_out = atoi(argv[2]);
-    if (argc >= 4)
-        rate = atof(argv[3]);
-
-    std::chrono::microseconds period_us((rate == 0) ? 0 : (int64_t) (1.0e6 / rate));
-
-    HDF5ReadOnlyFile data(argv[1]);
+    HDF5ReadOnlyFile data(argv[optind]);
     HDF5DataSet dataset(data, "/entry/data/data");
     HDF5DataSpace file_space(dataset);
 
@@ -53,6 +75,8 @@ int main(int argc, char **argv) {
     x.BeamX_pxl(1090).BeamY_pxl(1136).DetectorDistance_mm(75).IncidentEnergy_keV(WVL_1A_IN_KEV);
     x.MaskModuleEdges(true);
     x.MaskChipEdges(true);
+    if (rotation && rotation.value() != 0.0)
+        x.Goniometer(GoniometerAxis("omega", 0, rotation.value(), Coord(-1,0,0), std::nullopt));
 
     if ((file_space.GetDimensions()[1] == 2164) && (file_space.GetDimensions()[2] == 2068)) {
         std::cout << "JF4M with gaps detected (2068 x 2164)" << std::endl;
@@ -64,16 +88,12 @@ int main(int argc, char **argv) {
     logger.Info("Number of images in the original dataset: " + std::to_string(nimages));
 
     // Set file name
-    if (std::getenv("HDF5DATASET_WRITE_TEST_PREFIX") == nullptr)
-        x.FilePrefix("writing_test");
-    else
-        x.FilePrefix(std::getenv("HDF5DATASET_WRITE_TEST_PREFIX"));
+    x.FilePrefix(prefix);
+    x.SetFileWriterFormat(format);
+    x.OverwriteExistingFiles(true);
 
-    if (std::getenv("HDF5MASTER_NEW_FORMAT") != nullptr) {
-        std::cout << "Using new format for HDF5 master file" << std::endl;
-        x.SetFileWriterFormat(FileWriterFormat::NXmxVDS);
-    } else
-        x.SetFileWriterFormat(FileWriterFormat::NXmxLegacy);
+    if (images_per_file.has_value())
+        x.ImagesPerFile(images_per_file.value());
 
     x.ImagesPerTrigger(nimages);
 
@@ -127,8 +147,6 @@ int main(int argc, char **argv) {
 
     size_t total_image_size = 0;
     for (int i = 0; i < nimages_out; i++) {
-        std::this_thread::sleep_until(start_time + i * period_us);
-
         DataMessage message{};
         message.image = CompressedImage(output[i % nimages], x.GetXPixelsNum(), x.GetYPixelsNum(),
                                         x.GetImageMode(), x.GetCompressionAlgorithm());

tools/jfjoch_process.cpp

@@ -2,7 +2,6 @@
 // SPDX-License-Identifier: GPL-3.0-only
 
 #include <iostream>
-#include <iomanip>
 #include <vector>
 #include <string>
 #include <unistd.h>
@@ -28,7 +27,6 @@
 #include "../compression/JFJochCompressor.h"
 #include "../image_analysis/scale_merge/FrenchWilson.h"
 #include "../image_analysis/scale_merge/SearchSpaceGroup.h"
-#include "../image_analysis/lattice_search/LatticeSearch.h"
 #include "../image_analysis/WriteMmcif.h"
 
 void print_usage(Logger &logger) {
@@ -46,15 +44,12 @@ void print_usage(Logger &logger) {
     logger.Info("   -d<num>    High resolution limit for spot finding (default: 1.5)");
     logger.Info("   -D<num>    High resolution limit for scaling/merging (default: 0.0; no limit)");
     logger.Info("   -S<num>    Space group number");
-    logger.Info("   -M[num]    Scale and merge (refine mosaicity) and write scaled.hkl + image.dat, unless indexing rate is below threshold (default: 0.5)");
+    logger.Info("   -M         Scale and merge (refine mosaicity) and write scaled.hkl + image.dat");
     logger.Info("   -P<txt>    Partiality refinement fixed|rot|unity (default: fixed)");
     logger.Info("   -A         Anomalous mode (don't merge Friedel pairs)");
     logger.Info("   -C<cell>   Fix reference unit cell: -C\"a,b,c,alpha,beta,gamma\" (comma-separated, no spaces; quotes optional)");
     logger.Info("   -c<num>    Max spot count (default: 250)");
-    logger.Info("   -W<txt>    HDF5 file with analysis results is written. 'l' or 'light' deactivates image-output");
-    logger.Info("   -T<num>    Noise sigma level for spot finding (default: 3.0)");
-    logger.Info("   -m<num>    Min unique reflections per bin used for scaling. Use all when no value is specified (default: 10000)");
-    logger.Info("   -w         Refine wedge in scaling (default: false)");
+    logger.Info("   -W         HDF5 file with analysis results is written");
 }
 
 void trim_in_place(std::string& t) {
@@ -65,50 +60,6 @@ void trim_in_place(std::string& t) {
     t = t.substr(b, e - b);
 };
 
-void print_statistics(Logger &logger, const MergeStatistics &stats) {
-    logger.Info("");
-    logger.Info("  {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>10s}",
-                "d_min", "N_obs", "N_uniq", "Rmeas", "<I/sig>", "Complete");
-    logger.Info("  {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->10s}",
-                "", "", "", "", "", "");
-    for (const auto &sh: stats.shells) {
-        if (sh.unique_reflections == 0)
-            continue;
-        std::string compl_str = (sh.completeness > 0.0)
-                                    ? fmt::format("{:8.1f}%", sh.completeness * 100.0)
-                                    : "      N/A";
-        logger.Info("  {:8.2f} {:8d} {:8d} {:8.3f}% {:8.1f} {:>10s}",
-                    sh.d_min, sh.total_observations, sh.unique_reflections,
-                    sh.rmeas * 100, sh.mean_i_over_sigma, compl_str);
-    }
-    {
-        const auto &ov = stats.overall;
-        logger.Info("  {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->10s}",
-                    "", "", "", "", "", "");
-        std::string compl_str = (ov.completeness > 0.0)
-                                    ? fmt::format("{:8.1f}%", ov.completeness * 100.0)
-                                    : "      N/A";
-        logger.Info("  {:>8s} {:8d} {:8d} {:8.3f}% {:8.1f} {:>10s}",
-                    "Overall", ov.total_observations, ov.unique_reflections,
-                    ov.rmeas * 100, ov.mean_i_over_sigma, compl_str);
-    }
-    logger.Info("");
-}
-
-void print_unitcell(Logger &logger, const CrystalLattice &lattice, bool print_vector) {
-    auto vec0 = lattice.Vec0();
-    auto vec1 = lattice.Vec1();
-    auto vec2 = lattice.Vec2();
-    auto uc = lattice.GetUnitCell();
-    if (print_vector) {
-        logger.Info("Lattice vec0: {:.3f}, {:.3f}, {:.3f}", vec0.x, vec0.y, vec0.z);
-        logger.Info("Lattice vec1: {:.3f}, {:.3f}, {:.3f}", vec1.x, vec1.y, vec1.z);
-        logger.Info("Lattice vec2: {:.3f}, {:.3f}, {:.3f}", vec2.x, vec2.y, vec2.z);
-    }
-    logger.Info("Lattice:         a={:.2f} b={:.2f} c={:.2f} alpha={:.2f} beta={:.2f} gamma={:.2f}",
-                uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
-}
-
 std::optional<UnitCell> parse_unit_cell_arg(const char* arg)  {
     if (!arg)
         return std::nullopt;
@@ -181,14 +132,8 @@ int main(int argc, char **argv) {
     std::optional<int> space_group_number;
     std::optional<UnitCell> fixed_reference_unit_cell;
     bool write_output = false;
-    bool write_output_noimage = false;
-    bool filtering = true;
-    std::optional<int16_t> filter_min_per_bin;
     std::optional<int64_t> max_spot_count_override;
-    float sigma_spot_finding = 3.0;
-    std::optional<float> merging_threshold;
     IndexingAlgorithmEnum indexing_algorithm = IndexingAlgorithmEnum::Auto;
-    bool refine_wedge = false;
 
     ScaleMergeOptions::PartialityModel partiality_model = ScaleMergeOptions::PartialityModel::Fixed;
 
@@ -201,7 +146,7 @@ int main(int argc, char **argv) {
     }
 
     int opt;
-    while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:M::P:AD:C:T:W:m::w")) != -1) {
+    while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:MP:AD:C:W")) != -1) {
         switch (opt) {
             case 'o':
                 output_prefix = optarg;
@@ -217,10 +162,6 @@ int main(int argc, char **argv) {
                 break;
             case 'W':
                 write_output = true;
-                if (strcmp(optarg, "light") == 0 || strcmp(optarg, "l") == 0) {
-                    write_output_noimage = true;
-                    logger.Warning("Image data will not be saved.");
-                }
                 break;
             case 'v':
                 verbose = true;
@@ -263,14 +204,6 @@ int main(int argc, char **argv) {
             case 'x':
                 refine_beam_center = false;
                 break;
-            case 'm':
-                filtering = false;
-                if (optarg)
-                    filter_min_per_bin = atoi(optarg);
-                break;
-            case 'w':
-                refine_wedge = true;
-                break;
             case 'D':
                 d_min_scale_merge = atof(optarg);
                 logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_min_spot_finding);
@@ -280,15 +213,10 @@ int main(int argc, char **argv) {
                 break;
             case 'M':
                 run_scaling = true;
-                if (optarg) merging_threshold = atof(optarg);
                 break;
             case 'A':
                 anomalous_mode = true;
                 break;
-            case 'T':
-                sigma_spot_finding = atof(optarg);
-                logger.Info("Noise threshold level for spot finding set to {:.2f} sigma", sigma_spot_finding);
-                break;
             case 'C': {
                 auto uc = parse_unit_cell_arg(optarg);
                 if (!uc.has_value()) {
@@ -368,11 +296,8 @@ int main(int argc, char **argv) {
     experiment.OverwriteExistingFiles(true);
     experiment.PolarizationFactor(0.99);
 
-    if (fixed_reference_unit_cell.has_value()) {
+    if (fixed_reference_unit_cell.has_value())
         experiment.SetUnitCell(*fixed_reference_unit_cell);
-    } else {
-        experiment.SetUnitCell({});
-    }
 
     if (max_spot_count_override.has_value()) {
         experiment.MaxSpotCount(max_spot_count_override.value());
@@ -383,8 +308,6 @@ int main(int argc, char **argv) {
     IndexingSettings indexing_settings;
     indexing_settings.Algorithm(indexing_algorithm);
     indexing_settings.RotationIndexing(rotation_indexing);
-    if (rotation_indexing)
-        logger.Info("Rotation indexing is activated.");
     if (rotation_indexing_range.has_value())
         indexing_settings.RotationIndexingMinAngularRange_deg(rotation_indexing_range.value());
 
@@ -394,27 +317,10 @@ int main(int argc, char **argv) {
         indexing_settings.GeomRefinementAlgorithm(GeomRefinementAlgorithmEnum::None);
     experiment.ImportIndexingSettings(indexing_settings);
 
-    switch (experiment.GetIndexingAlgorithm()) {
-        case IndexingAlgorithmEnum::FFBIDX:
-            logger.Info("Indexer used: FFBIDX");
-            break;
-        case IndexingAlgorithmEnum::FFTW:
-            logger.Info("Indexer used: FFTW");
-            break;
-        case IndexingAlgorithmEnum::FFT:
-            logger.Info("Indexer used: FFT (CUDA)");
-            break;
-        case IndexingAlgorithmEnum::None:
-            logger.Warning("Indexer not defined!");
-            return 0;
-        default: ;
-    }
-
     SpotFindingSettings spot_settings;
     spot_settings.enable = true;
     spot_settings.indexing = true;
     spot_settings.high_resolution_limit = d_min_spot_finding;
-    spot_settings.signal_to_noise_threshold = sigma_spot_finding;
     if (d_min_scale_merge > 0)
         spot_settings.high_resolution_limit = d_min_spot_finding;
 
@@ -490,10 +396,6 @@ int main(int argc, char **argv) {
         compressed_buffer.resize(MaxCompressedSize(experiment.GetCompressionAlgorithm(),
                                                    experiment.GetPixelsNum(),
                                                    experiment.GetByteDepthImage()));
-        auto size = compressor.Compress(compressed_buffer.data(),
-                            compressed_buffer.data(),
-                            experiment.GetPixelsNum(),
-                            sizeof(uint8_t));
 
         // Thread-local analysis resources
         MXAnalysisWithoutFPGA analysis(experiment, mapping, pixel_mask, indexer);
@@ -544,11 +446,10 @@ int main(int argc, char **argv) {
             auto image_end_time = std::chrono::high_resolution_clock::now();
             std::chrono::duration<float> image_duration = image_end_time - image_start_time;
 
-            if (!write_output_noimage)
-                size = compressor.Compress(compressed_buffer.data(),
-                                                img->Image().data(),
-                                                experiment.GetPixelsNum(),
-                                                sizeof(int32_t));
+            auto size = compressor.Compress(compressed_buffer.data(),
+                                            img->Image().data(),
+                                            experiment.GetPixelsNum(),
+                                            sizeof(int32_t));
 
             msg.image = CompressedImage(compressed_buffer.data(),
                                         size, experiment.GetXPixelsNum(),
@@ -581,7 +482,7 @@ int main(int argc, char **argv) {
             }
 
             // Progress log
-            if ((current_idx_offset > 0 && (current_idx_offset+1) % 100 == 0) || image_idx == end_image - 1) {
+            if (current_idx_offset > 0 && current_idx_offset % 100 == 0) {
                 std::optional<float> indexing_rate;
                 {
                     std::lock_guard<std::mutex> lock(plots_mutex);
@@ -590,21 +491,11 @@ int main(int argc, char **argv) {
 
                 if (indexing_rate.has_value()) {
                     logger.Info("Processed {} / {} images (indexing rate {:.1f}%)",
-                                current_idx_offset+1, images_to_process,
+                                current_idx_offset, images_to_process,
                                 indexing_rate.value() * 100.0f);
                 } else {
                     logger.Info("Processed {} / {} images (indexing rate N/A)",
-                                current_idx_offset+1, images_to_process);
-                }
-                if (image_idx == end_image - 1) {
-                    if (!merging_threshold.has_value()) merging_threshold = 0.5f;
-                    if (!indexing_rate.has_value()) {
-                        run_scaling = false;
-                    } else if (indexing_rate.value() < merging_threshold) {
-                        run_scaling = false;
-                        logger.Warning("Not proceed to scale and merge with lower indexing rate: {:.1f}%",
-                            indexing_rate.value()*100.0f);
-                    }
+                                current_idx_offset, images_to_process);
                 }
             }
         }
@@ -653,11 +544,10 @@ int main(int argc, char **argv) {
             .niggli_class = rotation_indexer_ret->search_result.niggli_class,
             .crystal_system = rotation_indexer_ret->search_result.system
         };
-        logger.Info("Rotation Indexing found lattice: {} {}",
-            gemmi::crystal_system_str(rotation_indexer_ret->search_result.system),
-            rotation_indexer_ret->search_result.centering);
+        logger.Info("Rotation Indexing found lattice");
     }
 
+    // --- Optional: run scaling (mosaicity refinement) on accumulated reflections ---
     // --- Optional: run scaling (mosaicity refinement) on accumulated reflections ---
     if (run_scaling) {
         logger.Info("Running scaling (mosaicity refinement) ...");
@@ -668,16 +558,6 @@ int main(int argc, char **argv) {
         scale_opts.max_solver_time_s = 240.0; // generous cutoff for now
         scale_opts.merge_friedel = !anomalous_mode;
         scale_opts.d_min_limit_A = d_min_scale_merge.value_or(0.0);
-        scale_opts.refine_wedge = refine_wedge;
-        // scale_opts.min_partiality_for_merge = 0.4;
-        if (filter_min_per_bin.has_value()) {
-            scale_opts.selection = true;
-            scale_opts.filter_min_per_bin = filter_min_per_bin.value();
-        } else {
-            scale_opts.selection = filtering;
-        }
-        if (rotation_indexing)
-            scale_opts.filter_min_per_bin = scale_opts.filter_min_per_bin * 0.25;
 
         const bool fixed_space_group = space_group || experiment.GetGemmiSpaceGroup().has_value();
 
@@ -685,48 +565,27 @@ int main(int argc, char **argv) {
             scale_opts.space_group = *space_group;
         else
             scale_opts.space_group = experiment.GetGemmiSpaceGroup();
-        if (scale_opts.space_group->number == 0) scale_opts.space_group = *gemmi::find_spacegroup_by_number(1);
 
-        logger.Info("Starting SG-no.: {}", scale_opts.space_group->number);
         auto scale_start = std::chrono::steady_clock::now();
         auto scale_result = indexer.ScaleRotationData(scale_opts);
         auto scale_end = std::chrono::steady_clock::now();
         double scale_time = std::chrono::duration<double>(scale_end - scale_start).count();
 
-        if (scale_result) print_statistics(logger, scale_result->statistics);
-//        if (scale_opts.wedge_deg.has_value()) logger.Info("Refined wedge: {:.3f}", scale_opts.wedge_deg.value());
-
         if (scale_result && !fixed_space_group) {
             logger.Info("Searching for space group from P1-merged reflections ...");
 
-            std::vector<gemmi::SpaceGroup> candidates_preset;
-            if (rotation_indexer_ret.has_value()) {
-                std::vector<LatticeSearchResult> lattice_candidates = LatticeCandidates(rotation_indexer_ret->lattice);
-                for (const auto &lc: lattice_candidates) {
-                    logger.Info("Possible lattice : {} {} ({})", gemmi::crystal_system_str(lc.system), lc.centering, lc.lowest_spacegroup_number);
-                    print_unitcell(logger, lc.conventional, false);
-                    candidates_preset.push_back(*gemmi::find_spacegroup_by_number(lc.lowest_spacegroup_number));
-                }
-                candidates_preset.push_back(*gemmi::find_spacegroup_by_number(1));
-            }
-
             SearchSpaceGroupOptions sg_opts;
             sg_opts.crystal_system.reset();
             sg_opts.centering = '\0';
             sg_opts.merge_friedel = !anomalous_mode;
             sg_opts.d_min_limit_A = d_min_scale_merge.value_or(0.0);
-            sg_opts.min_i_over_sigma = 4.0; // 0.0; follows the DIALS's default
+            sg_opts.min_i_over_sigma = 0.0;
             sg_opts.min_operator_cc = 0.80;
             sg_opts.min_pairs_per_operator = 20;
             sg_opts.min_total_compared = 100;
             sg_opts.test_systematic_absences = true;
-            if (~candidates_preset.empty())
-                sg_opts.candidate_space_groups = candidates_preset;
 
-            auto sg_search_start = std::chrono::steady_clock::now();
             const auto sg_search = SearchSpaceGroup(scale_result->merged, sg_opts);
-            auto sg_search_end = std::chrono::steady_clock::now();
-            double sg_search_time = std::chrono::duration<double>(sg_search_end - sg_search_start).count();
 
             logger.Info("");
             {
@@ -740,24 +599,17 @@ int main(int argc, char **argv) {
             logger.Info("");
 
             if (sg_search.best_space_group.has_value()) {
-                logger.Info("SG-search wall-clock time: {:.2f} s", sg_search_time);
-                if (sg_search.best_space_group->number != 0) {
-                    if (sg_search.best_space_group->number != scale_opts.space_group->number) {
-                        logger.Info("Re-running scaling in detected space group {}", sg_search.best_space_group->short_name());
+                logger.Info("Re-running scaling in detected space group {}", sg_search.best_space_group->short_name());
 
-                        scale_opts.space_group = *sg_search.best_space_group;
+                scale_opts.space_group = *sg_search.best_space_group;
 
-                        auto rescale_start = std::chrono::steady_clock::now();
-                        auto refined_scale_result = indexer.ScaleRotationData(scale_opts);
-                        auto rescale_end = std::chrono::steady_clock::now();
+                auto rescale_start = std::chrono::steady_clock::now();
+                auto refined_scale_result = indexer.ScaleRotationData(scale_opts);
+                auto rescale_end = std::chrono::steady_clock::now();
 
-                        if (refined_scale_result) {
-                            scale_result = std::move(refined_scale_result);
-                            scale_time += std::chrono::duration<double>(rescale_end - rescale_start).count();
-                        }
-                    } else {
-                        logger.Info("No space group update indicated.");
-                    }
+                if (refined_scale_result) {
+                    scale_result = std::move(refined_scale_result);
+                    scale_time += std::chrono::duration<double>(rescale_end - rescale_start).count();
                 }
             } else {
                 logger.Warning("No space group accepted; keeping P1-merged result");
@@ -771,7 +623,36 @@ int main(int argc, char **argv) {
                         scale_time, scale_result->merged.size());
 
             // Print resolution-shell statistics table
-            print_statistics(logger, scale_result->statistics);
+            {
+                const auto &stats = scale_result->statistics;
+                logger.Info("");
+                logger.Info("  {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>10s}",
+                            "d_min", "N_obs", "N_uniq", "Rmeas", "<I/sig>", "Complete");
+                logger.Info("  {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->10s}",
+                            "", "", "", "", "", "");
+                for (const auto &sh: stats.shells) {
+                    if (sh.unique_reflections == 0)
+                        continue;
+                    std::string compl_str = (sh.completeness > 0.0)
+                                                ? fmt::format("{:8.1f}%", sh.completeness * 100.0)
+                                                : "      N/A";
+                    logger.Info("  {:8.2f} {:8d} {:8d} {:8.3f}% {:8.1f} {:>10s}",
+                                sh.d_min, sh.total_observations, sh.unique_reflections,
+                                sh.rmeas * 100, sh.mean_i_over_sigma, compl_str);
+                }
+                {
+                    const auto &ov = stats.overall;
+                    logger.Info("  {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->10s}",
+                                "", "", "", "", "", "");
+                    std::string compl_str = (ov.completeness > 0.0)
+                                                ? fmt::format("{:8.1f}%", ov.completeness * 100.0)
+                                                : "      N/A";
+                    logger.Info("  {:>8s} {:8d} {:8d} {:8.3f}% {:8.1f} {:>10s}",
+                                "Overall", ov.total_observations, ov.unique_reflections,
+                                ov.rmeas * 100, ov.mean_i_over_sigma, compl_str);
+                }
+                logger.Info("");
+            }
 
             {
                 const std::string img_path = output_prefix + "_image.dat";
@@ -817,8 +698,6 @@ int main(int argc, char **argv) {
                         cif_meta.unit_cell = rotation_indexer_ret->lattice.GetUnitCell();
                     } else if (experiment.GetUnitCell().has_value()) {
                         cif_meta.unit_cell = experiment.GetUnitCell().value();
-                    } else {
-                        logger.Warning("No UnitCell output");
                     }
 
                     if (scale_opts.space_group.has_value()) {
@@ -866,11 +745,12 @@ int main(int argc, char **argv) {
     double frame_rate = static_cast<double>(images_to_process) / processing_time;
 
     logger.Info("");
-    logger.Info("Processing time (excl. scaling): {:8.2f} s", processing_time);
-    logger.Info("Frame rate:                      {:8.2f} Hz", frame_rate);
-    logger.Info("Total throughput:                {:8.2f} MB/s", throughput_MBs);
+    logger.Info("Processing time: {:.2f} s", processing_time);
+    logger.Info("Frame rate:      {:.2f} Hz", frame_rate);
+    logger.Info("Total throughput:{:.2f} MB/s", throughput_MBs);
+
     // Print extended stats similar to Receiver
-    if (end_msg.indexing_rate.has_value()) {
+    if (!end_msg.indexing_rate.has_value()) {
         logger.Info("Indexing rate:   {:.2f}%", end_msg.indexing_rate.value() * 100.0);
     }
 
@@ -891,6 +771,14 @@ int main(int argc, char **argv) {
             latt = latt.Multiply(rot);
         }
 
-        print_unitcell(logger, rotation_indexer_ret->lattice, true);
+        auto vec0 = rotation_indexer_ret->lattice.Vec0();
+        auto vec1 = rotation_indexer_ret->lattice.Vec1();
+        auto vec2 = rotation_indexer_ret->lattice.Vec2();
+        auto uc = rotation_indexer_ret->lattice.GetUnitCell();
+        logger.Info("Lattice vec0: {:.3f}, {:.3f}, {:.3f}", vec0.x, vec0.y, vec0.z);
+        logger.Info("Lattice vec1: {:.3f}, {:.3f}, {:.3f}", vec1.x, vec1.y, vec1.z);
+        logger.Info("Lattice vec2: {:.3f}, {:.3f}, {:.3f}", vec2.x, vec2.y, vec2.z);
+        logger.Info("Lattice:         a={:.2f} b={:.2f} c={:.2f} alpha={:.2f} beta={:.2f} gamma={:.2f}",
+                    uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
     }
 }