mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases 72 Wiki Activity
Files
5cc2f811b33df7391574a6c41640e187d31bae8b
Jungfraujoch/tests/JFJochReaderTest.cpp
T
leonarski_f 54f2de31bb grid scan: say how the stationary spindle angle is stated, and pin it 
A grid scan is a set of stills at a stationary spindle, and the angle it stood
at is what relates one grid to another taken elsewhere on the circle. Users were
finding an all-zero omega in the file and concluding the angle could not be
recorded at all.

It already can, and has since the goniometer and the grid scan stopped being
alternatives: send the axis with step 0 and its start angle, and that angle is
written per image into the NXmx sample chain, read back by reader/, and taken by
dials.import as a set of stills. Measured on a generated 12-image grid: the
placeholder file carries omega = 0 x 12, the same file with the axis sent at
step 0 carries omega = 90 x 12, and dials.import reports "still: 1, sweep: 0"
for both. Nothing in the code needed changing, so nothing was; what was missing
was that nobody could tell, and that no test held the behaviour down.

So: the API and the HDF5 documentation now say it in as many words, and three
tests pin the three legs the value crosses - the OpenAPI request (which used to
drop the grid scan whenever an axis was present, unpinned until now), the CBOR
start message, and the file round trip.

Also corrects a claim two comments and the HDF5 page were making. NXmx can
express "no rotation" perfectly well - a sample may depend_on "." - so the
placeholder is not there for the standard's sake. It is there because dxtbx
cannot read a sample chain of translations alone: strip the rotation axis from a
grid scan master and dials.import dies in get_dxtbx_goniometer with a matmul
dimension mismatch. Recorded so nobody removes the placeholder on the strength
of the standard.

One thing the change does not fix, because it cannot: a stationary angle is
invisible to DIALS when a grid scan is present. dxtbx picks the first varying
axis as the scan axis, which is a grid translation, so the oscillation reads
(0, 0); and with exactly one rotation axis in the chain it builds a single-axis
goniometer whose fixed rotation is the identity, never consulting the angle. The
same angle IS visible when it is the only candidate (oscillation reads (90, 0))
or when a Smargon head puts a second rotation axis in the chain (the setting
rotation then carries it). The value is in the file and correct either way.
2026-09-02 11:04:29 +02:00

