From db908c6ec160a7e3fedd32273822869c76b5ece0 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 25 Aug 2026 21:21:55 +0200 Subject: [PATCH] _process.h5: record the change of basis, and the space group A _process.h5 was internally inconsistent whenever the space group re-seated the lattice. The per-image reflections and lattices go to file as each image is processed, but the conventional setting is only chosen afterwards, so the file kept pre-reindex indices beside a post-reindex cell. Measured on an affected dataset, max|d_file - d(hkl, stored cell)| was 18.6 A. --mode scale then compared each frame's lattice against the stored cell, found none of 1800 acceptable, rejected every observation and died in the merge with "resolution calculation failed". /entry/MX/reindexMatrix now carries M with hkl_cell = M . hkl_written, and the reader applies it, so everything it hands out is in the setting of /entry/sample/unit_cell. Absent means the identity, so a file written before this reads exactly as before. On the affected dataset M comes out [[1,1,0],[0,1,1],[1,0,1]], det 2 - the primitive-to-body-centred basis its volume ratio implied - and the same measure falls from 18.6 A to 1.8e-5. Writing the reflections in the final setting instead was rejected: the per-image writer is shared with the broker, which streams and cannot buffer a run; and h,k,l, predicted_x/y and the per-image lattice are one consistent statement about one image, which retro-editing the indices would silently break. Two things turned up while fixing it. There are three re-seat sites, not one - the space-group search's own centred-lattice test re-seats too, and logged nothing - so the matrix composes over all of them. And the space group itself was never written: it was set only on the arm that searches, while a two-pass rotation run reuses pass 1's group and takes the other arm, so the canonical file carried a cell but no group and --mode scale merged in P1. --mode scale now reproduces --mode mx on the affected dataset: same space group and cell, 39329 unique reflections both, 100% of reflections common, CC 0.99985, sum|dI|/sum|I| = 0.0081. The residual is three reflections in 11.09 M crossing an ice-band edge, because mx carries the integrator's d and scale recomputes it from the cell. An unaffected dataset is byte-identical in .mtz, .hkl and .cif. An older affected file still cannot merge - M is not recoverable from it - but now says so in 51 s, naming both cells and the -S/-C override to use, instead of failing inside the merge. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P --- common/JFJochMessages.h | 8 +++ docs/CHANGELOG.md | 1 + docs/HDF5.md | 14 +++++ docs/RUGNUX.md | 4 +- reader/HDF5MetadataSource.cpp | 53 +++++++++++++--- reader/JFJochReaderDataset.h | 10 +++ rugnux/Rugnux.cpp | 23 +++++++ rugnux/rugnux_cli.cpp | 44 ++++++++++++++ tests/JFJochReaderTest.cpp | 111 ++++++++++++++++++++++++++++++++++ writer/HDF5NXmx.cpp | 8 +++ 10 files changed, 266 insertions(+), 10 deletions(-) diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index 081de0f1..59fb283e 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -4,6 +4,7 @@ #pragma once #include +#include #include #include #include @@ -368,6 +369,13 @@ struct EndMessage { // 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; + // Change of basis (3x3 integers, row major) from the setting the PER-IMAGE reflections and + // lattices were written in to the setting of unit_cell / rotation_lattice above: + // hkl_cell = M . hkl_written, and the same M takes each per-image lattice across. The per-image + // data go to file as each image is processed, but the space group - and with it the conventional + // setting the merge works in - is only chosen afterwards, so the two settings are not always the + // same one. Absent means they are, i.e. the identity. + std::optional> reindex_matrix; // Geometry refined by the offline analysis. When present these overwrite the master-file geometry // at finalization (rugnux writes the refined values; the broker leaves them empty and keeps the diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 20fc56be..0f9e144c 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -8,6 +8,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * rugnux: where the space-group search reports alternatives that are not on the same lattice as the group it chose (`C2 or P21 or P2`), it now names the ones the reported unit cell does not describe and says how their cell differs from it. * rugnux: a lattice centring the data could not test - the crystal was integrated on the primitive sub-cell, so the reflections the centring extinguishes were never measured - is marked `UNTESTED` in the space-group candidate table and, where it is adopted, is warned about as coming from the lattice metric rather than from the intensities. * rugnux: the space-group search prints the twin-law disagreement H for every operator it tested and the H ratio of the point group it adopted, on every run, instead of only when that ratio refuses a promotion. +* rugnux: `--mode scale` works on a `_process.h5` whose space group came from re-seating the lattice; a file written before this stops with a message naming the two cells instead of failing inside the merge. * rugnux reports diffraction anisotropy: the anisotropic deltaB and the diffraction limit along each * rugnux: `--mode scale` reports the Wilson B-factor estimate instead of `WILSON_B= nan`. * rugnux: `--export-unmerged` writes the integrated observations as `_unmerged.mtz`, an unmerged MTZ readable by aimless, pointless, careless and `iotbx.merging_statistics`, in `--mode mx` and `--mode scale` alike; each rotation reflection's partials are summed into one full, and `--export-unmerged-partials` writes one row per image instead. Lattice-centring absences are not written; screw and glide absences are. diff --git a/docs/HDF5.md b/docs/HDF5.md index 5209f75d..e9966edd 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -327,6 +327,20 @@ variants. | `bkgEstimateMean` | photons | mean background over the run | | `iceRingScoreMean` | ratio | mean `iceRingScore` over the run — the single "how icy was this dataset" number (1 = no ice) | | `indexedLatticeCount` | | per-image lattice count summary (master). *Note: data files use `indexingLatticeCount`; readers accept either.* | +| `reindexMatrix` | | change of basis from the setting the per-image data are in to the setting of `/entry/sample/unit_cell` (`[9]`, `int32`, flattened 3×3, row major) — see below | + +**Reindex matrix.** The per-image `h`, `k`, `l` and `latticeIndexed` are written as each image is +processed, in the setting that image was *indexed* in. The space group is only settled afterwards, by +the merge, and settling it can re-seat the lattice into the group's conventional setting — so +`/entry/sample/unit_cell`, `/entry/sample/space_group_number` and `rotationLatticeIndexed` can be in a +different setting from the per-image data beside them. `reindexMatrix` **M** is the integral change of +basis between the two: `hkl_cell = M · hkl_written`, and the same **M** takes each per-image lattice +across (`latticeIndexedExtra` is not re-seated and stays as indexed). It is **absent** when the two +settings are the same one, which is the identity — as it is on every file written before rugnux +recorded it. The Jungfraujoch reader applies it, so everything it hands out is already in the cell's +setting; a third-party reader that ignores it will index the reflections in the wrong frame whenever +the dataset is present. Written by the offline `rugnux` path only — the broker never re-seats a +lattice — and not carried on the CBOR stream, in the same way as the other offline-only fields. **Sweep quality.** `sweepQuality` `[n_images]` (`uint8`) says why the stretch of the sweep this image belongs to was flagged as delivering much less than the rest of the run: **0** means it was diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 86fa063c..bd1135ce 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -237,7 +237,9 @@ The `scale` mode re-scales and merges the *already-integrated* reflections store `_process.h5` file, without re-running spot finding or integration. Use it to re-merge quickly with a different space group, resolution limit, anomalous setting or reference MTZ. It reuses the same `-o/-N/-s/-e/-S/-A/-B/-z/--scaling-*` options as the full run, and (unlike the full pipeline) does -not run a space-group search, so pass `-S` for the correct symmetry. +not run a space-group search: it merges in the space group and unit cell the file records, and `-S` +/ `-C` override them. A `_process.h5` written before the group was stored carries none, and merges +in P1 unless `-S` says otherwise. ## Detector calibration from powder rings (`rugnux --mode calibration`) diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 93755c1a..1ea3a69e 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -132,9 +132,22 @@ std::string dataset_name(const std::string &path) { return file; } +// Per-image reflections and lattices are written in the setting the images were INDEXED in; the +// unit cell, the run lattice and the space group beside them are in the setting the merge settled +// on, which the space-group search can re-seat to. /entry/MX/reindexMatrix is the integral change of +// basis between the two, so applying it here is what makes the file read as one consistent dataset. +// No matrix means the two settings are the same one. +CrystalLattice ApplyReindex(const CrystalLattice &latt, const std::optional> &m) { + if (!m) + return latt; + const auto &v = *m; + return latt.Multiply(gemmi::Mat33(v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8])); +} + bool ReadReflectionsFromGroup(HDF5Object &file, const std::string &image_group_name, - std::vector &reflections) { + std::vector &reflections, + const std::optional> &reindex) { if (!file.Exists("/entry/reflections") || !file.Exists(image_group_name)) return false; @@ -165,6 +178,15 @@ bool ReadReflectionsFromGroup(HDF5Object &file, throw JFJochException(JFJochExceptionCategory::HDF5, "Wrong size of reflections dataset"); for (size_t i = 0; i < h.size(); i++) { + int32_t hh = h.at(i), kk = k.at(i), ll = l.at(i); + if (reindex) { + const auto &m = *reindex; + const int32_t h0 = hh, k0 = kk, l0 = ll; + hh = m[0] * h0 + m[1] * k0 + m[2] * l0; + kk = m[3] * h0 + m[4] * k0 + m[5] * l0; + ll = m[6] * h0 + m[7] * k0 + m[8] * l0; + } + float lp_val = 0.0; if (lp.size() > i && lp[i] != 0.0f) lp_val = 1.0f / lp[i]; @@ -201,9 +223,9 @@ bool ReadReflectionsFromGroup(HDF5Object &file, } Reflection r{ - .h = h.at(i), - .k = k.at(i), - .l = l.at(i), + .h = hh, + .k = kk, + .l = ll, .image_number = image_number.at(i), .delta_phi_deg = delta_phi_val, .predicted_x = predicted_x.at(i), @@ -687,6 +709,13 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen .gamma = tmp[5] }); dataset->experiment.SpaceGroupNumber(master_file->GetOptInt("/entry/sample/space_group_number")); + // The setting the cell and space group just read are in, relative to the setting the per-image + // reflections and lattices were written in. Absent on every file written before the offline + // analysis started recording it, and on every run that never re-seated its lattice - both mean + // the identity, and both are read as such. + if (const auto m = master_file->ReadOptVector("/entry/MX/reindexMatrix"); m.size() == 9) + dataset->reindex_matrix = std::array{m[0], m[1], m[2], m[3], m[4], + m[5], m[6], m[7], m[8]}; dataset->experiment.SampleName(master_file->GetString("/entry/sample/name")); @@ -1070,7 +1099,7 @@ void HDF5MetadataSource::FillPerImage(DataMessage &message, int64_t requested_im ); if (tmp.size() == 9) - message.indexing_lattice = CrystalLattice(tmp); + message.indexing_lattice = ApplyReindex(CrystalLattice(tmp), dataset->reindex_matrix); std::optional lattice; if (master_file->Exists("/entry/MX/bravaisLattice")) @@ -1102,8 +1131,10 @@ void HDF5MetadataSource::FillPerImage(DataMessage &message, int64_t requested_im const std::string master_reflection_group_name = fmt::format("/entry/reflections/image_{:06d}", image_number); const std::string source_reflection_group_name = fmt::format("/entry/reflections/image_{:06d}", image_id); - if (!ReadReflectionsFromGroup(*master_file, master_reflection_group_name, message.reflections)) - ReadReflectionsFromGroup(*source_file, source_reflection_group_name, message.reflections); + if (!ReadReflectionsFromGroup(*master_file, master_reflection_group_name, message.reflections, + dataset->reindex_matrix)) + ReadReflectionsFromGroup(*source_file, source_reflection_group_name, message.reflections, + dataset->reindex_matrix); if (!message.reflections.empty()) { CalcISigma(message); CalcWilsonBFactor(message, !message.b_factor.has_value()); @@ -1214,6 +1245,10 @@ std::vector HDF5MetadataSource::ReadReflections(size_t start // reflections lazily from the source data files instead. const bool master_reflections_authoritative = master_file->Exists("/entry/reflections"); + // Everything below comes out in the setting of the dataset's unit cell and space group, not in the + // setting it was written in (see ApplyReindex). + const auto &reindex = dataset_->reindex_matrix; + for (size_t img = start_image; img <= end_image_val; img++) { IntegrationOutcome outcome; @@ -1235,7 +1270,7 @@ std::vector HDF5MetadataSource::ReadReflections(size_t start } // ── reflections ────────────────────────────────────────────────────── - ReadReflectionsFromGroup(*meta_file, refl_group, outcome.reflections); + ReadReflectionsFromGroup(*meta_file, refl_group, outcome.reflections, reindex); // ── per-image mosaicity ─────────────────────────────────────────────── if (meta_file->Exists("/entry/MX/mosaicity")) { @@ -1252,7 +1287,7 @@ std::vector HDF5MetadataSource::ReadReflections(size_t start auto lattice_vec = meta_file->ReadOptVector( "/entry/MX/latticeIndexed", {meta_image_id, 0}, {1, 9}); if (lattice_vec.size() == 9) - outcome.latt = CrystalLattice(lattice_vec); + outcome.latt = ApplyReindex(CrystalLattice(lattice_vec), reindex); } catch (...) { } } diff --git a/reader/JFJochReaderDataset.h b/reader/JFJochReaderDataset.h index b03c04e1..ab6d36e8 100644 --- a/reader/JFJochReaderDataset.h +++ b/reader/JFJochReaderDataset.h @@ -4,6 +4,7 @@ #pragma once #include +#include #include #include #include @@ -32,6 +33,15 @@ struct JFJochReaderDataset { std::string jfjoch_release; + // Change of basis (3x3 integers, row major) from the setting the per-image reflections and + // lattices in this file were written in to the setting of experiment's unit cell and space group: + // hkl_cell = M . hkl_written. The per-image data are written as each image is processed, but the + // space group - and the conventional setting that comes with it - is only settled after the merge, + // so the two are not always the same setting. Empty means they are. Read straight from + // /entry/MX/reindexMatrix; the reader applies it, so everything this class hands out is already in + // the cell's setting. + std::optional> reindex_matrix; + std::vector az_int_bin_to_q; std::vector az_int_bin_to_phi; size_t azimuthal_bins = 0; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 6bbde29e..77421591 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2059,6 +2059,21 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // new cell's i-th axis, and indices transform with the axes, so it is also the matrix that takes // hkl across. Rounding it to integers - and declining when the residual is large - keeps it exact // even though the refined cell is not exactly LatticeSearch's Niggli cell. + // + // The per-image reflections and lattices are already in the _process.h5 by now - they go to file + // as each image is processed, long before the group is known - so re-seating them here leaves the + // stored ones in the OLD setting beside a cell in the new one. Accumulate the change of basis + // (composed, in case more than one fires) and hand it to the writer, so the file says which setting + // its per-image data are in and a reader can bring them across. + gemmi::Mat33 stored_to_final; // identity until something below re-seats + const auto record_reindex = [&](const gemmi::Mat33 &m) { + stored_to_final = m.multiply(stored_to_final); + std::array v{}; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + v[3 * i + j] = static_cast(std::lround(stored_to_final.a[i][j])); + end_msg.reindex_matrix = v; + }; const auto reindex_into = [&](const LatticeSearchResult &cand) { const Coord cv[3] = {cand.conventional.Vec0(), cand.conventional.Vec1(), cand.conventional.Vec2()}; const Coord rs[3] = {end_msg.rotation_lattice->Astar(), end_msg.rotation_lattice->Bstar(), @@ -2087,6 +2102,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b end_msg.rotation_lattice = end_msg.rotation_lattice->Multiply(reindex); end_msg.rotation_lattice_type = LatticeMessage{ .centering = cand.centering, .niggli_class = end_msg.rotation_lattice_type->niggli_class, .crystal_system = cand.system }; + record_reindex(reindex); return true; }; @@ -2680,6 +2696,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // vectors) in the same conventional setting as the cell, reflections and space group. end_msg.rotation_lattice = commit_lattice; end_msg.rotation_lattice_type = commit_lattice_type; + record_reindex(*commit_reindex); if (rsm) { rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell, static_cast(config_.scaling_iter), @@ -2697,6 +2714,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } else { // A space group was fixed by the user; surface it so the viewer/CLI can still show it. result.space_group_number = experiment_.GetSpaceGroupNumber(); + // ...and record it, exactly as the search arm above does with the group it determined. + // This arm also carries the two-pass second pass, which REUSES the first pass's group + // rather than re-searching, so without this the _process.h5 of every two-pass rotation run + // came out with a cell but no group - and re-merging it defaulted to P1. + if (const auto sg = experiment_.GetSpaceGroupNumber(); sg.has_value()) + end_msg.space_group_number = static_cast(*sg); } // Rotation two-pass geometry pre-pass: the determined space group is now fixed on experiment_, so the diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 70bebdc3..9d4af119 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1375,6 +1375,50 @@ static int RunRugnux(int argc, char **argv) { logger.Error("Experiment unit cell not found, cannot update reflection resolution"); exit(EXIT_FAILURE); } + // The stored reflections have to be in the setting of the cell they are about to be merged in, + // or the merge's per-frame cell gate rejects every observation and the run dies far downstream + // complaining about resolution shells. The reader has already applied /entry/MX/reindexMatrix, + // so a file that carries one arrives consistent; one written before that dataset existed, by a + // run whose space-group search re-seated the lattice after the images had gone to file, does + // not - and the change of basis cannot be recovered from it, because the stored lattices carry + // each frame's own goniometer rotation. Say so here, where both cells can be named, instead of + // merging indices in a frame they do not belong to. + { + const auto &merge_cell = experiment.GetUnitCell().value(); + const float dist_tol = experiment.GetIndexingSettings().GetUnitCellDistTolerance(); + const float ang_tol = experiment.GetIndexingSettings().GetUnitCellAngleTolerance_deg(); + size_t n_indexed = 0, n_agree = 0; + std::optional file_cell; + for (const auto &o : reflections) + if (!o.reflections.empty()) { + ++n_indexed; + const auto cell = o.latt.GetUnitCell(); + if (cell.is_close(merge_cell, dist_tol, ang_tol)) + ++n_agree; + else if (!file_cell) + file_cell = cell; + } + if (n_indexed > 0 && n_agree == 0) { + logger.Error("The per-image lattices stored in this file are not in the setting of the " + "cell it is being merged in, so no observation can be used. Merging in " + "a={:.3f} b={:.3f} c={:.3f} alpha={:.2f} beta={:.2f} gamma={:.2f}, but all " + "{} indexed images are in a={:.3f} b={:.3f} c={:.3f} alpha={:.2f} " + "beta={:.2f} gamma={:.2f}{}. Merge in the setting the file is in with " + "-S 1 -C {:.4f},{:.4f},{:.4f},{:.4f},{:.4f},{:.4f}, or reprocess the " + "images with --mode mx.", + merge_cell.a, merge_cell.b, merge_cell.c, merge_cell.alpha, + merge_cell.beta, merge_cell.gamma, n_indexed, + file_cell->a, file_cell->b, file_cell->c, file_cell->alpha, + file_cell->beta, file_cell->gamma, + dataset->reindex_matrix + ? "" : " (the file carries no /entry/MX/reindexMatrix, so it was " + "written before that was recorded)", + file_cell->a, file_cell->b, file_cell->c, file_cell->alpha, + file_cell->beta, file_cell->gamma); + exit(EXIT_FAILURE); + } + } + auto refl_stats = UpdateReflectionResolution(experiment.GetUnitCell().value(), reflections); logger.Info("Read {} reflections from {} images", refl_stats.n_reflections, refl_stats.n_images); experiment.ImagesPerTrigger(refl_stats.n_images); diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index c9c554e9..11c56c63 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -2479,6 +2479,117 @@ TEST_CASE("JFJochReader_ReadReflections_VDS", "[HDF5][Full]") { REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); } +// The per-image reflections and lattices are written in the setting the images were indexed in, but +// the space group - and with it the conventional setting the cell beside them is in - is only settled +// after the merge, so the two can differ by an integral change of basis. /entry/MX/reindexMatrix +// carries it, and the reader applies it, so what comes out is in the cell's setting. A file without +// the dataset (every file written before it existed) is read as the identity, which is what the +// round-trip tests above check. +TEST_CASE("JFJochReader_ReadReflections_Reindex", "[HDF5][Full]") { + DiffractionExperiment x(DetJF(1)); + + x.FilePrefix("read_reflections_reindex") + .ImagesPerTrigger(2) + .ImagesPerFile(1) + .OverwriteExistingFiles(true) + .BitDepthImage(16) + .PixelSigned(true) + .SetFileWriterFormat(FileWriterFormat::NXmxVDS) + .IndexingAlgorithm(IndexingAlgorithmEnum::FFT) + .Compression(CompressionAlgorithm::NO_COMPRESSION); + + // hkl_cell = M . hkl_written, det 2 - the size of step a primitive-to-centred re-seat takes. + const std::array M = {1, 1, 0, + 0, 1, 1, + 1, 0, 1}; + + std::vector image(x.GetPixelsNum(), 0); + + RegisterHDF5Filter(); + + { + StartMessage start_message; + x.FillMessage(start_message); + + FileWriter writer(start_message); + ScanResultGenerator scan_result(x); + + for (int i = 0; i < x.GetImageNum(); i++) { + DataMessage message{}; + message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum()); + message.number = i; + + if (i == 1) { + message.integrated_reflections = 2; + message.reflections = {MakeTestReflection(i, 0), MakeTestReflection(i, 1)}; + message.indexing_result = true; + message.indexing_lattice = CrystalLattice({100,0,0}, {0,50,0}, {0,0,30}); + } + + REQUIRE_NOTHROW(writer.WriteHDF5(message)); + scan_result.Add(message); + } + + EndMessage end_message; + end_message.max_image_number = x.GetImageNum(); + end_message.reindex_matrix = M; + scan_result.FillEndMessage(end_message); + + writer.WriteHDF5(end_message); + writer.Finalize(); + } + + // hkl and the lattice come back in the cell's setting; every other field is untouched. + const auto check = [&](const Reflection &got, int j) { + const Reflection want = MakeTestReflection(1, j); + INFO("reflection " << j); + CHECK(got.h == M[0] * want.h + M[1] * want.k + M[2] * want.l); + CHECK(got.k == M[3] * want.h + M[4] * want.k + M[5] * want.l); + CHECK(got.l == M[6] * want.h + M[7] * want.k + M[8] * want.l); + CHECK(got.I == Catch::Approx(want.I)); + CHECK(got.d == Catch::Approx(want.d)); + CHECK(got.image_number == Catch::Approx(want.image_number)); + }; + + { + JFJochHDF5Reader reader; + REQUIRE_NOTHROW(reader.ReadFile("read_reflections_reindex_master.h5")); + + REQUIRE(reader.GetDataset()->reindex_matrix.has_value()); + CHECK(reader.GetDataset()->reindex_matrix.value() == M); + + auto reflections = reader.ReadReflections(); + REQUIRE(reflections.size() == 2); + REQUIRE(reflections[1].reflections.size() == 2); + check(reflections[1].reflections[0], 0); + check(reflections[1].reflections[1], 1); + + // latt = M . latt_written, row by row: (100,0,0)+(0,50,0), (0,50,0)+(0,0,30), (100,0,0)+(0,0,30). + CHECK(reflections[1].latt.Vec0().x == Catch::Approx(100.0f)); + CHECK(reflections[1].latt.Vec0().y == Catch::Approx(50.0f)); + CHECK(reflections[1].latt.Vec1().y == Catch::Approx(50.0f)); + CHECK(reflections[1].latt.Vec1().z == Catch::Approx(30.0f)); + CHECK(reflections[1].latt.Vec2().x == Catch::Approx(100.0f)); + CHECK(reflections[1].latt.Vec2().z == Catch::Approx(30.0f)); + CHECK(reflections[1].latt.CalcVolume() == Catch::Approx(2.0 * 100 * 50 * 30)); + + // The per-image message path (the viewer's) is re-seated the same way. + auto reader_image = reader.LoadImage(1); + REQUIRE(reader_image); + REQUIRE(reader_image->ImageData().reflections.size() == 2); + check(reader_image->ImageData().reflections[0], 0); + REQUIRE(reader_image->ImageData().indexing_lattice); + CHECK(reader_image->ImageData().indexing_lattice->CalcVolume() + == Catch::Approx(2.0 * 100 * 50 * 30)); + } + + remove("read_reflections_reindex_master.h5"); + remove("read_reflections_reindex_data_000001.h5"); + remove("read_reflections_reindex_data_000002.h5"); + + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); +} + static std::vector MakeTestSpots(int i) { return { SpotToSave{ diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index 3e3aa720..52a2ae81 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -1103,6 +1103,14 @@ void NXmx::Finalize(const EndMessage &end) { if (end.rotation_lattice_type) SaveScalar(*hdf5_file, "/entry/MX/rotationLatticeNiggliClass", end.rotation_lattice_type->niggli_class); + // The setting the per-image reflections and lattices were written in is the one they were + // indexed in, which is not always the setting of the cell above: the space group is chosen + // after the images have gone to file, and choosing it can re-seat the lattice. Write the + // change of basis between the two so the file is self-describing - absent means they agree. + if (end.reindex_matrix) + SaveVector(*hdf5_file, "/entry/MX/reindexMatrix", + std::vector(end.reindex_matrix->begin(), end.reindex_matrix->end())); + if (end.indexing_rate) { SaveScalar(*hdf5_file, "/entry/MX/imageIndexedMean", end.indexing_rate.value()); }