diff --git a/common/DetectorSetup.cpp b/common/DetectorSetup.cpp index 04d2edbea..43873e1c4 100644 --- a/common/DetectorSetup.cpp +++ b/common/DetectorSetup.cpp @@ -386,6 +386,42 @@ std::optional DetectorSetup::GetSaturationLimit() const { return saturation_limit; } +DetectorSetup &DetectorSetup::CountRateCorrectionApplied(bool input) { + countrate_correction_applied = input; + return *this; +} + +bool DetectorSetup::IsCountRateCorrectionApplied() const { + return countrate_correction_applied; +} + +DetectorSetup &DetectorSetup::FlatfieldApplied(bool input) { + flatfield_applied = input; + return *this; +} + +bool DetectorSetup::IsFlatfieldApplied() const { + return flatfield_applied; +} + +DetectorSetup &DetectorSetup::CountRateCorrectionLookupTable(const std::vector &input) { + countrate_correction_lookup_table = input; + return *this; +} + +const std::vector &DetectorSetup::GetCountRateCorrectionLookupTable() const { + return countrate_correction_lookup_table; +} + +DetectorSetup &DetectorSetup::VirtualPixelInterpolationApplied(std::optional input) { + virtual_pixel_interpolation_applied = input; + return *this; +} + +std::optional DetectorSetup::IsVirtualPixelInterpolationApplied() const { + return virtual_pixel_interpolation_applied; +} + DetectorSetup &DetectorSetup::DECTRISROI(const std::string &input) { dectris_roi = input; return *this; diff --git a/common/DetectorSetup.h b/common/DetectorSetup.h index a7095ffc4..9bef69623 100644 --- a/common/DetectorSetup.h +++ b/common/DetectorSetup.h @@ -59,6 +59,13 @@ class DetectorSetup { std::optional saturation_limit; std::optional settings; + // Corrections the detector itself applied to the pixels. Not configured here: a DECTRIS + // detector reports them in its stream2 start message. + bool countrate_correction_applied = false; + std::vector countrate_correction_lookup_table; + bool flatfield_applied = false; + std::optional virtual_pixel_interpolation_applied; + DetectorSetup(std::shared_ptr geom, DetectorType detector_type, const std::string &description = "Detector", @@ -93,6 +100,10 @@ public: DetectorSetup& MinCountTime(std::chrono::nanoseconds input); DetectorSetup& MinThreshold_keV(float input); DetectorSetup& SaturationLimit(std::optional input); + DetectorSetup& CountRateCorrectionApplied(bool input); + DetectorSetup& CountRateCorrectionLookupTable(const std::vector &input); + DetectorSetup& FlatfieldApplied(bool input); + DetectorSetup& VirtualPixelInterpolationApplied(std::optional input); DetectorSetup& Description(const std::string &input); DetectorSetup& DECTRISROI(const std::string &input); DetectorSetup& DefaultSettings(const std::optional &input); @@ -128,6 +139,10 @@ public: [[nodiscard]] std::string GetDECTRISStream2Addr() const; [[nodiscard]] float GetMinThreshold_keV() const; [[nodiscard]] std::optional GetSaturationLimit() const; + [[nodiscard]] bool IsCountRateCorrectionApplied() const; + [[nodiscard]] const std::vector &GetCountRateCorrectionLookupTable() const; + [[nodiscard]] bool IsFlatfieldApplied() const; + [[nodiscard]] std::optional IsVirtualPixelInterpolationApplied() const; [[nodiscard]] std::string GetDECTRISROI() const; [[nodiscard]] std::optional GetDefaultSettings() const; [[nodiscard]] int32_t GetTempThreshold_degC() const; diff --git a/common/DiffractionExperiment.cpp b/common/DiffractionExperiment.cpp index 19cc77b5c..4e1bf3d9b 100644 --- a/common/DiffractionExperiment.cpp +++ b/common/DiffractionExperiment.cpp @@ -736,8 +736,10 @@ void DiffractionExperiment::FillMessage(StartMessage &message) const { message.summation = GetSummation(); message.user_data = GetHeaderAppendix(); - message.countrate_correction_enabled = false; - message.flatfield_enabled = false; + message.countrate_correction_enabled = detector.IsCountRateCorrectionApplied(); + message.countrate_correction_lookup_table = detector.GetCountRateCorrectionLookupTable(); + message.flatfield_enabled = detector.IsFlatfieldApplied(); + message.virtual_pixel_interpolation_enabled = detector.IsVirtualPixelInterpolationApplied(); message.goniometer = dataset.GetGoniometer(); message.grid_scan = dataset.GetGridScan(); diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index d6b970040..dd68afaee 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -212,6 +212,8 @@ struct HDF5DataSourceMessage { uint64_t source_first_image = 0; uint64_t virtual_first_image = 0; uint64_t image_count = 0; + // Set when the source is 4D, [image, channel, y, x]: the channel linked to + std::optional source_channel; }; struct StartMessage { @@ -248,6 +250,9 @@ struct StartMessage { bool pixel_signed; // user data bool countrate_correction_enabled; + // Maps a measured count c to its corrected value [c]; sent by a DECTRIS detector + std::vector countrate_correction_lookup_table; + std::optional virtual_pixel_interpolation_enabled; float incident_energy; float incident_wavelength; diff --git a/docs/CBOR.md b/docs/CBOR.md index a6316e71b..9bcb8c441 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -23,7 +23,9 @@ There are minor differences at the moment: | direct_beam_x | float (optional) | Where the undeflected beam lands on the detector, X \[pixels\]. Not the same point as `beam_center_x`, which is the PONI - the foot of the perpendicular from the sample - and separates from the beam position as soon as the detector is tilted. This is the number a program that asks for "the beam centre" (XDS `ORGX`, for one) wants | | | direct_beam_y | float (optional) | Where the undeflected beam lands on the detector, Y \[pixels\] (XDS `ORGY`) | | | countrate_correction_enabled | bool | Countrate correction enabled | X | +| countrate_correction_lookup_table | uint32 array (optional) | Maps a measured count c to its corrected value \[c\], as sent by a DECTRIS detector | X | | flatfield_enabled | bool | Flatfield enabled | X | +| virtual_pixel_interpolation_enabled | bool (optional) | Virtual pixel interpolation enabled, as reported by a DECTRIS detector | X | | number_of_images | uint64 | Number of images in the series | X | | image_size_x | uint64 | Image width \[pixels\] | X | | image_size_y | uint64 | Image height \[pixels\] | X | diff --git a/docs/HDF5.md b/docs/HDF5.md index 9259e68cf..2cdc4c73b 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -181,6 +181,8 @@ File-level HDF5 attributes `file_name`, `file_time`, `HDF5_Version` are also set | `flatfield_applied` | NXmx | | | `pixel_mask`, `pixel_mask_applied` | NXmx | `pixel_mask` is `[y, x]`, hard-linked from `detectorSpecific/pixel_mask` | | `countrate_correction_applied` | NXmx | | +| `countrate_correction_lookup_table` | NXmx | only when the detector sent one (DECTRIS) | +| `virtual_pixel_interpolation_applied` | NXmx | only when the detector reported it (DECTRIS) | | `number_of_cycles` | base | frame-summation factor | #### Why `bit_depth_readout` is the image depth diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index 55572f9a0..a1f5ea29d 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -244,6 +244,43 @@ namespace { memcpy(v.data(), ptr, len); } + // DECTRIS may send this one compressed, like its images + void GetCBORUInt32Array(CborValue &value, std::vector &v) { + if (GetCBORTag(value) != TagUnsignedInt32BitLE) + throw JFJochException(JFJochExceptionCategory::CBORError, "Incorrect array type tag"); + + if (!cbor_value_is_tag(&value)) { + auto [ptr, len] = GetCBORByteString(value); + if (len % sizeof(uint32_t)) + throw JFJochException(JFJochExceptionCategory::CBORError, "Size mismatch"); + v.resize(len / sizeof(uint32_t)); + memcpy(v.data(), ptr, len); + return; + } + + if (GetCBORTag(value) != TagDECTRISCompression) + throw JFJochException(JFJochExceptionCategory::CBORError, "Unsupported tag"); + CborValue array_value; + cborErr(cbor_value_enter_container(&value, &array_value)); + auto algorithm_text = GetCBORString(array_value); + CompressionAlgorithm algorithm; + if (algorithm_text == "bslz4") + algorithm = CompressionAlgorithm::BSHUF_LZ4; + else if (algorithm_text == "bszstd") + algorithm = CompressionAlgorithm::BSHUF_ZSTD; + else + throw JFJochException(JFJochExceptionCategory::CBORError, "Unsupported compression algorithm"); + GetCBORUInt(array_value); // element size, known from the type tag + auto [ptr, len] = GetCBORByteString(array_value); + cborErr(cbor_value_leave_container(&value, &array_value)); + + // The bitshuffle header starts with the uncompressed size in bytes + if (len < 12) + throw JFJochException(JFJochExceptionCategory::CBORError, "Compressed array too short"); + size_t nelements = bshuf_read_uint64_BE(const_cast(ptr)) / sizeof(uint32_t); + JFJochDecompress(v, algorithm, ptr, len, nelements); + } + void GetCBORInt32Array(CborValue &value, std::vector &v) { if (GetCBORTag(value) != TagSignedInt32BitLE) throw JFJochException(JFJochExceptionCategory::CBORError, "Incorrect array type tag"); @@ -1270,8 +1307,12 @@ namespace { message.number_of_images = GetCBORUInt(value); else if (key == "countrate_correction_enabled") message.countrate_correction_enabled = GetCBORBool(value); + else if (key == "countrate_correction_lookup_table") + GetCBORUInt32Array(value, message.countrate_correction_lookup_table); else if (key == "flatfield_enabled") message.flatfield_enabled = GetCBORBool(value); + else if (key == "virtual_pixel_interpolation_enabled") + message.virtual_pixel_interpolation_enabled = GetCBORBool(value); else if (key == "image_size_x") message.image_size_x = GetCBORUInt(value); else if (key == "image_size_y") diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index f8675a247..385723e54 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -165,6 +165,12 @@ inline void CBOR_ENC(CborEncoder &encoder, const char* key, const std::vector& v) { + cborErr(cbor_encode_text_stringz(&encoder, key)); + cborErr(cbor_encode_tag(&encoder, TagUnsignedInt32BitLE)); + cborErr(cbor_encode_byte_string(&encoder, reinterpret_cast(v.data()), v.size() * sizeof(uint32_t))); +} + inline void CBOR_ENC(CborEncoder &encoder, const char* key, const std::vector& v) { cborErr(cbor_encode_text_stringz(&encoder, key)); cborErr(cbor_encode_tag(&encoder, TagSignedInt32BitLE)); @@ -696,7 +702,10 @@ void CBORStream2Serializer::SerializeSequenceStart(const StartMessage& message) CBOR_ENC(mapEncoder, "direct_beam_x", message.direct_beam_x); CBOR_ENC(mapEncoder, "direct_beam_y", message.direct_beam_y); CBOR_ENC(mapEncoder, "countrate_correction_enabled", message.countrate_correction_enabled); + if (!message.countrate_correction_lookup_table.empty()) + CBOR_ENC(mapEncoder, "countrate_correction_lookup_table", message.countrate_correction_lookup_table); CBOR_ENC(mapEncoder, "flatfield_enabled", message.flatfield_enabled); + CBOR_ENC(mapEncoder, "virtual_pixel_interpolation_enabled", message.virtual_pixel_interpolation_enabled); CBOR_ENC(mapEncoder, "number_of_images", message.number_of_images); CBOR_ENC(mapEncoder, "image_size_x", message.image_size_x); CBOR_ENC(mapEncoder, "image_size_y", message.image_size_y); diff --git a/reader/HDF5ImageLocator.cpp b/reader/HDF5ImageLocator.cpp index 3284d4970..05d6f908e 100644 --- a/reader/HDF5ImageLocator.cpp +++ b/reader/HDF5ImageLocator.cpp @@ -12,7 +12,8 @@ namespace { const std::string &dataset, uint64_t source_first_image, uint64_t virtual_first_image, - uint64_t image_count) { + uint64_t image_count, + std::optional source_channel = {}) { if (image_count == 0) return; @@ -20,6 +21,7 @@ namespace { auto &last = ret.back(); if (last.filename == filename && last.dataset == dataset + && last.source_channel == source_channel && last.source_first_image + last.image_count == source_first_image && last.virtual_first_image + last.image_count == virtual_first_image) { last.image_count += image_count; @@ -32,9 +34,22 @@ namespace { .dataset = dataset, .source_first_image = source_first_image, .virtual_first_image = virtual_first_image, - .image_count = image_count + .image_count = image_count, + .source_channel = source_channel }); } + + // A 4D mapping, [image, channel, y, x], may cover only some of the channels; the reader wants + // the first one. + bool CoversFirstChannel(const HDF5VirtualDatasetMapping &mapping) { + return mapping.virtual_start.size() != 4 || mapping.virtual_start[1] == 0; + } + + std::optional SourceChannel(const HDF5VirtualDatasetMapping &mapping) { + if (mapping.source_start.size() == 4) + return mapping.source_start[1]; + return {}; + } } void HDF5ImageLocator::Configure(Layout layout) { @@ -71,10 +86,10 @@ HDF5ImageLocator::Location HDF5ImageLocator::Resolve(int64_t global_image) const && layout_.data_layout == HDF5DataSetLayout::VIRTUAL) { const auto image = static_cast(global_image); for (const auto &mapping: layout_.vds_mappings) { - if (!mapping.ContainsVirtualImage(image)) + if (!CoversFirstChannel(mapping) || !mapping.ContainsVirtualImage(image)) continue; return {OpenCached(mapping.filename), static_cast(mapping.SourceImage(image)), - mapping.filename, mapping.dataset}; + mapping.filename, mapping.dataset, SourceChannel(mapping).value_or(0)}; } throw JFJochException(JFJochExceptionCategory::HDF5, "Image not covered by /entry/data/data VDS mappings"); @@ -127,13 +142,14 @@ std::vector HDF5ImageLocator::GetSourceMapping(uint64_t f bool found = false; for (const auto &mapping: layout_.vds_mappings) { - if (!mapping.ContainsVirtualImage(virtual_image)) + if (!CoversFirstChannel(mapping) || !mapping.ContainsVirtualImage(virtual_image)) continue; const uint64_t source_image = mapping.SourceImage(virtual_image); const std::string dataset = mapping.dataset.empty() ? "/entry/data/data" : mapping.dataset; - AppendOrExtendSourceMapping(ret, mapping.filename, dataset, source_image, local_image, 1); + AppendOrExtendSourceMapping(ret, mapping.filename, dataset, source_image, local_image, 1, + SourceChannel(mapping)); found = true; break; } diff --git a/reader/HDF5ImageLocator.h b/reader/HDF5ImageLocator.h index 3e9585c5d..bc5bf6579 100644 --- a/reader/HDF5ImageLocator.h +++ b/reader/HDF5ImageLocator.h @@ -31,6 +31,10 @@ public: // Where the images sit INSIDE that file. /entry/data/data everywhere DECTRIS writes, but a // VDS names its source dataset and is free to name another one, so take it at its word. std::string dataset = "/entry/data/data"; + // Channel to read when the dataset is 4D, [image, channel, y, x], as the DECTRIS "hdf5 nexus + // v2024.2 nxmx" format writes it (one channel per threshold). Only the first channel of the + // master is read. + hsize_t channel = 0; }; // One data file of a legacy multi-file dataset, with the dataset the master's link names diff --git a/reader/HDF5ImageSource.cpp b/reader/HDF5ImageSource.cpp index 9fc2d1a64..4d01896af 100644 --- a/reader/HDF5ImageSource.cpp +++ b/reader/HDF5ImageSource.cpp @@ -4,6 +4,8 @@ #include "HDF5ImageSource.h" #include "../common/JFJochException.h" +#include + #ifdef _WIN32 #include #else @@ -97,22 +99,26 @@ HDF5ImageSource::GetDataset(const HDF5ImageLocator::Location &loc) const { HDF5DataType datatype(*entry.dataset); HDF5Dcpl dcpl(*entry.dataset); - if (dataspace.GetNumOfDimensions() != 3) + const auto rank = dataspace.GetNumOfDimensions(); + if (rank != 3 && rank != 4) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - loc.dataset + " dataset must be 3D"); + loc.dataset + " dataset must be 3D or 4D"); + entry.multichannel = (rank == 4); if (datatype.IsFloat()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Float datasets not supported at this time"); auto dim = dataspace.GetDimensions(); - entry.height = dim[1]; - entry.width = dim[2]; + entry.height = dim[rank - 2]; + entry.width = dim[rank - 1]; entry.mode = CalcImageMode(datatype.GetElemSize(), datatype.IsFloat(), datatype.IsSigned()); + // One chunk per image: [1, h, w], or [1, 1, h, w] for a multichannel dataset auto chunk_size = dcpl.GetChunking(); - entry.direct_chunk = (chunk_size.size() == 3) && (chunk_size[0] == 1) - && (chunk_size[1] == dim[1]) && (chunk_size[2] == dim[2]); + entry.direct_chunk = (chunk_size.size() == rank) + && std::all_of(chunk_size.begin(), chunk_size.end() - 2, [](hsize_t c) { return c == 1; }) + && (chunk_size[rank - 2] == entry.height) && (chunk_size[rank - 1] == entry.width); if (entry.direct_chunk) entry.algorithm = dcpl.GetCompression(); @@ -155,7 +161,9 @@ HDF5ImageSource::PrepareDirectRead(const HDF5ImageLocator::Location &loc) const if (!ds.raw) return {}; - const hsize_t coord[3] = {static_cast(loc.local_index), 0, 0}; + const hsize_t coord_3d[3] = {loc.local_index, 0, 0}; + const hsize_t coord_4d[4] = {loc.local_index, loc.channel, 0, 0}; + const hsize_t *coord = ds.multichannel ? coord_4d : coord_3d; unsigned filter_mask = 0; haddr_t address = HADDR_UNDEF; hsize_t size = 0; diff --git a/reader/HDF5ImageSource.h b/reader/HDF5ImageSource.h index c4cda0ced..b3bc77dff 100644 --- a/reader/HDF5ImageSource.h +++ b/reader/HDF5ImageSource.h @@ -76,12 +76,17 @@ public: template CompressedImage ReadImageAt(std::vector &buffer, const HDF5ImageLocator::Location &loc) const { const auto &ds = GetDataset(loc); - const std::vector start = {static_cast(loc.local_index), 0, 0}; + std::vector start = {static_cast(loc.local_index), 0, 0}; + std::vector size = {1, ds.height, ds.width}; + if (ds.multichannel) { + start.insert(start.begin() + 1, loc.channel); + size.insert(size.begin() + 1, 1); + } if (ds.direct_chunk) ds.dataset->ReadDirectChunk(buffer, start); else - ds.dataset->ReadVectorToU8(buffer, start, {1, ds.height, ds.width}); + ds.dataset->ReadVectorToU8(buffer, start, size); return {buffer.data(), buffer.size(), ds.width, ds.height, ds.mode, ds.algorithm}; } @@ -127,6 +132,7 @@ private: CompressedImageMode mode{}; CompressionAlgorithm algorithm = CompressionAlgorithm::NO_COMPRESSION; bool direct_chunk = false; + bool multichannel = false; // 4D: [image, channel, y, x] }; // Keyed by file AND dataset path: a master whose VDS sources are datasets in itself serves // several of them out of one file, and keying by file alone would hand back the wrong one. diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 7c7de1169..0444a668f 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -153,9 +153,10 @@ inline std::pair parse_bravais_lattice(const std::st return {cs, centering}; } +// [image, y, x], or [image, channel, y, x] as the DECTRIS "hdf5 nexus v2024.2 nxmx" format writes it std::vector GetDimension(HDF5Object &object, const std::string &path) { const auto dim = object.GetDimension(path); - if (dim.size() != 3) + if (dim.size() != 3 && dim.size() != 4) throw JFJochException(JFJochExceptionCategory::HDF5, "Wrong dimension of " + path); return dim; } @@ -174,9 +175,9 @@ std::vector ReadVDSImageMappings(HDF5Object &file, if (mapping.dataset.empty()) throw JFJochException(JFJochExceptionCategory::HDF5, "VDS mapping has empty source dataset name"); - if (mapping.virtual_start.size() != 3) + if (mapping.virtual_start.size() != 3 && mapping.virtual_start.size() != 4) throw JFJochException(JFJochExceptionCategory::HDF5, - "Only 3D image VDS mappings are supported"); + "Only 3D or 4D image VDS mappings are supported"); } return mappings; @@ -529,6 +530,8 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen size_t image_size_x = 0; size_t image_size_y = 0; + // Of a multichannel master; only its first channel is read + std::string first_channel; if (master_file->Exists("/entry/data/data")) { HDF5DataSet data_dataset(*master_file, "/entry/data/data"); @@ -537,8 +540,10 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen auto dim = GetDimension(*master_file, "/entry/data/data"); number_of_images = dim[0]; - image_size_y = dim[1]; - image_size_x = dim[2]; + image_size_y = dim[dim.size() - 2]; + image_size_x = dim[dim.size() - 1]; + if (dim.size() == 4) + first_channel = master_file->ReadElement("/entry/data/channel", 0).value_or(""); images_per_file = number_of_images; if (data_layout == HDF5DataSetLayout::VIRTUAL) @@ -646,8 +651,8 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen const auto dim = GetDimension(data_file, data_path); fimages = dim[0]; if (nfiles == 0) { - image_size_y = dim[1]; - image_size_x = dim[2]; + image_size_y = dim[dim.size() - 2]; + image_size_x = dim[dim.size() - 1]; } legacy_format_files.push_back({fname, data_path}); @@ -1112,8 +1117,11 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen // it, which a deposition shipping no mask at all leaves behind. The link is there, so // the name exists; only opening it says whether the array does. std::vector mask_tmp; - for (const char *name: {"/entry/instrument/detector/pixel_mask", - "/entry/instrument/detector/detectorSpecific/pixel_mask"}) { + // A multichannel master keeps the mask per channel + const std::string channel_mask = "/entry/instrument/detector/" + first_channel + "_channel/pixel_mask"; + for (const std::string &name: {channel_mask, + std::string("/entry/instrument/detector/pixel_mask"), + std::string("/entry/instrument/detector/detectorSpecific/pixel_mask")}) { if (mask_tmp.empty() && master_file->IsDataSet(name)) mask_tmp = master_file->ReadVector(name, {0, 0}, {image_size_y, image_size_x}); } diff --git a/receiver/JFJochReceiverLite.cpp b/receiver/JFJochReceiverLite.cpp index a83d63156..bd550ffa3 100644 --- a/receiver/JFJochReceiverLite.cpp +++ b/receiver/JFJochReceiverLite.cpp @@ -217,6 +217,14 @@ void JFJochReceiverLite::Configure(const StartMessage &msg) { experiment.Detector().SensorMaterial(msg.sensor_material); experiment.Detector().SensorThickness_um(msg.sensor_thickness * 1e6); experiment.Detector().SaturationLimit(SaturationLimitFromValue(msg.saturation_value)); + // The detector decides the corrections applied to the pixels (DECTRIS enables them by default); + // the outgoing start message is rebuilt from the experiment, so without copying them NXmx would + // record them as not applied. + experiment.Detector().CountRateCorrectionApplied(msg.countrate_correction_enabled); + experiment.Detector().CountRateCorrectionLookupTable(msg.countrate_correction_lookup_table); + experiment.Detector().FlatfieldApplied(msg.flatfield_enabled); + experiment.Detector().VirtualPixelInterpolationApplied(msg.virtual_pixel_interpolation_enabled); + experiment.ApplyPixelMask(msg.pixel_mask_enabled); // Images are forwarded byte-for-byte, so the stream's own image_dtype - not anything configured // locally - decides both the width and the sign the outgoing metadata must declare. Taking only // the width used to leave a signed stream declared unsigned in NXmx. diff --git a/tests/CBORTest.cpp b/tests/CBORTest.cpp index 3fc35b980..abd43b243 100644 --- a/tests/CBORTest.cpp +++ b/tests/CBORTest.cpp @@ -26,6 +26,8 @@ TEST_CASE("CBORSerialize_Start", "[CBOR]") { .bit_depth_readout = 16, .pixel_signed = true, .countrate_correction_enabled = true, + .countrate_correction_lookup_table = {0, 1, 3, 7}, + .virtual_pixel_interpolation_enabled = true, .incident_energy = 12400, .incident_wavelength = 0.988, .beam_size_x = 8e-5, @@ -161,6 +163,8 @@ TEST_CASE("CBORSerialize_Start", "[CBOR]") { CHECK(output_message.gain_file_names == message.gain_file_names); CHECK(output_message.countrate_correction_enabled == message.countrate_correction_enabled); CHECK(output_message.flatfield_enabled == message.flatfield_enabled); + CHECK(output_message.countrate_correction_lookup_table == message.countrate_correction_lookup_table); + CHECK(output_message.virtual_pixel_interpolation_enabled == message.virtual_pixel_interpolation_enabled); CHECK(output_message.write_master_file == message.write_master_file); CHECK(output_message.data_reduction_factor_serialmx == message.data_reduction_factor_serialmx); CHECK(output_message.experiment_group == message.experiment_group); @@ -1471,3 +1475,35 @@ TEST_CASE("CBORSerialize_End_Transformations", "[CBOR]") { CHECK(!chain[1].IsConstant()); CHECK(chain[2].IsConstant()); } + +TEST_CASE("CBORDeserialize_Start_CompressedLookupTable", "[CBOR]") { + // A DECTRIS detector may send the count rate correction table compressed, like its images + std::vector lut(1000); + for (int i = 0; i < lut.size(); i++) + lut[i] = 2 * i + 1; + + JFJochBitShuffleCompressor compressor(CompressionAlgorithm::BSHUF_LZ4); + auto compressed = compressor.Compress(lut); + + std::vector buffer(64 * 1024); + CborEncoder encoder, map_encoder, array_encoder; + cbor_encoder_init(&encoder, buffer.data(), buffer.size(), 0); + REQUIRE(cbor_encode_tag(&encoder, CborSignatureTag) == CborNoError); + REQUIRE(cbor_encoder_create_map(&encoder, &map_encoder, CborIndefiniteLength) == CborNoError); + REQUIRE(cbor_encode_text_stringz(&map_encoder, "type") == CborNoError); + REQUIRE(cbor_encode_text_stringz(&map_encoder, "start") == CborNoError); + REQUIRE(cbor_encode_text_stringz(&map_encoder, "countrate_correction_lookup_table") == CborNoError); + REQUIRE(cbor_encode_tag(&map_encoder, TagUnsignedInt32BitLE) == CborNoError); + REQUIRE(cbor_encode_tag(&map_encoder, TagDECTRISCompression) == CborNoError); + REQUIRE(cbor_encoder_create_array(&map_encoder, &array_encoder, 3) == CborNoError); + REQUIRE(cbor_encode_text_stringz(&array_encoder, "bslz4") == CborNoError); + REQUIRE(cbor_encode_uint(&array_encoder, sizeof(uint32_t)) == CborNoError); + REQUIRE(cbor_encode_byte_string(&array_encoder, compressed.data(), compressed.size()) == CborNoError); + REQUIRE(cbor_encoder_close_container(&map_encoder, &array_encoder) == CborNoError); + REQUIRE(cbor_encoder_close_container(&encoder, &map_encoder) == CborNoError); + buffer.resize(cbor_encoder_get_buffer_size(&encoder, buffer.data())); + + auto deserialized = CBORStream2Deserialize(buffer); + REQUIRE(deserialized->start_message); + CHECK(deserialized->start_message->countrate_correction_lookup_table == lut); +} diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index 46b4c411d..5403c9764 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -4490,3 +4490,146 @@ TEST_CASE("JFJochCBFReader_GzippedSweep", "[HDF5][Full]") { remove(name.str().c_str()); } } + +// The DECTRIS "hdf5 nexus v2024.2 nxmx" format: /entry/data/data is 4D, [image, channel, y, x], a +// VDS over the data files, with the pixel mask kept per channel. Only the first channel is read. +TEST_CASE("JFJochReader_DECTRISMultichannel", "[HDF5][Full]") { + const hsize_t nimages = 3, nchannels = 2, ny = 6, nx = 8; + std::vector data(nimages * nchannels * ny * nx); + for (size_t i = 0; i < data.size(); i++) + data[i] = static_cast(i); + + { + const hsize_t dims[4] = {nimages, nchannels, ny, nx}; + const hsize_t chunk[4] = {1, 1, ny, nx}; + hid_t file = H5Fcreate("dectris_fw2_data_000001.h5", H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + hid_t space = H5Screate_simple(4, dims, nullptr); + hid_t dcpl = H5Pcreate(H5P_DATASET_CREATE); + H5Pset_chunk(dcpl, 4, chunk); + hid_t lcpl = H5Pcreate(H5P_LINK_CREATE); + H5Pset_create_intermediate_group(lcpl, 1); + hid_t dset = H5Dcreate2(file, "/entry/data/data", H5T_STD_U16LE, space, lcpl, dcpl, H5P_DEFAULT); + REQUIRE(H5Dwrite(dset, H5T_NATIVE_UINT16, H5S_ALL, H5S_ALL, H5P_DEFAULT, data.data()) >= 0); + H5Dclose(dset); + H5Pclose(lcpl); + H5Pclose(dcpl); + H5Sclose(space); + H5Fclose(file); + } + + std::vector mask_threshold_1(ny * nx, 0), mask_threshold_2(ny * nx, 0); + mask_threshold_1[5] = 2; + mask_threshold_2[7] = 2; + + { + HDF5File master("dectris_fw2_master.h5"); + HDF5Group entry(master, "entry"); + entry.SaveScalar("definition", "NXmx"); + HDF5Group instrument(entry, "instrument"); + HDF5Group beam(instrument, "beam"); + beam.SaveScalar("incident_wavelength", 1.0)->Units("angstrom"); + HDF5Group detector(instrument, "detector"); + detector.SaveScalar("description", "Dectris EIGER2"); + detector.SaveScalar("beam_center_x", 4.0)->Units("pixels"); + detector.SaveScalar("beam_center_y", 3.0)->Units("pixels"); + detector.SaveScalar("distance", 0.1)->Units("m"); + detector.SaveScalar("count_time", 0.01); + detector.SaveScalar("saturation_value", static_cast(65535)); + detector.SaveScalar("x_pixel_size", 75e-6)->Units("m"); + detector.SaveScalar("y_pixel_size", 75e-6)->Units("m"); + detector.SaveScalar("sensor_thickness", 450e-6)->Units("m"); + HDF5Group channel_1(detector, "threshold_1_channel"); + channel_1.SaveVector("pixel_mask", mask_threshold_1, {ny, nx}); + HDF5Group channel_2(detector, "threshold_2_channel"); + channel_2.SaveVector("pixel_mask", mask_threshold_2, {ny, nx}); + HDF5Group data_group(entry, "data"); + data_group.SaveVector("channel", std::vector{"threshold_1", "threshold_2"}); + + // One mapping per channel, the second channel's first, so a reader that takes the first + // mapping covering an image reads the wrong threshold + const hsize_t dims[4] = {nimages, nchannels, ny, nx}; + const hsize_t count[4] = {nimages, 1, ny, nx}; + hid_t vspace = H5Screate_simple(4, dims, nullptr); + hid_t sspace = H5Screate_simple(4, dims, nullptr); + hid_t dcpl = H5Pcreate(H5P_DATASET_CREATE); + for (hsize_t c: {hsize_t(1), hsize_t(0)}) { + const hsize_t start[4] = {0, c, 0, 0}; + H5Sselect_hyperslab(vspace, H5S_SELECT_SET, start, nullptr, count, nullptr); + H5Sselect_hyperslab(sspace, H5S_SELECT_SET, start, nullptr, count, nullptr); + REQUIRE(H5Pset_virtual(dcpl, vspace, "dectris_fw2_data_000001.h5", "/entry/data/data", sspace) >= 0); + } + H5Sselect_all(vspace); + hid_t dset = H5Dcreate2(data_group.GetID(), "data", H5T_STD_U16LE, vspace, H5P_DEFAULT, dcpl, H5P_DEFAULT); + REQUIRE(dset >= 0); + H5Dclose(dset); + H5Pclose(dcpl); + H5Sclose(sspace); + H5Sclose(vspace); + } + + std::vector source_data; + { + JFJochHDF5Reader reader; + REQUIRE_NOTHROW(reader.ReadFile("dectris_fw2_master.h5")); + auto dataset = reader.GetDataset(); + CHECK(dataset->experiment.GetImageNum() == nimages); + CHECK(dataset->experiment.GetXPixelsNum() == nx); + CHECK(dataset->experiment.GetYPixelsNum() == ny); + REQUIRE(dataset->pixel_mask); + CHECK(dataset->pixel_mask->GetMask() == mask_threshold_1); + + std::shared_ptr reader_image; + for (hsize_t i = 0; i < nimages; i++) { + REQUIRE_NOTHROW(reader_image = reader.GetRawImage(i)); + CHECK(reader_image->image.GetWidth() == nx); + CHECK(reader_image->image.GetHeight() == ny); + REQUIRE(reader_image->image_buffer.size() == ny * nx * sizeof(uint16_t)); + // Channel 0 of image i + CHECK(memcmp(reader_image->image_buffer.data(), data.data() + i * nchannels * ny * nx, + ny * nx * sizeof(uint16_t)) == 0); + } + + // A _process.h5 links its pictures to the source channel the reader used + source_data = reader.GetHDF5DataSource(0, nimages); + REQUIRE(source_data.size() == 1); + CHECK(source_data[0].source_channel == 0); + } + + { + DiffractionExperiment x(DetDECTRIS(nx, ny, "Test", {})); + x.FilePrefix("dectris_fw2_process").ImagesPerTrigger(nimages) + .SetFileWriterFormat(FileWriterFormat::NXmxIntegrated).OverwriteExistingFiles(true); + StartMessage start_message; + x.FillMessage(start_message); + start_message.bit_depth_image = 16; + start_message.pixel_signed = false; + start_message.write_images = false; + start_message.write_master_file = true; + start_message.hdf5_source_data = source_data; + + FileWriter writer(start_message); + std::vector image(ny * nx); + for (hsize_t i = 0; i < nimages; i++) { + DataMessage message{}; + message.number = i; + message.image = CompressedImage(image, nx, ny); + REQUIRE_NOTHROW(writer.WriteHDF5(message)); + } + EndMessage end_message; + end_message.max_image_number = nimages; + writer.WriteHDF5(end_message); + writer.Finalize(); + } + + { + HDF5ReadOnlyFile file("dectris_fw2_process_master.h5"); + auto image_1 = file.ReadVector("/entry/data/data", {1, 0, 0}, {1, ny, nx}); + CHECK(std::equal(image_1.begin(), image_1.end(), data.begin() + nchannels * ny * nx)); + } + + remove("dectris_fw2_process_master.h5"); + remove("dectris_fw2_master.h5"); + remove("dectris_fw2_data_000001.h5"); + // No leftover HDF5 objects + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); +} diff --git a/tests/JFJochReceiverLiteTest.cpp b/tests/JFJochReceiverLiteTest.cpp index 3d19fb506..f03001877 100644 --- a/tests/JFJochReceiverLiteTest.cpp +++ b/tests/JFJochReceiverLiteTest.cpp @@ -50,6 +50,13 @@ TEST_CASE("JFJochReceiverLite", "[JFJochReceiver]") { StartMessage start_msg; experiment.FillMessage(start_msg); + // As a DECTRIS detector sends it by default; the flags must reach the master file. + start_msg.countrate_correction_enabled = true; + start_msg.flatfield_enabled = true; + start_msg.countrate_correction_lookup_table = {0, 1, 3, 7}; + start_msg.virtual_pixel_interpolation_enabled = true; + // What the detector did to the pixels, not the local setting + start_msg.pixel_mask_enabled = !experiment.IsApplyPixelMask(); puller->Put(ImagePullerOutput{ .cbor = std::make_shared(start_msg) @@ -100,6 +107,14 @@ TEST_CASE("JFJochReceiverLite", "[JFJochReceiver]") { // No progress value at the end of the measurement REQUIRE(!service.GetProgress().has_value()); + + HDF5ReadOnlyFile master("crystal_test_lite_master.h5"); + CHECK(master.GetBool("/entry/instrument/detector/countrate_correction_applied")); + CHECK(master.GetBool("/entry/instrument/detector/flatfield_applied")); + CHECK(master.GetBool("/entry/instrument/detector/virtual_pixel_interpolation_applied")); + CHECK(master.GetBool("/entry/instrument/detector/pixel_mask_applied") == start_msg.pixel_mask_enabled); + CHECK(master.ReadVector("/entry/instrument/detector/countrate_correction_lookup_table") + == std::vector{0, 1, 3, 7}); } TEST_CASE("JFJochReceiverLite_Cancel", "[JFJochReceiver]") { diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index 5416cb5e8..6f99af71f 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -264,11 +264,17 @@ void NXmx::LinkToData_ProcessingVDS(const StartMessage &start, const EndMessage ); const hsize_t source_extent_images = mapping.source_first_image + image_count; - HDF5DataSpace source_data_space({source_extent_images, height, width}); - source_data_space.SelectHyperslab( - {static_cast(mapping.source_first_image), 0, 0}, - {image_count, height, width} - ); + std::vector source_dims = {source_extent_images, height, width}; + std::vector source_start = {static_cast(mapping.source_first_image), 0, 0}; + std::vector source_size = {image_count, height, width}; + if (mapping.source_channel) { + // One channel of a 4D source, [image, channel, y, x] + source_dims.insert(source_dims.begin() + 1, mapping.source_channel.value() + 1); + source_start.insert(source_start.begin() + 1, mapping.source_channel.value()); + source_size.insert(source_size.begin() + 1, 1); + } + HDF5DataSpace source_data_space(source_dims); + source_data_space.SelectHyperslab(source_start, source_size); dcpl.SetVirtual(mapping.filename, source_dataset, @@ -422,6 +428,10 @@ void NXmx::Detector(const StartMessage &start) { SaveScalar(group, "acquisition_type", "triggered"); SaveScalar(group, "countrate_correction_applied", start.countrate_correction_enabled); + if (!start.countrate_correction_lookup_table.empty()) + group.SaveVector("countrate_correction_lookup_table", start.countrate_correction_lookup_table); + if (start.virtual_pixel_interpolation_enabled) + SaveScalar(group, "virtual_pixel_interpolation_applied", start.virtual_pixel_interpolation_enabled.value()); SaveScalar(group, "number_of_cycles", start.summation); HDF5Group det_specific(group, "detectorSpecific");