4101 lines
169 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <catch2/catch_all.hpp>
#include "../common/DiffractionExperiment.h"
#include "../common/ScanResultGenerator.h"
#include "../writer/FileWriter.h"
#include "../reader/JFJochHDF5Reader.h"
#include "../reader/JFJochCBFReader.h"
#include "../reader/MiniCBF.h"
#include "../compression/JFJochCompressor.h"
#include <fstream>
#include <future>
#include <iomanip>
#include <sstream>
TEST_CASE("HDF5DataType_Sign","[HDF5]") {
HDF5DataType type_u8((uint8_t)0), type_fl(0.0f), type_i32((int32_t) 0), type_u32((uint32_t) 0);
CHECK(!type_u8.IsSigned());
CHECK(type_fl.IsSigned());
CHECK(type_i32.IsSigned());
CHECK(!type_u32.IsSigned());
}
TEST_CASE("HDF5DataType_ElemSize","[HDF5]") {
HDF5DataType type_u8((uint8_t)0), type_fl(0.0f), type_i32((int32_t) 0), type_u32((uint32_t) 0);
CHECK(type_u8.GetElemSize() == 1);
CHECK(type_fl.GetElemSize() == 4);
CHECK(type_i32.GetElemSize() == 4);
CHECK(type_u32.GetElemSize() == 4);
}
TEST_CASE("HDF5DataType_ElemType","[HDF5]") {
HDF5DataType type_u8((uint8_t)0), type_fl(0.0f), type_i32((int32_t) 0), type_u32((uint32_t) 0);
CHECK(type_u8.IsInteger());
CHECK(!type_fl.IsInteger());
CHECK(type_fl.IsFloat());
CHECK(type_i32.IsInteger());
CHECK(type_u32.IsInteger());
}
TEST_CASE("JFJochReader_SpaceGroupSetting", "[HDF5][Full]") {
// A space group is carried through the master file as its name, not its number, because a number
// only ever names the reference setting. Both groups here are non-reference settings that a
// number destroys: "P 1 1 2" comes back from 3 as "P 1 2 1", and "R 3:R" from 146 as "R 3:H"
// (short_name() loses that one too - only xhm() is faithful).
const auto setting = GENERATE(std::string("P 1 1 2"), std::string("R 3:R"));
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(setting);
REQUIRE(sg != nullptr);
// The number is not a carrier for it, which is the whole reason the name is written.
DiffractionExperiment by_number(DetJF(1));
by_number.SpaceGroupNumber(sg->number);
CHECK(by_number.GetGemmiSpaceGroup()->xhm() != setting);
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test_sg_setting").ImagesPerTrigger(1).OverwriteExistingFiles(true);
x.SetSpaceGroup(*sg);
CHECK(x.GetGemmiSpaceGroup()->xhm() == setting);
CHECK(x.GetSpaceGroupNumber() == sg->number);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
EndMessage end_message;
end_message.max_image_number = 0;
end_message.space_group_name = sg->xhm();
std::unique_ptr<NXmx> master = std::make_unique<NXmx>(start_message);
master->Finalize(end_message);
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("test_sg_setting_master.h5"));
auto dataset = reader.GetDataset();
REQUIRE(dataset->experiment.GetGemmiSpaceGroup().has_value());
CHECK(dataset->experiment.GetGemmiSpaceGroup()->xhm() == setting);
}
remove("test_sg_setting_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_SpaceGroupNumberOnly", "[HDF5][Full]") {
// A file written before the name was recorded carries only the number; it still reads back, as
// the reference setting the number names.
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test_sg_number").ImagesPerTrigger(1).OverwriteExistingFiles(true);
x.SpaceGroupNumber(96);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
EndMessage end_message;
end_message.max_image_number = 0;
std::unique_ptr<NXmx> master = std::make_unique<NXmx>(start_message);
master->Finalize(end_message);
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("test_sg_number_master.h5"));
auto dataset = reader.GetDataset();
REQUIRE(dataset->experiment.GetGemmiSpaceGroup().has_value());
CHECK(dataset->experiment.GetGemmiSpaceGroup()->xhm() == "P 43 21 2");
}
remove("test_sg_number_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_MasterFile", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test08").ImagesPerTrigger(950).OverwriteExistingFiles(true);
x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
.IncidentEnergy_keV(WVL_1A_IN_KEV)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10))
.SetUnitCell(UnitCell{.a= 10, .b= 20, .c= 30, .alpha= 90, .beta= 101, .gamma = 90});
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
EndMessage end_message;
end_message.max_image_number = 0;
std::unique_ptr<NXmx> master = std::make_unique<NXmx>(start_message);
master->Finalize(end_message);
master.reset();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("test08_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetBeamX_pxl() == Catch::Approx(x.GetBeamX_pxl()));
CHECK(dataset->experiment.GetBeamY_pxl() == Catch::Approx(x.GetBeamY_pxl()));
CHECK(dataset->experiment.GetDetectorDistance_mm() == Catch::Approx(x.GetDetectorDistance_mm()));
CHECK(dataset->experiment.GetFrameTime() == x.GetFrameTime());
CHECK(dataset->experiment.GetFrameCountTime() == x.GetFrameCountTime());
CHECK(dataset->experiment.GetWavelength_A() == Catch::Approx(x.GetWavelength_A()));
CHECK(dataset->experiment.GetImageNum() == 0);
REQUIRE(dataset->experiment.GetUnitCell().has_value());
CHECK(dataset->experiment.GetUnitCell()->b == 20.0);
CHECK(dataset->experiment.GetUnitCell()->beta == 101.0);
CHECK(dataset->calibration_data.empty());
}
remove("test08_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_MasterFile_Calibration", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test_reader_calibration").ImagesPerTrigger(1).OverwriteExistingFiles(true);
RegisterHDF5Filter();
std::vector<uint16_t> calib_1(200*300, 10);
std::vector<int32_t> calib_2(100*400, 55);
std::vector<float> calib_f(100*400, 1234.56f);
JFJochBitShuffleCompressor compressor(CompressionAlgorithm::BSHUF_LZ4);
auto calib_1_compressed = compressor.Compress(calib_1);
{
StartMessage start_message;
x.FillMessage(start_message);
CompressedImage calibration_01(calib_1, 200, 300);
CompressedImage calibration_02(calib_2, 100, 400);
CompressedImage calibration_f(calib_f, 100, 400);
CompressedImage calibration_01_lz4(
calib_1_compressed.data(), calib_1_compressed.size(),
200, 300, CompressedImageMode::Uint16, CompressionAlgorithm::BSHUF_LZ4
);
calibration_01.Channel("c1");
calibration_02.Channel("c2");
calibration_f.Channel("cf");
calibration_01_lz4.Channel("c1_lz4");
EndMessage end_message;
end_message.max_image_number = 0;
std::unique_ptr<NXmx> master = std::make_unique<NXmx>(start_message);
master->WriteCalibration(calibration_01);
master->WriteCalibration(calibration_01_lz4);
master->WriteCalibration(calibration_02);
master->WriteCalibration(calibration_f);
master->Finalize(end_message);
master.reset();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("test_reader_calibration_master.h5"));
auto dataset = reader.GetDataset();
REQUIRE(dataset->calibration_data.size() == 4);
CHECK(dataset->calibration_data[0] == "c1");
CHECK(dataset->calibration_data[1] == "c1_lz4");
CHECK(dataset->calibration_data[2] == "c2");
CHECK(dataset->calibration_data[3] == "cf");
std::vector<uint8_t> buffer;
std::vector<uint8_t> buff_2;
REQUIRE_THROWS(reader.ReadCalibration(buffer, "c3"));
CompressedImage test;
REQUIRE_NOTHROW(test = reader.ReadCalibration(buffer, "c1"));
CHECK(test.GetByteDepth() == 2);
CHECK(test.GetHeight() == 300);
CHECK(test.GetWidth() == 200);
CHECK(test.GetMode() == CompressedImageMode::Uint16);
CHECK(reinterpret_cast<const uint16_t *>(test.GetUncompressedPtr(buff_2))[76] == 10);
REQUIRE_NOTHROW(test = reader.ReadCalibration(buffer, "c1_lz4"));
CHECK(test.GetByteDepth() == 2);
CHECK(test.GetHeight() == 300);
CHECK(test.GetWidth() == 200);
CHECK(test.GetMode() == CompressedImageMode::Uint16);
CHECK(reinterpret_cast<const uint16_t *>(test.GetUncompressedPtr(buff_2))[76] == 10);
REQUIRE_NOTHROW(test = reader.ReadCalibration(buffer, "c2"));
CHECK(test.GetByteDepth() == 4);
CHECK(test.GetHeight() == 400);
CHECK(test.GetWidth() == 100);
CHECK(test.GetMode() == CompressedImageMode::Int32);
CHECK(reinterpret_cast<const int32_t *>(test.GetUncompressedPtr(buff_2))[76] == 55);
REQUIRE_NOTHROW(test = reader.ReadCalibration(buffer, "cf"));
CHECK(test.GetByteDepth() == 4);
CHECK(test.GetHeight() == 400);
CHECK(test.GetWidth() == 100);
CHECK(test.GetMode() == CompressedImageMode::Float32);
CHECK(reinterpret_cast<const float *>(test.GetUncompressedPtr(buff_2))[76] == Catch::Approx(1234.56f));
}
remove("test_reader_calibration_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_DefaultExperiment", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test_def").OverwriteExistingFiles(true);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
EndMessage end_message;
end_message.max_image_number = 0;
std::unique_ptr<NXmx> master = std::make_unique<NXmx>(start_message);
master->Finalize(end_message);
master.reset();
}
{
JFJochHDF5Reader reader;
DiffractionExperiment x1;
IndexingSettings is;
is.FFT_NumVectors(1024);
x1.ImportIndexingSettings(is);
reader.Experiment(x1);
REQUIRE_NOTHROW(reader.ReadFile("test_def_master.h5"));
auto dataset = reader.GetDataset();
REQUIRE(x1.GetIndexingSettings().GetFFT_NumVectors() == 1024);
}
remove("test_def_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_PixelMask", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test16").ImagesPerTrigger(950).OverwriteExistingFiles(true);
x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
.IncidentEnergy_keV(WVL_1A_IN_KEV).PixelSigned(false).BitDepthImage(16)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
RegisterHDF5Filter();
std::vector<uint32_t> pixel_mask(x.GetPixelsNum(), 0);
pixel_mask[5767] = 1;
pixel_mask[x.GetPixelsNum() - 1] = 4;
pixel_mask[0] = 256;
pixel_mask[3] = 1u << PixelMask::BeamStopPixelBit;
ScanResultGenerator generator(x);
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
{
StartMessage start_message;
x.FillMessage(start_message);
start_message.pixel_mask["default"] = pixel_mask;
FileWriter file_set(start_message);
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = 0;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
REQUIRE_NOTHROW(generator.Add(message));
EndMessage end_message;
end_message.max_image_number = 1;
generator.FillEndMessage(end_message);
file_set.WriteHDF5(end_message);
file_set.Finalize();
}
{
JFJochHDF5Reader reader;
reader.ReadFile("test16_master.h5");
auto dataset = reader.GetDataset();
REQUIRE(dataset->pixel_mask->GetMask().size() == x.GetPixelsNum());
CHECK(dataset->pixel_mask->GetMask() == pixel_mask);
std::shared_ptr<JFJochReaderImage> reader_image;
REQUIRE_NOTHROW(reader_image = reader.LoadImage(0));
REQUIRE(reader_image);
CHECK(reader_image->Image().at(5767) == GAP_PXL_VALUE);
CHECK(reader_image->Image().at(0) == ERROR_PXL_VALUE);
CHECK(reader_image->Image().at(1) == 0);
CHECK(reader_image->Image().at(2) == 0);
// The beam-stop shadow reads back as its own marker, not as a bad pixel
CHECK(reader_image->Image().at(3) == BEAM_STOP_PXL_VALUE);
CHECK(reader_image->Image().at(x.GetPixelsNum() - 1) == ERROR_PXL_VALUE);
}
remove("test16_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_ROIDefinitions", "[HDF5][Full]") {
RegisterHDF5Filter();
// ROI definitions and the bitmap live in the master file for every format.
auto format = GENERATE(FileWriterFormat::NXmxVDS, FileWriterFormat::NXmxIntegrated);
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test_roi").ImagesPerTrigger(950).OverwriteExistingFiles(true)
.SetFileWriterFormat(format);
x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
.IncidentEnergy_keV(WVL_1A_IN_KEV).PixelSigned(false).BitDepthImage(16)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
ROIDefinition defs;
defs.boxes.emplace_back("mybox", 10, 20, 30, 40);
defs.circles.emplace_back("mycircle", 100, 200, 15);
defs.azimuthal.emplace_back("mywedge", 2.0f, 4.0f, 30.0f, 90.0f);
x.ROI().SetROI(defs);
ScanResultGenerator generator(x);
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
{
StartMessage start_message;
x.FillMessage(start_message);
start_message.rois = x.ROI().ExportMetadata();
start_message.roi_map = x.ExportROIMap();
FileWriter file_set(start_message);
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = 0;
for (const auto &name : {"mybox", "mycircle", "mywedge"})
message.roi[name] = ROIMessage{.sum = 100, .sum_square = 1000, .max_count = 50,
.pixels = 10, .x_weighted = 500, .y_weighted = 600};
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
EndMessage end_message;
end_message.max_image_number = 1;
generator.FillEndMessage(end_message);
file_set.WriteHDF5(end_message);
file_set.Finalize();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("test_roi_master.h5"));
auto dataset = reader.GetDataset();
const auto &rd = dataset->experiment.ROI().GetROIDefinition();
REQUIRE(rd.boxes.size() == 1);
REQUIRE(rd.circles.size() == 1);
REQUIRE(rd.azimuthal.size() == 1);
CHECK(rd.boxes[0].GetName() == "mybox");
CHECK(rd.boxes[0].GetXMin() == 10);
CHECK(rd.boxes[0].GetXMax() == 20);
CHECK(rd.circles[0].GetName() == "mycircle");
CHECK(rd.circles[0].GetRadius_pxl() == 15.0f);
CHECK(rd.azimuthal[0].GetName() == "mywedge");
CHECK(rd.azimuthal[0].HasPhi());
CHECK(rd.azimuthal[0].GetPhiMin_deg() == 30.0f);
// bitmap read back with the per-pixel footprint and the name->bit index
CHECK(dataset->roi_map.size() == x.GetXPixelsNumConv() * x.GetYPixelsNumConv());
CHECK(dataset->roi_bit_index.size() == 3);
CHECK(dataset->roi_bit_index.at("mybox") == 0);
// per-image ROI results surface from the master (VDS-linked for VDS format)
REQUIRE(dataset->roi.size() == 3);
auto it = std::find(dataset->roi.begin(), dataset->roi.end(), "mybox");
REQUIRE(it != dataset->roi.end());
const size_t idx = std::distance(dataset->roi.begin(), it);
CHECK(dataset->roi_sum.at(idx).at(0) == 100);
}
remove("test_roi_master.h5");
remove("test_roi_data_000001.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_Goniometer", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test17").ImagesPerTrigger(950).OverwriteExistingFiles(true);
x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
.IncidentEnergy_keV(WVL_1A_IN_KEV).PixelSigned(false).BitDepthImage(16)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
x.Goniometer(GoniometerAxis("omega", 95, 0.1f, Coord(0,-1,0),{}).ScreeningWedge(0.01f));
RegisterHDF5Filter();
ScanResultGenerator generator(x);
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
DataMessage message{};
for (int i = 0; i < 5; i++) {
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = i;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
REQUIRE_NOTHROW(generator.Add(message));
}
EndMessage end_message;
end_message.max_image_number = 5;
generator.FillEndMessage(end_message);
file_set.WriteHDF5(end_message);
file_set.Finalize();
}
{
JFJochHDF5Reader reader;
reader.ReadFile("test17_master.h5");
auto dataset = reader.GetDataset();
REQUIRE(!dataset->experiment.GetGridScan().has_value());
REQUIRE(dataset->experiment.GetGoniometer().has_value());
CHECK(dataset->experiment.GetGoniometer()->GetStart_deg() == 95.0);
CHECK(dataset->experiment.GetGoniometer()->GetIncrement_deg() == Catch::Approx(0.1f).margin(0.00001f));
CHECK(dataset->experiment.GetGoniometer()->GetWedge_deg() == Catch::Approx(0.01f).margin(0.00001f));
CHECK(dataset->experiment.GetGoniometer()->GetName() == "omega");
CHECK(dataset->experiment.GetGoniometer()->GetAxis().x == 0);
CHECK(dataset->experiment.GetGoniometer()->GetAxis().y == -1);
CHECK(dataset->experiment.GetGoniometer()->GetAxis().z == 0);
}
remove("test17_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// The axis name is free-form in the API, on the wire and in the writer - tests/CBORTest.cpp round
// trips one literally called "z". The reader used to look only for "omega", so a sweep recorded
// under any other name came back as stills, with nothing to indicate it. This is that case.
TEST_CASE("JFJochReader_Goniometer_NonOmegaName", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test17b").ImagesPerTrigger(950).OverwriteExistingFiles(true);
x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
.IncidentEnergy_keV(WVL_1A_IN_KEV).PixelSigned(false).BitDepthImage(16)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
x.Goniometer(GoniometerAxis("phi", 12, 0.2f, Coord(-1,0,0),{}));
RegisterHDF5Filter();
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
DataMessage message{};
for (int i = 0; i < 5; i++) {
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = i;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
}
EndMessage end_message;
end_message.max_image_number = 5;
file_set.WriteHDF5(end_message);
file_set.Finalize();
}
{
JFJochHDF5Reader reader;
reader.ReadFile("test17b_master.h5");
auto dataset = reader.GetDataset();
REQUIRE(dataset->experiment.GetGoniometer().has_value());
CHECK(dataset->experiment.GetGoniometer()->GetName() == "phi");
CHECK(dataset->experiment.GetGoniometer()->GetStart_deg() == 12.0);
CHECK(dataset->experiment.GetGoniometer()->GetIncrement_deg() == Catch::Approx(0.2f).margin(0.00001f));
CHECK(dataset->experiment.GetGoniometer()->IsScanning());
CHECK(dataset->experiment.GetGoniometer()->GetAxis().x == -1);
}
remove("test17b_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_GridScan", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.PixelSigned(false).BitDepthImage(16).OverwriteExistingFiles(true);
x.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
DatasetSettings d;
d.FilePrefix("test_reader_grid_scan").ImagesPerTrigger(5);
d.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
.PhotonEnergy_keV(WVL_1A_IN_KEV)
.GridScan(GridScanSettings(3, -7.5, 8.0, true, true));
x.ImportDatasetSettings(d);
RegisterHDF5Filter();
ScanResultGenerator generator(x);
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
DataMessage message{};
for (int i = 0; i < 5; i++) {
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = i;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
}
EndMessage end_message;
end_message.max_image_number = 5;
generator.FillEndMessage(end_message);
file_set.WriteHDF5(end_message);
file_set.Finalize();
}
{
JFJochHDF5Reader reader;
reader.ReadFile("test_reader_grid_scan_master.h5");
auto dataset = reader.GetDataset();
// A grid scan carries a stationary spindle: NXmx cannot say "no rotation", and a chain of
// translations alone is not readable (dxtbx raises on it). It must not read back as a sweep.
REQUIRE(dataset->experiment.GetGoniometer().has_value());
CHECK(!dataset->experiment.GetGoniometer()->IsScanning());
CHECK(dataset->experiment.GetGoniometer()->GetIncrement_deg() == 0.0f);
REQUIRE(dataset->experiment.GetGridScan().has_value());
CHECK(dataset->experiment.GetGridScan()->IsSnakeScan());
CHECK(dataset->experiment.GetGridScan()->IsVerticalScan());
CHECK(dataset->experiment.GetGridScan()->GetNFast() == 3);
CHECK(dataset->experiment.GetGridScan()->GetNSlow() == 2);
CHECK(dataset->experiment.GetGridScan()->GetNElem() == 6);
CHECK(dataset->experiment.GetGridScan()->GetGridStepX_um() == Catch::Approx(-7.5));
CHECK(dataset->experiment.GetGridScan()->GetGridStepY_um() == Catch::Approx(8.0));
}
{
// That placeholder spindle must carry one entry per image, not a scalar. Our own reader
// copes with either, so the check has to be on the stored shape: a third-party reader takes
// the image count from the innermost axis of the sample chain, and grid_scan_x/y are
// translations and are passed over - so with a scalar here the whole scan reads as one image.
hid_t file = H5Fopen("test_reader_grid_scan_master.h5", H5F_ACC_RDONLY, H5P_DEFAULT);
REQUIRE(file >= 0);
hid_t omega = H5Dopen2(file, "/entry/sample/transformations/omega", H5P_DEFAULT);
REQUIRE(omega >= 0);
hid_t space = H5Dget_space(omega);
CHECK(H5Sget_simple_extent_ndims(space) == 1);
hsize_t dim = 0;
H5Sget_simple_extent_dims(space, &dim, nullptr);
CHECK(dim == 5);
H5Sclose(space);
H5Dclose(omega);
H5Fclose(file);
}
remove("test_reader_grid_scan_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_DataI16", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test09").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());
image[0] = INT16_MAX;
image[1] = INT16_MIN;
image[2] = 456;
image[3] = -3456;
ScanResultGenerator generator(x);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
for (int i = 0; i < x.GetImageNum(); i++) {
std::vector<SpotToSave> spots;
image[5678] = i;
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.spots = spots;
message.indexing_result = (i % 2 == 0);
message.bkg_estimate = i * 345.6;
message.number = i;
message.profile_radius = 123.09;
generator.Add(message);
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
}
EndMessage end_message;
end_message.max_image_number = x.GetImageNum();
generator.FillEndMessage(end_message);
file_set.WriteHDF5(end_message);
file_set.Finalize();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("test09_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 4);
REQUIRE(dataset->spot_count.size() == 4);
REQUIRE(dataset->bkg_estimate.size() == 4);
REQUIRE(dataset->profile_radius.size() == 4);
REQUIRE_THROWS(reader.LoadImage(4));
std::shared_ptr<JFJochReaderImage> reader_image;
for (int i = 0; i < 4; i++) {
REQUIRE_NOTHROW(reader_image = reader.LoadImage(i));
REQUIRE(reader_image);
CHECK(reader_image->Image()[0] == SATURATED_PXL_VALUE);
CHECK(reader_image->Image()[1] == ERROR_PXL_VALUE);
CHECK(reader_image->Image()[2] == image[2]);
CHECK(reader_image->Image()[3] == image[3]);
CHECK(reader_image->Image()[5678] == i);
CHECK(dataset->indexing_result[i] == (i % 2 == 0));
CHECK(dataset->bkg_estimate[i] == Catch::Approx(i * 345.6));
CHECK(dataset->profile_radius[i] == Catch::Approx(123.09));
}
}
remove("test09_master.h5");
remove("test09_data_000001.h5");
remove("test09_data_000002.h5");
remove("test09_data_000003.h5");
remove("test09_data_000004.h5");
// No leftover HDF5 objects
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_DataI16_OldMasterFormat", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test15").ImagesPerTrigger(4).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(1).SetFileWriterFormat(FileWriterFormat::NXmxLegacy).PixelSigned(true)
.IndexingAlgorithm(IndexingAlgorithmEnum::FFT);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
std::vector<int16_t> image(x.GetPixelsNum());
image[0] = INT16_MAX;
image[1] = INT16_MIN;
image[2] = 456;
image[3] = -3456;
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
for (int i = 0; i < x.GetImageNum(); i++) {
std::vector<SpotToSave> spots;
image[5678] = i;
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.spots = spots;
message.indexing_result = (i % 2 == 0);
message.bkg_estimate = i * 345.6;
message.number = i;
message.profile_radius = 1.64;
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("test15_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 4);
REQUIRE(dataset->spot_count.size() == 4);
REQUIRE(dataset->bkg_estimate.size() == 4);
REQUIRE(dataset->profile_radius.size() == 4);
REQUIRE_THROWS(reader.LoadImage(4));
std::shared_ptr<JFJochReaderImage> reader_image;
for (int i = 0; i < 4; i++) {
REQUIRE_NOTHROW(reader_image = reader.LoadImage(i));
REQUIRE(reader_image);
CHECK(reader_image->Image()[0] == SATURATED_PXL_VALUE);
CHECK(reader_image->Image()[1] == ERROR_PXL_VALUE);
CHECK(reader_image->Image()[2] == image[2]);
CHECK(reader_image->Image()[3] == image[3]);
CHECK(reader_image->Image()[5678] == i);
CHECK(dataset->profile_radius[i] == Catch::Approx(1.64));
CHECK(dataset->indexing_result[i] == (i % 2 == 0));
CHECK(dataset->bkg_estimate[i] == Catch::Approx(i * 345.6));
}
}
remove("test15_master.h5");
remove("test15_data_000001.h5");
remove("test15_data_000002.h5");
remove("test15_data_000003.h5");
remove("test15_data_000004.h5");
// No leftover HDF5 objects
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_DataU16", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test10").ImagesPerTrigger(4).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(1).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(false)
.IndexingAlgorithm(IndexingAlgorithmEnum::FFT);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
std::vector<uint16_t> image(x.GetPixelsNum());
image[0] = UINT16_MAX;
image[1] = INT16_MAX;
image[2] = 456;
ScanResultGenerator generator(x);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
for (int i = 0; i < x.GetImageNum(); i++) {
std::vector<SpotToSave> spots;
image[5678] = i;
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.spots = spots;
message.indexing_result = (i % 2 == 0);
message.bkg_estimate = i * 345.6;
message.number = i;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
}
EndMessage end_message;
end_message.max_image_number = x.GetImageNum();
generator.FillEndMessage(end_message);
file_set.WriteHDF5(end_message);
file_set.Finalize();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("test10_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 4);
REQUIRE_THROWS(reader.LoadImage(4));
std::shared_ptr<JFJochReaderImage> reader_image;
for (int i = 0; i < 4; i++) {
REQUIRE_NOTHROW(reader_image = reader.LoadImage(i));
REQUIRE(reader_image);
CHECK(reader_image->Image()[0] == SATURATED_PXL_VALUE);
CHECK(reader_image->Image()[1] == INT16_MAX);
CHECK(reader_image->Image()[2] == 456);
CHECK(reader_image->Image()[5678] == i);
CHECK(dataset->indexing_result[i] == (i % 2 == 0));
CHECK(dataset->bkg_estimate[i] == Catch::Approx(i * 345.6));
}
}
remove("test10_master.h5");
remove("test10_data_000001.h5");
remove("test10_data_000002.h5");
remove("test10_data_000003.h5");
remove("test10_data_000004.h5");
// No leftover HDF5 objects
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_DataI32", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test11").ImagesPerTrigger(4).OverwriteExistingFiles(true);
x.BitDepthImage(32).ImagesPerFile(1).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
std::vector<int32_t> image(x.GetPixelsNum());
image[0] = INT32_MAX;
image[1] = INT32_MIN;
image[2] = 456;
ScanResultGenerator generator(x);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
for (int i = 0; i < x.GetImageNum(); i++) {
std::vector<SpotToSave> spots;
image[5678] = i;
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.spots = spots;
message.number = i;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
}
EndMessage end_message;
end_message.max_image_number = x.GetImageNum();
generator.FillEndMessage(end_message);
file_set.WriteHDF5(end_message);
file_set.Finalize();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("test11_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 4);
REQUIRE_THROWS(reader.LoadImage(4));
std::shared_ptr<JFJochReaderImage> reader_image;
for (int i = 0; i < 4; i++) {
REQUIRE_NOTHROW(reader_image = reader.LoadImage(i));
REQUIRE(reader_image);
CHECK(reader_image->Image()[0] == SATURATED_PXL_VALUE);
CHECK(reader_image->Image()[1] == ERROR_PXL_VALUE);
CHECK(reader_image->Image()[2] == 456);
CHECK(reader_image->Image()[5678] == i);
}
}
remove("test11_master.h5");
remove("test11_data_000001.h5");
remove("test11_data_000002.h5");
remove("test11_data_000003.h5");
remove("test11_data_000004.h5");
// No leftover HDF5 objects
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_DataU32", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test12").ImagesPerTrigger(4).OverwriteExistingFiles(true);
x.BitDepthImage(32).ImagesPerFile(1).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(false);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
std::vector<uint32_t> image(x.GetPixelsNum());
image[0] = UINT32_MAX;
image[1] = static_cast<uint32_t>(INT32_MAX) + 50;
image[2] = 456;
image[3] = INT32_MAX;
image[4] = INT32_MAX - 1;
ScanResultGenerator generator(x);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
for (int i = 0; i < x.GetImageNum(); i++) {
std::vector<SpotToSave> spots;
image[5678] = i;
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.spots = spots;
message.number = i;
generator.Add(message);
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
}
EndMessage end_message;
end_message.max_image_number = x.GetImageNum();
generator.FillEndMessage(end_message);
file_set.WriteHDF5(end_message);
file_set.Finalize();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("test12_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 4);
REQUIRE_THROWS(reader.LoadImage(4));
std::shared_ptr<JFJochReaderImage> reader_image;
for (int i = 0; i < 4; i++) {
REQUIRE_NOTHROW(reader_image = reader.LoadImage(i));
REQUIRE(reader_image);
CHECK(reader_image->Image()[0] == INT32_MAX);
CHECK(reader_image->Image()[1] == INT32_MAX);
CHECK(reader_image->Image()[2] == 456);
CHECK(reader_image->Image()[3] == INT32_MAX);
CHECK(reader_image->Image()[4] == INT32_MAX - 1);
CHECK(reader_image->Image()[5678] == i);
}
}
remove("test12_master.h5");
remove("test12_data_000001.h5");
remove("test12_data_000002.h5");
remove("test12_data_000003.h5");
remove("test12_data_000004.h5");
// No leftover HDF5 objects
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_Summation", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test30").ImagesPerTrigger(3).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(3).SetFileWriterFormat(FileWriterFormat::NXmxVDS)
.PixelSigned(true);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
std::vector<int16_t> image_1(x.GetPixelsNum(),1);
std::vector<int16_t> image_2(x.GetPixelsNum(),2);
std::vector<int16_t> image_3(x.GetPixelsNum(),3);
image_3[0] = INT16_MAX;
image_2[1] = INT16_MIN;
ScanResultGenerator generator(x);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
std::vector<SpotToSave> spots;
DataMessage message{};
message.spots = spots;
message.image = CompressedImage(image_1, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = 0;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
message.image = CompressedImage(image_2, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = 1;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
message.image = CompressedImage(image_3, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = 2;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
EndMessage end_message;
end_message.max_image_number = x.GetImageNum();
generator.FillEndMessage(end_message);
file_set.WriteHDF5(end_message);
file_set.Finalize();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("test30_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 3);
std::shared_ptr<JFJochReaderImage> reader_image;
REQUIRE_NOTHROW(reader_image = reader.LoadImage(0, 3));
REQUIRE(reader_image);
CHECK(reader_image->Image()[0] == SATURATED_PXL_VALUE);
CHECK(reader_image->Image()[1] == ERROR_PXL_VALUE);
CHECK(reader_image->Image()[2] == 1 + 2 +3);
CHECK(reader_image->Image()[5678] == 1 + 2 +3);
CHECK(reader_image->Image()[x.GetPixelsNum() - 1] == 1 + 2 +3);
REQUIRE_THROWS(reader.LoadImage(1, 3));
}
remove("test30_master.h5");
remove("test30_data_000001.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_Summation_5", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test31").ImagesPerTrigger(5).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(5).SetFileWriterFormat(FileWriterFormat::NXmxVDS)
.PixelSigned(true);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
std::vector<int16_t> image_1(x.GetPixelsNum(),1);
std::vector<int16_t> image_2(x.GetPixelsNum(),2);
std::vector<int16_t> image_3(x.GetPixelsNum(),3);
std::vector<int16_t> image_4(x.GetPixelsNum(),4);
std::vector<int16_t> image_5(x.GetPixelsNum(),5);
image_3[0] = INT16_MAX;
image_2[1] = INT16_MIN;
ScanResultGenerator generator(x);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
std::vector<SpotToSave> spots;
DataMessage message{};
message.spots = spots;
message.image = CompressedImage(image_1, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = 0;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
message.image = CompressedImage(image_2, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = 1;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
message.image = CompressedImage(image_3, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = 2;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
message.image = CompressedImage(image_4, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = 3;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
message.image = CompressedImage(image_5, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = 4;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
EndMessage end_message;
end_message.max_image_number = x.GetImageNum();
generator.FillEndMessage(end_message);
file_set.WriteHDF5(end_message);
file_set.Finalize();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("test31_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 5);
std::shared_ptr<JFJochReaderImage> reader_image;
REQUIRE_NOTHROW(reader_image = reader.LoadImage(0, 5));
REQUIRE(reader_image);
CHECK(reader_image->Image()[0] == SATURATED_PXL_VALUE);
CHECK(reader_image->Image()[1] == ERROR_PXL_VALUE);
CHECK(reader_image->Image()[2] == 1 + 2 + 3 + 4 + 5);
CHECK(reader_image->Image()[5678] == 1 + 2 + 3 + 4 + 5);
CHECK(reader_image->Image()[x.GetPixelsNum() - 1] == 1 + 2 + 3 + 4 + 5);
REQUIRE_THROWS(reader.LoadImage(1, 6));
}
remove("test31_master.h5");
remove("test31_data_000001.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_Azint", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test27").ImagesPerTrigger(4).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(1).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(false);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
AzimuthalIntegrationSettings azint_settings;
azint_settings.AzimuthalBinCount(4);
x.ImportAzimuthalIntegrationSettings(azint_settings);
// The high-q limit is unset, i.e. "as far as the detector reaches", so read the settings back from
// the experiment, where that has been resolved against the geometry - that is what the bins are.
azint_settings = x.GetAzimuthalIntegrationSettings();
std::vector<uint16_t> image(x.GetPixelsNum());
AzimuthalIntegrationMapping azint(x, PixelMask(x));
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
start_message.az_int_bin_to_q = azint.GetBinToQ();
start_message.az_int_bin_to_phi = azint.GetBinToPhi();
start_message.az_int_q_bin_count = azint.GetQBinCount();
start_message.az_int_phi_bin_count = azint.GetAzimuthalBinCount();
FileWriter file_set(start_message);
for (int i = 0; i < x.GetImageNum(); i++) {
std::vector<SpotToSave> spots;
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.spots = spots;
message.number = i;
message.az_int_profile = std::vector<float>(azint_settings.GetBinCount(), 57);
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("test27_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 4);
std::shared_ptr<JFJochReaderImage> reader_image;
REQUIRE_NOTHROW(reader_image = reader.LoadImage(0));
REQUIRE(reader_image);
CHECK(reader_image->Dataset().az_int_bin_to_q.size() == azint_settings.GetBinCount());
CHECK(reader_image->Dataset().azimuthal_bins == azint_settings.GetAzimuthalBinCount());
CHECK(reader_image->Dataset().q_bins == azint_settings.GetQBinCount());
REQUIRE(reader_image->ImageData().az_int_profile.size() == azint_settings.GetBinCount());
CHECK(reader_image->ImageData().az_int_profile[23] == 57);
CHECK(reader_image->GetAzInt1D_BinToQ().size() == azint_settings.GetQBinCount());
REQUIRE(reader_image->GetAzInt1D().size() == azint_settings.GetQBinCount());
CHECK(reader_image->GetAzInt1D()[23] == 4 * 57);
}
remove("test27_master.h5");
remove("test27_data_000001.h5");
remove("test27_data_000002.h5");
remove("test27_data_000003.h5");
remove("test27_data_000004.h5");
// No leftover HDF5 objects
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_NiggliClass", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test95").ImagesPerTrigger(4).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(1).SetFileWriterFormat(FileWriterFormat::NXmxLegacy).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_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);
x.BitDepthImage(16).ImagesPerFile(10).SetFileWriterFormat(FileWriterFormat::NXmxLegacy).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);
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.spot_count_indexed = 56;
message.spot_count = 85;
message.b_factor = 123.45;
message.spots = {SpotToSave{.x = 10, .y=50, .intensity = 80}};
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
message.number = 1;
message.indexing_result = false;
message.indexing_lattice = std::nullopt;
message.spot_count_indexed = std::nullopt;
message.spot_count = 70;
message.b_factor = std::nullopt;
message.spots = {SpotToSave{.x = 10, .y=50, .intensity = 80}};
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("test96_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 2);
REQUIRE(dataset->b_factor.size() == 2);
REQUIRE(dataset->spot_count_indexed.size() == 2);
CHECK(dataset->b_factor[0] == Catch::Approx(123.45));
CHECK(std::isnan(dataset->b_factor[1]));
CHECK(dataset->spot_count_indexed[0] == 56);
CHECK(dataset->spot_count_indexed[1] == 0);
std::shared_ptr<JFJochReaderImage> reader_image, reader_image_2;
REQUIRE_NOTHROW(reader_image = reader.LoadImage(0));
REQUIRE(reader_image);
REQUIRE(reader_image->ImageData().b_factor.has_value());
CHECK(reader_image->ImageData().b_factor.value() == Catch::Approx(123.45));
REQUIRE(reader_image->ImageData().spot_count_indexed.has_value());
CHECK(reader_image->ImageData().spot_count_indexed.value() == 56);
REQUIRE_NOTHROW(reader_image_2 = reader.LoadImage(1));
REQUIRE(reader_image_2);
CHECK(reader_image_2->ImageData().spot_count_indexed.has_value());
CHECK(reader_image_2->ImageData().spot_count_indexed.value() == 0);
}
remove("test96_master.h5");
remove("test96_data_000001.h5");
// No leftover HDF5 objects
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_Spots_OldMasterFormat", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test40").ImagesPerTrigger(4).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(1).SetFileWriterFormat(FileWriterFormat::NXmxLegacy).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);
for (int i = 0; i < x.GetImageNum(); i++) {
std::vector<SpotToSave> spots;
spots.push_back(SpotToSave{
.x = 1, .y = 2, .intensity = 376,
.h = 11, .k = -3, .l = -5,
.dist_ewald_sphere = 0.1234f,
.ice_ring = true,
.indexed = true
});
spots.push_back(SpotToSave{
.x = 7, .y = -3, .intensity = 0.156f,
.ice_ring = false,
.indexed = false,
});
image[5678] = i;
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.spots = spots;
message.indexing_result = (i % 2 == 0);
message.number = 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("test40_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 4);
CHECK(dataset->spot_count[1] == 72 + 1);
CHECK(dataset->spot_count_ice_rings[2] == 45 + 2 * 2);
CHECK(dataset->spot_count_low_res[3] == 12 + 3 * 3);
CHECK(dataset->spot_count_indexed[0] == 15);
REQUIRE_THROWS(reader.LoadImage(4));
std::shared_ptr<JFJochReaderImage> reader_image;
for (int i = 0; i < 4; i++) {
REQUIRE_NOTHROW(reader_image = reader.LoadImage(i));
REQUIRE(reader_image);
CHECK(reader_image->ImageData().spot_count == 72 + i);
CHECK(reader_image->ImageData().spot_count_ice_rings == 45 + 2 * i);
CHECK(reader_image->ImageData().spot_count_low_res == 12 + 3 * i);
CHECK(reader_image->ImageData().spot_count_indexed == 15 + 4 * i);
REQUIRE(reader_image->ImageData().spots.size() == 2);
CHECK(reader_image->ImageData().spots[0].x == 1);
CHECK(reader_image->ImageData().spots[0].y == 2);
CHECK(reader_image->ImageData().spots[0].intensity == 376);
CHECK(reader_image->ImageData().spots[0].ice_ring == true);
CHECK(reader_image->ImageData().spots[0].indexed == true);
CHECK(reader_image->ImageData().spots[0].h == 11);
CHECK(reader_image->ImageData().spots[0].k == -3);
CHECK(reader_image->ImageData().spots[0].l == -5);
CHECK(reader_image->ImageData().spots[0].dist_ewald_sphere == Catch::Approx(0.1234f));
CHECK(reader_image->ImageData().spots[1].x == 7);
CHECK(reader_image->ImageData().spots[1].y == -3);
CHECK(reader_image->ImageData().spots[1].intensity == Catch::Approx(0.156f));
CHECK(reader_image->ImageData().spots[1].ice_ring == false);
CHECK(reader_image->ImageData().spots[1].indexed == false);
}
}
remove("test40_master.h5");
remove("test40_data_000001.h5");
remove("test40_data_000002.h5");
remove("test40_data_000003.h5");
remove("test40_data_000004.h5");
// No leftover HDF5 objects
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_Spots_VDS", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test41").ImagesPerTrigger(4).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(1).SetFileWriterFormat(FileWriterFormat::NXmxLegacy).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);
for (int i = 0; i < x.GetImageNum(); i++) {
std::vector<SpotToSave> spots;
spots.push_back(SpotToSave{
.x = 1, .y = 2, .intensity = 376,
.h = 11, .k = -3, .l = -5,
.dist_ewald_sphere = 0.1234f,
.ice_ring = true,
.indexed = true
});
spots.push_back(SpotToSave{
.x = 7, .y = -3, .intensity = 0.156f,
.ice_ring = false,
.indexed = false,
});
image[5678] = i;
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.spots = spots;
message.indexing_result = (i % 2 == 0);
message.number = 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("test41_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 4);
CHECK(dataset->spot_count[1] == 72 + 1);
CHECK(dataset->spot_count_ice_rings[2] == 45 + 2 * 2);
CHECK(dataset->spot_count_low_res[3] == 12 + 3 * 3);
CHECK(dataset->spot_count_indexed[0] == 15);
REQUIRE_THROWS(reader.LoadImage(4));
std::shared_ptr<JFJochReaderImage> reader_image;
for (int i = 0; i < 4; i++) {
REQUIRE_NOTHROW(reader_image = reader.LoadImage(i));
REQUIRE(reader_image);
CHECK(reader_image->ImageData().spot_count == 72 + i);
CHECK(reader_image->ImageData().spot_count_ice_rings == 45 + 2 * i);
CHECK(reader_image->ImageData().spot_count_low_res == 12 + 3 * i);
CHECK(reader_image->ImageData().spot_count_indexed == 15 + 4 * i);
REQUIRE(reader_image->ImageData().spots.size() == 2);
CHECK(reader_image->ImageData().spots[0].x == 1);
CHECK(reader_image->ImageData().spots[0].y == 2);
CHECK(reader_image->ImageData().spots[0].intensity == 376);
CHECK(reader_image->ImageData().spots[0].ice_ring == true);
CHECK(reader_image->ImageData().spots[0].indexed == true);
CHECK(reader_image->ImageData().spots[0].h == 11);
CHECK(reader_image->ImageData().spots[0].k == -3);
CHECK(reader_image->ImageData().spots[0].l == -5);
CHECK(reader_image->ImageData().spots[0].dist_ewald_sphere == Catch::Approx(0.1234f));
CHECK(reader_image->ImageData().spots[1].x == 7);
CHECK(reader_image->ImageData().spots[1].y == -3);
CHECK(reader_image->ImageData().spots[1].intensity == Catch::Approx(0.156f));
CHECK(reader_image->ImageData().spots[1].ice_ring == false);
CHECK(reader_image->ImageData().spots[1].indexed == false);
}
}
remove("test41_master.h5");
remove("test41_data_000001.h5");
remove("test41_data_000002.h5");
remove("test41_data_000003.h5");
remove("test41_data_000004.h5");
// No leftover HDF5 objects
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_InstrumentMetadata_Sample_RingCurrent", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
// Set identifying names and ring current (mA in API; writer stores A)
x.FilePrefix("test_meta").ImagesPerTrigger(0).OverwriteExistingFiles(true);
InstrumentMetadata metadata;
metadata.InstrumentName("PXI").SourceName("SLS");
x.ImportInstrumentMetadata(metadata).SampleName("test_sample").RingCurrent_mA(399.5); // 0.3995 A
x.TotalFlux(1e7).AttenuatorTransmission(0.56);
x.DetectIceRings(false);
// Minimal other required fields
x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
.IncidentEnergy_keV(WVL_1A_IN_KEV)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
// also set fluorescence spectrum
x.FluorescenceSpectrum(XrayFluorescenceSpectrum({1.0f, 2.0f, 3.0f}, {5.0f, 7.0f, 6.0f}));
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
EndMessage end_message;
end_message.max_image_number = 0;
std::unique_ptr<NXmx> master = std::make_unique<NXmx>(start_message);
master->Finalize(end_message);
master.reset();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("test_meta_master.h5"));
auto dataset = reader.GetDataset();
auto meta = dataset->experiment.GetInstrumentMetadata();
CHECK(meta.GetInstrumentName() == "PXI");
CHECK(meta.GetSourceName() == "SLS");
CHECK(dataset->experiment.GetAttenuatorTransmission() == Catch::Approx(0.56));
CHECK(dataset->experiment.GetTotalFlux() == Catch::Approx(1e7));
// Sample name
CHECK(dataset->experiment.GetSampleName() == "test_sample");
// Ring current read back in mA; allow small fp tolerance
CHECK(dataset->experiment.GetRingCurrent_mA().has_value());
CHECK(dataset->experiment.GetRingCurrent_mA().value() == Catch::Approx(399.5));
CHECK(!dataset->experiment.IsDetectIceRings());
// Fluorescence spectrum presence and values
REQUIRE(!dataset->experiment.GetFluorescenceSpectrum().empty());
const auto& fl = dataset->experiment.GetFluorescenceSpectrum();
CHECK(fl.GetEnergy_eV().size() == 3);
CHECK(fl.GetData().size() == 3);
CHECK(fl.GetEnergy_eV()[0] == Catch::Approx(1.0f));
CHECK(fl.GetData()[1] == Catch::Approx(7.0f));
}
remove("test_meta_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_NXmxIntegrated", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test_reader_integrated").ImagesPerTrigger(3).OverwriteExistingFiles(true);
x.BitDepthImage(16).PixelSigned(false).SetFileWriterFormat(FileWriterFormat::NXmxIntegrated);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
.IncidentEnergy_keV(WVL_1A_IN_KEV)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
AzimuthalIntegrationSettings azint_settings;
azint_settings.AzimuthalBinCount(4);
x.ImportAzimuthalIntegrationSettings(azint_settings);
// The high-q limit is unset, i.e. "as far as the detector reaches", so read the settings back from
// the experiment, where that has been resolved against the geometry - that is what the bins are.
azint_settings = x.GetAzimuthalIntegrationSettings();
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
image[0] = UINT16_MAX;
image[1] = 123;
image[5678] = 321;
AzimuthalIntegrationMapping azint(x, PixelMask(x));
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
start_message.az_int_bin_to_q = azint.GetBinToQ();
start_message.az_int_bin_to_phi = azint.GetBinToPhi();
start_message.az_int_q_bin_count = azint.GetQBinCount();
start_message.az_int_phi_bin_count = azint.GetAzimuthalBinCount();
FileWriter file_set(start_message);
for (int i = 0; i < x.GetImageNum(); i++) {
DataMessage message{};
image[5678] = 321 + i;
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = i;
message.image_collection_efficiency = 0.9f + 0.01f * i;
message.az_int_profile = std::vector<float>(azint_settings.GetBinCount(), static_cast<float>(50 + 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("test_reader_integrated_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 3);
REQUIRE(dataset->efficiency.size() == 3);
CHECK(dataset->efficiency[0] == Catch::Approx(0.90f));
CHECK(dataset->efficiency[1] == Catch::Approx(0.91f));
CHECK(dataset->efficiency[2] == Catch::Approx(0.92f));
CHECK(dataset->az_int_bin_to_q.size() == azint_settings.GetBinCount());
CHECK(dataset->azimuthal_bins == azint_settings.GetAzimuthalBinCount());
CHECK(dataset->q_bins == azint_settings.GetQBinCount());
std::shared_ptr<JFJochReaderImage> reader_image;
REQUIRE_NOTHROW(reader_image = reader.LoadImage(1));
REQUIRE(reader_image);
CHECK(reader_image->Image()[0] == SATURATED_PXL_VALUE);
CHECK(reader_image->Image()[1] == 123);
CHECK(reader_image->Image()[5678] == 322);
REQUIRE(reader_image->ImageData().image_collection_efficiency.has_value());
CHECK(reader_image->ImageData().image_collection_efficiency.value() == Catch::Approx(0.91f));
REQUIRE(reader_image->ImageData().az_int_profile.size() == azint_settings.GetBinCount());
CHECK(reader_image->ImageData().az_int_profile[0] == Catch::Approx(51.0f));
CHECK(reader_image->ImageData().az_int_profile[23] == Catch::Approx(51.0f));
}
remove("test_reader_integrated_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_GetRawImage_NXmxLegacy", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test_read_raw_image").ImagesPerTrigger(4).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(2).SetFileWriterFormat(FileWriterFormat::NXmxLegacy).PixelSigned(true)
.IndexingAlgorithm(IndexingAlgorithmEnum::FFT);
x.Compression(CompressionAlgorithm::BSHUF_ZSTD);
std::vector<int16_t> image(x.GetPixelsNum());
for (int i = 0; i < image.size(); i++)
image[i] = static_cast<int16_t>((i * 7 + 33) % UINT16_MAX);
RegisterHDF5Filter();
JFJochBitShuffleCompressor compressor(CompressionAlgorithm::BSHUF_ZSTD);
auto compressed_image = compressor.Compress(image);
{
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(compressed_image, x.GetXPixelsNum(), x.GetYPixelsNum(),
CompressedImageMode::Int16, CompressionAlgorithm::BSHUF_ZSTD);
message.number = 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("test_read_raw_image_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 4);
std::shared_ptr<JFJochReaderRawImage> reader_image;
for (int i = 0; i < 4; i++) {
REQUIRE_NOTHROW(reader_image = reader.GetRawImage(i));
CHECK(reader_image->image.GetMode() == CompressedImageMode::Int16);
CHECK(reader_image->image.GetCompressionAlgorithm() == CompressionAlgorithm::BSHUF_ZSTD);
CHECK(reader_image->image.GetWidth() == x.GetXPixelsNum());
CHECK(reader_image->image.GetHeight() == x.GetYPixelsNum());
CHECK(reader_image->image.GetCompressedSize() == compressed_image.size());
CHECK(reader_image->image.GetCompressed() == reader_image->image_buffer.data());
REQUIRE(reader_image->image_buffer.size() == compressed_image.size());
CHECK(memcmp(reader_image->image_buffer.data(), compressed_image.data(), compressed_image.size()) == 0);
}
}
remove("test_read_raw_image_master.h5");
remove("test_read_raw_image_data_000001.h5");
remove("test_read_raw_image_data_000002.h5");
// No leftover HDF5 objects
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// GetRawImage takes the chunk address under the HDF5 lock, then reads the bytes outside it, so this
// is the one path where several workers are inside the reader at once - which is how rugnux uses it.
// The per-image cases above are all single-threaded and would not notice the file being pulled from
// under a read, nor a cache entry racing its own creation.
TEST_CASE("JFJochReader_GetRawImage_Concurrent", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test_raw_concurrent").ImagesPerTrigger(16).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(4).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true);
x.Compression(CompressionAlgorithm::BSHUF_ZSTD);
std::vector<int16_t> image(x.GetPixelsNum());
for (size_t i = 0; i < image.size(); i++)
image[i] = static_cast<int16_t>((i * 11 + 5) % UINT16_MAX);
RegisterHDF5Filter();
JFJochBitShuffleCompressor compressor(CompressionAlgorithm::BSHUF_ZSTD);
const auto compressed_image = compressor.Compress(image);
{
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(compressed_image, x.GetXPixelsNum(), x.GetYPixelsNum(),
CompressedImageMode::Int16, CompressionAlgorithm::BSHUF_ZSTD);
message.number = 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("test_raw_concurrent_master.h5"));
std::vector<std::future<bool>> workers;
workers.reserve(8);
for (int w = 0; w < 8; w++) {
workers.push_back(std::async(std::launch::async, [&reader, &compressed_image, &x]() {
for (int i = 0; i < x.GetImageNum(); i++) {
auto raw = reader.GetRawImage(i);
if (raw->image_buffer.size() != compressed_image.size())
return false;
if (memcmp(raw->image_buffer.data(), compressed_image.data(),
compressed_image.size()) != 0)
return false;
}
return true;
}));
}
for (auto &worker: workers)
CHECK(worker.get());
}
remove("test_raw_concurrent_master.h5");
for (int f = 1; f <= 4; f++)
remove(("test_raw_concurrent_data_00000" + std::to_string(f) + ".h5").c_str());
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_GetRawImage_VDS", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test_read_raw_image").ImagesPerTrigger(4).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(2).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true)
.IndexingAlgorithm(IndexingAlgorithmEnum::FFT);
x.Compression(CompressionAlgorithm::BSHUF_ZSTD);
std::vector<int16_t> image(x.GetPixelsNum());
for (int i = 0; i < image.size(); i++)
image[i] = static_cast<int16_t>((i * 7 + 33) % UINT16_MAX);
RegisterHDF5Filter();
JFJochBitShuffleCompressor compressor(CompressionAlgorithm::BSHUF_ZSTD);
auto compressed_image = compressor.Compress(image);
{
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(compressed_image, x.GetXPixelsNum(), x.GetYPixelsNum(),
CompressedImageMode::Int16, CompressionAlgorithm::BSHUF_ZSTD);
message.number = 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("test_read_raw_image_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 4);
std::shared_ptr<JFJochReaderRawImage> reader_image;
for (int i = 0; i < 4; i++) {
REQUIRE_NOTHROW(reader_image = reader.GetRawImage(i));
CHECK(reader_image->image.GetMode() == CompressedImageMode::Int16);
CHECK(reader_image->image.GetCompressionAlgorithm() == CompressionAlgorithm::BSHUF_ZSTD);
CHECK(reader_image->image.GetWidth() == x.GetXPixelsNum());
CHECK(reader_image->image.GetHeight() == x.GetYPixelsNum());
CHECK(reader_image->image.GetCompressedSize() == compressed_image.size());
CHECK(reader_image->image.GetCompressed() == reader_image->image_buffer.data());
REQUIRE(reader_image->image_buffer.size() == compressed_image.size());
CHECK(memcmp(reader_image->image_buffer.data(), compressed_image.data(), compressed_image.size()) == 0);
}
}
remove("test_read_raw_image_master.h5");
remove("test_read_raw_image_data_000001.h5");
remove("test_read_raw_image_data_000002.h5");
// No leftover HDF5 objects
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_GetRawImage_Integrated", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test_read_raw_image").ImagesPerTrigger(4).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(2).SetFileWriterFormat(FileWriterFormat::NXmxIntegrated).PixelSigned(true)
.IndexingAlgorithm(IndexingAlgorithmEnum::FFT);
x.Compression(CompressionAlgorithm::BSHUF_ZSTD);
std::vector<int16_t> image(x.GetPixelsNum());
for (int i = 0; i < image.size(); i++)
image[i] = static_cast<int16_t>((i * 7 + 33) % UINT16_MAX);
RegisterHDF5Filter();
JFJochBitShuffleCompressor compressor(CompressionAlgorithm::BSHUF_ZSTD);
auto compressed_image = compressor.Compress(image);
{
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(compressed_image, x.GetXPixelsNum(), x.GetYPixelsNum(),
CompressedImageMode::Int16, CompressionAlgorithm::BSHUF_ZSTD);
message.number = 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("test_read_raw_image_master.h5"));
auto dataset = reader.GetDataset();
CHECK(dataset->experiment.GetImageNum() == 4);
std::shared_ptr<JFJochReaderRawImage> reader_image;
for (int i = 0; i < 4; i++) {
REQUIRE_NOTHROW(reader_image = reader.GetRawImage(i));
CHECK(reader_image->image.GetMode() == CompressedImageMode::Int16);
CHECK(reader_image->image.GetCompressionAlgorithm() == CompressionAlgorithm::BSHUF_ZSTD);
CHECK(reader_image->image.GetWidth() == x.GetXPixelsNum());
CHECK(reader_image->image.GetHeight() == x.GetYPixelsNum());
CHECK(reader_image->image.GetCompressedSize() == compressed_image.size());
CHECK(reader_image->image.GetCompressed() == reader_image->image_buffer.data());
REQUIRE(reader_image->image_buffer.size() == compressed_image.size());
CHECK(memcmp(reader_image->image_buffer.data(), compressed_image.data(), compressed_image.size()) == 0);
}
}
remove("test_read_raw_image_master.h5");
// No leftover HDF5 objects
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_HDF5DataSource_Integrated", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("source_integrated").ImagesPerTrigger(5).OverwriteExistingFiles(true);
x.BitDepthImage(16).SetFileWriterFormat(FileWriterFormat::NXmxIntegrated).PixelSigned(true);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
std::vector<int16_t> image(x.GetPixelsNum(), 17);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter writer(start_message);
for (int i = 0; i < x.GetImageNum(); i++) {
image[5678] = static_cast<int16_t>(100 + i);
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = i;
REQUIRE_NOTHROW(writer.WriteHDF5(message));
}
EndMessage end_message;
end_message.max_image_number = x.GetImageNum();
writer.WriteHDF5(end_message);
writer.Finalize();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("source_integrated_master.h5"));
auto source = reader.GetHDF5DataSource(1, 3);
REQUIRE(source.size() == 1);
CHECK(source[0].filename == "source_integrated_master.h5");
CHECK(source[0].dataset == "/entry/data/data");
CHECK(source[0].source_first_image == 1);
CHECK(source[0].virtual_first_image == 0);
CHECK(source[0].image_count == 3);
}
remove("source_integrated_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_HDF5DataSource_VDS", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("source_vds_mapping").ImagesPerTrigger(5).ImagesPerFile(2).OverwriteExistingFiles(true);
x.BitDepthImage(16).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
std::vector<int16_t> image(x.GetPixelsNum(), 21);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter writer(start_message);
for (int i = 0; i < x.GetImageNum(); i++) {
image[5678] = static_cast<int16_t>(200 + i);
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = i;
REQUIRE_NOTHROW(writer.WriteHDF5(message));
}
EndMessage end_message;
end_message.max_image_number = x.GetImageNum();
writer.WriteHDF5(end_message);
writer.Finalize();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("source_vds_mapping_master.h5"));
// Range crosses file boundary:
// global images 1,2,3 map to:
// data_000001 image 1
// data_000002 images 0,1
auto source = reader.GetHDF5DataSource(1, 3);
REQUIRE(source.size() == 2);
CHECK(source[0].filename == "source_vds_mapping_data_000001.h5");
CHECK(source[0].dataset == "/entry/data/data");
CHECK(source[0].source_first_image == 1);
CHECK(source[0].virtual_first_image == 0);
CHECK(source[0].image_count == 1);
CHECK(source[1].filename == "source_vds_mapping_data_000002.h5");
CHECK(source[1].dataset == "/entry/data/data");
CHECK(source[1].source_first_image == 0);
CHECK(source[1].virtual_first_image == 1);
CHECK(source[1].image_count == 2);
}
remove("source_vds_mapping_master.h5");
remove("source_vds_mapping_data_000001.h5");
remove("source_vds_mapping_data_000002.h5");
remove("source_vds_mapping_data_000003.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_HDF5DataSource_Legacy", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("source_legacy_mapping").ImagesPerTrigger(5).ImagesPerFile(2).OverwriteExistingFiles(true);
x.BitDepthImage(16).SetFileWriterFormat(FileWriterFormat::NXmxLegacy).PixelSigned(true);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
std::vector<int16_t> image(x.GetPixelsNum(), 31);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter writer(start_message);
for (int i = 0; i < x.GetImageNum(); i++) {
image[5678] = static_cast<int16_t>(300 + i);
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = i;
REQUIRE_NOTHROW(writer.WriteHDF5(message));
}
EndMessage end_message;
end_message.max_image_number = x.GetImageNum();
writer.WriteHDF5(end_message);
writer.Finalize();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("source_legacy_mapping_master.h5"));
auto source = reader.GetHDF5DataSource(1, 3);
REQUIRE(source.size() == 2);
CHECK(source[0].filename == "source_legacy_mapping_data_000001.h5");
CHECK(source[0].dataset == "/entry/data/data");
CHECK(source[0].source_first_image == 1);
CHECK(source[0].virtual_first_image == 0);
CHECK(source[0].image_count == 1);
CHECK(source[1].filename == "source_legacy_mapping_data_000002.h5");
CHECK(source[1].dataset == "/entry/data/data");
CHECK(source[1].source_first_image == 0);
CHECK(source[1].virtual_first_image == 1);
CHECK(source[1].image_count == 2);
}
remove("source_legacy_mapping_master.h5");
remove("source_legacy_mapping_data_000001.h5");
remove("source_legacy_mapping_data_000002.h5");
remove("source_legacy_mapping_data_000003.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
TEST_CASE("JFJochReader_ProcessingHDF5_FromVDS_MapsToDataFiles", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("proc_source_vds").ImagesPerTrigger(5).ImagesPerFile(2).OverwriteExistingFiles(true);
x.BitDepthImage(16).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
std::vector<int16_t> image(x.GetPixelsNum(), 51);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter writer(start_message);
for (int i = 0; i < x.GetImageNum(); i++) {
image[5678] = static_cast<int16_t>(500 + i);
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = i;
REQUIRE_NOTHROW(writer.WriteHDF5(message));
}
EndMessage end_message;
end_message.max_image_number = x.GetImageNum();
writer.WriteHDF5(end_message);
writer.Finalize();
}
std::vector<HDF5DataSourceMessage> source_data;
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("proc_source_vds_master.h5"));
source_data = reader.GetHDF5DataSource(1, 3);
REQUIRE(source_data.size() == 2);
CHECK(source_data[0].filename == "proc_source_vds_data_000001.h5");
CHECK(source_data[1].filename == "proc_source_vds_data_000002.h5");
}
{
DiffractionExperiment proc_x = x;
proc_x.FilePrefix("proc_from_vds")
.ImagesPerTrigger(3)
.SetFileWriterFormat(FileWriterFormat::NXmxIntegrated)
.OverwriteExistingFiles(true);
StartMessage start_message;
proc_x.FillMessage(start_message);
start_message.number_of_images = 3;
start_message.images_per_file = 3;
start_message.write_images = false;
start_message.write_master_file = true;
start_message.hdf5_source_data = source_data;
FileWriter writer(start_message);
for (int i = 0; i < 3; i++) {
DataMessage message{};
message.number = i;
message.original_number = i + 1;
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.spot_count = 200 + i;
REQUIRE_NOTHROW(writer.WriteHDF5(message));
}
EndMessage end_message;
end_message.max_image_number = 3;
writer.WriteHDF5(end_message);
writer.Finalize();
}
{
HDF5ReadOnlyFile file("proc_from_vds_master.h5");
HDF5DataSet data(file, "/entry/data/data");
HDF5Dcpl dcpl(data);
REQUIRE(dcpl.GetLayout() == HDF5DataSetLayout::VIRTUAL);
auto mappings = dcpl.GetVirtualMappings();
REQUIRE(mappings.size() == 2);
CHECK(mappings[0].filename == "proc_source_vds_data_000001.h5");
CHECK(mappings[1].filename == "proc_source_vds_data_000002.h5");
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("proc_from_vds_master.h5"));
auto img0 = reader.LoadImage(0);
REQUIRE(img0);
CHECK(img0->Image()[5678] == 501);
auto img2 = reader.LoadImage(2);
REQUIRE(img2);
CHECK(img2->Image()[5678] == 503);
}
remove("proc_source_vds_master.h5");
remove("proc_source_vds_data_000001.h5");
remove("proc_source_vds_data_000002.h5");
remove("proc_source_vds_data_000003.h5");
remove("proc_from_vds_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// rugnux --stride N processes every Nth image, so image i of the _process.h5 IS source image
// start + i*N. The mapping used to be built without the stride, linking the first N images instead,
// which put each frame's picture next to a different frame's analysis. Read the pixels back through
// the written VDS rather than only inspecting the mapping: that is what a user opening the file sees.
TEST_CASE("JFJochReader_ProcessingHDF5_Strided_LinksTheProcessedImages", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("proc_stride_src").ImagesPerTrigger(10).ImagesPerFile(2).OverwriteExistingFiles(true);
x.BitDepthImage(16).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
std::vector<int16_t> image(x.GetPixelsNum(), 51);
RegisterHDF5Filter();
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter writer(start_message);
for (int i = 0; i < x.GetImageNum(); i++) {
image[5678] = static_cast<int16_t>(500 + i); // per-image tag: which source frame is this?
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = i;
REQUIRE_NOTHROW(writer.WriteHDF5(message));
}
EndMessage end_message;
end_message.max_image_number = x.GetImageNum();
writer.WriteHDF5(end_message);
writer.Finalize();
}
// start 1, 3 images, stride 3 -> source images 1, 4, 7. Two images per data file, so those sit in
// data_000001 (holds 0,1), data_000003 (holds 4,5) and data_000004 (holds 6,7).
std::vector<HDF5DataSourceMessage> source_data;
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("proc_stride_src_master.h5"));
source_data = reader.GetHDF5DataSource(1, 3, 3);
// Non-adjacent source images cannot be merged into one run, so one mapping per image.
REQUIRE(source_data.size() == 3);
CHECK(source_data[0].filename == "proc_stride_src_data_000001.h5");
CHECK(source_data[0].source_first_image == 1);
CHECK(source_data[0].virtual_first_image == 0);
CHECK(source_data[0].image_count == 1);
CHECK(source_data[1].filename == "proc_stride_src_data_000003.h5");
CHECK(source_data[1].source_first_image == 0);
CHECK(source_data[1].virtual_first_image == 1);
CHECK(source_data[1].image_count == 1);
CHECK(source_data[2].filename == "proc_stride_src_data_000004.h5");
CHECK(source_data[2].source_first_image == 1);
CHECK(source_data[2].virtual_first_image == 2);
CHECK(source_data[2].image_count == 1);
}
{
DiffractionExperiment proc_x = x;
proc_x.FilePrefix("proc_from_stride")
.ImagesPerTrigger(3)
.SetFileWriterFormat(FileWriterFormat::NXmxIntegrated)
.OverwriteExistingFiles(true);
StartMessage start_message;
proc_x.FillMessage(start_message);
start_message.number_of_images = 3;
start_message.images_per_file = 3;
start_message.write_images = false;
start_message.write_master_file = true;
start_message.hdf5_source_data = source_data;
FileWriter writer(start_message);
for (int i = 0; i < 3; i++) {
DataMessage message{};
message.number = i;
message.original_number = 1 + i * 3;
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.spot_count = 200 + i;
REQUIRE_NOTHROW(writer.WriteHDF5(message));
}
EndMessage end_message;
end_message.max_image_number = 3;
writer.WriteHDF5(end_message);
writer.Finalize();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("proc_from_stride_master.h5"));
// The whole point: frame i of the process file carries source frame 1 + 3i, not 1 + i.
auto img0 = reader.LoadImage(0);
REQUIRE(img0);
CHECK(img0->Image()[5678] == 501);
auto img1 = reader.LoadImage(1);
REQUIRE(img1);
CHECK(img1->Image()[5678] == 504);
auto img2 = reader.LoadImage(2);
REQUIRE(img2);
CHECK(img2->Image()[5678] == 507);
}
remove("proc_stride_src_master.h5");
for (int i = 1; i <= 5; i++)
remove(("proc_stride_src_data_00000" + std::to_string(i) + ".h5").c_str());
remove("proc_from_stride_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// One reflection with a distinct value in EVERY field that is meant to survive a write/read cycle,
// keyed on the image and the reflection index so no two are alike and a field read back from the
// wrong place shows up. Fractional image_number on purpose: a 3D-integrated reflection has one, and
// it is the field the offline --scale path silently lost when it was not read back.
static Reflection MakeTestReflection(int i, int j) {
const auto f = static_cast<float>(i * 10 + j);
return Reflection{
.h = 10 + i + 30 * j,
.k = 20 + j,
.l = 30 + j,
.image_number = static_cast<float>(i) + 0.25f * static_cast<float>(j + 1),
.delta_phi_deg = 0.1f + 0.01f * f,
.predicted_x = 100.0f + f,
.predicted_y = 200.0f + f,
.observed_x = 100.5f + f,
.observed_y = 200.5f + f,
.d = 1.5f + 0.1f * f,
.I = 1000.0f + f,
.bkg = 10.0f + f,
.sigma = 2.0f + 0.5f * f,
.rlp = 1.0f + 0.125f * f,
.partiality = 0.5f + 0.01f * f,
.zeta = 0.01f + 0.001f * f,
.image_scale_corr = 1.0f + 0.25f * f
};
}
// Every field of the round trip, against the reflection that was written. Deliberately NOT checked,
// because they are not part of it: dist_ewald and observed are prediction/integration scratch that
// is never written, and on_ice_ring is recomputed from the resolution by whoever scales (see the
// ice-ring handling in Rugnux and in the --scale path).
static void CheckReflectionRoundTrip(const Reflection &got, int i, int j) {
INFO("image " << i << " reflection " << j);
const Reflection want = MakeTestReflection(i, j);
CHECK(got.h == want.h);
CHECK(got.k == want.k);
CHECK(got.l == want.l);
CHECK(got.image_number == Catch::Approx(want.image_number));
CHECK(got.delta_phi_deg == Catch::Approx(want.delta_phi_deg));
CHECK(got.predicted_x == Catch::Approx(want.predicted_x));
CHECK(got.predicted_y == Catch::Approx(want.predicted_y));
CHECK(got.observed_x == Catch::Approx(want.observed_x));
CHECK(got.observed_y == Catch::Approx(want.observed_y));
CHECK(got.d == Catch::Approx(want.d));
CHECK(got.I == Catch::Approx(want.I));
CHECK(got.bkg == Catch::Approx(want.bkg));
CHECK(got.sigma == Catch::Approx(want.sigma));
CHECK(got.rlp == Catch::Approx(want.rlp)); // stored as 1/rlp, inverted again on read
CHECK(got.partiality == Catch::Approx(want.partiality));
CHECK(got.zeta == Catch::Approx(want.zeta));
CHECK(got.image_scale_corr == Catch::Approx(want.image_scale_corr));
}
TEST_CASE("JFJochReader_ReadReflections_VDS", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.FilePrefix("read_reflections_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 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 || i == 3) {
message.integrated_reflections = 2;
message.reflections = {MakeTestReflection(i, 0), MakeTestReflection(i, 1)};
message.mosaicity_deg = i*0.15f;
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();
scan_result.FillEndMessage(end_message);
writer.WriteHDF5(end_message);
writer.Finalize();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("read_reflections_vds_master.h5"));
auto reflections = reader.ReadReflections();
REQUIRE(reflections.size() == 4);
CHECK(reflections[0].reflections.empty());
REQUIRE(reflections[1].reflections.size() == 2);
CheckReflectionRoundTrip(reflections[1].reflections[0], 1, 0);
CheckReflectionRoundTrip(reflections[1].reflections[1], 1, 1);
CHECK(reflections[1].mosaicity_deg == Catch::Approx(0.15f));
CHECK(reflections[1].latt.CalcVolume() == Catch::Approx(100*50*30));
CHECK(reflections[2].reflections.empty());
REQUIRE(reflections[3].reflections.size() == 2);
CheckReflectionRoundTrip(reflections[3].reflections[0], 3, 0);
CheckReflectionRoundTrip(reflections[3].reflections[1], 3, 1);
CHECK(reflections[3].mosaicity_deg == Catch::Approx(0.45f));
CHECK(reflections[3].latt.Vec0().x == Catch::Approx(100.0f));
CHECK(reflections[3].latt.Vec1().y == Catch::Approx(50.0f));
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("read_reflections_vds_master.h5"));
auto reflections = reader.ReadReflections(1, 3);
REQUIRE(reflections.size() == 3);
REQUIRE(reflections[0].reflections.size() == 2); // original image 1
CheckReflectionRoundTrip(reflections[0].reflections[0], 1, 0);
CHECK(reflections[1].reflections.empty()); // original image 2
REQUIRE(reflections[2].reflections.size() == 2); // original image 3
CheckReflectionRoundTrip(reflections[2].reflections[0], 3, 0);
}
remove("read_reflections_vds_master.h5");
remove("read_reflections_vds_data_000001.h5");
remove("read_reflections_vds_data_000002.h5");
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);
}
// 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<int32_t, 9> M = {1, 1, 0,
0, 1, 1,
1, 0, 1};
std::vector<int16_t> 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<SpotToSave> MakeTestSpots(int i) {
return {
SpotToSave{
.x = 1, .y = 2, .intensity = 376,
.h = 11, .k = -3, .l = -5,
.dist_ewald_sphere = 0.1234f,
.ice_ring = true,
.indexed = true
},
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);
}
TEST_CASE("JFJochReader_Snapshots", "[HDF5][Full]") {
RegisterHDF5Filter();
DiffractionExperiment x(DetJF(1));
x.FilePrefix("test_snap").ImagesPerTrigger(4).OverwriteExistingFiles(true);
x.BitDepthImage(16).ImagesPerFile(1).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true);
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
// 1. Original dataset: distinct pixels + "original" MX metadata (not indexed, bkg = 10 + i).
std::vector<int16_t> image(x.GetPixelsNum());
{
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
ScanResultGenerator generator(x);
for (int i = 0; i < 4; i++) {
image[5678] = 100 + i;
DataMessage message{};
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = i;
message.indexing_result = false;
message.bkg_estimate = 10.0 + i;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
}
EndMessage end_message;
end_message.max_image_number = 4;
generator.FillEndMessage(end_message);
file_set.WriteHDF5(end_message);
file_set.Finalize();
}
// 2. A reprocessing result over the same 4 images: integrated master, "reprocessed" MX
// metadata (all indexed, bkg = 99 + i) and deliberately wrong pixels that must NOT surface
// (snapshot pixels still come from the original image source).
{
DiffractionExperiment px(x);
px.FilePrefix("test_snap_proc").SetFileWriterFormat(FileWriterFormat::NXmxIntegrated);
StartMessage start_message;
px.FillMessage(start_message);
FileWriter file_set(start_message);
ScanResultGenerator generator(px);
std::vector<int16_t> proc_image(x.GetPixelsNum(), 7);
for (int i = 0; i < 4; i++) {
DataMessage message{};
message.image = CompressedImage(proc_image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = i;
message.indexing_result = true;
message.bkg_estimate = 99.0 + i;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
generator.Add(message);
}
EndMessage end_message;
end_message.max_image_number = 4;
generator.FillEndMessage(end_message);
file_set.WriteHDF5(end_message);
file_set.Finalize();
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("test_snap_master.h5"));
REQUIRE(reader.GetNumberOfImages() == 4);
CHECK(reader.ActiveSnapshot() == "Original");
CHECK(reader.SnapshotNames() == std::vector<std::string>{"Original"});
// Dataset-level plot arrays come from the active (original) metadata source.
REQUIRE(reader.GetDataset()->bkg_estimate.size() == 4);
CHECK(reader.GetDataset()->bkg_estimate[2] == Catch::Approx(12.0));
auto orig0 = reader.LoadImage(0);
REQUIRE(orig0);
CHECK(orig0->Image()[5678] == 100);
REQUIRE(orig0->ImageData().indexing_result.has_value());
CHECK(orig0->ImageData().indexing_result.value() == false);
CHECK(orig0->ImageData().bkg_estimate.value() == Catch::Approx(10.0));
// Register the reprocessing result as a second metadata source over the same images.
REQUIRE_NOTHROW(reader.RegisterSnapshot("Reprocess", "test_snap_proc_master.h5"));
{
auto names = reader.SnapshotNames();
CHECK(std::find(names.begin(), names.end(), "Original") != names.end());
CHECK(std::find(names.begin(), names.end(), "Reprocess") != names.end());
}
REQUIRE_NOTHROW(reader.SetActiveSnapshot("Reprocess"));
CHECK(reader.ActiveSnapshot() == "Reprocess");
// Plots now come from the reprocessing master.
CHECK(reader.GetDataset()->bkg_estimate[2] == Catch::Approx(101.0));
auto repro0 = reader.LoadImage(0);
REQUIRE(repro0);
// Pixels still from the original image source, not the 7's stored in the process file.
CHECK(repro0->Image()[5678] == 100);
// Metadata from the reprocessing snapshot.
CHECK(repro0->ImageData().indexing_result.value() == true);
CHECK(repro0->ImageData().bkg_estimate.value() == Catch::Approx(99.0));
// Switch back to the original metadata.
REQUIRE_NOTHROW(reader.SetActiveSnapshot("Original"));
auto orig0b = reader.LoadImage(0);
REQUIRE(orig0b);
CHECK(orig0b->Image()[5678] == 100);
CHECK(orig0b->ImageData().indexing_result.value() == false);
REQUIRE_THROWS(reader.SetActiveSnapshot("Nonexistent"));
}
remove("test_snap_master.h5");
remove("test_snap_data_000001.h5");
remove("test_snap_data_000002.h5");
remove("test_snap_data_000003.h5");
remove("test_snap_data_000004.h5");
remove("test_snap_proc_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// The chain may be sent in the END message or left to the writer to build. Both must produce the
// same file - the sent one is written verbatim, which is what will later allow measured positions to
// be reported, and the built one is what a producer that does not send it gets.
TEST_CASE("JFJochReader_TransformationChain_SentAndBuilt", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(5).OverwriteExistingFiles(true);
x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
.IncidentEnergy_keV(WVL_1A_IN_KEV).PixelSigned(false).BitDepthImage(16)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
x.Goniometer(GoniometerAxis("omega", 95, 0.1f, Coord(0,-1,0), {}));
x.Smargon(SmargonPosition{.phi_deg = -7.25f, .chi_deg = 12.5f});
RegisterHDF5Filter();
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
const auto write = [&](const std::string &prefix, bool send_chain) {
DiffractionExperiment local = x;
local.FilePrefix(prefix);
StartMessage start_message;
local.FillMessage(start_message);
FileWriter file_set(start_message);
DataMessage message{};
for (int i = 0; i < local.GetImageNum(); i++) {
message.image = CompressedImage(image, local.GetXPixelsNum(), local.GetYPixelsNum());
message.number = i;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
}
EndMessage end_message;
end_message.max_image_number = local.GetImageNum();
if (send_chain)
end_message.transformations = local.BuildTransformationChain(local.GetImageNum());
file_set.WriteHDF5(end_message);
file_set.Finalize();
};
write("test_chain_built", false);
write("test_chain_sent", true);
const auto read = [](const std::string &prefix) {
JFJochHDF5Reader reader;
reader.ReadFile(prefix + "_master.h5");
return reader.GetDataset()->experiment;
};
const auto built = read("test_chain_built");
const auto sent = read("test_chain_sent");
REQUIRE(built.GetGoniometer().has_value());
REQUIRE(sent.GetGoniometer().has_value());
CHECK(sent.GetGoniometer()->GetName() == built.GetGoniometer()->GetName());
CHECK(sent.GetGoniometer()->GetStart_deg()
== Catch::Approx(built.GetGoniometer()->GetStart_deg()).margin(1e-3));
CHECK(sent.GetGoniometer()->GetIncrement_deg()
== Catch::Approx(built.GetGoniometer()->GetIncrement_deg()).margin(1e-4));
// chi/phi survive both routes, which they did not before they became ordinary axes.
REQUIRE(built.GetDatasetSettings().GetSmargonPosition().has_value());
REQUIRE(sent.GetDatasetSettings().GetSmargonPosition().has_value());
CHECK(sent.GetDatasetSettings().GetSmargonPosition()->chi_deg
== Catch::Approx(12.5f).margin(1e-3));
CHECK(sent.GetDatasetSettings().GetSmargonPosition()->phi_deg
== Catch::Approx(-7.25f).margin(1e-3));
// Compared on the files, not through the reader: the reader reads neither AXISNAME_end nor the
// rotation width, so it cannot see the two routes diverge - and it did, until the writer started
// deriving them for a chain it was handed.
{
HDF5ReadOnlyFile built_file("test_chain_built_master.h5");
HDF5ReadOnlyFile sent_file("test_chain_sent_master.h5");
CHECK(built_file.FindLeafs("/entry/sample/transformations")
== sent_file.FindLeafs("/entry/sample/transformations"));
CHECK(sent_file.ReadVector<double>("/entry/sample/transformations/omega_end").size()
== static_cast<size_t>(x.GetImageNum()));
CHECK(sent_file.ReadVector<double>("/entry/sample/transformations/omega_range_average").at(0)
== Catch::Approx(0.1).margin(1e-4));
}
remove("test_chain_built_master.h5");
remove("test_chain_sent_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// Recovering the sample axes from a written file, across the configurations the writer produces.
// The reader searches every leaf of /entry/sample/transformations, so on the way it meets the
// writer's own AXISNAME_end and rotation-width datasets, which are not axes and carry no
// transformation_type. It used to throw on them: a master whose axis did not turn never stopped the
// search early, walked into omega_end and could not be opened at all - which took out rugnux's own
// output for a grid scan.
TEST_CASE("JFJochReader_AxisRecovery", "[HDF5][Full]") {
RegisterHDF5Filter();
const auto round_trip = [](DiffractionExperiment x, const std::string &prefix) {
x.FilePrefix(prefix).OverwriteExistingFiles(true)
.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
.IncidentEnergy_keV(WVL_1A_IN_KEV).PixelSigned(false).BitDepthImage(16)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
StartMessage start_message;
x.FillMessage(start_message);
{
FileWriter file_set(start_message);
DataMessage message{};
for (int i = 0; i < x.GetImageNum(); i++) {
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = 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(prefix + "_master.h5"));
return reader.GetDataset()->experiment;
};
const auto cleanup = [](const std::string &prefix) {
remove((prefix + "_master.h5").c_str());
remove((prefix + "_data_000001.h5").c_str());
};
SECTION("a sweep") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(5).Goniometer(GoniometerAxis("omega", 95, 0.1f, Coord(0,-1,0), {}));
const auto out = round_trip(x, "test_ax_sweep");
REQUIRE(out.GetGoniometer().has_value());
CHECK(out.GetGoniometer()->GetName() == "omega");
CHECK(out.GetGoniometer()->IsScanning());
CHECK(out.GetGoniometer()->GetStart_deg() == Catch::Approx(95).margin(1e-3));
CHECK(out.GetGoniometer()->GetIncrement_deg() == Catch::Approx(0.1).margin(1e-4));
CHECK(out.GetGoniometer()->GetAxis() == Coord(0,-1,0));
cleanup("test_ax_sweep");
}
SECTION("a sweep about an axis that is not called omega") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(5).Goniometer(GoniometerAxis("kappa", 10, 0.5f, Coord(0,-1,0), {}));
const auto out = round_trip(x, "test_ax_kappa");
REQUIRE(out.GetGoniometer().has_value());
CHECK(out.GetGoniometer()->GetName() == "kappa");
CHECK(out.GetGoniometer()->IsScanning());
CHECK(out.GetGoniometer()->GetIncrement_deg() == Catch::Approx(0.5).margin(1e-4));
cleanup("test_ax_kappa");
}
SECTION("a spindle that does not turn") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(5).Goniometer(GoniometerAxis("omega", 12.5f, 0.0f, Coord(0,-1,0), {}));
const auto out = round_trip(x, "test_ax_still");
REQUIRE(out.GetGoniometer().has_value());
CHECK(out.GetGoniometer()->GetName() == "omega");
CHECK(!out.GetGoniometer()->IsScanning());
CHECK(out.GetGoniometer()->GetStart_deg() == Catch::Approx(12.5).margin(1e-3));
cleanup("test_ax_still");
}
SECTION("a grid scan, which sits on a spindle that does not turn") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(6).GridScan(GridScanSettings(3, 10.0f, 20.0f, false, false).ImageNum(6));
const auto out = round_trip(x, "test_ax_grid");
REQUIRE(out.GetGridScan().has_value());
CHECK(out.GetGridScan()->GetNFast() == 3);
CHECK(out.GetGridScan()->GetGridStepX_um() == Catch::Approx(10.0).margin(1e-3));
CHECK(out.GetGridScan()->GetGridStepY_um() == Catch::Approx(20.0).margin(1e-3));
// NXmx cannot say "no rotation", so the writer records the spindle standing still.
REQUIRE(out.GetGoniometer().has_value());
CHECK(!out.GetGoniometer()->IsScanning());
cleanup("test_ax_grid");
}
// A grid scan is taken at a stationary spindle, and the angle it stood at is what relates one
// grid to another taken elsewhere on the circle. It is stated by sending the axis with step 0;
// send nothing and the spindle is recorded at 0, which says only that nobody told us.
SECTION("a grid scan at a stationary head position") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(6).GridScan(GridScanSettings(3, 10.0f, 20.0f, false, false).ImageNum(6))
.Goniometer(GoniometerAxis("omega", 90.0f, 0.0f, Coord(-1,0,0), {}));
const auto out = round_trip(x, "test_ax_gridstill");
REQUIRE(out.GetGridScan().has_value());
CHECK(out.GetGridScan()->GetNFast() == 3);
REQUIRE(out.GetGoniometer().has_value());
CHECK(out.GetGoniometer()->GetName() == "omega");
CHECK(!out.GetGoniometer()->IsScanning());
CHECK(out.GetGoniometer()->GetStart_deg() == Catch::Approx(90).margin(1e-3));
CHECK(out.GetGoniometer()->GetAxis() == Coord(-1,0,0));
cleanup("test_ax_gridstill");
}
SECTION("a grid scan under a turning spindle") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(6).GridScan(GridScanSettings(3, 10.0f, 20.0f, false, false).ImageNum(6))
.Goniometer(GoniometerAxis("omega", 0, 0.2f, Coord(0,-1,0), {}));
const auto out = round_trip(x, "test_ax_gridsweep");
REQUIRE(out.GetGridScan().has_value());
CHECK(out.GetGridScan()->GetNFast() == 3);
REQUIRE(out.GetGoniometer().has_value());
CHECK(out.GetGoniometer()->IsScanning());
CHECK(out.GetGoniometer()->GetIncrement_deg() == Catch::Approx(0.2).margin(1e-4));
cleanup("test_ax_gridsweep");
}
SECTION("a sweep with the head at a Smargon position") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(5).Goniometer(GoniometerAxis("omega", 95, 0.1f, Coord(0,-1,0), {}))
.Smargon(SmargonPosition{.phi_deg = -7.25f, .chi_deg = 12.5f});
const auto out = round_trip(x, "test_ax_smargon");
REQUIRE(out.GetGoniometer().has_value());
CHECK(out.GetGoniometer()->GetName() == "omega");
CHECK(out.GetGoniometer()->IsScanning());
REQUIRE(out.GetDatasetSettings().GetSmargonPosition().has_value());
CHECK(out.GetDatasetSettings().GetSmargonPosition()->chi_deg == Catch::Approx(12.5).margin(1e-3));
CHECK(out.GetDatasetSettings().GetSmargonPosition()->phi_deg == Catch::Approx(-7.25).margin(1e-3));
cleanup("test_ax_smargon");
}
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// saturation_value is written inclusive and used exclusive, so a read has to add the count back. It
// did not, and the value fell by one on every write-read-write cycle - unbounded, and compounding
// whenever a _process.h5 was reprocessed. Nothing caught it: no test asserted the read-back limit.
TEST_CASE("JFJochReader_SaturationSurvivesRoundTrip", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(2).OverwriteExistingFiles(true).FilePrefix("test_satrt");
x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
.IncidentEnergy_keV(WVL_1A_IN_KEV).PixelSigned(false).BitDepthImage(16)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
const int64_t original_limit = x.GetSaturationLimit();
RegisterHDF5Filter();
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
// The second pass writes metadata only: reading pins the experiment to signed 32-bit, the
// container images are handed out in, so feeding it the uint16 frames again would - rightly -
// be refused by the writer's pixel-format check.
const auto write = [&](const DiffractionExperiment &src, const std::string &prefix,
bool with_images) {
DiffractionExperiment local = src;
local.FilePrefix(prefix).OverwriteExistingFiles(true);
StartMessage start_message;
local.FillMessage(start_message);
FileWriter file_set(start_message);
DataMessage message{};
if (with_images) {
for (int i = 0; i < 2; i++) {
message.image = CompressedImage(image, local.GetXPixelsNum(), local.GetYPixelsNum());
message.number = i;
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
}
}
EndMessage end_message;
end_message.max_image_number = 2;
file_set.WriteHDF5(end_message);
file_set.Finalize();
return start_message.saturation_value;
};
const auto read = [](const std::string &prefix) {
JFJochHDF5Reader reader;
reader.ReadFile(prefix + "_master.h5");
return reader.GetDataset()->experiment;
};
const int64_t declared_once = write(x, "test_satrt", true);
CHECK(declared_once == SaturationValueFromLimit(original_limit));
const auto once = read("test_satrt");
CHECK(once.GetSaturationLimit() == original_limit);
// The cycle that used to lose a count: read a file, write what was read, read it again.
const int64_t declared_twice = write(once, "test_satrt2", false);
CHECK(declared_twice == declared_once);
CHECK(read("test_satrt2").GetSaturationLimit() == original_limit);
remove("test_satrt_master.h5");
remove("test_satrt_data_000001.h5");
remove("test_satrt2_master.h5");
remove("test_satrt2_data_000001.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// A head position is not a sweep, and it is not the spindle either. Both properties are carried by
// the file itself - the axis length says how many images there are, the equipment_component tag says
// what the axis is - so both are checked here on the file, not through the reader: the reader alone
// cannot see a shape or an attribute it never looks at.
TEST_CASE("JFJochReader_Smargon_StillIsNotOneImage", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(5).OverwriteExistingFiles(true).FilePrefix("test_smargon");
x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
.IncidentEnergy_keV(WVL_1A_IN_KEV).PixelSigned(false).BitDepthImage(16)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
x.Smargon(SmargonPosition{.phi_deg = -7.25f, .chi_deg = 12.5f});
REQUIRE(!x.GetGoniometer().has_value());
RegisterHDF5Filter();
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
DataMessage message{};
for (int i = 0; i < x.GetImageNum(); i++) {
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = 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();
{
HDF5ReadOnlyFile master("test_smargon_master.h5");
// One entry per image. A reader takes the image count from the innermost axis of the sample
// chain when no axis varies; as scalars these read back as a single image.
CHECK(master.GetDimension("/entry/sample/transformations/chi")
== std::vector<hsize_t>{static_cast<hsize_t>(x.GetImageNum())});
CHECK(master.GetDimension("/entry/sample/transformations/phi")
== std::vector<hsize_t>{static_cast<hsize_t>(x.GetImageNum())});
CHECK(master.ReadVector<double>("/entry/sample/transformations/phi")
== std::vector<double>(x.GetImageNum(), -7.25));
// Tagged, so neither is mistaken for the spindle - and so a phi from anywhere else is not
// mistaken for a head position.
HDF5DataSet chi(master, "/entry/sample/transformations/chi");
HDF5DataSet phi(master, "/entry/sample/transformations/phi");
REQUIRE(chi.AttrExists("equipment_component"));
REQUIRE(phi.AttrExists("equipment_component"));
CHECK(chi.ReadAttrStr("equipment_component") == "smargon");
CHECK(phi.ReadAttrStr("equipment_component") == "smargon");
}
const auto read = [](const std::string &prefix) {
JFJochHDF5Reader reader;
reader.ReadFile(prefix + "_master.h5");
return reader.GetDataset()->experiment;
};
const auto read_back = read("test_smargon");
// chi is the alphabetically first stationary axis in the file; it must not become the spindle.
CHECK(!read_back.GetGoniometer().has_value());
REQUIRE(read_back.GetDatasetSettings().GetSmargonPosition().has_value());
CHECK(read_back.GetDatasetSettings().GetSmargonPosition()->chi_deg == Catch::Approx(12.5f).margin(1e-3));
CHECK(read_back.GetDatasetSettings().GetSmargonPosition()->phi_deg == Catch::Approx(-7.25f).margin(1e-3));
remove("test_smargon_master.h5");
remove("test_smargon_data_000001.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// phi is an ordinary spindle name in MX. A file whose rotation axis is called phi carries no
// equipment_component, so it stays the spindle and no head position is invented from it - which also
// means the file can be written back out, instead of colliding on a second dataset called phi.
TEST_CASE("JFJochReader_Goniometer_NamedPhiIsNotSmargon", "[HDF5][Full]") {
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(5).OverwriteExistingFiles(true).FilePrefix("test_phispindle");
x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
.IncidentEnergy_keV(WVL_1A_IN_KEV).PixelSigned(false).BitDepthImage(16)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
x.Goniometer(GoniometerAxis("phi", 30, 0.2f, Coord(0,-1,0), {}));
RegisterHDF5Filter();
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
DataMessage message{};
for (int i = 0; i < x.GetImageNum(); i++) {
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = 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();
const auto read = [](const std::string &prefix) {
JFJochHDF5Reader reader;
reader.ReadFile(prefix + "_master.h5");
return reader.GetDataset()->experiment;
};
const auto read_back = read("test_phispindle");
REQUIRE(read_back.GetGoniometer().has_value());
CHECK(read_back.GetGoniometer()->GetName() == "phi");
CHECK(read_back.GetGoniometer()->GetStart_deg() == Catch::Approx(30).margin(1e-3));
CHECK(!read_back.GetDatasetSettings().GetSmargonPosition().has_value());
// Writing what was read must not try to create phi a second time.
DiffractionExperiment rewrite = read_back;
rewrite.FilePrefix("test_phispindle_out").OverwriteExistingFiles(true);
StartMessage out_start;
rewrite.FillMessage(out_start);
FileWriter out(out_start);
EndMessage out_end;
out_end.max_image_number = rewrite.GetImageNum();
REQUIRE_NOTHROW(out.WriteHDF5(out_end));
REQUIRE_NOTHROW(out.Finalize());
remove("test_phispindle_master.h5");
remove("test_phispindle_data_000001.h5");
remove("test_phispindle_out_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// A valid NXmx master written outside the DECTRIS toolchain, in the shape a Diamond-written one
// takes: lengths in millimetres, no detectorSpecific, the distance one level up in NXinstrument, a
// pixel_mask that is an external link into a file not holding it, and per-file links naming a plain
// /data rather than /entry/data/data. Every one of those was enough on its own to stop the file
// opening, and the last one did it without an error - the run reported no images and succeeded.
namespace {
void WriteThirdPartyDataFile(const std::string &filename, const std::vector<uint16_t> &image,
hsize_t nimages, hsize_t ny, hsize_t nx) {
std::vector<uint16_t> block;
for (hsize_t i = 0; i < nimages; i++)
block.insert(block.end(), image.begin(), image.end());
HDF5File file(filename);
file.SaveVector("/data", block, {nimages, ny, nx});
}
}
TEST_CASE("JFJochReader_ThirdPartyNXmxMaster", "[HDF5][Full]") {
const hsize_t nx = 8, ny = 6, per_file = 2;
std::vector<uint16_t> image(nx * ny);
for (size_t i = 0; i < image.size(); i++)
image[i] = static_cast<uint16_t>(i * 3 + 1);
WriteThirdPartyDataFile("third_party_000001.h5", image, per_file, ny, nx);
WriteThirdPartyDataFile("third_party_000002.h5", image, per_file, ny, nx);
{
HDF5File master("third_party_master.h5");
HDF5Group entry(master, "entry");
entry.SaveScalar("definition", "NXmx");
HDF5Group instrument(entry, "instrument");
// The distance NXdetector does not carry, in millimetres
instrument.SaveScalar("detector_distance", 287.5)->Units("mm");
HDF5Group beam(instrument, "beam");
beam.SaveScalar("incident_wavelength", 0.9794)->Units("angstrom");
HDF5Group detector(instrument, "detector");
detector.SaveScalar("description", "Eiger 16M");
detector.SaveScalar("beam_center_x", 4.0)->Units("pixels");
detector.SaveScalar("beam_center_y", 3.0)->Units("pixels");
detector.SaveScalar("count_time", 0.01);
detector.SaveScalar("saturation_value", static_cast<int64_t>(65535));
detector.SaveScalar("x_pixel_size", 0.075)->Units("mm");
detector.SaveScalar("y_pixel_size", 0.075)->Units("mm");
detector.SaveScalar("sensor_thickness", 0.45)->Units("mm");
// Links into a file that does not exist at all, so neither can be dereferenced
detector.ExternalLink("third_party_no_such_meta.h5", "/mask", "pixel_mask");
HDF5Group data(entry, "data");
data.ExternalLink("third_party_000001.h5", "/data", "data_000001");
data.ExternalLink("third_party_000002.h5", "/data", "data_000002");
}
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("third_party_master.h5"));
auto dataset = reader.GetDataset();
// Images found through the link's own target path, and sized from the array itself
CHECK(dataset->experiment.GetImageNum() == 2 * per_file);
CHECK(dataset->experiment.GetXPixelsNum() == nx);
CHECK(dataset->experiment.GetYPixelsNum() == ny);
// Millimetres read as millimetres
CHECK(dataset->experiment.GetDetectorDistance_mm() == Catch::Approx(287.5));
CHECK(dataset->experiment.GetDetectorSetup().GetPixelSize_mm() == Catch::Approx(0.075));
CHECK(dataset->experiment.GetDetectorSetup().GetSensorThickness_um() == Catch::Approx(450.0));
// Both mask links dangle; the reader must fall back to an empty mask, not throw
REQUIRE(dataset->pixel_mask);
std::shared_ptr<JFJochReaderRawImage> reader_image;
for (int i = 0; i < 2 * static_cast<int>(per_file); i++) {
REQUIRE_NOTHROW(reader_image = reader.GetRawImage(i));
CHECK(reader_image->image.GetWidth() == nx);
CHECK(reader_image->image.GetHeight() == ny);
}
}
// A master that names data files nothing can be read from is an error, not an empty data set
remove("third_party_000001.h5");
remove("third_party_000002.h5");
{
JFJochHDF5Reader reader;
REQUIRE_THROWS(reader.ReadFile("third_party_master.h5"));
}
remove("third_party_master.h5");
// No leftover HDF5 objects
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// A detector swung out on a 2theta arm. NXmx has no field for it: the swing is one rotation in the
// depends_on chain the detector's position is stated as, and "two_theta" is only one beamline's name
// for that dataset. So the chain is what the reader follows, and the chain here carries two rotations
// about different axes, outboard of the translation that sets the distance - a file that stated only
// the innermost one, or composed them the other way round, gives a different plane.
//
// Both axes are stated in McStas, which is the internal frame turned half a turn about z: a reader
// that takes the vector as it stands swings the detector the wrong way, which is twice the error of
// not reading it at all.
TEST_CASE("JFJochReader_DetectorTwoThetaArm", "[HDF5][Full]") {
const hsize_t nx = 8, ny = 6;
const double two_theta_deg = 20.0, tilt_deg = 7.0;
std::vector<uint16_t> image(nx * ny, 5);
WriteThirdPartyDataFile("two_theta_000001.h5", image, 2, ny, nx);
{
HDF5File master("two_theta_master.h5");
HDF5Group entry(master, "entry");
entry.SaveScalar("definition", "NXmx");
HDF5Group instrument(entry, "instrument");
HDF5Group beam(instrument, "beam");
beam.SaveScalar("incident_wavelength", 0.6889)->Units("angstrom");
HDF5Group transformations(instrument, "transformations");
// Outermost first in the file, innermost first along the chain: det_z -> two_theta -> tilt
transformations.SaveVector("tilt", std::vector<double>{tilt_deg})
->Transformation("deg", ".", "detector", "", "rotation", {0, 1, 0});
transformations.SaveVector("two_theta", std::vector<double>{two_theta_deg})
->Transformation("deg", "/entry/instrument/transformations/tilt",
"detector", "", "rotation", {-1, 0, 0});
transformations.SaveVector("det_z", std::vector<double>{160.0})
->Transformation("mm", "/entry/instrument/transformations/two_theta",
"detector", "", "translation", {0, 0, 1});
HDF5Group detector(instrument, "detector");
detector.SaveScalar("depends_on", "/entry/instrument/transformations/det_z");
detector.SaveScalar("description", "PILATUS 2M");
detector.SaveScalar("beam_center_x", 4.0)->Units("pixels");
detector.SaveScalar("beam_center_y", 3.0)->Units("pixels");
detector.SaveScalar("distance", 0.160)->Units("m");
detector.SaveScalar("x_pixel_size", 0.172)->Units("mm");
detector.SaveScalar("y_pixel_size", 0.172)->Units("mm");
detector.SaveScalar("sensor_thickness", 0.32)->Units("mm");
detector.SaveScalar("count_time", 0.2);
detector.SaveScalar("saturation_value", static_cast<int64_t>(65535));
HDF5Group data(entry, "data");
data.ExternalLink("two_theta_000001.h5", "/data", "data_000001");
}
DiffractionGeometry geom;
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("two_theta_master.h5"));
geom = reader.GetDataset()->experiment.GetDiffractionGeometry();
}
// The chain as it stands in the internal frame: McStas (-1,0,0) is internal (1,0,0) and McStas
// (0,1,0) is internal (0,-1,0), and the outer rotation multiplies on the left.
const auto to_rad = [](double deg) { return static_cast<float>(deg * PI / 180.0); };
const RotMatrix expected = RotMatrix(to_rad(tilt_deg), {0, -1, 0})
* RotMatrix(to_rad(two_theta_deg), {1, 0, 0});
for (int64_t column = 0; column < 3; column++)
CHECK((geom.GetDetectorMatrix().Column(column) - expected.Column(column)).Length() < 1e-5f);
// Distance and beam centre are the ones the file states: the arm turns the detector about the
// sample and moves neither.
CHECK(geom.GetDetectorDistance_mm() == Catch::Approx(160.0));
CHECK(geom.GetBeamX_pxl() == Catch::Approx(4.0));
CHECK(geom.GetBeamY_pxl() == Catch::Approx(3.0));
// And the beam centre pixel is now that far from the beam - the whole point of a 2theta arm.
CHECK(geom.TwoTheta_rad(4.0f, 3.0f) * 180.0f / PI
== Catch::Approx(angle_deg(expected * Coord(0, 0, 1), Coord(0, 0, 1))));
CHECK(geom.TwoTheta_rad(4.0f, 3.0f) * 180.0f / PI > two_theta_deg);
remove("two_theta_000001.h5");
remove("two_theta_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// The same file with the arm parked at zero: the chain is there, nothing in it turns, and the
// geometry must be exactly the square-on one. This is nearly every file, so it has to cost nothing.
TEST_CASE("JFJochReader_DetectorTwoThetaZeroIsSquareOn", "[HDF5][Full]") {
const hsize_t nx = 8, ny = 6;
std::vector<uint16_t> image(nx * ny, 5);
WriteThirdPartyDataFile("two_theta_zero_000001.h5", image, 2, ny, nx);
{
HDF5File master("two_theta_zero_master.h5");
HDF5Group entry(master, "entry");
entry.SaveScalar("definition", "NXmx");
HDF5Group instrument(entry, "instrument");
HDF5Group beam(instrument, "beam");
beam.SaveScalar("incident_wavelength", 0.6889)->Units("angstrom");
HDF5Group transformations(instrument, "transformations");
transformations.SaveVector("two_theta", std::vector<double>{0.0})
->Transformation("deg", ".", "detector", "", "rotation", {-1, 0, 0});
transformations.SaveVector("det_z", std::vector<double>{160.0})
->Transformation("mm", "/entry/instrument/transformations/two_theta",
"detector", "", "translation", {0, 0, 1});
HDF5Group detector(instrument, "detector");
detector.SaveScalar("depends_on", "/entry/instrument/transformations/det_z");
detector.SaveScalar("description", "PILATUS 2M");
detector.SaveScalar("beam_center_x", 4.0)->Units("pixels");
detector.SaveScalar("beam_center_y", 3.0)->Units("pixels");
detector.SaveScalar("distance", 0.160)->Units("m");
detector.SaveScalar("x_pixel_size", 0.172)->Units("mm");
detector.SaveScalar("y_pixel_size", 0.172)->Units("mm");
detector.SaveScalar("sensor_thickness", 0.32)->Units("mm");
detector.SaveScalar("count_time", 0.2);
detector.SaveScalar("saturation_value", static_cast<int64_t>(65535));
HDF5Group data(entry, "data");
data.ExternalLink("two_theta_zero_000001.h5", "/data", "data_000001");
}
DiffractionGeometry geom;
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("two_theta_zero_master.h5"));
geom = reader.GetDataset()->experiment.GetDiffractionGeometry();
}
CHECK(geom.GetPoniRot1_rad() == 0.0f);
CHECK(geom.GetPoniRot2_rad() == 0.0f);
CHECK(geom.GetPoniRot3_rad() == 0.0f);
CHECK(geom.TwoTheta_rad(4.0f, 3.0f) == 0.0f);
remove("two_theta_zero_000001.h5");
remove("two_theta_zero_master.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// A file this system wrote states its own PONI angles twice: as the three scalars the reader takes
// them from, and as three rotations in the detector's depends_on chain. Following the chain must
// therefore skip them - applied on top of the scalars they would tilt the detector twice, which is
// how a correct 2theta reader breaks every tilted file this system has ever written. A test that only
// wrote an untilted detector could not see it.
TEST_CASE("JFJochReader_DetectorChainDoesNotDoubleTheTilt", "[HDF5][Full]") {
const float rot1 = 0.031f, rot2 = -0.047f, rot3 = 0.019f;
DiffractionExperiment x(DetJF(1));
x.ImagesPerTrigger(2).OverwriteExistingFiles(true).FilePrefix("test_ponichain");
x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
.IncidentEnergy_keV(WVL_1A_IN_KEV).PixelSigned(false).BitDepthImage(16)
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10));
x.PoniRot1_rad(rot1).PoniRot2_rad(rot2).PoniRot3_rad(rot3);
RegisterHDF5Filter();
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
StartMessage start_message;
x.FillMessage(start_message);
FileWriter file_set(start_message);
DataMessage message{};
for (int i = 0; i < x.GetImageNum(); i++) {
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
message.number = 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();
DiffractionGeometry geom;
{
JFJochHDF5Reader reader;
REQUIRE_NOTHROW(reader.ReadFile("test_ponichain_master.h5"));
geom = reader.GetDataset()->experiment.GetDiffractionGeometry();
}
CHECK(geom.GetPoniRot1_rad() == Catch::Approx(rot1).margin(1e-6));
CHECK(geom.GetPoniRot2_rad() == Catch::Approx(rot2).margin(1e-6));
CHECK(geom.GetPoniRot3_rad() == Catch::Approx(rot3).margin(1e-6));
for (int64_t column = 0; column < 3; column++)
CHECK((geom.GetDetectorMatrix().Column(column)
- PoniRotMatrix(rot1, rot2, rot3).Column(column)).Length() < 1e-5f);
remove("test_ponichain_master.h5");
remove("test_ponichain_data_000001.h5");
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
}
// A miniCBF header states more about the instrument than the "# " lines do: the CBF template block
// some beamlines write carries a full imgCIF axis table, saying which laboratory direction the image's
// columns and rows run along and which the spindle turns about. The reader assumed all three, and two
// instruments in the corpus are not what it assumed - one stores its image a quarter turn round, the
// other turns about the VERTICAL. Either way the spindle ends up 90 degrees from the image, which is
// not a sign and so is beyond the axis-sign rescue; both indexed nothing at all.
namespace {
// Two frames of a sweep whose pixels are all zero. Every delta of a zero image is zero, so the
// byte-offset stream is one 0x00 per pixel - which is a valid stream and enough to open a sweep.
void WriteMiniCBFSweep(const std::string &prefix, const std::string &header_body,
int64_t nx, int64_t ny) {
for (int frame = 1; frame <= 2; frame++) {
std::ostringstream head;
head << "###CBF: VERSION 1.5\n_array_data.header_convention \"PILATUS_1.2\"\n"
<< "_array_data.header_contents\n"
<< "# Detector: PILATUS3 6M, S/N 60-0119\n"
<< "# Pixel_size 172e-6 m x 172e-6 m\n"
<< "# Silicon sensor, thickness 0.000450 m\n"
<< "# Exposure_time 0.1 s\n# Exposure_period 0.1 s\n# Count_cutoff 768595 counts\n"
<< "# Wavelength 0.96864 A\n# Detector_distance 0.33161 m\n"
<< "# Beam_xy (12.00, 8.00) pixels\n"
<< "# Start_angle " << (frame - 1) * 0.1 << " deg.\n# Angle_increment 0.1000 deg.\n"
<< "# Omega " << (frame - 1) * 0.1 << " deg.\n# Omega_increment 0.1000 deg.\n"
<< "# Phi 0.0000 deg.\n# Phi_increment 0.0000 deg.\n"
<< "# Chi 0.0000 deg.\n# Chi_increment 0.0000 deg.\n"
<< header_body
<< "_array_data.data\n--CIF-BINARY-FORMAT-SECTION--\n"
<< "Content-Type: application/octet-stream;\n"
<< " conversions=\"x-CBF_BYTE_OFFSET\"\n"
<< "Content-Transfer-Encoding: BINARY\n"
<< "X-Binary-Size: " << nx * ny << "\n"
<< "X-Binary-Element-Type: \"signed 32-bit integer\"\n"
<< "X-Binary-Number-of-Elements: " << nx * ny << "\n"
<< "X-Binary-Size-Fastest-Dimension: " << nx << "\n"
<< "X-Binary-Size-Second-Dimension: " << ny << "\n\n";
std::ostringstream name;
name << prefix << "_" << std::setfill('0') << std::setw(4) << frame << ".cbf";
std::ofstream f(name.str(), std::ios::binary);
const std::string text = head.str();
f.write(text.data(), static_cast<std::streamsize>(text.size()));
f.write(reinterpret_cast<const char *>(minicbf::BINARY_SEPARATOR),
sizeof(minicbf::BINARY_SEPARATOR));
const std::vector<char> zeros(static_cast<size_t>(nx * ny), 0);
f.write(zeros.data(), static_cast<std::streamsize>(zeros.size()));
}
}
void RemoveMiniCBFSweep(const std::string &prefix) {
for (int frame = 1; frame <= 2; frame++) {
std::ostringstream name;
name << prefix << "_" << std::setfill('0') << std::setw(4) << frame << ".cbf";
remove(name.str().c_str());
}
}
// The axis table in the form these headers write it, several tags to a line.
std::string AxisTable(const std::string &rows, int64_t nx, int64_t ny) {
return "loop_\n_axis.id\n_axis.type\n_axis.equipment\n_axis.depends_on\n"
"_axis.vector[1] _axis.vector[2] _axis.vector[3]\n"
"_axis.offset[1] _axis.offset[2] _axis.offset[3]\n"
+ rows +
"loop_\n_array_structure_list.array_id\n_array_structure_list.index\n"
"_array_structure_list.dimension\n_array_structure_list.precedence\n"
"_array_structure_list.direction\n_array_structure_list.axis_set_id\n"
"ARRAY1 1 " + std::to_string(nx) + " 1 increasing ELEMENT_X\n"
"ARRAY1 2 " + std::to_string(ny) + " 2 increasing ELEMENT_Y\n"
"loop_\n_array_structure_list_axis.axis_set_id\n_array_structure_list_axis.axis_id\n"
"_array_structure_list_axis.displacement\n_array_structure_list_axis.displacement_increment\n"
"ELEMENT_X ELEMENT_X 0.0 0.1720\nELEMENT_Y ELEMENT_Y 0.0 0.1720\n";
}
}
TEST_CASE("JFJochCBFReader_AxisTableStatesTheMounting", "[HDF5][Full]") {
const int64_t nx = 24, ny = 16;
// A header that states nothing: the assumption, and the behaviour of nearly every file there is.
SECTION("no table, no hint - the assumption stands") {
WriteMiniCBFSweep("cbfaxis_plain", "# Detector_2theta 0.0000 deg.\n# Oscillation_axis OMEGA\n",
nx, ny);
JFJochCBFReader reader;
REQUIRE_NOTHROW(reader.ReadFiles("cbfaxis_plain_0001.cbf"));
const auto x = reader.GetDataset()->experiment;
REQUIRE(x.GetGoniometer().has_value());
CHECK((x.GetGoniometer()->GetAxis() - Coord(-1, 0, 0)).Length() < 1e-6f);
CHECK(x.GetDetectorSetup().GetImageOrientation().IsIdentity());
CHECK(x.GetDiffractionGeometry().GetPoniRot2_rad() == 0.0f);
reader.Close();
RemoveMiniCBFSweep("cbfaxis_plain");
}
// A spindle that turns about the VERTICAL, with the image mounted the usual way round. imgCIF Y is
// up and the internal frame's y is down, so the stated (0,1,0) is internal (0,-1,0) - and NOT the
// (-1,0,0) that was assumed, which is 90 degrees away and indexes nothing.
SECTION("vertical spindle, standard image") {
WriteMiniCBFSweep("cbfaxis_vert",
"# Detector_2theta 0.0000 deg.\n# Oscillation_axis X.CW +SLOW\n"
+ AxisTable("GON_OMEGA rotation goniometer . 0 1 0 . . .\n"
"DET_Z translation detector . 0 0 -1 0 0 0\n"
"ELEMENT_X translation detector DET_Z 1 0 0 -1 1 0\n"
"ELEMENT_Y translation detector ELEMENT_X 0 -1 0 0 0 0\n", nx, ny),
nx, ny);
JFJochCBFReader reader;
REQUIRE_NOTHROW(reader.ReadFiles("cbfaxis_vert_0001.cbf"));
const auto x = reader.GetDataset()->experiment;
REQUIRE(x.GetGoniometer().has_value());
CHECK((x.GetGoniometer()->GetAxis() - Coord(0, -1, 0)).Length() < 1e-6f);
// The image itself is standard, so nothing about it is turned - the axis was the whole error.
CHECK(x.GetDetectorSetup().GetImageOrientation().IsIdentity());
reader.Close();
RemoveMiniCBFSweep("cbfaxis_vert");
}
// The same vertical spindle, stated only by the "+SLOW" token, which is all a header with no axis
// table says. Two datasets from that instrument are in this state.
SECTION("vertical spindle from the +SLOW token alone") {
WriteMiniCBFSweep("cbfaxis_slow", "# Detector_2theta 0.0000 deg.\n# Oscillation_axis X.CW +SLOW\n",
nx, ny);
JFJochCBFReader reader;
REQUIRE_NOTHROW(reader.ReadFiles("cbfaxis_slow_0001.cbf"));
const auto x = reader.GetDataset()->experiment;
REQUIRE(x.GetGoniometer().has_value());
CHECK((x.GetGoniometer()->GetAxis() - Coord(0, -1, 0)).Length() < 1e-6f);
reader.Close();
RemoveMiniCBFSweep("cbfaxis_slow");
}
// An image stored a quarter turn round, on a detector swung out to 30 degrees. The two are read
// together or not at all: the arm turns about a laboratory axis, and which way that runs across
// the stored image is exactly what the mounting says.
SECTION("quarter-turned image on a swung arm") {
WriteMiniCBFSweep("cbfaxis_turn",
"# Detector_2theta 30.0000 deg.\n# Oscillation_axis OMEGA\n"
+ AxisTable("GON_OMEGA rotation goniometer . 1 0 0 . . .\n"
"DET_2THETA rotation detector . 1 0 0 . . .\n"
"DET_Z translation detector DET_2THETA 0 0 -1 0 0 0\n"
"ELEMENT_X translation detector DET_Z 0 1 0 -1 1 0\n"
"ELEMENT_Y translation detector ELEMENT_X 1 0 0 0 0 0\n", nx, ny),
nx, ny);
JFJochCBFReader reader;
REQUIRE_NOTHROW(reader.ReadFiles("cbfaxis_turn_0001.cbf"));
const auto x = reader.GetDataset()->experiment;
REQUIRE(x.GetGoniometer().has_value());
CHECK((x.GetGoniometer()->GetAxis() - Coord(1, 0, 0)).Length() < 1e-6f);
// fast = imgCIF (0,1,0) = internal (0,-1,0), slow = imgCIF (1,0,0) = internal (1,0,0)
CHECK(x.GetDetectorSetup().GetImageOrientation() == DetectorOrientation(false, 3));
// and the arm turns about its own stated axis, internal +x, by the stated 30 degrees
const auto geom = x.GetDiffractionGeometry();
const RotMatrix expected = RotMatrix(static_cast<float>(30.0 * PI / 180.0), {1, 0, 0})
* DetectorOrientation(false, 3).Matrix();
for (int64_t column = 0; column < 3; column++)
CHECK((geom.GetDetectorMatrix().Column(column) - expected.Column(column)).Length() < 1e-5f);
reader.Close();
RemoveMiniCBFSweep("cbfaxis_turn");
}
}
namespace {
// A byte-offset CBF whose header lines the caller chooses, so a header that is MISSING something
// can be built. WriteMiniCBFSweep above always writes a complete PILATUS head.
void WriteRawMiniCBF(const std::string &name, const std::string &head_lines,
int64_t nx, int64_t ny) {
std::ostringstream head;
head << "###CBF: VERSION 1.5\n_array_data.header_contents\n" << head_lines
<< "_array_data.data\n--CIF-BINARY-FORMAT-SECTION--\n"
<< "Content-Type: application/octet-stream;\n"
<< " conversions=\"x-CBF_BYTE_OFFSET\"\n"
<< "Content-Transfer-Encoding: BINARY\n"
<< "X-Binary-Size: " << nx * ny << "\n"
<< "X-Binary-Element-Type: \"signed 32-bit integer\"\n"
<< "X-Binary-Number-of-Elements: " << nx * ny << "\n"
<< "X-Binary-Size-Fastest-Dimension: " << nx << "\n"
<< "X-Binary-Size-Second-Dimension: " << ny << "\n\n";
std::ofstream f(name, std::ios::binary);
const std::string text = head.str();
f.write(text.data(), static_cast<std::streamsize>(text.size()));
f.write(reinterpret_cast<const char *>(minicbf::BINARY_SEPARATOR),
sizeof(minicbf::BINARY_SEPARATOR));
const std::vector<char> zeros(static_cast<size_t>(nx * ny), 0);
f.write(zeros.data(), static_cast<std::streamsize>(zeros.size()));
}
}
// Three ways a CBF that is not a detector image, or is one with a hole in its head, used to be
// opened anyway - each of them silently, which is the failure this project cares most about.
TEST_CASE("JFJochCBFReader_incomplete_header_is_refused_not_misread", "[HDF5][Full]") {
const int64_t nx = 24, ny = 16;
SECTION("a byte-offset CBF with no PILATUS header is not ours to read") {
// XDS writes its correction files in exactly this shape: a real byte-offset binary section
// and not one '#' line. Claiming it opened it with a pixel size of zero, which collapses
// every resolution and scattering vector the run computes.
WriteRawMiniCBF("cbfbare_0001.cbf", "", nx, ny);
CHECK_FALSE(JFJochCBFReader::CanRead("cbfbare_0001.cbf"));
JFJochCBFReader reader;
CHECK_THROWS_AS(reader.ReadFiles("cbfbare_0001.cbf"), JFJochException);
remove("cbfbare_0001.cbf");
}
SECTION("a header number that does not parse is a malformed header, not a raw std throw") {
// The captures are character classes, not number grammars: "[\\d.eE+-]+" matches a bare ".".
// std::stod answers that with std::invalid_argument, which CanRead does not catch - so merely
// LOOKING at the file threw out of the format probe.
WriteRawMiniCBF("cbfbadnum_0001.cbf",
"# Pixel_size 172e-6 m x 172e-6 m\n# Wavelength . A\n"
"# Detector_distance 0.3 m\n# Count_cutoff 1000 counts\n", nx, ny);
CHECK_NOTHROW(JFJochCBFReader::CanRead("cbfbadnum_0001.cbf"));
CHECK_FALSE(JFJochCBFReader::CanRead("cbfbadnum_0001.cbf"));
remove("cbfbadnum_0001.cbf");
}
SECTION("no Count_cutoff does not mean every pixel is saturated") {
// SaturationLimitFromValue(0) is 1, so an absent line marked every pixel at or above one
// count as an overload and the integration accept gate then dropped the whole reflection.
WriteMiniCBFSweep("cbfnocut", "# Oscillation_axis OMEGA\n", nx, ny);
// ...rewrite frame 1 without the Count_cutoff line, keeping everything else.
WriteRawMiniCBF("cbfnocut_0001.cbf",
"# Detector: PILATUS3 6M, S/N 60-0119\n"
"# Pixel_size 172e-6 m x 172e-6 m\n"
"# Silicon sensor, thickness 0.000450 m\n"
"# Wavelength 0.96864 A\n# Detector_distance 0.33161 m\n"
"# Beam_xy (12.00, 8.00) pixels\n"
"# Start_angle 0 deg.\n# Angle_increment 0.1000 deg.\n", nx, ny);
JFJochCBFReader reader;
REQUIRE_NOTHROW(reader.ReadFiles("cbfnocut_0001.cbf"));
const auto x = reader.GetDataset()->experiment;
CHECK(x.GetSaturationLimit() > 1);
reader.Close();
RemoveMiniCBFSweep("cbfnocut");
}
}
Reference in New Issue View Git Blame Copy Permalink