accurate sigma estimation on merging with partiality-weighting

.gitea/workflows/build_and_test.yml

@@ -76,16 +76,27 @@ jobs:
             upload_url: https://gitea.psi.ch/api/packages/mx/debian/pool/noble/nocuda/upload
     steps:
       - uses: actions/checkout@v4
-      - name: Setup build (cmake)
+      - 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
         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
@@ -122,32 +133,6 @@ 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 ci
+            COMMAND npm install
             COMMAND npm run build
             COMMAND npm run redocly
             WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/frontend)

broker/JFJochServices.cpp

@@ -18,12 +18,9 @@ void JFJochServices::Start(const DiffractionExperiment& experiment,
     if (receiver != nullptr) {
         logger.Info("   ... receiver start");
         receiver->Start(experiment, pixel_mask, calibration);
-        {
-            std::shared_lock ul(detector_mutex);
-            if (detector && !experiment.IsUsingInternalPacketGen()) {
-                logger.Info("   ... detector start");
-                detector->Start(experiment);
-            }
+        if (detector && !experiment.IsUsingInternalPacketGen()) {
+            logger.Info("   ... detector start");
+            detector->Start(experiment);
         }
     }
 
@@ -31,38 +28,27 @@ void JFJochServices::Start(const DiffractionExperiment& experiment,
 }
 
 void JFJochServices::Off() {
-    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
+    if (detector) {
+        detector->Deactivate();
+        detector.reset();
     }
 }
 
 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:
-                new_detector = std::make_unique<SLSDetectorWrapper>();
+                detector = std::make_unique<SLSDetectorWrapper>();
                 break;
             case DetectorType::DECTRIS:
-                new_detector = std::make_unique<DectrisDetectorWrapper>();
+                detector = std::make_unique<DectrisDetectorWrapper>();
                 break;
         }
-        new_detector->Initialize(x, receiver->GetNetworkConfig());
-        {
-            std::unique_lock ul(detector_mutex);
-            detector = std::move(new_detector);
-        }
+        detector->Initialize(x, receiver->GetNetworkConfig());
         logger.Info("   ... done");
     }
 }
@@ -76,28 +62,24 @@ JFJochServicesOutput JFJochServices::Stop() {
 
     if (receiver != nullptr) {
         try {
-            {
-                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;
-                    }
+            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();
 
@@ -135,7 +117,6 @@ 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 {
@@ -182,16 +163,15 @@ 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 {
-    std::shared_lock ul(detector_mutex, std::defer_lock);
-    if (ul.try_lock_for(std::chrono::milliseconds(500)) && detector)
+    if (detector)
         return detector->GetStatus();
-    return {};
+    else
+        return {};
 }
 
 std::string JFJochServices::GetPreviewJPEG(const PreviewImageSettings &settings, int64_t image_number) const {
@@ -209,7 +189,6 @@ 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);
 }
@@ -283,7 +262,6 @@ void JFJochServices::ClearImageBuffer() const {
 }
 
 void JFJochServices::LoadDetectorPixelMask(PixelMask &mask) {
-    std::shared_lock ul(detector_mutex);
     if (detector)
         detector->LoadPixelMask(mask);
 }

broker/JFJochServices.h

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

docker/rocky9/Dockerfile

@@ -8,9 +8,6 @@ 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 && \
@@ -68,8 +65,6 @@ 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
 
@@ -190,23 +185,9 @@ 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,25 +78,10 @@ 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; \
@@ -181,6 +166,11 @@ 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/scale_merge/ScaleAndMerge.cpp

@@ -241,6 +241,7 @@ namespace {
         int hkl_slot;
         double I_corr;
         double sigma_corr;
+        double weight;
     };
 
     void scale(const ScaleMergeOptions &opt,
@@ -408,18 +409,16 @@ namespace {
         struct HKLAccum {
             double sum_wI = 0.0;
             double sum_w = 0.0;
-            double sum_wI2 = 0.0;
-            int count = 0;
+            double sum_wsigma2 = 0.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);
             auto &a = accum[co.hkl_slot];
             a.sum_wI += w * co.I_corr;
             a.sum_w += w;
-            a.sum_wI2 += w * co.I_corr * co.I_corr;
-            a.count++;
+            a.sum_wsigma2 += w * w * co.sigma_corr * co.sigma_corr;
         }
 
         for (int h = 0; h < nhkl; ++h) {
@@ -428,15 +427,7 @@ namespace {
                 continue;
 
             out.merged[h].I = a.sum_wI / 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;
-            }
+            out.merged[h].sigma = std::sqrt(a.sum_wsigma2) / a.sum_w;
         }
     }
 
@@ -643,7 +634,8 @@ namespace {
             corr_obs.push_back({
                 o.hkl_slot,
                 static_cast<double>(r.I) / correction,
-                o.sigma / correction
+                o.sigma / correction,
+                1 / correction
             });
         }
     }

make_doc.sh

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

tools/jfjoch_hdf5_test.cpp

@@ -18,48 +18,26 @@ int main(int argc, char **argv) {
 
     RegisterHDF5Filter();
 
-    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;
+    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;
         exit(EXIT_FAILURE);
     }
 
-    HDF5ReadOnlyFile data(argv[optind]);
+    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]);
     HDF5DataSet dataset(data, "/entry/data/data");
     HDF5DataSpace file_space(dataset);
 
@@ -75,8 +53,6 @@ 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;
@@ -88,12 +64,16 @@ int main(int argc, char **argv) {
     logger.Info("Number of images in the original dataset: " + std::to_string(nimages));
 
     // Set file name
-    x.FilePrefix(prefix);
-    x.SetFileWriterFormat(format);
-    x.OverwriteExistingFiles(true);
+    if (std::getenv("HDF5DATASET_WRITE_TEST_PREFIX") == nullptr)
+        x.FilePrefix("writing_test");
+    else
+        x.FilePrefix(std::getenv("HDF5DATASET_WRITE_TEST_PREFIX"));
 
-    if (images_per_file.has_value())
-        x.ImagesPerFile(images_per_file.value());
+    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);
 
     x.ImagesPerTrigger(nimages);
 
@@ -147,6 +127,8 @@ 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());