diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index e5abba6e..da24ebe8 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -24,7 +24,7 @@ #include "XrayFluorescenceSpectrum.h" #include "../symmetry/gemmi/symmetry.hpp" -constexpr const uint64_t user_data_release = 5; +constexpr const uint64_t user_data_release = 6; constexpr const uint64_t user_data_magic_number = 0x52320000UL | user_data_release; enum class CBORImageType {START, END, IMAGE, CALIBRATION, METADATA, NONE}; @@ -168,6 +168,15 @@ struct DataMessage { std::optional lattice_type; }; +struct HDF5DataSourceMessage { + std::string filename; + std::string dataset = "/entry/data/data"; + + uint64_t source_first_image = 0; + uint64_t virtual_first_image = 0; + uint64_t image_count = 0; +}; + struct StartMessage { float detector_distance; float beam_center_x; @@ -283,6 +292,8 @@ struct StartMessage { XrayFluorescenceSpectrum fluorescence_spectrum; std::optional detect_ice_rings; + + std::vector hdf5_source_data; }; struct EndMessage { @@ -306,7 +317,26 @@ struct EndMessage { std::optional rotation_lattice_type; std::optional rotation_lattice; - std::vector scale_factor; + // Vectors with end result: + std::vector data_collection_efficiency; + std::vector spot_count; + std::vector spot_count_ice_ring; + std::vector spot_count_low_res; + std::vector spot_count_indexed; + std::vector image_indexed; + std::vector v_bkg_estimate; + std::vector profile_radius; + std::vector mosaicity; + std::vector bFactor; + std::vector resolution_estimate; + std::vector min_viable_pixel_value; + std::vector max_viable_pixel_value; + std::vector saturated_pixel_count; + std::vector error_pixel_count; + std::vector image_scale_factor; + std::vector integrated_reflections; + std::vector niggli_class; + std::vector pixel_sum; }; struct MetadataMessage { diff --git a/common/ScanResult.h b/common/ScanResult.h index 62453123..28686df9 100644 --- a/common/ScanResult.h +++ b/common/ScanResult.h @@ -22,6 +22,7 @@ struct ScanResultElem { std::optional angle_deg; std::optional pixel_sum; + std::optional min_viable_pixel; std::optional max_viable_pixel; std::optional err_pixels; std::optional sat_pixels; @@ -36,6 +37,10 @@ struct ScanResultElem { std::optional res; std::optional uc; std::optional xfel_pulse_id; + std::optional mosaicity; + std::optional image_scale_factor; + std::optional niggli_class; + std::optional integrated_reflections; }; struct ScanResult { diff --git a/common/ScanResultGenerator.cpp b/common/ScanResultGenerator.cpp index 5c2c9df9..202bf068 100644 --- a/common/ScanResultGenerator.cpp +++ b/common/ScanResultGenerator.cpp @@ -1,8 +1,18 @@ // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include +#include + #include "ScanResultGenerator.h" +namespace { + template + T value_or_zero(const std::optional& v) { + return v.value_or(static_cast(0)); + } +} + ScanResultGenerator::ScanResultGenerator(const DiffractionExperiment &experiment) { grid_scan = experiment.GetGridScan(); goniometer_axis = experiment.GetGoniometer(); @@ -21,7 +31,7 @@ void ScanResultGenerator::Add(const DataMessage &message) { if (grid_scan) image_number = grid_scan->Rearrange(image_number); - if (image_number < v.size()) { + if (image_number >= 0 && static_cast(image_number) < v.size()) { if (grid_scan) { v[image_number].x = grid_scan->GetElementPosX_step(message.number); v[image_number].y = grid_scan->GetElementPosY_step(message.number); @@ -36,16 +46,22 @@ void ScanResultGenerator::Add(const DataMessage &message) { v[image_number].spot_count = message.spot_count; v[image_number].indexing_solution = message.indexing_result; v[image_number].profile_radius = message.profile_radius; + v[image_number].mosaicity = message.mosaicity_deg; v[image_number].b_factor = message.b_factor; v[image_number].uc = message.indexing_unit_cell; v[image_number].xfel_pulse_id = message.xfel_pulse_id; v[image_number].err_pixels = message.error_pixel_count; + v[image_number].min_viable_pixel = message.min_viable_pixel_value; v[image_number].max_viable_pixel = message.max_viable_pixel_value; v[image_number].sat_pixels = message.saturated_pixel_count; v[image_number].spot_count_ice = message.spot_count_ice_rings; v[image_number].spot_count_low_res = message.spot_count_low_res; v[image_number].spot_count_indexed = message.spot_count_indexed; v[image_number].res = message.resolution_estimate; + v[image_number].integrated_reflections = message.integrated_reflections; + + if (message.lattice_type) + v[image_number].niggli_class = message.lattice_type->niggli_class; } } @@ -56,3 +72,64 @@ ScanResult ScanResultGenerator::GetResult() const { ret.images = v; return ret; } + +void ScanResultGenerator::FillEndMessage(EndMessage &message) const { + std::unique_lock ul(m); + + size_t n = 0; + for (const auto &e: v) { + if (e.number >= 0) + n = std::max(n, static_cast(e.number) + 1); + } + if (n == 0) + return; + + message.data_collection_efficiency.resize(n); + message.spot_count.resize(n); + message.spot_count_ice_ring.resize(n); + message.spot_count_low_res.resize(n); + message.spot_count_indexed.resize(n); + message.image_indexed.resize(n); + message.v_bkg_estimate.resize(n); + message.profile_radius.resize(n); + message.mosaicity.resize(n); + message.bFactor.resize(n); + message.resolution_estimate.resize(n); + message.min_viable_pixel_value.resize(n); + message.max_viable_pixel_value.resize(n); + message.saturated_pixel_count.resize(n); + message.error_pixel_count.resize(n); + message.image_scale_factor.resize(n); + message.integrated_reflections.resize(n); + message.niggli_class.resize(n); + message.pixel_sum.resize(n); + + for (const auto &e: v) { + if (e.number < 0) + continue; + + const auto number = static_cast(e.number); + if (number >= n) + continue; + + message.data_collection_efficiency[number] = e.collection_efficiency; + message.spot_count[number] = static_cast(value_or_zero(e.spot_count)); + message.spot_count_ice_ring[number] = static_cast(value_or_zero(e.spot_count_ice)); + message.spot_count_low_res[number] = static_cast(value_or_zero(e.spot_count_low_res)); + message.spot_count_indexed[number] = static_cast(value_or_zero(e.spot_count_indexed)); + message.image_indexed[number] = static_cast(e.indexing_solution.value_or(0)); + message.v_bkg_estimate[number] = e.bkg.value_or(NAN); + message.profile_radius[number] = e.profile_radius.value_or(NAN); + message.mosaicity[number] = e.mosaicity.value_or(NAN); + message.bFactor[number] = e.b_factor.value_or(NAN); + message.resolution_estimate[number] = e.res.value_or(NAN); + message.min_viable_pixel_value[number] = value_or_zero(e.min_viable_pixel); + message.max_viable_pixel_value[number] = value_or_zero(e.max_viable_pixel); + message.saturated_pixel_count[number] = static_cast(value_or_zero(e.sat_pixels)); + message.error_pixel_count[number] = static_cast(value_or_zero(e.err_pixels)); + message.image_scale_factor[number] = e.image_scale_factor.value_or(NAN); + message.integrated_reflections[number] = static_cast(value_or_zero(e.integrated_reflections)); + message.niggli_class[number] = static_cast(value_or_zero(e.niggli_class)); + message.pixel_sum[number] = static_cast(value_or_zero(e.pixel_sum)); + } +} diff --git a/common/ScanResultGenerator.h b/common/ScanResultGenerator.h index 76be34e0..c7e59c99 100644 --- a/common/ScanResultGenerator.h +++ b/common/ScanResultGenerator.h @@ -25,6 +25,8 @@ public: explicit ScanResultGenerator(const DiffractionExperiment& experiment); void Add(const DataMessage& message); ScanResult GetResult() const; + + void FillEndMessage(EndMessage &message) const; }; diff --git a/docs/CBOR.md b/docs/CBOR.md index b8a785a6..cafb7292 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -12,106 +12,106 @@ There are minor differences at the moment: ## Start message -| Field name | Type | Description | Present in DECTRIS format | -|----------------------------------|----------------------|-----------------------------------------------------------------------------------------------------------------------------------|:-------------------------:| -| type | String | value "start" | X | -| magic_number | uint64 | Number used to describe version of the Jungfraujoch data interface - to allow to detect inconsistency between sender and receiver | | -| detector_distance | float | Detector distance \[m\] | | -| detector_translation | Array(float) | Detector translation vector \[m\] | X | -| beam_center_x | float | Beam center in X direction \[pixels\] | X | -| beam_center_y | float | Beam center in Y direction \[pixels\] | X | -| countrate_correction_enabled | bool | Countrate correction enabled | X | -| flatfield_enabled | bool | Flatfield enabled | 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 | -| incident_energy | float | X-ray energy \[eV\] | X | -| incident_wavelength | float | X-ray wavelength \[Angstrom\] | X | -| frame_time | float | Frame time, if multiple frames per trigger \[s\] | X | -| count_time | float | Exposure time \[s\] | X | -| saturation_value | int64 | Maximum valid sample value | X | -| error_value | int64 (optional) | Value used in images to describe pixels that are in error state or missing | | -| pixel_size_x | float | Pixel width \[m\] | X | -| pixel_size_y | float | Pixel height \[m\] | X | -| sensor_thickness | float | Sensor thickness \[m\] | X | -| sensor_material | string | Sensor material | X | -| arm_date | date | Approximate date of arming | X | -| pixel_mask_enabled | bool | Pixel mask applied on images | X | -| detector_description | string | Name of the detector | X | -| detector_serial_number | string | Detector serial number | X | -| series_unique_id | string | Unique text ID of the series (run_name parameter) | X | -| series_id | uint64 | Unique numeric ID of the series (run_number parameter) | X | -| fluorescence | object (optional) | X-ray fluorescence spectrum collected at start | | -| - energy | Array(float) | Energy of measuring point \[eV\] | | -| - data | Array(float) | Fluorescence scan result `data` \[arbitrary units\]; must be strictly the same length as energy | | -| goniometer | Map | Definition of rotation axis (optional) | X | -| - `AXIS` | string | Rotation axis name (e.g. omega) - only one axis is supported in Jungfraujoch | X | -| - - increment | float | Rotation axis increment (per image) in degree \[deg\] | X | -| - - start | float | Rotation axis start angle \[deg\] | X | -| - - axis | Array(float) | Vector for the rotation axis | | -| - - helical_step | Array(float) | Translation for helical scan for 1 image \[m\] | | -| - - screening_wedge | Array(float) | Wedge for screening \[deg\] (increment would correspond to difference between screening points) | | -| grid_scan | object | Grid scan definition (optional and exclusive with rotation axis) | | -| - n_fast | uint64 | Number of elements along fast axis | | -| - n_slow | uint64 | Number of elements along slow axis | | -| - step_x_axis | float | Step along X axis, can be negative \[m\] | | -| - step_y_axis | float | Step along Y axis, can be negative \[m\] | | -| - snake_scan | bool | Snake scan (rows alternate direction) | | -| - vertical_scan | bool | Vertical scan (enabled: fast direction = Y, disabled: fast direction = X) | | -| jungfrau_conversion_enabled | bool (optional) | Applying JUNGFRAU pixel conversion (to photons or keV) | | -| jungfrau_conversion_factor | float (optional) | Factor used for JUNGFRAU conversion \[eV\] | | -| geometry_transformation_enabled | bool (optional) | Transformation from detector module geometry (512x1024) to full detector geometry | | -| pixel_mask | Map(string -> Image) | Pixel mask - multiple in case of storage cells | X | -| channels | Array(string) | List of image channels | X | -| max_spot_count | uint64 | Maximum number of spots identified in spot finding | | -| storage_cell_number | uint64 (optional) | Number of storage cells used by JUNGFRAU | | -| storage_cell_delay | Rational | Delay of storage cells in JUNGFRAU | | -| threshold_energy | float | Threshold energy for EIGER detector \[eV\] | | -| image_dtype | string | Pixel bit type (e.g. uint16) | X | -| unit_cell | object (optional) | Unit cell of the system: a, b, c \[angstrom\] and alpha, beta, gamma \[degree\] | | -| az_int_q_bin_count | uint64 | Number of azimuthal integration bins in the radial direction | | -| az_int_phi_bin_count | uint64 | Number of azimuthal integration bins in the phi angle direction | | -| az_int_bin_to_q | Array(float) | Q value for each azimuthal integration bin \[angstrom^-1\] | | -| az_int_bin_to_two_theta | Array(float) | Two theta angle value for each azimuthal integration bin \[deg\] | | -| az_int_bin_to_phi | Array(float) | Phi value for each azimuthal integration bin \[deg\] | | -| summation | uint64 | Factor of frame summation | | -| user_data | string | JSON serialized to string that can contain the following fields (all fields are optional): | X | -| - file_prefix | string | File prefix | | -| - images_per_file | uint64 | Number of images written per file | | -| - images_per_trigger | uint64 | Number of images collected per trigger | | -| - source_name | string | Facility name | | -| - source_type | string | Type of X-ray source (use NXsource/type values, for example "Synchrotron X-ray Source" or "Free-Electron Laser") | | -| - instrument_name | string | Instrument name | | -| - sample_name | string | Name of the sample | | -| - user | any valid JSON | Value of header_appendix provided at collection start to Jungfraujoch | | -| - attenuator_transmission | float | Attenuator transmission \[\] | | -| - total_flux | float | Total flux \[ph/s\] | | -| - space_group_number | uint64 | Space group number | | -| - summation_mode | string | Summation mode (internal\|fpga\|cpu) | | -| - overwrite | bool | Overwrite existing HDF5 files | | -| - file_format | int | File writer format: 0 = no master file, 1 = soft links, 2 = virtual dataset, 3 = CBF, 4 = TIFF | | -| - roi | Array(object) | ROI configurations; each element is one of: | | -| | | type "box": xmin, xmax, ymin, ymax (numbers) | | -| | | type "circle": r, x, y (numbers) | | -| | | type "azim": qmin, qmax (numbers) | | -| - gain_file_names | Array(string) | Names of JUNGFRAU gain files used for the current detector | | -| - write_master_file | bool | With multiple sockets, it selects which socket will provide master file | | -| - write_images | bool | Write images in the HDF5 file (if false, will only write metadata) | | -| - data_reduction_factor_serialmx | uint64 | Data reduction factor for serial MX | | -| - experiment_group | string | ID of instrument user, e.g., p-group (SLS/SwissFEL) or proposal number | | -| - jfjoch_release | string | Jungfraujoch release number | | -| - socket_number | uint64 | Number of ZeroMQ socket (on `jfjoch_broker` side) used for transmission | | -| - bit_depth_readout | uint64 | Bit depth of the detector readout | | -| - writer_notification_zmq_addr | string | ZeroMQ address to inform `jfjoch_broker` about writers that finished operation | | -| - xfel_pulse_id | uint64 | Pulse IDs are recorded for images | | -| - ring_current_mA | float | Ring current at the start of the measurement | | -| - sample_temperature_K | float | Sample temperature \[K\] | | -| - detect_ice_rings | bool | Ice ring detection feature is enabled | | -| - indexing_algorithm | string | Indexing algorithm used on-the-fly; allowed values: ffbidx, fft, fftw, none | | -| - geom_refinement_algorithm | string | Post-indexing detector geometry refinement algorithm; allowed values: none, beam_center, beam_center_tetragonal | | -| - poni_rot1 | float | Tilt of the detector rot1 according to PyFAI PONI convention \[rad\] | | -| - poni_rot2 | float | Tilt of the detector rot2 according to PyFAI PONI convention \[rad\] | | -| - poni_rot3 | float | Tilt of the detector rot3 according to PyFAI PONI convention \[rad\] | | +| Field name | Type | Description | Present in DECTRIS format | +|----------------------------------|----------------------|------------------------------------------------------------------------------------------------------------------------------------------------|:-------------------------:| +| type | String | value "start" | X | +| magic_number | uint64 | Number used to describe version of the Jungfraujoch data interface - to allow to detect inconsistency between sender and receiver | | +| detector_distance | float | Detector distance \[m\] | | +| detector_translation | Array(float) | Detector translation vector \[m\] | X | +| beam_center_x | float | Beam center in X direction \[pixels\] | X | +| beam_center_y | float | Beam center in Y direction \[pixels\] | X | +| countrate_correction_enabled | bool | Countrate correction enabled | X | +| flatfield_enabled | bool | Flatfield enabled | 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 | +| incident_energy | float | X-ray energy \[eV\] | X | +| incident_wavelength | float | X-ray wavelength \[Angstrom\] | X | +| frame_time | float | Frame time, if multiple frames per trigger \[s\] | X | +| count_time | float | Exposure time \[s\] | X | +| saturation_value | int64 | Maximum valid sample value | X | +| error_value | int64 (optional) | Value used in images to describe pixels that are in error state or missing | | +| pixel_size_x | float | Pixel width \[m\] | X | +| pixel_size_y | float | Pixel height \[m\] | X | +| sensor_thickness | float | Sensor thickness \[m\] | X | +| sensor_material | string | Sensor material | X | +| arm_date | date | Approximate date of arming | X | +| pixel_mask_enabled | bool | Pixel mask applied on images | X | +| detector_description | string | Name of the detector | X | +| detector_serial_number | string | Detector serial number | X | +| series_unique_id | string | Unique text ID of the series (run_name parameter) | X | +| series_id | uint64 | Unique numeric ID of the series (run_number parameter) | X | +| fluorescence | object (optional) | X-ray fluorescence spectrum collected at start | | +| - energy | Array(float) | Energy of measuring point \[eV\] | | +| - data | Array(float) | Fluorescence scan result `data` \[arbitrary units\]; must be strictly the same length as energy | | +| goniometer | Map | Definition of rotation axis (optional) | X | +| - `AXIS` | string | Rotation axis name (e.g. omega) - only one axis is supported in Jungfraujoch | X | +| - - increment | float | Rotation axis increment (per image) in degree \[deg\] | X | +| - - start | float | Rotation axis start angle \[deg\] | X | +| - - axis | Array(float) | Vector for the rotation axis | | +| - - helical_step | Array(float) | Translation for helical scan for 1 image \[m\] | | +| - - screening_wedge | Array(float) | Wedge for screening \[deg\] (increment would correspond to difference between screening points) | | +| grid_scan | object | Grid scan definition (optional and exclusive with rotation axis) | | +| - n_fast | uint64 | Number of elements along fast axis | | +| - n_slow | uint64 | Number of elements along slow axis | | +| - step_x_axis | float | Step along X axis, can be negative \[m\] | | +| - step_y_axis | float | Step along Y axis, can be negative \[m\] | | +| - snake_scan | bool | Snake scan (rows alternate direction) | | +| - vertical_scan | bool | Vertical scan (enabled: fast direction = Y, disabled: fast direction = X) | | +| jungfrau_conversion_enabled | bool (optional) | Applying JUNGFRAU pixel conversion (to photons or keV) | | +| jungfrau_conversion_factor | float (optional) | Factor used for JUNGFRAU conversion \[eV\] | | +| geometry_transformation_enabled | bool (optional) | Transformation from detector module geometry (512x1024) to full detector geometry | | +| pixel_mask | Map(string -> Image) | Pixel mask - multiple in case of storage cells | X | +| channels | Array(string) | List of image channels | X | +| max_spot_count | uint64 | Maximum number of spots identified in spot finding | | +| storage_cell_number | uint64 (optional) | Number of storage cells used by JUNGFRAU | | +| storage_cell_delay | Rational | Delay of storage cells in JUNGFRAU | | +| threshold_energy | float | Threshold energy for EIGER detector \[eV\] | | +| image_dtype | string | Pixel bit type (e.g. uint16) | X | +| unit_cell | object (optional) | Unit cell of the system: a, b, c \[angstrom\] and alpha, beta, gamma \[degree\] | | +| az_int_q_bin_count | uint64 | Number of azimuthal integration bins in the radial direction | | +| az_int_phi_bin_count | uint64 | Number of azimuthal integration bins in the phi angle direction | | +| az_int_bin_to_q | Array(float) | Q value for each azimuthal integration bin \[angstrom^-1\] | | +| az_int_bin_to_two_theta | Array(float) | Two theta angle value for each azimuthal integration bin \[deg\] | | +| az_int_bin_to_phi | Array(float) | Phi value for each azimuthal integration bin \[deg\] | | +| summation | uint64 | Factor of frame summation | | +| user_data | string | JSON serialized to string that can contain the following fields (all fields are optional): | X | +| - file_prefix | string | File prefix | | +| - images_per_file | uint64 | Number of images written per file | | +| - images_per_trigger | uint64 | Number of images collected per trigger | | +| - source_name | string | Facility name | | +| - source_type | string | Type of X-ray source (use NXsource/type values, for example "Synchrotron X-ray Source" or "Free-Electron Laser") | | +| - instrument_name | string | Instrument name | | +| - sample_name | string | Name of the sample | | +| - user | any valid JSON | Value of header_appendix provided at collection start to Jungfraujoch | | +| - attenuator_transmission | float | Attenuator transmission \[\] | | +| - total_flux | float | Total flux \[ph/s\] | | +| - space_group_number | uint64 | Space group number | | +| - summation_mode | string | Summation mode (internal\|fpga\|cpu) | | +| - overwrite | bool | Overwrite existing HDF5 files | | +| - file_format | int | File writer format: 0 = only data files, 1 = NXmx legacy soft links, 2 = NXmx VDS, 3 = NXmx integrated, 4 = CBF, 5 = TIFF, 6 = no file written | | +| - roi | Array(object) | ROI configurations; each element is one of: | | +| | | type "box": xmin, xmax, ymin, ymax (numbers) | | +| | | type "circle": r, x, y (numbers) | | +| | | type "azim": qmin, qmax (numbers) | | +| - gain_file_names | Array(string) | Names of JUNGFRAU gain files used for the current detector | | +| - write_master_file | bool | With multiple sockets, it selects which socket will provide master file | | +| - write_images | bool | Write images in the HDF5 file (if false, will only write metadata) | | +| - data_reduction_factor_serialmx | uint64 | Data reduction factor for serial MX | | +| - experiment_group | string | ID of instrument user, e.g., p-group (SLS/SwissFEL) or proposal number | | +| - jfjoch_release | string | Jungfraujoch release number | | +| - socket_number | uint64 | Number of ZeroMQ socket (on `jfjoch_broker` side) used for transmission | | +| - bit_depth_readout | uint64 | Bit depth of the detector readout | | +| - writer_notification_zmq_addr | string | ZeroMQ address to inform `jfjoch_broker` about writers that finished operation | | +| - xfel_pulse_id | uint64 | Pulse IDs are recorded for images | | +| - ring_current_mA | float | Ring current at the start of the measurement | | +| - sample_temperature_K | float | Sample temperature \[K\] | | +| - detect_ice_rings | bool | Ice ring detection feature is enabled | | +| - indexing_algorithm | string | Indexing algorithm used on-the-fly; allowed values: ffbidx, fft, fftw, none | | +| - geom_refinement_algorithm | string | Post-indexing detector geometry refinement algorithm; allowed values: none, beam_center, beam_center_tetragonal | | +| - poni_rot1 | float | Tilt of the detector rot1 according to PyFAI PONI convention \[rad\] | | +| - poni_rot2 | float | Tilt of the detector rot2 according to PyFAI PONI convention \[rad\] | | +| - poni_rot3 | float | Tilt of the detector rot3 according to PyFAI PONI convention \[rad\] | | See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/stream_v2) for definition of Image as MultiDimArray with optional compression. @@ -253,28 +253,49 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s ## End message -| Field name | Type | Description | Present in DECTRIS format | -|----------------------------|--------------------------|-----------------------------------------------------------------------------------------------------------------------------------|:-------------------------:| -| type | String | value "end" | X | -| magic_number | uint64 | Number used to describe version of the Jungfraujoch data interface - to allow to detect inconsistency between sender and receiver | | -| series_unique_id | string | Unique text ID of the series (run_name parameter) | X | -| series_id | uint64 | Unique numeric ID of the series (run_number parameter) | X | -| end_date | date | Approximate end date | | -| max_image_number | uint64 | Number of image with the highest number (this is counted from 1 - to distinguish zero images and one image) | | -| images_collected | uint64 | Number of image collected | | -| images_sent_to_write | uint64 | Number of image sent to writer; if writer queues were full, it is possible this is less than images collected | | -| data_collection_efficiency | float | Network packets collected / Network packets expected \[\] | | -| az_int_result | Map(text->Array(float)) | Azimuthal integration results, use az_int_bin_to_q from start message for legend | | -| adu_histogram | Map(text->Array(uint64)) | ADU values histogram | | -| adu_histogram_bin_width | uint64 | Width of bins in the above histogram \[ADU\] | | -| max_receiver_delay | uint64 | Internal performance of Jungfraujoch | | -| bkg_estimate | float | Mean background estimate for the whole run | | -| indexing_rate | float | Mean indexing rate for the whole run | | -| rotation_lattice_type | object | Bravais lattice classification of the total rotation solution over the run (if available); same schema as `lattice_type` | | -| - centering | string | One-letter centering code: P, A, B, C, I, F, or R | | -| - niggli_class | int64 | Integer identifier for the Niggli-reduced Bravais class | | -| - system | string | Crystal system: triclinic, monoclinic, orthorhombic, tetragonal, trigonal, hexagonal, cubic | | -| rotation_lattice | Array(9 * float) | Real-space lattice basis (flattened 3x3 in row-major), corresponding to the rotation indexing result | | +| Field name | Type | Description | Present in DECTRIS format | +|------------------------------------|--------------------------|-----------------------------------------------------------------------------------------------------------------------------------|:-------------------------:| +| type | String | value "end" | X | +| magic_number | uint64 | Number used to describe version of the Jungfraujoch data interface - to allow to detect inconsistency between sender and receiver | | +| series_unique_id | string | Unique text ID of the series (run_name parameter) | X | +| series_id | uint64 | Unique numeric ID of the series (run_number parameter) | X | +| end_date | date | Approximate end date | | +| max_image_number | uint64 | Number of image with the highest number; counted from 1 to distinguish zero images and one image | | +| images_collected | uint64 | Number of images collected | | +| images_sent_to_write | uint64 | Number of images sent to writer; if writer queues were full, it is possible this is less than images collected | | +| data_collection_efficiency | float | Overall network packets collected / network packets expected | | +| az_int_result | Map(text->Array(float)) | Azimuthal integration results, use az_int_bin_to_q from start message for legend | | +| adu_histogram | Map(text->Array(uint64)) | ADU values histogram | | +| adu_histogram_bin_width | uint64 | Width of bins in the above histogram \[ADU\] | | +| max_receiver_delay | uint64 | Internal performance of Jungfraujoch | | +| bkg_estimate | float | Mean background estimate for the whole run | | +| indexing_rate | float | Mean indexing rate for the whole run | | +| rotation_lattice_type | object | Bravais lattice classification of the total rotation solution over the run, if available; same schema as `lattice_type` | | +| - centering | string | One-letter centering code: P, A, B, C, I, F, or R | | +| - niggli_class | int64 | Integer identifier for the Niggli-reduced Bravais class | | +| - system | string | Crystal system: triclinic, monoclinic, orthorhombic, tetragonal, trigonal, hexagonal, cubic | | +| rotation_lattice | Array(9 * float) | Real-space lattice basis, flattened 3x3 in row-major order | | +| data_collection_efficiency_image | Array(float) | Per-image data collection efficiency. Missing values are encoded as 0 or 1 depending on producer context | | +| spot_count | Array(int32) | Per-image spot count | | +| spot_count_ice_ring | Array(int32) | Per-image number of spots within identified ice-ring resolution ranges | | +| spot_count_low_res | Array(int32) | Per-image number of low-resolution spots | | +| spot_count_indexed | Array(int32) | Per-image number of spots fitting indexing solution | | +| image_indexed | Array(uint8) | Per-image indexing result; 0 = not indexed, nonzero = indexed | | +| v_bkg_estimate | Array(float) | Per-image background estimate | | +| profile_radius | Array(float) | Per-image profile radius \[Angstrom^-1\] | | +| mosaicity | Array(float) | Per-image mosaicity \[degree\] | | +| bFactor | Array(float) | Per-image estimated B-factor \[Angstrom^2\] | | +| resolution_estimate | Array(float) | Per-image diffraction resolution estimate \[Angstrom\] | | +| min_viable_pixel_value | Array(int64) | Per-image minimum valid pixel value, excluding error/saturated pixels | | +| max_viable_pixel_value | Array(int64) | Per-image maximum valid pixel value, excluding error/saturated pixels | | +| saturated_pixel_count | Array(int32) | Per-image saturated pixel count | | +| error_pixel_count | Array(int32) | Per-image error pixel count | | +| image_scale_factor | Array(float) | Per-image scale factor, if scaling/merging was performed | | +| integrated_reflections | Array(int32) | Per-image count of integrated reflections | | +| niggli_class | Array(uint8) | Per-image Niggli class identifier for indexed images; 0 if unavailable | | +| pixel_sum | Array(int64) | Per-image sum of all valid pixels, excluding error/saturated pixels | | + +End-message vector fields are optional. When present, they provide master-file summary data so readers can inspect scan-level and per-image analysis results without opening every linked data file. Missing optional per-image values are encoded by the producer as zero unless otherwise noted. ## Calibration message @@ -301,4 +322,16 @@ Therefore `user_data` is serialized by Jungfraujoch as CBOR object. There is mem - Compression: - Uncompressed: raw CBOR byte string - Bitshuffle+LZ4: tag with \["bslz4", elem_size, bytes\] - - Bitshuffle+Zstandard: tag with \["bszstd", elem_size, bytes\] \ No newline at end of file + - Bitshuffle+Zstandard: tag with \["bszstd", elem_size, bytes\] + +### Notes on typed arrays + +Jungfraujoch uses RFC 8746-style typed byte-string tags for compact numeric arrays. + +Common tags used in this protocol include: + +- float32 little-endian arrays for `Array(float)` +- uint8 arrays for compact boolean/integer flags such as `image_indexed` +- int32 little-endian arrays for per-image counts +- int64 little-endian arrays for large per-image integer values +- uint64 little-endian arrays for histograms \ No newline at end of file diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index e60188d4..6b4d9fd6 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -228,6 +228,42 @@ namespace { return {ptr, len}; } + void GetCBORUInt8Array(CborValue &value, std::vector &v) { + if (GetCBORTag(value) != TagUnsignedInt8Bit) + throw JFJochException(JFJochExceptionCategory::CBORError, "Incorrect array type tag"); + + auto [ptr, len] = GetCBORByteString(value); + + v.resize(len); + memcpy(v.data(), ptr, len); + } + + void GetCBORInt32Array(CborValue &value, std::vector &v) { + if (GetCBORTag(value) != TagSignedInt32BitLE) + throw JFJochException(JFJochExceptionCategory::CBORError, "Incorrect array type tag"); + + auto [ptr, len] = GetCBORByteString(value); + + if (len % sizeof(int32_t)) + throw JFJochException(JFJochExceptionCategory::CBORError, "Size mismatch"); + + v.resize(len / sizeof(int32_t)); + memcpy(v.data(), ptr, len); + } + + void GetCBORInt64Array(CborValue &value, std::vector &v) { + if (GetCBORTag(value) != TagSignedInt64BitLE) + throw JFJochException(JFJochExceptionCategory::CBORError, "Incorrect array type tag"); + + auto [ptr, len] = GetCBORByteString(value); + + if (len % sizeof(int64_t)) + throw JFJochException(JFJochExceptionCategory::CBORError, "Size mismatch"); + + v.resize(len / sizeof(int64_t)); + memcpy(v.data(), ptr, len); + } + void GetCBORFloatArray(CborValue &value, std::vector &v) { if (GetCBORTag(value) != TagFloatLE) throw JFJochException(JFJochExceptionCategory::CBORError, "Incorrect array type tag"); @@ -1284,8 +1320,44 @@ namespace { message.bkg_estimate = GetCBORFloat(value); else if (key == "indexing_rate") message.indexing_rate = GetCBORFloat(value); + else if (key == "data_collection_efficiency_image") + GetCBORFloatArray(value, message.data_collection_efficiency); + else if (key == "spot_count") + GetCBORInt32Array(value, message.spot_count); + else if (key == "spot_count_ice_ring") + GetCBORInt32Array(value, message.spot_count_ice_ring); + else if (key == "spot_count_low_res") + GetCBORInt32Array(value, message.spot_count_low_res); + else if (key == "spot_count_indexed") + GetCBORInt32Array(value, message.spot_count_indexed); + else if (key == "image_indexed") + GetCBORUInt8Array(value, message.image_indexed); + else if (key == "v_bkg_estimate") + GetCBORFloatArray(value, message.v_bkg_estimate); + else if (key == "profile_radius") + GetCBORFloatArray(value, message.profile_radius); + else if (key == "mosaicity") + GetCBORFloatArray(value, message.mosaicity); + else if (key == "bFactor") + GetCBORFloatArray(value, message.bFactor); + else if (key == "resolution_estimate") + GetCBORFloatArray(value, message.resolution_estimate); + else if (key == "min_viable_pixel_value") + GetCBORInt64Array(value, message.min_viable_pixel_value); + else if (key == "max_viable_pixel_value") + GetCBORInt64Array(value, message.max_viable_pixel_value); + else if (key == "saturated_pixel_count") + GetCBORInt32Array(value, message.saturated_pixel_count); + else if (key == "error_pixel_count") + GetCBORInt32Array(value, message.error_pixel_count); else if (key == "image_scale_factor") - GetCBORFloatArray(value, message.scale_factor); + GetCBORFloatArray(value, message.image_scale_factor); + else if (key == "integrated_reflections") + GetCBORInt32Array(value, message.integrated_reflections); + else if (key == "niggli_class") + GetCBORUInt8Array(value, message.niggli_class); + else if (key == "pixel_sum") + GetCBORInt64Array(value, message.pixel_sum); else if (key == "rotation_lattice") { std::vector tmp; GetCBORFloatArray(value, tmp); diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index 195e5e05..fcd39f18 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -153,6 +153,23 @@ 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, TagUnsignedInt8Bit)); + cborErr(cbor_encode_byte_string(&encoder, v.data(), v.size() * sizeof(uint8_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)); + cborErr(cbor_encode_byte_string(&encoder, reinterpret_cast(v.data()), v.size() * sizeof(int32_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, TagSignedInt64BitLE)); + cborErr(cbor_encode_byte_string(&encoder, reinterpret_cast(v.data()), v.size() * sizeof(int64_t))); +} inline void CBOR_ENC(CborEncoder &encoder, const char* key, const std::vector& v) { cborErr(cbor_encode_text_stringz(&encoder, key)); @@ -679,7 +696,27 @@ void CBORStream2Serializer::SerializeSequenceEnd(const EndMessage& message) { CBOR_ENC(mapEncoder, "rotation_lattice_type", message.rotation_lattice_type); if (message.rotation_lattice.has_value()) CBOR_ENC(mapEncoder, "rotation_lattice", message.rotation_lattice->GetVector()); - CBOR_ENC(mapEncoder, "image_scale_factor", message.scale_factor); + + CBOR_ENC(mapEncoder, "data_collection_efficiency_image", message.data_collection_efficiency); + CBOR_ENC(mapEncoder, "spot_count", message.spot_count); + CBOR_ENC(mapEncoder, "spot_count_ice_ring", message.spot_count_ice_ring); + CBOR_ENC(mapEncoder, "spot_count_low_res", message.spot_count_low_res); + CBOR_ENC(mapEncoder, "spot_count_indexed", message.spot_count_indexed); + CBOR_ENC(mapEncoder, "image_indexed", message.image_indexed); + CBOR_ENC(mapEncoder, "v_bkg_estimate", message.v_bkg_estimate); + CBOR_ENC(mapEncoder, "profile_radius", message.profile_radius); + CBOR_ENC(mapEncoder, "mosaicity", message.mosaicity); + CBOR_ENC(mapEncoder, "bFactor", message.bFactor); + CBOR_ENC(mapEncoder, "resolution_estimate", message.resolution_estimate); + CBOR_ENC(mapEncoder, "min_viable_pixel_value", message.min_viable_pixel_value); + CBOR_ENC(mapEncoder, "max_viable_pixel_value", message.max_viable_pixel_value); + CBOR_ENC(mapEncoder, "saturated_pixel_count", message.saturated_pixel_count); + CBOR_ENC(mapEncoder, "error_pixel_count", message.error_pixel_count); + CBOR_ENC(mapEncoder, "image_scale_factor", message.image_scale_factor); + CBOR_ENC(mapEncoder, "integrated_reflections", message.integrated_reflections); + CBOR_ENC(mapEncoder, "niggli_class", message.niggli_class); + CBOR_ENC(mapEncoder, "pixel_sum", message.pixel_sum); + cborErr(cbor_encoder_close_container(&encoder, &mapEncoder)); curr_size = cbor_encoder_get_buffer_size(&encoder, buffer); diff --git a/frame_serialize/CborUtil.h b/frame_serialize/CborUtil.h index 8d56bec7..2f2f3ea6 100644 --- a/frame_serialize/CborUtil.h +++ b/frame_serialize/CborUtil.h @@ -24,5 +24,6 @@ constexpr const CborTag TagUnsignedInt32BitLE = 0b01000110; constexpr const CborTag TagSignedInt32BitLE = 0b01001110; constexpr const CborTag TagUnsignedInt64BitLE = 0b01000111; +constexpr const CborTag TagSignedInt64BitLE = 0b01001111; #endif //JUNGFRAUJOCH_CBORUTIL_H diff --git a/image_pusher/CBORFilePusher.cpp b/image_pusher/CBORFilePusher.cpp index 5e39c481..f09a761e 100644 --- a/image_pusher/CBORFilePusher.cpp +++ b/image_pusher/CBORFilePusher.cpp @@ -11,7 +11,7 @@ void CBORFilePusher::StartDataCollection(StartMessage &message) { dataset_number++; img_number = 0; - size_t approx_size = 1024*1024; + size_t approx_size = 256*1024*1024; for (const auto &[key, value] : message.pixel_mask) approx_size += value.size() * sizeof(uint32_t); @@ -25,7 +25,7 @@ void CBORFilePusher::StartDataCollection(StartMessage &message) { bool CBORFilePusher::EndDataCollection(const EndMessage &message) { std::unique_lock ul(m); - size_t approx_size = 1024*1024; + size_t approx_size = 256*1024*1024; std::vector serialization_buffer(approx_size); CBORStream2Serializer serializer(serialization_buffer.data(), serialization_buffer.size()); diff --git a/reader/JFJochHDF5Reader.cpp b/reader/JFJochHDF5Reader.cpp index 46bace0b..4f4d3209 100644 --- a/reader/JFJochHDF5Reader.cpp +++ b/reader/JFJochHDF5Reader.cpp @@ -1,6 +1,7 @@ // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include #include "JFJochHDF5Reader.h" #include "spdlog/fmt/fmt.h" #include "../image_analysis/bragg_integration/CalcISigma.h" @@ -55,7 +56,7 @@ inline std::pair parse_niggli_class(int64_t val) { } } -std::vector GetDimension(HDF5Object& object, const std::string& path) { +std::vector GetDimension(HDF5Object &object, const std::string &path) { const auto dim = object.GetDimension(path); if (dim.size() != 3) throw JFJochException(JFJochExceptionCategory::HDF5, "Wrong dimension of /entry/data/data"); @@ -93,7 +94,7 @@ std::string ResolveRelativeToMaster(const std::string &directory, return (std::filesystem::path(directory) / path).string(); } -template +template void JFJochHDF5Reader::ReadVector(std::vector &v, HDF5Object &file, const std::string &dataset_name, @@ -106,33 +107,125 @@ void JFJochHDF5Reader::ReadVector(std::vector &v, for (int i = 0; i < tmp.size(); i++) v[image0 + i] = tmp[i]; } - } catch (JFJochException &e) {} + } catch (JFJochException &e) { + } } -std::string removeSuffix(const std::string& s, const std::string& suffix) -{ +std::string removeSuffix(const std::string &s, const std::string &suffix) { if (s.ends_with(suffix)) return s.substr(0, s.size() - suffix.size()); return s; } -std::string dataset_name(const std::string& path) { +std::string dataset_name(const std::string &path) { std::string file = path; int pos = file.rfind('/'); if (pos != std::string::npos) - file = file.substr(pos+1); + file = file.substr(pos + 1); file = removeSuffix(file, "_master.h5"); -// If previous suffix was not found, try removing this one + // If previous suffix was not found, try removing this one file = removeSuffix(file, ".h5"); return file; } -void JFJochHDF5Reader::ReadFile(const std::string& filename) { +bool ReadReflectionsFromGroup(HDF5Object &file, + const std::string &image_group_name, + DataMessage &message) { + if (!file.Exists("/entry/reflections") || !file.Exists(image_group_name)) + return false; + + auto h = file.ReadOptVector(image_group_name + "/h"); + auto k = file.ReadOptVector(image_group_name + "/k"); + auto l = file.ReadOptVector(image_group_name + "/l"); + auto predicted_x = file.ReadOptVector(image_group_name + "/predicted_x"); + auto predicted_y = file.ReadOptVector(image_group_name + "/predicted_y"); + + auto d = file.ReadOptVector(image_group_name + "/d"); + auto int_sum = file.ReadOptVector(image_group_name + "/int_sum"); + auto int_err = file.ReadOptVector(image_group_name + "/int_err"); + auto bkg = file.ReadOptVector(image_group_name + "/background_mean"); + auto lp = file.ReadOptVector(image_group_name + "/lp"); + auto partiality = file.ReadOptVector(image_group_name + "/partiality"); + auto phi = file.ReadOptVector(image_group_name + "/delta_phi"); + auto zeta = file.ReadOptVector(image_group_name + "/zeta"); + + if (h.size() != l.size() || h.size() != k.size() || h.size() != d.size() + || h.size() != predicted_x.size() || h.size() != predicted_y.size() + || h.size() != int_sum.size() || h.size() != int_err.size() || h.size() != bkg.size()) + throw JFJochException(JFJochExceptionCategory::HDF5, "Wrong size of reflections dataset"); + + for (size_t i = 0; i < h.size(); i++) { + float lp_val = 0.0; + if (lp.size() > i && lp[i] != 0.0f) + lp_val = 1.0f / lp[i]; + + float partiality_val = -1.0f; + if (partiality.size() > i && partiality[i] >= 0.0f) + partiality_val = partiality[i]; + float delta_phi_val = NAN; + if (phi.size() > i) + delta_phi_val = phi[i]; + float zeta_val = NAN; + if (zeta.size() > i) + zeta_val = zeta[i]; + + Reflection r{ + .h = h.at(i), + .k = k.at(i), + .l = l.at(i), + .delta_phi_deg = delta_phi_val, + .predicted_x = predicted_x.at(i), + .predicted_y = predicted_y.at(i), + .d = d.at(i), + .I = int_sum.at(i), + .bkg = bkg.at(i), + .sigma = int_err.at(i), + .rlp = lp_val, + .partiality = partiality_val, + .zeta = zeta_val + }; + message.reflections.emplace_back(r); + } + + CalcISigma(message); + CalcWilsonBFactor(message, !message.b_factor.has_value()); + return true; +} + +template +std::optional ReadElementMasterFirst(HDF5Object &master_file, + HDF5Object &source_file, + const std::string &path, + hsize_t master_image, + hsize_t source_image) { + if (master_file.Exists(path)) + return master_file.ReadElement(path, master_image); + if (source_file.Exists(path)) + return source_file.ReadElement(path, source_image); + return {}; +} + +template +std::vector ReadVectorMasterFirst(HDF5Object &master_file, + HDF5Object &source_file, + const std::string &path, + const std::vector &master_start, + const std::vector &source_start, + const std::vector &size) { + if (master_file.Exists(path)) + return master_file.ReadOptVector(path, master_start, size); + if (source_file.Exists(path)) + return source_file.ReadOptVector(path, source_start, size); + return {}; +} + +void JFJochHDF5Reader::ReadFile(const std::string &filename) { std::unique_lock ul(hdf5_mutex); try { auto dataset = std::make_shared(); master_file = std::make_shared(filename); + master_filename = filename; dataset->experiment = default_experiment; std::filesystem::path master_path(filename); @@ -165,7 +258,7 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) { vds_data_mappings = ReadVDSImageMappings(*master_file, "/entry/data/data"); if (master_file->Exists("/entry/instrument/detector/detectorSpecific/data_collection_efficiency_image")) - dataset->efficiency = master_file->ReadVector( + dataset->efficiency = master_file->ReadVector( "/entry/instrument/detector/detectorSpecific/data_collection_efficiency_image"); else dataset->efficiency = std::vector(number_of_images, 1.0); @@ -198,6 +291,8 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) { dataset->profile_radius = master_file->ReadOptVector("/entry/MX/profileRadius"); dataset->mosaicity_deg = master_file->ReadOptVector("/entry/MX/mosaicity"); dataset->b_factor = master_file->ReadOptVector("/entry/MX/bFactor"); + dataset->scale_factor = master_file->ReadOptVector("/entry/MX/imageScaleFactor"); + dataset->integrated_reflections = master_file->ReadOptVector("/entry/MX/integratedReflections"); } if (master_file->Exists("/entry/image")) dataset->max_value = master_file->ReadOptVector("/entry/image/max_value"); @@ -317,8 +412,8 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) { number_of_images, fimages); ReadVector(dataset->mosaicity_deg, - data_file, "/entry/MX/mosaicity", - number_of_images, fimages); + data_file, "/entry/MX/mosaicity", + number_of_images, fimages); ReadVector(dataset->b_factor, data_file, "/entry/MX/bFactor", @@ -334,7 +429,8 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) { data_file, "/entry/image/max_value", number_of_images, fimages); } - } catch (JFJochException &e) {} + } catch (JFJochException &e) { + } if (nfiles == 0) images_per_file = fimages; @@ -353,17 +449,19 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) { dataset->experiment.IndexingAlgorithm(IndexingAlgorithmEnum::FFT); else if (indexing == "ffbidx") dataset->experiment.IndexingAlgorithm(IndexingAlgorithmEnum::FFBIDX); - - dataset->scale_factor = master_file->ReadOptVector("/entry/MX/imageScaleFactor"); } auto ring_current_A = master_file->GetOptFloat("/entry/source/current"); if (ring_current_A) dataset->experiment.RingCurrent_mA(ring_current_A.value() * 1000.0); - dataset->experiment.DetectIceRings(master_file->GetOptBool("/entry/instrument/detector/detectorSpecific/detect_ice_rings").value_or(false)); - dataset->experiment.PoniRot1_rad(master_file->GetOptFloat("/entry/instrument/detector/transformations/rot1").value_or(0.0)); - dataset->experiment.PoniRot2_rad(master_file->GetOptFloat("/entry/instrument/detector/transformations/rot2").value_or(0.0)); - dataset->experiment.PoniRot3_rad(master_file->GetOptFloat("/entry/instrument/detector/transformations/rot3").value_or(0.0)); + dataset->experiment.DetectIceRings( + master_file->GetOptBool("/entry/instrument/detector/detectorSpecific/detect_ice_rings").value_or(false)); + dataset->experiment.PoniRot1_rad( + master_file->GetOptFloat("/entry/instrument/detector/transformations/rot1").value_or(0.0)); + dataset->experiment.PoniRot2_rad( + master_file->GetOptFloat("/entry/instrument/detector/transformations/rot2").value_or(0.0)); + dataset->experiment.PoniRot3_rad( + master_file->GetOptFloat("/entry/instrument/detector/transformations/rot3").value_or(0.0)); dataset->experiment.SampleTemperature_K(master_file->GetOptFloat("/entry/sample/temperature")); dataset->experiment.BeamX_pxl(master_file->GetFloat("/entry/instrument/detector/beam_center_x")); @@ -374,7 +472,8 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) { det_distance = 0.1; // Set to 100 mm, if det distance is less than 1 mm dataset->experiment.DetectorDistance_mm(det_distance * 1000.0); - dataset->experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / master_file->GetFloat("/entry/instrument/beam/incident_wavelength")); + dataset->experiment.IncidentEnergy_keV( + WVL_1A_IN_KEV / master_file->GetFloat("/entry/instrument/beam/incident_wavelength")); dataset->error_value = master_file->GetOptInt("/entry/instrument/detector/error_value"); @@ -391,12 +490,12 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) { dataset->experiment.Goniometer(omega); } else if (master_file->Exists("/entry/sample/grid_scan")) { GridScanSettings grid( - master_file->GetInt("/entry/sample/grid_scan/n_fast"), - master_file->GetFloat("/entry/sample/grid_scan/step_x") * 1e6f, - master_file->GetFloat("/entry/sample/grid_scan/step_y") * 1e6f, - master_file->GetOptBool("/entry/sample/grid_scan/snake_scan").value_or(false), - master_file->GetOptBool("/entry/sample/grid_scan/vertical_scan").value_or(false) - ); + master_file->GetInt("/entry/sample/grid_scan/n_fast"), + master_file->GetFloat("/entry/sample/grid_scan/step_x") * 1e6f, + master_file->GetFloat("/entry/sample/grid_scan/step_y") * 1e6f, + master_file->GetOptBool("/entry/sample/grid_scan/snake_scan").value_or(false), + master_file->GetOptBool("/entry/sample/grid_scan/vertical_scan").value_or(false) + ); grid.ImageNum(number_of_images); dataset->experiment.GridScan(grid); } @@ -405,18 +504,20 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) { auto tmp = master_file->ReadOptVector("/entry/sample/unit_cell"); if (tmp.size() == 6) dataset->experiment.SetUnitCell(UnitCell{ - .a = tmp[0], - .b = tmp[1], - .c = tmp[2], - .alpha = tmp[3], - .beta = tmp[4], - .gamma = tmp[5]}); + .a = tmp[0], + .b = tmp[1], + .c = tmp[2], + .alpha = tmp[3], + .beta = tmp[4], + .gamma = tmp[5] + }); dataset->experiment.SpaceGroupNumber(master_file->GetOptInt("/entry/sample/space_group_number")); dataset->experiment.SampleName(master_file->GetString("/entry/sample/name")); if (master_file->Exists("/entry/instrument/attenuator")) - dataset->experiment.AttenuatorTransmission(master_file->GetOptFloat("/entry/instrument/attenuator/attenuator_transmission")); + dataset->experiment.AttenuatorTransmission( + master_file->GetOptFloat("/entry/instrument/attenuator/attenuator_transmission")); dataset->experiment.TotalFlux(master_file->GetOptFloat("/entry/instrument/beam/total_flux")); if (master_file->Exists("/entry/azint") && master_file->Exists("/entry/azint/bin_to_q")) { @@ -432,14 +533,14 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) { dataset->azimuthal_bins = dim[0]; dataset->q_bins = dim[1]; dataset->az_int_bin_to_q.resize(dim[0] * dim[1]); - bin_to_q_dataset.ReadVector(dataset->az_int_bin_to_q, {0,0}, dim); + bin_to_q_dataset.ReadVector(dataset->az_int_bin_to_q, {0, 0}, dim); } else throw JFJochException(JFJochExceptionCategory::HDF5, "Wrong dimension of /entry/azint/image dataset"); if (master_file->Exists("/entry/azint/bin_to_phi")) { HDF5DataSet bin_to_phi_dataset(*master_file, "/entry/azint/bin_to_phi"); if (dataset->q_bins > 0) { dataset->az_int_bin_to_phi.resize(dim[0] * dim[1]); - bin_to_phi_dataset.ReadVector(dataset->az_int_bin_to_phi, {0,0}, dim); + bin_to_phi_dataset.ReadVector(dataset->az_int_bin_to_phi, {0, 0}, dim); } else { bin_to_phi_dataset.ReadVector(dataset->az_int_bin_to_phi); } @@ -461,7 +562,7 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) { detector.SaturationLimit(master_file->GetInt("/entry/instrument/detector/saturation_value")); detector.MinFrameTime(std::chrono::microseconds(0)); detector.MinCountTime(std::chrono::microseconds(0)); - detector.ReadOutTime(std::chrono::nanoseconds (0)); + detector.ReadOutTime(std::chrono::nanoseconds(0)); dataset->experiment.Detector(detector); dataset->experiment.FrameTime( @@ -480,9 +581,9 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) { if (image_size_x * image_size_y > 0) { auto mask_tmp = master_file->ReadOptVector( - "/entry/instrument/detector/pixel_mask", - {0, 0}, - {image_size_y, image_size_x} + "/entry/instrument/detector/pixel_mask", + {0, 0}, + {image_size_y, image_size_x} ); if (mask_tmp.empty()) mask_tmp = std::vector(image_size_x * image_size_y); @@ -490,7 +591,7 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) { } dataset->experiment.ImagesPerTrigger(number_of_images); SetStartMessage(dataset); - } catch (const std::exception& e) { + } catch (const std::exception &e) { master_file = {}; number_of_images = 0; SetStartMessage({}); @@ -531,11 +632,14 @@ CompressedImage JFJochHDF5Reader::LoadImageDataset(std::vector &tmp, HD } if (datatype.IsFloat()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,"Float datasets not supported at this time"); + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Float datasets not supported at this time"); - return {tmp, dim[2], dim[1], - CalcImageMode(datatype.GetElemSize(), datatype.IsFloat(), datatype.IsSigned()), - algorithm}; + return { + tmp, dim[2], dim[1], + CalcImageMode(datatype.GetElemSize(), datatype.IsFloat(), datatype.IsSigned()), + algorithm + }; } std::pair, uint32_t> JFJochHDF5Reader::GetImageLocation(int64_t image_number) { @@ -558,7 +662,7 @@ std::pair, uint32_t> JFJochHDF5Reader::GetImag image_id = static_cast(mapping.SourceImage(image)); data_file = std::make_shared( - ResolveRelativeToMaster(master_file_directory, mapping.filename) + ResolveRelativeToMaster(master_file_directory, mapping.filename) ); return {std::move(data_file), image_id}; } @@ -601,69 +705,106 @@ bool JFJochHDF5Reader::LoadImage_i(std::shared_ptr &dataset "Cannot load image if file not loaded"); auto [source_file, image_id] = GetImageLocation(image_number); - + message.image = LoadImageDataset(buffer, *source_file, image_id); message.number = image_number; - auto spot_count_opt = source_file->ReadElement("/entry/MX/nPeaks", image_id); + const auto master_image = static_cast(image_number); + const auto source_image = static_cast(image_id); + + auto spot_count_opt = ReadElementMasterFirst(*master_file, + *source_file, + "/entry/MX/nPeaks", + master_image, + source_image); if (spot_count_opt.has_value() && spot_count_opt.value() > 0) { size_t spot_count = spot_count_opt.value(); - auto spot_x = source_file->ReadVector( - "/entry/MX/peakXPosRaw", - {(hsize_t) image_id, 0}, - {1, spot_count} + auto spot_x = ReadVectorMasterFirst( + *master_file, + *source_file, + "/entry/MX/peakXPosRaw", + {master_image, 0}, + {source_image, 0}, + {1, spot_count} ); - auto spot_y = source_file->ReadVector( - "/entry/MX/peakYPosRaw", - {(hsize_t) image_id, 0}, - {1, spot_count} + auto spot_y = ReadVectorMasterFirst( + *master_file, + *source_file, + "/entry/MX/peakYPosRaw", + {master_image, 0}, + {source_image, 0}, + {1, spot_count} ); - auto spot_intensity = source_file->ReadVector( - "/entry/MX/peakTotalIntensity", - {(hsize_t) image_id, 0}, - {1, spot_count} - ); - auto spot_indexed = source_file->ReadOptVector( - "/entry/MX/peakIndexed", - {(hsize_t) image_id, 0}, - {1, spot_count} + auto spot_intensity = ReadVectorMasterFirst( + *master_file, + *source_file, + "/entry/MX/peakTotalIntensity", + {master_image, 0}, + {source_image, 0}, + {1, spot_count} ); - auto spot_ice = source_file->ReadOptVector( - "/entry/MX/peakIceRingRes", - {(hsize_t) image_id, 0}, - {1, spot_count} + if (spot_x.size() < spot_count || spot_y.size() < spot_count || spot_intensity.size() < spot_count) + throw JFJochException(JFJochExceptionCategory::HDF5, "Wrong size of spot dataset"); + + auto spot_indexed = ReadVectorMasterFirst( + *master_file, + *source_file, + "/entry/MX/peakIndexed", + {master_image, 0}, + {source_image, 0}, + {1, spot_count} ); - auto spot_h = source_file->ReadOptVector( - "/entry/MX/peakH", - {(hsize_t) image_id, 0}, - {1, spot_count} + auto spot_ice = ReadVectorMasterFirst( + *master_file, + *source_file, + "/entry/MX/peakIceRingRes", + {master_image, 0}, + {source_image, 0}, + {1, spot_count} ); - auto spot_k = source_file->ReadOptVector( - "/entry/MX/peakK", - {(hsize_t) image_id, 0}, - {1, spot_count} + auto spot_h = ReadVectorMasterFirst( + *master_file, + *source_file, + "/entry/MX/peakH", + {master_image, 0}, + {source_image, 0}, + {1, spot_count} ); - auto spot_l = source_file->ReadOptVector( - "/entry/MX/peakL", - {(hsize_t) image_id, 0}, - {1, spot_count} + auto spot_k = ReadVectorMasterFirst( + *master_file, + *source_file, + "/entry/MX/peakK", + {master_image, 0}, + {source_image, 0}, + {1, spot_count} ); - auto spot_dist_ewald_sphere = source_file->ReadOptVector( - "/entry/MX/peakDistEwaldSphere", - {(hsize_t) image_id, 0}, - {1, spot_count} + auto spot_l = ReadVectorMasterFirst( + *master_file, + *source_file, + "/entry/MX/peakL", + {master_image, 0}, + {source_image, 0}, + {1, spot_count} + ); + + auto spot_dist_ewald_sphere = ReadVectorMasterFirst( + *master_file, + *source_file, + "/entry/MX/peakDistEwaldSphere", + {master_image, 0}, + {source_image, 0}, + {1, spot_count} ); auto geom = dataset->experiment.GetDiffractionGeometry(); for (int i = 0; i < spot_count; i++) { - auto x = spot_x.at(i); auto y = spot_y.at(i); auto d = geom.PxlToRes(x, y); @@ -689,40 +830,47 @@ bool JFJochHDF5Reader::LoadImage_i(std::shared_ptr &dataset s.ice_ring = (spot_ice.at(i) != 0); message.spots.emplace_back(s); } - if (source_file->Exists("/entry/MX/peakCountUnfiltered")) - message.spot_count = source_file->ReadElement("/entry/MX/peakCountUnfiltered", image_id); + + if (auto v = ReadElementMasterFirst(*master_file, *source_file, + "/entry/MX/peakCountUnfiltered", + master_image, source_image); v) + message.spot_count = v; else - message.spot_count = source_file->ReadElement("/entry/MX/nPeaks", image_id); - if (source_file->Exists("/entry/MX/peakCountIceRingRes")) - message.spot_count_ice_rings = source_file->ReadElement("/entry/MX/peakCountIceRingRes", image_id); - if (source_file->Exists("/entry/MX/peakCountLowRes")) - message.spot_count_low_res = source_file->ReadElement("/entry/MX/peakCountLowRes", image_id); - if (source_file->Exists("/entry/MX/peakCountIndexed")) - message.spot_count_indexed = source_file->ReadElement("/entry/MX/peakCountIndexed", image_id); + message.spot_count = spot_count_opt; + + message.spot_count_ice_rings = ReadElementMasterFirst( + *master_file, *source_file, "/entry/MX/peakCountIceRingRes", master_image, source_image); + message.spot_count_low_res = ReadElementMasterFirst( + *master_file, *source_file, "/entry/MX/peakCountLowRes", master_image, source_image); + message.spot_count_indexed = ReadElementMasterFirst( + *master_file, *source_file, "/entry/MX/peakCountIndexed", master_image, source_image); GenerateSpotPlot(message, 1.5); } - if (source_file->Exists("/entry/MX/integratedReflections")) - message.integrated_reflections = source_file->ReadElement("/entry/MX/integratedReflections", image_id); - if (!dataset->az_int_bin_to_q.empty()) { - if (dataset->azimuthal_bins == 0) - message.az_int_profile = source_file->ReadOptVector( + if (dataset->azimuthal_bins == 0) { + message.az_int_profile = ReadVectorMasterFirst( + *master_file, + *source_file, "/entry/azint/image", - {(hsize_t) image_id, 0}, + {master_image, 0}, + {source_image, 0}, {1, dataset->az_int_bin_to_q.size()} ); - else { - message.az_int_profile.resize(dataset->azimuthal_bins * dataset->q_bins, 0); - message.az_int_profile = source_file->ReadOptVector( + } else { + message.az_int_profile = ReadVectorMasterFirst( + *master_file, + *source_file, "/entry/azint/image", - {(hsize_t) image_id, 0, 0}, + {master_image, 0, 0}, + {source_image, 0, 0}, {1, dataset->azimuthal_bins, dataset->q_bins} ); } } - + if (dataset->integrated_reflections.size() > image_number) + message.integrated_reflections = static_cast(std::lround(dataset->integrated_reflections.at(image_number))); if (dataset->resolution_estimate.size() > image_number) message.resolution_estimate = dataset->resolution_estimate[image_number]; if (dataset->indexing_result.size() > image_number) @@ -740,17 +888,30 @@ bool JFJochHDF5Reader::LoadImage_i(std::shared_ptr &dataset if (dataset->indexing_result.size() > image_number && dataset->indexing_result[image_number] != 0 - && source_file->Exists("/entry/MX") - && source_file->Exists("/entry/MX/latticeIndexed")) { + && (master_file->Exists("/entry/MX/latticeIndexed") || + source_file->Exists("/entry/MX/latticeIndexed"))) { + std::vector tmp = ReadVectorMasterFirst( + *master_file, + *source_file, + "/entry/MX/latticeIndexed", + {master_image, 0}, + {source_image, 0}, + {1, 9} + ); - std::vector tmp = source_file->ReadVector( - "/entry/MX/latticeIndexed", - {(hsize_t) image_id, 0}, - {1, 9} - ); - message.indexing_lattice = CrystalLattice(tmp); + if (tmp.size() == 9) + message.indexing_lattice = CrystalLattice(tmp); + + std::optional niggli_opt; + if (master_file->Exists("/entry/MX/niggli_class")) + niggli_opt = master_file->ReadElement("/entry/MX/niggli_class", image_number); + else if (master_file->Exists("/entry/MX/niggliClass")) + niggli_opt = master_file->ReadElement("/entry/MX/niggliClass", image_number); + else if (source_file->Exists("/entry/MX/niggli_class")) + niggli_opt = source_file->ReadElement("/entry/MX/niggli_class", image_id); + else if (source_file->Exists("/entry/MX/niggliClass")) + niggli_opt = source_file->ReadElement("/entry/MX/niggliClass", image_id); - std::optional niggli_opt = source_file->ReadElement("/entry/MX/niggli_class", image_id); if (niggli_opt) { auto symm_info = parse_niggli_class(niggli_opt.value()); @@ -762,65 +923,11 @@ bool JFJochHDF5Reader::LoadImage_i(std::shared_ptr &dataset } } - std::string image_group_name = fmt::format("/entry/reflections/image_{:06d}", image_id); + const std::string master_reflection_group_name = fmt::format("/entry/reflections/image_{:06d}", image_number); + const std::string source_reflection_group_name = fmt::format("/entry/reflections/image_{:06d}", image_id); - if (source_file->Exists("/entry/reflections") && source_file->Exists(image_group_name)) { - auto h = source_file->ReadOptVector(image_group_name + "/h"); - auto k = source_file->ReadOptVector(image_group_name + "/k"); - auto l = source_file->ReadOptVector(image_group_name + "/l"); - auto predicted_x = source_file->ReadOptVector(image_group_name + "/predicted_x"); - auto predicted_y = source_file->ReadOptVector(image_group_name + "/predicted_y"); - - auto d = source_file->ReadOptVector(image_group_name + "/d"); - auto int_sum = source_file->ReadOptVector(image_group_name + "/int_sum"); - auto int_err = source_file->ReadOptVector(image_group_name + "/int_err"); - auto bkg = source_file->ReadOptVector(image_group_name + "/background_mean"); - auto lp = source_file->ReadOptVector(image_group_name + "/lp"); - auto partiality = source_file->ReadOptVector(image_group_name + "/partiality"); - auto phi = source_file->ReadOptVector(image_group_name + "/delta_phi"); - auto zeta = source_file->ReadOptVector(image_group_name + "/zeta"); - - if (h.size() != l.size() || h.size() != k.size() || h.size() != d.size() - || h.size() != predicted_x.size() || h.size() != predicted_y.size() - || h.size() != int_sum.size() || h.size() != int_err.size() || h.size() != bkg.size()) - throw JFJochException(JFJochExceptionCategory::HDF5, "Wrong size of reflections dataset"); - - for (size_t i = 0; i < h.size(); i++) { - float lp_val = 0.0; - if (lp.size() > i && lp[i] != 0.0f) - lp_val = 1.0f / lp[i]; - - float partiality_val = -1.0f; - if (partiality.size() > i && partiality[i] >= 0.0f) - partiality_val = partiality[i]; - float delta_phi_val = NAN; - if (phi.size() > i) - delta_phi_val = phi[i]; - float zeta_val = NAN; - if (zeta.size() > i) - zeta_val = zeta[i]; - - Reflection r{ - .h = h.at(i), - .k = k.at(i), - .l = l.at(i), - .delta_phi_deg = delta_phi_val, - .predicted_x = predicted_x.at(i), - .predicted_y = predicted_y.at(i), - .d = d.at(i), - .I = int_sum.at(i), - .bkg = bkg.at(i), - .sigma = int_err.at(i), - .rlp = lp_val, - .partiality = partiality_val, - .zeta = zeta_val - }; - message.reflections.emplace_back(r); - } - - CalcISigma(message); - CalcWilsonBFactor(message, !message.b_factor.has_value()); - } + if (!ReadReflectionsFromGroup(*master_file, master_reflection_group_name, message)) + ReadReflectionsFromGroup(*source_file, source_reflection_group_name, message); return true; } @@ -832,6 +939,7 @@ void JFJochHDF5Reader::Close() { legacy_format_files.clear(); vds_data_mappings.clear(); master_file_directory.clear(); + master_filename.clear(); SetStartMessage({}); } @@ -892,7 +1000,148 @@ CompressedImage JFJochHDF5Reader::ReadCalibration(std::vector &tmp, con dataset.ReadVectorToU8(tmp, start, {dim[0], dim[1]}); algorithm = CompressionAlgorithm::NO_COMPRESSION; - return {tmp, dim[1], dim[0], - CalcImageMode(datatype.GetElemSize(), datatype.IsFloat(), datatype.IsSigned()), - algorithm}; + return { + tmp, dim[1], dim[0], + CalcImageMode(datatype.GetElemSize(), datatype.IsFloat(), datatype.IsSigned()), + algorithm + }; } + +template +std::vector ReadVectorMasterFirst(HDF5Object &master_file, + HDF5Object &source_file, + const std::string &path, + const std::vector &master_start, + const std::vector &source_start, + const std::vector &size) { + if (master_file.Exists(path)) + return master_file.ReadOptVector(path, master_start, size); + if (source_file.Exists(path)) + return source_file.ReadOptVector(path, source_start, size); + return {}; +} + +void AppendOrExtendSourceMapping(std::vector &ret, + const std::string &filename, + const std::string &dataset, + uint64_t source_first_image, + uint64_t virtual_first_image, + uint64_t image_count) { + if (image_count == 0) + return; + + if (!ret.empty()) { + auto &last = ret.back(); + if (last.filename == filename + && last.dataset == dataset + && 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; + return; + } + } + + ret.push_back(HDF5DataSourceMessage{ + .filename = filename, + .dataset = dataset, + .source_first_image = source_first_image, + .virtual_first_image = virtual_first_image, + .image_count = image_count + }); +}; + +std::vector JFJochHDF5Reader::GetHDF5DataSource( + uint64_t first_image, + std::optional image_count +) const { + std::unique_lock ul(hdf5_mutex); + + if (!master_file) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Cannot generate HDF5 source mapping if file not loaded"); + + if (first_image > number_of_images) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "First image outside dataset range"); + + const uint64_t requested_count = image_count.value_or(number_of_images - first_image); + if (first_image + requested_count > number_of_images) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Requested image range outside dataset range"); + + std::vector ret; + if (requested_count == 0) + return ret; + + // Integrated / contiguous source: link directly to original master file. + if (format == FileWriterFormat::NXmxVDS && data_layout != HDF5DataSetLayout::VIRTUAL) { + AppendOrExtendSourceMapping(ret, + master_filename, + "/entry/data/data", + first_image, + 0, + requested_count); + return ret; + } + + // VDS source: expand VDS mappings to original source files, not to the VDS master. + if (format == FileWriterFormat::NXmxVDS && data_layout == HDF5DataSetLayout::VIRTUAL) { + for (uint64_t local_image = 0; local_image < requested_count; ++local_image) { + const hsize_t virtual_image = first_image + local_image; + + bool found = false; + for (const auto &mapping: vds_data_mappings) { + if (!mapping.ContainsVirtualImage(virtual_image)) + continue; + + const uint64_t source_image = mapping.SourceImage(virtual_image); + const std::string filename = ResolveRelativeToMaster(master_file_directory, mapping.filename); + const std::string dataset = mapping.dataset.empty() ? "/entry/data/data" : mapping.dataset; + + AppendOrExtendSourceMapping(ret, + filename, + dataset, + source_image, + local_image, + 1); + found = true; + break; + } + + if (!found) + throw JFJochException(JFJochExceptionCategory::HDF5, + "Image not covered by /entry/data/data VDS mappings"); + } + + return ret; + } + + // Legacy source: link directly to linked data files. + if (format == FileWriterFormat::NXmxLegacy) { + if (images_per_file == 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Cannot generate HDF5 source mapping: images_per_file is zero"); + + for (uint64_t local_image = 0; local_image < requested_count; ++local_image) { + const uint64_t source_global_image = first_image + local_image; + const uint64_t file_id = source_global_image / images_per_file; + const uint64_t source_image = source_global_image % images_per_file; + + if (file_id >= legacy_format_files.size()) + throw JFJochException(JFJochExceptionCategory::HDF5, + "Legacy image source file missing"); + + AppendOrExtendSourceMapping(ret, + legacy_format_files.at(file_id), + "/entry/data/data", + source_image, + local_image, + 1); + } + + return ret; + } + + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Unsupported HDF5 file layout for source mapping"); +} \ No newline at end of file diff --git a/reader/JFJochHDF5Reader.h b/reader/JFJochHDF5Reader.h index 49e7e137..8f3781e4 100644 --- a/reader/JFJochHDF5Reader.h +++ b/reader/JFJochHDF5Reader.h @@ -16,6 +16,7 @@ class JFJochHDF5Reader : public JFJochReader { std::vector legacy_format_files; std::vector vds_data_mappings; std::string master_file_directory; + std::string master_filename; size_t images_per_file = 1; size_t number_of_images = 0; @@ -43,8 +44,14 @@ public: void ReadFile(const std::string& filename); uint64_t GetNumberOfImages() const override; + void Close() override; + std::vector GetHDF5DataSource( + uint64_t first_image = 0, + std::optional image_count = {} + ) const; + CompressedImage ReadCalibration(std::vector &tmp, const std::string &name) const; std::shared_ptr GetRawImage(int64_t image_number) override; diff --git a/reader/JFJochHttpReader.cpp b/reader/JFJochHttpReader.cpp index 08935fc5..1f453037 100644 --- a/reader/JFJochHttpReader.cpp +++ b/reader/JFJochHttpReader.cpp @@ -182,6 +182,7 @@ std::shared_ptr JFJochHttpReader::UpdateDataset_i() { dataset->b_factor = GetPlot_i("b_factor"); dataset->resolution_estimate = GetPlot_i("resolution_estimate"); dataset->efficiency = GetPlot_i("image_collection_efficiency"); + dataset->integrated_reflections = GetPlot_i("integrated_reflections"); if (msg->start_message->goniometer) dataset->experiment.Goniometer(msg->start_message->goniometer); diff --git a/reader/JFJochReaderDataset.h b/reader/JFJochReaderDataset.h index e13c100b..d3d48f15 100644 --- a/reader/JFJochReaderDataset.h +++ b/reader/JFJochReaderDataset.h @@ -40,6 +40,7 @@ struct JFJochReaderDataset { std::vector profile_radius; std::vector mosaicity_deg; std::vector b_factor; + std::vector integrated_reflections; std::vector scale_factor; diff --git a/receiver/JFJochReceiver.cpp b/receiver/JFJochReceiver.cpp index 40c5b5a7..46bbd53b 100644 --- a/receiver/JFJochReceiver.cpp +++ b/receiver/JFJochReceiver.cpp @@ -175,6 +175,8 @@ void JFJochReceiver::SendEndMessage() { for (int i = 0; i < adu_histogram_module.size(); i++) message.adu_histogram["module" + std::to_string(i)] = adu_histogram_module[i]->GetHistogram(); + scan_result.FillEndMessage(message); + if (push_images_to_writer) { if (!image_pusher.EndDataCollection(message)) logger.Error("End message not sent via ZeroMQ (time-out)"); diff --git a/tests/HDF5WritingTest.cpp b/tests/HDF5WritingTest.cpp index 1269b116..a1ec0525 100644 --- a/tests/HDF5WritingTest.cpp +++ b/tests/HDF5WritingTest.cpp @@ -468,7 +468,9 @@ TEST_CASE("HDF5Writer", "[HDF5][Full]") { DiffractionExperiment x(DetJF4M()); std::vector spots; - x.FilePrefix("test02_1p10").ImagesPerTrigger(5).ImagesPerFile(2).Compression(CompressionAlgorithm::NO_COMPRESSION); + x.FilePrefix("test02_1p10").ImagesPerTrigger(5).ImagesPerFile(2).Compression(CompressionAlgorithm::NO_COMPRESSION) + .OverwriteExistingFiles(true); + StartMessage start_message; x.FillMessage(start_message); diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index 3807a470..b08b008a 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -1670,4 +1670,278 @@ TEST_CASE("JFJochReader_GetRawImage_Integrated", "[HDF5][Full]") { 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 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(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 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(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 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(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 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(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 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); } \ No newline at end of file diff --git a/writer/HDF5DataFilePluginDetector.cpp b/writer/HDF5DataFilePluginDetector.cpp index a7ab1efc..775c5bab 100644 --- a/writer/HDF5DataFilePluginDetector.cpp +++ b/writer/HDF5DataFilePluginDetector.cpp @@ -15,8 +15,6 @@ void HDF5DataFilePluginDetector::OpenFile(HDF5File &in_data_file, const DataMess efficiency.reserve(images_per_file); packets_received.reserve(images_per_file); packets_expected.reserve(images_per_file); - pixel_sum.reserve(images_per_file); - processing_time.reserve(images_per_file); } void HDF5DataFilePluginDetector::Write(const DataMessage &msg, uint64_t image_number) { @@ -33,8 +31,6 @@ void HDF5DataFilePluginDetector::Write(const DataMessage &msg, uint64_t image_nu packets_received[image_number] = msg.packets_received.value(); if (msg.packets_expected.has_value()) packets_expected[image_number] = msg.packets_expected.value(); - if (msg.pixel_sum.has_value()) - pixel_sum[image_number] = msg.pixel_sum.value(); if (msg.image_collection_efficiency.has_value()) efficiency[image_number] = msg.image_collection_efficiency.value(); @@ -57,7 +53,5 @@ void HDF5DataFilePluginDetector::WriteFinal(HDF5File &data_file) { data_file.SaveVector(prefix + "/packets_received", packets_received.vec()); if (!packets_expected.empty()) data_file.SaveVector(prefix + "/packets_expected", packets_expected.vec()); - if (!pixel_sum.empty()) - data_file.SaveVector(prefix + "/pixel_sum", pixel_sum.vec()); data_file.SaveVector(prefix + "/data_collection_efficiency_image", efficiency.vec()); } diff --git a/writer/HDF5DataFilePluginDetector.h b/writer/HDF5DataFilePluginDetector.h index c61d7d51..99a8f6b6 100644 --- a/writer/HDF5DataFilePluginDetector.h +++ b/writer/HDF5DataFilePluginDetector.h @@ -16,8 +16,6 @@ class HDF5DataFilePluginDetector : public HDF5DataFilePlugin { AutoIncrVector efficiency; AutoIncrVector packets_received; AutoIncrVector packets_expected; - AutoIncrVector pixel_sum; - AutoIncrVector processing_time{NAN}; public: HDF5DataFilePluginDetector(const StartMessage& msg); void OpenFile(HDF5File &data_file, const DataMessage& msg, size_t images_per_file) override; diff --git a/writer/HDF5DataFilePluginImageStats.cpp b/writer/HDF5DataFilePluginImageStats.cpp index 74c8982e..ac2cd790 100644 --- a/writer/HDF5DataFilePluginImageStats.cpp +++ b/writer/HDF5DataFilePluginImageStats.cpp @@ -8,6 +8,7 @@ void HDF5DataFilePluginImageStats::OpenFile(HDF5File &data_file, const DataMessa min_value.reserve(images_per_file); error_pixels.reserve(images_per_file); saturated_pixels.reserve(images_per_file); + pixel_sum.reserve(images_per_file); } void HDF5DataFilePluginImageStats::Write(const DataMessage &msg, uint64_t image_number) { @@ -19,6 +20,8 @@ void HDF5DataFilePluginImageStats::Write(const DataMessage &msg, uint64_t image_ error_pixels[image_number] = msg.error_pixel_count.value(); if (msg.saturated_pixel_count) saturated_pixels[image_number] = msg.saturated_pixel_count.value(); + if (msg.pixel_sum) + pixel_sum[image_number] = msg.pixel_sum.value(); } void HDF5DataFilePluginImageStats::WriteFinal(HDF5File &data_file) { @@ -31,4 +34,6 @@ void HDF5DataFilePluginImageStats::WriteFinal(HDF5File &data_file) { data_file.SaveVector("/entry/image/error_pixels", error_pixels.vec()); if (!saturated_pixels.empty()) data_file.SaveVector("/entry/image/saturated_pixels", saturated_pixels.vec()); + if (!pixel_sum.empty()) + data_file.SaveVector("/entry/image/pixel_sum", pixel_sum.vec()); } diff --git a/writer/HDF5DataFilePluginImageStats.h b/writer/HDF5DataFilePluginImageStats.h index 0d4dc87e..b23e0ae5 100644 --- a/writer/HDF5DataFilePluginImageStats.h +++ b/writer/HDF5DataFilePluginImageStats.h @@ -13,6 +13,7 @@ class HDF5DataFilePluginImageStats : public HDF5DataFilePlugin { AutoIncrVector min_value; AutoIncrVector error_pixels; AutoIncrVector saturated_pixels; + AutoIncrVector pixel_sum; public: void OpenFile(HDF5File &data_file, const DataMessage& msg, size_t images_per_file) override; void Write(const DataMessage& msg, uint64_t image_number) override; diff --git a/writer/HDF5DataFilePluginMX.cpp b/writer/HDF5DataFilePluginMX.cpp index fa3951bc..e8f372a3 100644 --- a/writer/HDF5DataFilePluginMX.cpp +++ b/writer/HDF5DataFilePluginMX.cpp @@ -225,9 +225,9 @@ void HDF5DataFilePluginMX::WriteFinal(HDF5File &data_file) { if (!beam_corr_y.empty()) data_file.SaveVector("/entry/MX/beam_corr_y", beam_corr_y.vec())->Units("pixel"); if (!niggli_class.empty()) - data_file.SaveVector("/entry/MX/niggli_class", niggli_class.vec()); + data_file.SaveVector("/entry/MX/niggliClass", niggli_class.vec()); if (!bravais_lattice.empty()) - data_file.SaveVector("/entry/MX/bravais_lattice", bravais_lattice.vec()); + data_file.SaveVector("/entry/MX/bravaisLattice", bravais_lattice.vec()); if (!resolution_estimate.empty()) data_file.SaveVector("/entry/MX/resolutionEstimate", resolution_estimate.vec())->Units("Angstrom"); if (!integrated_reflections.empty()) diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index 29f2a757..fdee855c 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -21,7 +21,7 @@ NXmx::NXmx(const StartMessage &start) uint64_t tmp_suffix; try { if (!start.arm_date.empty()) - tmp_suffix = parse_UTC_to_ms(start.arm_date); + tmp_suffix = parse_UTC_to_ms(start.arm_date); } catch (...) { tmp_suffix = std::chrono::system_clock::now().time_since_epoch().count(); } @@ -33,7 +33,8 @@ NXmx::NXmx(const StartMessage &start) MakeDirectory(filename); - bool v1_10 = (start.file_format == FileWriterFormat::NXmxVDS); + bool v1_10 = (start.file_format == FileWriterFormat::NXmxVDS) + || !start.hdf5_source_data.empty(); hdf5_file = std::make_shared(tmp_filename, v1_10); hdf5_file->Attr("file_name", filename); @@ -119,54 +120,11 @@ void NXmx::LinkToData_VDS(const StartMessage &start, const EndMessage &end) { data_dataset->Attr("image_nr_low", (int32_t) 1) .Attr("image_nr_high",(int32_t) total_images); - VDS(start, - "/entry/detector/data_collection_efficiency_image", - "/entry/instrument/detector/detectorSpecific/data_collection_efficiency_image", - {total_images}, - HDF5DataType(0.0f)); - - VDS(start, - "/entry/detector/pixel_sum", - "/entry/instrument/detector/detectorSpecific/pixel_sum", - {total_images}, - HDF5DataType((int64_t) 0)); - - HDF5Group(*hdf5_file, "/entry/image").NXClass("NXCollection"); - - VDS(start, - "/entry/image/max_value", - {total_images}, - HDF5DataType((int64_t)0)); - - VDS(start, - "/entry/image/min_value", - {total_images}, - HDF5DataType((int64_t)0)); - - VDS(start, - "/entry/image/error_pixels", - {total_images}, - HDF5DataType((int64_t)0)); - - VDS(start, - "/entry/image/saturated_pixels", - {total_images}, - HDF5DataType((int64_t)0)); - if (start.max_spot_count > 0) { VDS(start, "/entry/MX/peakXPosRaw",{total_images, start.max_spot_count}, HDF5DataType(0.0f)); VDS(start, "/entry/MX/peakYPosRaw",{total_images, start.max_spot_count}, HDF5DataType(0.0f)); VDS(start, "/entry/MX/peakTotalIntensity",{total_images, start.max_spot_count}, HDF5DataType(0.0f)); VDS(start, "/entry/MX/peakIceRingRes", {total_images, start.max_spot_count}, HDF5DataType((uint8_t) 0)); - VDS(start, "/entry/MX/nPeaks",{total_images}, HDF5DataType((uint32_t) 0)); - VDS(start, "/entry/MX/strongPixels", {total_images}, HDF5DataType((uint32_t) 0)); - VDS(start, "/entry/MX/bkgEstimate", {total_images}, HDF5DataType(0.0f)); - VDS(start, "/entry/MX/resolutionEstimate",{total_images}, HDF5DataType(0.0f)); - - VDS(start, "/entry/MX/peakCountUnfiltered",{total_images}, HDF5DataType((uint32_t) 0)); - VDS(start, "/entry/MX/peakCountIceRingRes",{total_images}, HDF5DataType((uint32_t) 0)); - VDS(start, "/entry/MX/peakCountLowRes",{total_images}, HDF5DataType((uint32_t) 0)); - VDS(start, "/entry/MX/peakCountIndexed",{total_images}, HDF5DataType((uint32_t) 0)); } if (start.indexing_algorithm != IndexingAlgorithmEnum::None) { @@ -175,11 +133,7 @@ void NXmx::LinkToData_VDS(const StartMessage &start, const EndMessage &end) { VDS(start, "/entry/MX/peakK", {total_images, start.max_spot_count}, HDF5DataType((int32_t) 0)); VDS(start, "/entry/MX/peakL", {total_images, start.max_spot_count}, HDF5DataType((int32_t) 0)); VDS(start, "/entry/MX/peakDistEwaldSphere", {total_images, start.max_spot_count}, HDF5DataType((float) 0)); - - VDS(start, "/entry/MX/imageIndexed", {total_images}, HDF5DataType((uint8_t) 0)); VDS(start, "/entry/MX/latticeIndexed", {total_images,9}, HDF5DataType((float) 0))->Units("Angstrom"); - VDS(start, "/entry/MX/profileRadius", {total_images}, HDF5DataType(0.0f))->Units("Angstrom^-1"); - VDS(start, "/entry/MX/bFactor", {total_images}, HDF5DataType(0.0f))->Units("Angstrom^2"); } if (start.geom_refinement_algorithm != GeomRefinementAlgorithmEnum::None) { @@ -215,20 +169,89 @@ void NXmx::LinkToData_VDS(const StartMessage &start, const EndMessage &end) { {total_images}, HDF5DataType((uint8_t) 0)); - if (!start.rois.empty()) { - HDF5Group(*hdf5_file, "/entry/roi").NXClass("NXcollection"); + LinkToReflections_VDS(start, end); + } +} - for (const auto &r: start.rois) { - std::string roi = r.name; - HDF5Group(*hdf5_file, "/entry/roi/" + roi); - VDS(start, "/entry/roi/" + roi + "/max", {total_images}, HDF5DataType((int64_t) 0)); - VDS(start, "/entry/roi/" + roi + "/sum", {total_images}, HDF5DataType((int64_t) 0)); - VDS(start, "/entry/roi/" + roi + "/sum_sq", {total_images}, HDF5DataType((int64_t) 0)); - VDS(start, "/entry/roi/" + roi + "/npixel", {total_images}, HDF5DataType((int64_t) 0)); - VDS(start, "/entry/roi/" + roi + "/x", {total_images}, HDF5DataType((float) 0)); - VDS(start, "/entry/roi/" + roi + "/y", {total_images}, HDF5DataType((float) 0)); - } - } +void NXmx::LinkToData_ProcessingVDS(const StartMessage &start, const EndMessage &end) { + if (start.hdf5_source_data.empty() || end.max_image_number == 0) + return; + + const hsize_t total_images = end.max_image_number; + const hsize_t width = start.image_size_x; + const hsize_t height = start.image_size_y; + + HDF5Group(*hdf5_file, "/entry/data").NXClass("NXdata"); + + HDF5DataSpace full_data_space({total_images, height, width}); + HDF5Dcpl dcpl; + dcpl.SetChunking({1, height, width}); + + for (const auto &mapping: start.hdf5_source_data) { + if (mapping.image_count == 0) + continue; + + if (mapping.virtual_first_image + mapping.image_count > total_images) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Processing VDS mapping exceeds output image count"); + + const std::string source_dataset = mapping.dataset.empty() + ? "/entry/data/data" + : mapping.dataset; + + HDF5DataSpace virtual_data_space({total_images, height, width}); + virtual_data_space.SelectHyperslab( + {static_cast(mapping.virtual_first_image), 0, 0}, + {static_cast(mapping.image_count), height, width} + ); + + const hsize_t source_extent_images = mapping.source_first_image + mapping.image_count; + HDF5DataSpace source_data_space({source_extent_images, height, width}); + source_data_space.SelectHyperslab( + {static_cast(mapping.source_first_image), 0, 0}, + {static_cast(mapping.image_count), height, width} + ); + + dcpl.SetVirtual(mapping.filename, + source_dataset, + source_data_space, + virtual_data_space); + } + + auto data_dataset = std::make_unique( + *hdf5_file, + "/entry/data/data", + HDF5DataType(start.bit_depth_image / 8, start.pixel_signed), + full_data_space, + dcpl + ); + + data_dataset->Attr("image_nr_low", static_cast(1)) + .Attr("image_nr_high", static_cast(total_images)); +} + +void NXmx::LinkToReflections_VDS(const StartMessage &start, const EndMessage &end) { + if (end.integrated_reflections.empty()) + return; + + HDF5Group(*hdf5_file, "/entry/reflections").NXClass("NXcollection"); + + for (size_t image = 0; image < end.integrated_reflections.size(); ++image) { + if (end.integrated_reflections[image] <= 0) + continue; + + if (start.images_per_file <= 0) + continue; + + const uint64_t file_id = image / static_cast(start.images_per_file); + const uint64_t image_in_file = image % static_cast(start.images_per_file); + + const std::string local_name = fmt::format("/entry/reflections/image_{:06d}", image); + const std::string source_name = fmt::format("/entry/reflections/image_{:06d}", image_in_file); + + hdf5_file->ExternalLink(HDF5Metadata::DataFileName(start, file_id), + source_name, + local_name); } } @@ -696,6 +719,19 @@ void NXmx::ADUHistogram(const EndMessage &end) { } } +template +void SaveVectorIfMissing(HDF5Object &object, + const std::string &path, + const std::vector &values, + const std::string &units = "") { + if (values.empty() || object.Exists(path)) + return; + + auto dataset = object.SaveVector(path, values); + if (!units.empty()) + dataset->Units(units); +} + void NXmx::Finalize(const EndMessage &end) { try { if (!hdf5_file) @@ -715,6 +751,7 @@ void NXmx::Finalize(const EndMessage &end) { Sample(start_message, end); AzimuthalIntegration(start_message, end); ADUHistogram(end); + EndResultVectors(end); switch (start_message.file_format.value_or(FileWriterFormat::NXmxLegacy)) { case FileWriterFormat::NXmxLegacy: @@ -724,6 +761,9 @@ void NXmx::Finalize(const EndMessage &end) { LinkToData_VDS(start_message, end); break; case FileWriterFormat::NXmxIntegrated: + if (!start_message.hdf5_source_data.empty()) + LinkToData_ProcessingVDS(start_message, end); + break; default: break; } @@ -742,9 +782,6 @@ void NXmx::Finalize(const EndMessage &end) { SaveScalar(*hdf5_file, "/entry/MX/bkgEstimateMean", end.bkg_estimate.value()); } - if (!end.scale_factor.empty()) - SaveVector(*hdf5_file, "/entry/MX/imageScaleFactor", end.scale_factor); - hdf5_file->Close(); hdf5_file.reset(); } catch (const JFJochException &e) { @@ -784,3 +821,49 @@ void NXmx::UserData(const StartMessage &start) { std::shared_ptr NXmx::GetFile() { return hdf5_file; } + +void NXmx::EndResultVectors(const EndMessage &end) { + if (!end.data_collection_efficiency.empty()) { + HDF5Group det_specific(*hdf5_file, "/entry/instrument/detector/detectorSpecific"); + det_specific.NXClass("NXcollection"); + SaveVectorIfMissing(*hdf5_file, + "/entry/instrument/detector/detectorSpecific/data_collection_efficiency_image", + end.data_collection_efficiency); + } + + if (!end.max_viable_pixel_value.empty() || + !end.min_viable_pixel_value.empty() || + !end.error_pixel_count.empty() || + !end.saturated_pixel_count.empty() || + !end.pixel_sum.empty()) { + HDF5Group image_group(*hdf5_file, "/entry/image"); + image_group.NXClass("NXcollection"); + + SaveVectorIfMissing(*hdf5_file, "/entry/image/max_value", end.max_viable_pixel_value); + SaveVectorIfMissing(*hdf5_file, "/entry/image/min_value", end.min_viable_pixel_value); + SaveVectorIfMissing(*hdf5_file, "/entry/image/error_pixels", end.error_pixel_count); + SaveVectorIfMissing(*hdf5_file, "/entry/image/saturated_pixels", end.saturated_pixel_count); + SaveVectorIfMissing(*hdf5_file, "/entry/image/pixel_sum", end.pixel_sum); + } + + HDF5Group mx_group(*hdf5_file, "/entry/MX"); + mx_group.NXClass("NXcollection"); + + SaveVectorIfMissing(*hdf5_file, "/entry/MX/nPeaks", end.spot_count); + SaveVectorIfMissing(*hdf5_file, "/entry/MX/peakCountIceRingRes", end.spot_count_ice_ring); + SaveVectorIfMissing(*hdf5_file, "/entry/MX/peakCountLowRes", end.spot_count_low_res); + SaveVectorIfMissing(*hdf5_file, "/entry/MX/peakCountIndexed", end.spot_count_indexed); + SaveVectorIfMissing(*hdf5_file, "/entry/MX/imageIndexed", end.image_indexed); + SaveVectorIfMissing(*hdf5_file, "/entry/MX/bkgEstimate", end.v_bkg_estimate); + SaveVectorIfMissing(*hdf5_file, "/entry/MX/profileRadius", end.profile_radius, "Angstrom^-1"); + SaveVectorIfMissing(*hdf5_file, "/entry/MX/mosaicity", end.mosaicity, "deg"); + SaveVectorIfMissing(*hdf5_file, "/entry/MX/bFactor", end.bFactor, "Angstrom^2"); + SaveVectorIfMissing(*hdf5_file, "/entry/MX/resolutionEstimate", end.resolution_estimate, "Angstrom"); + SaveVectorIfMissing(*hdf5_file, "/entry/MX/imageScaleFactor", end.image_scale_factor); + SaveVectorIfMissing(*hdf5_file, "/entry/MX/integratedReflections", end.integrated_reflections); + + if (!end.niggli_class.empty()) { + SaveVectorIfMissing(*hdf5_file, "/entry/MX/niggliClass", end.niggli_class); + SaveVectorIfMissing(*hdf5_file, "/entry/MX/niggli_class", end.niggli_class); + } +} diff --git a/writer/HDF5NXmx.h b/writer/HDF5NXmx.h index 8bc0b36d..30612dc0 100644 --- a/writer/HDF5NXmx.h +++ b/writer/HDF5NXmx.h @@ -1,8 +1,7 @@ // SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only -#ifndef JUNGFRAUJOCH_HDF5NXMX_H -#define JUNGFRAUJOCH_HDF5NXMX_H +#pragma once #include "../common/JFJochMessages.h" @@ -23,6 +22,9 @@ class NXmx { void LinkToData(const StartMessage &start, const EndMessage &end); void LinkToData_VDS(const StartMessage &start, const EndMessage &end); + void LinkToData_ProcessingVDS(const StartMessage &start, const EndMessage &end); + void LinkToReflections_VDS(const StartMessage &start, const EndMessage &end); + std::unique_ptr VDS(const StartMessage &start, const std::string& name, const std::vector &dim, @@ -49,9 +51,11 @@ class NXmx { void SaveCBORImage(const std::string& hdf5_path, const CompressedImage &image); void AzimuthalIntegration(const StartMessage &start, const EndMessage &end); void ADUHistogram(const EndMessage &end); + void EndResultVectors(const EndMessage &end); void UserData(const StartMessage &start); void MX(const StartMessage &start); void Fluorescence(const StartMessage &start); + public: NXmx(const StartMessage& start); ~NXmx(); @@ -62,5 +66,3 @@ public: std::shared_ptr GetFile(); }; - -#endif //JUNGFRAUJOCH_HDF5NXMX_H