v1.0.0-rc.147 (#57)

This is an UNSTABLE release. The release has significant modifications for data processing - in case of troubles go back to 1.0.0-rc.144.

* jfjoch_viewer: Add reciprocal space viewer
* jfjoch_process: Two pass algorithm that does spot finding/indexing + integration of full dataset
* jfjoch_process: Improve logic for rotation indexer, to make execution more deterministic (still work in progress)

Reviewed-on: #57
Co-authored-by: Filip Leonarski <filip.leonarski@psi.ch>
Co-committed-by: Filip Leonarski <filip.leonarski@psi.ch>

tests/JFJochReaderTest.cpp

@@ -1038,6 +1038,84 @@ 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_0{};
+        message_0.number = 0;
+        message_0.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
+        message_0.indexing_result = true;
+        message_0.indexing_lattice = CrystalLattice(40, 50, 60, 90, 90, 90);
+        message_0.lattice_type = lm;
+        REQUIRE_NOTHROW(file_set.WriteHDF5(message_0));
+
+        DataMessage message_1{};
+        message_1.number = 1;
+        message_1.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
+        message_1.indexing_result = false;
+        message_1.indexing_lattice = std::nullopt;
+        REQUIRE_NOTHROW(file_set.WriteHDF5(message_1));
+
+        EndMessage end_message;
+        end_message.max_image_number = 2;
+        end_message.image_indexed = {true, false};
+
+        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);
+        REQUIRE(reader_image->ImageData().indexing_result.has_value());
+        CHECK(reader_image->ImageData().indexing_result.value() == true);
+        REQUIRE(reader_image->ImageData().indexing_lattice.has_value());
+        REQUIRE(reader_image->ImageData().lattice_type.has_value());
+        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);
@@ -2042,5 +2120,416 @@ TEST_CASE("JFJochReader_ReadReflections_VDS", "[HDF5][Full]") {
     remove("read_reflections_vds_data_000003.h5");
     remove("read_reflections_vds_data_000004.h5");
 
+    REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
+}
+
+static std::vector<SpotToSave> MakeTestSpots(int i) {
+    return {
+        SpotToSave{
+            .x = 1, .y = 2, .intensity = 376,
+            .ice_ring = true,
+            .indexed = true,
+            .h = 11, .k = -3, .l = -5,
+            .dist_ewald_sphere = 0.1234f
+        },
+        SpotToSave{
+            .x = 7, .y = static_cast<float>(-3 - i), .intensity = 0.156f,
+            .ice_ring = false,
+            .indexed = false,
+        }
+    };
+}
+
+// Assert the full field set on spots[0] and the per-image variation on
+// spots[1].y, which is the only field that differs across images.
+static void CheckSpotFields(const SpotToSave &s0, const SpotToSave &s1, int i) {
+    CHECK(s0.x == 1);
+    CHECK(s0.y == 2);
+    CHECK(s0.intensity == Catch::Approx(376));
+    CHECK(s0.ice_ring == true);
+    CHECK(s0.indexed == true);
+    CHECK(s0.h == 11);
+    CHECK(s0.k == -3);
+    CHECK(s0.l == -5);
+    CHECK(s0.dist_ewald_sphere == Catch::Approx(0.1234f));
+    CHECK(s0.image == i);
+
+    CHECK(s1.x == Catch::Approx(7));
+    CHECK(s1.y == Catch::Approx(static_cast<float>(-3 - i)));
+    CHECK(s1.intensity == Catch::Approx(0.156f));
+    CHECK(s1.ice_ring == false);
+    CHECK(s1.indexed == false);
+    CHECK(s1.image == i);
+}
+
+TEST_CASE("JFJochReader_ReadSpots_Legacy", "[HDF5][Full]") {
+    DiffractionExperiment x(DetJF(1));
+
+    x.FilePrefix("read_spots_legacy")
+     .ImagesPerTrigger(4)
+     .ImagesPerFile(1)
+     .OverwriteExistingFiles(true)
+     .BitDepthImage(16)
+     .PixelSigned(true)
+     .SetFileWriterFormat(FileWriterFormat::NXmxLegacy)
+     .IndexingAlgorithm(IndexingAlgorithmEnum::FFT)
+     .Compression(CompressionAlgorithm::NO_COMPRESSION);
+
+    std::vector<int16_t> image(x.GetPixelsNum(), 0);
+
+    RegisterHDF5Filter();
+    {
+        StartMessage start_message;
+        x.FillMessage(start_message);
+        FileWriter file_set(start_message);
+
+        for (int i = 0; i < x.GetImageNum(); i++) {
+            DataMessage message{};
+            message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
+            message.number = i;
+            message.spots = MakeTestSpots(i);
+            message.spot_count = 72 + i;
+            message.spot_count_ice_rings = 45 + 2 * i;
+            message.spot_count_low_res = 12 + 3 * i;
+            message.spot_count_indexed = 15 + 4 * i;
+            REQUIRE_NOTHROW(file_set.WriteHDF5(message));
+        }
+
+        EndMessage end_message;
+        end_message.max_image_number = x.GetImageNum();
+        file_set.WriteHDF5(end_message);
+        file_set.Finalize();
+    }
+
+    // All images, one at a time
+    {
+        JFJochHDF5Reader reader;
+        REQUIRE_NOTHROW(reader.ReadFile("read_spots_legacy_master.h5"));
+
+        for (int i = 0; i < 4; i++) {
+            std::vector<SpotToSave> spots;
+            REQUIRE_NOTHROW(spots = reader.ReadSpots(i));
+            REQUIRE(spots.size() == 2);
+            CheckSpotFields(spots[0], spots[1], i);
+        }
+    }
+
+    // Out-of-range must throw
+    {
+        JFJochHDF5Reader reader;
+        REQUIRE_NOTHROW(reader.ReadFile("read_spots_legacy_master.h5"));
+        REQUIRE_THROWS(reader.ReadSpots(4));
+    }
+
+    remove("read_spots_legacy_master.h5");
+    remove("read_spots_legacy_data_000001.h5");
+    remove("read_spots_legacy_data_000002.h5");
+    remove("read_spots_legacy_data_000003.h5");
+    remove("read_spots_legacy_data_000004.h5");
+
+    REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
+}
+
+TEST_CASE("JFJochReader_ReadSpots_VDS", "[HDF5][Full]") {
+    DiffractionExperiment x(DetJF(1));
+
+    x.FilePrefix("read_spots_vds")
+     .ImagesPerTrigger(4)
+     .ImagesPerFile(1)
+     .OverwriteExistingFiles(true)
+     .BitDepthImage(16)
+     .PixelSigned(true)
+     .SetFileWriterFormat(FileWriterFormat::NXmxVDS)
+     .IndexingAlgorithm(IndexingAlgorithmEnum::FFT)
+     .Compression(CompressionAlgorithm::NO_COMPRESSION);
+
+    std::vector<int16_t> image(x.GetPixelsNum(), 0);
+
+    RegisterHDF5Filter();
+    {
+        StartMessage start_message;
+        x.FillMessage(start_message);
+        FileWriter file_set(start_message);
+
+        for (int i = 0; i < x.GetImageNum(); i++) {
+            DataMessage message{};
+            message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
+            message.number = i;
+            message.spots = MakeTestSpots(i);
+            message.spot_count = 72 + i;
+            message.spot_count_ice_rings = 45 + 2 * i;
+            message.spot_count_low_res = 12 + 3 * i;
+            message.spot_count_indexed = 15 + 4 * i;
+            REQUIRE_NOTHROW(file_set.WriteHDF5(message));
+        }
+
+        EndMessage end_message;
+        end_message.max_image_number = x.GetImageNum();
+        file_set.WriteHDF5(end_message);
+        file_set.Finalize();
+    }
+
+    // All images, one at a time — also verifies that .image carries the
+    // correct global index across the virtual-to-source remapping.
+    {
+        JFJochHDF5Reader reader;
+        REQUIRE_NOTHROW(reader.ReadFile("read_spots_vds_master.h5"));
+
+        for (int i = 0; i < 4; i++) {
+            std::vector<SpotToSave> spots;
+            REQUIRE_NOTHROW(spots = reader.ReadSpots(i));
+            REQUIRE(spots.size() == 2);
+            CheckSpotFields(spots[0], spots[1], i);
+        }
+    }
+
+    // Out-of-range must throw
+    {
+        JFJochHDF5Reader reader;
+        REQUIRE_NOTHROW(reader.ReadFile("read_spots_vds_master.h5"));
+        REQUIRE_THROWS(reader.ReadSpots(4));
+    }
+
+    // Image with no spots returns an empty vector
+    {
+        JFJochHDF5Reader reader;
+        REQUIRE_NOTHROW(reader.ReadFile("read_spots_vds_master.h5"));
+
+        // Write a separate 2-image VDS file where only image 0 has spots.
+        DiffractionExperiment y(DetJF(1));
+        y.FilePrefix("read_spots_vds_sparse")
+         .ImagesPerTrigger(2)
+         .ImagesPerFile(1)
+         .OverwriteExistingFiles(true)
+         .BitDepthImage(16)
+         .PixelSigned(true)
+         .SetFileWriterFormat(FileWriterFormat::NXmxVDS)
+         .IndexingAlgorithm(IndexingAlgorithmEnum::FFT)
+         .Compression(CompressionAlgorithm::NO_COMPRESSION);
+
+        {
+            StartMessage start_message;
+            y.FillMessage(start_message);
+            FileWriter file_set(start_message);
+
+            for (int i = 0; i < 2; i++) {
+                DataMessage message{};
+                message.image = CompressedImage(image, y.GetXPixelsNum(), y.GetYPixelsNum());
+                message.number = i;
+                if (i == 0)
+                    message.spots = MakeTestSpots(i);
+                REQUIRE_NOTHROW(file_set.WriteHDF5(message));
+            }
+
+            EndMessage end_message;
+            end_message.max_image_number = 2;
+            file_set.WriteHDF5(end_message);
+            file_set.Finalize();
+        }
+
+        JFJochHDF5Reader sparse_reader;
+        REQUIRE_NOTHROW(sparse_reader.ReadFile("read_spots_vds_sparse_master.h5"));
+
+        std::vector<SpotToSave> spots_0, spots_1;
+        REQUIRE_NOTHROW(spots_0 = sparse_reader.ReadSpots(0));
+        REQUIRE_NOTHROW(spots_1 = sparse_reader.ReadSpots(1));
+
+        REQUIRE(spots_0.size() == 2);
+        CHECK(spots_0[0].image == 0);
+        CHECK(spots_1.empty());
+
+        remove("read_spots_vds_sparse_master.h5");
+        remove("read_spots_vds_sparse_data_000001.h5");
+        remove("read_spots_vds_sparse_data_000002.h5");
+    }
+
+    remove("read_spots_vds_master.h5");
+    remove("read_spots_vds_data_000001.h5");
+    remove("read_spots_vds_data_000002.h5");
+    remove("read_spots_vds_data_000003.h5");
+    remove("read_spots_vds_data_000004.h5");
+
+    REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
+}
+
+TEST_CASE("JFJochReader_ReadAllSpots_VDS", "[HDF5][Full]") {
+    DiffractionExperiment x(DetJF(1));
+
+    x.FilePrefix("read_spots_vds")
+     .ImagesPerTrigger(20)
+     .ImagesPerFile(3)
+     .OverwriteExistingFiles(true)
+     .BitDepthImage(16)
+     .PixelSigned(true)
+     .SetFileWriterFormat(FileWriterFormat::NXmxVDS)
+     .IndexingAlgorithm(IndexingAlgorithmEnum::FFT)
+     .Compression(CompressionAlgorithm::NO_COMPRESSION);
+
+    std::vector<int16_t> image(x.GetPixelsNum(), 0);
+
+    RegisterHDF5Filter();
+    {
+        StartMessage start_message;
+        x.FillMessage(start_message);
+        FileWriter file_set(start_message);
+
+        for (int i = 0; i < x.GetImageNum(); i++) {
+            DataMessage message{};
+            message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
+            message.number = i;
+            message.spots = MakeTestSpots(i);
+            message.spot_count = 72 + i;
+            message.spot_count_ice_rings = 45 + 2 * i;
+            message.spot_count_low_res = 12 + 3 * i;
+            message.spot_count_indexed = 15 + 4 * i;
+            REQUIRE_NOTHROW(file_set.WriteHDF5(message));
+        }
+
+        EndMessage end_message;
+        end_message.max_image_number = x.GetImageNum();
+        file_set.WriteHDF5(end_message);
+        file_set.Finalize();
+    }
+
+    {
+        JFJochHDF5Reader reader;
+        REQUIRE_NOTHROW(reader.ReadFile("read_spots_vds_master.h5"));
+
+        std::shared_ptr<JFJochReaderSpots> ret;
+        REQUIRE_NOTHROW(ret = reader.ReadAllSpots(1, 15, 2));
+
+        // 1,3,5,7,9,11,13,15
+        REQUIRE(ret);
+        REQUIRE(ret->start_image == 1);
+        REQUIRE(ret->stride == 2);
+        REQUIRE(ret->spots.size() == 8);
+
+        for (int i = 0; i < ret->spots.size(); i++) {
+            REQUIRE(ret->spots[i].size() == 2);
+            CheckSpotFields(ret->spots[i][0], ret->spots[i][1], 2 * i + 1);
+        }
+    }
+
+    {
+        JFJochHDF5Reader reader;
+        REQUIRE_NOTHROW(reader.ReadFile("read_spots_vds_master.h5"));
+        REQUIRE_THROWS(reader.ReadAllSpots(-5,0));
+        REQUIRE_THROWS(reader.ReadAllSpots(5,0));
+    }
+
+    remove("read_spots_vds_master.h5");
+    remove("read_spots_vds_data_000001.h5");
+    remove("read_spots_vds_data_000002.h5");
+    remove("read_spots_vds_data_000003.h5");
+    remove("read_spots_vds_data_000004.h5");
+
+    REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
+}
+
+
+
+TEST_CASE("JFJochReader_ReadSpots_Integrated", "[HDF5][Full]") {
+    DiffractionExperiment x(DetJF(1));
+
+    x.FilePrefix("read_spots_integrated")
+     .ImagesPerTrigger(4)
+     .OverwriteExistingFiles(true)
+     .BitDepthImage(16)
+     .PixelSigned(true)
+     .SetFileWriterFormat(FileWriterFormat::NXmxIntegrated)
+     .IndexingAlgorithm(IndexingAlgorithmEnum::FFT)
+     .Compression(CompressionAlgorithm::NO_COMPRESSION);
+
+    std::vector<int16_t> image(x.GetPixelsNum(), 0);
+
+    RegisterHDF5Filter();
+    {
+        StartMessage start_message;
+        x.FillMessage(start_message);
+        FileWriter file_set(start_message);
+
+        for (int i = 0; i < x.GetImageNum(); i++) {
+            DataMessage message{};
+            message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
+            message.number = i;
+            message.spots = MakeTestSpots(i);
+            message.spot_count = 72 + i;
+            message.spot_count_ice_rings = 45 + 2 * i;
+            message.spot_count_low_res = 12 + 3 * i;
+            message.spot_count_indexed = 15 + 4 * i;
+            REQUIRE_NOTHROW(file_set.WriteHDF5(message));
+        }
+
+        EndMessage end_message;
+        end_message.max_image_number = x.GetImageNum();
+        file_set.WriteHDF5(end_message);
+        file_set.Finalize();
+    }
+
+    // All images, one at a time
+    {
+        JFJochHDF5Reader reader;
+        REQUIRE_NOTHROW(reader.ReadFile("read_spots_integrated_master.h5"));
+
+        for (int i = 0; i < 4; i++) {
+            std::vector<SpotToSave> spots;
+            REQUIRE_NOTHROW(spots = reader.ReadSpots(i));
+            REQUIRE(spots.size() == 2);
+            CheckSpotFields(spots[0], spots[1], i);
+        }
+    }
+
+    // Out-of-range must throw
+    {
+        JFJochHDF5Reader reader;
+        REQUIRE_NOTHROW(reader.ReadFile("read_spots_integrated_master.h5"));
+        REQUIRE_THROWS(reader.ReadSpots(4));
+    }
+
+    // Image with no spots returns an empty vector
+    {
+        DiffractionExperiment y(DetJF(1));
+        y.FilePrefix("read_spots_integrated_sparse")
+         .ImagesPerTrigger(3)
+         .OverwriteExistingFiles(true)
+         .BitDepthImage(16)
+         .PixelSigned(true)
+         .SetFileWriterFormat(FileWriterFormat::NXmxIntegrated)
+         .IndexingAlgorithm(IndexingAlgorithmEnum::FFT)
+         .Compression(CompressionAlgorithm::NO_COMPRESSION);
+
+        {
+            StartMessage start_message;
+            y.FillMessage(start_message);
+            FileWriter file_set(start_message);
+
+            for (int i = 0; i < 3; i++) {
+                DataMessage message{};
+                message.image = CompressedImage(image, y.GetXPixelsNum(), y.GetYPixelsNum());
+                message.number = i;
+                if (i == 1)
+                    message.spots = MakeTestSpots(i);
+                REQUIRE_NOTHROW(file_set.WriteHDF5(message));
+            }
+
+            EndMessage end_message;
+            end_message.max_image_number = 3;
+            file_set.WriteHDF5(end_message);
+            file_set.Finalize();
+        }
+
+        JFJochHDF5Reader reader;
+        REQUIRE_NOTHROW(reader.ReadFile("read_spots_integrated_sparse_master.h5"));
+
+        CHECK(reader.ReadSpots(0).empty());
+        REQUIRE(reader.ReadSpots(1).size() == 2);
+        CheckSpotFields(reader.ReadSpots(1)[0], reader.ReadSpots(1)[1], 1);
+        CHECK(reader.ReadSpots(2).empty());
+
+        remove("read_spots_integrated_sparse_master.h5");
+    }
+
+    remove("read_spots_integrated_master.h5");
+
     REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
 }