diff --git a/reader/JFJochHDF5Reader.cpp b/reader/JFJochHDF5Reader.cpp index 66cb9a2a..723a70a0 100644 --- a/reader/JFJochHDF5Reader.cpp +++ b/reader/JFJochHDF5Reader.cpp @@ -2,58 +2,62 @@ // SPDX-License-Identifier: GPL-3.0-only #include +#include + #include "JFJochHDF5Reader.h" #include "spdlog/fmt/fmt.h" #include "../image_analysis/bragg_integration/CalcISigma.h" #include "../image_analysis/spot_finding/SpotUtils.h" #include "../common/GridScanSettings.h" -inline std::pair parse_niggli_class(int64_t val) { - switch (val) { - case 1: return {gemmi::CrystalSystem::Cubic, 'F'}; - case 2: return {gemmi::CrystalSystem::Trigonal, 'R'}; - case 3: return {gemmi::CrystalSystem::Cubic, 'P'}; - case 5: return {gemmi::CrystalSystem::Cubic, 'I'}; - case 4: return {gemmi::CrystalSystem::Trigonal, 'R'}; - case 6: return {gemmi::CrystalSystem::Tetragonal, 'I'}; - case 7: return {gemmi::CrystalSystem::Tetragonal, 'I'}; - case 9: return {gemmi::CrystalSystem::Trigonal, 'R'}; - case 10: return {gemmi::CrystalSystem::Monoclinic, 'C'}; - case 11: return {gemmi::CrystalSystem::Tetragonal, 'P'}; - case 12: return {gemmi::CrystalSystem::Hexagonal, 'P'}; - case 13: return {gemmi::CrystalSystem::Orthorhombic, 'C'}; - case 15: return {gemmi::CrystalSystem::Tetragonal, 'I'}; - case 16: return {gemmi::CrystalSystem::Orthorhombic, 'F'}; - case 14: return {gemmi::CrystalSystem::Monoclinic, 'C'}; - case 17: return {gemmi::CrystalSystem::Monoclinic, 'C'}; - case 18: return {gemmi::CrystalSystem::Tetragonal, 'I'}; - case 19: return {gemmi::CrystalSystem::Orthorhombic, 'I'}; - case 20: return {gemmi::CrystalSystem::Monoclinic, 'C'}; - case 21: return {gemmi::CrystalSystem::Tetragonal, 'P'}; - case 22: return {gemmi::CrystalSystem::Hexagonal, 'P'}; - case 23: return {gemmi::CrystalSystem::Orthorhombic, 'C'}; - case 24: return {gemmi::CrystalSystem::Trigonal, 'R'}; - case 25: return {gemmi::CrystalSystem::Monoclinic, 'C'}; - case 26: return {gemmi::CrystalSystem::Orthorhombic, 'F'}; - case 27: return {gemmi::CrystalSystem::Monoclinic, 'C'}; - case 28: return {gemmi::CrystalSystem::Monoclinic, 'C'}; - case 29: return {gemmi::CrystalSystem::Monoclinic, 'C'}; - case 30: return {gemmi::CrystalSystem::Monoclinic, 'C'}; - case 31: return {gemmi::CrystalSystem::Triclinic, 'P'}; - case 32: return {gemmi::CrystalSystem::Orthorhombic, 'P'}; - case 40: return {gemmi::CrystalSystem::Orthorhombic, 'C'}; - case 35: return {gemmi::CrystalSystem::Monoclinic, 'P'}; - case 36: return {gemmi::CrystalSystem::Orthorhombic, 'C'}; - case 33: return {gemmi::CrystalSystem::Monoclinic, 'P'}; - case 38: return {gemmi::CrystalSystem::Orthorhombic, 'C'}; - case 34: return {gemmi::CrystalSystem::Monoclinic, 'P'}; - case 42: return {gemmi::CrystalSystem::Orthorhombic, 'I'}; - case 41: return {gemmi::CrystalSystem::Monoclinic, 'C'}; - case 37: return {gemmi::CrystalSystem::Monoclinic, 'C'}; - case 39: return {gemmi::CrystalSystem::Monoclinic, 'C'}; - case 44: return {gemmi::CrystalSystem::Triclinic, 'P'}; - default: return {gemmi::CrystalSystem::Triclinic, 'P'}; +inline std::pair parse_bravais_lattice(const std::string &val) { + if (val.empty()) + return {gemmi::CrystalSystem::Triclinic, 'P'}; + + if (val.size() != 2) + throw JFJochException(JFJochExceptionCategory::HDF5, "Wrong Bravais lattice encoding"); + + gemmi::CrystalSystem cs; + char centering = val[1]; + std::set allowed_centering; + + switch (val[0]) { + case 'a': + cs = gemmi::CrystalSystem::Triclinic; + allowed_centering = {'P'}; + break; + case 'm': + cs = gemmi::CrystalSystem::Monoclinic; + allowed_centering = {'P', 'A', 'B', 'C'}; + break; + case 'o': + cs = gemmi::CrystalSystem::Orthorhombic; + allowed_centering = {'P', 'A', 'B', 'C', 'I', 'F'}; + break; + case 't': + cs = gemmi::CrystalSystem::Tetragonal; + allowed_centering = {'P', 'I'}; + break; + case 'h': + if (centering == 'P') + cs = gemmi::CrystalSystem::Hexagonal; + else if (centering == 'R') + cs = gemmi::CrystalSystem::Trigonal; + allowed_centering = {'P', 'R'}; + break; + case 'c': + cs = gemmi::CrystalSystem::Cubic; + allowed_centering = {'P', 'F', 'I'}; + break; + default: + // allowed_centering is empty and exception will be always thrown + break; } + + if (!allowed_centering.contains(centering)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Invalid lattice encoding " + val); + + return {cs, centering}; } std::vector GetDimension(HDF5Object &object, const std::string &path) { @@ -914,6 +918,12 @@ bool JFJochHDF5Reader::LoadImage_i(std::shared_ptr &dataset if (tmp.size() == 9) message.indexing_lattice = CrystalLattice(tmp); + std::optional lattice; + if (master_file->Exists("/entry/MX/bravaisLattice")) + lattice = master_file->ReadElement("/entry/MX/bravaisLattice", image_number); + else + lattice = source_file->ReadElement("/entry/MX/bravaisLattice", image_id); + std::optional niggli_opt; if (master_file->Exists("/entry/MX/niggli_class")) niggli_opt = master_file->ReadElement("/entry/MX/niggli_class", image_number); @@ -924,12 +934,12 @@ bool JFJochHDF5Reader::LoadImage_i(std::shared_ptr &dataset else if (source_file->Exists("/entry/MX/niggliClass")) niggli_opt = source_file->ReadElement("/entry/MX/niggliClass", image_id); - if (niggli_opt) { - auto symm_info = parse_niggli_class(niggli_opt.value()); + if (lattice && !lattice->empty()) { + auto symm_info = parse_bravais_lattice(lattice.value()); message.lattice_type = LatticeMessage{ .centering = symm_info.second, - .niggli_class = static_cast(niggli_opt.value()), + .niggli_class = static_cast(niggli_opt.value_or(0)), .crystal_system = symm_info.first, }; } diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index cd7e556f..d7635940 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -1038,6 +1038,79 @@ TEST_CASE("JFJochReader_NiggliClass", "[HDF5][Full]") { REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); } +TEST_CASE("JFJochReader_NiggliClass_VDS", "[HDF5][Full]") { + DiffractionExperiment x(DetJF(1)); + x.FilePrefix("test95").ImagesPerTrigger(4).OverwriteExistingFiles(true); + x.BitDepthImage(16).ImagesPerFile(1).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true) + .IndexingAlgorithm(IndexingAlgorithmEnum::FFT); + x.Compression(CompressionAlgorithm::NO_COMPRESSION); + + std::vector image(x.GetPixelsNum()); + + RegisterHDF5Filter(); + { + StartMessage start_message; + x.FillMessage(start_message); + FileWriter file_set(start_message); + + LatticeMessage lm{ + .centering = 'F', + .niggli_class = 1, + .crystal_system = gemmi::CrystalSystem::Cubic, + }; + + DataMessage message{}; + message.number = 0; + message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum()); + message.indexing_result = true; + message.indexing_lattice = CrystalLattice(40, 50, 60, 90, 90, 90); + message.lattice_type = lm; + REQUIRE_NOTHROW(file_set.WriteHDF5(message)); + + message.number = 1; + message.indexing_result = false; + message.indexing_lattice = std::nullopt; + REQUIRE_NOTHROW(file_set.WriteHDF5(message)); + + EndMessage end_message; + end_message.max_image_number = 2; + file_set.WriteHDF5(end_message); + file_set.Finalize(); + } + { + JFJochHDF5Reader reader; + REQUIRE_NOTHROW(reader.ReadFile("test95_master.h5")); + auto dataset = reader.GetDataset(); + CHECK(dataset->experiment.GetImageNum() == 2); + + std::shared_ptr reader_image, reader_image_2; + + REQUIRE_NOTHROW(reader_image = reader.LoadImage(0)); + + REQUIRE(reader_image); + CHECK(reader_image->ImageData().indexing_result.value() == true); + REQUIRE(reader_image->ImageData().indexing_lattice); + REQUIRE(reader_image->ImageData().lattice_type); + CHECK(reader_image->ImageData().lattice_type->centering == 'F'); + CHECK(reader_image->ImageData().lattice_type->niggli_class == 1); + CHECK(reader_image->ImageData().lattice_type->crystal_system == gemmi::CrystalSystem::Cubic); + + REQUIRE_NOTHROW(reader_image_2 = reader.LoadImage(1)); + + REQUIRE(reader_image_2); + CHECK(!reader_image_2->ImageData().indexing_result.value()); + REQUIRE(!reader_image_2->ImageData().indexing_lattice); + REQUIRE(!reader_image_2->ImageData().lattice_type); + } + remove("test95_master.h5"); + remove("test95_data_000001.h5"); + remove("test95_data_000002.h5"); + + // No leftover HDF5 objects + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); +} + + TEST_CASE("JFJochReader_MissingEntries", "[HDF5][Full]") { DiffractionExperiment x(DetJF(1)); x.FilePrefix("test96").ImagesPerTrigger(4).OverwriteExistingFiles(true); diff --git a/writer/HDF5Objects.h b/writer/HDF5Objects.h index bd0b170d..c2fe512b 100644 --- a/writer/HDF5Objects.h +++ b/writer/HDF5Objects.h @@ -329,27 +329,7 @@ public: } template - std::optional ReadElement(size_t n) const { - HDF5DataSpace file_space(*this); - auto dims = file_space.GetDimensions(); - // Allow scalar (0-D) to be read with n==0 as a convenience. - if (file_space.GetNumOfDimensions() == 0) { - if (n != 0) return std::nullopt; - return ReadScalar(); - } - if (file_space.GetNumOfDimensions() != 1) - throw JFJochException(JFJochExceptionCategory::HDF5, "ReadNth requires a 1D dataset"); - if (n >= dims[0]) - return std::nullopt; - - // Select a single-element hyperslab at position n. - file_space.SelectHyperslab({static_cast(n)}, {1}); - HDF5DataSpace mem_space({1}); - T out{}; - if (H5Dread(id, HDF5DataType(out).GetID(), mem_space.GetID(), file_space.GetID(), H5P_DEFAULT, &out) < 0) - throw JFJochException(JFJochExceptionCategory::HDF5, "ReadNth unsuccessful"); - return out; - } + std::optional ReadElement(size_t n) const; std::string ReadString() const; void Close(); @@ -437,8 +417,81 @@ std::optional HDF5Object::ReadElement(const std::string &name, size_t n) cons return std::nullopt; } +template +std::optional HDF5DataSet::ReadElement(size_t n) const { + HDF5DataSpace file_space(*this); + auto dims = file_space.GetDimensions(); + // Allow scalar (0-D) to be read with n==0 as a convenience. + if (file_space.GetNumOfDimensions() == 0) { + if (n != 0) return std::nullopt; + return ReadScalar(); + } + if (file_space.GetNumOfDimensions() != 1) + throw JFJochException(JFJochExceptionCategory::HDF5, "ReadElement requires a 1D dataset"); + if (n >= dims[0]) + return std::nullopt; -template std::unique_ptr SaveVector(const HDF5Object& parent, const std::string &name, const std::vector &val, + // Select a single-element hyperslab at position n. + file_space.SelectHyperslab({static_cast(n)}, {1}); + HDF5DataSpace mem_space({1}); + T out{}; + if (H5Dread(id, HDF5DataType(out).GetID(), mem_space.GetID(), file_space.GetID(), H5P_DEFAULT, &out) < 0) + throw JFJochException(JFJochExceptionCategory::HDF5, "ReadElement unsuccessful"); + return out; +} + +template<> +inline std::optional HDF5DataSet::ReadElement(size_t n) const { + HDF5DataSpace file_space(*this); + auto dims = file_space.GetDimensions(); + + if (file_space.GetNumOfDimensions() == 0) { + if (n != 0) return std::nullopt; + return ReadString(); + } + + if (file_space.GetNumOfDimensions() != 1) + throw JFJochException(JFJochExceptionCategory::HDF5, "ReadElement requires a 1D dataset"); + + if (n >= dims[0]) + return std::nullopt; + + file_space.SelectHyperslab({static_cast(n)}, {1}); + + HDF5DataType file_type(*this); + if (H5Tget_class(file_type.GetID()) != H5T_STRING) + throw JFJochException(JFJochExceptionCategory::HDF5, "ReadElement: dataset is not a string type"); + + // Variable-length string + if (H5Tis_variable_str(file_type.GetID()) > 0) { + HDF5DataSpace mem_space({1}); + char *tmp = nullptr; + + if (H5Dread(id, file_type.GetID(), mem_space.GetID(), file_space.GetID(), H5P_DEFAULT, &tmp) < 0) + throw JFJochException(JFJochExceptionCategory::HDF5, "ReadElement unsuccessful"); + + std::string out = (tmp != nullptr) ? std::string(tmp) : std::string{}; + if (tmp != nullptr) + H5free_memory(tmp); + return out; + } + + // Fixed-length string + const size_t elem_size = file_type.GetElemSize(); + if (elem_size == 0) + return std::string{}; + + std::vector buffer(elem_size + 1, '\0'); + HDF5DataSpace mem_space({1}); + + if (H5Dread(id, file_type.GetID(), mem_space.GetID(), file_space.GetID(), H5P_DEFAULT, buffer.data()) < 0) + throw JFJochException(JFJochExceptionCategory::HDF5, "ReadElement unsuccessful"); + + return std::string(buffer.data()); +} + +template +std::unique_ptr SaveVector(const HDF5Object& parent, const std::string &name, const std::vector &val, std::vector dims = {}, CompressionAlgorithm algorithm = CompressionAlgorithm::NO_COMPRESSION) { if (dims.empty()) dims = {val.size()}; @@ -459,7 +512,8 @@ template std::unique_ptr SaveVector(const HDF5Object& par return dataset; } -template std::unique_ptr HDF5Object::SaveVector(const std::string &name, const std::vector &val, +template +std::unique_ptr HDF5Object::SaveVector(const std::string &name, const std::vector &val, std::vector dims, CompressionAlgorithm algorithm) { if (val.empty()) throw JFJochException(JFJochExceptionCategory::HDF5, "Cannot write empty vector");