mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 78 Wiki Activity
Files
rc171
Jungfraujoch/tests/JFJochReaderTest.cpp
T
leonarski_fandClaude Opus 5 2ab8c55dfa rugnux and the viewer read SMV and gzipped miniCBF 
Two more of the formats deposited data actually arrives in, found by processing a
corpus of it: every PETRA III EMBL set is .cbf.gz, and NSRRC and the whole ADSC
Quantum era are SMV. Both were previously "no native input".

SMV is an ASCII "KEY=value;" block between braces, then the pixels - no container,
no compression, nothing to decode by offset - so reader/SMV.{h,cpp} and
JFJochSMVReader are a smaller job than the marCCD pair they sit beside, and need no
new dependency at all. Two things the format does not give us, both said out loud
rather than papered over:

* It states no saturation value, so overloads are judged on the 16-bit container
  alone. That can only fail to call a pixel saturated, never condemn a good one,
  but a CCD at the top of its range does saturate, so the reader warns once.
* Its beam centre is in MILLIMETRES and which of X/Y is the fast direction is a
  convention rather than a rule. Measured on one ALS ADSC sweep the file's value is
  TRANSPOSED: as stated it indexes 2/60 frames, and the run's own beam-centre
  measurement (which adopts the right one automatically) indexes 60/60. Swapping it
  here would fit that writer and might break another, so the header is read as the
  format defines it and the measurement stays the arbiter. Revisit with a second
  vendor's SMV in hand.

.cbf.gz needed only Slurp() in MiniCBF.cpp, through which every read already passes:
it sniffs the two-byte gzip magic - not the file name - and takes a zlib path when it
is there, leaving the plain path free of zlib's buffer copy. zlib-ng is already in the
build, so this is a link line, not a dependency. The sweep template grew a suffix,
because ".cbf" and ".cbf.gz" are separate sweeps and std::filesystem cannot split the
double extension on its own.

The viewer's single cbf_reader becomes three, dispatched by CanRead() in the same
order as rugnux. Dispatch is by CONTENT in both: ".img" is used by miniCBF, marCCD
AND SMV depending on the writer, and a PDB detector label has now been wrong about
the format four times, so an extension decides nothing.

Measured, de novo, no flags: 9fcg (1800 gzipped frames) gives P4 and a cell 0.06%
from the deposited one at 1.37 A against a deposited 1.54; 6oel (ADSC SMV) gives
F4132 - 96 operations, the most a protein space group can have - and a cell 0.05%
out, 100% indexed. Tests cover both formats and the transposed-beam-centre case with
fixtures written byte for byte, so they need no external data.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-09-16 22:35:05 +02:00

