From d2159bde92fa00900d6ca35e469e6d0a57c49c32 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 26 Jun 2026 11:42:44 +0200 Subject: [PATCH 01/66] jfjoch_viewer: Grey out "Analyze dataset" on a live HTTP connection Dataset re-processing reads a stored HDF5 file, so it is unavailable for the live HTTP stream. Disable the "Analyze dataset" hero button while a source is connected (with an explanatory tooltip) instead of letting the user click through to a "open a file first" dialog afterwards. Co-Authored-By: Claude Opus 4.8 (1M context) --- viewer/JFJochViewerWindow.cpp | 3 +++ viewer/widgets/JFJochViewerSettingsDock.cpp | 17 ++++++++++++----- viewer/widgets/JFJochViewerSettingsDock.h | 6 ++++++ 3 files changed, 21 insertions(+), 5 deletions(-) diff --git a/viewer/JFJochViewerWindow.cpp b/viewer/JFJochViewerWindow.cpp index b0f9ee4b..a6547f5f 100644 --- a/viewer/JFJochViewerWindow.cpp +++ b/viewer/JFJochViewerWindow.cpp @@ -481,6 +481,9 @@ JFJochViewerWindow::JFJochViewerWindow(QWidget *parent, bool dbus, const QString reading_worker, &JFJochImageReadingWorker::ReanalyzeImages); connect(settingsPanel, &JFJochViewerSettingsDock::analyzeDataset, processingJobsWindow, &JFJochProcessingJobsWindow::newJob); + // Grey out "Analyze dataset" on a live HTTP source (re-processing needs a stored file). + connect(reading_worker, &JFJochImageReadingWorker::httpConnectionChanged, + settingsPanel, &JFJochViewerSettingsDock::setHttpConnection); connect(settingsPanel, &JFJochViewerSettingsDock::experimentChanged, reading_worker, &JFJochImageReadingWorker::UpdateDataset); connect(settingsPanel, &JFJochViewerSettingsDock::findBeamCenter, diff --git a/viewer/widgets/JFJochViewerSettingsDock.cpp b/viewer/widgets/JFJochViewerSettingsDock.cpp index da00c622..3c8a5c9f 100644 --- a/viewer/widgets/JFJochViewerSettingsDock.cpp +++ b/viewer/widgets/JFJochViewerSettingsDock.cpp @@ -71,16 +71,16 @@ JFJochViewerSettingsDock::JFJochViewerSettingsDock(const SpotFindingSettings &sp analyzeImageBtn->setStyleSheet(heroStyle); analyzeImageBtn->setToolTip("Re-analyse the current image now, and keep re-analysing on every" " image / settings change while active"); - auto *analyzeDatasetBtn = new QPushButton(FramesIcon(3), " Analyze dataset", this); - analyzeDatasetBtn->setStyleSheet(heroStyle); - analyzeDatasetBtn->setToolTip("Process the whole dataset (MX or azimuthal, per the toggle below)"); + analyzeDatasetBtn_ = new QPushButton(FramesIcon(3), " Analyze dataset", this); + analyzeDatasetBtn_->setStyleSheet(heroStyle); + analyzeDatasetBtn_->setToolTip("Process the whole dataset (MX or azimuthal, per the toggle below)"); auto *analyzeRow = new QHBoxLayout(); analyzeRow->addWidget(analyzeImageBtn); - analyzeRow->addWidget(analyzeDatasetBtn); + analyzeRow->addWidget(analyzeDatasetBtn_); layout->addLayout(analyzeRow); layout->addSpacing(10); connect(analyzeImageBtn, &QPushButton::toggled, this, &JFJochViewerSettingsDock::reanalyzeImage); - connect(analyzeDatasetBtn, &QPushButton::clicked, this, [this] { emit analyzeDataset(azint_mode_); }); + connect(analyzeDatasetBtn_, &QPushButton::clicked, this, [this] { emit analyzeDataset(azint_mode_); }); // Segmented MX / AzInt toggle: the two communities pick their page; pages never share a screen. auto *mxButton = new QPushButton("MX", this); @@ -115,6 +115,13 @@ JFJochViewerSettingsDock::JFJochViewerSettingsDock(const SpotFindingSettings &sp layout->addStretch(); } +void JFJochViewerSettingsDock::setHttpConnection(bool connected, QString) { + analyzeDatasetBtn_->setEnabled(!connected); + analyzeDatasetBtn_->setToolTip(connected + ? "Dataset re-processing is only available for an open file, not a live HTTP stream" + : "Process the whole dataset (MX or azimuthal, per the toggle below)"); +} + QWidget *JFJochViewerSettingsDock::BuildGeometrySection() { auto *section = new CollapsibleSection("Geometry", this); auto *geom = new QFormLayout(); diff --git a/viewer/widgets/JFJochViewerSettingsDock.h b/viewer/widgets/JFJochViewerSettingsDock.h index 27e22e21..d4e36b2a 100644 --- a/viewer/widgets/JFJochViewerSettingsDock.h +++ b/viewer/widgets/JFJochViewerSettingsDock.h @@ -40,6 +40,9 @@ public slots: void datasetLoaded(std::shared_ptr dataset); void loadImage(std::shared_ptr image); void referenceLoaded(ReferenceMtzInfo info); // worker reports a (un)loaded reference MTZ + // Dataset re-processing reads a stored HDF5 file, so it is unavailable on a live HTTP stream: + // grey out "Analyze dataset" while connected rather than failing with a dialog afterwards. + void setHttpConnection(bool connected, QString addr); signals: void spotFindingChanged(const SpotFindingSettings &spot, const IndexingSettings &indexing, int64_t max_spots); @@ -64,6 +67,9 @@ private: int64_t max_spots_ = 1000; bool azint_mode_ = false; // false = MX page, true = AzInt page (drives "Analyze dataset") + // "Analyze dataset" hero button, disabled while a live HTTP source is connected. + QPushButton *analyzeDatasetBtn_ = nullptr; + // Indexing-algorithm combo + the description line under it (kept so the Auto resolution refreshes). QComboBox *algo_ = nullptr; QLabel *algoDesc_ = nullptr; -- 2.54.0 From 786af96b3be3b3da966d78a8011c5950ad870798 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 26 Jun 2026 16:11:28 +0200 Subject: [PATCH 02/66] SearchSpaceGroup: POINTLESS-style rewrite, pipeline integration, twinning test Space-group search (image_analysis/scale_merge/SearchSpaceGroup): - Two-stage POINTLESS-style determination. Stage A scores each distinct rotation operator once (was once per candidate space group, ~34x faster on lysozyme: ~26s -> <1s) and picks the largest point group all of whose operators confirm. Stage B picks the maximal space group whose predicted absences are confirmed weak, fixing the prototype's default to the symmorphic group (it returned P422 instead of P4(3)2(1)2). Enantiomorphic / origin-ambiguous pairs (P4(1) vs P4(3), I222 vs I2(1)2(1)2(1)) are reported as indistinguishable. - Constrain candidates to subgroups of the lattice (metric) holohedry and weigh centering only P-vs-metric, fed from rotation indexing's LatticeSearch result. Integration / pipeline: - With no user-fixed space group, predict in P (IndexAndRefine) so the centering-absent reflections are integrated and the search can confirm/deny centering (catching pseudo-centering / a missed superstructure) instead of trusting the metric; a user-fixed group still rejects absences in integration. - JFJochProcess: scale+merge in P1 -> determine the space group -> set it and re-scale+merge in it (statistics then come out in the right symmetry) -> write it to /entry/sample/space_group_number (new EndMessage.space_group_number, preferred by NXmx::Sample). jfjoch_scale no longer searches; it consumes the file's space group (and no longer clobbers it with an empty -S). Twinning (new image_analysis/scale_merge/TwinningAnalysis): Padilla-Yeates L-test (<|L|>, ; acentric-only, positive intensities so L is bounded) plus a shell-normalised /^2 second moment and a twin-fraction estimate. Reported after the final merge in jfjoch_process and jfjoch_scale, and surfaced in the jfjoch_viewer merge-statistics window with a red outline when twinning is suspected. Co-Authored-By: Claude Opus 4.8 (1M context) --- common/JFJochMessages.h | 3 + image_analysis/IndexAndRefine.cpp | 6 +- image_analysis/scale_merge/CMakeLists.txt | 2 + .../scale_merge/SearchSpaceGroup.cpp | 834 +++++++----------- image_analysis/scale_merge/SearchSpaceGroup.h | 128 +-- .../scale_merge/TwinningAnalysis.cpp | 156 ++++ image_analysis/scale_merge/TwinningAnalysis.h | 44 + process/JFJochProcess.cpp | 153 ++-- process/JFJochProcess.h | 8 + tests/SearchSpaceGroupTest.cpp | 25 +- tools/jfjoch_scale.cpp | 35 +- .../image_viewer/JFJochDiffractionImage.cpp | 8 + viewer/windows/JFJochMergeStatsWindow.cpp | 24 +- viewer/windows/JFJochMergeStatsWindow.h | 4 +- viewer/windows/JFJochProcessingJobsWindow.cpp | 4 +- viewer/windows/JFJochProcessingJobsWindow.h | 1 + writer/HDF5NXmx.cpp | 10 +- 17 files changed, 789 insertions(+), 656 deletions(-) create mode 100644 image_analysis/scale_merge/TwinningAnalysis.cpp create mode 100644 image_analysis/scale_merge/TwinningAnalysis.h diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index 11a29806..2e04885e 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -338,6 +338,9 @@ struct EndMessage { std::optional rotation_lattice; std::vector rotation_extra_lattices; std::optional unit_cell; + // Space group determined by the offline analysis (overrides the start message when writing the + // master, since it is only known after merging). + std::optional space_group_number; // Vectors with end result: std::vector data_collection_efficiency; diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 64a394f0..c1b0ac7e 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -312,7 +312,11 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(), .ewald_dist_cutoff = ewald_dist_cutoff, .max_hkl = 100, - .centering = outcome.symmetry.centering, + // Centering is a hypothesis to confirm, not assume: with no user-fixed space group, predict + // in P so the centering-absent reflections are integrated and the space-group search can + // confirm or disprove centering (and catch a missed superstructure). A user-fixed space + // group is trusted, so reject its absences here. + .centering = experiment.GetGemmiSpaceGroup().has_value() ? outcome.symmetry.centering : 'P', .wedge_deg = std::fabs(wedge_deg), .mosaicity_deg = std::fabs(mos_deg), // FWHM -> sigma; 0 when monochromatic, leaving the prediction unchanged. diff --git a/image_analysis/scale_merge/CMakeLists.txt b/image_analysis/scale_merge/CMakeLists.txt index 190040c7..ffa3020a 100644 --- a/image_analysis/scale_merge/CMakeLists.txt +++ b/image_analysis/scale_merge/CMakeLists.txt @@ -1,6 +1,8 @@ ADD_LIBRARY(JFJochScaleMerge SearchSpaceGroup.cpp SearchSpaceGroup.h + TwinningAnalysis.cpp + TwinningAnalysis.h Merge.cpp Merge.h ScaleOnTheFly.cpp diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index add88fcb..de71f6a7 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -6,382 +6,156 @@ #include #include #include -#include +#include #include #include +#include #include #include #include -#include #include namespace { - struct MergedHKLKey { - int h = 0; - int k = 0; - int l = 0; - bool plus = true; - - bool operator==(const MergedHKLKey& o) const noexcept { - return h == o.h && k == o.k && l == o.l && plus == o.plus; - } + // A merged reflection, folded onto the +/- Friedel-equivalent it represents, used as a + // hash key to match symmetry-related reflections. + struct HKLKey { + int h = 0, k = 0, l = 0; + bool operator==(const HKLKey& o) const noexcept { return h == o.h && k == o.k && l == o.l; } }; - struct MergedHKLKeyHash { - size_t operator()(const MergedHKLKey& key) const noexcept { + struct HKLKeyHash { + size_t operator()(const HKLKey& key) const noexcept { auto mix = [](uint64_t x) { - x ^= x >> 33; - x *= 0xff51afd7ed558ccdULL; - x ^= x >> 33; - x *= 0xc4ceb9fe1a85ec53ULL; - x ^= x >> 33; - return x; + x ^= x >> 33; x *= 0xff51afd7ed558ccdULL; + x ^= x >> 33; x *= 0xc4ceb9fe1a85ec53ULL; + x ^= x >> 33; return x; }; - return static_cast( - mix(static_cast(key.h)) ^ - (mix(static_cast(key.k)) << 1) ^ - (mix(static_cast(key.l)) << 2) ^ - (mix(static_cast(key.plus ? 1 : 0)) << 3)); + return static_cast(mix(static_cast(key.h)) ^ + (mix(static_cast(key.k)) << 1) ^ + (mix(static_cast(key.l)) << 2)); } }; - bool IsIdentityOp(const gemmi::Op& op) { - return op.rot == gemmi::Op::identity().rot; - } - - bool IsParityChanging(const gemmi::Op& op) { - return op.det_rot() < 0; - } - - MergedHKLKey CanonicalizeMergedHKL(int h, int k, int l, bool merge_friedel) { - MergedHKLKey key{h, k, l, true}; - if (!merge_friedel) - return key; - - const std::tuple pos{h, k, l}; - const std::tuple neg{-h, -k, -l}; - if (neg < pos) { - key.h = -h; - key.k = -k; - key.l = -l; - key.plus = false; - } - return key; - } - - bool ReflectionPassesFilters(const MergedReflection& r, const SearchSpaceGroupOptions& opt) { - if (!std::isfinite(r.I)) - return false; - if (!std::isfinite(r.sigma) || r.sigma <= 0.0) - return false; - if (!std::isfinite(r.d) || r.d <= 0.0) - return false; - if (opt.d_min_limit_A > 0.0 && r.d < opt.d_min_limit_A) - return false; - if (opt.min_i_over_sigma > 0.0 && r.I / r.sigma < opt.min_i_over_sigma) - return false; - return true; - } - - std::unordered_map - BuildMergedLookup(const std::vector& merged, const SearchSpaceGroupOptions& opt) { - std::unordered_map out; - out.reserve(merged.size()); - - for (const auto& r : merged) { - if (!ReflectionPassesFilters(r, opt)) - continue; - const auto key = CanonicalizeMergedHKL(r.h, r.k, r.l, opt.merge_friedel); - out.emplace(key, &r); - } - - return out; + HKLKey Canonicalize(int h, int k, int l, bool merge_friedel) { + if (merge_friedel && std::make_tuple(-h, -k, -l) < std::make_tuple(h, k, l)) + return {-h, -k, -l}; + return {h, k, l}; } double PearsonCC(const std::vector& x, const std::vector& y) { - if (x.size() != y.size() || x.size() < 2) + if (x.size() < 2) return std::numeric_limits::quiet_NaN(); - double sx = 0.0, sy = 0.0; - double sxx = 0.0, syy = 0.0, sxy = 0.0; - + double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0; for (size_t i = 0; i < x.size(); ++i) { - sx += x[i]; - sy += y[i]; - sxx += x[i] * x[i]; - syy += y[i] * y[i]; - sxy += x[i] * y[i]; + sx += x[i]; sy += y[i]; + sxx += x[i] * x[i]; syy += y[i] * y[i]; sxy += x[i] * y[i]; } - const double n = static_cast(x.size()); - const double cov = sxy - sx * sy / n; const double vx = sxx - sx * sx / n; const double vy = syy - sy * sy / n; - - if (vx <= 0.0 || vy <= 0.0) + if (vx <= 0 || vy <= 0) return std::numeric_limits::quiet_NaN(); - - return cov / std::sqrt(vx * vy); + return (sxy - sx * sy / n) / std::sqrt(vx * vy); } - std::optional ScoreOperator( - const gemmi::Op& op, - const std::vector& merged, - const std::unordered_map& lookup, - const SearchSpaceGroupOptions& opt) { - - std::vector x; - std::vector y; - x.reserve(merged.size() / 2); - y.reserve(merged.size() / 2); - - std::unordered_set used; - used.reserve(merged.size() / 2); - - for (const auto& r : merged) { - if (!ReflectionPassesFilters(r, opt)) - continue; - - const auto key1 = CanonicalizeMergedHKL(r.h, r.k, r.l, opt.merge_friedel); - if (used.find(key1) != used.end()) - continue; - - const gemmi::Op::Miller hkl{{r.h, r.k, r.l}}; - const auto hkl2 = op.apply_to_hkl(hkl); - const auto key2 = CanonicalizeMergedHKL(hkl2[0], hkl2[1], hkl2[2], opt.merge_friedel); - - if (key1.h == key2.h && key1.k == key2.k && key1.l == key2.l && key1.plus == key2.plus) - continue; - - auto it = lookup.find(key2); - if (it == lookup.end()) - continue; - - const auto* mate = it->second; - if (mate == nullptr) - continue; - - x.push_back(r.I); - y.push_back(mate->I); - - used.insert(key1); - used.insert(key2); + // A reflection is extinct from lattice centering alone (independent of any screw/glide) when a + // centering translation makes its structure factor cancel. Mirrors the centering half of + // gemmi::GroupOps::is_systematically_absent, so screw absences can be judged separately. + bool CenteringAbsent(const gemmi::GroupOps& gops, const gemmi::Op::Miller& hkl) { + for (size_t i = 1; i < gops.cen_ops.size(); ++i) { + const auto& t = gops.cen_ops[i]; + if ((t[0] * hkl[0] + t[1] * hkl[1] + t[2] * hkl[2]) % gemmi::Op::DEN != 0) + return true; } + return false; + } - if (x.empty()) - return std::nullopt; - - SpaceGroupOperatorScore out; - out.op_triplet_hkl = op.as_hkl().triplet('h'); - out.compared = static_cast(x.size()); - out.cc = PearsonCC(x, y); - out.accepted = (out.compared >= opt.min_pairs_per_operator && - std::isfinite(out.cc) && - out.cc >= opt.min_operator_cc); + std::array RotKey(const gemmi::Op& op) { + std::array out{}; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + out[i * 3 + j] = op.rot[i][j]; return out; } - struct ResolutionBinAccumulator { - double d_min_A = std::numeric_limits::infinity(); - double d_max_A = 0.0; - double absent_sum = 0.0; - int absent_count = 0; - double allowed_sum = 0.0; - int allowed_count = 0; + // The rotation part of a space group in the reference setting (identity included), as a + // sorted list of matrices - the key that groups space groups into a candidate point group. + // It must be the rotation SET, not gemmi's PointGroup enum: P321 and P312 are both "32" yet + // have their 2-folds along different directions, and only the matrices tell them apart. + using RotationSet = std::vector>; + + RotationSet RotationSetOf(const gemmi::SpaceGroup& sg) { + RotationSet out; + for (const auto& op : sg.operations().derive_symmorphic().sym_ops) + out.push_back(RotKey(op)); + std::sort(out.begin(), out.end()); + return out; + } + + // Proper rotations of a crystal system's holohedry (the highest lattice symmetry it can host), + // in the reference setting. Any candidate point group must be a subgroup of this. + RotationSet HolohedryRotationSet(gemmi::CrystalSystem system) { + int number = 0; + switch (system) { + case gemmi::CrystalSystem::Triclinic: number = 1; break; // P1 + case gemmi::CrystalSystem::Monoclinic: number = 3; break; // P2 (unique axis b) + case gemmi::CrystalSystem::Orthorhombic: number = 16; break; // P222 + case gemmi::CrystalSystem::Tetragonal: number = 89; break; // P422 + case gemmi::CrystalSystem::Trigonal: number = 155; break; // R32 + case gemmi::CrystalSystem::Hexagonal: number = 177; break; // P622 + case gemmi::CrystalSystem::Cubic: number = 207; break; // P432 + } + const auto* sg = gemmi::find_spacegroup_by_number(number); + return sg ? RotationSetOf(*sg) : RotationSet{}; + } + + // A candidate point group: its proper rotations (reference setting) and a representative + // symmorphic space group (used when only the point group is wanted, or for display). + struct PointGroupInfo { + RotationSet rotation_set; + std::vector rotations; // non-identity proper rotations + const gemmi::SpaceGroup* representative = nullptr; }; - SpaceGroupAbsenceScore ScoreSystematicAbsences( - const gemmi::GroupOps& gops, - const std::vector& merged, - const SearchSpaceGroupOptions& opt) { - - SpaceGroupAbsenceScore out; - - if (!opt.test_systematic_absences) - return out; - - int n_bins = opt.absence_resolution_bins; - if (n_bins <= 0) - n_bins = 1; - - double min_inv_d2 = std::numeric_limits::infinity(); - double max_inv_d2 = -std::numeric_limits::infinity(); - - for (const auto& r : merged) { - if (!ReflectionPassesFilters(r, opt)) - continue; - - const double inv_d2 = 1.0 / (r.d * r.d); - min_inv_d2 = std::min(min_inv_d2, inv_d2); - max_inv_d2 = std::max(max_inv_d2, inv_d2); - } - - if (!std::isfinite(min_inv_d2) || !std::isfinite(max_inv_d2)) - return out; - - if (max_inv_d2 < min_inv_d2) - std::swap(max_inv_d2, min_inv_d2); - - std::vector bins(static_cast(n_bins)); - - auto bin_index = [&](double d) { - if (n_bins == 1 || max_inv_d2 <= min_inv_d2) - return 0; - - const double inv_d2 = 1.0 / (d * d); - const double t = (inv_d2 - min_inv_d2) / (max_inv_d2 - min_inv_d2); - int idx = static_cast(t * n_bins); - if (idx < 0) - idx = 0; - if (idx >= n_bins) - idx = n_bins - 1; - return idx; - }; - - for (const auto& r : merged) { - if (!ReflectionPassesFilters(r, opt)) - continue; - - const int idx = bin_index(r.d); - auto& bin = bins[static_cast(idx)]; - bin.d_min_A = std::min(bin.d_min_A, r.d); - bin.d_max_A = std::max(bin.d_max_A, r.d); - - const double i_over_sigma = std::max(0.0, r.I / r.sigma); - const gemmi::Op::Miller hkl{{r.h, r.k, r.l}}; - - if (gops.is_systematically_absent(hkl)) { - bin.absent_sum += i_over_sigma; - bin.absent_count += 1; - out.absent_reflections += 1; - } else { - bin.allowed_sum += i_over_sigma; - bin.allowed_count += 1; - out.allowed_reflections += 1; - } - } - - double global_absent_sum = 0.0; - double global_allowed_sum = 0.0; - double weighted_ratio_sum = 0.0; - int weighted_ratio_weight = 0; - double worst_ratio = 0.0; - bool any_decision_bin_failed = false; - - out.bins.reserve(bins.size()); - - for (const auto& bin : bins) { - SpaceGroupAbsenceBinScore bin_score; - if (std::isfinite(bin.d_min_A)) - bin_score.d_min_A = bin.d_min_A; - if (bin.d_max_A > 0.0) - bin_score.d_max_A = bin.d_max_A; - - bin_score.absent_reflections = bin.absent_count; - bin_score.allowed_reflections = bin.allowed_count; - - if (bin.absent_count > 0) - bin_score.mean_absent_i_over_sigma = bin.absent_sum / static_cast(bin.absent_count); - if (bin.allowed_count > 0) - bin_score.mean_allowed_i_over_sigma = bin.allowed_sum / static_cast(bin.allowed_count); - - global_absent_sum += bin.absent_sum; - global_allowed_sum += bin.allowed_sum; - - if (bin.absent_count >= opt.min_absent_reflections_per_bin && - bin.allowed_count >= opt.min_allowed_reflections_per_bin && - bin_score.mean_allowed_i_over_sigma > 0.0) { - - bin_score.used_for_decision = true; - bin_score.absent_to_allowed_ratio = - bin_score.mean_absent_i_over_sigma / bin_score.mean_allowed_i_over_sigma; - bin_score.accepted = - bin_score.absent_to_allowed_ratio <= opt.max_absent_to_allowed_i_over_sigma_ratio_in_any_bin; - - out.compared_bins += 1; - if (bin_score.accepted) - out.accepted_bins += 1; - else - any_decision_bin_failed = true; - - weighted_ratio_sum += bin_score.absent_to_allowed_ratio * bin.absent_count; - weighted_ratio_weight += bin.absent_count; - worst_ratio = std::max(worst_ratio, bin_score.absent_to_allowed_ratio); - } - - out.bins.push_back(bin_score); - } - - if (out.absent_reflections > 0) - out.mean_absent_i_over_sigma = global_absent_sum / static_cast(out.absent_reflections); - if (out.allowed_reflections > 0) - out.mean_allowed_i_over_sigma = global_allowed_sum / static_cast(out.allowed_reflections); - - if (weighted_ratio_weight > 0) { - out.weighted_absent_to_allowed_ratio = weighted_ratio_sum / static_cast(weighted_ratio_weight); - out.worst_absent_to_allowed_ratio = worst_ratio; - out.accepted = - !any_decision_bin_failed && - out.weighted_absent_to_allowed_ratio <= opt.max_absent_to_allowed_i_over_sigma_ratio; - } else if (out.absent_reflections == 0) { - out.weighted_absent_to_allowed_ratio = 0.0; - out.worst_absent_to_allowed_ratio = 0.0; - out.accepted = true; - } else if (out.mean_allowed_i_over_sigma > 0.0) { - const double global_ratio = out.mean_absent_i_over_sigma / out.mean_allowed_i_over_sigma; - out.weighted_absent_to_allowed_ratio = global_ratio; - out.worst_absent_to_allowed_ratio = global_ratio; - out.accepted = global_ratio <= opt.max_absent_to_allowed_i_over_sigma_ratio_in_any_bin; - } else { - out.weighted_absent_to_allowed_ratio = std::numeric_limits::infinity(); - out.worst_absent_to_allowed_ratio = std::numeric_limits::infinity(); - out.accepted = false; - } - - return out; - } - - bool IsCenteringCompatible(char requested, char candidate) { - if (requested == '\0') - return true; - return std::toupper(static_cast(requested)) == - std::toupper(static_cast(candidate)); - } - - bool IsCandidateSpaceGroup(const gemmi::SpaceGroup& sg, - const std::optional& crystal_system, - char centering) { - if (!sg.is_reference_setting()) - return false; - if (!sg.is_sohncke()) - return false; - if (crystal_system.has_value() && sg.crystal_system() != crystal_system.value()) - return false; - if (!IsCenteringCompatible(centering, sg.centring_type())) - return false; - return true; - } - - std::vector EnumerateCandidateSpaceGroups( - const std::optional& crystal_system, - char centering) { - - std::vector out; + // Enumerate candidate point groups. When a holohedry is given (from the lattice metric), keep + // only its subgroups - this both skips operators the lattice forbids and avoids accepting a + // coincidental higher symmetry; all subgroups down to P1 are still candidates. + std::vector EnumeratePointGroups(const std::optional& holohedry) { + std::vector out; + std::map index; for (const auto& sg : gemmi::spacegroup_tables::main) { - if (!IsCandidateSpaceGroup(sg, crystal_system, centering)) + if (!sg.is_sohncke() || !sg.is_reference_setting()) continue; - out.push_back(sg); - } - std::sort(out.begin(), out.end(), - [](const gemmi::SpaceGroup& a, const gemmi::SpaceGroup& b) { - const int order_a = a.operations().derive_symmorphic().order(); - const int order_b = b.operations().derive_symmorphic().order(); - if (order_a != order_b) - return order_a < order_b; - return a.number < b.number; - }); + RotationSet rs = RotationSetOf(sg); + if (holohedry.has_value() && + !std::includes(holohedry->begin(), holohedry->end(), rs.begin(), rs.end())) + continue; + + auto it = index.find(rs); + size_t pos; + if (it == index.end()) { + PointGroupInfo info; + for (const auto& op : sg.operations().derive_symmorphic().sym_ops) { + if (op.rot == gemmi::Op::identity().rot) + continue; + info.rotations.push_back(gemmi::Op{op.rot, {0, 0, 0}, op.notation}); + } + info.rotation_set = rs; + pos = out.size(); + index[rs] = pos; + out.push_back(std::move(info)); + } else { + pos = it->second; + } + // Prefer a symmorphic representative (the plain point-group setting). + auto& info = out[pos]; + if (info.representative == nullptr || + (!info.representative->is_symmorphic() && sg.is_symmorphic())) + info.representative = &sg; + } return out; } } @@ -394,97 +168,219 @@ SearchSpaceGroupResult SearchSpaceGroup( if (merged.empty()) return result; - const auto lookup = BuildMergedLookup(merged, opt); - const auto candidates = EnumerateCandidateSpaceGroups(opt.crystal_system, opt.centering); + const size_t n = merged.size(); - for (const auto& sg : candidates) { - SpaceGroupCandidateScore score{.space_group = sg}; + // Flatten the reflections and mark which ones each stage may use. The correlation stage drops + // weak reflections; the absence stage must keep them - that is where the screw-axis signal is. + std::vector H(n), K(n), L(n); + std::vector I(n), IoverSigma(n); + std::vector key(n); + std::vector pass_absence(n, 0), pass_cc(n, 0); - const gemmi::GroupOps gops_full = sg.operations(); - const gemmi::GroupOps gops_rot = gops_full.derive_symmorphic(); + for (size_t i = 0; i < n; ++i) { + const auto& r = merged[i]; + H[i] = r.h; K[i] = r.k; L[i] = r.l; + I[i] = r.I; + key[i] = Canonicalize(r.h, r.k, r.l, opt.merge_friedel); - double cc_sum = 0.0; - int cc_count = 0; - int compared_total = 0; - double min_cc = std::numeric_limits::infinity(); + const bool finite = std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0 && + std::isfinite(r.d) && r.d > 0; + const bool in_range = finite && (opt.d_min_limit_A <= 0 || r.d >= opt.d_min_limit_A); + IoverSigma[i] = finite ? r.I / r.sigma : 0.0; + pass_absence[i] = in_range; + // The correlation stage uses only genuinely-present reflections. Near-zero (systematically + // absent) reflections would otherwise form a second cluster at the origin and fake a high + // correlation for false operators - fatal on centered lattices, where half the reflections + // are extinct. + pass_cc[i] = in_range && IoverSigma[i] >= opt.present_i_over_sigma && + (opt.min_i_over_sigma <= 0 || IoverSigma[i] >= opt.min_i_over_sigma); + } - for (const auto& op : gops_rot.sym_ops) { - if (IsIdentityOp(op)) + std::unordered_map key_to_index; + key_to_index.reserve(n * 2); + for (size_t i = 0; i < n; ++i) + if (pass_absence[i]) + key_to_index.emplace(key[i], static_cast(i)); + + // --- Stage A: score each distinct rotation operator once --- + std::vector visited(n, 0); + uint32_t epoch = 0; + + auto score_operator = [&](const gemmi::Op& op) -> SpaceGroupOperatorScore { + ++epoch; + std::vector x, y; + for (size_t i = 0; i < n; ++i) { + if (!pass_cc[i] || visited[i] == epoch) continue; - if (IsParityChanging(op)) + const auto m2 = op.apply_to_hkl(gemmi::Op::Miller{{H[i], K[i], L[i]}}); + const HKLKey k2 = Canonicalize(m2[0], m2[1], m2[2], opt.merge_friedel); + if (k2 == key[i]) + continue; // reflection lies on this rotation axis + const auto it = key_to_index.find(k2); + if (it == key_to_index.end()) continue; - - auto op_score = ScoreOperator(op, merged, lookup, opt); - if (!op_score.has_value()) + const int j = it->second; + if (!pass_cc[j]) continue; + x.push_back(I[i]); + y.push_back(I[j]); + visited[i] = epoch; + visited[j] = epoch; + } + SpaceGroupOperatorScore s; + s.op_triplet_hkl = op.as_hkl().triplet('h'); + s.n_pairs = static_cast(x.size()); + s.cc = PearsonCC(x, y); + s.present = s.n_pairs >= opt.min_pairs_per_operator && std::isfinite(s.cc) && + s.cc >= opt.min_operator_cc; + return s; + }; - compared_total += op_score->compared; - score.operator_scores.push_back(*op_score); + std::map, SpaceGroupOperatorScore> op_cache; + auto operator_score = [&](const gemmi::Op& op) -> const SpaceGroupOperatorScore& { + const auto rk = RotKey(op); + auto it = op_cache.find(rk); + if (it != op_cache.end()) + return it->second; + return op_cache.emplace(rk, score_operator(op)).first->second; + }; - if (op_score->compared >= opt.min_pairs_per_operator && std::isfinite(op_score->cc)) { - cc_sum += op_score->cc; - min_cc = std::min(min_cc, op_score->cc); - cc_count += 1; - if (op_score->accepted) - score.accepted_operators += 1; + std::optional holohedry; + if (opt.lattice_system.has_value()) + holohedry = HolohedryRotationSet(opt.lattice_system.value()); + const auto point_groups = EnumeratePointGroups(holohedry); + + // Choose the largest point group all of whose operators are confirmed (ties -> higher min CC). + const PointGroupInfo* best_pg = nullptr; + int best_pg_order = 0; + double best_pg_min_cc = -2.0; + + for (const auto& pg : point_groups) { + bool all_present = true; + double min_cc = pg.rotations.empty() ? 1.0 : std::numeric_limits::infinity(); + for (const auto& op : pg.rotations) { + const auto& s = operator_score(op); + all_present = all_present && s.present; + min_cc = std::min(min_cc, s.cc); + } + if (!all_present) + continue; + + const int order = static_cast(pg.rotations.size()) + 1; + if (order > best_pg_order || (order == best_pg_order && min_cc > best_pg_min_cc)) { + best_pg = &pg; + best_pg_order = order; + best_pg_min_cc = min_cc; + } + } + + for (const auto& [rk, s] : op_cache) + result.operator_scores.push_back(s); + std::sort(result.operator_scores.begin(), result.operator_scores.end(), + [](const auto& a, const auto& b) { return a.cc > b.cc; }); + + if (best_pg == nullptr) // should not happen (C1 always qualifies) + return result; + + if (best_pg->representative) + result.point_group_hm = best_pg->representative->point_group_hm(); + + // --- Stage B: pick the space group within the point group --- + // Without screw/centering determination, return the symmorphic representative. + if (!opt.determine_space_group || best_pg->rotations.empty()) { + if (best_pg->representative) + result.best_space_group = *best_pg->representative; + return result; + } + + for (const auto& sg : gemmi::spacegroup_tables::main) { + if (!sg.is_sohncke() || !sg.is_reference_setting() || RotationSetOf(sg) != best_pg->rotation_set) + continue; + // Centering is a metric property: when the lattice centering is known, only weigh P against + // that centering, never a different one. + if (opt.lattice_centering != '\0' && + sg.centring_type() != 'P' && sg.centring_type() != opt.lattice_centering) + continue; + + const gemmi::GroupOps gops = sg.operations(); + SpaceGroupCandidateScore s{.space_group = sg}; + double absent_sum = 0, present_sum = 0; + int present_n = 0; + + // Judge centering and screw/glide absences on separate reflection sets. Lumping them lets + // a large, correct centering-absent set hide a few strong screw violations and over-claim + // screw axes (e.g. I4_132 on I432 data). + int centering_absent = 0, centering_violations = 0; + int screw_absent = 0, screw_violations = 0; + + for (size_t i = 0; i < n; ++i) { + if (!pass_absence[i]) + continue; + const gemmi::Op::Miller hkl{{H[i], K[i], L[i]}}; + const bool present = IoverSigma[i] > opt.present_i_over_sigma; + + if (CenteringAbsent(gops, hkl)) { + s.absent_observed += 1; + absent_sum += IoverSigma[i]; + centering_absent += 1; + if (present) { s.absent_violations += 1; centering_violations += 1; } + } else if (gops.is_systematically_absent(hkl)) { + s.absent_observed += 1; + absent_sum += IoverSigma[i]; + screw_absent += 1; + if (present) { s.absent_violations += 1; screw_violations += 1; } + } else { + present_n += 1; + present_sum += IoverSigma[i]; } } - score.absence_score = ScoreSystematicAbsences(gops_full, merged, opt); - score.tested_operators = static_cast(score.operator_scores.size()); - score.compared_total = compared_total; - score.mean_cc = (cc_count > 0) ? (cc_sum / cc_count) : 0.0; - score.min_cc = std::isfinite(min_cc) ? min_cc : 0.0; + if (s.absent_observed > 0) + s.absent_mean_i_over_sigma = absent_sum / s.absent_observed; + if (present_n > 0) + s.present_mean_i_over_sigma = present_sum / present_n; - const bool trivial_group = (gops_rot.order() <= 1); - const bool rotationally_accepted = - trivial_group || - ((score.tested_operators > 0) && - (score.accepted_operators == score.tested_operators) && - (score.compared_total >= opt.min_total_compared)); - - score.accepted = rotationally_accepted && score.absence_score.accepted; - - result.candidates.push_back(std::move(score)); + const bool centering_ok = centering_absent == 0 || + centering_violations <= opt.max_absent_violation_fraction * centering_absent; + const bool screw_ok = screw_absent == 0 || + screw_violations <= opt.max_absent_violation_fraction * screw_absent; + s.consistent = centering_ok && screw_ok; + result.candidates.push_back(std::move(s)); } + // A candidate is eligible when its absences are confirmed and there are enough of them to + // trust (the symmorphic group, with no absences, is always eligible as the fallback). Rank + // eligible candidates by how many absences they explain - the screw/centering content that is + // both real and maximal wins, instead of defaulting to the symmorphic group. + auto eligible = [&](const SpaceGroupCandidateScore& s) { + return s.consistent && (s.absent_observed == 0 || s.absent_observed >= opt.min_absent_observed); + }; std::sort(result.candidates.begin(), result.candidates.end(), - [](const SpaceGroupCandidateScore& a, const SpaceGroupCandidateScore& b) { - if (a.accepted != b.accepted) - return a.accepted > b.accepted; - - if (a.absence_score.weighted_absent_to_allowed_ratio != - b.absence_score.weighted_absent_to_allowed_ratio) - return a.absence_score.weighted_absent_to_allowed_ratio < - b.absence_score.weighted_absent_to_allowed_ratio; - - if (a.absence_score.worst_absent_to_allowed_ratio != - b.absence_score.worst_absent_to_allowed_ratio) - return a.absence_score.worst_absent_to_allowed_ratio < - b.absence_score.worst_absent_to_allowed_ratio; - - const int order_a = a.space_group.operations().derive_symmorphic().order(); - const int order_b = b.space_group.operations().derive_symmorphic().order(); - if (order_a != order_b) - return order_a > order_b; - - if (a.absence_score.absent_reflections != b.absence_score.absent_reflections) - return a.absence_score.absent_reflections > b.absence_score.absent_reflections; - - if (a.accepted_operators != b.accepted_operators) - return a.accepted_operators > b.accepted_operators; - if (a.min_cc != b.min_cc) - return a.min_cc > b.min_cc; - if (a.mean_cc != b.mean_cc) - return a.mean_cc > b.mean_cc; - - return a.space_group.number > b.space_group.number; + [&](const SpaceGroupCandidateScore& a, const SpaceGroupCandidateScore& b) { + if (eligible(a) != eligible(b)) + return eligible(a); + if (a.absent_observed != b.absent_observed) + return a.absent_observed > b.absent_observed; + // Tie: prefer the metric centering, then the lower space-group number. + if (opt.lattice_centering != '\0') { + const bool am = a.space_group.centring_type() == opt.lattice_centering; + const bool bm = b.space_group.centring_type() == opt.lattice_centering; + if (am != bm) + return am; + } + return a.space_group.number < b.space_group.number; }); - for (const auto& cand : result.candidates) { - if (cand.accepted) { - result.best_space_group = cand.space_group; - break; + if (!result.candidates.empty() && eligible(result.candidates.front())) { + const int best_absent = result.candidates.front().absent_observed; + for (auto& s : result.candidates) { + if (!eligible(s) || s.absent_observed != best_absent) + continue; + s.selected = true; + if (!result.best_space_group.has_value()) + result.best_space_group = s.space_group; // representative (lowest number) + else + result.alternatives.push_back(s.space_group); } } @@ -495,89 +391,45 @@ std::string SearchSpaceGroupResultToText(const SearchSpaceGroupResult& result, size_t max_candidates_to_print) { std::ostringstream os; - os << "Space-group candidates\n"; - os << " " - << std::setw(10) << "SG" - << " " - << std::setw(4) << "Acc" - << " " - << std::setw(8) << "" - << " " - << std::setw(8) << "minCC" - << " " - << std::setw(9) << "compared" - << " " - << std::setw(7) << "ops" - << " " - << std::setw(5) << "AbsOK" - << " " - << std::setw(8) << "Nabs" - << " " - << std::setw(8) << "Nallow" - << " " - << std::setw(8) << "Abs/All" - << " " - << std::setw(8) << "worst" - << "\n"; + os << "Point group: " << (result.point_group_hm.empty() ? "?" : result.point_group_hm) + << " (from intensity correlations)\n"; - os << " " - << std::setw(10) << "----------" - << " " - << std::setw(4) << "----" - << " " - << std::setw(8) << "--------" - << " " - << std::setw(8) << "--------" - << " " - << std::setw(9) << "---------" - << " " - << std::setw(7) << "-------" - << " " - << std::setw(5) << "-----" - << " " - << std::setw(8) << "--------" - << " " - << std::setw(8) << "--------" - << " " - << std::setw(8) << "--------" - << " " - << std::setw(8) << "--------" - << "\n"; - - const size_t n = std::min(max_candidates_to_print, result.candidates.size()); - for (size_t i = 0; i < n; ++i) { - const auto& c = result.candidates[i]; - os << " " - << std::setw(10) << c.space_group.short_name() - << " " - << std::setw(4) << (c.accepted ? "yes" : "no") - << " " - << std::setw(8) << std::fixed << std::setprecision(3) << c.mean_cc - << " " - << std::setw(8) << std::fixed << std::setprecision(3) << c.min_cc - << " " - << std::setw(9) << c.compared_total - << " " - << std::setw(3) << c.accepted_operators << "/" << std::setw(3) << c.tested_operators - << " " - << std::setw(5) << (c.absence_score.accepted ? "yes" : "no") - << " " - << std::setw(8) << c.absence_score.absent_reflections - << " " - << std::setw(8) << c.absence_score.allowed_reflections - << " " - << std::setw(8) << std::fixed << std::setprecision(3) - << c.absence_score.weighted_absent_to_allowed_ratio - << " " - << std::setw(8) << std::fixed << std::setprecision(3) - << c.absence_score.worst_absent_to_allowed_ratio - << "\n"; + os << " " << std::setw(14) << std::left << "operator" << std::right + << std::setw(9) << "CC" << std::setw(10) << "pairs" << std::setw(9) << "symm" << "\n"; + for (const auto& s : result.operator_scores) { + os << " " << std::setw(14) << std::left << s.op_triplet_hkl << std::right + << std::setw(9) << std::fixed << std::setprecision(3) << s.cc + << std::setw(10) << s.n_pairs + << std::setw(9) << (s.present ? "yes" : "no") << "\n"; } - if (result.best_space_group.has_value()) - os << "Best space group: " << result.best_space_group->short_name() << "\n"; - else - os << "Best space group: none accepted\n"; + os << "\nSpace-group candidates\n"; + os << " " << std::setw(10) << std::left << "SG" << std::right + << std::setw(9) << "absent" << std::setw(7) << "viol" + << std::setw(11) << "abs" << std::setw(11) << "pres" + << std::setw(6) << "OK" << "\n"; + + const size_t count = std::min(max_candidates_to_print, result.candidates.size()); + for (size_t i = 0; i < count; ++i) { + const auto& c = result.candidates[i]; + os << (c.selected ? "* " : " ") + << std::setw(10) << std::left << c.space_group.short_name() << std::right + << std::setw(9) << c.absent_observed << std::setw(7) << c.absent_violations + << std::setw(11) << std::fixed << std::setprecision(2) << c.absent_mean_i_over_sigma + << std::setw(11) << std::fixed << std::setprecision(2) << c.present_mean_i_over_sigma + << std::setw(6) << (c.consistent ? "yes" : "no") << "\n"; + } + + if (result.best_space_group.has_value()) { + os << "Best space group: " << result.best_space_group->short_name(); + for (const auto& alt : result.alternatives) + os << " or " << alt.short_name(); + if (!result.alternatives.empty()) + os << " (indistinguishable from these data)"; + os << "\n"; + } else { + os << "Best space group: none determined\n"; + } return os.str(); -} \ No newline at end of file +} diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index 00c56ba9..e65c7e24 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -7,95 +7,97 @@ #include #include -#include "Merge.h" +#include "../../common/Reflection.h" #include "gemmi/symmetry.hpp" +// Determine the likely space group of a dataset from its P1-merged intensities, in the spirit +// of POINTLESS (Evans 2006): +// +// Stage A - point group (Laue) symmetry. Every candidate rotation operator is scored once by +// the correlation of I(h) with I(Rh). The chosen point group is the largest one all +// of whose operators are confirmed (high CC). A wrong operator scores ~0, so this is +// self-pruning - the unit cell is not needed. +// +// Stage B - space group within that point group. Each Sohncke space group of the point group +// predicts a set of systematically absent reflections (centering + screw axes). The +// chosen space group is the one that explains the MOST absences while every reflection +// it predicts absent is in fact weak. The symmorphic group (no absences) is the +// fallback when no screw/centering is supported. Enantiomorphic pairs (e.g. P4_1 vs +// P4_3) are indistinguishable from intensities and are reported as a pair. + struct SpaceGroupOperatorScore { - std::string op_triplet_hkl; - double cc = 0.0; - int compared = 0; - bool accepted = false; -}; - -struct SpaceGroupAbsenceBinScore { - double d_min_A = 0.0; - double d_max_A = 0.0; - int absent_reflections = 0; - int allowed_reflections = 0; - double mean_absent_i_over_sigma = 0.0; - double mean_allowed_i_over_sigma = 0.0; - double absent_to_allowed_ratio = 0.0; - bool used_for_decision = false; - bool accepted = true; -}; - -struct SpaceGroupAbsenceScore { - int absent_reflections = 0; - int allowed_reflections = 0; - int compared_bins = 0; - int accepted_bins = 0; - double mean_absent_i_over_sigma = 0.0; - double mean_allowed_i_over_sigma = 0.0; - double weighted_absent_to_allowed_ratio = 0.0; - double worst_absent_to_allowed_ratio = 0.0; - bool accepted = true; - std::vector bins; + std::string op_triplet_hkl; // reciprocal-space triplet of the rotation, e.g. "-h,-k,l" + double cc = 0.0; // correlation of I(h) with I(Rh) + int n_pairs = 0; // independent reflection pairs the CC was computed from + bool present = false; // operator confirmed as a real symmetry of the intensities }; struct SpaceGroupCandidateScore { gemmi::SpaceGroup space_group; - double mean_cc = 0.0; - double min_cc = 0.0; - int compared_total = 0; - int accepted_operators = 0; - int tested_operators = 0; - bool accepted = false; - SpaceGroupAbsenceScore absence_score; - std::vector operator_scores; + int absent_observed = 0; // observed reflections this SG predicts systematically absent + int absent_violations = 0; // of those, how many are nonetheless strongly present + double absent_mean_i_over_sigma = 0.0; + double present_mean_i_over_sigma = 0.0; + bool consistent = false; // absent class confirmed weak (few violations) + bool selected = false; // chosen result (or its enantiomorph) }; struct SearchSpaceGroupOptions { - // If not set, search all crystal systems. - std::optional crystal_system; + // Lattice (metric) symmetry from LatticeSearch. When set, the point-group search is limited to + // the Sohncke subgroups of this system's holohedry - the metric is an upper bound on the + // intensity symmetry, so e.g. a tetragonal metric tests {1,2,222,4,422} and never 3/6/23. All + // subgroups are still tested (down to P1), so a pseudo-symmetric metric never forces a higher + // symmetry than the intensities support. Unset = search every Sohncke system. + std::optional lattice_system; - // '\0' means: search all centerings compatible with the candidate system. - char centering = '\0'; + // Lattice centering from LatticeSearch ('P','C','I','F','R'; '\0' = unknown). Centering is a + // metric property of the lattice, so when known the space-group search only weighs P vs this + // centering and prefers the metric one when the data cannot decide between them. + char lattice_centering = '\0'; - // If true, Friedel mates are treated as equivalent when matching HKLs. + // Friedel mates are treated as equivalent when matching HKLs (i.e. anomalous signal ignored). bool merge_friedel = true; - // Ignore reflections beyond this resolution (higher-resolution means smaller d). - // Set to 0 to disable the cutoff. + // Ignore reflections beyond this resolution (smaller d = higher resolution). 0 disables. double d_min_limit_A = 0.0; - // Optional I/sigma filter on merged reflections. + // Drop reflections weaker than this from the correlation stage only (the absence stage must + // keep weak reflections - that is where the screw-axis signal lives). double min_i_over_sigma = 0.0; - // Acceptance thresholds for each tested rotational operator. - double min_operator_cc = 0.80; + // --- 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. + double min_operator_cc = 0.5; int min_pairs_per_operator = 20; - // Minimum total number of compared reflection pairs for a candidate SG. - int min_total_compared = 100; + // --- Stage B: space group (screw axes / centering) --- + bool determine_space_group = true; // false: stop at the symmorphic representative - // Test systematic absences from centering / screws / glides. - bool test_systematic_absences = true; + // Signed I/sigma above which a reflection counts as genuinely present. Signed, so negative + // noise is never mistaken for a real reflection. Used both to spot reflections that violate a + // wrongly-assumed absence, and to keep the correlation stage from pairing near-zero reflections. + double present_i_over_sigma = 3.0; - // Number of resolution bins used for absence-vs-allowed comparison. - int absence_resolution_bins = 10; + // A candidate's absence conditions are accepted when at most this fraction of the reflections + // it predicts absent are in fact strongly present. + double max_absent_violation_fraction = 0.10; - // Minimum counts in a bin before it contributes to the decision. - int min_absent_reflections_per_bin = 5; - int min_allowed_reflections_per_bin = 10; - - // Acceptance thresholds for absent/allowed 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; + // Need at least this many observed reflections in the predicted-absent class before a + // screw/centering is claimed (guards against deciding from a handful of reflections). + int min_absent_observed = 8; }; struct SearchSpaceGroupResult { std::optional best_space_group; - std::vector candidates; + // Other space groups that fit the data equally well (same systematic absences): enantiomorphic + // partners (P4_1 vs P4_3), origin-ambiguous pairs (I222 vs I2_12_12_1), or groups left + // undetermined by incomplete data. The data cannot choose between best_space_group and these. + std::vector alternatives; + + std::string point_group_hm; // chosen point group, e.g. "422" + std::vector operator_scores; // Stage A, all distinct operators tested + std::vector candidates; // Stage B, ranked }; SearchSpaceGroupResult SearchSpaceGroup( @@ -104,4 +106,4 @@ SearchSpaceGroupResult SearchSpaceGroup( std::string SearchSpaceGroupResultToText( const SearchSpaceGroupResult& result, - size_t max_candidates_to_print = 20); \ No newline at end of file + size_t max_candidates_to_print = 20); diff --git a/image_analysis/scale_merge/TwinningAnalysis.cpp b/image_analysis/scale_merge/TwinningAnalysis.cpp new file mode 100644 index 00000000..97b70c30 --- /dev/null +++ b/image_analysis/scale_merge/TwinningAnalysis.cpp @@ -0,0 +1,156 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "TwinningAnalysis.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + int64_t PackHKL(int h, int k, int l) { + constexpr int64_t bias = 1 << 20; // indices assumed within +/- 2^20 + return ((h + bias) << 42) | ((k + bias) << 21) | (l + bias); + } + + bool UsableIntensity(const MergedReflection& r) { + return std::isfinite(r.I) && std::isfinite(r.d) && r.d > 0.0; + } +} + +TwinningAnalysisResult AnalyzeTwinning(const std::vector& merged, + const gemmi::SpaceGroup* space_group, int resolution_shells) { + TwinningAnalysisResult result; + if (merged.empty()) + return result; + + // Centric reflections follow different statistics and must be excluded. In P1 (no space group) + // none are centric. + const gemmi::GroupOps gops = space_group ? space_group->operations() : gemmi::GroupOps{}; + auto acentric = [&](const MergedReflection& r) { + return !space_group || !gops.is_reflection_centric(gemmi::Op::Miller{{r.h, r.k, r.l}}); + }; + + // --- L-test --- + // Pair each reflection with a symmetry-independent neighbour two steps away along an axis (the + // step of 2 keeps the partner local in resolution while avoiding the reflection itself). The + // merged reflections are unique in the asymmetric unit, so any other merged reflection is + // genuinely non-equivalent - exactly the pairing the L-test wants. Only acentric reflections + // with positive intensity enter, which also keeps L = (I1-I2)/(I1+I2) bounded in [-1, 1]. + std::unordered_map intensity; + intensity.reserve(merged.size() * 2); + for (const auto& r : merged) + if (UsableIntensity(r) && r.I > 0.0 && acentric(r)) + intensity.emplace(PackHKL(r.h, r.k, r.l), r.I); + + const std::array, 3> steps{{{2, 0, 0}, {0, 2, 0}, {0, 0, 2}}}; + + double sum_abs_l = 0.0, sum_l2 = 0.0; + int n_pairs = 0; + for (const auto& [key, i1] : intensity) { + const int h = static_cast((key >> 42) & 0x1FFFFF) - (1 << 20); + const int k = static_cast((key >> 21) & 0x1FFFFF) - (1 << 20); + const int l = static_cast(key & 0x1FFFFF) - (1 << 20); + for (const auto& s : steps) { + const auto it = intensity.find(PackHKL(h + s[0], k + s[1], l + s[2])); + if (it == intensity.end()) + continue; + const double lstat = (i1 - it->second) / (i1 + it->second); + sum_abs_l += std::fabs(lstat); + sum_l2 += lstat * lstat; + ++n_pairs; + break; // one neighbour per reflection + } + } + result.l_test_pairs = n_pairs; + if (n_pairs > 0) { + result.mean_abs_l = sum_abs_l / n_pairs; + result.mean_l_squared = sum_l2 / n_pairs; + } + + // --- Second moment /^2 of acentric intensities, normalised per resolution shell --- + // Binning by 1/d^2 removes the resolution fall-off, so the moment is 2.0 (untwinned) or 1.5 + // (perfect twin) regardless of the overall B-factor. + int n_shells = std::max(1, resolution_shells); + double min_s = std::numeric_limits::infinity(); + double max_s = -std::numeric_limits::infinity(); + for (const auto& r : merged) { + if (!UsableIntensity(r) || !acentric(r)) + continue; + const double s = 1.0 / (r.d * r.d); + min_s = std::min(min_s, s); + max_s = std::max(max_s, s); + } + if (std::isfinite(min_s) && max_s > min_s) { + auto shell_of = [&](double d) { + const double t = (1.0 / (d * d) - min_s) / (max_s - min_s); + return std::min(n_shells - 1, std::max(0, static_cast(t * n_shells))); + }; + std::vector shell_sum(n_shells, 0.0); + std::vector shell_n(n_shells, 0); + for (const auto& r : merged) { + if (!UsableIntensity(r) || !acentric(r)) + continue; + const int b = shell_of(r.d); + shell_sum[b] += r.I; + shell_n[b] += 1; + } + std::vector shell_mean(n_shells, 0.0); + for (int b = 0; b < n_shells; ++b) + if (shell_n[b] > 0) + shell_mean[b] = shell_sum[b] / shell_n[b]; + + double sum_e4 = 0.0; + int n_moment = 0; + for (const auto& r : merged) { + if (!UsableIntensity(r) || !acentric(r)) + continue; + const double mean = shell_mean[shell_of(r.d)]; + if (mean <= 0.0) + continue; + const double e2 = r.I / mean; + sum_e4 += e2 * e2; + ++n_moment; + } + result.moment_reflections = n_moment; + if (n_moment > 0) + result.second_moment = sum_e4 / n_moment; + } + + // Twin fraction from the second moment M = 2(1 - a + a^2): a = (1 - sqrt(2M-3))/2. + if (result.second_moment > 0.0) { + const double m = std::clamp(result.second_moment, 1.5, 2.0); + result.estimated_twin_fraction = (1.0 - std::sqrt(std::max(0.0, 2.0 * m - 3.0))) / 2.0; + } + // Either indicator dropping clearly below its untwinned value is suspicious. + result.twinning_suspected = (result.l_test_pairs > 0 && result.mean_abs_l < 0.44) || + (result.moment_reflections > 0 && result.second_moment < 1.85); + return result; +} + +std::string TwinningAnalysisToText(const TwinningAnalysisResult& result) { + std::ostringstream os; + os << std::fixed << std::setprecision(3); + os << "Twinning analysis\n"; + if (result.l_test_pairs > 0) + os << " L-test (Padilla-Yeates): <|L|> = " << result.mean_abs_l + << ", = " << result.mean_l_squared + << " [untwinned 0.500 / 0.333, perfect twin 0.375 / 0.200; " << result.l_test_pairs + << " pairs]\n"; + if (result.moment_reflections > 0) + os << " Second moment /^2 = " << result.second_moment + << " [untwinned 2.00, perfect twin 1.50]\n"; + if (result.twinning_suspected) + os << " => Twinning suspected (estimated twin fraction ~" + << result.estimated_twin_fraction << "). Statistics flag the presence of twinning, not\n" + << " the twin law; confirm with a dedicated twin-law analysis.\n"; + else + os << " => No twinning indicated.\n"; + return os.str(); +} diff --git a/image_analysis/scale_merge/TwinningAnalysis.h b/image_analysis/scale_merge/TwinningAnalysis.h new file mode 100644 index 00000000..6c0ec7d8 --- /dev/null +++ b/image_analysis/scale_merge/TwinningAnalysis.h @@ -0,0 +1,44 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "../../common/Reflection.h" +#include "gemmi/symmetry.hpp" + +// Intensity-statistics test for merohedral twinning, in the spirit of phenix.xtriage. Two robust, +// normalisation-free indicators on the merged intensities: +// +// - Padilla & Yeates L-test (Acta Cryst D59, 2003): for pairs of nearby but symmetry-independent +// reflections, L = (I1 - I2)/(I1 + I2). Untwinned acentric data give <|L|> = 0.5, = 1/3; +// a perfect merohedral twin narrows this to <|L|> = 0.375, = 0.2. Insensitive to +// anisotropy and overall scale because it is a local ratio. +// +// - Second moment /^2 of acentric intensities, normalised per resolution shell so it is +// 2.0 for untwinned data and 1.5 for a perfect twin (= 2(1 - a + a^2) for twin fraction a). +// +// These flag the *presence* of twinning; they do not by themselves identify the twin law. + +struct TwinningAnalysisResult { + int l_test_pairs = 0; + double mean_abs_l = 0.0; // <|L|>: 0.5 untwinned, 0.375 perfect twin + double mean_l_squared = 0.0; // : 0.333 untwinned, 0.2 perfect twin + + int moment_reflections = 0; + double second_moment = 0.0; // /^2 (shell-normalised): 2.0 untwinned, 1.5 perfect twin + + // Twin fraction estimated from the second moment, a = (1 - sqrt(2M-3))/2, clamped to [0, 0.5]. + double estimated_twin_fraction = 0.0; + bool twinning_suspected = false; +}; + +// The space group (when known) is used to drop centric reflections, which follow different +// statistics (<|L|>, /^2 = 3 not 2); pass nullptr for P1 data, where none are centric. +TwinningAnalysisResult AnalyzeTwinning(const std::vector& merged, + const gemmi::SpaceGroup* space_group = nullptr, + int resolution_shells = 20); + +std::string TwinningAnalysisToText(const TwinningAnalysisResult& result); diff --git a/process/JFJochProcess.cpp b/process/JFJochProcess.cpp index ab7e9d47..e6926784 100644 --- a/process/JFJochProcess.cpp +++ b/process/JFJochProcess.cpp @@ -29,6 +29,7 @@ #include "../image_analysis/scale_merge/Merge.h" #include "../image_analysis/scale_merge/ScaleOnTheFly.h" #include "../image_analysis/scale_merge/SearchSpaceGroup.h" +#include "../image_analysis/scale_merge/TwinningAnalysis.h" #include "../image_analysis/scale_merge/Combine3D.h" #include "../image_analysis/WriteReflections.h" @@ -409,85 +410,111 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { }; phase("Scaling and merging"); - // ScaleOnTheFly self-scaling is only for the no-reference path; with a reference each image - // is already scaled against it during the per-image pass, so we merge directly. - if (config_.reference_data.empty()) { - logger.Info("Running scaling ..."); - ScalingResult scale_result(0); - double t_merge_all = 0.0, t_scale = 0.0; - for (int i = 0; i < config_.scaling_iter; i++) { - phase("Scaling images (iteration " + std::to_string(i + 1) + "/" - + std::to_string(config_.scaling_iter) + ")"); - const auto a = std::chrono::steady_clock::now(); - auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome(), false); - const auto b = std::chrono::steady_clock::now(); - scale_result = indexer->ScaleAllImages(merge_result); - const auto c = std::chrono::steady_clock::now(); - t_merge_all += std::chrono::duration(b - a).count(); - t_scale += std::chrono::duration(c - b).count(); + // Scale the images and merge. Factored so it can run twice: first in P1 to give the + // space-group search a merged dataset, then again in the determined space group so the + // scaling sees symmetry equivalents and the final statistics are in the right symmetry. + // (With a user-fixed space group it simply runs once, already in that symmetry.) + struct ScaleMergeResult { + std::vector merged; + MergeStatistics statistics; + }; + auto scale_and_merge = [&](const std::string &label) -> ScaleMergeResult { + // ScaleOnTheFly self-scaling is only for the no-reference path; with a reference each + // image is already scaled against it during the per-image pass, so we merge directly. + if (config_.reference_data.empty()) { + logger.Info("Running scaling ({}) ...", label); + 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); + } } - logger.Info("[timing] scaling loop ({} iter): MergeAll(serial) {:.2f} s, ScaleAllImages(parallel) {:.2f} s", - config_.scaling_iter, t_merge_all, t_scale); - } - // -P rot3d: weight-sum each reflection's per-frame partials into one full before merging, so - // the error model sees counting statistics (high ISa) instead of rocking-curve slicing scatter. - const bool rot3d = experiment_.GetScalingSettings().GetCombine3D(); - std::vector combined; - if (rot3d) { - phase("Combining 3D partials"); - combined = CombineRotationObservations(indexer->GetIntegrationOutcome(), experiment_, &logger, - config_.observation_dump_path); - } - if (rot3d && experiment_.GetScalingSettings().GetScaleFulls()) { - phase("Scaling fulls (XDS order)"); - ScaleFulls(experiment_, combined, static_cast(config_.scaling_iter), config_.nthreads, logger); - } - const std::vector &merge_input = - rot3d ? combined : indexer->GetIntegrationOutcome(); + // -P rot3d: weight-sum each reflection's per-frame partials into one full before + // merging, so the error model sees counting statistics (high ISa) instead of + // rocking-curve slicing scatter. + const bool rot3d = experiment_.GetScalingSettings().GetCombine3D(); + std::vector combined; + if (rot3d) { + phase("Combining 3D partials"); + combined = CombineRotationObservations(indexer->GetIntegrationOutcome(), experiment_, &logger, + config_.observation_dump_path); + } + if (rot3d && experiment_.GetScalingSettings().GetScaleFulls()) { + phase("Scaling fulls (XDS order)"); + ScaleFulls(experiment_, combined, static_cast(config_.scaling_iter), config_.nthreads, logger); + } + const std::vector &merge_input = + rot3d ? combined : indexer->GetIntegrationOutcome(); - phase("Merging"); - MergeOnTheFly merge_engine(experiment_); - if (result.consensus_cell.has_value()) - merge_engine.ReferenceCell(*result.consensus_cell); - merge_engine.RefineErrorModel(merge_input); - if (merge_engine.ErrorModelActive()) - logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f}", merge_engine.ErrorModelA(), - merge_engine.ErrorModelB(), - merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0); - for (const auto &outcome: merge_input) - merge_engine.AddImage(outcome); + phase("Merging"); + MergeOnTheFly merge_engine(experiment_); + if (result.consensus_cell.has_value()) + merge_engine.ReferenceCell(*result.consensus_cell); + merge_engine.RefineErrorModel(merge_input); + if (merge_engine.ErrorModelActive()) + logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f}", merge_engine.ErrorModelA(), + merge_engine.ErrorModelB(), + merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0); + for (const auto &outcome: merge_input) + merge_engine.AddImage(outcome); - auto merged_reflections = merge_engine.ExportReflections(); - phase("Computing statistics"); - auto merged_statistics = merge_engine.MergeStats(merged_reflections, merge_input, - config_.reference_data); - logger.Info("Merge complete ({} unique reflections)", merged_reflections.size()); + ScaleMergeResult out; + out.merged = merge_engine.ExportReflections(); + phase("Computing statistics"); + out.statistics = merge_engine.MergeStats(out.merged, merge_input, config_.reference_data); + result.error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0; + logger.Info("Merge complete ({} unique reflections, {})", out.merged.size(), label); + return out; + }; + + // First pass: P1 when searching, or directly the user-fixed space group. + const auto initial_sg = experiment_.GetGemmiSpaceGroup(); + auto sm = scale_and_merge(initial_sg ? initial_sg->short_name() : "P1"); std::ostringstream stats_text; if (!experiment_.GetGemmiSpaceGroup().has_value()) { SearchSpaceGroupOptions sg_opts; - sg_opts.crystal_system.reset(); - sg_opts.centering = '\0'; sg_opts.merge_friedel = experiment_.GetScalingSettings().GetMergeFriedel(); sg_opts.d_min_limit_A = experiment_.GetScalingSettings().GetHighResolutionLimit_A().value_or(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; - stats_text << SearchSpaceGroupResultToText(SearchSpaceGroup(merged_reflections, sg_opts)) << "\n\n"; - } - stats_text << merged_statistics; - result.merge_statistics_text = stats_text.str(); + // Constrain the search to subgroups of the lattice (metric) symmetry found by rotation + // indexing, and weigh centering only against the metric one. + if (end_msg.rotation_lattice_type.has_value()) { + sg_opts.lattice_system = end_msg.rotation_lattice_type->crystal_system; + sg_opts.lattice_centering = end_msg.rotation_lattice_type->centering; + } + const auto sg_search = SearchSpaceGroup(sm.merged, sg_opts); + stats_text << SearchSpaceGroupResultToText(sg_search) << "\n\n"; + // Adopt the determined space group and re-scale/merge in it, so scaling uses symmetry + // equivalents and the statistics come out in the right symmetry. P1 stands when nothing + // is determined. + if (sg_search.best_space_group.has_value()) { + const auto &sg = *sg_search.best_space_group; + logger.Info("Adopting space group {} (number {})", sg.short_name(), sg.number); + experiment_.SpaceGroupNumber(sg.number); + end_msg.space_group_number = sg.number; + result.space_group_number = sg.number; + phase("Re-scaling in space group " + sg.short_name()); + sm = scale_and_merge(sg.short_name()); + } + } + + const auto twin_sg_number = experiment_.GetSpaceGroupNumber(); + const gemmi::SpaceGroup *twin_sg = twin_sg_number + ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr; + result.twinning = AnalyzeTwinning(sm.merged, twin_sg); + stats_text << TwinningAnalysisToText(result.twinning) << "\n"; + stats_text << sm.statistics; + result.merge_statistics_text = stats_text.str(); result.has_merge_statistics = true; - result.merge_statistics = merged_statistics; - result.error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0; + result.merge_statistics = sm.statistics; result.has_reference = !config_.reference_data.empty(); if (result.consensus_cell && write_files) { phase("Writing reflections"); - WriteReflections(merged_reflections, *result.consensus_cell, experiment_, config_.output_prefix); + WriteReflections(sm.merged, *result.consensus_cell, experiment_, config_.output_prefix); } phase(""); // flush the last phase's timing } diff --git a/process/JFJochProcess.h b/process/JFJochProcess.h index 8cb89bee..71b2e153 100644 --- a/process/JFJochProcess.h +++ b/process/JFJochProcess.h @@ -16,6 +16,7 @@ #include "../common/JFJochReceiverPlots.h" // MeanProcessingTime #include "../image_analysis/spot_finding/SpotFindingSettings.h" #include "../image_analysis/scale_merge/Merge.h" // MergeStatistics +#include "../image_analysis/scale_merge/TwinningAnalysis.h" // TwinningAnalysisResult class JFJochHDF5Reader; @@ -79,6 +80,13 @@ struct ProcessResult { MergeStatistics merge_statistics; double error_model_isa = 0.0; bool has_reference = false; + + // Twinning analysis of the final merged intensities (l_test_pairs == 0 when not computed). + TwinningAnalysisResult twinning; + + // Space group determined by the search (when the user did not fix one), used for the final + // re-scale/merge and written to the master file. + std::optional space_group_number; }; // Callbacks for progress and live results. Methods may be called from worker threads, so an diff --git a/tests/SearchSpaceGroupTest.cpp b/tests/SearchSpaceGroupTest.cpp index 6e19baf5..c330f5a5 100644 --- a/tests/SearchSpaceGroupTest.cpp +++ b/tests/SearchSpaceGroupTest.cpp @@ -3,6 +3,7 @@ #include "../image_analysis/scale_merge/SearchSpaceGroup.h" #include "gemmi/symmetry.hpp" +#include #include #include #include @@ -126,25 +127,23 @@ TEST_CASE("SearchSpaceGroup detects synthetic space groups") { const auto merged = GenerateMergedReflectionsForSpaceGroup(sg); SearchSpaceGroupOptions opt; - opt.crystal_system.reset(); - opt.centering = '\0'; opt.merge_friedel = true; - opt.d_min_limit_A = 2.0; - opt.min_operator_cc = 0.99; - opt.min_pairs_per_operator = 20; - opt.min_total_compared = 80; - opt.test_systematic_absences = true; - opt.absence_resolution_bins = 10; - opt.min_absent_reflections_per_bin = 5; - opt.min_allowed_reflections_per_bin = 10; - opt.max_absent_to_allowed_i_over_sigma_ratio = 0.05; - opt.max_absent_to_allowed_i_over_sigma_ratio_in_any_bin = 0.10; const auto result = SearchSpaceGroup(merged, opt); + // Several inputs cannot be told apart from intensities alone: enantiomorphic partners + // (P4_3 vs P4_1) and origin-ambiguous pairs (I2_12_12_1 vs I222, I2_13 vs I2_3) share + // the same systematic absences. The search reports those as alternatives, so the + // expected group must appear among the best group and its alternatives. + std::vector accepted; + if (result.best_space_group.has_value()) + accepted.push_back(result.best_space_group->short_name()); + for (const auto& alt : result.alternatives) + accepted.push_back(alt.short_name()); + INFO(SearchSpaceGroupResultToText(result)); REQUIRE(result.best_space_group.has_value()); - CHECK(result.best_space_group->short_name() == tc.expected_short_name); + CHECK(std::find(accepted.begin(), accepted.end(), tc.expected_short_name) != accepted.end()); } } } \ No newline at end of file diff --git a/tools/jfjoch_scale.cpp b/tools/jfjoch_scale.cpp index e7de4ed6..4ce1e338 100644 --- a/tools/jfjoch_scale.cpp +++ b/tools/jfjoch_scale.cpp @@ -21,7 +21,7 @@ #include "../compression/JFJochCompressor.h" #include "../image_analysis/LoadFCalcFromMtz.h" #include "../image_analysis/scale_merge/Merge.h" -#include "../image_analysis/scale_merge/SearchSpaceGroup.h" +#include "../image_analysis/scale_merge/TwinningAnalysis.h" #include "../image_analysis/scale_merge/Combine3D.h" #include "../image_analysis/WriteReflections.h" #include "../image_analysis/UpdateReflectionResolution.h" @@ -310,7 +310,10 @@ int main(int argc, char **argv) { experiment.OverwriteExistingFiles(true); experiment.PolarizationFactor(0.99); experiment.SetFileWriterFormat(FileWriterFormat::NXmxLegacy); - experiment.SpaceGroupNumber(space_group_number); + // Keep the space group read from the input file (written by the full pipeline) unless -S + // overrides it. + if (space_group_number.has_value()) + experiment.SpaceGroupNumber(space_group_number); experiment.NumTriggers(1); ScalingSettings scaling_settings; @@ -388,25 +391,19 @@ int main(int argc, char **argv) { // Print resolution-shell statistics table std::cout << merged_statistics; + // Space-group determination lives in the full jfjoch_process pipeline, where the lattice search + // and a consistent integration are available; this re-scaling tool only consumes a space group + // (from the file's /entry/sample/space_group_number or -S) and merges in it. const bool fixed_space_group = space_group || experiment.GetGemmiSpaceGroup().has_value(); + if (!fixed_space_group) + logger.Warning("No space group in the input file or on the command line - merged in P1. " + "Re-run jfjoch_process (which determines and stores the space group) or pass " + "-S to scale and merge in the correct symmetry."); - if (!fixed_space_group) { - logger.Info("Searching for space group from P1-merged reflections ..."); - - SearchSpaceGroupOptions sg_opts; - sg_opts.crystal_system.reset(); - sg_opts.centering = '\0'; - sg_opts.merge_friedel = experiment.GetScalingSettings().GetMergeFriedel(); - sg_opts.d_min_limit_A = experiment.GetScalingSettings().GetHighResolutionLimit_A().value_or(0.0); - 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; - - const auto sg_search = SearchSpaceGroup(merged_reflections, sg_opts); - std::cout << std::endl << SearchSpaceGroupResultToText(sg_search) << std::endl; - } + const auto twin_sg_number = experiment.GetSpaceGroupNumber(); + const gemmi::SpaceGroup *twin_sg = twin_sg_number + ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr; + std::cout << std::endl << TwinningAnalysisToText(AnalyzeTwinning(merged_reflections, twin_sg)) << std::endl; if (!output_prefix.empty()) WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, output_prefix); diff --git a/viewer/image_viewer/JFJochDiffractionImage.cpp b/viewer/image_viewer/JFJochDiffractionImage.cpp index 9864a268..d00ca443 100644 --- a/viewer/image_viewer/JFJochDiffractionImage.cpp +++ b/viewer/image_viewer/JFJochDiffractionImage.cpp @@ -7,6 +7,7 @@ #include "../../common/DiffractionGeometry.h" #include "../../common/JFJochMath.h" #include "../../common/ROIAzimuthal.h" +#include "../../image_analysis/bragg_integration/SystematicAbsence.h" #include #include @@ -187,7 +188,14 @@ void JFJochDiffractionImage::DrawPredictions() { // Compute current visible area in scene coordinates const QRectF visibleRect = mapToScene(viewport()->geometry()).boundingRect(); + // In space-group search mode the centering-absent reflections are integrated (to confirm the + // centering), but they are not real predictions - keep them out of the overlay. + const char centering = image->ImageData().lattice_type.has_value() + ? image->ImageData().lattice_type->centering : 'P'; + for (const auto &s: image->ImageData().reflections) { + if (systematic_absence(s.h, s.k, s.l, centering)) + continue; // Skip reflections outside the viewport if (!visibleRect.contains(QPointF{s.predicted_x, s.predicted_y})) continue; diff --git a/viewer/windows/JFJochMergeStatsWindow.cpp b/viewer/windows/JFJochMergeStatsWindow.cpp index bf4910b5..78091240 100644 --- a/viewer/windows/JFJochMergeStatsWindow.cpp +++ b/viewer/windows/JFJochMergeStatsWindow.cpp @@ -61,7 +61,8 @@ namespace { } JFJochMergeStatsWindow::JFJochMergeStatsWindow(const QString &title, const MergeStatistics &stats, - double isa, bool has_reference, QWidget *parent) + double isa, bool has_reference, + const TwinningAnalysisResult &twinning, QWidget *parent) : QWidget(parent, Qt::Window), stats_(stats), has_reference_(has_reference) { setWindowTitle("Merge statistics — " + title); setAttribute(Qt::WA_DeleteOnClose); @@ -88,6 +89,27 @@ JFJochMergeStatsWindow::JFJochMergeStatsWindow(const QString &title, const Merge hero->addWidget(MakeCard(pct(o.cc_ref * 100.0), "CCref", this)); layout->addLayout(hero); + // --- Twinning test (Padilla-Yeates L-test + second moment) --- + // A red-outlined banner flags suspected twinning; otherwise a quiet grey "all clear". + if (twinning.l_test_pairs > 0) { + QString text = QString("Twinning: ⟨|L|⟩ = %1 ⟨I²⟩/⟨I⟩² = %2") + .arg(twinning.mean_abs_l, 0, 'f', 3) + .arg(twinning.second_moment, 0, 'f', 2); + if (twinning.twinning_suspected) + text += QString(" — twinning suspected (estimated fraction ≈ %1)") + .arg(twinning.estimated_twin_fraction, 0, 'f', 2); + else + text += " — no twinning indicated"; + + auto *banner = new QLabel(text, this); + banner->setContentsMargins(10, 6, 10, 6); + if (twinning.twinning_suspected) + banner->setStyleSheet("color: rgb(160,0,0); border: 1px solid red; border-radius: 4px;"); + else + banner->setStyleSheet("color: gray;"); + layout->addWidget(banner); + } + // --- Controls: plot/table view toggle (icon buttons) + per-resolution metric selector --- auto *row = new QHBoxLayout(); auto *plotBtn = new QToolButton(this); diff --git a/viewer/windows/JFJochMergeStatsWindow.h b/viewer/windows/JFJochMergeStatsWindow.h index eae1e115..a95eb4fe 100644 --- a/viewer/windows/JFJochMergeStatsWindow.h +++ b/viewer/windows/JFJochMergeStatsWindow.h @@ -6,6 +6,7 @@ #include #include "../../image_analysis/scale_merge/Merge.h" // MergeStatistics +#include "../../image_analysis/scale_merge/TwinningAnalysis.h" // TwinningAnalysisResult class QComboBox; class QLabel; @@ -21,7 +22,8 @@ class JFJochMergeStatsWindow : public QWidget { Q_OBJECT public: JFJochMergeStatsWindow(const QString &title, const MergeStatistics &stats, - double isa, bool has_reference, QWidget *parent = nullptr); + double isa, bool has_reference, + const TwinningAnalysisResult &twinning, QWidget *parent = nullptr); private: MergeStatistics stats_; diff --git a/viewer/windows/JFJochProcessingJobsWindow.cpp b/viewer/windows/JFJochProcessingJobsWindow.cpp index 3fefd5cd..648919c5 100644 --- a/viewer/windows/JFJochProcessingJobsWindow.cpp +++ b/viewer/windows/JFJochProcessingJobsWindow.cpp @@ -343,7 +343,8 @@ void JFJochProcessingJobsWindow::addRowActions(int row, const QString &id) { void JFJochProcessingJobsWindow::showStats(const QString &id) { for (const auto &j: jobs_) { if (j.id == id && j.has_merge_stats) { - auto *win = new JFJochMergeStatsWindow(j.label, j.merge_stats, j.isa, j.merge_has_reference, window()); + auto *win = new JFJochMergeStatsWindow(j.label, j.merge_stats, j.isa, j.merge_has_reference, + j.twinning, window()); win->show(); return; } @@ -496,6 +497,7 @@ void JFJochProcessingJobsWindow::onFinished(ProcessResult result) { jobs_[row].merge_stats = result.merge_statistics; jobs_[row].isa = result.error_model_isa; jobs_[row].merge_has_reference = result.has_reference; + jobs_[row].twinning = result.twinning; if (jobs_[row].graph_btn) jobs_[row].graph_btn->setEnabled(true); showStats(jobs_[row].id); diff --git a/viewer/windows/JFJochProcessingJobsWindow.h b/viewer/windows/JFJochProcessingJobsWindow.h index 346cf6a7..1994b916 100644 --- a/viewer/windows/JFJochProcessingJobsWindow.h +++ b/viewer/windows/JFJochProcessingJobsWindow.h @@ -63,6 +63,7 @@ private: MergeStatistics merge_stats; double isa = 0.0; bool merge_has_reference = false; + TwinningAnalysisResult twinning; // twinning test of the merged intensities QToolButton *graph_btn = nullptr; // per-row "show statistics" button (enabled once stats exist) }; struct JobSpec { diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index 1da5057f..39164292 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -637,9 +637,13 @@ void NXmx::Sample(const StartMessage &start, const EndMessage &end) { if (!start.sample_name.empty()) group.SaveScalar("name", start.sample_name); - if (start.space_group_number) { - group.SaveScalar("space_group_number", start.space_group_number.value()); - auto *sg = gemmi::find_spacegroup_by_number(start.space_group_number.value()); + // The offline analysis determines the space group only after merging, so it arrives on the end + // message; prefer it over the (usually empty) start-message value. + const auto space_group_number = end.space_group_number ? end.space_group_number + : start.space_group_number; + if (space_group_number) { + group.SaveScalar("space_group_number", space_group_number.value()); + auto *sg = gemmi::find_spacegroup_by_number(space_group_number.value()); if (sg != nullptr) group.SaveScalar("space_group", sg->short_name()); } -- 2.54.0 From 2398330e521743ae9232859cb9e08a7c68937587 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 26 Jun 2026 17:11:41 +0200 Subject: [PATCH 03/66] Drop oversized-metadata frames instead of aborting the collection A serial-crystallography run on a detector with a large converted geometry (JF17T16, modules tiled vertically + horizontally) aborted with "Array out of bounds (Not enough memory to save image)". An indexed still on such a detector predicts/integrates close to the kMaxReflections (10000) cap; at ~170 B per serialized Reflection that is ~1.7 MB of per-image CBOR metadata, which overflowed the fixed 1 MB the buffer slot reserved on top of the image. The serialization guard then threw and cancelled the whole run. - Raise the per-image metadata headroom from 1 MB to 4 MB (GetImageBufferLocationSize). The worst case - 10000 reflections + 2000 spots (API max) + 65534 azimuthal bins - serializes to 2.78 MB, leaving margin while staying negligible next to the multi-MB image slot. - When metadata still does not fit, drop just that frame (log metadata/image/slot sizes + recycle the slot) instead of aborting, in both the FPGA and Lite receivers. - Add a regression test asserting the worst-case metadata fits the headroom. Co-Authored-By: Claude Opus 4.8 (1M context) --- common/DiffractionExperiment.cpp | 8 ++- receiver/JFJochReceiverFPGA.cpp | 77 ++++++++++++++++------------- receiver/JFJochReceiverLite.cpp | 83 +++++++++++++++++++++++--------- tests/CBORTest.cpp | 58 ++++++++++++++++++++++ 4 files changed, 169 insertions(+), 57 deletions(-) diff --git a/common/DiffractionExperiment.cpp b/common/DiffractionExperiment.cpp index 3fac1a9c..51740eed 100644 --- a/common/DiffractionExperiment.cpp +++ b/common/DiffractionExperiment.cpp @@ -1107,7 +1107,13 @@ int64_t DiffractionExperiment::GetImagesPerFile() const { } int64_t DiffractionExperiment::GetImageBufferLocationSize() const { - return GetMaxCompressedSize() + 1024 * 1024; + // A buffer slot holds the compressed image plus the per-image CBOR metadata + // (spot list, reflection list, azimuthal profile, ...). 4 MB of headroom covers + // the worst case: up to kMaxReflections (10000) reflections serialize to ~1.7 MB + // (~170 B each), plus the capped spot list and the azimuthal profile. A frame + // whose metadata still does not fit is dropped individually by the receiver + // rather than aborting the collection. + return GetMaxCompressedSize() + 4 * 1024 * 1024; } float DiffractionExperiment::GetLossyCompressionSerialMX() const { diff --git a/receiver/JFJochReceiverFPGA.cpp b/receiver/JFJochReceiverFPGA.cpp index 470a9096..e4f3749f 100644 --- a/receiver/JFJochReceiverFPGA.cpp +++ b/receiver/JFJochReceiverFPGA.cpp @@ -473,42 +473,51 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) { experiment.GetImageMode(), experiment.GetCompressionAlgorithm()); serializer.SerializeImage(message); - if (experiment.GetImageBufferLocationSize() - serializer.GetImageAppendOffset() - < experiment.GetMaxCompressedSize()) - throw JFJochException(JFJochExceptionCategory::ArrayOutOfBounds, - "Not enough memory to save image"); - const auto compression_start_time = std::chrono::steady_clock::now(); - size_t image_size = transformation.CompressImage(writer_buffer + serializer.GetImageAppendOffset()); - const auto compression_end_time = std::chrono::steady_clock::now(); - message.compression_time_s = std::chrono::duration(compression_end_time - compression_start_time).count(); - - serializer.AppendImage(image_size); - compressed_size += image_size; - if (image_size > 0) - message.compression_ratio = - static_cast(experiment.GetPixelsNum() * experiment.GetByteDepthImage()) - / static_cast(image_size); - - loc->SetImageNumber(image_number); - loc->SetImageSize(serializer.GetBufferSize()); - loc->SetIndexed(message.indexing_result.value_or(false)); - loc->ReadyToSend(); - - if (zmq_preview_socket != nullptr) - zmq_preview_socket->SendImage(writer_buffer, serializer.GetBufferSize()); - - if (zmq_metadata_socket != nullptr) - zmq_metadata_socket->AddDataMessage(message); - - if (push_images_to_writer) { - // Only count images the pusher accepted; TCP can drop on a - // broken/absent connection or an expired enqueue deadline. - if (image_pusher.SendImage(*loc)) - ++images_sent; - } else + const size_t metadata_size = serializer.GetImageAppendOffset(); + const size_t max_image_size = experiment.GetMaxCompressedSize(); + const size_t slot_size = experiment.GetImageBufferLocationSize(); + if (slot_size - metadata_size < max_image_size) { + // The per-image CBOR metadata (spot/reflection lists, ...) leaves + // too little room for the compressed image in the buffer slot. Drop + // just this frame instead of aborting the whole data collection. + logger.Error("Dropping image {} - CBOR metadata too large to store image: " + "metadata {} B + image {} B do not fit in buffer slot {} B", + message.number, metadata_size, max_image_size, slot_size); loc->release(); - UpdateMaxImageSent(message.number); + } else { + const auto compression_start_time = std::chrono::steady_clock::now(); + size_t image_size = transformation.CompressImage(writer_buffer + serializer.GetImageAppendOffset()); + const auto compression_end_time = std::chrono::steady_clock::now(); + message.compression_time_s = std::chrono::duration(compression_end_time - compression_start_time).count(); + + serializer.AppendImage(image_size); + compressed_size += image_size; + if (image_size > 0) + message.compression_ratio = + static_cast(experiment.GetPixelsNum() * experiment.GetByteDepthImage()) + / static_cast(image_size); + + loc->SetImageNumber(image_number); + loc->SetImageSize(serializer.GetBufferSize()); + loc->SetIndexed(message.indexing_result.value_or(false)); + loc->ReadyToSend(); + + if (zmq_preview_socket != nullptr) + zmq_preview_socket->SendImage(writer_buffer, serializer.GetBufferSize()); + + if (zmq_metadata_socket != nullptr) + zmq_metadata_socket->AddDataMessage(message); + + if (push_images_to_writer) { + // Only count images the pusher accepted; TCP can drop on a + // broken/absent connection or an expired enqueue deadline. + if (image_pusher.SendImage(*loc)) + ++images_sent; + } else + loc->release(); + UpdateMaxImageSent(message.number); + } } } diff --git a/receiver/JFJochReceiverLite.cpp b/receiver/JFJochReceiverLite.cpp index b88972e2..cbf004c2 100644 --- a/receiver/JFJochReceiverLite.cpp +++ b/receiver/JFJochReceiverLite.cpp @@ -7,6 +7,41 @@ using namespace std::chrono_literals; +namespace { + // Serialize a processed image into its buffer slot. The slot holds the compressed + // image plus the per-image CBOR metadata (spot list, reflection list, ...); if an + // unusually rich frame makes that metadata too large to fit alongside the image, + // drop just this frame (log + recycle the slot) instead of letting the serializer + // throw, which would abort the whole data collection. Returns true when the image + // was serialized and the slot is ready for the caller to commit. + bool SerializeImageOrDrop(CBORStream2Serializer &serializer, DataMessage &msg, + ZeroCopyReturnValue &loc, size_t slot_size, Logger &logger) { + try { + serializer.SerializeImage(msg); + return true; + } catch (const JFJochException &) { + // Confirm this is a slot-overflow (not some other CBOR error) by re-serializing + // the metadata alone with an empty image placeholder: that always fits when the + // image is what overflowed, and gives the exact metadata size for the log. If it + // throws too, it is not a simple overflow and the exception propagates. + const auto image = msg.image; + msg.image = CompressedImage(nullptr, 0, image.GetWidth(), image.GetHeight(), + image.GetMode(), image.GetCompressionAlgorithm(), + image.GetChannel()); + CBORStream2Serializer meas(static_cast(loc.GetImage()), slot_size); + meas.SerializeImage(msg); + const size_t metadata_size = meas.GetImageAppendOffset(); + msg.image = image; + + logger.Error("Dropping image {} - CBOR metadata too large to store image: " + "metadata {} B + image {} B do not fit in buffer slot {} B", + msg.number, metadata_size, image.GetCompressedSize(), slot_size); + loc.release(); + return false; + } + } +} + int64_t JFJochReceiverLite::NumberOfDataAnalysisThreads(int64_t requested_thread_number, const DiffractionExperiment& in_experiment) { int64_t number_of_images = in_experiment.GetImageNum(); @@ -208,11 +243,13 @@ void JFJochReceiverLite::MaskThread(uint32_t id) { else { auto writer_buffer = (uint8_t *) loc->GetImage(); CBORStream2Serializer serializer(writer_buffer, experiment.GetImageBufferLocationSize()); - serializer.SerializeImage(data_msg); - loc->SetImageNumber(data_msg.number); - loc->SetImageSize(serializer.GetBufferSize()); - loc->SetIndexed(false); - loc->release(); + if (SerializeImageOrDrop(serializer, data_msg, *loc, + experiment.GetImageBufferLocationSize(), logger)) { + loc->SetImageNumber(data_msg.number); + loc->SetImageSize(serializer.GetBufferSize()); + loc->SetIndexed(false); + loc->release(); + } } } current_status.SetProgress(GetProgress()); @@ -315,26 +352,28 @@ void JFJochReceiverLite::DataAnalysisThread(uint32_t id) { else { auto writer_buffer = (uint8_t *) loc->GetImage(); CBORStream2Serializer serializer(writer_buffer, experiment.GetImageBufferLocationSize()); - serializer.SerializeImage(data_msg); - loc->SetImageNumber(data_msg.number); - loc->SetImageSize(serializer.GetBufferSize()); - loc->SetIndexed(data_msg.indexing_result.value_or(false)); - loc->ReadyToSend(); + if (SerializeImageOrDrop(serializer, data_msg, *loc, + experiment.GetImageBufferLocationSize(), logger)) { + loc->SetImageNumber(data_msg.number); + loc->SetImageSize(serializer.GetBufferSize()); + loc->SetIndexed(data_msg.indexing_result.value_or(false)); + loc->ReadyToSend(); - if (zmq_preview_socket != nullptr) - zmq_preview_socket->SendImage(writer_buffer, serializer.GetBufferSize()); + if (zmq_preview_socket != nullptr) + zmq_preview_socket->SendImage(writer_buffer, serializer.GetBufferSize()); - if (zmq_metadata_socket != nullptr) - zmq_metadata_socket->AddDataMessage(data_msg); + if (zmq_metadata_socket != nullptr) + zmq_metadata_socket->AddDataMessage(data_msg); - if (push_images_to_writer) { - // Only count images the pusher accepted; TCP can drop on a - // broken/absent connection or an expired enqueue deadline. - if (image_pusher.SendImage(*loc)) - ++images_sent; - } else - loc->release(); - UpdateMaxImageSent(data_msg.number); + if (push_images_to_writer) { + // Only count images the pusher accepted; TCP can drop on a + // broken/absent connection or an expired enqueue deadline. + if (image_pusher.SendImage(*loc)) + ++images_sent; + } else + loc->release(); + UpdateMaxImageSent(data_msg.number); + } } } diff --git a/tests/CBORTest.cpp b/tests/CBORTest.cpp index 347aef50..9c5e89dc 100644 --- a/tests/CBORTest.cpp +++ b/tests/CBORTest.cpp @@ -997,6 +997,64 @@ TEST_CASE("CBORSerialize_Image_Reflections") { } } +TEST_CASE("CBORSerialize_Image_MetadataHeadroom", "[CBOR]") { + // The receiver serializes the per-image CBOR metadata (spots, reflections, + // azimuthal profile, ...) into the same buffer slot as the compressed image, and + // DiffractionExperiment::GetImageBufferLocationSize() reserves 4 MB of headroom + // on top of GetMaxCompressedSize() for it. A frame whose metadata exceeds that + // headroom is dropped by the receiver. This test guards that even the largest + // possible metadata still fits, so ordinary frames are never dropped: the + // kMaxReflections (10000, IndexAndRefine.cpp) reflection cap, the API maximum of + // 2000 spots (jfjoch_api.yaml), and the 65534-bin azimuthal cap + // (AzimuthalIntegrationMapping, UINT16_MAX - 1), all serialized together. + constexpr size_t max_reflections = 10000; + constexpr size_t max_spots = 2000; + constexpr size_t max_az_bins = 65534; + constexpr size_t headroom = 4 * 1024 * 1024; + + std::vector reflections(max_reflections); + for (size_t i = 0; i < reflections.size(); ++i) + reflections[i] = Reflection{ + .h = static_cast(i % 200) - 100, + .k = static_cast(i % 200) - 100, + .l = static_cast(i % 200) - 100, + .image_number = 789.0f, .delta_phi_deg = 0.123f, + .predicted_x = 1234.5f, .predicted_y = 678.9f, + .observed_x = 1234.0f, .observed_y = 678.0f, + .d = 2.345f, .I = 1234.5f, .bkg = 12.3f, .sigma = 35.6f, + .dist_ewald = 0.001f, .rlp = 0.5f, .partiality = 0.9f, + .zeta = 0.3f, .image_scale_corr = 1.1f}; + + std::vector spots(max_spots); + for (size_t i = 0; i < spots.size(); ++i) + spots[i] = SpotToSave{ + .x = 1234.5f, .y = 678.9f, .intensity = 1234.5f, .maxc = 5000, + .lattice = 0, .image = 789, + .h = static_cast(i % 200) - 100, + .k = static_cast(i % 200) - 100, + .l = static_cast(i % 200) - 100, + .dist_ewald_sphere = 0.01f, .indexed = true}; + + DataMessage message{ + .number = 789, + .image = CompressedImage(nullptr, 0, 4000, 4000, + CompressedImageMode::Uint32, CompressionAlgorithm::BSHUF_LZ4), + .spots = spots, + .az_int_profile = std::vector(max_az_bins, 1.0f), + .az_int_profile_std = std::vector(max_az_bins, 1.0f), + .az_int_profile_count = std::vector(max_az_bins, 1), + .reflections = reflections}; + + std::vector buffer(16 * 1024 * 1024); + CBORStream2Serializer serializer(buffer.data(), buffer.size()); + REQUIRE_NOTHROW(serializer.SerializeImage(message)); + + // GetImageAppendOffset() is where the compressed image is appended, i.e. the size + // of everything-but-the-image; it must stay below the 4 MB metadata headroom. + INFO("worst-case metadata size = " << serializer.GetImageAppendOffset() << " B"); + CHECK(serializer.GetImageAppendOffset() < headroom); +} + TEST_CASE("CBORSerialize_Image_ROI", "[CBOR]") { std::vector buffer(8 * 1024 * 1024); -- 2.54.0 From 8e40822e5b45295a1a07e343e5aa2a3c0f950ef0 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Fri, 26 Jun 2026 15:44:42 +0200 Subject: [PATCH 04/66] Build viewer on Windows --- .gitea/workflows/build_and_test.yml | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/.gitea/workflows/build_and_test.yml b/.gitea/workflows/build_and_test.yml index fa13a3ce..8e33668d 100644 --- a/.gitea/workflows/build_and_test.yml +++ b/.gitea/workflows/build_and_test.yml @@ -64,6 +64,27 @@ jobs: run: | cd build/tools ./jfjoch_hdf5_test ../../tests/test_data/compression_benchmark.h5 + build-windows: + name: build:windows + runs-on: windows-11-cuda-qt + timeout-minutes: 120 + steps: + - uses: actions/checkout@v4 + - name: Configure viewer build + shell: cmd + run: | + for /f "usebackq tokens=*" %%i in (`"%ProgramFiles(x86)%\Microsoft Visual Studio\Installer\vswhere.exe" -latest -property installationPath`) do set "VSPATH=%%i" + call "%VSPATH%\VC\Auxiliary\Build\vcvars64.bat" + cmake -G Ninja -B build -DCMAKE_BUILD_TYPE=Release -DCMAKE_PREFIX_PATH="C:/deps;C:/Qt/6.11.1/msvc2022_64" + - name: Build viewer + shell: cmd + run: | + cmake --build build --target jfjoch_viewer + - name: Build installer (NSIS) + shell: cmd + run: | + cd build + cpack build-rpm: name: build:rpm (${{ matrix.distro }}) if: github.ref_type != 'workflow_dispatch' -- 2.54.0 From 62a7fb3ab65a6348a7e6ca2f044cfc271a599a22 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Fri, 26 Jun 2026 15:49:34 +0200 Subject: [PATCH 05/66] Use vcvars64.bat for all build steps --- .gitea/workflows/build_and_test.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.gitea/workflows/build_and_test.yml b/.gitea/workflows/build_and_test.yml index 8e33668d..b08e9e07 100644 --- a/.gitea/workflows/build_and_test.yml +++ b/.gitea/workflows/build_and_test.yml @@ -79,10 +79,14 @@ jobs: - name: Build viewer shell: cmd run: | + for /f "usebackq tokens=*" %%i in (`"%ProgramFiles(x86)%\Microsoft Visual Studio\Installer\vswhere.exe" -latest -property installationPath`) do set "VSPATH=%%i" + call "%VSPATH%\VC\Auxiliary\Build\vcvars64.bat" cmake --build build --target jfjoch_viewer - name: Build installer (NSIS) shell: cmd run: | + for /f "usebackq tokens=*" %%i in (`"%ProgramFiles(x86)%\Microsoft Visual Studio\Installer\vswhere.exe" -latest -property installationPath`) do set "VSPATH=%%i" + call "%VSPATH%\VC\Auxiliary\Build\vcvars64.bat" cd build cpack build-rpm: -- 2.54.0 From 5aad009cd7d059dfe1cb53b927372badf12d8d05 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Fri, 26 Jun 2026 18:27:33 +0200 Subject: [PATCH 06/66] CI: Build all on Windows (not only jfjoch_viewer target) --- .gitea/workflows/build_and_test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitea/workflows/build_and_test.yml b/.gitea/workflows/build_and_test.yml index b08e9e07..28feda56 100644 --- a/.gitea/workflows/build_and_test.yml +++ b/.gitea/workflows/build_and_test.yml @@ -81,7 +81,7 @@ jobs: run: | for /f "usebackq tokens=*" %%i in (`"%ProgramFiles(x86)%\Microsoft Visual Studio\Installer\vswhere.exe" -latest -property installationPath`) do set "VSPATH=%%i" call "%VSPATH%\VC\Auxiliary\Build\vcvars64.bat" - cmake --build build --target jfjoch_viewer + cmake --build build - name: Build installer (NSIS) shell: cmd run: | -- 2.54.0 From c332e45a54f12c80445ab36147eb00c3b1bd65a7 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Fri, 26 Jun 2026 20:10:52 +0200 Subject: [PATCH 07/66] Compressor: make Compress size-aware, drop frames that don't fit JFJochBitShuffleCompressor::Compress now takes a dest_size and returns a negative value when the compressed output would not fit, instead of writing past the destination buffer. The check is lazy: before each block it verifies the remaining space still covers that block's worst case (mirrored by the new MaxCompressedBlockSize helper, consistent with MaxCompressedSize so a dest sized to MaxCompressedSize never fails). On overflow the dest content is undefined - no rescue. The receiver uses this to compress directly into the writer buffer slot and drop just the oversized frame instead of pre-reserving the full worst-case image size next to the per-image CBOR metadata. Co-Authored-By: Claude Opus 4.8 --- compression/JFJochCompressor.cpp | 38 +++++++++++++++++++++++++++----- compression/JFJochCompressor.h | 6 ++--- receiver/FrameTransformation.cpp | 10 ++++----- receiver/FrameTransformation.h | 2 +- receiver/JFJochReceiverFPGA.cpp | 16 ++++++++------ tests/ZSTDCompressorTest.cpp | 4 ++-- 6 files changed, 53 insertions(+), 23 deletions(-) diff --git a/compression/JFJochCompressor.cpp b/compression/JFJochCompressor.cpp index aa4d06b5..1afa0177 100644 --- a/compression/JFJochCompressor.cpp +++ b/compression/JFJochCompressor.cpp @@ -16,6 +16,20 @@ extern "C" { void bshuf_write_uint64_BE(void* buf, uint64_t num); } +// Worst-case size of one compressed block, including its 4-byte length prefix. Mirrors the +// per-block term of MaxCompressedSize(), so a dest sized to MaxCompressedSize() never fails. +static size_t MaxCompressedBlockSize(CompressionAlgorithm algorithm, size_t src_size) { + switch (algorithm) { + case CompressionAlgorithm::BSHUF_LZ4: + return LZ4_compressBound(src_size) + 4; + case CompressionAlgorithm::BSHUF_ZSTD: + case CompressionAlgorithm::BSHUF_ZSTD_RLE: + return ZSTD_compressBound(src_size) + 4; + default: + return src_size + 4; + } +} + JFJochBitShuffleCompressor::JFJochBitShuffleCompressor(CompressionAlgorithm in_algorithm) { algorithm = in_algorithm; } @@ -59,12 +73,12 @@ size_t JFJochBitShuffleCompressor::CompressBlock(char *dest, const char *source, std::vector JFJochBitShuffleCompressor::Compress(const void *source, size_t nelements, size_t elem_size) { std::vector tmp(MaxCompressedSize(algorithm, nelements, elem_size)); - size_t tmp_size = Compress(tmp.data(), source, nelements, elem_size); + size_t tmp_size = Compress(tmp.data(), tmp.size(), source, nelements, elem_size); tmp.resize(tmp_size); return tmp; } -size_t JFJochBitShuffleCompressor::Compress(void *dest, const void *source, size_t nelements, size_t elem_size) { +int64_t JFJochBitShuffleCompressor::Compress(void *dest, size_t dest_size, const void *source, size_t nelements, size_t elem_size) { auto c_dest = (char *) dest; auto c_source = (char *) source; @@ -72,10 +86,15 @@ size_t JFJochBitShuffleCompressor::Compress(void *dest, const void *source, size if (algorithm == CompressionAlgorithm::NO_COMPRESSION) { // Trivial case if no compression - copy content + if (nelements * elem_size > dest_size) + return -1; memcpy(dest, source, nelements * elem_size); return nelements * elem_size; } + if (dest_size < 12) + return -1; + bshuf_write_uint64_BE(c_dest, nelements * elem_size); bshuf_write_uint32_BE(c_dest + 8, DefaultBlockSize * elem_size); @@ -86,18 +105,27 @@ size_t JFJochBitShuffleCompressor::Compress(void *dest, const void *source, size size_t reminder_size = nelements - num_full_blocks * DefaultBlockSize; size_t compressed_size = 12; - - for (int i = 0; i < num_full_blocks; i++) + // Blocks are small relative to the image, so before each one we just check that the + // remaining space still covers that block's worst case, and bail out (-1) if not. + for (int i = 0; i < num_full_blocks; i++) { + if (compressed_size + MaxCompressedBlockSize(algorithm, DefaultBlockSize * elem_size) > dest_size) + return -1; compressed_size += CompressBlock(c_dest + compressed_size, c_source + i * DefaultBlockSize * elem_size, DefaultBlockSize, elem_size); + } size_t last_block_size = reminder_size - reminder_size % BSHUF_BLOCKED_MULT; - if (last_block_size > 0) + if (last_block_size > 0) { + if (compressed_size + MaxCompressedBlockSize(algorithm, last_block_size * elem_size) > dest_size) + return -1; compressed_size += CompressBlock(c_dest + compressed_size, c_source + num_full_blocks * DefaultBlockSize * elem_size, last_block_size, elem_size); + } size_t leftover_bytes = (reminder_size % BSHUF_BLOCKED_MULT) * elem_size; if (leftover_bytes > 0) { + if (compressed_size + leftover_bytes > dest_size) + return -1; memcpy(c_dest + compressed_size, c_source + (num_full_blocks * DefaultBlockSize + last_block_size) * elem_size, leftover_bytes); compressed_size += leftover_bytes; } diff --git a/compression/JFJochCompressor.h b/compression/JFJochCompressor.h index 0da07ce6..255a49cb 100644 --- a/compression/JFJochCompressor.h +++ b/compression/JFJochCompressor.h @@ -25,8 +25,8 @@ public: explicit JFJochBitShuffleCompressor(CompressionAlgorithm algorithm); template - size_t Compress(void *dest, const std::vector &src) { - return Compress((char *) dest, (char *) src.data(), src.size(), sizeof(T)); + int64_t Compress(void *dest, size_t dest_size, const std::vector &src) { + return Compress(dest, dest_size, src.data(), src.size(), sizeof(T)); }; template @@ -34,7 +34,7 @@ public: return Compress(src.data(), src.size(), sizeof(T)); } std::vector Compress(const void* source, size_t nelements, size_t elem_size); - size_t Compress(void *dest, const void* source, size_t nelements, size_t elem_size); + int64_t Compress(void *dest, size_t dest_size, const void* source, size_t nelements, size_t elem_size); private: char scratch[DefaultBlockSize * sizeof(uint64_t)]; }; diff --git a/receiver/FrameTransformation.cpp b/receiver/FrameTransformation.cpp index aa8f7cfb..82ec2d8c 100644 --- a/receiver/FrameTransformation.cpp +++ b/receiver/FrameTransformation.cpp @@ -54,8 +54,8 @@ FrameTransformation::FrameTransformation(const DiffractionExperiment &in_experim } } -size_t FrameTransformation::CompressImage(void *output) { - return compressor.Compress(output, precompression_buffer.data(), experiment.GetPixelsNum(), +int64_t FrameTransformation::CompressImage(void *output, size_t output_size) { + return compressor.Compress(output, output_size, precompression_buffer.data(), experiment.GetPixelsNum(), experiment.GetByteDepthImage()); } @@ -63,14 +63,14 @@ CompressedImage FrameTransformation::GetCompressedImage() { const uint8_t *data; size_t size; if (experiment.GetCompressionAlgorithm() == CompressionAlgorithm::NO_COMPRESSION) { - data = (uint8_t *) precompression_buffer.data(); + data = reinterpret_cast(precompression_buffer.data()); size = experiment.GetPixelsNum() * experiment.GetByteDepthImage(); } else { compressed_buffer.resize(MaxCompressedSize(experiment.GetCompressionAlgorithm(), experiment.GetPixelsNum(), experiment.GetByteDepthImage())); - data = (uint8_t *) compressed_buffer.data(); - size = CompressImage(compressed_buffer.data()); + data = reinterpret_cast(compressed_buffer.data()); + size = static_cast(CompressImage(compressed_buffer.data(), compressed_buffer.size())); } return CompressedImage(data, size,experiment.GetXPixelsNum(), diff --git a/receiver/FrameTransformation.h b/receiver/FrameTransformation.h index e155f013..870f52e5 100644 --- a/receiver/FrameTransformation.h +++ b/receiver/FrameTransformation.h @@ -24,6 +24,6 @@ public: void ProcessModule(const void *input, uint16_t module_number, int data_stream); void FillNotCollectedModule(uint16_t module_number, int data_stream); CompressedImage GetCompressedImage(); - size_t CompressImage(void *output); + int64_t CompressImage(void *output, size_t output_size); // < 0 - error in Compression const void *GetImage() const; }; diff --git a/receiver/JFJochReceiverFPGA.cpp b/receiver/JFJochReceiverFPGA.cpp index e4f3749f..6fe2a27e 100644 --- a/receiver/JFJochReceiverFPGA.cpp +++ b/receiver/JFJochReceiverFPGA.cpp @@ -475,21 +475,23 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) { serializer.SerializeImage(message); const size_t metadata_size = serializer.GetImageAppendOffset(); - const size_t max_image_size = experiment.GetMaxCompressedSize(); const size_t slot_size = experiment.GetImageBufferLocationSize(); - if (slot_size - metadata_size < max_image_size) { + + const auto compression_start_time = std::chrono::steady_clock::now(); + int64_t image_size = transformation.CompressImage(writer_buffer + serializer.GetImageAppendOffset(), + slot_size - (metadata_size + 32)); // keep 32 byte extra for close, etc. + if (image_size < 0) { // The per-image CBOR metadata (spot/reflection lists, ...) leaves // too little room for the compressed image in the buffer slot. Drop // just this frame instead of aborting the whole data collection. logger.Error("Dropping image {} - CBOR metadata too large to store image: " - "metadata {} B + image {} B do not fit in buffer slot {} B", - message.number, metadata_size, max_image_size, slot_size); + "metadata {} do not fit in buffer slot {} B", + message.number, metadata_size, slot_size); loc->release(); } else { - const auto compression_start_time = std::chrono::steady_clock::now(); - size_t image_size = transformation.CompressImage(writer_buffer + serializer.GetImageAppendOffset()); const auto compression_end_time = std::chrono::steady_clock::now(); - message.compression_time_s = std::chrono::duration(compression_end_time - compression_start_time).count(); + message.compression_time_s = std::chrono::duration( + compression_end_time - compression_start_time).count(); serializer.AppendImage(image_size); compressed_size += image_size; diff --git a/tests/ZSTDCompressorTest.cpp b/tests/ZSTDCompressorTest.cpp index 699f7ce5..0d0d4262 100644 --- a/tests/ZSTDCompressorTest.cpp +++ b/tests/ZSTDCompressorTest.cpp @@ -240,7 +240,7 @@ TEST_CASE("JFJochCompressor_JFJochDecompressor_ZSTD","[ZSTD]") { std::vector tmp(x.GetPixelsNum() * sizeof(int32_t) * 4 + 12); - auto tmp_size = compressor.Compress(tmp.data(), image); + auto tmp_size = compressor.Compress(tmp.data(), tmp.size(), image); tmp.resize(tmp_size); std::vector output; @@ -262,7 +262,7 @@ TEST_CASE("JFJochCompressor_JFJochDecompressor_LZ4","[ZSTD]") { std::vector tmp(x.GetPixelsNum() * sizeof(int32_t) * 4 + 12); - auto tmp_size = compressor.Compress(tmp.data(), image); + auto tmp_size = compressor.Compress(tmp.data(), tmp.size(), image); tmp.resize(tmp_size); std::vector output; -- 2.54.0 From aadba5b34304429ad9eaf57df3cdeb5180426476 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Fri, 26 Jun 2026 20:35:32 +0200 Subject: [PATCH 08/66] Compressor: bump bitshuffle block size 4096 -> 16384 elements On sparse lyso frames the larger block improves compression ratio across all bshuf algorithms (16-bit data): ZSTD 8.58 -> 9.30, LZ4 7.38 -> 7.58, RLE 6.82 -> 6.90. 16384 captures most of the gain available from even larger blocks (ZSTD tops out ~9.55 at 65536) while staying close to the cache sweet spot: the cheap codecs (LZ4, RLE) peak in throughput once a block's working set fits L1d (~4096 elem here), so very large blocks trade real throughput for diminishing ratio - and that penalty is worse on the Xeon Gold/Platinum production hosts (smaller private L2, shared-L3 contention under many parallel compression threads). The block size is stored per-dataset in the bitshuffle HDF5 filter params, so existing readers (XDS/Neggia/Durin/CrystFEL) stay compatible. Move the per-block bitshuffle scratch off the inline member array onto a lazily-sized heap vector, like tmp_space, so the block size no longer bloats every stack-allocated compressor (incl. the transient ones in CBORStream2Serializer). Co-Authored-By: Claude Opus 4.8 --- compression/JFJochCompressor.cpp | 4 +++- compression/JFJochCompressor.h | 5 ++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/compression/JFJochCompressor.cpp b/compression/JFJochCompressor.cpp index 1afa0177..c144abb1 100644 --- a/compression/JFJochCompressor.cpp +++ b/compression/JFJochCompressor.cpp @@ -37,7 +37,7 @@ JFJochBitShuffleCompressor::JFJochBitShuffleCompressor(CompressionAlgorithm in_a size_t JFJochBitShuffleCompressor::CompressBlock(char *dest, const char *source, size_t nelements, size_t elem_size) { // Assert nelements < block_size const char *src_ptr; - int64_t bshuf_ret = bitshuf_encode_block(tmp_space.data(), source, scratch, nelements, elem_size); + int64_t bshuf_ret = bitshuf_encode_block(tmp_space.data(), source, scratch.data(), nelements, elem_size); if (bshuf_ret < 0) throw JFJochException(JFJochExceptionCategory::Compression, "bshuf_trans_bit_elem error"); src_ptr = tmp_space.data(); @@ -100,6 +100,8 @@ int64_t JFJochBitShuffleCompressor::Compress(void *dest, size_t dest_size, const if (tmp_space.size() < DefaultBlockSize * elem_size) tmp_space.resize(DefaultBlockSize * elem_size); + if (scratch.size() < DefaultBlockSize * elem_size) + scratch.resize(DefaultBlockSize * elem_size); size_t num_full_blocks = nelements / DefaultBlockSize; size_t reminder_size = nelements - num_full_blocks * DefaultBlockSize; diff --git a/compression/JFJochCompressor.h b/compression/JFJochCompressor.h index 255a49cb..30bcbcc7 100644 --- a/compression/JFJochCompressor.h +++ b/compression/JFJochCompressor.h @@ -17,10 +17,11 @@ class JFJochBitShuffleCompressor { JFJochZstdCompressor zstd_compressor; CompressionAlgorithm algorithm; std::vector tmp_space; + std::vector scratch; size_t CompressBlock(char *dest, const char * source, size_t nelements, size_t elem_size); public: - constexpr static const size_t DefaultBlockSize = 4096; + constexpr static const size_t DefaultBlockSize = 16384; explicit JFJochBitShuffleCompressor(CompressionAlgorithm algorithm); @@ -35,8 +36,6 @@ public: } std::vector Compress(const void* source, size_t nelements, size_t elem_size); int64_t Compress(void *dest, size_t dest_size, const void* source, size_t nelements, size_t elem_size); -private: - char scratch[DefaultBlockSize * sizeof(uint64_t)]; }; template std::vector bitshuffle(const std::vector &input, size_t block_size) { -- 2.54.0 From 74584c23ac290b3720dcd459abcc8336b03c6949 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 27 Jun 2026 07:56:08 +0200 Subject: [PATCH 09/66] Compressor: scale block size to a per-algorithm byte target Replace the fixed-element DefaultBlockSize with a byte target divided by elem_size to get the block element count, so the per-block working set (and thus cache behaviour) stays constant across pixel bit depths instead of halving from 8- to 16- to 32-bit. The target is per-algorithm, following the measured sweet spots on sparse data: LZ4 wants a small, cache-resident block for throughput (16 kB), ZSTD/RLE want a large block for ratio (128 kB). The gap is widest on extreme-sparsity inputs such as the uint32 pixel_mask, where large-block ZSTD reaches 100-1800x vs ~160x for LZ4. The block size is read back per-dataset from the bitshuffle stream header (block_size = header_bytes / elem_size) and the HDF5 filter params, so the decompressor and external readers (XDS/Neggia/Durin/CrystFEL) need no change. Co-Authored-By: Claude Opus 4.8 --- compression/JFJochCompressor.cpp | 33 +++++++++++++++++++------------ compression/JFJochCompressor.h | 13 +++++++++++- compression/MaxCompressedSize.cpp | 4 ++-- tests/FrameTransformationTest.cpp | 21 +++++++++++++------- writer/HDF5DataFile.cpp | 4 +++- 5 files changed, 51 insertions(+), 24 deletions(-) diff --git a/compression/JFJochCompressor.cpp b/compression/JFJochCompressor.cpp index c144abb1..4540db96 100644 --- a/compression/JFJochCompressor.cpp +++ b/compression/JFJochCompressor.cpp @@ -16,6 +16,13 @@ extern "C" { void bshuf_write_uint64_BE(void* buf, uint64_t num); } +// Necessary condition for BlockSize() to be a valid bitshuffle block: with elem_size a power of +// two, a byte target that is a multiple of BSHUF_BLOCKED_MULT keeps the element count a multiple +// of BSHUF_BLOCKED_MULT too. +static_assert(JFJochBitShuffleCompressor::DefaultBlockSizeBytes(CompressionAlgorithm::BSHUF_LZ4) % BSHUF_BLOCKED_MULT == 0 + && JFJochBitShuffleCompressor::DefaultBlockSizeBytes(CompressionAlgorithm::BSHUF_ZSTD) % BSHUF_BLOCKED_MULT == 0, + "block byte target must be a multiple of the bitshuffle block multiple"); + // Worst-case size of one compressed block, including its 4-byte length prefix. Mirrors the // per-block term of MaxCompressedSize(), so a dest sized to MaxCompressedSize() never fails. static size_t MaxCompressedBlockSize(CompressionAlgorithm algorithm, size_t src_size) { @@ -82,8 +89,6 @@ int64_t JFJochBitShuffleCompressor::Compress(void *dest, size_t dest_size, const auto c_dest = (char *) dest; auto c_source = (char *) source; - static_assert(DefaultBlockSize % BSHUF_BLOCKED_MULT == 0, "Block size must be multiple of 8"); - if (algorithm == CompressionAlgorithm::NO_COMPRESSION) { // Trivial case if no compression - copy content if (nelements * elem_size > dest_size) @@ -95,25 +100,27 @@ int64_t JFJochBitShuffleCompressor::Compress(void *dest, size_t dest_size, const if (dest_size < 12) return -1; + const size_t block_size = BlockSize(algorithm, elem_size); + bshuf_write_uint64_BE(c_dest, nelements * elem_size); - bshuf_write_uint32_BE(c_dest + 8, DefaultBlockSize * elem_size); + bshuf_write_uint32_BE(c_dest + 8, block_size * elem_size); - if (tmp_space.size() < DefaultBlockSize * elem_size) - tmp_space.resize(DefaultBlockSize * elem_size); - if (scratch.size() < DefaultBlockSize * elem_size) - scratch.resize(DefaultBlockSize * elem_size); + if (tmp_space.size() < block_size * elem_size) + tmp_space.resize(block_size * elem_size); + if (scratch.size() < block_size * elem_size) + scratch.resize(block_size * elem_size); - size_t num_full_blocks = nelements / DefaultBlockSize; - size_t reminder_size = nelements - num_full_blocks * DefaultBlockSize; + size_t num_full_blocks = nelements / block_size; + size_t reminder_size = nelements - num_full_blocks * block_size; size_t compressed_size = 12; // Blocks are small relative to the image, so before each one we just check that the // remaining space still covers that block's worst case, and bail out (-1) if not. for (int i = 0; i < num_full_blocks; i++) { - if (compressed_size + MaxCompressedBlockSize(algorithm, DefaultBlockSize * elem_size) > dest_size) + if (compressed_size + MaxCompressedBlockSize(algorithm, block_size * elem_size) > dest_size) return -1; compressed_size += CompressBlock(c_dest + compressed_size, - c_source + i * DefaultBlockSize * elem_size, DefaultBlockSize, elem_size); + c_source + i * block_size * elem_size, block_size, elem_size); } size_t last_block_size = reminder_size - reminder_size % BSHUF_BLOCKED_MULT; @@ -121,14 +128,14 @@ int64_t JFJochBitShuffleCompressor::Compress(void *dest, size_t dest_size, const if (compressed_size + MaxCompressedBlockSize(algorithm, last_block_size * elem_size) > dest_size) return -1; compressed_size += CompressBlock(c_dest + compressed_size, - c_source + num_full_blocks * DefaultBlockSize * elem_size, last_block_size, elem_size); + c_source + num_full_blocks * block_size * elem_size, last_block_size, elem_size); } size_t leftover_bytes = (reminder_size % BSHUF_BLOCKED_MULT) * elem_size; if (leftover_bytes > 0) { if (compressed_size + leftover_bytes > dest_size) return -1; - memcpy(c_dest + compressed_size, c_source + (num_full_blocks * DefaultBlockSize + last_block_size) * elem_size, leftover_bytes); + memcpy(c_dest + compressed_size, c_source + (num_full_blocks * block_size + last_block_size) * elem_size, leftover_bytes); compressed_size += leftover_bytes; } return compressed_size; diff --git a/compression/JFJochCompressor.h b/compression/JFJochCompressor.h index 30bcbcc7..a6a3f97c 100644 --- a/compression/JFJochCompressor.h +++ b/compression/JFJochCompressor.h @@ -21,7 +21,18 @@ class JFJochBitShuffleCompressor { size_t CompressBlock(char *dest, const char * source, size_t nelements, size_t elem_size); public: - constexpr static const size_t DefaultBlockSize = 16384; + // The bitshuffle block size is chosen as a byte target rather than a fixed element count, so + // the per-block working set - and thus the cache behaviour - stays constant across pixel bit + // depths. The target is per-algorithm: LZ4 favours a small, cache-resident block (throughput), + // ZSTD/RLE a large block (ratio). The element block size scales inversely with elem_size, + // which is always a power of two here, so the result is a power of two >= 8 (a valid bitshuffle + // block, multiple of 8). E.g. ZSTD 128 kB -> 65536 elem at 16-bit, 32768 at 32-bit. + constexpr static size_t DefaultBlockSizeBytes(CompressionAlgorithm algorithm) { + return algorithm == CompressionAlgorithm::BSHUF_LZ4 ? 16384 : 131072; + } + constexpr static size_t BlockSize(CompressionAlgorithm algorithm, size_t elem_size) { + return DefaultBlockSizeBytes(algorithm) / elem_size; + } explicit JFJochBitShuffleCompressor(CompressionAlgorithm algorithm); diff --git a/compression/MaxCompressedSize.cpp b/compression/MaxCompressedSize.cpp index 2d2b9243..613a3480 100644 --- a/compression/MaxCompressedSize.cpp +++ b/compression/MaxCompressedSize.cpp @@ -9,10 +9,10 @@ int64_t MaxCompressedSize(CompressionAlgorithm algorithm, int64_t pixels_number, uint16_t pixel_depth) { switch (algorithm) { case CompressionAlgorithm::BSHUF_LZ4: - return bshuf_compress_lz4_bound(pixels_number, pixel_depth, JFJochBitShuffleCompressor::DefaultBlockSize) + 12; + return bshuf_compress_lz4_bound(pixels_number, pixel_depth, JFJochBitShuffleCompressor::BlockSize(algorithm, pixel_depth)) + 12; case CompressionAlgorithm::BSHUF_ZSTD: case CompressionAlgorithm::BSHUF_ZSTD_RLE: - return bshuf_compress_zstd_bound(pixels_number, pixel_depth,JFJochBitShuffleCompressor::DefaultBlockSize) + 12; + return bshuf_compress_zstd_bound(pixels_number, pixel_depth, JFJochBitShuffleCompressor::BlockSize(algorithm, pixel_depth)) + 12; default: return pixels_number * pixel_depth; } diff --git a/tests/FrameTransformationTest.cpp b/tests/FrameTransformationTest.cpp index 1ac5d701..22f99752 100644 --- a/tests/FrameTransformationTest.cpp +++ b/tests/FrameTransformationTest.cpp @@ -97,7 +97,8 @@ TEST_CASE("FrameTransformation_Raw_NoCompression_bshuf_lz4" ,"") { auto image = transformation.GetCompressedImage(); REQUIRE(read_be64(image.GetCompressed()) == experiment.GetPixelsNum() * experiment.GetByteDepthImage()); - REQUIRE(read_be32(image.GetCompressed() + 8) == JFJochBitShuffleCompressor::DefaultBlockSize * + REQUIRE(read_be32(image.GetCompressed() + 8) == JFJochBitShuffleCompressor::BlockSize(experiment.GetCompressionAlgorithm(), + experiment.GetByteDepthImage()) * experiment.GetByteDepthImage()); std::vector output; @@ -238,7 +239,8 @@ TEST_CASE("FrameTransformation_Converted_bshuf_lz4" ,"") { auto image = transformation.GetCompressedImage(); REQUIRE(read_be64(image.GetCompressed()) == experiment.GetPixelsNum() * experiment.GetByteDepthImage()); - REQUIRE(read_be32(image.GetCompressed() + 8) == JFJochBitShuffleCompressor::DefaultBlockSize * + REQUIRE(read_be32(image.GetCompressed() + 8) == JFJochBitShuffleCompressor::BlockSize(experiment.GetCompressionAlgorithm(), + experiment.GetByteDepthImage()) * experiment.GetByteDepthImage()); std::vector output; @@ -297,7 +299,8 @@ TEST_CASE("FrameTransformation_Converted_bshuf_zstd" ,"") { auto image = transformation.GetCompressedImage(); REQUIRE(read_be64(image.GetCompressed()) == experiment.GetPixelsNum() * experiment.GetByteDepthImage()); - REQUIRE(read_be32(image.GetCompressed() + 8) == JFJochBitShuffleCompressor::DefaultBlockSize * + REQUIRE(read_be32(image.GetCompressed() + 8) == JFJochBitShuffleCompressor::BlockSize(experiment.GetCompressionAlgorithm(), + experiment.GetByteDepthImage()) * experiment.GetByteDepthImage()); std::vector output; @@ -356,7 +359,8 @@ TEST_CASE("FrameTransformation_Converted_bshuf_zstd_rle" ,"") { auto image = transformation.GetCompressedImage(); REQUIRE(read_be64(image.GetCompressed()) == experiment.GetPixelsNum() * experiment.GetByteDepthImage()); - REQUIRE(read_be32(image.GetCompressed() + 8) == JFJochBitShuffleCompressor::DefaultBlockSize * + REQUIRE(read_be32(image.GetCompressed() + 8) == JFJochBitShuffleCompressor::BlockSize(experiment.GetCompressionAlgorithm(), + experiment.GetByteDepthImage()) * experiment.GetByteDepthImage()); std::vector output; @@ -418,7 +422,8 @@ TEST_CASE("FrameTransformation_Converted_bshuf_zstd_32bit" ,"") { auto image = transformation.GetCompressedImage(); REQUIRE(read_be64(image.GetCompressed()) == experiment.GetPixelsNum() * experiment.GetByteDepthImage()); - REQUIRE(read_be32(image.GetCompressed() + 8) == JFJochBitShuffleCompressor::DefaultBlockSize * + REQUIRE(read_be32(image.GetCompressed() + 8) == JFJochBitShuffleCompressor::BlockSize(experiment.GetCompressionAlgorithm(), + experiment.GetByteDepthImage()) * experiment.GetByteDepthImage()); std::vector output; @@ -481,7 +486,8 @@ TEST_CASE("FrameTransformation_Converted_bshuf_zstd_8bit" ,"") { auto image = transformation.GetCompressedImage(); REQUIRE(read_be64(image.GetCompressed()) == experiment.GetPixelsNum() * experiment.GetByteDepthImage()); - REQUIRE(read_be32(image.GetCompressed() + 8) == JFJochBitShuffleCompressor::DefaultBlockSize * + REQUIRE(read_be32(image.GetCompressed() + 8) == JFJochBitShuffleCompressor::BlockSize(experiment.GetCompressionAlgorithm(), + experiment.GetByteDepthImage()) * experiment.GetByteDepthImage()); std::vector output; @@ -548,7 +554,8 @@ TEST_CASE("FrameTransformation_Converted_bshuf_zstd_unsigned_16bit" ,"") { auto image = transformation.GetCompressedImage(); REQUIRE(read_be64(image.GetCompressed()) == experiment.GetPixelsNum() * experiment.GetByteDepthImage()); - REQUIRE(read_be32(image.GetCompressed() + 8) == JFJochBitShuffleCompressor::DefaultBlockSize * + REQUIRE(read_be32(image.GetCompressed() + 8) == JFJochBitShuffleCompressor::BlockSize(experiment.GetCompressionAlgorithm(), + experiment.GetByteDepthImage()) * experiment.GetByteDepthImage()); std::vector output; diff --git a/writer/HDF5DataFile.cpp b/writer/HDF5DataFile.cpp index 378b1ca3..88fc38ed 100644 --- a/writer/HDF5DataFile.cpp +++ b/writer/HDF5DataFile.cpp @@ -140,7 +140,9 @@ void HDF5DataFile::CreateFile(const DataMessage& msg, std::shared_ptr xpixel = msg.image.GetWidth(); ypixel = msg.image.GetHeight(); - dcpl.SetCompression(msg.image.GetCompressionAlgorithm(), JFJochBitShuffleCompressor::DefaultBlockSize); + dcpl.SetCompression(msg.image.GetCompressionAlgorithm(), + JFJochBitShuffleCompressor::BlockSize(msg.image.GetCompressionAlgorithm(), + msg.image.GetByteDepth())); dcpl.SetChunking( {1, ypixel, xpixel}); H5Pset_fill_time(dcpl.GetID(), H5D_FILL_TIME_NEVER); -- 2.54.0 From 02514bb6b733d8738fc3014d75eb26712dc90200 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 27 Jun 2026 08:17:21 +0200 Subject: [PATCH 10/66] Compressor: throw on overflow instead of returning a negative size Compress() and FrameTransformation::CompressImage() returned int64_t with a negative value meaning "did not fit". That is a footgun: the negative result silently converts to a huge size_t if a caller forgets to check it. Return size_t and instead throw a named CompressionBufferTooSmallException (deriving from JFJochException, Compression category) when the output would not fit the destination buffer. The receiver catches it explicitly and drops just that frame, as before; the offline/GetCompressedImage path uses a worst-case buffer so it never throws. Add a test that a too-small destination throws and a worst-case buffer does not. Co-Authored-By: Claude Opus 4.8 --- common/JFJochException.h | 8 ++++++++ compression/JFJochCompressor.cpp | 14 +++++++------- compression/JFJochCompressor.h | 5 +++-- receiver/FrameTransformation.cpp | 4 ++-- receiver/FrameTransformation.h | 3 ++- receiver/JFJochReceiverFPGA.cpp | 23 ++++++++++++----------- tests/ZSTDCompressorTest.cpp | 18 ++++++++++++++++++ 7 files changed, 52 insertions(+), 23 deletions(-) diff --git a/common/JFJochException.h b/common/JFJochException.h index d2ecdb0f..4c132a5a 100644 --- a/common/JFJochException.h +++ b/common/JFJochException.h @@ -162,3 +162,11 @@ public: msg += " (" + std::string(strerror(errno)) + ")"; } }; + +// Thrown by JFJochBitShuffleCompressor::Compress when the compressed output does not fit in the +// caller-provided destination buffer. Lets callers (e.g. the receiver) drop just that frame. +class CompressionBufferTooSmallException : public JFJochException { +public: + explicit CompressionBufferTooSmallException(const std::string &description) + : JFJochException(JFJochExceptionCategory::Compression, description) {} +}; diff --git a/compression/JFJochCompressor.cpp b/compression/JFJochCompressor.cpp index 4540db96..3470e728 100644 --- a/compression/JFJochCompressor.cpp +++ b/compression/JFJochCompressor.cpp @@ -85,20 +85,20 @@ std::vector JFJochBitShuffleCompressor::Compress(const void *source, si return tmp; } -int64_t JFJochBitShuffleCompressor::Compress(void *dest, size_t dest_size, const void *source, size_t nelements, size_t elem_size) { +size_t JFJochBitShuffleCompressor::Compress(void *dest, size_t dest_size, const void *source, size_t nelements, size_t elem_size) { auto c_dest = (char *) dest; auto c_source = (char *) source; if (algorithm == CompressionAlgorithm::NO_COMPRESSION) { // Trivial case if no compression - copy content if (nelements * elem_size > dest_size) - return -1; + throw CompressionBufferTooSmallException("compressed output exceeds destination buffer"); memcpy(dest, source, nelements * elem_size); return nelements * elem_size; } if (dest_size < 12) - return -1; + throw CompressionBufferTooSmallException("compressed output exceeds destination buffer"); const size_t block_size = BlockSize(algorithm, elem_size); @@ -115,10 +115,10 @@ int64_t JFJochBitShuffleCompressor::Compress(void *dest, size_t dest_size, const size_t compressed_size = 12; // Blocks are small relative to the image, so before each one we just check that the - // remaining space still covers that block's worst case, and bail out (-1) if not. + // remaining space still covers that block's worst case, and throw if not. for (int i = 0; i < num_full_blocks; i++) { if (compressed_size + MaxCompressedBlockSize(algorithm, block_size * elem_size) > dest_size) - return -1; + throw CompressionBufferTooSmallException("compressed output exceeds destination buffer"); compressed_size += CompressBlock(c_dest + compressed_size, c_source + i * block_size * elem_size, block_size, elem_size); } @@ -126,7 +126,7 @@ int64_t JFJochBitShuffleCompressor::Compress(void *dest, size_t dest_size, const size_t last_block_size = reminder_size - reminder_size % BSHUF_BLOCKED_MULT; if (last_block_size > 0) { if (compressed_size + MaxCompressedBlockSize(algorithm, last_block_size * elem_size) > dest_size) - return -1; + throw CompressionBufferTooSmallException("compressed output exceeds destination buffer"); compressed_size += CompressBlock(c_dest + compressed_size, c_source + num_full_blocks * block_size * elem_size, last_block_size, elem_size); } @@ -134,7 +134,7 @@ int64_t JFJochBitShuffleCompressor::Compress(void *dest, size_t dest_size, const size_t leftover_bytes = (reminder_size % BSHUF_BLOCKED_MULT) * elem_size; if (leftover_bytes > 0) { if (compressed_size + leftover_bytes > dest_size) - return -1; + throw CompressionBufferTooSmallException("compressed output exceeds destination buffer"); memcpy(c_dest + compressed_size, c_source + (num_full_blocks * block_size + last_block_size) * elem_size, leftover_bytes); compressed_size += leftover_bytes; } diff --git a/compression/JFJochCompressor.h b/compression/JFJochCompressor.h index a6a3f97c..644900d2 100644 --- a/compression/JFJochCompressor.h +++ b/compression/JFJochCompressor.h @@ -37,7 +37,7 @@ public: explicit JFJochBitShuffleCompressor(CompressionAlgorithm algorithm); template - int64_t Compress(void *dest, size_t dest_size, const std::vector &src) { + size_t Compress(void *dest, size_t dest_size, const std::vector &src) { return Compress(dest, dest_size, src.data(), src.size(), sizeof(T)); }; @@ -46,7 +46,8 @@ public: return Compress(src.data(), src.size(), sizeof(T)); } std::vector Compress(const void* source, size_t nelements, size_t elem_size); - int64_t Compress(void *dest, size_t dest_size, const void* source, size_t nelements, size_t elem_size); + // Throws CompressionBufferTooSmallException if the compressed output would not fit dest_size. + size_t Compress(void *dest, size_t dest_size, const void* source, size_t nelements, size_t elem_size); }; template std::vector bitshuffle(const std::vector &input, size_t block_size) { diff --git a/receiver/FrameTransformation.cpp b/receiver/FrameTransformation.cpp index 82ec2d8c..1f824217 100644 --- a/receiver/FrameTransformation.cpp +++ b/receiver/FrameTransformation.cpp @@ -54,7 +54,7 @@ FrameTransformation::FrameTransformation(const DiffractionExperiment &in_experim } } -int64_t FrameTransformation::CompressImage(void *output, size_t output_size) { +size_t FrameTransformation::CompressImage(void *output, size_t output_size) { return compressor.Compress(output, output_size, precompression_buffer.data(), experiment.GetPixelsNum(), experiment.GetByteDepthImage()); } @@ -70,7 +70,7 @@ CompressedImage FrameTransformation::GetCompressedImage() { experiment.GetPixelsNum(), experiment.GetByteDepthImage())); data = reinterpret_cast(compressed_buffer.data()); - size = static_cast(CompressImage(compressed_buffer.data(), compressed_buffer.size())); + size = CompressImage(compressed_buffer.data(), compressed_buffer.size()); } return CompressedImage(data, size,experiment.GetXPixelsNum(), diff --git a/receiver/FrameTransformation.h b/receiver/FrameTransformation.h index 870f52e5..586a37ad 100644 --- a/receiver/FrameTransformation.h +++ b/receiver/FrameTransformation.h @@ -24,6 +24,7 @@ public: void ProcessModule(const void *input, uint16_t module_number, int data_stream); void FillNotCollectedModule(uint16_t module_number, int data_stream); CompressedImage GetCompressedImage(); - int64_t CompressImage(void *output, size_t output_size); // < 0 - error in Compression + // Throws CompressionBufferTooSmallException if the image does not fit output_size. + size_t CompressImage(void *output, size_t output_size); const void *GetImage() const; }; diff --git a/receiver/JFJochReceiverFPGA.cpp b/receiver/JFJochReceiverFPGA.cpp index 6fe2a27e..5ba633b2 100644 --- a/receiver/JFJochReceiverFPGA.cpp +++ b/receiver/JFJochReceiverFPGA.cpp @@ -478,17 +478,10 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) { const size_t slot_size = experiment.GetImageBufferLocationSize(); const auto compression_start_time = std::chrono::steady_clock::now(); - int64_t image_size = transformation.CompressImage(writer_buffer + serializer.GetImageAppendOffset(), - slot_size - (metadata_size + 32)); // keep 32 byte extra for close, etc. - if (image_size < 0) { - // The per-image CBOR metadata (spot/reflection lists, ...) leaves - // too little room for the compressed image in the buffer slot. Drop - // just this frame instead of aborting the whole data collection. - logger.Error("Dropping image {} - CBOR metadata too large to store image: " - "metadata {} do not fit in buffer slot {} B", - message.number, metadata_size, slot_size); - loc->release(); - } else { + try { + size_t image_size = transformation.CompressImage( + writer_buffer + serializer.GetImageAppendOffset(), + slot_size - (metadata_size + 32)); // keep 32 byte extra for close, etc. const auto compression_end_time = std::chrono::steady_clock::now(); message.compression_time_s = std::chrono::duration( compression_end_time - compression_start_time).count(); @@ -519,6 +512,14 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) { } else loc->release(); UpdateMaxImageSent(message.number); + } catch (const CompressionBufferTooSmallException &) { + // The per-image CBOR metadata (spot/reflection lists, ...) leaves + // too little room for the compressed image in the buffer slot. Drop + // just this frame instead of aborting the whole data collection. + logger.Error("Dropping image {} - CBOR metadata too large to store image: " + "metadata {} do not fit in buffer slot {} B", + message.number, metadata_size, slot_size); + loc->release(); } } } diff --git a/tests/ZSTDCompressorTest.cpp b/tests/ZSTDCompressorTest.cpp index 0d0d4262..ba2fe8a4 100644 --- a/tests/ZSTDCompressorTest.cpp +++ b/tests/ZSTDCompressorTest.cpp @@ -249,6 +249,24 @@ TEST_CASE("JFJochCompressor_JFJochDecompressor_ZSTD","[ZSTD]") { REQUIRE(memcmp(image.data(), output.data(), x.GetPixelsNum() * sizeof(int32_t)) == 0); } +TEST_CASE("JFJochCompressor_DestTooSmall_Throws","[ZSTD]") { + DiffractionExperiment x(DetJF4M()); + x.Compression(CompressionAlgorithm::BSHUF_ZSTD).BitDepthImage(32).PixelSigned(true); + + std::vector image(x.GetPixelsNum(), 345); + + JFJochBitShuffleCompressor compressor(x.GetCompressionAlgorithm()); + + // A destination far too small for the compressed output must throw, not overflow it. + std::vector tiny(64); + REQUIRE_THROWS_AS(compressor.Compress(tiny.data(), tiny.size(), image), + CompressionBufferTooSmallException); + + // A buffer sized to the worst case never throws. + std::vector big(MaxCompressedSize(x.GetCompressionAlgorithm(), x.GetPixelsNum(), sizeof(int32_t))); + REQUIRE_NOTHROW(compressor.Compress(big.data(), big.size(), image)); +} + TEST_CASE("JFJochCompressor_JFJochDecompressor_LZ4","[ZSTD]") { DiffractionExperiment x(DetJF4M()); x.Compression(CompressionAlgorithm::BSHUF_LZ4).BitDepthImage(32).PixelSigned(true); -- 2.54.0 From 7e7a73062c16862df78d44521ecd9f97c168605c Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 27 Jun 2026 14:41:46 +0200 Subject: [PATCH 11/66] Compression: add BSHUF_ZSTD_RLE_HUFF (RLE runs + Huffman literals) New CompressionAlgorithm that emits a standard Zstandard frame: zero/0xFF runs become RLE_Blocks (like BSHUF_ZSTD_RLE) and literal regions become Compressed_Blocks with per-block adaptive Huffman literals and no sequences (Number_of_Sequences=0). Short runs are absorbed into the literal stream; incompressible literals fall back to Raw_Blocks so the worst case stays within ZSTD_compressBound. The Huffman tree + bitstream are produced by zstd's own HUF_compress{1,4}X_repeat (the same calls ZSTD_compressLiterals uses); only the frame/block/literals-section framing is hand-written, with comments citing zstd_compression_format.md so it can be checked clause by clause. Output decodes with stock ZSTD_decompress, so no reader changes are needed (decode routes like BSHUF_ZSTD). On sparse diffraction this gives ~12% smaller files than bitshuffle/LZ4 at about the same end-to-end speed, sitting between LZ4 and full ZSTD; for maximum ratio use BSHUF_ZSTD. Robust on any input: tests round-trip pure zeros, Poisson(10), Mersenne-Twister noise (checked against the size bound), an extreme-sparsity mask, and a real lyso image through stock ZSTD_decompress. API: exposed as "bszstd_rlehuf"; regenerate the Python/TS clients (update_version.sh) to surface the new value there. Co-Authored-By: Claude Opus 4.8 --- broker/OpenAPIConvert.cpp | 2 + broker/jfjoch_api.yaml | 1 + common/DatasetSettings.cpp | 1 + compression/CMakeLists.txt | 2 + compression/CompressionAlgorithmEnum.h | 2 +- compression/JFJochCompressor.cpp | 4 + compression/JFJochCompressor.h | 2 + compression/JFJochDecompress.h | 7 +- compression/JFJochZstdHuffCompressor.cpp | 164 ++++++++++++++++++++++ compression/JFJochZstdHuffCompressor.h | 39 +++++ compression/MaxCompressedSize.cpp | 1 + frame_serialize/CBORStream2Serializer.cpp | 1 + tests/ZSTDCompressorTest.cpp | 63 ++++++++- writer/HDF5Objects.cpp | 1 + 14 files changed, 286 insertions(+), 4 deletions(-) create mode 100644 compression/JFJochZstdHuffCompressor.cpp create mode 100644 compression/JFJochZstdHuffCompressor.h diff --git a/broker/OpenAPIConvert.cpp b/broker/OpenAPIConvert.cpp index d137c2e1..8da187f9 100644 --- a/broker/OpenAPIConvert.cpp +++ b/broker/OpenAPIConvert.cpp @@ -625,6 +625,8 @@ DatasetSettings Convert(const org::openapitools::server::model::Dataset_settings ret.Compression(CompressionAlgorithm::BSHUF_ZSTD); else if (compr == "bszstd_rle") ret.Compression(CompressionAlgorithm::BSHUF_ZSTD_RLE); + else if (compr == "bszstd_rlehuf") + ret.Compression(CompressionAlgorithm::BSHUF_ZSTD_RLE_HUFF); else if (compr == "none") ret.Compression(CompressionAlgorithm::NO_COMPRESSION); else diff --git a/broker/jfjoch_api.yaml b/broker/jfjoch_api.yaml index 95845dfc..8e78be09 100644 --- a/broker/jfjoch_api.yaml +++ b/broker/jfjoch_api.yaml @@ -428,6 +428,7 @@ components: "bslz4", "bszstd", "bszstd_rle", + "bszstd_rlehuf", "none" ] default: "bslz4" diff --git a/common/DatasetSettings.cpp b/common/DatasetSettings.cpp index 2bcb144b..db185921 100644 --- a/common/DatasetSettings.cpp +++ b/common/DatasetSettings.cpp @@ -93,6 +93,7 @@ DatasetSettings &DatasetSettings::Compression(CompressionAlgorithm input) { case CompressionAlgorithm::BSHUF_LZ4: case CompressionAlgorithm::BSHUF_ZSTD: case CompressionAlgorithm::BSHUF_ZSTD_RLE: + case CompressionAlgorithm::BSHUF_ZSTD_RLE_HUFF: compression = input; break; default: diff --git a/compression/CMakeLists.txt b/compression/CMakeLists.txt index aacc889e..73a875d3 100644 --- a/compression/CMakeLists.txt +++ b/compression/CMakeLists.txt @@ -8,6 +8,8 @@ ADD_LIBRARY(Compression STATIC bitshuffle_hperf/bitshuffle.h JFJochZstdCompressor.cpp JFJochZstdCompressor.h + JFJochZstdHuffCompressor.cpp + JFJochZstdHuffCompressor.h JFJochCompressor.cpp JFJochCompressor.h JFJochDecompress.h diff --git a/compression/CompressionAlgorithmEnum.h b/compression/CompressionAlgorithmEnum.h index 783ddb14..e6559140 100644 --- a/compression/CompressionAlgorithmEnum.h +++ b/compression/CompressionAlgorithmEnum.h @@ -3,4 +3,4 @@ #pragma once -enum class CompressionAlgorithm {BSHUF_LZ4, BSHUF_ZSTD, BSHUF_ZSTD_RLE, NO_COMPRESSION}; +enum class CompressionAlgorithm {BSHUF_LZ4, BSHUF_ZSTD, BSHUF_ZSTD_RLE, BSHUF_ZSTD_RLE_HUFF, NO_COMPRESSION}; diff --git a/compression/JFJochCompressor.cpp b/compression/JFJochCompressor.cpp index 3470e728..b837fa5b 100644 --- a/compression/JFJochCompressor.cpp +++ b/compression/JFJochCompressor.cpp @@ -31,6 +31,7 @@ static size_t MaxCompressedBlockSize(CompressionAlgorithm algorithm, size_t src_ return LZ4_compressBound(src_size) + 4; case CompressionAlgorithm::BSHUF_ZSTD: case CompressionAlgorithm::BSHUF_ZSTD_RLE: + case CompressionAlgorithm::BSHUF_ZSTD_RLE_HUFF: return ZSTD_compressBound(src_size) + 4; default: return src_size + 4; @@ -69,6 +70,9 @@ size_t JFJochBitShuffleCompressor::CompressBlock(char *dest, const char *source, throw JFJochException(JFJochExceptionCategory::ZSTDCompressionError, e.what()); } break; + case CompressionAlgorithm::BSHUF_ZSTD_RLE_HUFF: + compressed_size = huff_compressor.Compress(((uint8_t *) dest) + 4, (const uint64_t *) src_ptr, src_size); + break; default: throw JFJochException(JFJochExceptionCategory::Compression, "Algorithm not supported"); } diff --git a/compression/JFJochCompressor.h b/compression/JFJochCompressor.h index 644900d2..c4591c71 100644 --- a/compression/JFJochCompressor.h +++ b/compression/JFJochCompressor.h @@ -12,9 +12,11 @@ #include "MaxCompressedSize.h" #include "JFJochZstdCompressor.h" +#include "JFJochZstdHuffCompressor.h" class JFJochBitShuffleCompressor { JFJochZstdCompressor zstd_compressor; + JFJochZstdHuffCompressor huff_compressor; CompressionAlgorithm algorithm; std::vector tmp_space; std::vector scratch; diff --git a/compression/JFJochDecompress.h b/compression/JFJochDecompress.h index a2479d04..0e54a3d0 100644 --- a/compression/JFJochDecompress.h +++ b/compression/JFJochDecompress.h @@ -30,7 +30,8 @@ inline size_t JFJochDecompressHperfPtr(uint8_t *output, size_t block_size) { if ((algorithm != CompressionAlgorithm::BSHUF_LZ4) && (algorithm != CompressionAlgorithm::BSHUF_ZSTD) && - (algorithm != CompressionAlgorithm::BSHUF_ZSTD_RLE)) + (algorithm != CompressionAlgorithm::BSHUF_ZSTD_RLE) && + (algorithm != CompressionAlgorithm::BSHUF_ZSTD_RLE_HUFF)) throw JFJochException(JFJochExceptionCategory::Compression, "Algorithm not supported by hperf decompressor"); if ((block_size % BSHUF_BLOCKED_MULT) != 0) @@ -65,7 +66,8 @@ inline size_t JFJochDecompressHperfPtr(uint8_t *output, break; } case CompressionAlgorithm::BSHUF_ZSTD: - case CompressionAlgorithm::BSHUF_ZSTD_RLE: { + case CompressionAlgorithm::BSHUF_ZSTD_RLE: + case CompressionAlgorithm::BSHUF_ZSTD_RLE_HUFF: { const size_t ret = ZSTD_decompress(decompressed_block.data(), expected_size, src_ptr, @@ -138,6 +140,7 @@ inline void JFJochDecompressPtr(uint8_t *output, } break; case CompressionAlgorithm::BSHUF_ZSTD_RLE: + case CompressionAlgorithm::BSHUF_ZSTD_RLE_HUFF: case CompressionAlgorithm::BSHUF_ZSTD: if (use_hperf) { if (JFJochDecompressHperfPtr(output, algorithm, source + 12, source_size - 12, diff --git a/compression/JFJochZstdHuffCompressor.cpp b/compression/JFJochZstdHuffCompressor.cpp new file mode 100644 index 00000000..25c449d5 --- /dev/null +++ b/compression/JFJochZstdHuffCompressor.cpp @@ -0,0 +1,164 @@ +// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +// This file hand-assembles a standard Zstandard frame so it can be decoded by stock +// ZSTD_decompress. Every structure below follows the Zstandard compression format spec: +// https://github.com/facebook/zstd/blob/dev/doc/zstd_compression_format.md +// Comments cite that document's section names so the framing can be checked clause by clause. +// Only the framing is hand-written; the Huffman tree + bitstream inside each Compressed_Block's +// Literals_Section are produced by zstd's own HUF_compress*X_repeat (so those bytes are trusted). + +#include "JFJochZstdHuffCompressor.h" + +#include +#include +#include + +extern "C" { +#include +size_t HUF_compress1X_repeat(void *dst, size_t dstSize, const void *src, size_t srcSize, + unsigned maxSymbolValue, unsigned tableLog, void *workSpace, size_t wkspSize, + HUF_CElt *hufTable, HUF_repeat *repeat, int flags); +} + +namespace { +// --- "Zstandard frames" / "Frame_Header" --- +constexpr uint32_t MAGIC_NUMBER = 0xFD2FB528u; +// Frame_Header_Descriptor: Frame_Content_Size_flag = 2 (4-byte size), Single_Segment_flag = 1 +// (=> no Window_Descriptor; Window_Size = Frame_Content_Size), no checksum, no dictionary. +constexpr uint8_t FRAME_DESCRIPTOR = 0xA0; + +// --- "Blocks": Block_Type field of Block_Header --- +constexpr uint32_t BLOCK_RAW = 0; +constexpr uint32_t BLOCK_RLE = 1; +constexpr uint32_t BLOCK_COMPRESSED = 2; + +// --- "Literals_Section_Header": Literals_Block_Type field --- +constexpr uint8_t LITERALS_COMPRESSED = 2; // carries a Huffman_Tree_Description +constexpr uint8_t LITERALS_TREELESS = 3; // reuses the tree from a previous Compressed block + +constexpr size_t CHUNK = 65536; // cap each block's regenerated size well under Block_Maximum_Size (128 KiB) +constexpr size_t RUN_MIN = 64; // runs shorter than this are absorbed into the literal stream +constexpr unsigned LIT_LOG = 11; // Huffman table-log limit for literals, per the format +} + +JFJochZstdHuffCompressor::JFJochZstdHuffCompressor() + : ctable(HUF_CTABLE_SIZE_ST(255)), entwksp(HUF_WORKSPACE_SIZE_U64) {} + +// append the low nbytes of v, little-endian (Zstd integer fields are little-endian) +void JFJochZstdHuffCompressor::put_le(uint64_t v, int nbytes) { + for (int i = 0; i < nbytes; i++) out.push_back((uint8_t)(v >> (8 * i))); +} + +// Block_Header (3 bytes, little-endian): [ Last_Block:1 | Block_Type:2 | Block_Size:21 ]. +// Last_Block (bit 0) is left 0 here and OR-ed onto the final block afterwards. Returns the +// offset of the header so that final block can be marked. +size_t JFJochZstdHuffCompressor::blk_hdr(uint32_t type, uint32_t size) { + size_t off = out.size(); + put_le((uint64_t)(type << 1) | ((uint64_t)size << 3), 3); + return off; +} + +// RLE_Block: Block_Size is the number of repeats (the regenerated size); Block_Content is the +// single repeated byte. Split into <=CHUNK pieces so Block_Size stays under Block_Maximum_Size. +void JFJochZstdHuffCompressor::emit_run(uint8_t value, size_t nbytes, size_t &last_off) { + for (size_t off = 0; off < nbytes; off += CHUNK) { + size_t c = std::min(CHUNK, nbytes - off); + last_off = blk_hdr(BLOCK_RLE, (uint32_t)c); + out.push_back(value); + } +} + +// Emit one Compressed_Block whose Literals_Section holds n Huffman-coded literal bytes and whose +// Sequences_Section is empty (Number_of_Sequences = 0). Falls back to a Raw_Block if Huffman does +// not help. The Huffman payload (tree + streams, or streams only when reusing) is produced by zstd. +void JFJochZstdHuffCompressor::emit_lit_chunk(const uint8_t *lits, size_t n, size_t &last_off) { + HUF_CElt *ct = reinterpret_cast(ctable.data()); + HUF_repeat rep = (HUF_repeat)repeat_state; + // Size_Format also selects the stream count: format 0 = 1 stream, formats 1/2/3 = 4 streams. + // Use a single stream only for small inputs (matches zstd's own ZSTD_compressLiterals heuristic). + const bool single = n < 256; + const size_t wkspSize = entwksp.size() * sizeof(uint64_t); + size_t hs = single + ? HUF_compress1X_repeat(hufbuf.data(), hufbuf.size(), lits, n, 255, LIT_LOG, entwksp.data(), wkspSize, ct, &rep, 0) + : HUF_compress4X_repeat(hufbuf.data(), hufbuf.size(), lits, n, 255, LIT_LOG, entwksp.data(), wkspSize, ct, &rep, 0); + + // Literals_Section_Header, read as a little-endian value, is laid out as + // [ Literals_Block_Type:2 | Size_Format:2 | Regenerated_Size | Compressed_Size ]. + // Size_Format (chosen from Regenerated_Size n) sets the header length and field widths: + // n < 1 KiB -> 3-byte header, 10-bit fields (format 0 single-stream, else 1) + // n < 16 KiB -> 4-byte header, 14-bit fields (format 2) + // else -> 5-byte header, 18-bit fields (format 3) + // Regenerated_Size starts at bit 4; Compressed_Size starts right after it, at bit (4 + width). + int size_format, lh_bytes, regen_width; + if (n < 1024) { size_format = single ? 0 : 1; lh_bytes = 3; regen_width = 10; } + else if (n < 16384) { size_format = 2; lh_bytes = 4; regen_width = 14; } + else { size_format = 3; lh_bytes = 5; regen_width = 18; } + + // hs == 1 is HUF's "single-symbol alphabet" signal (not a usable Huffman section); only keep + // the Huffman block when it is well-formed and actually smaller than a Raw_Block. + if (hs > 1 && !HUF_isError(hs) && (size_t)(lh_bytes + 1) + hs < n) { + repeat_state = HUF_repeat_check; // the table now in ct is valid -> reuse it next time + // HUF_compress*_repeat leaves *repeat != none only when it reused the previous table. + uint8_t lit_type = (rep != HUF_repeat_none) ? LITERALS_TREELESS : LITERALS_COMPRESSED; + uint64_t hdr = (uint64_t)lit_type | ((uint64_t)size_format << 2) + | ((uint64_t)n << 4) | ((uint64_t)hs << (4 + regen_width)); + last_off = blk_hdr(BLOCK_COMPRESSED, (uint32_t)(lh_bytes + hs + 1)); + put_le(hdr, lh_bytes); // Literals_Section_Header + out.insert(out.end(), hufbuf.data(), hufbuf.data() + hs); // Huffman tree + stream(s) + out.push_back(0x00); // Sequences_Section: Number_of_Sequences = 0 + } else { + repeat_state = HUF_repeat_none; // discard any table HUF rebuilt but we did not emit + last_off = blk_hdr(BLOCK_RAW, (uint32_t)n); + out.insert(out.end(), lits, lits + n); + } +} + +// src = one bitshuffled block (src_size bytes, a multiple of 8). Emits a complete Zstd frame. +size_t JFJochZstdHuffCompressor::Compress(uint8_t *dst, const uint64_t *src, size_t src_size) { + const size_t W = src_size / 8; + out.clear(); literals.clear(); segs.clear(); + repeat_state = HUF_repeat_none; + if (hufbuf.size() < src_size + 1024) hufbuf.resize(src_size + 1024); + + // Pass 1: classify 8-byte words into runs (0x00 / 0xFF) and literals. Runs >= RUN_MIN become + // RLE_Blocks; shorter runs are folded into the literal stream (where 0x00 costs ~1 Huffman bit). + size_t lit_start = 0; + auto close_lit = [&]() { if (literals.size() > lit_start) { + segs.push_back({2, literals.size() - lit_start, lit_start}); lit_start = literals.size(); } }; + for (size_t i = 0; i < W; ) { + if (src[i] == 0 || src[i] == UINT64_MAX) { + uint64_t val = src[i]; size_t j = i; while (j < W && src[j] == val) j++; + size_t bytes = (j - i) * 8; uint8_t v = val ? 0xFF : 0x00; + if (bytes >= RUN_MIN) { close_lit(); segs.push_back({(uint8_t)(val ? 1 : 0), bytes, 0}); } + else literals.insert(literals.end(), bytes, v); + i = j; + } else { + const uint8_t *wb = reinterpret_cast(&src[i]); + literals.insert(literals.end(), wb, wb + 8); + i++; + } + } + close_lit(); + + // Frame_Header: Magic_Number + Frame_Header_Descriptor + Frame_Content_Size (4 bytes). + put_le(MAGIC_NUMBER, 4); + out.push_back(FRAME_DESCRIPTOR); + put_le(src_size, 4); + + // Pass 2: emit the blocks in order. + size_t last_off = 0; bool emitted = false; + for (const Seg &s : segs) { + emitted = true; + if (s.type == 2) + for (size_t off = 0; off < s.bytes; off += CHUNK) + emit_lit_chunk(literals.data() + s.lit_off + off, std::min(CHUNK, s.bytes - off), last_off); + else + emit_run(s.type == 1 ? 0xFF : 0x00, s.bytes, last_off); + } + if (!emitted) last_off = blk_hdr(BLOCK_RAW, 0); // empty input -> a single empty block + out[last_off] |= 1; // set Last_Block on the final block + + memcpy(dst, out.data(), out.size()); + return out.size(); +} diff --git a/compression/JFJochZstdHuffCompressor.h b/compression/JFJochZstdHuffCompressor.h new file mode 100644 index 00000000..4d1ac84e --- /dev/null +++ b/compression/JFJochZstdHuffCompressor.h @@ -0,0 +1,39 @@ +// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include + +// Produces a STANDARD Zstandard frame from bitshuffled data, decodable by stock ZSTD_decompress: +// - zero / 0xFF runs -> RLE_Blocks (cheap, like JFJochZstdCompressor) +// - literal regions -> Compressed_Blocks with Huffman literals (no sequences); short +// runs are absorbed into the literal stream +// - incompressible literals -> Raw_Blocks (bounded worst case) +// Faster than full ZSTD (no match search) and better ratio than the plain RLE compressor (it +// entropy-codes the literals). The Huffman table is built/reused per block from that block's own +// literals via zstd's HUF_compress*X_repeat, so it is robust on any input (random, Poisson, masks, +// zeros) with no trained tables. +class JFJochZstdHuffCompressor { + std::vector out; // assembled frame + std::vector literals; // literal bytes in stream order (incl. absorbed short runs) + std::vector hufbuf; // scratch for one Huffman-coded literal chunk + std::vector ctable; // HUF_CElt[] (size_t-aligned) + std::vector entwksp; // HUF compression workspace + unsigned repeat_state = 0; // HUF_repeat across literal blocks within the current frame + + struct Seg { uint8_t type; size_t bytes; size_t lit_off; }; // type 0=run0, 1=runFF, 2=literals + std::vector segs; + + void put_le(uint64_t v, int nbytes); + size_t blk_hdr(uint32_t type, uint32_t size); + void emit_run(uint8_t value, size_t nbytes, size_t &last_off); + void emit_lit_chunk(const uint8_t *lits, size_t n, size_t &last_off); +public: + JFJochZstdHuffCompressor(); + // src = bitshuffled block (src_size bytes, a multiple of 8). Writes one zstd frame to dst and + // returns its size. dst must hold at least ZSTD_compressBound(src_size) + 12 bytes. + size_t Compress(uint8_t *dst, const uint64_t *src, size_t src_size); +}; diff --git a/compression/MaxCompressedSize.cpp b/compression/MaxCompressedSize.cpp index 613a3480..1a28e0c8 100644 --- a/compression/MaxCompressedSize.cpp +++ b/compression/MaxCompressedSize.cpp @@ -12,6 +12,7 @@ int64_t MaxCompressedSize(CompressionAlgorithm algorithm, int64_t pixels_number, return bshuf_compress_lz4_bound(pixels_number, pixel_depth, JFJochBitShuffleCompressor::BlockSize(algorithm, pixel_depth)) + 12; case CompressionAlgorithm::BSHUF_ZSTD: case CompressionAlgorithm::BSHUF_ZSTD_RLE: + case CompressionAlgorithm::BSHUF_ZSTD_RLE_HUFF: return bshuf_compress_zstd_bound(pixels_number, pixel_depth, JFJochBitShuffleCompressor::BlockSize(algorithm, pixel_depth)) + 12; default: return pixels_number * pixel_depth; diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index 0cf19abd..73cea1c7 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -68,6 +68,7 @@ void CBOR_ENC_COMPRESSED(CborEncoder &encoder, break; case CompressionAlgorithm::BSHUF_ZSTD: case CompressionAlgorithm::BSHUF_ZSTD_RLE: + case CompressionAlgorithm::BSHUF_ZSTD_RLE_HUFF: cborErr(cbor_encode_text_stringz(&arrayEncoder, "bszstd")); break; default: diff --git a/tests/ZSTDCompressorTest.cpp b/tests/ZSTDCompressorTest.cpp index ba2fe8a4..29d8fb25 100644 --- a/tests/ZSTDCompressorTest.cpp +++ b/tests/ZSTDCompressorTest.cpp @@ -10,6 +10,7 @@ #include "../compression/JFJochCompressor.h" #include "../compression/JFJochDecompress.h" #include "../common/DiffractionExperiment.h" +#include "../writer/HDF5Objects.h" TEST_CASE("JFjochZstdCompressor_Raw_Block","[ZSTD]") { uint8_t src[128*8]; @@ -315,4 +316,64 @@ TEST_CASE("Bitshuffle_ZSTD","[ZSTD]") { REQUIRE(bshuf_decompress_zstd(compressed.data(), decompressed.data(), RAW_MODULE_SIZE, 4, 0) == out_size); REQUIRE(memcmp(image.data(), decompressed.data(), RAW_MODULE_SIZE*sizeof(uint32_t)) == 0); -} \ No newline at end of file +} +// ---- BSHUF_ZSTD_RLE_HUFF (RLE runs + adaptive Huffman literals, standard zstd frame) ---- +// Round-trips through JFJochDecompress (which calls stock ZSTD_decompress per block) and checks +// the output never exceeds the advertised worst-case size. Run on good and adversarial data to +// confirm the algorithm always produces valid output, even when the ratio is poor. +template +static void RequireHuffRoundTrip(const std::vector &image) { + JFJochBitShuffleCompressor compressor(CompressionAlgorithm::BSHUF_ZSTD_RLE_HUFF); + auto compressed = compressor.Compress(image); + + REQUIRE(compressed.size() <= (size_t) MaxCompressedSize(CompressionAlgorithm::BSHUF_ZSTD_RLE_HUFF, + image.size(), sizeof(T))); + std::vector output; + REQUIRE_NOTHROW(JFJochDecompress(output, CompressionAlgorithm::BSHUF_ZSTD_RLE_HUFF, + compressed, image.size())); + REQUIRE(output.size() == image.size()); + REQUIRE(memcmp(image.data(), output.data(), image.size() * sizeof(T)) == 0); +} + +TEST_CASE("ZstdHuff_PureZeros", "[ZSTD]") { + std::vector image(200000, 0); // all runs, no literals + RequireHuffRoundTrip(image); +} + +TEST_CASE("ZstdHuff_Poisson10", "[ZSTD]") { + std::vector image(200000); + std::mt19937 gen(12345); + std::poisson_distribution dist(10); + for (auto &v : image) v = (uint16_t) dist(gen); + RequireHuffRoundTrip(image); +} + +TEST_CASE("ZstdHuff_RandomIncompressible", "[ZSTD]") { + std::vector image(200000); // Mersenne-Twister noise: ~zero compression + std::mt19937 gen(98765); + for (auto &v : image) v = (uint16_t) gen(); + RequireHuffRoundTrip(image); // must still be valid and within bound +} + +TEST_CASE("ZstdHuff_MaskLike", "[ZSTD]") { + std::vector image(200000, 0); // extreme sparsity like a pixel_mask + std::mt19937 gen(555); + for (auto &v : image) { + uint32_t r = gen() % 100; + if (r < 5) v = 1; + else if (r == 5) v = 0x80000000u; // module-gap style flag + } + RequireHuffRoundTrip(image); +} + +TEST_CASE("ZstdHuff_LysoImage", "[ZSTD]") { + RegisterHDF5Filter(); // bitshuffle filter, needed to read the compressed benchmark dataset + HDF5ReadOnlyFile data("../../tests/test_data/compression_benchmark.h5"); + HDF5DataSet dataset(data, "/entry/data/data"); + HDF5DataSpace file_space(dataset); + auto dims = file_space.GetDimensions(); + std::vector image(dims[1] * dims[2]); + std::vector start = {0, 0, 0}, size = {1, dims[1], dims[2]}; + dataset.ReadVector(image, start, size); + RequireHuffRoundTrip(image); +} diff --git a/writer/HDF5Objects.cpp b/writer/HDF5Objects.cpp index dc8329ec..66aa2fe2 100644 --- a/writer/HDF5Objects.cpp +++ b/writer/HDF5Objects.cpp @@ -972,6 +972,7 @@ void HDF5Dcpl::SetCompression(CompressionAlgorithm c, size_t block_size) { break; case CompressionAlgorithm::BSHUF_ZSTD: case CompressionAlgorithm::BSHUF_ZSTD_RLE: + case CompressionAlgorithm::BSHUF_ZSTD_RLE_HUFF: #ifdef USE_ZSTD params[0] = block_size; params[1] = BSHUF_H5_COMPRESS_ZSTD; -- 2.54.0 From ac75177c5a1dd71105580c8609c63a1d2d4a8853 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 27 Jun 2026 14:47:16 +0200 Subject: [PATCH 12/66] Test: prove Huffman-literals Compressed_Block with Number_of_Sequences=0 is valid Builds a single Compressed_Block (Huffman-coded Literals_Section, empty Sequences_Section) and checks: the block type is Compressed, its trailing Number_of_Sequences byte is 0, and stock ZSTD_decompress reconstructs the literals exactly. This is the format guarantee from zstd_compression_format.md ("if Number_of_Sequences == 0 ... Block's decompressed content is defined solely by the Literals Section content"), locked into the test suite. Co-Authored-By: Claude Opus 4.8 --- tests/ZSTDCompressorTest.cpp | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/tests/ZSTDCompressorTest.cpp b/tests/ZSTDCompressorTest.cpp index 29d8fb25..3c840a04 100644 --- a/tests/ZSTDCompressorTest.cpp +++ b/tests/ZSTDCompressorTest.cpp @@ -377,3 +377,33 @@ TEST_CASE("ZstdHuff_LysoImage", "[ZSTD]") { dataset.ReadVector(image, start, size); RequireHuffRoundTrip(image); } + +// On-the-record proof that a Compressed_Block with Huffman literals and Number_of_Sequences = 0 is +// a valid Zstandard block that outputs its literals. Per zstd_compression_format.md: "if +// (Number_of_Sequences == 0) ... Block's decompressed content is defined solely by the Literals +// Section content." We build exactly such a block and confirm stock ZSTD_decompress returns it. +TEST_CASE("ZstdHuff_LiteralsOnly_NbSeqZero", "[ZSTD]") { + std::vector lits(4096); + std::mt19937 gen(7); + for (auto &b : lits) b = (uint8_t)(gen() % 6 + 1); // 6-symbol alphabet: Huffman-compressible, no 0x00/0xFF runs + + JFJochZstdHuffCompressor huff; + std::vector frame(lits.size() * 2 + 4096); + size_t fsz = huff.Compress(frame.data(), (const uint64_t *) lits.data(), lits.size()); + + // Frame = Magic(4) + Frame_Header_Descriptor(1) + Frame_Content_Size(4) = 9 bytes, then one Block. + // Block_Header (3 bytes, little-endian): [ Last_Block:1 | Block_Type:2 | Block_Size:21 ]. + uint32_t bh = frame[9] | (frame[10] << 8) | (frame[11] << 16); + bool last_block = bh & 1; + uint32_t block_type = (bh >> 1) & 3; + uint32_t block_size = bh >> 3; + REQUIRE(last_block); // a single block for this all-literal input + REQUIRE(block_type == 2); // Compressed_Block (Huffman literals path, not Raw/RLE) + REQUIRE(frame[12 + block_size - 1] == 0x00); // its Sequences_Section is Number_of_Sequences = 0 + + std::vector out(lits.size()); + size_t r = ZSTD_decompress(out.data(), out.size(), frame.data(), fsz); + REQUIRE_FALSE(ZSTD_isError(r)); + REQUIRE(r == lits.size()); + REQUIRE(memcmp(lits.data(), out.data(), lits.size()) == 0); +} -- 2.54.0 From 809f6f92a333b642464df9b8f37b7e1be4985723 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 27 Jun 2026 14:52:12 +0200 Subject: [PATCH 13/66] Regenerate OpenAPI clients for bszstd_rlehuf compression value update_version.sh regenerated the TypeScript client (types.gen.ts, zod.gen.ts) and Redoc docs to include the new "bszstd_rlehuf" compression enum value added to jfjoch_api.yaml; the package-lock version field follows VERSION. The C++ server model treats compression as a plain string, so it needed no change. Co-Authored-By: Claude Opus 4.8 --- broker/redoc-static.html | 4 ++-- frontend/package-lock.json | 4 ++-- frontend/src/client/types.gen.ts | 2 +- frontend/src/client/zod.gen.ts | 1 + 4 files changed, 6 insertions(+), 5 deletions(-) diff --git a/broker/redoc-static.html b/broker/redoc-static.html index 76176298..d7262869 100644 --- a/broker/redoc-static.html +++ b/broker/redoc-static.html @@ -457,7 +457,7 @@ Incident particle (photon, electron) energy in keV

space_group_number
integer <int64> [ 1 .. 194 ]

Number of space group for the crystal. Currently used solely as metadata, not relevant for image processing done in Jungfraujoch.

sample_name
string
Default: ""

/entry/sample/name in NXmx Sample name

-
compression
string
Default: "bslz4"
Enum: "bslz4" "bszstd" "bszstd_rle" "none"

Compression type for the images transferred over ZeroMQ and saved to HDF5 file.

+
compression
string
Default: "bslz4"
Enum: "bslz4" "bszstd" "bszstd_rle" "bszstd_rlehuf" "none"

Compression type for the images transferred over ZeroMQ and saved to HDF5 file.

total_flux
number <float>

/entry/beam/total_flux in NXmx Flux incident on beam plane in photons per second. In other words this is the flux integrated over area. [photons/s]

transmission
number <float> [ 0 .. 1 ]

/entry/instrument/attenuator/attenuator_transmission @@ -937,7 +937,7 @@ then image might be replaced in the buffer between calling /images and /image.cb