JFJochHDF5Reader: Use lattice, and not Niggli class, to decode symmetry (if symmetry is imposed, then lattice search returns niggli class 0

reader/JFJochHDF5Reader.cpp

@@ -2,58 +2,62 @@
 // SPDX-License-Identifier: GPL-3.0-only
 
 #include <cmath>
+#include <set>
+
 #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<gemmi::CrystalSystem, char> 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<gemmi::CrystalSystem, char> 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<char> 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<hsize_t> GetDimension(HDF5Object &object, const std::string &path) {
@@ -914,6 +918,12 @@ bool JFJochHDF5Reader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &dataset
         if (tmp.size() == 9)
             message.indexing_lattice = CrystalLattice(tmp);
 
+        std::optional<std::string> lattice;
+        if (master_file->Exists("/entry/MX/bravaisLattice"))
+            lattice = master_file->ReadElement<std::string>("/entry/MX/bravaisLattice", image_number);
+        else
+            lattice = source_file->ReadElement<std::string>("/entry/MX/bravaisLattice", image_id);
+
         std::optional<uint32_t> niggli_opt;
         if (master_file->Exists("/entry/MX/niggli_class"))
             niggli_opt = master_file->ReadElement<uint32_t>("/entry/MX/niggli_class", image_number);
@@ -924,12 +934,12 @@ bool JFJochHDF5Reader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &dataset
         else if (source_file->Exists("/entry/MX/niggliClass"))
             niggli_opt = source_file->ReadElement<uint32_t>("/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<int64_t>(niggli_opt.value()),
+                .niggli_class = static_cast<int64_t>(niggli_opt.value_or(0)),
                 .crystal_system = symm_info.first,
             };
         }

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<int16_t> 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<JFJochReaderImage> 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);

writer/HDF5Objects.h

@@ -329,27 +329,7 @@ public:
     }
 
     template<class T>
-    std::optional<T> 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<T>();
-        }
-        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<hsize_t>(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<T> ReadElement(size_t n) const;
 
     std::string ReadString() const;
     void Close();
@@ -437,8 +417,81 @@ std::optional<T> HDF5Object::ReadElement(const std::string &name, size_t n) cons
     return std::nullopt;
 }
 
+template<class T>
+std::optional<T> 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<T>();
+    }
+    if (file_space.GetNumOfDimensions() != 1)
+        throw JFJochException(JFJochExceptionCategory::HDF5, "ReadElement requires a 1D dataset");
+    if (n >= dims[0])
+        return std::nullopt;
 
-template <class T> std::unique_ptr<HDF5DataSet> SaveVector(const HDF5Object& parent, const std::string &name, const std::vector<T> &val,
+    // Select a single-element hyperslab at position n.
+    file_space.SelectHyperslab({static_cast<hsize_t>(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<std::string> HDF5DataSet::ReadElement<std::string>(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<string> requires a 1D dataset");
+
+    if (n >= dims[0])
+        return std::nullopt;
+
+    file_space.SelectHyperslab({static_cast<hsize_t>(n)}, {1});
+
+    HDF5DataType file_type(*this);
+    if (H5Tget_class(file_type.GetID()) != H5T_STRING)
+        throw JFJochException(JFJochExceptionCategory::HDF5, "ReadElement<string>: 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<string> 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<char> 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<string> unsuccessful");
+
+    return std::string(buffer.data());
+}
+
+template <class T>
+std::unique_ptr<HDF5DataSet> SaveVector(const HDF5Object& parent, const std::string &name, const std::vector<T> &val,
                                                            std::vector<hsize_t> dims = {},
                                                            CompressionAlgorithm algorithm = CompressionAlgorithm::NO_COMPRESSION) {
     if (dims.empty()) dims = {val.size()};
@@ -459,7 +512,8 @@ template <class T> std::unique_ptr<HDF5DataSet> SaveVector(const HDF5Object& par
     return dataset;
 }
 
-template <class T> std::unique_ptr<HDF5DataSet> HDF5Object::SaveVector(const std::string &name, const std::vector<T> &val,
+template <class T>
+std::unique_ptr<HDF5DataSet> HDF5Object::SaveVector(const std::string &name, const std::vector<T> &val,
                                                                        std::vector<hsize_t> dims, CompressionAlgorithm algorithm) {
     if (val.empty())
         throw JFJochException(JFJochExceptionCategory::HDF5, "Cannot write empty vector");