4410 lines
183 KiB
C++
Raw Permalink 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/JFJochMarCCDReader.h"
#include "../reader/JFJochSMVReader.h"
#include "../reader/MiniCBF.h"
#include "../compression/JFJochCompressor.h"
#include <cstring>
#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;
// Only one frame has a spindle severity: the others must come back ABSENT (stored as
// NaN), because no value is the CANNOT-SAY trigger state and a zero is not.
if (i == 2)
message.spindle_blind_fraction = 0.75f;
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(dataset->spindle_blind_fraction.size() == 4);
for (int i = 0; i < 4; i++) {
if (i == 2)
CHECK(dataset->spindle_blind_fraction[i] == Catch::Approx(0.75f));
else
CHECK(std::isnan(dataset->spindle_blind_fraction[i]));
}
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,
.prescaling_corr = 1.0f + 0.125f * f,
.qe_corr = 1.0f / (1.0f + 0.01f * f), // always in (0, 1], as the correction itself is
.flight_corr = 1.0f + 0.02f * f, // always >= 1, as the correction itself is
.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.prescaling_corr == Catch::Approx(want.prescaling_corr)); // stored as 1/prescaling_corr, inverted again on read
CHECK(got.qe_corr == Catch::Approx(want.qe_corr)); // likewise stored as its reciprocal
CHECK(got.flight_corr == Catch::Approx(want.flight_corr)); // likewise
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");
}
}
// marCCD: a TIFF whose instrument header sits in the gap between the TIFF header and the pixels.
// The fixtures below are written byte for byte rather than through libtiff, because the layout IS
// what the reader has to cope with - the private tag that states where the instrument header
// begins, and pixels that start well past the end of the IFD.
namespace {
constexpr size_t MARCCD_HEADER_OFF = 1024;
constexpr size_t MARCCD_DATA_OFF = 4096;
void PutU16(std::vector<uint8_t> &b, size_t off, uint16_t v) {
b[off] = v & 0xff; b[off + 1] = (v >> 8) & 0xff;
}
void PutU32(std::vector<uint8_t> &b, size_t off, uint32_t v) {
for (int i = 0; i < 4; i++) b[off + i] = (v >> (8 * i)) & 0xff;
}
void PutI32(std::vector<uint8_t> &b, size_t off, int32_t v) {
PutU32(b, off, static_cast<uint32_t>(v));
}
struct MarCCDFields {
int32_t start_phi_mdeg = 0;
int32_t end_phi_mdeg = 200;
int32_t rotation_range_mdeg = 200;
int32_t two_theta_mdeg = 0;
int32_t distance_um = 300000;
int32_t beam_x_mpx = 12000; // thousandths of a pixel
int32_t beam_y_mpx = 8000;
int32_t pixel_nm = 73242;
int32_t wavelength_e5A = 97872;
int32_t saturated = 65535;
int32_t depth = 2;
bool write_instrument_header = true;
};
// One frame. The image is all zeros: the reader is being asked about geometry, and a pixel
// pattern would say nothing extra about that.
void WriteMarCCDFrame(const std::string &name, int64_t nx, int64_t ny, const MarCCDFields &f) {
const size_t npixel = static_cast<size_t>(nx) * static_cast<size_t>(ny);
std::vector<uint8_t> file(MARCCD_DATA_OFF + npixel * 2, 0);
// TIFF header: little-endian, IFD at byte 8.
file[0] = 'I'; file[1] = 'I'; PutU16(file, 2, 42); PutU32(file, 4, 8);
struct Entry { uint16_t tag; uint16_t type; uint32_t count; uint32_t value; };
const std::vector<Entry> ifd = {
{256, 4, 1, static_cast<uint32_t>(nx)}, // ImageWidth
{257, 4, 1, static_cast<uint32_t>(ny)}, // ImageLength
{258, 3, 1, 16}, // BitsPerSample
{259, 3, 1, 1}, // Compression: none
{262, 3, 1, 1}, // Photometric: min-is-black
{273, 4, 1, MARCCD_DATA_OFF}, // StripOffsets
{277, 3, 1, 1}, // SamplesPerPixel
{278, 4, 1, static_cast<uint32_t>(ny)}, // RowsPerStrip: one strip
{279, 4, 1, static_cast<uint32_t>(npixel * 2)}, // StripByteCounts
{34710, 4, 1, MARCCD_HEADER_OFF}, // where the instrument header begins
};
PutU16(file, 8, static_cast<uint16_t>(ifd.size()));
size_t e = 10;
for (const auto &t : ifd) {
PutU16(file, e, t.tag); PutU16(file, e + 2, t.type); PutU32(file, e + 4, t.count);
// A SHORT value sits in the first two bytes of the value field, a LONG fills it.
if (t.type == 3) PutU16(file, e + 8, static_cast<uint16_t>(t.value));
else PutU32(file, e + 8, t.value);
e += 12;
}
PutU32(file, e, 0); // no next IFD
if (f.write_instrument_header) {
const size_t h = MARCCD_HEADER_OFF;
std::memcpy(file.data() + h + 4, "MMX", 3);
PutI32(file, h + 80, static_cast<int32_t>(nx)); // nfast
PutI32(file, h + 84, static_cast<int32_t>(ny)); // nslow
PutI32(file, h + 88, f.depth);
PutI32(file, h + 104, f.saturated);
PutI32(file, h + 640, f.distance_um);
PutI32(file, h + 644, f.beam_x_mpx);
PutI32(file, h + 648, f.beam_y_mpx);
PutI32(file, h + 656, 1000); // exposure, ms
PutI32(file, h + 668, f.two_theta_mdeg); // start_twotheta
PutI32(file, h + 684, f.start_phi_mdeg); // start_phi
PutI32(file, h + 700, f.two_theta_mdeg); // end_twotheta
PutI32(file, h + 716, f.end_phi_mdeg); // end_phi
PutI32(file, h + 736, f.rotation_range_mdeg);
PutI32(file, h + 772, f.pixel_nm);
PutI32(file, h + 776, f.pixel_nm);
PutI32(file, h + 908, f.wavelength_e5A);
}
std::ofstream out(name, std::ios::binary);
out.write(reinterpret_cast<const char *>(file.data()),
static_cast<std::streamsize>(file.size()));
}
// Two frames 0.2 deg apart, under whatever naming scheme the caller asks for.
void WriteMarCCDSweep(const std::string &prefix, const std::string &suffix,
int64_t nx, int64_t ny, MarCCDFields f = {}) {
for (int frame = 1; frame <= 2; frame++) {
f.start_phi_mdeg = (frame - 1) * 200;
f.end_phi_mdeg = frame * 200;
std::ostringstream name;
name << prefix << std::setfill('0') << std::setw(3) << frame << suffix;
WriteMarCCDFrame(name.str(), nx, ny, f);
}
}
void RemoveMarCCDSweep(const std::string &prefix, const std::string &suffix) {
for (int frame = 1; frame <= 2; frame++) {
std::ostringstream name;
name << prefix << std::setfill('0') << std::setw(3) << frame << suffix;
remove(name.str().c_str());
}
}
}
TEST_CASE("JFJochMarCCDReader_Geometry", "[HDF5][Full]") {
const int64_t nx = 24, ny = 16;
SECTION("the header's geometry reaches the experiment in the units the rest of the code uses") {
WriteMarCCDSweep("marccd_", ".mccd", nx, ny);
REQUIRE(JFJochMarCCDReader::CanRead("marccd_001.mccd"));
JFJochMarCCDReader reader;
REQUIRE_NOTHROW(reader.ReadFiles("marccd_001.mccd"));
CHECK(reader.GetNumberOfImages() == 2);
const auto x = reader.GetDataset()->experiment;
CHECK(x.GetDetectorDistance_mm() == Catch::Approx(300.0));
CHECK(x.GetBeamX_pxl() == Catch::Approx(12.0));
CHECK(x.GetBeamY_pxl() == Catch::Approx(8.0));
CHECK(x.GetPixelSize_mm() == Catch::Approx(0.073242).epsilon(1e-4));
CHECK(x.GetWavelength_A() == Catch::Approx(0.97872).epsilon(1e-5));
REQUIRE(x.GetGoniometer().has_value());
// The circle that moved between start and end is the scanned one, and the step comes from
// the two frames rather than from the header's nominal range.
CHECK(x.GetGoniometer()->GetName() == "phi");
CHECK(x.GetGoniometer()->GetIncrement_deg() == Catch::Approx(0.2));
CHECK(x.GetGoniometer()->GetStart_deg() == Catch::Approx(0.0));
reader.Close();
RemoveMarCCDSweep("marccd_", ".mccd");
}
SECTION("a frame number as the file extension names a sweep just as well") {
// What the mar software writes by default: D1.001, D1.002, ... The frame number is the last
// number in the name either way, which is the one rule the sweep template uses.
WriteMarCCDSweep("marccdnum.", "", nx, ny);
REQUIRE(JFJochMarCCDReader::CanRead("marccdnum.001"));
JFJochMarCCDReader reader;
REQUIRE_NOTHROW(reader.ReadFiles("marccdnum.001"));
CHECK(reader.GetNumberOfImages() == 2);
reader.Close();
RemoveMarCCDSweep("marccdnum.", "");
}
SECTION("the pixels come back as the int32 the rest of the code reads") {
WriteMarCCDSweep("marccdpix_", ".mccd", nx, ny);
JFJochMarCCDReader reader;
REQUIRE_NOTHROW(reader.ReadFiles("marccdpix_001.mccd"));
JFJochReaderRawImage image;
REQUIRE(reader.ReadRawImage(0, image));
CHECK(image.image.GetWidth() == nx);
CHECK(image.image.GetHeight() == ny);
CHECK(image.image.GetMode() == CompressedImageMode::Int32);
reader.Close();
RemoveMarCCDSweep("marccdpix_", ".mccd");
}
SECTION("a TIFF that is not marCCD is refused, not misread") {
// The name rule alone matches any numbered TIFF in the directory, so a file with no
// instrument header must not be opened with a pixel size and a distance of zero.
MarCCDFields f;
f.write_instrument_header = false;
WriteMarCCDFrame("marccdplain_001.tif", nx, ny, f);
CHECK_FALSE(JFJochMarCCDReader::CanRead("marccdplain_001.tif"));
JFJochMarCCDReader reader;
CHECK_THROWS(reader.ReadFiles("marccdplain_001.tif"));
remove("marccdplain_001.tif");
}
}
// SMV: an ASCII "KEY=value;" block between braces, then the pixels. The fixtures are written by
// hand here, as for marCCD, because the layout is what the reader has to cope with.
namespace {
struct SMVFields {
int nx = 24, ny = 16;
double pixel_mm = 0.1;
double distance_mm = 250.0;
double osc_start = 5.0, osc_range = 0.2;
double wavelength = 0.9998;
double beam_x_mm = 1.2, beam_y_mm = 0.8;
std::string byte_order = "little_endian";
bool write_header = true;
};
void WriteSMVFrame(const std::string &name, const SMVFields &f) {
std::ostringstream h;
h << "{\n";
if (f.write_header) {
h << "HEADER_BYTES= 512;\nDIM=2;\nBYTE_ORDER=" << f.byte_order << ";\n"
<< "TYPE=unsigned_short;\nSIZE1=" << f.nx << ";\nSIZE2=" << f.ny << ";\n"
<< "PIXEL_SIZE=" << f.pixel_mm << ";\nDISTANCE=" << f.distance_mm << ";\n"
<< "OSC_START=" << f.osc_start << ";\nOSC_RANGE=" << f.osc_range << ";\n"
<< "WAVELENGTH=" << f.wavelength << ";\n"
<< "BEAM_CENTER_X=" << f.beam_x_mm << ";\nBEAM_CENTER_Y=" << f.beam_y_mm << ";\n";
} else {
h << "SOMETHING_ELSE=1;\n"; // a brace block that is not an image header
}
h << "}\n";
std::string head = h.str();
head.resize(512, '\f'); // pad to HEADER_BYTES, as every writer does
std::ofstream out(name, std::ios::binary);
out.write(head.data(), static_cast<std::streamsize>(head.size()));
const std::vector<char> pixels(static_cast<size_t>(f.nx) * f.ny * 2, 0);
out.write(pixels.data(), static_cast<std::streamsize>(pixels.size()));
}
void WriteSMVSweep(const std::string &prefix, const std::string &suffix, SMVFields f = {}) {
for (int frame = 1; frame <= 2; frame++) {
f.osc_start = 5.0 + (frame - 1) * 0.2;
std::ostringstream name;
name << prefix << std::setfill('0') << std::setw(3) << frame << suffix;
WriteSMVFrame(name.str(), f);
}
}
void RemoveSMVSweep(const std::string &prefix, const std::string &suffix) {
for (int frame = 1; frame <= 2; frame++) {
std::ostringstream name;
name << prefix << std::setfill('0') << std::setw(3) << frame << suffix;
remove(name.str().c_str());
}
}
}
TEST_CASE("JFJochSMVReader_Geometry", "[HDF5][Full]") {
SECTION("the header's geometry reaches the experiment in the units the rest of the code uses") {
WriteSMVSweep("smv_", ".img");
REQUIRE(JFJochSMVReader::CanRead("smv_001.img"));
JFJochSMVReader reader;
REQUIRE_NOTHROW(reader.ReadFiles("smv_001.img"));
CHECK(reader.GetNumberOfImages() == 2);
const auto x = reader.GetDataset()->experiment;
CHECK(x.GetDetectorDistance_mm() == Catch::Approx(250.0));
CHECK(x.GetPixelSize_mm() == Catch::Approx(0.1));
CHECK(x.GetWavelength_A() == Catch::Approx(0.9998).epsilon(1e-4));
// The file states the beam centre in MILLIMETRES; the rest of the code wants pixels.
CHECK(x.GetBeamX_pxl() == Catch::Approx(12.0));
CHECK(x.GetBeamY_pxl() == Catch::Approx(8.0));
REQUIRE(x.GetGoniometer().has_value());
CHECK(x.GetGoniometer()->GetStart_deg() == Catch::Approx(5.0));
CHECK(x.GetGoniometer()->GetIncrement_deg() == Catch::Approx(0.2));
reader.Close();
RemoveSMVSweep("smv_", ".img");
}
SECTION("a brace block that is not an image header is refused") {
SMVFields f; f.write_header = false;
WriteSMVFrame("smvbare_001.img", f);
CHECK_FALSE(JFJochSMVReader::CanRead("smvbare_001.img"));
remove("smvbare_001.img");
}
}
TEST_CASE("JFJochCBFReader_GzippedSweep", "[HDF5][Full]") {
// EMBL Hamburg writes .cbf.gz by default. The reader decompresses in place, so a gzipped sweep
// opens exactly as a plain one does and needs no conversion step.
const int64_t nx = 24, ny = 16;
WriteMiniCBFSweep("cbfgz", "# Oscillation_axis OMEGA\n", nx, ny);
for (int frame = 1; frame <= 2; frame++) {
std::ostringstream name;
name << "cbfgz_" << std::setfill('0') << std::setw(4) << frame << ".cbf";
// gzip the frame in place, so the directory holds only the .cbf.gz form
const std::string cmd = "gzip -f " + name.str();
REQUIRE(std::system(cmd.c_str()) == 0);
}
REQUIRE(JFJochCBFReader::CanRead("cbfgz_0001.cbf.gz"));
JFJochCBFReader reader;
REQUIRE_NOTHROW(reader.ReadFiles("cbfgz_0001.cbf.gz"));
CHECK(reader.GetNumberOfImages() == 2);
const auto x = reader.GetDataset()->experiment;
CHECK(x.GetDetectorDistance_mm() == Catch::Approx(331.61));
CHECK(x.GetBeamX_pxl() == Catch::Approx(12.0));
CHECK(x.GetWavelength_A() == Catch::Approx(0.96864).epsilon(1e-5));
JFJochReaderRawImage image;
REQUIRE(reader.ReadRawImage(0, image));
CHECK(image.image.GetWidth() == nx);
CHECK(image.image.GetHeight() == ny);
reader.Close();
for (int frame = 1; frame <= 2; frame++) {
std::ostringstream name;
name << "cbfgz_" << std::setfill('0') << std::setw(4) << frame << ".cbf.gz";
remove(name.str().c_str());
}
}
Reference in New Issue View Git Blame Copy Permalink