Compare commits
24 Commits
main
..
2605-scaling
| Author | SHA1 | Date | |
|---|---|---|---|
 leonarski_f
|
4b6a042e4a | ||
|
leonarski_f
|
b9af590ff5 | ||
|
leonarski_f
|
a441e96b50 | ||
|
leonarski_f
|
a90ec13418 | ||
|
leonarski_f
|
81af116b4d | ||
|
leonarski_f
|
3e4e6019c4 | ||
|
leonarski_f
|
e392a3ae1b | ||
|
leonarski_f
|
cc42ae6bf6 | ||
|
leonarski_f
|
cf2ca90fb7 | ||
|
leonarski_f
|
38c2826c09 | ||
|
leonarski_f
|
d718999335 | ||
|
leonarski_f
|
93914e1fb9 | ||
|
leonarski_f
|
0503ef0271 | ||
|
leonarski_f
|
1b8c528ed2 | ||
|
leonarski_f
|
8b5eb5a208 | ||
|
leonarski_f
|
da15714080 | ||
|
leonarski_f
|
91c714f1a4 | ||
|
leonarski_f
|
025c9b3aee | ||
|
leonarski_f
|
27bcb19328 | ||
|
leonarski_f
|
f5176b56a9 | ||
|
leonarski_f
|
cd5d97aa55 | ||
|
leonarski_f
|
6c8c953c92 | ||
|
leonarski_f
|
173198be40 | ||
|
leonarski_f
|
930cfb0b35 |
@@ -17,6 +17,8 @@ constexpr size_t CONVERTED_MODULE_COLS = 1030;
|
||||
constexpr size_t CONVERTED_MODULE_SIZE = CONVERTED_MODULE_LINES * CONVERTED_MODULE_COLS;
|
||||
constexpr size_t JUNGFRAU_PACKET_SIZE_BYTES = 8192;
|
||||
|
||||
constexpr int MAX_IMAGE_NUMBER = 2*1024*1024;
|
||||
|
||||
constexpr std::chrono::nanoseconds MIN_COUNT_TIME = std::chrono::microseconds(3);
|
||||
constexpr std::chrono::nanoseconds MIN_STORAGE_CELL_DELAY = std::chrono::nanoseconds(2100);
|
||||
constexpr std::chrono::nanoseconds MIN_FRAME_TIME_JUNGFRAU_HALF_SPEED = std::chrono::microseconds(1000);
|
||||
|
||||
@@ -1060,6 +1060,12 @@ DiffractionExperiment &DiffractionExperiment::ImportDatasetSettings(const Datase
|
||||
+ std::to_string(MAX_FRAMES));
|
||||
}
|
||||
|
||||
if (GetImageNum() > MAX_IMAGE_NUMBER) {
|
||||
dataset = tmp;
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"Number of images cannot exceed " + std::to_string(MAX_IMAGE_NUMBER));
|
||||
}
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
|
||||
+34
-2
@@ -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<LatticeMessage> 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;
|
||||
@@ -251,6 +260,7 @@ struct StartMessage {
|
||||
std::optional<float> attenuator_transmission;
|
||||
|
||||
std::optional<bool> write_master_file;
|
||||
std::optional<bool> write_images;
|
||||
|
||||
nlohmann::json user_data;
|
||||
|
||||
@@ -282,6 +292,9 @@ struct StartMessage {
|
||||
XrayFluorescenceSpectrum fluorescence_spectrum;
|
||||
|
||||
std::optional<bool> detect_ice_rings;
|
||||
|
||||
std::vector<HDF5DataSourceMessage> hdf5_source_data;
|
||||
std::optional<std::string> master_suffix;
|
||||
};
|
||||
|
||||
struct EndMessage {
|
||||
@@ -305,7 +318,26 @@ struct EndMessage {
|
||||
std::optional<LatticeMessage> rotation_lattice_type;
|
||||
std::optional<CrystalLattice> rotation_lattice;
|
||||
|
||||
std::vector<float> scale_factor;
|
||||
// Vectors with end result:
|
||||
std::vector<float> data_collection_efficiency;
|
||||
std::vector<int32_t> spot_count;
|
||||
std::vector<int32_t> spot_count_ice_ring;
|
||||
std::vector<int32_t> spot_count_low_res;
|
||||
std::vector<int32_t> spot_count_indexed;
|
||||
std::vector<uint8_t> image_indexed;
|
||||
std::vector<float> v_bkg_estimate;
|
||||
std::vector<float> profile_radius;
|
||||
std::vector<float> mosaicity;
|
||||
std::vector<float> bFactor;
|
||||
std::vector<float> resolution_estimate;
|
||||
std::vector<int64_t> min_viable_pixel_value;
|
||||
std::vector<int64_t> max_viable_pixel_value;
|
||||
std::vector<int32_t> saturated_pixel_count;
|
||||
std::vector<int32_t> error_pixel_count;
|
||||
std::vector<float> image_scale_factor;
|
||||
std::vector<int32_t> integrated_reflections;
|
||||
std::vector<uint8_t> niggli_class;
|
||||
std::vector<int64_t> pixel_sum;
|
||||
};
|
||||
|
||||
struct MetadataMessage {
|
||||
|
||||
@@ -24,7 +24,17 @@ struct Reflection {
|
||||
float rlp;
|
||||
float partiality;
|
||||
float zeta;
|
||||
float scaling_correction; // I_true = scaling_correction * I; scaling_correction = rlp / (partiality * image_scale)
|
||||
bool observed = false;
|
||||
};
|
||||
|
||||
struct MergedReflection {
|
||||
int32_t h;
|
||||
int32_t k;
|
||||
int32_t l;
|
||||
float I;
|
||||
float sigma;
|
||||
float d = 0.0;
|
||||
};
|
||||
|
||||
#endif //JFJOCH_REFLECTION_H
|
||||
|
||||
@@ -22,6 +22,7 @@ struct ScanResultElem {
|
||||
std::optional<float> angle_deg;
|
||||
|
||||
std::optional<uint64_t> pixel_sum;
|
||||
std::optional<int64_t> min_viable_pixel;
|
||||
std::optional<int64_t> max_viable_pixel;
|
||||
std::optional<int64_t> err_pixels;
|
||||
std::optional<int64_t> sat_pixels;
|
||||
@@ -36,6 +37,10 @@ struct ScanResultElem {
|
||||
std::optional<float> res;
|
||||
std::optional<UnitCell> uc;
|
||||
std::optional<uint64_t> xfel_pulse_id;
|
||||
std::optional<float> mosaicity;
|
||||
std::optional<float> image_scale_factor;
|
||||
std::optional<int64_t> niggli_class;
|
||||
std::optional<int64_t> integrated_reflections;
|
||||
};
|
||||
|
||||
struct ScanResult {
|
||||
|
||||
@@ -1,8 +1,18 @@
|
||||
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
|
||||
#include "ScanResultGenerator.h"
|
||||
|
||||
namespace {
|
||||
template<class T>
|
||||
T value_or_zero(const std::optional<T>& v) {
|
||||
return v.value_or(static_cast<T>(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<size_t>(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<size_t>(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<size_t>(e.number);
|
||||
if (number >= n)
|
||||
continue;
|
||||
|
||||
message.data_collection_efficiency[number] = e.collection_efficiency;
|
||||
message.spot_count[number] = static_cast<int32_t>(value_or_zero(e.spot_count));
|
||||
message.spot_count_ice_ring[number] = static_cast<int32_t>(value_or_zero(e.spot_count_ice));
|
||||
message.spot_count_low_res[number] = static_cast<int32_t>(value_or_zero(e.spot_count_low_res));
|
||||
message.spot_count_indexed[number] = static_cast<int32_t>(value_or_zero(e.spot_count_indexed));
|
||||
message.image_indexed[number] = static_cast<uint8_t>(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<int32_t>(value_or_zero(e.sat_pixels));
|
||||
message.error_pixel_count[number] = static_cast<int32_t>(value_or_zero(e.err_pixels));
|
||||
message.image_scale_factor[number] = e.image_scale_factor.value_or(NAN);
|
||||
message.integrated_reflections[number] = static_cast<int32_t>(value_or_zero(e.integrated_reflections));
|
||||
message.niggli_class[number] = static_cast<uint8_t>(value_or_zero(e.niggli_class));
|
||||
message.pixel_sum[number] = static_cast<int64_t>(value_or_zero(e.pixel_sum));
|
||||
}
|
||||
}
|
||||
|
||||
@@ -25,6 +25,8 @@ public:
|
||||
explicit ScanResultGenerator(const DiffractionExperiment& experiment);
|
||||
void Add(const DataMessage& message);
|
||||
ScanResult GetResult() const;
|
||||
|
||||
void FillEndMessage(EndMessage &message) const;
|
||||
};
|
||||
|
||||
|
||||
|
||||
+156
-122
@@ -12,105 +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 | |
|
||||
| - 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.
|
||||
|
||||
@@ -252,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
|
||||
|
||||
@@ -300,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\]
|
||||
- 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
|
||||
@@ -228,6 +228,42 @@ namespace {
|
||||
return {ptr, len};
|
||||
}
|
||||
|
||||
void GetCBORUInt8Array(CborValue &value, std::vector<uint8_t> &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<int32_t> &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<int64_t> &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<float> &v) {
|
||||
if (GetCBORTag(value) != TagFloatLE)
|
||||
throw JFJochException(JFJochExceptionCategory::CBORError, "Incorrect array type tag");
|
||||
@@ -1012,6 +1048,8 @@ namespace {
|
||||
ProcessROIConfig(message, j["roi"]);
|
||||
if (j.contains("gain_file_names"))
|
||||
message.gain_file_names = j["gain_file_names"];
|
||||
if (j.contains("write_images"))
|
||||
message.write_images = j["write_images"];
|
||||
if (j.contains("write_master_file"))
|
||||
message.write_master_file = j["write_master_file"];
|
||||
if (j.contains("data_reduction_factor_serialmx"))
|
||||
@@ -1282,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<float> tmp;
|
||||
GetCBORFloatArray(value, tmp);
|
||||
|
||||
@@ -153,6 +153,23 @@ inline void CBOR_ENC(CborEncoder &encoder, const char* key, const std::vector<fl
|
||||
cborErr(cbor_encode_byte_string(&encoder, (uint8_t *) v.data(), v.size() * sizeof(float)));
|
||||
}
|
||||
|
||||
inline void CBOR_ENC(CborEncoder &encoder, const char* key, const std::vector<uint8_t>& 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<int32_t>& v) {
|
||||
cborErr(cbor_encode_text_stringz(&encoder, key));
|
||||
cborErr(cbor_encode_tag(&encoder, TagSignedInt32BitLE));
|
||||
cborErr(cbor_encode_byte_string(&encoder, reinterpret_cast<const uint8_t *>(v.data()), v.size() * sizeof(int32_t)));
|
||||
}
|
||||
|
||||
inline void CBOR_ENC(CborEncoder &encoder, const char* key, const std::vector<int64_t>& v) {
|
||||
cborErr(cbor_encode_text_stringz(&encoder, key));
|
||||
cborErr(cbor_encode_tag(&encoder, TagSignedInt64BitLE));
|
||||
cborErr(cbor_encode_byte_string(&encoder, reinterpret_cast<const uint8_t *>(v.data()), v.size() * sizeof(int64_t)));
|
||||
}
|
||||
|
||||
inline void CBOR_ENC(CborEncoder &encoder, const char* key, const std::vector<uint64_t>& v) {
|
||||
cborErr(cbor_encode_text_stringz(&encoder, key));
|
||||
@@ -493,6 +510,8 @@ inline void CBOR_ENC_START_USER_DATA(CborEncoder& encoder, const char* key,
|
||||
j["gain_file_names"] = message.gain_file_names;
|
||||
if (message.write_master_file)
|
||||
j["write_master_file"] = message.write_master_file.value();
|
||||
if (message.write_images)
|
||||
j["write_images"] = message.write_images.value();
|
||||
if (message.data_reduction_factor_serialmx)
|
||||
j["data_reduction_factor_serialmx"] = message.data_reduction_factor_serialmx.value();
|
||||
j["experiment_group"] = message.experiment_group;
|
||||
@@ -677,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);
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -9,6 +9,7 @@
|
||||
#include "indexing/AnalyzeIndexing.h"
|
||||
#include "indexing/FFTIndexer.h"
|
||||
#include "lattice_search/LatticeSearch.h"
|
||||
#include "scale_merge/ScaleAll.h"
|
||||
|
||||
IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer)
|
||||
: index_ice_rings(x.GetIndexingSettings().GetIndexIceRings()),
|
||||
@@ -275,7 +276,7 @@ std::optional<RotationIndexerResult> IndexAndRefine::Finalize() {
|
||||
return {};
|
||||
}
|
||||
|
||||
std::optional<ScaleMergeResult> IndexAndRefine::ScaleRotationData(const ScaleMergeOptions &opts) const {
|
||||
std::optional<ScaleMergeResult> IndexAndRefine::ScaleRotationData(const ScaleMergeOptions &opts) {
|
||||
size_t nrefl = 0;
|
||||
for (const auto &i: reflections)
|
||||
nrefl += i.size();
|
||||
|
||||
@@ -12,7 +12,7 @@
|
||||
#include "bragg_prediction/BraggPrediction.h"
|
||||
#include "indexing/IndexerThreadPool.h"
|
||||
#include "lattice_search/LatticeSearch.h"
|
||||
#include "scale_merge/ScaleAndMerge.h"
|
||||
#include "scale_merge/Merge.h"
|
||||
#include "RotationIndexer.h"
|
||||
#include "RotationParameters.h"
|
||||
|
||||
@@ -63,7 +63,7 @@ public:
|
||||
/// Run scale-and-merge on accumulated reflections to refine per-image
|
||||
/// mosaicity (and optionally B-factors / scale factors).
|
||||
/// Returns std::nullopt if there are too few reflections to be meaningful.
|
||||
std::optional<ScaleMergeResult> ScaleRotationData(const ScaleMergeOptions &opts = {}) const;
|
||||
std::optional<ScaleMergeResult> ScaleRotationData(const ScaleMergeOptions &opts = {});
|
||||
|
||||
std::optional<RotationIndexerResult> Finalize();
|
||||
};
|
||||
|
||||
@@ -167,6 +167,7 @@ std::vector<Reflection> IntegrateInternal(const DiffractionExperiment &experimen
|
||||
if (r.observed) {
|
||||
if (experiment.GetPolarizationFactor())
|
||||
r.rlp /= geom.CalcAzIntPolarizationCorr(r.predicted_x, r.predicted_y, experiment.GetPolarizationFactor().value());
|
||||
r.scaling_correction = r.rlp / r.partiality;
|
||||
r.image_number = static_cast<float>(image_number);
|
||||
ret.emplace_back(r);
|
||||
}
|
||||
|
||||
@@ -131,7 +131,8 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
|
||||
.dist_ewald = dist_ewald_sphere,
|
||||
.rlp = 1.0,
|
||||
.partiality = 1.0,
|
||||
.zeta = 1.0
|
||||
.zeta = 1.0,
|
||||
.scaling_correction = 1.0
|
||||
};
|
||||
++i;
|
||||
}
|
||||
|
||||
@@ -115,6 +115,7 @@ namespace {
|
||||
out.rlp = 1.0f;
|
||||
out.partiality = 1.0f;
|
||||
out.zeta = 1.0f;
|
||||
out.scaling_correction = 1.0f;
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
@@ -141,7 +141,8 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys
|
||||
.dist_ewald = dist_ewald_sphere,
|
||||
.rlp = lorentz_reciprocal,
|
||||
.partiality = partiality,
|
||||
.zeta = zeta_abs
|
||||
.zeta = zeta_abs,
|
||||
.scaling_correction = lorentz_reciprocal / partiality,
|
||||
};
|
||||
i++;
|
||||
}
|
||||
|
||||
@@ -158,6 +158,7 @@ namespace {
|
||||
out[count].rlp = lorentz;
|
||||
out[count].partiality = partiality;
|
||||
out[count].zeta = zeta_abs;
|
||||
out[count].scaling_correction = lorentz / partiality;
|
||||
count++;
|
||||
}
|
||||
return count;
|
||||
|
||||
@@ -1,4 +1,10 @@
|
||||
ADD_LIBRARY(JFJochScaleMerge ScaleAndMerge.cpp ScaleAndMerge.h FrenchWilson.cpp FrenchWilson.h
|
||||
ADD_LIBRARY(JFJochScaleMerge ScaleAll.cpp ScaleAll.h FrenchWilson.cpp FrenchWilson.h
|
||||
SearchSpaceGroup.cpp
|
||||
SearchSpaceGroup.h)
|
||||
SearchSpaceGroup.h
|
||||
Merge.cpp
|
||||
Merge.h
|
||||
ScaleOnTheFly.cpp
|
||||
ScaleOnTheFly.h
|
||||
HKLKey.cpp
|
||||
HKLKey.h)
|
||||
TARGET_LINK_LIBRARIES(JFJochScaleMerge Ceres::ceres Eigen3::Eigen JFJochCommon)
|
||||
@@ -112,7 +112,7 @@ FrenchWilson(const std::vector<MergedReflection>& merged,
|
||||
out[i].k = merged[i].k;
|
||||
out[i].l = merged[i].l;
|
||||
out[i].sigmaI = merged[i].sigma;
|
||||
const double I_pos = std::max(merged[i].I, 0.0);
|
||||
const double I_pos = std::max(merged[i].I, 0.0f);
|
||||
out[i].I = I_pos;
|
||||
out[i].F = std::sqrt(I_pos);
|
||||
out[i].sigmaF = 0.0;
|
||||
|
||||
@@ -6,7 +6,7 @@
|
||||
#include <vector>
|
||||
#include <cstdint>
|
||||
|
||||
#include "ScaleAndMerge.h"
|
||||
#include "Merge.h"
|
||||
|
||||
/// Result of the French-Wilson procedure for a single reflection
|
||||
struct FrenchWilsonReflection {
|
||||
|
||||
@@ -0,0 +1,56 @@
|
||||
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include <cmath>
|
||||
|
||||
#include "HKLKey.h"
|
||||
#include "gemmi/symmetry.hpp"
|
||||
|
||||
HKLKey CanonicalHKL(int32_t h, int32_t k, int32_t l, bool merge_friedel, const std::optional<gemmi::SpaceGroup> &in_sg) {
|
||||
HKLKey key{h, k, l, true};
|
||||
|
||||
if (!in_sg.has_value()) {
|
||||
if (!merge_friedel) {
|
||||
const HKLKey neg{-h, -k, -l, true};
|
||||
if (std::tie(key.h, key.k, key.l) < std::tie(neg.h, neg.k, neg.l)) {
|
||||
key.h = -key.h;
|
||||
key.k = -key.k;
|
||||
key.l = -key.l;
|
||||
key.plus = false;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
gemmi::SpaceGroup sg = in_sg.value();
|
||||
const auto ops = sg.operations();
|
||||
const gemmi::ReciprocalAsu asu(&sg);
|
||||
|
||||
const gemmi::Op::Miller in{h, k, l};
|
||||
const auto [hkl, sign_plus] = asu.to_asu_sign(in, ops);
|
||||
|
||||
key.h = hkl[0];
|
||||
key.k = hkl[1];
|
||||
key.l = hkl[2];
|
||||
key.plus = merge_friedel ? true : sign_plus;
|
||||
}
|
||||
return key;
|
||||
}
|
||||
|
||||
HKLKey CanonicalHKL(const Reflection &r, bool merge_friedel, const std::optional<gemmi::SpaceGroup> &sg) {
|
||||
return CanonicalHKL(r.h, r.k, r.l, merge_friedel, sg);
|
||||
}
|
||||
|
||||
HKLKey CanonicalHKL(const MergedReflection &r, bool merge_friedel, const std::optional<gemmi::SpaceGroup> &sg) {
|
||||
return CanonicalHKL(r.h, r.k, r.l, merge_friedel, sg);
|
||||
}
|
||||
|
||||
bool AcceptReflection(const Reflection &r, double d_min_limit) {
|
||||
if (!std::isfinite(r.I))
|
||||
return false;
|
||||
if (!std::isfinite(r.d) || r.d <= 0.0f)
|
||||
return false;
|
||||
if (d_min_limit > 0.0 && r.d < d_min_limit)
|
||||
return false;
|
||||
if (!std::isfinite(r.rlp) || r.rlp == 0.0f)
|
||||
return false;
|
||||
return true;
|
||||
}
|
||||
@@ -0,0 +1,25 @@
|
||||
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <tuple>
|
||||
#include "../../common/Reflection.h"
|
||||
#include "gemmi/symmetry.hpp"
|
||||
|
||||
struct HKLKey {
|
||||
int h = 0;
|
||||
int k = 0;
|
||||
int l = 0;
|
||||
bool plus = true;
|
||||
|
||||
bool operator<(const HKLKey &o) const {
|
||||
return std::tie(h, k, l, plus) < std::tie(o.h, o.k, o.l, o.plus);
|
||||
}
|
||||
};
|
||||
|
||||
HKLKey CanonicalHKL(const Reflection &r, bool merge_friedel, const std::optional<gemmi::SpaceGroup> &sg);
|
||||
HKLKey CanonicalHKL(const MergedReflection &r, bool merge_friedel, const std::optional<gemmi::SpaceGroup> &sg);
|
||||
HKLKey CanonicalHKL(int32_t h, int32_t k, int32_t l, bool merge_friedel, const std::optional<gemmi::SpaceGroup> &sg);
|
||||
|
||||
bool AcceptReflection(const Reflection &r, double d_min_limit);
|
||||
@@ -0,0 +1,300 @@
|
||||
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include "Merge.h"
|
||||
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
#include <map>
|
||||
#include <stdexcept>
|
||||
#include <tuple>
|
||||
#include <unordered_set>
|
||||
|
||||
#include "../../common/ResolutionShells.h"
|
||||
#include "HKLKey.h"
|
||||
|
||||
namespace {
|
||||
struct Obs {
|
||||
const Reflection *r = nullptr;
|
||||
int hkl = -1;
|
||||
double sigma = 1.0;
|
||||
};
|
||||
|
||||
double SafeSigma(double sigma, double min_sigma) {
|
||||
if (!std::isfinite(sigma) || sigma <= 0.0)
|
||||
return min_sigma;
|
||||
return std::max(sigma, min_sigma);
|
||||
}
|
||||
|
||||
double SafeInv(double x, double fallback) {
|
||||
if (!std::isfinite(x) || x == 0.0)
|
||||
return fallback;
|
||||
return 1.0 / x;
|
||||
}
|
||||
|
||||
std::vector<Obs> BuildObservations(const std::vector<std::vector<Reflection>> &observations,
|
||||
const ScaleMergeOptions &opt,
|
||||
std::vector<HKLKey> &slot_to_hkl) {
|
||||
std::map<HKLKey, int> hkl_to_slot;
|
||||
std::vector<Obs> out;
|
||||
|
||||
size_t nrefl = 0;
|
||||
for (const auto &image: observations)
|
||||
nrefl += image.size();
|
||||
out.reserve(nrefl);
|
||||
|
||||
for (const auto &image: observations) {
|
||||
for (const auto &r: image) {
|
||||
|
||||
if (r.scaling_correction <= 0.0 || !std::isfinite(r.scaling_correction))
|
||||
continue;
|
||||
if (!AcceptReflection(r, opt.d_min_limit_A))
|
||||
continue;
|
||||
|
||||
HKLKey key;
|
||||
try {
|
||||
key = CanonicalHKL(r, opt.merge_friedel, opt.space_group);
|
||||
} catch (...) {
|
||||
continue;
|
||||
}
|
||||
|
||||
auto it = hkl_to_slot.find(key);
|
||||
if (it == hkl_to_slot.end()) {
|
||||
const int slot = static_cast<int>(slot_to_hkl.size());
|
||||
it = hkl_to_slot.emplace(key, slot).first;
|
||||
slot_to_hkl.push_back(key);
|
||||
}
|
||||
|
||||
out.push_back({
|
||||
.r = &r,
|
||||
.hkl = it->second,
|
||||
.sigma = SafeSigma(r.sigma, opt.min_sigma)
|
||||
});
|
||||
}
|
||||
}
|
||||
|
||||
return out;
|
||||
}
|
||||
|
||||
ScaleMergeResult InitResult(const std::vector<HKLKey> &slot_to_hkl,
|
||||
const std::vector<Obs> &obs) {
|
||||
ScaleMergeResult out;
|
||||
out.merged.resize(slot_to_hkl.size());
|
||||
|
||||
for (int i = 0; i < static_cast<int>(slot_to_hkl.size()); ++i) {
|
||||
out.merged[i].h = slot_to_hkl[i].h;
|
||||
out.merged[i].k = slot_to_hkl[i].k;
|
||||
out.merged[i].l = slot_to_hkl[i].l;
|
||||
out.merged[i].I = 0.0;
|
||||
out.merged[i].sigma = 0.0;
|
||||
out.merged[i].d = 0.0;
|
||||
}
|
||||
|
||||
std::vector<std::vector<double>> d_values(slot_to_hkl.size());
|
||||
for (const auto &o: obs) {
|
||||
if (std::isfinite(o.r->d) && o.r->d > 0.0f)
|
||||
d_values[o.hkl].push_back(o.r->d);
|
||||
}
|
||||
|
||||
for (int h = 0; h < static_cast<int>(d_values.size()); ++h) {
|
||||
auto &v = d_values[h];
|
||||
if (v.empty())
|
||||
continue;
|
||||
|
||||
std::nth_element(v.begin(), v.begin() + static_cast<long>(v.size() / 2), v.end());
|
||||
out.merged[h].d = v[v.size() / 2];
|
||||
}
|
||||
|
||||
return out;
|
||||
}
|
||||
|
||||
void Merge(size_t nhkl, ScaleMergeResult &out, const std::vector<Obs> &obs) {
|
||||
struct Accum {
|
||||
double sum_wI = 0.0;
|
||||
double sum_w = 0.0;
|
||||
double sum_wsigma2 = 0.0;
|
||||
};
|
||||
|
||||
std::vector<Accum> acc(nhkl);
|
||||
|
||||
for (const auto &o: obs) {
|
||||
const double I_corr = static_cast<double>(o.r->I) * o.r->scaling_correction;
|
||||
const double sigma_corr = o.sigma * o.r->scaling_correction;
|
||||
|
||||
if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0)
|
||||
continue;
|
||||
|
||||
// Extra factor o.r->scaling_correction down-weights weak images / low partiality observations.
|
||||
const double w = o.r->scaling_correction / (sigma_corr * sigma_corr);
|
||||
|
||||
auto &a = acc[o.hkl];
|
||||
a.sum_wI += w * I_corr;
|
||||
a.sum_w += w;
|
||||
a.sum_wsigma2 += w * w * sigma_corr * sigma_corr;
|
||||
}
|
||||
|
||||
for (int h = 0; h < static_cast<int>(nhkl); ++h) {
|
||||
const auto &a = acc[h];
|
||||
if (a.sum_w <= 0.0)
|
||||
continue;
|
||||
|
||||
out.merged[h].I = a.sum_wI / a.sum_w;
|
||||
out.merged[h].sigma = std::sqrt(a.sum_wsigma2) / a.sum_w;
|
||||
}
|
||||
}
|
||||
|
||||
void Stats(const ScaleMergeOptions &opt, ScaleMergeResult &out, const std::vector<Obs> &obs) {
|
||||
constexpr int n_shells = 10;
|
||||
|
||||
float d_min = std::numeric_limits<float>::max();
|
||||
float d_max = 0.0f;
|
||||
|
||||
for (const auto &m: out.merged) {
|
||||
const auto d = static_cast<float>(m.d);
|
||||
if (!std::isfinite(d) || d <= 0.0f)
|
||||
continue;
|
||||
if (opt.d_min_limit_A > 0.0 && d < static_cast<float>(opt.d_min_limit_A))
|
||||
continue;
|
||||
|
||||
d_min = std::min(d_min, d);
|
||||
d_max = std::max(d_max, d);
|
||||
}
|
||||
|
||||
if (!(d_min < d_max && d_min > 0.0f))
|
||||
return;
|
||||
|
||||
const float d_min_pad = d_min * 0.999f;
|
||||
const float d_max_pad = d_max * 1.001f;
|
||||
|
||||
ResolutionShells shells(d_min_pad, d_max_pad, n_shells);
|
||||
const auto shell_mean_1_d2 = shells.GetShellMeanOneOverResSq();
|
||||
const auto shell_min_res = shells.GetShellMinRes();
|
||||
|
||||
std::vector<int> hkl_shell(out.merged.size(), -1);
|
||||
for (int h = 0; h < static_cast<int>(out.merged.size()); ++h) {
|
||||
auto s = shells.GetShell(out.merged[h].d);
|
||||
if (s)
|
||||
hkl_shell[h] = *s;
|
||||
}
|
||||
|
||||
struct PerHKL {
|
||||
double sum_I = 0.0;
|
||||
std::vector<double> I;
|
||||
};
|
||||
|
||||
std::vector<PerHKL> per_hkl(out.merged.size());
|
||||
|
||||
for (const auto &o: obs) {
|
||||
if (o.hkl < 0 || o.hkl >= static_cast<int>(per_hkl.size()))
|
||||
continue;
|
||||
if (hkl_shell[o.hkl] < 0)
|
||||
continue;
|
||||
|
||||
const double I_corr = static_cast<double>(o.r->I) * o.r->scaling_correction;
|
||||
if (!std::isfinite(I_corr))
|
||||
continue;
|
||||
|
||||
per_hkl[o.hkl].sum_I += I_corr;
|
||||
per_hkl[o.hkl].I.push_back(I_corr);
|
||||
}
|
||||
|
||||
struct ShellAccum {
|
||||
int total_obs = 0;
|
||||
std::unordered_set<int> unique;
|
||||
double rmeas_num = 0.0;
|
||||
double rmeas_den = 0.0;
|
||||
double sum_i_over_sigma = 0.0;
|
||||
int n_i_over_sigma = 0;
|
||||
};
|
||||
|
||||
std::vector<ShellAccum> acc(n_shells);
|
||||
|
||||
for (int h = 0; h < static_cast<int>(per_hkl.size()); ++h) {
|
||||
const int s = hkl_shell[h];
|
||||
if (s < 0 || per_hkl[h].I.empty())
|
||||
continue;
|
||||
|
||||
auto &sa = acc[s];
|
||||
const auto &ph = per_hkl[h];
|
||||
const int n = static_cast<int>(ph.I.size());
|
||||
const double mean_I = ph.sum_I / n;
|
||||
|
||||
sa.unique.insert(h);
|
||||
sa.total_obs += n;
|
||||
|
||||
if (n >= 2) {
|
||||
double sum_abs_dev = 0.0;
|
||||
for (double I: ph.I)
|
||||
sum_abs_dev += std::abs(I - mean_I);
|
||||
|
||||
sa.rmeas_num += std::sqrt(static_cast<double>(n) / (n - 1.0)) * sum_abs_dev;
|
||||
}
|
||||
|
||||
for (double I: ph.I)
|
||||
sa.rmeas_den += std::abs(I);
|
||||
|
||||
if (out.merged[h].sigma > 0.0) {
|
||||
sa.sum_i_over_sigma += out.merged[h].I / out.merged[h].sigma;
|
||||
++sa.n_i_over_sigma;
|
||||
}
|
||||
}
|
||||
|
||||
out.statistics.shells.resize(n_shells);
|
||||
|
||||
for (int s = 0; s < n_shells; ++s) {
|
||||
const auto &sa = acc[s];
|
||||
auto &ss = out.statistics.shells[s];
|
||||
|
||||
ss.mean_one_over_d2 = shell_mean_1_d2[s];
|
||||
ss.d_min = shell_min_res[s];
|
||||
ss.d_max = s == 0 ? d_max_pad : shell_min_res[s - 1];
|
||||
ss.total_observations = sa.total_obs;
|
||||
ss.unique_reflections = static_cast<int>(sa.unique.size());
|
||||
ss.rmeas = sa.rmeas_den > 0.0 ? sa.rmeas_num / sa.rmeas_den : 0.0;
|
||||
ss.mean_i_over_sigma = sa.n_i_over_sigma > 0
|
||||
? sa.sum_i_over_sigma / sa.n_i_over_sigma
|
||||
: 0.0;
|
||||
ss.completeness = 0.0;
|
||||
ss.possible_reflections = 0;
|
||||
}
|
||||
|
||||
auto &overall = out.statistics.overall;
|
||||
overall.d_min = d_min;
|
||||
overall.d_max = d_max;
|
||||
|
||||
std::unordered_set<int> all_unique;
|
||||
double rmeas_num = 0.0;
|
||||
double rmeas_den = 0.0;
|
||||
double sum_i_over_sigma = 0.0;
|
||||
int n_i_over_sigma = 0;
|
||||
|
||||
for (const auto &sa: acc) {
|
||||
overall.total_observations += sa.total_obs;
|
||||
all_unique.insert(sa.unique.begin(), sa.unique.end());
|
||||
rmeas_num += sa.rmeas_num;
|
||||
rmeas_den += sa.rmeas_den;
|
||||
sum_i_over_sigma += sa.sum_i_over_sigma;
|
||||
n_i_over_sigma += sa.n_i_over_sigma;
|
||||
}
|
||||
|
||||
overall.unique_reflections = static_cast<int>(all_unique.size());
|
||||
overall.rmeas = rmeas_den > 0.0 ? rmeas_num / rmeas_den : 0.0;
|
||||
overall.mean_i_over_sigma = n_i_over_sigma > 0 ? sum_i_over_sigma / n_i_over_sigma : 0.0;
|
||||
overall.completeness = 0.0;
|
||||
overall.possible_reflections = 0;
|
||||
}
|
||||
}
|
||||
|
||||
ScaleMergeResult MergeReflections(const std::vector<std::vector<Reflection>> &observations,
|
||||
const ScaleMergeOptions &opt) {
|
||||
std::vector<HKLKey> slot_to_hkl;
|
||||
auto obs = BuildObservations(observations, opt, slot_to_hkl);
|
||||
|
||||
auto out = InitResult(slot_to_hkl, obs);
|
||||
|
||||
Merge(slot_to_hkl.size(), out, obs);
|
||||
Stats(opt, out, obs);
|
||||
|
||||
return out;
|
||||
}
|
||||
@@ -0,0 +1,70 @@
|
||||
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <cstdint>
|
||||
#include <optional>
|
||||
#include <vector>
|
||||
|
||||
#include "../../common/Reflection.h"
|
||||
#include "gemmi/symmetry.hpp"
|
||||
|
||||
struct ScaleMergeOptions {
|
||||
int max_num_iterations = 100;
|
||||
double max_solver_time_s = 1.0;
|
||||
|
||||
double image_number_rounding = 1.0;
|
||||
double min_sigma = 1e-3;
|
||||
|
||||
std::optional<gemmi::SpaceGroup> space_group;
|
||||
bool merge_friedel = true;
|
||||
|
||||
std::optional<double> wedge_deg;
|
||||
|
||||
double mosaicity_init_deg = 0.17;
|
||||
double mosaicity_min_deg = 1e-3;
|
||||
double mosaicity_max_deg = 2.0;
|
||||
std::vector<double> mosaicity_init_deg_vec;
|
||||
|
||||
bool regularize_scale_to_one = true;
|
||||
double scale_regularization_sigma = 0.05;
|
||||
|
||||
bool smoothen_g = true;
|
||||
bool smoothen_mos = true;
|
||||
|
||||
double d_min_limit_A = 0.0;
|
||||
|
||||
int64_t image_cluster = 1;
|
||||
|
||||
bool refine_wedge = false;
|
||||
|
||||
enum class PartialityModel { Fixed, Rotation, Unity, Still } partiality_model = PartialityModel::Fixed;
|
||||
};
|
||||
|
||||
struct MergeStatisticsShell {
|
||||
float d_min = 0.0f;
|
||||
float d_max = 0.0f;
|
||||
float mean_one_over_d2 = 0;
|
||||
int total_observations = 0;
|
||||
int unique_reflections = 0;
|
||||
double rmeas = 0.0;
|
||||
double mean_i_over_sigma = 0.0;
|
||||
double completeness = 0.0;
|
||||
int possible_reflections = 0;
|
||||
};
|
||||
|
||||
struct MergeStatistics {
|
||||
std::vector<MergeStatisticsShell> shells;
|
||||
MergeStatisticsShell overall;
|
||||
};
|
||||
|
||||
struct ScaleMergeResult {
|
||||
std::vector<MergedReflection> merged;
|
||||
std::vector<float> image_scale_g;
|
||||
std::vector<float> mosaicity_deg;
|
||||
MergeStatistics statistics;
|
||||
};
|
||||
|
||||
ScaleMergeResult MergeReflections(const std::vector<std::vector<Reflection>> &observations,
|
||||
const ScaleMergeOptions &opt = {});
|
||||
@@ -0,0 +1,461 @@
|
||||
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include "ScaleAll.h"
|
||||
|
||||
#include <ceres/ceres.h>
|
||||
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <iostream>
|
||||
#include <limits>
|
||||
#include <stdexcept>
|
||||
#include <thread>
|
||||
#include "HKLKey.h"
|
||||
|
||||
namespace {
|
||||
struct ScaleObs {
|
||||
const Reflection *r = nullptr;
|
||||
int image = 0;
|
||||
int hkl = -1;
|
||||
double sigma = 1.0;
|
||||
};
|
||||
|
||||
double SafeSigma(double sigma, double min_sigma) {
|
||||
if (!std::isfinite(sigma) || sigma <= 0.0)
|
||||
return min_sigma;
|
||||
return std::max(sigma, min_sigma);
|
||||
}
|
||||
|
||||
double SafeInv(double x, double fallback) {
|
||||
if (!std::isfinite(x) || x == 0.0)
|
||||
return fallback;
|
||||
return 1.0 / x;
|
||||
}
|
||||
|
||||
bool AcceptReflectionForScaling(const Reflection &r, const ScaleMergeOptions &opt) {
|
||||
if (!AcceptReflection(r, opt.d_min_limit_A))
|
||||
return false;
|
||||
|
||||
switch (opt.partiality_model) {
|
||||
case ScaleMergeOptions::PartialityModel::Rotation:
|
||||
return std::isfinite(r.zeta) && r.zeta > 0.0f;
|
||||
|
||||
case ScaleMergeOptions::PartialityModel::Still:
|
||||
return std::isfinite(r.dist_ewald);
|
||||
|
||||
case ScaleMergeOptions::PartialityModel::Fixed:
|
||||
case ScaleMergeOptions::PartialityModel::Unity:
|
||||
return true;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
struct IntensityFixedResidual {
|
||||
IntensityFixedResidual(const Reflection &r, double sigma, double partiality)
|
||||
: Iobs(r.I),
|
||||
weight(SafeInv(sigma, 1.0)),
|
||||
correction(partiality * SafeInv(r.rlp, 1.0)) {
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
bool operator()(const T *const G, const T *const Itrue, T *residual) const {
|
||||
residual[0] = (T(correction) * G[0] * Itrue[0] - T(Iobs)) * T(weight);
|
||||
return true;
|
||||
}
|
||||
|
||||
double Iobs;
|
||||
double weight;
|
||||
double correction;
|
||||
};
|
||||
|
||||
struct IntensityRotationResidual {
|
||||
IntensityRotationResidual(const Reflection &r, double sigma)
|
||||
: Iobs(r.I),
|
||||
weight(SafeInv(sigma, 1.0)),
|
||||
delta_phi_deg(r.delta_phi_deg),
|
||||
lp(SafeInv(r.rlp, 1.0)),
|
||||
c1(r.zeta / std::sqrt(2.0)) {
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
bool operator()(const T *const G,
|
||||
const T *const mosaicity,
|
||||
const T *const Itrue,
|
||||
const T *const wedge,
|
||||
T *residual) const {
|
||||
const T half_wedge = wedge[0] / T(2.0);
|
||||
const T arg_plus = T(delta_phi_deg + half_wedge) * T(c1) / mosaicity[0];
|
||||
const T arg_minus = T(delta_phi_deg - half_wedge) * T(c1) / mosaicity[0];
|
||||
const T partiality = (ceres::erf(arg_plus) - ceres::erf(arg_minus)) / T(2.0);
|
||||
|
||||
residual[0] = (G[0] * partiality * T(lp) * Itrue[0] - T(Iobs)) * T(weight);
|
||||
return true;
|
||||
}
|
||||
|
||||
double Iobs;
|
||||
double weight;
|
||||
double delta_phi_deg;
|
||||
double lp;
|
||||
double c1;
|
||||
};
|
||||
|
||||
struct IntensityStillResidual {
|
||||
IntensityStillResidual(const Reflection &r, double sigma)
|
||||
: Iobs(r.I),
|
||||
weight(SafeInv(sigma, 1.0)),
|
||||
lp(SafeInv(r.rlp, 1.0)),
|
||||
dist_ewald_sq(r.dist_ewald * r.dist_ewald) {
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
bool operator()(const T *const G,
|
||||
const T *const R_sq,
|
||||
const T *const Itrue,
|
||||
T *residual) const {
|
||||
const T partiality = ceres::exp(-T(dist_ewald_sq) / R_sq[0]);
|
||||
residual[0] = (G[0] * partiality * T(lp) * Itrue[0] - T(Iobs)) * T(weight);
|
||||
return true;
|
||||
}
|
||||
|
||||
double Iobs;
|
||||
double weight;
|
||||
double lp;
|
||||
double dist_ewald_sq;
|
||||
};
|
||||
|
||||
struct ScaleRegularizationResidual {
|
||||
explicit ScaleRegularizationResidual(double sigma)
|
||||
: inv_sigma(SafeInv(sigma, 1.0)) {
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
bool operator()(const T *const k, T *residual) const {
|
||||
residual[0] = (k[0] - T(1.0)) * T(inv_sigma);
|
||||
return true;
|
||||
}
|
||||
|
||||
double inv_sigma;
|
||||
};
|
||||
|
||||
struct SmoothnessResidual {
|
||||
explicit SmoothnessResidual(double sigma)
|
||||
: inv_sigma(SafeInv(sigma, 1.0)) {
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
bool operator()(const T *const x0,
|
||||
const T *const x1,
|
||||
const T *const x2,
|
||||
T *residual) const {
|
||||
residual[0] = (ceres::log(x0[0]) + ceres::log(x2[0]) - T(2.0) * ceres::log(x1[0])) * T(inv_sigma);
|
||||
return true;
|
||||
}
|
||||
|
||||
double inv_sigma;
|
||||
};
|
||||
|
||||
std::vector<ScaleObs> BuildScaleObs(const std::vector<std::vector<Reflection>> &observations,
|
||||
const ScaleMergeOptions &opt,
|
||||
std::vector<uint8_t> &image_used,
|
||||
int &nhkl) {
|
||||
std::map<HKLKey, int> hkl_to_slot;
|
||||
std::vector<ScaleObs> obs;
|
||||
|
||||
size_t nrefl = 0;
|
||||
for (const auto &image: observations)
|
||||
nrefl += image.size();
|
||||
obs.reserve(nrefl);
|
||||
|
||||
for (int image = 0; image < static_cast<int>(observations.size()); ++image) {
|
||||
const int image_slot = image / static_cast<int>(opt.image_cluster);
|
||||
|
||||
for (const auto &r: observations[image]) {
|
||||
if (!AcceptReflectionForScaling(r, opt))
|
||||
continue;
|
||||
|
||||
HKLKey key;
|
||||
try {
|
||||
key = CanonicalHKL(r, opt.merge_friedel, opt.space_group);
|
||||
} catch (...) {
|
||||
continue;
|
||||
}
|
||||
auto it = hkl_to_slot.find(key);
|
||||
if (it == hkl_to_slot.end()) {
|
||||
const int slot = static_cast<int>(hkl_to_slot.size());
|
||||
it = hkl_to_slot.emplace(key, slot).first;
|
||||
}
|
||||
|
||||
image_used[image_slot] = 1;
|
||||
|
||||
obs.push_back({
|
||||
.r = &r,
|
||||
.image = image_slot,
|
||||
.hkl = it->second,
|
||||
.sigma = SafeSigma(r.sigma, opt.min_sigma)
|
||||
});
|
||||
}
|
||||
}
|
||||
|
||||
nhkl = static_cast<int>(hkl_to_slot.size());
|
||||
return obs;
|
||||
}
|
||||
|
||||
std::vector<double> InitialIntensities(int nhkl,
|
||||
const ScaleMergeOptions &opt,
|
||||
const std::vector<ScaleObs> &obs) {
|
||||
std::vector<std::vector<double>> values(nhkl);
|
||||
|
||||
for (const auto &o: obs)
|
||||
values[o.hkl].push_back(o.r->I);
|
||||
|
||||
std::vector<double> Itrue(nhkl, opt.min_sigma);
|
||||
|
||||
for (int h = 0; h < nhkl; ++h) {
|
||||
auto &v = values[h];
|
||||
if (v.empty())
|
||||
continue;
|
||||
|
||||
std::nth_element(v.begin(), v.begin() + static_cast<long>(v.size() / 2), v.end());
|
||||
|
||||
Itrue[h] = v[v.size() / 2];
|
||||
if (!std::isfinite(Itrue[h]) || Itrue[h] <= opt.min_sigma)
|
||||
Itrue[h] = opt.min_sigma;
|
||||
}
|
||||
|
||||
return Itrue;
|
||||
}
|
||||
|
||||
void Scale(const ScaleMergeOptions &opt,
|
||||
const std::vector<ScaleObs> &obs,
|
||||
const std::vector<uint8_t> &image_used,
|
||||
int nhkl,
|
||||
std::vector<double> &G,
|
||||
std::vector<double> &mosaicity,
|
||||
std::vector<double> &R_sq) {
|
||||
ceres::Problem problem;
|
||||
|
||||
auto Itrue = InitialIntensities(nhkl, opt, obs);
|
||||
double wedge = opt.wedge_deg.value_or(0.0);
|
||||
|
||||
for (const auto &o: obs) {
|
||||
switch (opt.partiality_model) {
|
||||
case ScaleMergeOptions::PartialityModel::Rotation: {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<IntensityRotationResidual, 1, 1, 1, 1, 1>(
|
||||
new IntensityRotationResidual(*o.r, o.sigma));
|
||||
problem.AddResidualBlock(cost, nullptr, &G[o.image], &mosaicity[o.image], &Itrue[o.hkl], &wedge);
|
||||
break;
|
||||
}
|
||||
|
||||
case ScaleMergeOptions::PartialityModel::Still: {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<IntensityStillResidual, 1, 1, 1, 1>(
|
||||
new IntensityStillResidual(*o.r, o.sigma));
|
||||
problem.AddResidualBlock(cost, nullptr, &G[o.image], &R_sq[o.image], &Itrue[o.hkl]);
|
||||
break;
|
||||
}
|
||||
|
||||
case ScaleMergeOptions::PartialityModel::Unity: {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<IntensityFixedResidual, 1, 1, 1>(
|
||||
new IntensityFixedResidual(*o.r, o.sigma, 1.0));
|
||||
problem.AddResidualBlock(cost, nullptr, &G[o.image], &Itrue[o.hkl]);
|
||||
break;
|
||||
}
|
||||
|
||||
case ScaleMergeOptions::PartialityModel::Fixed: {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<IntensityFixedResidual, 1, 1, 1>(
|
||||
new IntensityFixedResidual(*o.r, o.sigma, o.r->partiality));
|
||||
problem.AddResidualBlock(cost, nullptr, &G[o.image], &Itrue[o.hkl]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = 0; i < static_cast<int>(G.size()); ++i) {
|
||||
if (!image_used[i])
|
||||
continue;
|
||||
|
||||
problem.SetParameterLowerBound(&G[i], 0, 1e-12);
|
||||
|
||||
if (opt.regularize_scale_to_one) {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<ScaleRegularizationResidual, 1, 1>(
|
||||
new ScaleRegularizationResidual(opt.scale_regularization_sigma));
|
||||
problem.AddResidualBlock(cost, nullptr, &G[i]);
|
||||
}
|
||||
}
|
||||
|
||||
if (opt.smoothen_g) {
|
||||
for (int i = 0; i + 2 < static_cast<int>(G.size()); ++i) {
|
||||
if (!(image_used[i] && image_used[i + 1] && image_used[i + 2]))
|
||||
continue;
|
||||
|
||||
auto *cost = new ceres::AutoDiffCostFunction<SmoothnessResidual, 1, 1, 1, 1>(
|
||||
new SmoothnessResidual(0.05));
|
||||
problem.AddResidualBlock(cost, nullptr, &G[i], &G[i + 1], &G[i + 2]);
|
||||
}
|
||||
}
|
||||
|
||||
if (opt.partiality_model == ScaleMergeOptions::PartialityModel::Rotation) {
|
||||
for (int i = 0; i < static_cast<int>(mosaicity.size()); ++i) {
|
||||
if (!image_used[i])
|
||||
continue;
|
||||
|
||||
problem.SetParameterLowerBound(&mosaicity[i], 0, opt.mosaicity_min_deg);
|
||||
problem.SetParameterUpperBound(&mosaicity[i], 0, opt.mosaicity_max_deg);
|
||||
}
|
||||
|
||||
if (opt.smoothen_mos) {
|
||||
for (int i = 0; i + 2 < static_cast<int>(mosaicity.size()); ++i) {
|
||||
if (!(image_used[i] && image_used[i + 1] && image_used[i + 2]))
|
||||
continue;
|
||||
|
||||
auto *cost = new ceres::AutoDiffCostFunction<SmoothnessResidual, 1, 1, 1, 1>(
|
||||
new SmoothnessResidual(0.05));
|
||||
problem.AddResidualBlock(cost, nullptr, &mosaicity[i], &mosaicity[i + 1], &mosaicity[i + 2]);
|
||||
}
|
||||
}
|
||||
|
||||
if (!opt.refine_wedge)
|
||||
problem.SetParameterBlockConstant(&wedge);
|
||||
else
|
||||
problem.SetParameterLowerBound(&wedge, 0, 0.0);
|
||||
}
|
||||
|
||||
if (opt.partiality_model == ScaleMergeOptions::PartialityModel::Still) {
|
||||
for (int i = 0; i < static_cast<int>(R_sq.size()); ++i) {
|
||||
if (!image_used[i])
|
||||
continue;
|
||||
|
||||
problem.SetParameterLowerBound(&R_sq[i], 0, 1e-9);
|
||||
problem.SetParameterUpperBound(&R_sq[i], 0, 1.0);
|
||||
}
|
||||
}
|
||||
|
||||
unsigned int hw = std::thread::hardware_concurrency();
|
||||
if (hw == 0)
|
||||
hw = 1;
|
||||
|
||||
ceres::Solver::Options options;
|
||||
options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
|
||||
options.minimizer_progress_to_stdout = true;
|
||||
options.max_num_iterations = opt.max_num_iterations;
|
||||
options.max_solver_time_in_seconds = opt.max_solver_time_s;
|
||||
options.num_threads = static_cast<int>(hw);
|
||||
options.function_tolerance = 1e-4;
|
||||
|
||||
ceres::Solver::Summary summary;
|
||||
ceres::Solve(options, &problem, &summary);
|
||||
|
||||
std::cout << summary.FullReport() << std::endl;
|
||||
}
|
||||
|
||||
double Partiality(const Reflection &r,
|
||||
const ScaleMergeOptions &opt,
|
||||
int image_slot,
|
||||
const std::vector<double> &mosaicity,
|
||||
const std::vector<double> &R_sq) {
|
||||
switch (opt.partiality_model) {
|
||||
case ScaleMergeOptions::PartialityModel::Fixed:
|
||||
return r.partiality;
|
||||
|
||||
case ScaleMergeOptions::PartialityModel::Unity:
|
||||
return 1.0;
|
||||
|
||||
case ScaleMergeOptions::PartialityModel::Rotation: {
|
||||
const double half_wedge = opt.wedge_deg.value_or(0.0) / 2.0;
|
||||
const double c1 = r.zeta / std::sqrt(2.0);
|
||||
const double arg_plus = (r.delta_phi_deg + half_wedge) * c1 / mosaicity[image_slot];
|
||||
const double arg_minus = (r.delta_phi_deg - half_wedge) * c1 / mosaicity[image_slot];
|
||||
|
||||
return (std::erf(arg_plus) - std::erf(arg_minus)) / 2.0;
|
||||
}
|
||||
|
||||
case ScaleMergeOptions::PartialityModel::Still:
|
||||
return std::exp(-r.dist_ewald * r.dist_ewald / R_sq[image_slot]);
|
||||
}
|
||||
|
||||
return 1.0;
|
||||
}
|
||||
|
||||
void CalcCorrections(std::vector<std::vector<Reflection>> &observations,
|
||||
const ScaleMergeOptions &opt,
|
||||
const std::vector<double> &G,
|
||||
const std::vector<double> &mosaicity,
|
||||
const std::vector<double> &R_sq) {
|
||||
|
||||
size_t nrefl = 0;
|
||||
for (const auto &image: observations)
|
||||
nrefl += image.size();
|
||||
|
||||
for (int image = 0; image < static_cast<int>(observations.size()); ++image) {
|
||||
const int image_slot = image / static_cast<int>(opt.image_cluster);
|
||||
|
||||
for (auto &r: observations[image]) {
|
||||
if (!AcceptReflectionForScaling(r, opt)) {
|
||||
r.scaling_correction = 0.0;
|
||||
continue;
|
||||
}
|
||||
|
||||
const double partiality = Partiality(r, opt, image_slot, mosaicity, R_sq);
|
||||
if (!std::isfinite(partiality) || partiality < 0.01) {
|
||||
r.partiality = 0.0;
|
||||
r.scaling_correction = 0.0;
|
||||
} else {
|
||||
const double denom = G[image_slot] * partiality;
|
||||
const double correction = denom > 0.0 ? r.rlp / denom : 0.0;
|
||||
r.partiality = partiality;
|
||||
r.scaling_correction = std::isfinite(correction) ? correction : 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
ScaleMergeResult ScaleAndMergeReflectionsCeres(std::vector<std::vector<Reflection>> &observations,
|
||||
const ScaleMergeOptions &opt) {
|
||||
if (opt.image_cluster <= 0)
|
||||
throw std::invalid_argument("image_cluster must be positive");
|
||||
|
||||
const size_t n_image_slots = observations.size() / opt.image_cluster +
|
||||
(observations.size() % opt.image_cluster > 0 ? 1 : 0);
|
||||
|
||||
std::vector<uint8_t> image_used(n_image_slots, 0);
|
||||
|
||||
int nhkl = 0;
|
||||
auto scale_obs = BuildScaleObs(observations, opt, image_used, nhkl);
|
||||
|
||||
std::vector<double> G(n_image_slots, 1.0);
|
||||
std::vector<double> mosaicity(n_image_slots, opt.mosaicity_init_deg);
|
||||
std::vector<double> R_sq(n_image_slots, 0.001 * 0.001);
|
||||
|
||||
for (int i = 0; i < static_cast<int>(n_image_slots); ++i) {
|
||||
if (!image_used[i]) {
|
||||
G[i] = NAN;
|
||||
mosaicity[i] = NAN;
|
||||
R_sq[i] = NAN;
|
||||
} else if (opt.mosaicity_init_deg_vec.size() > static_cast<size_t>(i) &&
|
||||
std::isfinite(opt.mosaicity_init_deg_vec[i])) {
|
||||
mosaicity[i] = opt.mosaicity_init_deg_vec[i];
|
||||
}
|
||||
}
|
||||
|
||||
Scale(opt, scale_obs, image_used, nhkl, G, mosaicity, R_sq);
|
||||
|
||||
CalcCorrections(observations, opt, G, mosaicity, R_sq);
|
||||
|
||||
auto out = MergeReflections(observations, opt);
|
||||
|
||||
out.image_scale_g.resize(observations.size(), NAN);
|
||||
out.mosaicity_deg.resize(observations.size(), NAN);
|
||||
|
||||
for (int image = 0; image < static_cast<int>(observations.size()); ++image) {
|
||||
const int image_slot = image / static_cast<int>(opt.image_cluster);
|
||||
|
||||
if (image_slot < static_cast<int>(image_used.size()) && image_used[image_slot]) {
|
||||
out.image_scale_g[image] = G[image_slot];
|
||||
out.mosaicity_deg[image] = mosaicity[image_slot];
|
||||
}
|
||||
}
|
||||
|
||||
return out;
|
||||
}
|
||||
@@ -0,0 +1,9 @@
|
||||
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#pragma once
|
||||
|
||||
#include "Merge.h"
|
||||
|
||||
ScaleMergeResult ScaleAndMergeReflectionsCeres(std::vector<std::vector<Reflection>>& observations,
|
||||
const ScaleMergeOptions& opt = {});
|
||||
@@ -1,789 +0,0 @@
|
||||
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include "ScaleAndMerge.h"
|
||||
|
||||
#include <ceres/ceres.h>
|
||||
|
||||
#include <thread>
|
||||
#include <iostream>
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
#include <stdexcept>
|
||||
#include <tuple>
|
||||
#include <unordered_map>
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
#include "../../common/ResolutionShells.h"
|
||||
|
||||
namespace {
|
||||
struct HKLKey {
|
||||
int64_t h = 0;
|
||||
int64_t k = 0;
|
||||
int64_t l = 0;
|
||||
bool is_positive = true; // only relevant if opt.merge_friedel == false
|
||||
|
||||
bool operator==(const HKLKey &o) const noexcept {
|
||||
return h == o.h && k == o.k && l == o.l && is_positive == o.is_positive;
|
||||
}
|
||||
};
|
||||
|
||||
struct HKLKeyHash {
|
||||
size_t operator()(const HKLKey &key) const noexcept {
|
||||
auto mix = [](uint64_t x) {
|
||||
x ^= x >> 33;
|
||||
x *= 0xff51afd7ed558ccdULL;
|
||||
x ^= x >> 33;
|
||||
x *= 0xc4ceb9fe1a85ec53ULL;
|
||||
x ^= x >> 33;
|
||||
return x;
|
||||
};
|
||||
const uint64_t a = static_cast<uint64_t>(key.h);
|
||||
const uint64_t b = static_cast<uint64_t>(key.k);
|
||||
const uint64_t c = static_cast<uint64_t>(key.l);
|
||||
const uint64_t d = static_cast<uint64_t>(key.is_positive ? 1 : 0);
|
||||
return static_cast<size_t>(mix(a) ^ (mix(b) << 1) ^ (mix(c) << 2) ^ (mix(d) << 3));
|
||||
}
|
||||
};
|
||||
|
||||
inline int RoundImageId(float image_number, double rounding_step) {
|
||||
if (!(rounding_step > 0.0))
|
||||
rounding_step = 1.0;
|
||||
const double x = static_cast<double>(image_number) / rounding_step;
|
||||
const double r = std::round(x) * rounding_step;
|
||||
return static_cast<int>(std::llround(r / rounding_step));
|
||||
}
|
||||
|
||||
inline double SafeSigma(double s, double min_sigma) {
|
||||
if (!std::isfinite(s) || s <= 0.0)
|
||||
return min_sigma;
|
||||
return std::max(s, min_sigma);
|
||||
}
|
||||
|
||||
inline double SafeD(double d) {
|
||||
if (!std::isfinite(d) || d <= 0.0)
|
||||
return std::numeric_limits<double>::quiet_NaN();
|
||||
return d;
|
||||
}
|
||||
|
||||
inline int SafeToInt(int64_t x) {
|
||||
if (x < std::numeric_limits<int>::min() || x > std::numeric_limits<int>::max())
|
||||
throw std::out_of_range("HKL index out of int range for Gemmi");
|
||||
return static_cast<int>(x);
|
||||
}
|
||||
|
||||
inline double SafeInv(double x, double fallback) {
|
||||
if (!std::isfinite(x) || x == 0.0)
|
||||
return fallback;
|
||||
return 1.0 / x;
|
||||
}
|
||||
|
||||
inline HKLKey CanonicalizeHKLKey(const Reflection &r, const ScaleMergeOptions &opt) {
|
||||
HKLKey key{};
|
||||
key.h = r.h;
|
||||
key.k = r.k;
|
||||
key.l = r.l;
|
||||
key.is_positive = true;
|
||||
|
||||
if (!opt.space_group.has_value()) {
|
||||
if (!opt.merge_friedel) {
|
||||
const HKLKey neg{-r.h, -r.k, -r.l, true};
|
||||
const bool pos = std::tie(key.h, key.k, key.l) >= std::tie(neg.h, neg.k, neg.l);
|
||||
if (!pos) {
|
||||
key.h = -key.h;
|
||||
key.k = -key.k;
|
||||
key.l = -key.l;
|
||||
key.is_positive = false;
|
||||
}
|
||||
}
|
||||
return key;
|
||||
}
|
||||
|
||||
const gemmi::SpaceGroup &sg = *opt.space_group;
|
||||
const gemmi::GroupOps gops = sg.operations();
|
||||
const gemmi::ReciprocalAsu rasu(&sg);
|
||||
|
||||
const gemmi::Op::Miller in{{SafeToInt(r.h), SafeToInt(r.k), SafeToInt(r.l)}};
|
||||
const auto [asu_hkl, sign_plus] = rasu.to_asu_sign(in, gops);
|
||||
|
||||
key.h = asu_hkl[0];
|
||||
key.k = asu_hkl[1];
|
||||
key.l = asu_hkl[2];
|
||||
key.is_positive = opt.merge_friedel ? true : sign_plus;
|
||||
return key;
|
||||
}
|
||||
|
||||
struct IntensityRotResidual {
|
||||
IntensityRotResidual(const Reflection &r, double sigma_obs, double wedge_deg)
|
||||
: Iobs_(static_cast<double>(r.I)),
|
||||
weight_(SafeInv(sigma_obs, 1.0)),
|
||||
delta_phi_(r.delta_phi_deg),
|
||||
lp_(SafeInv(r.rlp, 1.0)),
|
||||
c1_(r.zeta / std::sqrt(2.0)),
|
||||
partiality_(r.partiality) {
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
bool operator()(const T *const G,
|
||||
const T *const mosaicity,
|
||||
const T *const Itrue,
|
||||
const T *const wedge,
|
||||
T *residual) const {
|
||||
T partiality;
|
||||
if (mosaicity[0] >= 0.0) {
|
||||
const T half_wedge = wedge[0] / T(2.0);
|
||||
const T arg_plus = T(delta_phi_ + half_wedge) * T(c1_) / mosaicity[0];
|
||||
const T arg_minus = T(delta_phi_ - half_wedge) * T(c1_) / mosaicity[0];
|
||||
partiality = (ceres::erf(arg_plus) - ceres::erf(arg_minus)) / T(2.0);
|
||||
} else
|
||||
partiality = T(1.0);
|
||||
|
||||
const T Ipred = G[0] * partiality * T(lp_) * Itrue[0];
|
||||
residual[0] = (Ipred - T(Iobs_)) * T(weight_);
|
||||
return true;
|
||||
}
|
||||
|
||||
double Iobs_;
|
||||
double weight_;
|
||||
double delta_phi_;
|
||||
double lp_;
|
||||
double c1_;
|
||||
double partiality_;
|
||||
};
|
||||
|
||||
struct IntensityStillResidual {
|
||||
IntensityStillResidual(const Reflection &r, double sigma_obs)
|
||||
: Iobs_(static_cast<double>(r.I)),
|
||||
weight_(SafeInv(sigma_obs, 1.0)),
|
||||
lp_(SafeInv(r.rlp, 1.0)),
|
||||
dist_ewald_sq_(r.dist_ewald * r.dist_ewald) {
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
bool operator()(const T *const G,
|
||||
const T *const R,
|
||||
const T *const Itrue,
|
||||
T *residual) const {
|
||||
const T partiality = ceres::exp(-T(dist_ewald_sq_)/R[0]);
|
||||
const T Ipred = G[0] * partiality * T(lp_) * Itrue[0];
|
||||
residual[0] = (Ipred - T(Iobs_)) * T(weight_);
|
||||
return true;
|
||||
}
|
||||
|
||||
double Iobs_;
|
||||
double weight_;
|
||||
double lp_;
|
||||
double dist_ewald_sq_;
|
||||
};
|
||||
|
||||
struct IntensityFixedResidual {
|
||||
IntensityFixedResidual(const Reflection &r, double sigma_obs, double partiality)
|
||||
: Iobs_(static_cast<double>(r.I)),
|
||||
weight_(SafeInv(sigma_obs, 1.0)),
|
||||
corr_(partiality * SafeInv(r.rlp, 1.0))
|
||||
{}
|
||||
|
||||
template<typename T>
|
||||
bool operator()(const T *const G,
|
||||
const T *const Itrue,
|
||||
T *residual) const {
|
||||
const T Ipred = T(corr_) * G[0] * Itrue[0];
|
||||
residual[0] = (Ipred - T(Iobs_)) * T(weight_);
|
||||
return true;
|
||||
}
|
||||
|
||||
double Iobs_;
|
||||
double weight_;
|
||||
double corr_;
|
||||
};
|
||||
|
||||
struct ScaleRegularizationResidual {
|
||||
explicit ScaleRegularizationResidual(double sigma_k)
|
||||
: inv_sigma_(SafeInv(sigma_k, 1.0)) {
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
bool operator()(const T *const k, T *residual) const {
|
||||
residual[0] = (k[0] - T(1.0)) * T(inv_sigma_);
|
||||
return true;
|
||||
}
|
||||
|
||||
double inv_sigma_;
|
||||
};
|
||||
|
||||
struct SmoothnessRegularizationResidual {
|
||||
explicit SmoothnessRegularizationResidual(double sigma)
|
||||
: inv_sigma_(SafeInv(sigma, 1.0)) {
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
bool operator()(const T *const k0,
|
||||
const T *const k1,
|
||||
const T *const k2,
|
||||
T *residual) const {
|
||||
residual[0] = (ceres::log(k0[0]) + ceres::log(k2[0]) - T(2.0) * ceres::log(k1[0])) * T(inv_sigma_);
|
||||
return true;
|
||||
}
|
||||
|
||||
double inv_sigma_;
|
||||
};
|
||||
|
||||
struct ObsRef {
|
||||
const Reflection *r = nullptr;
|
||||
int img_id = 0;
|
||||
int hkl_slot = -1;
|
||||
double sigma = 0.0;
|
||||
};
|
||||
|
||||
struct CorrectedObs {
|
||||
int hkl_slot;
|
||||
double I_corr;
|
||||
double sigma_corr;
|
||||
double weight;
|
||||
};
|
||||
|
||||
void scale(const ScaleMergeOptions &opt,
|
||||
std::vector<double> &g,
|
||||
std::vector<double> &mosaicity,
|
||||
std::vector<double> &R_sq,
|
||||
const std::vector<uint8_t> &image_slot_used,
|
||||
bool rotation_crystallography,
|
||||
size_t nhkl,
|
||||
const std::vector<ObsRef> &obs) {
|
||||
ceres::Problem problem;
|
||||
|
||||
std::vector<double> Itrue(nhkl, 0.0);
|
||||
|
||||
// Initialize Itrue from per-HKL median of observed intensities
|
||||
{
|
||||
std::vector<std::vector<double> > per_hkl_I(nhkl);
|
||||
for (const auto &o: obs) {
|
||||
per_hkl_I[o.hkl_slot].push_back(static_cast<double>(o.r->I));
|
||||
}
|
||||
for (int h = 0; h < nhkl; ++h) {
|
||||
auto &v = per_hkl_I[h];
|
||||
if (v.empty()) {
|
||||
Itrue[h] = std::max(opt.min_sigma, 1e-6);
|
||||
continue;
|
||||
}
|
||||
std::nth_element(v.begin(), v.begin() + static_cast<long>(v.size() / 2), v.end());
|
||||
double med = v[v.size() / 2];
|
||||
if (!std::isfinite(med) || med <= opt.min_sigma)
|
||||
med = opt.min_sigma;
|
||||
Itrue[h] = med;
|
||||
}
|
||||
}
|
||||
|
||||
double wedge = opt.wedge_deg.value_or(0.0);
|
||||
|
||||
std::vector<bool> is_valid_hkl_slot(nhkl, false);
|
||||
for (const auto &o: obs) {
|
||||
switch (opt.partiality_model) {
|
||||
case ScaleMergeOptions::PartialityModel::Rotation: {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<IntensityRotResidual, 1, 1, 1, 1, 1>(
|
||||
new IntensityRotResidual(*o.r, o.sigma, opt.wedge_deg.value_or(0.0)));
|
||||
problem.AddResidualBlock(cost,
|
||||
nullptr,
|
||||
&g[o.img_id],
|
||||
&mosaicity[o.img_id],
|
||||
&Itrue[o.hkl_slot],
|
||||
&wedge);
|
||||
}
|
||||
break;
|
||||
case ScaleMergeOptions::PartialityModel::Still: {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<IntensityStillResidual, 1, 1, 1, 1>(
|
||||
new IntensityStillResidual(*o.r, o.sigma));
|
||||
problem.AddResidualBlock(cost,
|
||||
nullptr,
|
||||
&g[o.img_id],
|
||||
&R_sq[o.img_id],
|
||||
&Itrue[o.hkl_slot]);
|
||||
}
|
||||
break;
|
||||
case ScaleMergeOptions::PartialityModel::Unity: {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<IntensityFixedResidual, 1, 1, 1>(
|
||||
new IntensityFixedResidual(*o.r, o.sigma, 1.0));
|
||||
problem.AddResidualBlock(cost,
|
||||
nullptr,
|
||||
&g[o.img_id],
|
||||
&Itrue[o.hkl_slot]);
|
||||
}
|
||||
break;
|
||||
case ScaleMergeOptions::PartialityModel::Fixed: {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<IntensityFixedResidual, 1, 1, 1>(
|
||||
new IntensityFixedResidual(*o.r, o.sigma, o.r->partiality));
|
||||
problem.AddResidualBlock(cost,
|
||||
nullptr,
|
||||
&g[o.img_id],
|
||||
&Itrue[o.hkl_slot]);
|
||||
}
|
||||
break;
|
||||
}
|
||||
is_valid_hkl_slot[o.hkl_slot] = true;
|
||||
}
|
||||
|
||||
for (int i = 0; i < g.size(); ++i) {
|
||||
if (image_slot_used[i]) {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<ScaleRegularizationResidual, 1, 1>(
|
||||
new ScaleRegularizationResidual(0.05));
|
||||
problem.AddResidualBlock(cost, nullptr, &g[i]);
|
||||
}
|
||||
}
|
||||
|
||||
if (rotation_crystallography) {
|
||||
if (opt.smoothen_g) {
|
||||
for (int i = 0; i < g.size() - 2; ++i) {
|
||||
if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<SmoothnessRegularizationResidual, 1, 1, 1, 1>(
|
||||
new SmoothnessRegularizationResidual(0.05));
|
||||
|
||||
problem.AddResidualBlock(cost, nullptr, &g[i], &g[i + 1], &g[i + 2]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (opt.smoothen_mos && opt.partiality_model == ScaleMergeOptions::PartialityModel::Rotation) {
|
||||
for (int i = 0; i < mosaicity.size() - 2; ++i) {
|
||||
if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<SmoothnessRegularizationResidual, 1, 1, 1, 1>(
|
||||
new SmoothnessRegularizationResidual(0.05));
|
||||
problem.AddResidualBlock(cost, nullptr, &mosaicity[i], &mosaicity[i + 1], &mosaicity[i + 2]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (opt.partiality_model == ScaleMergeOptions::PartialityModel::Still) {
|
||||
for (int i = 0; i < R_sq.size(); ++i) {
|
||||
if (image_slot_used[i]) {
|
||||
problem.SetParameterLowerBound(&R_sq[i], 0, 1e-9);
|
||||
problem.SetParameterUpperBound(&R_sq[i], 0, 1.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Scaling factors must be always positive
|
||||
for (int i = 0; i < g.size(); i++) {
|
||||
if (image_slot_used[i])
|
||||
problem.SetParameterLowerBound(&g[i], 0, 1e-12);
|
||||
}
|
||||
|
||||
// Mosaicity refinement + bounds
|
||||
if (opt.partiality_model == ScaleMergeOptions::PartialityModel::Rotation) {
|
||||
for (int i = 0; i < mosaicity.size(); ++i) {
|
||||
if (image_slot_used[i]) {
|
||||
problem.SetParameterLowerBound(&mosaicity[i], 0, opt.mosaicity_min_deg);
|
||||
problem.SetParameterUpperBound(&mosaicity[i], 0, opt.mosaicity_max_deg);
|
||||
}
|
||||
}
|
||||
if (!opt.refine_wedge)
|
||||
problem.SetParameterBlockConstant(&wedge);
|
||||
else
|
||||
problem.SetParameterLowerBound(&wedge, 0, 0.0);
|
||||
}
|
||||
|
||||
// use all available threads
|
||||
unsigned int hw = std::thread::hardware_concurrency();
|
||||
if (hw == 0)
|
||||
hw = 1; // fallback
|
||||
|
||||
ceres::Solver::Options options;
|
||||
|
||||
options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
|
||||
options.minimizer_progress_to_stdout = true;
|
||||
options.max_num_iterations = opt.max_num_iterations;
|
||||
options.max_solver_time_in_seconds = opt.max_solver_time_s;
|
||||
options.num_threads = static_cast<int>(hw);
|
||||
options.function_tolerance = 1e-4;
|
||||
|
||||
ceres::Solver::Summary summary;
|
||||
ceres::Solve(options, &problem, &summary);
|
||||
std::cout << summary.FullReport() << std::endl;
|
||||
}
|
||||
|
||||
void merge(size_t nhkl, ScaleMergeResult &out, const std::vector<CorrectedObs> &corr_obs) {
|
||||
// ---- Merge (XDS/XSCALE style: inverse-variance weighted mean) ----
|
||||
|
||||
struct HKLAccum {
|
||||
double sum_wI = 0.0;
|
||||
double sum_w = 0.0;
|
||||
double sum_wsigma2 = 0.0;
|
||||
};
|
||||
std::vector<HKLAccum> accum(nhkl);
|
||||
|
||||
for (const auto &co: corr_obs) {
|
||||
const double w = co.weight / (co.sigma_corr * co.sigma_corr);
|
||||
auto &a = accum[co.hkl_slot];
|
||||
a.sum_wI += w * co.I_corr;
|
||||
a.sum_w += w;
|
||||
a.sum_wsigma2 += w * w * co.sigma_corr * co.sigma_corr;
|
||||
}
|
||||
|
||||
for (int h = 0; h < nhkl; ++h) {
|
||||
const auto &a = accum[h];
|
||||
if (a.sum_w <= 0.0)
|
||||
continue;
|
||||
|
||||
out.merged[h].I = a.sum_wI / a.sum_w;
|
||||
out.merged[h].sigma = std::sqrt(a.sum_wsigma2) / a.sum_w;
|
||||
}
|
||||
}
|
||||
|
||||
void stats(const ScaleMergeOptions &opt, size_t nhkl, ScaleMergeResult &out, const std::vector<CorrectedObs> &corr_obs)
|
||||
// ---- Compute per-shell merging statistics ----
|
||||
{
|
||||
constexpr int kStatShells = 10;
|
||||
|
||||
float stat_d_min = std::numeric_limits<float>::max();
|
||||
float stat_d_max = 0.0f;
|
||||
for (int h = 0; h < nhkl; ++h) {
|
||||
const auto d = static_cast<float>(out.merged[h].d);
|
||||
if (std::isfinite(d) && d > 0.0f) {
|
||||
if (opt.d_min_limit_A > 0.0 && d < static_cast<float>(opt.d_min_limit_A))
|
||||
continue;
|
||||
stat_d_min = std::min(stat_d_min, d);
|
||||
stat_d_max = std::max(stat_d_max, d);
|
||||
}
|
||||
}
|
||||
|
||||
if (stat_d_min < stat_d_max && stat_d_min > 0.0f) {
|
||||
const float d_min_pad = stat_d_min * 0.999f;
|
||||
const float d_max_pad = stat_d_max * 1.001f;
|
||||
ResolutionShells stat_shells(d_min_pad, d_max_pad, kStatShells);
|
||||
|
||||
const auto shell_mean_1_d2 = stat_shells.GetShellMeanOneOverResSq();
|
||||
const auto shell_min_res = stat_shells.GetShellMinRes();
|
||||
|
||||
// Assign each unique reflection to a shell
|
||||
std::vector<int32_t> hkl_shell(nhkl, -1);
|
||||
for (int h = 0; h < nhkl; ++h) {
|
||||
const auto d = static_cast<float>(out.merged[h].d);
|
||||
if (std::isfinite(d) && d > 0.0f) {
|
||||
if (opt.d_min_limit_A > 0.0 && d < static_cast<float>(opt.d_min_limit_A))
|
||||
continue;
|
||||
auto s = stat_shells.GetShell(d);
|
||||
if (s.has_value())
|
||||
hkl_shell[h] = s.value();
|
||||
}
|
||||
}
|
||||
|
||||
// Per-HKL: collect corrected intensities for Rmeas
|
||||
struct HKLStats {
|
||||
double sum_I = 0.0;
|
||||
int n = 0;
|
||||
std::vector<double> I_list;
|
||||
};
|
||||
std::vector<HKLStats> per_hkl(nhkl);
|
||||
|
||||
for (const auto &co: corr_obs) {
|
||||
if (hkl_shell[co.hkl_slot] < 0)
|
||||
continue;
|
||||
auto &hs = per_hkl[co.hkl_slot];
|
||||
hs.sum_I += co.I_corr;
|
||||
hs.n += 1;
|
||||
hs.I_list.push_back(co.I_corr);
|
||||
}
|
||||
|
||||
// Accumulators per shell
|
||||
struct ShellAccum {
|
||||
int total_obs = 0;
|
||||
std::unordered_set<int> unique_hkls;
|
||||
double rmeas_num = 0.0;
|
||||
double rmeas_den = 0.0;
|
||||
double sum_i_over_sig = 0.0;
|
||||
int n_merged_with_sigma = 0;
|
||||
};
|
||||
std::vector<ShellAccum> shell_acc(kStatShells);
|
||||
|
||||
for (int h = 0; h < nhkl; ++h) {
|
||||
if (hkl_shell[h] < 0)
|
||||
continue;
|
||||
const int s = hkl_shell[h];
|
||||
auto &sa = shell_acc[s];
|
||||
const auto &hs = per_hkl[h];
|
||||
if (hs.n == 0)
|
||||
continue;
|
||||
|
||||
sa.unique_hkls.insert(h);
|
||||
sa.total_obs += hs.n;
|
||||
|
||||
const double mean_I = hs.sum_I / static_cast<double>(hs.n);
|
||||
|
||||
if (hs.n >= 2) {
|
||||
double sum_abs_dev = 0.0;
|
||||
for (double Ii: hs.I_list)
|
||||
sum_abs_dev += std::abs(Ii - mean_I);
|
||||
sa.rmeas_num += std::sqrt(static_cast<double>(hs.n) / (hs.n - 1.0)) * sum_abs_dev;
|
||||
}
|
||||
|
||||
for (double Ii: hs.I_list)
|
||||
sa.rmeas_den += std::abs(Ii);
|
||||
|
||||
if (out.merged[h].sigma > 0.0) {
|
||||
sa.sum_i_over_sig += out.merged[h].I / out.merged[h].sigma;
|
||||
sa.n_merged_with_sigma += 1;
|
||||
}
|
||||
}
|
||||
|
||||
// Completeness (not yet available without unit cell)
|
||||
std::vector<int> possible_per_shell(kStatShells, 0);
|
||||
int total_possible = 0;
|
||||
bool have_completeness = false;
|
||||
|
||||
// Fill output statistics
|
||||
out.statistics.shells.resize(kStatShells);
|
||||
for (int s = 0; s < kStatShells; ++s) {
|
||||
auto &ss = out.statistics.shells[s];
|
||||
const auto &sa = shell_acc[s];
|
||||
|
||||
ss.mean_one_over_d2 = shell_mean_1_d2[s];
|
||||
ss.d_min = shell_min_res[s];
|
||||
ss.d_max = (s == 0) ? d_max_pad : shell_min_res[s - 1];
|
||||
|
||||
ss.total_observations = sa.total_obs;
|
||||
ss.unique_reflections = static_cast<int>(sa.unique_hkls.size());
|
||||
ss.rmeas = (sa.rmeas_den > 0.0) ? (sa.rmeas_num / sa.rmeas_den) : 0.0;
|
||||
ss.mean_i_over_sigma = (sa.n_merged_with_sigma > 0)
|
||||
? (sa.sum_i_over_sig / sa.n_merged_with_sigma)
|
||||
: 0.0;
|
||||
ss.possible_reflections = possible_per_shell[s];
|
||||
ss.completeness = (have_completeness && possible_per_shell[s] > 0)
|
||||
? static_cast<double>(sa.unique_hkls.size()) / possible_per_shell[s]
|
||||
: 0.0;
|
||||
}
|
||||
|
||||
// Overall statistics
|
||||
{
|
||||
auto &ov = out.statistics.overall;
|
||||
ov.d_min = stat_d_min;
|
||||
ov.d_max = stat_d_max;
|
||||
ov.mean_one_over_d2 = 0.0f;
|
||||
|
||||
int total_obs_all = 0;
|
||||
std::unordered_set<int> all_unique;
|
||||
double rmeas_num_all = 0.0, rmeas_den_all = 0.0;
|
||||
double sum_isig_all = 0.0;
|
||||
int n_isig_all = 0;
|
||||
|
||||
for (int s = 0; s < kStatShells; ++s) {
|
||||
const auto &sa = shell_acc[s];
|
||||
total_obs_all += sa.total_obs;
|
||||
all_unique.insert(sa.unique_hkls.begin(), sa.unique_hkls.end());
|
||||
rmeas_num_all += sa.rmeas_num;
|
||||
rmeas_den_all += sa.rmeas_den;
|
||||
sum_isig_all += sa.sum_i_over_sig;
|
||||
n_isig_all += sa.n_merged_with_sigma;
|
||||
}
|
||||
|
||||
ov.total_observations = total_obs_all;
|
||||
ov.unique_reflections = static_cast<int>(all_unique.size());
|
||||
ov.rmeas = (rmeas_den_all > 0.0) ? (rmeas_num_all / rmeas_den_all) : 0.0;
|
||||
ov.mean_i_over_sigma = (n_isig_all > 0) ? (sum_isig_all / n_isig_all) : 0.0;
|
||||
ov.possible_reflections = total_possible;
|
||||
ov.completeness = (have_completeness && total_possible > 0)
|
||||
? static_cast<double>(all_unique.size()) / total_possible
|
||||
: 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
void calc_obs(const ScaleMergeOptions &opt,
|
||||
std::vector<double> &g,
|
||||
std::vector<double> &mosaicity,
|
||||
std::vector<double> &R_sq,
|
||||
const std::vector<ObsRef> &obs,
|
||||
std::vector<CorrectedObs> &corr_obs) {
|
||||
|
||||
// ---- Compute corrected observations once (used for both merging and statistics) ----
|
||||
const double half_wedge = opt.wedge_deg.value_or(0.0) / 2.0;
|
||||
|
||||
for (const auto &o: obs) {
|
||||
const Reflection &r = *o.r;
|
||||
const double lp = SafeInv(static_cast<double>(r.rlp), 1.0);
|
||||
const double G_i = g[o.img_id];
|
||||
|
||||
// Compute partiality with refined mosaicity
|
||||
double partiality = 1.0;
|
||||
|
||||
switch (opt.partiality_model) {
|
||||
case ScaleMergeOptions::PartialityModel::Fixed:
|
||||
partiality = r.partiality;
|
||||
break;
|
||||
case ScaleMergeOptions::PartialityModel::Rotation: {
|
||||
const double c1 = r.zeta / std::sqrt(2.0);
|
||||
const double arg_plus = (r.delta_phi_deg + half_wedge) * c1 / mosaicity[o.img_id];
|
||||
const double arg_minus = (r.delta_phi_deg - half_wedge) * c1 / mosaicity[o.img_id];
|
||||
partiality = (std::erf(arg_plus) - std::erf(arg_minus)) / 2.0;
|
||||
}
|
||||
break;
|
||||
case ScaleMergeOptions::PartialityModel::Still:
|
||||
partiality = std::exp(-r.dist_ewald * r.dist_ewald / R_sq[o.img_id]);
|
||||
break;
|
||||
case ScaleMergeOptions::PartialityModel::Unity:
|
||||
break;
|
||||
}
|
||||
|
||||
if (partiality <= opt.min_partiality_for_merge)
|
||||
continue;
|
||||
|
||||
const double correction = G_i * partiality * lp;
|
||||
if (correction <= 0.0)
|
||||
continue;
|
||||
|
||||
corr_obs.push_back({
|
||||
o.hkl_slot,
|
||||
static_cast<double>(r.I) / correction,
|
||||
o.sigma / correction,
|
||||
1 / correction
|
||||
});
|
||||
}
|
||||
}
|
||||
|
||||
void proc_obs(const std::vector<std::vector<Reflection> > &observations,
|
||||
const ScaleMergeOptions &opt,
|
||||
std::vector<uint8_t> &image_slot_used,
|
||||
std::vector<ObsRef> &obs,
|
||||
std::unordered_map<HKLKey, int, HKLKeyHash> &hklToSlot
|
||||
) {
|
||||
for (int i = 0; i < observations.size(); i++) {
|
||||
for (const auto &r: observations[i]) {
|
||||
const double d = SafeD(r.d);
|
||||
if (!std::isfinite(d))
|
||||
continue;
|
||||
if (!std::isfinite(r.I))
|
||||
continue;
|
||||
|
||||
if (opt.d_min_limit_A > 0.0 && d < opt.d_min_limit_A)
|
||||
continue;
|
||||
|
||||
if (!std::isfinite(r.zeta) || r.zeta <= 0.0f)
|
||||
continue;
|
||||
if (!std::isfinite(r.rlp) || r.rlp == 0.0f)
|
||||
continue;
|
||||
|
||||
const double sigma = SafeSigma(r.sigma, opt.min_sigma);
|
||||
|
||||
const int img_id = i / opt.image_cluster;
|
||||
image_slot_used[img_id] = 1;
|
||||
|
||||
int hkl_slot;
|
||||
try {
|
||||
const HKLKey key = CanonicalizeHKLKey(r, opt);
|
||||
auto it = hklToSlot.find(key);
|
||||
if (it == hklToSlot.end()) {
|
||||
hkl_slot = static_cast<int>(hklToSlot.size());
|
||||
hklToSlot.emplace(key, hkl_slot);
|
||||
} else {
|
||||
hkl_slot = it->second;
|
||||
}
|
||||
} catch (...) {
|
||||
continue;
|
||||
}
|
||||
|
||||
ObsRef o;
|
||||
o.r = &r;
|
||||
o.img_id = img_id;
|
||||
o.hkl_slot = hkl_slot;
|
||||
o.sigma = sigma;
|
||||
obs.push_back(o);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
} // namespace
|
||||
|
||||
ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Reflection> > &observations,
|
||||
const ScaleMergeOptions &opt) {
|
||||
if (opt.image_cluster <= 0)
|
||||
throw std::invalid_argument("image_cluster must be positive");
|
||||
|
||||
const bool rotation_crystallography = opt.wedge_deg.has_value();
|
||||
|
||||
size_t nrefl = 0;
|
||||
for (const auto &i: observations)
|
||||
nrefl += i.size();
|
||||
|
||||
std::vector<ObsRef> obs;
|
||||
std::vector<CorrectedObs> corr_obs;
|
||||
|
||||
obs.reserve(nrefl);
|
||||
corr_obs.reserve(nrefl);
|
||||
|
||||
std::unordered_map<HKLKey, int, HKLKeyHash> hklToSlot;
|
||||
hklToSlot.reserve(nrefl);
|
||||
|
||||
size_t n_image_slots = observations.size() / opt.image_cluster +
|
||||
(observations.size() % opt.image_cluster > 0 ? 1 : 0);
|
||||
|
||||
std::vector<uint8_t> image_slot_used(n_image_slots, 0);
|
||||
|
||||
proc_obs(observations, opt, image_slot_used, obs, hklToSlot);
|
||||
|
||||
const int nhkl = static_cast<int>(hklToSlot.size());
|
||||
|
||||
std::vector<double> g(n_image_slots, 1.0);
|
||||
std::vector<double> mosaicity(n_image_slots, opt.mosaicity_init_deg);
|
||||
std::vector<double> R_sq(n_image_slots, 0.001 * 0.001);
|
||||
|
||||
for (int i = 0; i < n_image_slots; i++) {
|
||||
if (!image_slot_used[i]) {
|
||||
mosaicity[i] = NAN;
|
||||
g[i] = NAN;
|
||||
R_sq[i] = NAN;
|
||||
} else if (opt.mosaicity_init_deg_vec.size() > i && std::isfinite(opt.mosaicity_init_deg_vec[i])) {
|
||||
mosaicity[i] = opt.mosaicity_init_deg_vec[i];
|
||||
}
|
||||
}
|
||||
|
||||
scale(opt, g, mosaicity, R_sq, image_slot_used, rotation_crystallography, nhkl, obs);
|
||||
|
||||
ScaleMergeResult out;
|
||||
|
||||
out.image_scale_g.resize(observations.size(), NAN);
|
||||
out.mosaicity_deg.resize(observations.size(), NAN);
|
||||
for (int i = 0; i < observations.size(); i++) {
|
||||
size_t img_slot = i / opt.image_cluster;
|
||||
if (image_slot_used[img_slot]) {
|
||||
out.image_scale_g[i] = g[img_slot];
|
||||
out.mosaicity_deg[i] = mosaicity[img_slot];
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<HKLKey> slotToHKL(nhkl);
|
||||
for (const auto &kv: hklToSlot)
|
||||
slotToHKL[kv.second] = kv.first;
|
||||
|
||||
out.merged.resize(nhkl);
|
||||
for (int h = 0; h < nhkl; ++h) {
|
||||
out.merged[h].h = slotToHKL[h].h;
|
||||
out.merged[h].k = slotToHKL[h].k;
|
||||
out.merged[h].l = slotToHKL[h].l;
|
||||
out.merged[h].I = 0.0;
|
||||
out.merged[h].sigma = 0.0;
|
||||
out.merged[h].d = 0.0;
|
||||
}
|
||||
|
||||
// Populate d from median of observations per HKL
|
||||
{
|
||||
std::vector<std::vector<double> > per_hkl_d(nhkl);
|
||||
for (const auto &o: obs) {
|
||||
const double d_val = static_cast<double>(o.r->d);
|
||||
if (std::isfinite(d_val) && d_val > 0.0)
|
||||
per_hkl_d[o.hkl_slot].push_back(d_val);
|
||||
}
|
||||
for (int h = 0; h < nhkl; ++h) {
|
||||
auto &v = per_hkl_d[h];
|
||||
if (!v.empty()) {
|
||||
std::nth_element(v.begin(), v.begin() + static_cast<long>(v.size() / 2), v.end());
|
||||
out.merged[h].d = v[v.size() / 2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
calc_obs(opt, g, mosaicity, R_sq, obs, corr_obs);
|
||||
merge(nhkl, out, corr_obs);
|
||||
stats(opt, nhkl, out, corr_obs);
|
||||
|
||||
return out;
|
||||
}
|
||||
@@ -1,94 +0,0 @@
|
||||
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <vector>
|
||||
#include <cstdint>
|
||||
#include <optional>
|
||||
|
||||
#include "../../common/Reflection.h"
|
||||
#include "gemmi/symmetry.hpp"
|
||||
|
||||
struct ScaleMergeOptions {
|
||||
int max_num_iterations = 100;
|
||||
double max_solver_time_s = 1.0;
|
||||
|
||||
double image_number_rounding = 1.0;
|
||||
double min_sigma = 1e-3;
|
||||
|
||||
// Symmetry canonicalization of HKL prior to merging/scaling.
|
||||
// If not set, the routine uses raw HKL as-is.
|
||||
std::optional<gemmi::SpaceGroup> space_group;
|
||||
|
||||
// If true, treat Friedel mates as equivalent (merge anomalous pairs).
|
||||
// If false, keep them separate by including a sign flag in the HKL key.
|
||||
bool merge_friedel = true;
|
||||
|
||||
// --- Kabsch(XDS)-style partiality model ---
|
||||
// Rotation range (wedge) used in partiality calculation.
|
||||
// Set to 0 to disable partiality correction.
|
||||
std::optional<double> wedge_deg;
|
||||
|
||||
// --- Mosaicity (user input in degrees; internally converted to radians) ---
|
||||
double mosaicity_init_deg = 0.17;
|
||||
double mosaicity_min_deg = 1e-3;
|
||||
double mosaicity_max_deg = 2.0;
|
||||
std::vector<double> mosaicity_init_deg_vec;
|
||||
|
||||
// --- Optional: regularize per-image scale k towards 1 (Kabsch-like) ---
|
||||
bool regularize_scale_to_one = true;
|
||||
double scale_regularization_sigma = 0.05;
|
||||
|
||||
double min_partiality_for_merge = 0.2;
|
||||
bool smoothen_g = true;
|
||||
bool smoothen_mos = true;
|
||||
|
||||
double d_min_limit_A = 0.0;
|
||||
|
||||
int64_t image_cluster = 1;
|
||||
|
||||
bool refine_wedge = false;
|
||||
|
||||
enum class PartialityModel {Fixed, Rotation, Unity, Still} partiality_model = PartialityModel::Fixed;
|
||||
};
|
||||
|
||||
struct MergedReflection {
|
||||
int h;
|
||||
int k;
|
||||
int l;
|
||||
double I;
|
||||
double sigma;
|
||||
double d = 0.0;
|
||||
};
|
||||
|
||||
/// Per-resolution-shell merging statistics
|
||||
struct MergeStatisticsShell {
|
||||
float d_min = 0.0f; ///< High-resolution limit of this shell (Å)
|
||||
float d_max = 0.0f; ///< Low-resolution limit of this shell (Å)
|
||||
float mean_one_over_d2 = 0; ///< Mean 1/d² in shell
|
||||
int total_observations = 0; ///< Total number of (corrected) observations
|
||||
int unique_reflections = 0; ///< Number of unique HKLs with observations
|
||||
double rmeas = 0.0; ///< Redundancy-independent merging R-factor
|
||||
double mean_i_over_sigma = 0.0; ///< Mean I/σ(I) of the merged reflections
|
||||
double completeness = 0.0; ///< Fraction of possible reflections observed (0 if unknown)
|
||||
int possible_reflections = 0;///< Theoretical number of reflections in this shell (0 if unknown)
|
||||
};
|
||||
|
||||
/// Overall + per-shell merging statistics
|
||||
struct MergeStatistics {
|
||||
std::vector<MergeStatisticsShell> shells;
|
||||
MergeStatisticsShell overall; ///< Statistics over all shells combined
|
||||
};
|
||||
|
||||
struct ScaleMergeResult {
|
||||
std::vector<MergedReflection> merged;
|
||||
std::vector<float> image_scale_g;
|
||||
std::vector<float> mosaicity_deg;
|
||||
|
||||
/// Per-shell and overall merging statistics (populated after merging)
|
||||
MergeStatistics statistics;
|
||||
};
|
||||
|
||||
ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Reflection>>& observations,
|
||||
const ScaleMergeOptions& opt = {});
|
||||
@@ -0,0 +1,13 @@
|
||||
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include "ScaleOnTheFly.h"
|
||||
|
||||
ScaleOnTheFly::ScaleOnTheFly(const std::vector<MergedReflection> &ref,
|
||||
std::optional<gemmi::SpaceGroup> sg,
|
||||
bool merge_friedel) : sg(sg), merge_friedel(merge_friedel) {
|
||||
for (const auto &r : ref) {
|
||||
const auto key = CanonicalHKL(r, merge_friedel, sg);
|
||||
reference_data[key] = r.I;
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,19 @@
|
||||
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#pragma once
|
||||
|
||||
#include "HKLKey.h"
|
||||
#include "Merge.h"
|
||||
|
||||
#include <map>
|
||||
|
||||
class ScaleOnTheFly {
|
||||
std::optional<gemmi::SpaceGroup> sg;
|
||||
bool merge_friedel;
|
||||
std::map<HKLKey, double> reference_data;
|
||||
public:
|
||||
ScaleOnTheFly(const std::vector<MergedReflection> &ref, std::optional<gemmi::SpaceGroup> sg, bool merge_friedel = true);
|
||||
void Scale(std::vector<Reflection> &reflections);
|
||||
};
|
||||
|
||||
@@ -244,10 +244,10 @@ namespace {
|
||||
|
||||
const int idx = bin_index(r.d);
|
||||
auto& bin = bins[static_cast<size_t>(idx)];
|
||||
bin.d_min_A = std::min(bin.d_min_A, r.d);
|
||||
bin.d_max_A = std::max(bin.d_max_A, r.d);
|
||||
bin.d_min_A = std::min<double>(bin.d_min_A, r.d);
|
||||
bin.d_max_A = std::max<double>(bin.d_max_A, r.d);
|
||||
|
||||
const double i_over_sigma = std::max(0.0, r.I / r.sigma);
|
||||
const double i_over_sigma = std::max<double>(0.0, r.I / r.sigma);
|
||||
const gemmi::Op::Miller hkl{{r.h, r.k, r.l}};
|
||||
|
||||
if (gops.is_systematically_absent(hkl)) {
|
||||
|
||||
@@ -7,7 +7,7 @@
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
#include "ScaleAndMerge.h"
|
||||
#include "Merge.h"
|
||||
#include "gemmi/symmetry.hpp"
|
||||
|
||||
struct SpaceGroupOperatorScore {
|
||||
|
||||
@@ -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<uint8_t> serialization_buffer(approx_size);
|
||||
CBORStream2Serializer serializer(serialization_buffer.data(), serialization_buffer.size());
|
||||
|
||||
@@ -27,12 +27,6 @@ void HDF5FilePusher::StartDataCollection(StartMessage &message) {
|
||||
StartMessage repub_message = message;
|
||||
repub_message.writer_notification_zmq_addr = "";
|
||||
|
||||
size_t approx_size = 1024 * 1024;
|
||||
for (const auto &[key, value] : repub_message.pixel_mask)
|
||||
approx_size += value.size() * sizeof(uint32_t);
|
||||
|
||||
std::vector<uint8_t> serialization_buffer(approx_size);
|
||||
CBORStream2Serializer serializer(serialization_buffer.data(), serialization_buffer.size());
|
||||
serializer.SerializeSequenceStart(repub_message);
|
||||
|
||||
repub_active = repub_socket->Send(serialization_buffer.data(), serializer.GetBufferSize(), true);
|
||||
@@ -59,9 +53,6 @@ bool HDF5FilePusher::EndDataCollection(const EndMessage &message) {
|
||||
|
||||
if (repub_socket) {
|
||||
try {
|
||||
size_t approx_size = 1024 * 1024;
|
||||
std::vector<uint8_t> serialization_buffer(approx_size);
|
||||
CBORStream2Serializer serializer(serialization_buffer.data(), serialization_buffer.size());
|
||||
serializer.SerializeSequenceEnd(message);
|
||||
|
||||
if (repub_active)
|
||||
|
||||
@@ -13,6 +13,10 @@ void PrepareCBORImage(DataMessage& message,
|
||||
experiment.GetCompressionAlgorithm());
|
||||
}
|
||||
|
||||
ImagePusher::ImagePusher()
|
||||
: serialization_buffer(MESSAGE_SIZE_FOR_START_END),
|
||||
serializer(serialization_buffer.data(), serialization_buffer.size()) {}
|
||||
|
||||
std::string ImagePusher::Finalize() {
|
||||
return "";
|
||||
}
|
||||
|
||||
@@ -10,6 +10,7 @@
|
||||
#include "../common/JFJochMessages.h"
|
||||
#include "../common/ZeroCopyReturnValue.h"
|
||||
#include "../common/Logger.h"
|
||||
#include "../frame_serialize/CBORStream2Serializer.h"
|
||||
|
||||
enum class ImagePusherType {HDF5, CBOR, TCP, ZMQ, Test, None};
|
||||
|
||||
@@ -33,6 +34,10 @@ void PrepareCBORImage(DataMessage& message,
|
||||
void *image, size_t image_size);
|
||||
|
||||
class ImagePusher {
|
||||
protected:
|
||||
std::vector<uint8_t> serialization_buffer;
|
||||
CBORStream2Serializer serializer;
|
||||
ImagePusher();
|
||||
public:
|
||||
virtual void StartDataCollection(StartMessage& message) = 0;
|
||||
virtual bool EndDataCollection(const EndMessage& message) = 0; // Non-blocking
|
||||
|
||||
@@ -108,9 +108,7 @@ void TCPStreamPusher::CloseFd(std::atomic<int>& fd) {
|
||||
TCPStreamPusher::TCPStreamPusher(const std::string& addr,
|
||||
size_t in_max_connections,
|
||||
std::optional<int32_t> in_send_buffer_size)
|
||||
: serialization_buffer(256 * 1024 * 1024),
|
||||
serializer(serialization_buffer.data(), serialization_buffer.size()),
|
||||
endpoint(addr),
|
||||
: endpoint(addr),
|
||||
max_connections(in_max_connections),
|
||||
send_buffer_size(in_send_buffer_size) {
|
||||
if (endpoint.empty())
|
||||
|
||||
@@ -97,9 +97,6 @@ class TCPStreamPusher : public ImagePusher {
|
||||
std::chrono::steady_clock::time_point last_keepalive_recv{};
|
||||
};
|
||||
|
||||
std::vector<uint8_t> serialization_buffer;
|
||||
CBORStream2Serializer serializer;
|
||||
|
||||
std::string endpoint;
|
||||
size_t max_connections;
|
||||
std::optional<int32_t> send_buffer_size;
|
||||
|
||||
@@ -7,8 +7,6 @@
|
||||
ZMQStream2Pusher::ZMQStream2Pusher(const std::vector<std::string> &addr,
|
||||
std::optional<int32_t> send_buffer_high_watermark,
|
||||
std::optional<int32_t> send_buffer_size)
|
||||
: serialization_buffer(256*1024*1024),
|
||||
serializer(serialization_buffer.data(), serialization_buffer.size())
|
||||
{
|
||||
if (addr.empty())
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No writer ZMQ address provided");
|
||||
|
||||
@@ -9,12 +9,8 @@
|
||||
#include "../preview/PreviewCounter.h"
|
||||
#include "ZMQWriterNotificationPuller.h"
|
||||
#include "ZMQStream2PusherSocket.h"
|
||||
#include "../frame_serialize/CBORStream2Serializer.h"
|
||||
|
||||
class ZMQStream2Pusher : public ImagePusher {
|
||||
std::vector<uint8_t> serialization_buffer;
|
||||
CBORStream2Serializer serializer;
|
||||
|
||||
std::vector<std::unique_ptr<ZMQStream2PusherSocket>> socket;
|
||||
|
||||
std::unique_ptr<ZMQWriterNotificationPuller> writer_notification_socket;
|
||||
|
||||
+440
-173
@@ -1,6 +1,7 @@
|
||||
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include <cmath>
|
||||
#include "JFJochHDF5Reader.h"
|
||||
#include "spdlog/fmt/fmt.h"
|
||||
#include "../image_analysis/bragg_integration/CalcISigma.h"
|
||||
@@ -55,7 +56,7 @@ inline std::pair<gemmi::CrystalSystem, char> parse_niggli_class(int64_t val) {
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<hsize_t> GetDimension(HDF5Object& object, const std::string& path) {
|
||||
std::vector<hsize_t> 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 <class T>
|
||||
template<class T>
|
||||
void JFJochHDF5Reader::ReadVector(std::vector<T> &v,
|
||||
HDF5Object &file,
|
||||
const std::string &dataset_name,
|
||||
@@ -106,33 +107,123 @@ void JFJochHDF5Reader::ReadVector(std::vector<T> &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,
|
||||
std::vector<Reflection> &reflections) {
|
||||
if (!file.Exists("/entry/reflections") || !file.Exists(image_group_name))
|
||||
return false;
|
||||
|
||||
auto h = file.ReadOptVector<int32_t>(image_group_name + "/h");
|
||||
auto k = file.ReadOptVector<int32_t>(image_group_name + "/k");
|
||||
auto l = file.ReadOptVector<int32_t>(image_group_name + "/l");
|
||||
auto predicted_x = file.ReadOptVector<float>(image_group_name + "/predicted_x");
|
||||
auto predicted_y = file.ReadOptVector<float>(image_group_name + "/predicted_y");
|
||||
|
||||
auto d = file.ReadOptVector<float>(image_group_name + "/d");
|
||||
auto int_sum = file.ReadOptVector<float>(image_group_name + "/int_sum");
|
||||
auto int_err = file.ReadOptVector<float>(image_group_name + "/int_err");
|
||||
auto bkg = file.ReadOptVector<float>(image_group_name + "/background_mean");
|
||||
auto lp = file.ReadOptVector<float>(image_group_name + "/lp");
|
||||
auto partiality = file.ReadOptVector<float>(image_group_name + "/partiality");
|
||||
auto phi = file.ReadOptVector<float>(image_group_name + "/delta_phi");
|
||||
auto zeta = file.ReadOptVector<float>(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
|
||||
};
|
||||
reflections.emplace_back(r);
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
template<class T>
|
||||
std::optional<T> 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<T>(path, master_image);
|
||||
if (source_file.Exists(path))
|
||||
return source_file.ReadElement<T>(path, source_image);
|
||||
return {};
|
||||
}
|
||||
|
||||
template<class T>
|
||||
std::vector<T> ReadVectorMasterFirst(HDF5Object &master_file,
|
||||
HDF5Object &source_file,
|
||||
const std::string &path,
|
||||
const std::vector<hsize_t> &master_start,
|
||||
const std::vector<hsize_t> &source_start,
|
||||
const std::vector<hsize_t> &size) {
|
||||
if (master_file.Exists(path))
|
||||
return master_file.ReadOptVector<T>(path, master_start, size);
|
||||
if (source_file.Exists(path))
|
||||
return source_file.ReadOptVector<T>(path, source_start, size);
|
||||
return {};
|
||||
}
|
||||
|
||||
void JFJochHDF5Reader::ReadFile(const std::string &filename) {
|
||||
std::unique_lock ul(hdf5_mutex);
|
||||
try {
|
||||
auto dataset = std::make_shared<JFJochReaderDataset>();
|
||||
master_file = std::make_shared<HDF5ReadOnlyFile>(filename);
|
||||
master_filename = filename;
|
||||
dataset->experiment = default_experiment;
|
||||
|
||||
std::filesystem::path master_path(filename);
|
||||
@@ -165,7 +256,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<float>(
|
||||
dataset->efficiency = master_file->ReadVector<float>(
|
||||
"/entry/instrument/detector/detectorSpecific/data_collection_efficiency_image");
|
||||
else
|
||||
dataset->efficiency = std::vector<float>(number_of_images, 1.0);
|
||||
@@ -198,6 +289,8 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) {
|
||||
dataset->profile_radius = master_file->ReadOptVector<float>("/entry/MX/profileRadius");
|
||||
dataset->mosaicity_deg = master_file->ReadOptVector<float>("/entry/MX/mosaicity");
|
||||
dataset->b_factor = master_file->ReadOptVector<float>("/entry/MX/bFactor");
|
||||
dataset->scale_factor = master_file->ReadOptVector<float>("/entry/MX/imageScaleFactor");
|
||||
dataset->integrated_reflections = master_file->ReadOptVector<float>("/entry/MX/integratedReflections");
|
||||
}
|
||||
if (master_file->Exists("/entry/image"))
|
||||
dataset->max_value = master_file->ReadOptVector<int64_t>("/entry/image/max_value");
|
||||
@@ -317,8 +410,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 +427,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 +447,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<float>("/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 +470,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 +488,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 +502,20 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) {
|
||||
auto tmp = master_file->ReadOptVector<float>("/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 +531,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 +560,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 +579,9 @@ void JFJochHDF5Reader::ReadFile(const std::string& filename) {
|
||||
|
||||
if (image_size_x * image_size_y > 0) {
|
||||
auto mask_tmp = master_file->ReadOptVector<uint32_t>(
|
||||
"/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<uint32_t>(image_size_x * image_size_y);
|
||||
@@ -490,7 +589,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 +630,14 @@ CompressedImage JFJochHDF5Reader::LoadImageDataset(std::vector<uint8_t> &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<std::shared_ptr<HDF5ReadOnlyFile>, uint32_t> JFJochHDF5Reader::GetImageLocation(int64_t image_number) {
|
||||
@@ -558,7 +660,7 @@ std::pair<std::shared_ptr<HDF5ReadOnlyFile>, uint32_t> JFJochHDF5Reader::GetImag
|
||||
|
||||
image_id = static_cast<uint32_t>(mapping.SourceImage(image));
|
||||
data_file = std::make_shared<HDF5ReadOnlyFile>(
|
||||
ResolveRelativeToMaster(master_file_directory, mapping.filename)
|
||||
ResolveRelativeToMaster(master_file_directory, mapping.filename)
|
||||
);
|
||||
return {std::move(data_file), image_id};
|
||||
}
|
||||
@@ -601,69 +703,106 @@ bool JFJochHDF5Reader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &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<uint32_t>("/entry/MX/nPeaks", image_id);
|
||||
const auto master_image = static_cast<hsize_t>(image_number);
|
||||
const auto source_image = static_cast<hsize_t>(image_id);
|
||||
|
||||
auto spot_count_opt = ReadElementMasterFirst<uint32_t>(*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<float>(
|
||||
"/entry/MX/peakXPosRaw",
|
||||
{(hsize_t) image_id, 0},
|
||||
{1, spot_count}
|
||||
auto spot_x = ReadVectorMasterFirst<float>(
|
||||
*master_file,
|
||||
*source_file,
|
||||
"/entry/MX/peakXPosRaw",
|
||||
{master_image, 0},
|
||||
{source_image, 0},
|
||||
{1, spot_count}
|
||||
);
|
||||
auto spot_y = source_file->ReadVector<float>(
|
||||
"/entry/MX/peakYPosRaw",
|
||||
{(hsize_t) image_id, 0},
|
||||
{1, spot_count}
|
||||
auto spot_y = ReadVectorMasterFirst<float>(
|
||||
*master_file,
|
||||
*source_file,
|
||||
"/entry/MX/peakYPosRaw",
|
||||
{master_image, 0},
|
||||
{source_image, 0},
|
||||
{1, spot_count}
|
||||
);
|
||||
auto spot_intensity = source_file->ReadVector<float>(
|
||||
"/entry/MX/peakTotalIntensity",
|
||||
{(hsize_t) image_id, 0},
|
||||
{1, spot_count}
|
||||
);
|
||||
auto spot_indexed = source_file->ReadOptVector<uint8_t>(
|
||||
"/entry/MX/peakIndexed",
|
||||
{(hsize_t) image_id, 0},
|
||||
{1, spot_count}
|
||||
auto spot_intensity = ReadVectorMasterFirst<float>(
|
||||
*master_file,
|
||||
*source_file,
|
||||
"/entry/MX/peakTotalIntensity",
|
||||
{master_image, 0},
|
||||
{source_image, 0},
|
||||
{1, spot_count}
|
||||
);
|
||||
|
||||
auto spot_ice = source_file->ReadOptVector<uint8_t>(
|
||||
"/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<uint8_t>(
|
||||
*master_file,
|
||||
*source_file,
|
||||
"/entry/MX/peakIndexed",
|
||||
{master_image, 0},
|
||||
{source_image, 0},
|
||||
{1, spot_count}
|
||||
);
|
||||
|
||||
auto spot_h = source_file->ReadOptVector<int32_t>(
|
||||
"/entry/MX/peakH",
|
||||
{(hsize_t) image_id, 0},
|
||||
{1, spot_count}
|
||||
auto spot_ice = ReadVectorMasterFirst<uint8_t>(
|
||||
*master_file,
|
||||
*source_file,
|
||||
"/entry/MX/peakIceRingRes",
|
||||
{master_image, 0},
|
||||
{source_image, 0},
|
||||
{1, spot_count}
|
||||
);
|
||||
|
||||
auto spot_k = source_file->ReadOptVector<int32_t>(
|
||||
"/entry/MX/peakK",
|
||||
{(hsize_t) image_id, 0},
|
||||
{1, spot_count}
|
||||
auto spot_h = ReadVectorMasterFirst<int32_t>(
|
||||
*master_file,
|
||||
*source_file,
|
||||
"/entry/MX/peakH",
|
||||
{master_image, 0},
|
||||
{source_image, 0},
|
||||
{1, spot_count}
|
||||
);
|
||||
|
||||
auto spot_l = source_file->ReadOptVector<int32_t>(
|
||||
"/entry/MX/peakL",
|
||||
{(hsize_t) image_id, 0},
|
||||
{1, spot_count}
|
||||
auto spot_k = ReadVectorMasterFirst<int32_t>(
|
||||
*master_file,
|
||||
*source_file,
|
||||
"/entry/MX/peakK",
|
||||
{master_image, 0},
|
||||
{source_image, 0},
|
||||
{1, spot_count}
|
||||
);
|
||||
|
||||
auto spot_dist_ewald_sphere = source_file->ReadOptVector<float>(
|
||||
"/entry/MX/peakDistEwaldSphere",
|
||||
{(hsize_t) image_id, 0},
|
||||
{1, spot_count}
|
||||
auto spot_l = ReadVectorMasterFirst<int32_t>(
|
||||
*master_file,
|
||||
*source_file,
|
||||
"/entry/MX/peakL",
|
||||
{master_image, 0},
|
||||
{source_image, 0},
|
||||
{1, spot_count}
|
||||
);
|
||||
|
||||
auto spot_dist_ewald_sphere = ReadVectorMasterFirst<float>(
|
||||
*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 +828,47 @@ bool JFJochHDF5Reader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &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<uint32_t>("/entry/MX/peakCountUnfiltered", image_id);
|
||||
|
||||
if (auto v = ReadElementMasterFirst<uint32_t>(*master_file, *source_file,
|
||||
"/entry/MX/peakCountUnfiltered",
|
||||
master_image, source_image); v)
|
||||
message.spot_count = v;
|
||||
else
|
||||
message.spot_count = source_file->ReadElement<uint32_t>("/entry/MX/nPeaks", image_id);
|
||||
if (source_file->Exists("/entry/MX/peakCountIceRingRes"))
|
||||
message.spot_count_ice_rings = source_file->ReadElement<uint32_t>("/entry/MX/peakCountIceRingRes", image_id);
|
||||
if (source_file->Exists("/entry/MX/peakCountLowRes"))
|
||||
message.spot_count_low_res = source_file->ReadElement<uint32_t>("/entry/MX/peakCountLowRes", image_id);
|
||||
if (source_file->Exists("/entry/MX/peakCountIndexed"))
|
||||
message.spot_count_indexed = source_file->ReadElement<uint32_t>("/entry/MX/peakCountIndexed", image_id);
|
||||
message.spot_count = spot_count_opt;
|
||||
|
||||
message.spot_count_ice_rings = ReadElementMasterFirst<uint32_t>(
|
||||
*master_file, *source_file, "/entry/MX/peakCountIceRingRes", master_image, source_image);
|
||||
message.spot_count_low_res = ReadElementMasterFirst<uint32_t>(
|
||||
*master_file, *source_file, "/entry/MX/peakCountLowRes", master_image, source_image);
|
||||
message.spot_count_indexed = ReadElementMasterFirst<uint32_t>(
|
||||
*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<float>("/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<float>(
|
||||
if (dataset->azimuthal_bins == 0) {
|
||||
message.az_int_profile = ReadVectorMasterFirst<float>(
|
||||
*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<float>(
|
||||
} else {
|
||||
message.az_int_profile = ReadVectorMasterFirst<float>(
|
||||
*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<int64_t>(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 +886,30 @@ bool JFJochHDF5Reader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &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<float> tmp = ReadVectorMasterFirst<float>(
|
||||
*master_file,
|
||||
*source_file,
|
||||
"/entry/MX/latticeIndexed",
|
||||
{master_image, 0},
|
||||
{source_image, 0},
|
||||
{1, 9}
|
||||
);
|
||||
|
||||
std::vector<float> tmp = source_file->ReadVector<float>(
|
||||
"/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<uint32_t> niggli_opt;
|
||||
if (master_file->Exists("/entry/MX/niggli_class"))
|
||||
niggli_opt = master_file->ReadElement<uint32_t>("/entry/MX/niggli_class", image_number);
|
||||
else if (master_file->Exists("/entry/MX/niggliClass"))
|
||||
niggli_opt = master_file->ReadElement<uint32_t>("/entry/MX/niggliClass", image_number);
|
||||
else if (source_file->Exists("/entry/MX/niggli_class"))
|
||||
niggli_opt = source_file->ReadElement<uint32_t>("/entry/MX/niggli_class", image_id);
|
||||
else if (source_file->Exists("/entry/MX/niggliClass"))
|
||||
niggli_opt = source_file->ReadElement<uint32_t>("/entry/MX/niggliClass", image_id);
|
||||
|
||||
std::optional<uint32_t> niggli_opt = source_file->ReadElement<uint32_t>("/entry/MX/niggli_class", image_id);
|
||||
if (niggli_opt) {
|
||||
auto symm_info = parse_niggli_class(niggli_opt.value());
|
||||
|
||||
@@ -762,66 +921,15 @@ bool JFJochHDF5Reader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &dataset
|
||||
}
|
||||
}
|
||||
|
||||
std::string image_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<int32_t>(image_group_name + "/h");
|
||||
auto k = source_file->ReadOptVector<int32_t>(image_group_name + "/k");
|
||||
auto l = source_file->ReadOptVector<int32_t>(image_group_name + "/l");
|
||||
auto predicted_x = source_file->ReadOptVector<float>(image_group_name + "/predicted_x");
|
||||
auto predicted_y = source_file->ReadOptVector<float>(image_group_name + "/predicted_y");
|
||||
|
||||
auto d = source_file->ReadOptVector<float>(image_group_name + "/d");
|
||||
auto int_sum = source_file->ReadOptVector<float>(image_group_name + "/int_sum");
|
||||
auto int_err = source_file->ReadOptVector<float>(image_group_name + "/int_err");
|
||||
auto bkg = source_file->ReadOptVector<float>(image_group_name + "/background_mean");
|
||||
auto lp = source_file->ReadOptVector<float>(image_group_name + "/lp");
|
||||
auto partiality = source_file->ReadOptVector<float>(image_group_name + "/partiality");
|
||||
auto phi = source_file->ReadOptVector<float>(image_group_name + "/delta_phi");
|
||||
auto zeta = source_file->ReadOptVector<float>(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);
|
||||
}
|
||||
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 (!ReadReflectionsFromGroup(*master_file, master_reflection_group_name, message.reflections))
|
||||
ReadReflectionsFromGroup(*source_file, source_reflection_group_name, message.reflections);
|
||||
if (!message.reflections.empty()) {
|
||||
CalcISigma(message);
|
||||
CalcWilsonBFactor(message, !message.b_factor.has_value());
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
@@ -832,6 +940,7 @@ void JFJochHDF5Reader::Close() {
|
||||
legacy_format_files.clear();
|
||||
vds_data_mappings.clear();
|
||||
master_file_directory.clear();
|
||||
master_filename.clear();
|
||||
SetStartMessage({});
|
||||
}
|
||||
|
||||
@@ -892,7 +1001,165 @@ CompressedImage JFJochHDF5Reader::ReadCalibration(std::vector<uint8_t> &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
|
||||
};
|
||||
}
|
||||
|
||||
void AppendOrExtendSourceMapping(std::vector<HDF5DataSourceMessage> &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<HDF5DataSourceMessage> JFJochHDF5Reader::GetHDF5DataSource(
|
||||
uint64_t first_image,
|
||||
std::optional<uint64_t> 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<HDF5DataSourceMessage> 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");
|
||||
}
|
||||
|
||||
std::vector<std::vector<Reflection> > JFJochHDF5Reader::ReadReflections(size_t start_image, std::optional<size_t> end_image) const
|
||||
{
|
||||
std::unique_lock ul(hdf5_mutex);
|
||||
if (start_image >= number_of_images)
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"start_image must be less than number_of_images");
|
||||
|
||||
if (end_image.has_value() && end_image.value() < start_image)
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"end_image must be greater than start_image if provided");
|
||||
|
||||
int end_image_val = end_image.value_or(number_of_images);
|
||||
|
||||
if (end_image_val > number_of_images)
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"end_image_val must be less than or equal to number_of_images");
|
||||
if (!master_file)
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"Cannot read all reflections if file not loaded");
|
||||
|
||||
std::vector<std::vector<Reflection> > ret(end_image_val - start_image);
|
||||
|
||||
for (int i = 0; i < end_image_val - start_image; ++i)
|
||||
ReadReflectionsFromGroup(*master_file,
|
||||
fmt::format("/entry/reflections/image_{:06d}", start_image + i),
|
||||
ret[i]);
|
||||
|
||||
return ret;
|
||||
|
||||
}
|
||||
|
||||
@@ -16,6 +16,7 @@ class JFJochHDF5Reader : public JFJochReader {
|
||||
std::vector<std::string> legacy_format_files;
|
||||
std::vector<HDF5VirtualDatasetMapping> 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,16 @@ public:
|
||||
void ReadFile(const std::string& filename);
|
||||
|
||||
uint64_t GetNumberOfImages() const override;
|
||||
|
||||
void Close() override;
|
||||
|
||||
std::vector<HDF5DataSourceMessage> GetHDF5DataSource(
|
||||
uint64_t first_image = 0,
|
||||
std::optional<uint64_t> image_count = {}
|
||||
) const;
|
||||
|
||||
std::vector<std::vector<Reflection>> ReadReflections(size_t start_image = 0, std::optional<size_t> end_image = {}) const;
|
||||
|
||||
CompressedImage ReadCalibration(std::vector<uint8_t> &tmp, const std::string &name) const;
|
||||
|
||||
std::shared_ptr<JFJochReaderRawImage> GetRawImage(int64_t image_number) override;
|
||||
|
||||
@@ -182,6 +182,7 @@ std::shared_ptr<JFJochReaderDataset> 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);
|
||||
|
||||
@@ -40,6 +40,7 @@ struct JFJochReaderDataset {
|
||||
std::vector<float> profile_radius;
|
||||
std::vector<float> mosaicity_deg;
|
||||
std::vector<float> b_factor;
|
||||
std::vector<float> integrated_reflections;
|
||||
|
||||
std::vector<float> scale_factor;
|
||||
|
||||
|
||||
@@ -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)");
|
||||
|
||||
@@ -468,7 +468,9 @@ TEST_CASE("HDF5Writer", "[HDF5][Full]") {
|
||||
DiffractionExperiment x(DetJF4M());
|
||||
std::vector<SpotToSave> 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);
|
||||
|
||||
|
||||
+447
-82
@@ -4,6 +4,7 @@
|
||||
#include <catch2/catch_all.hpp>
|
||||
|
||||
#include "../common/DiffractionExperiment.h"
|
||||
#include "../common/ScanResultGenerator.h"
|
||||
#include "../writer/FileWriter.h"
|
||||
#include "../reader/JFJochHDF5Reader.h"
|
||||
#include "../compression/JFJochCompressor.h"
|
||||
@@ -39,13 +40,13 @@ TEST_CASE("HDF5DataType_ElemType","[HDF5]") {
|
||||
TEST_CASE("JFJochReader_MasterFile", "[HDF5][Full]") {
|
||||
DiffractionExperiment x(DetJF(1));
|
||||
|
||||
|
||||
x.FilePrefix("test08").ImagesPerTrigger(950).OverwriteExistingFiles(true);
|
||||
x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150)
|
||||
.IncidentEnergy_keV(WVL_1A_IN_KEV)
|
||||
.FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10))
|
||||
.SetUnitCell(UnitCell{.a= 10, .b= 20, .c= 30, .alpha= 90, .beta= 101, .gamma = 90});
|
||||
RegisterHDF5Filter();
|
||||
|
||||
{
|
||||
StartMessage start_message;
|
||||
x.FillMessage(start_message);
|
||||
@@ -210,6 +211,8 @@ TEST_CASE("JFJochReader_PixelMask", "[HDF5][Full]") {
|
||||
pixel_mask[x.GetPixelsNum() - 1] = 4;
|
||||
pixel_mask[0] = 256;
|
||||
|
||||
ScanResultGenerator generator(x);
|
||||
|
||||
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
|
||||
{
|
||||
StartMessage start_message;
|
||||
@@ -223,9 +226,11 @@ TEST_CASE("JFJochReader_PixelMask", "[HDF5][Full]") {
|
||||
message.number = 0;
|
||||
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
|
||||
REQUIRE_NOTHROW(generator.Add(message));
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = 1;
|
||||
generator.FillEndMessage(end_message);
|
||||
|
||||
file_set.WriteHDF5(end_message);
|
||||
|
||||
file_set.Finalize();
|
||||
@@ -263,6 +268,7 @@ TEST_CASE("JFJochReader_Goniometer", "[HDF5][Full]") {
|
||||
|
||||
RegisterHDF5Filter();
|
||||
|
||||
ScanResultGenerator generator(x);
|
||||
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
|
||||
{
|
||||
StartMessage start_message;
|
||||
@@ -276,10 +282,12 @@ TEST_CASE("JFJochReader_Goniometer", "[HDF5][Full]") {
|
||||
message.number = i;
|
||||
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
REQUIRE_NOTHROW(generator.Add(message));
|
||||
}
|
||||
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = 5;
|
||||
generator.FillEndMessage(end_message);
|
||||
file_set.WriteHDF5(end_message);
|
||||
|
||||
file_set.Finalize();
|
||||
@@ -320,6 +328,7 @@ TEST_CASE("JFJochReader_GridScan", "[HDF5][Full]") {
|
||||
x.ImportDatasetSettings(d);
|
||||
|
||||
RegisterHDF5Filter();
|
||||
ScanResultGenerator generator(x);
|
||||
|
||||
std::vector<uint16_t> image(x.GetPixelsNum(), 0);
|
||||
{
|
||||
@@ -334,10 +343,12 @@ TEST_CASE("JFJochReader_GridScan", "[HDF5][Full]") {
|
||||
message.number = i;
|
||||
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
generator.Add(message);
|
||||
}
|
||||
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = 5;
|
||||
generator.FillEndMessage(end_message);
|
||||
file_set.WriteHDF5(end_message);
|
||||
|
||||
file_set.Finalize();
|
||||
@@ -378,6 +389,8 @@ TEST_CASE("JFJochReader_DataI16", "[HDF5][Full]") {
|
||||
image[2] = 456;
|
||||
image[3] = -3456;
|
||||
|
||||
ScanResultGenerator generator(x);
|
||||
|
||||
RegisterHDF5Filter();
|
||||
{
|
||||
StartMessage start_message;
|
||||
@@ -395,12 +408,13 @@ TEST_CASE("JFJochReader_DataI16", "[HDF5][Full]") {
|
||||
message.bkg_estimate = i * 345.6;
|
||||
message.number = i;
|
||||
message.profile_radius = 123.09;
|
||||
|
||||
generator.Add(message);
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
|
||||
}
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = x.GetImageNum();
|
||||
generator.FillEndMessage(end_message);
|
||||
file_set.WriteHDF5(end_message);
|
||||
file_set.Finalize();
|
||||
}
|
||||
@@ -529,6 +543,8 @@ TEST_CASE("JFJochReader_DataU16", "[HDF5][Full]") {
|
||||
image[1] = INT16_MAX;
|
||||
image[2] = 456;
|
||||
|
||||
ScanResultGenerator generator(x);
|
||||
|
||||
RegisterHDF5Filter();
|
||||
{
|
||||
StartMessage start_message;
|
||||
@@ -547,10 +563,11 @@ TEST_CASE("JFJochReader_DataU16", "[HDF5][Full]") {
|
||||
message.number = i;
|
||||
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
|
||||
generator.Add(message);
|
||||
}
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = x.GetImageNum();
|
||||
generator.FillEndMessage(end_message);
|
||||
file_set.WriteHDF5(end_message);
|
||||
file_set.Finalize();
|
||||
}
|
||||
@@ -597,6 +614,8 @@ TEST_CASE("JFJochReader_DataI32", "[HDF5][Full]") {
|
||||
image[1] = INT32_MIN;
|
||||
image[2] = 456;
|
||||
|
||||
ScanResultGenerator generator(x);
|
||||
|
||||
RegisterHDF5Filter();
|
||||
{
|
||||
StartMessage start_message;
|
||||
@@ -613,10 +632,11 @@ TEST_CASE("JFJochReader_DataI32", "[HDF5][Full]") {
|
||||
message.number = i;
|
||||
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
|
||||
generator.Add(message);
|
||||
}
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = x.GetImageNum();
|
||||
generator.FillEndMessage(end_message);
|
||||
file_set.WriteHDF5(end_message);
|
||||
file_set.Finalize();
|
||||
}
|
||||
@@ -663,6 +683,8 @@ TEST_CASE("JFJochReader_DataU32", "[HDF5][Full]") {
|
||||
image[3] = INT32_MAX;
|
||||
image[4] = INT32_MAX - 1;
|
||||
|
||||
ScanResultGenerator generator(x);
|
||||
|
||||
RegisterHDF5Filter();
|
||||
{
|
||||
StartMessage start_message;
|
||||
@@ -678,11 +700,12 @@ TEST_CASE("JFJochReader_DataU32", "[HDF5][Full]") {
|
||||
message.spots = spots;
|
||||
message.number = i;
|
||||
|
||||
generator.Add(message);
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
|
||||
}
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = x.GetImageNum();
|
||||
generator.FillEndMessage(end_message);
|
||||
file_set.WriteHDF5(end_message);
|
||||
file_set.Finalize();
|
||||
}
|
||||
@@ -731,6 +754,8 @@ TEST_CASE("JFJochReader_Summation", "[HDF5][Full]") {
|
||||
image_3[0] = INT16_MAX;
|
||||
image_2[1] = INT16_MIN;
|
||||
|
||||
ScanResultGenerator generator(x);
|
||||
|
||||
RegisterHDF5Filter();
|
||||
{
|
||||
StartMessage start_message;
|
||||
@@ -745,17 +770,21 @@ TEST_CASE("JFJochReader_Summation", "[HDF5][Full]") {
|
||||
message.image = CompressedImage(image_1, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.number = 0;
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
generator.Add(message);
|
||||
|
||||
message.image = CompressedImage(image_2, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.number = 1;
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
generator.Add(message);
|
||||
|
||||
message.image = CompressedImage(image_3, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.number = 2;
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
generator.Add(message);
|
||||
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = x.GetImageNum();
|
||||
generator.FillEndMessage(end_message);
|
||||
file_set.WriteHDF5(end_message);
|
||||
file_set.Finalize();
|
||||
}
|
||||
@@ -800,6 +829,8 @@ TEST_CASE("JFJochReader_Summation_5", "[HDF5][Full]") {
|
||||
image_3[0] = INT16_MAX;
|
||||
image_2[1] = INT16_MIN;
|
||||
|
||||
ScanResultGenerator generator(x);
|
||||
|
||||
RegisterHDF5Filter();
|
||||
{
|
||||
StartMessage start_message;
|
||||
@@ -814,25 +845,31 @@ TEST_CASE("JFJochReader_Summation_5", "[HDF5][Full]") {
|
||||
message.image = CompressedImage(image_1, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.number = 0;
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
generator.Add(message);
|
||||
|
||||
message.image = CompressedImage(image_2, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.number = 1;
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
generator.Add(message);
|
||||
|
||||
message.image = CompressedImage(image_3, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.number = 2;
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
generator.Add(message);
|
||||
|
||||
message.image = CompressedImage(image_4, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.number = 3;
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
generator.Add(message);
|
||||
|
||||
message.image = CompressedImage(image_5, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.number = 4;
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
generator.Add(message);
|
||||
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = x.GetImageNum();
|
||||
generator.FillEndMessage(end_message);
|
||||
file_set.WriteHDF5(end_message);
|
||||
file_set.Finalize();
|
||||
}
|
||||
@@ -860,82 +897,6 @@ TEST_CASE("JFJochReader_Summation_5", "[HDF5][Full]") {
|
||||
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
|
||||
}
|
||||
|
||||
TEST_CASE("JFJochReader_ROI", "[HDF5][Full]") {
|
||||
DiffractionExperiment x(DetJF(1));
|
||||
|
||||
x.FilePrefix("test25").ImagesPerTrigger(4).OverwriteExistingFiles(true);
|
||||
x.BitDepthImage(16).ImagesPerFile(1).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(false);
|
||||
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
|
||||
|
||||
x.ROI().SetROI(ROIDefinition{
|
||||
.boxes = {ROIBox("beam", 100, 120, 20, 30)},
|
||||
.circles = {ROICircle("roi1", 500, 800, 10)}
|
||||
});
|
||||
|
||||
std::vector<uint16_t> image(x.GetPixelsNum());
|
||||
|
||||
RegisterHDF5Filter();
|
||||
{
|
||||
StartMessage start_message;
|
||||
x.FillMessage(start_message);
|
||||
FileWriter file_set(start_message);
|
||||
|
||||
for (int i = 0; i < x.GetImageNum(); i++) {
|
||||
std::vector<SpotToSave> spots;
|
||||
|
||||
DataMessage message{};
|
||||
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.spots = spots;
|
||||
message.number = i;
|
||||
message.roi["beam"] = ROIMessage{
|
||||
.sum = 12 + i,
|
||||
.sum_square = (uint64_t) 123 * i,
|
||||
.max_count = 123 - i,
|
||||
.pixels = (uint64_t) 189 + i
|
||||
};
|
||||
|
||||
message.roi["roi1"] = ROIMessage{
|
||||
.sum = 25 + i,
|
||||
.sum_square = (uint64_t) 15 * i,
|
||||
.max_count = 67 - i,
|
||||
.pixels = (uint64_t) 95 + i
|
||||
};
|
||||
REQUIRE_NOTHROW(file_set.WriteHDF5(message));
|
||||
|
||||
}
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = x.GetImageNum();
|
||||
file_set.WriteHDF5(end_message);
|
||||
file_set.Finalize();
|
||||
}
|
||||
{
|
||||
JFJochHDF5Reader reader;
|
||||
REQUIRE_NOTHROW(reader.ReadFile("test25_master.h5"));
|
||||
auto dataset = reader.GetDataset();
|
||||
CHECK(dataset->experiment.GetImageNum() == 4);
|
||||
|
||||
CHECK(dataset->roi.size() == 2);
|
||||
REQUIRE(dataset->roi_max.size() == 2);
|
||||
int index = 0;
|
||||
if (dataset->roi[1] == "beam")
|
||||
index = 1;
|
||||
|
||||
REQUIRE(dataset->roi_max[index].size() == 4);
|
||||
CHECK(dataset->roi_max[index][2] == 123 - 2);
|
||||
CHECK(dataset->roi_sum_sq[index][0] == 0);
|
||||
CHECK(dataset->roi_sum_sq[index][1] == 123);
|
||||
CHECK(dataset->roi_npixel[index][3] == 189 + 3);
|
||||
CHECK(dataset->roi_sum[1 - index][3] == 25 + 3);
|
||||
}
|
||||
remove("test25_master.h5");
|
||||
remove("test25_data_000001.h5");
|
||||
remove("test25_data_000002.h5");
|
||||
remove("test25_data_000003.h5");
|
||||
remove("test25_data_000004.h5");
|
||||
// No leftover HDF5 objects
|
||||
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
|
||||
}
|
||||
|
||||
TEST_CASE("JFJochReader_Azint", "[HDF5][Full]") {
|
||||
DiffractionExperiment x(DetJF(1));
|
||||
|
||||
@@ -1669,5 +1630,409 @@ 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<int16_t> image(x.GetPixelsNum(), 17);
|
||||
|
||||
RegisterHDF5Filter();
|
||||
{
|
||||
StartMessage start_message;
|
||||
x.FillMessage(start_message);
|
||||
FileWriter writer(start_message);
|
||||
|
||||
for (int i = 0; i < x.GetImageNum(); i++) {
|
||||
image[5678] = static_cast<int16_t>(100 + i);
|
||||
|
||||
DataMessage message{};
|
||||
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.number = i;
|
||||
REQUIRE_NOTHROW(writer.WriteHDF5(message));
|
||||
}
|
||||
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = x.GetImageNum();
|
||||
writer.WriteHDF5(end_message);
|
||||
writer.Finalize();
|
||||
}
|
||||
|
||||
{
|
||||
JFJochHDF5Reader reader;
|
||||
REQUIRE_NOTHROW(reader.ReadFile("source_integrated_master.h5"));
|
||||
|
||||
auto source = reader.GetHDF5DataSource(1, 3);
|
||||
REQUIRE(source.size() == 1);
|
||||
|
||||
CHECK(source[0].filename == "source_integrated_master.h5");
|
||||
CHECK(source[0].dataset == "/entry/data/data");
|
||||
CHECK(source[0].source_first_image == 1);
|
||||
CHECK(source[0].virtual_first_image == 0);
|
||||
CHECK(source[0].image_count == 3);
|
||||
}
|
||||
|
||||
remove("source_integrated_master.h5");
|
||||
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
|
||||
}
|
||||
|
||||
TEST_CASE("JFJochReader_HDF5DataSource_VDS", "[HDF5][Full]") {
|
||||
DiffractionExperiment x(DetJF(1));
|
||||
|
||||
x.FilePrefix("source_vds_mapping").ImagesPerTrigger(5).ImagesPerFile(2).OverwriteExistingFiles(true);
|
||||
x.BitDepthImage(16).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true);
|
||||
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
|
||||
|
||||
std::vector<int16_t> image(x.GetPixelsNum(), 21);
|
||||
|
||||
RegisterHDF5Filter();
|
||||
{
|
||||
StartMessage start_message;
|
||||
x.FillMessage(start_message);
|
||||
FileWriter writer(start_message);
|
||||
|
||||
for (int i = 0; i < x.GetImageNum(); i++) {
|
||||
image[5678] = static_cast<int16_t>(200 + i);
|
||||
|
||||
DataMessage message{};
|
||||
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.number = i;
|
||||
REQUIRE_NOTHROW(writer.WriteHDF5(message));
|
||||
}
|
||||
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = x.GetImageNum();
|
||||
writer.WriteHDF5(end_message);
|
||||
writer.Finalize();
|
||||
}
|
||||
|
||||
{
|
||||
JFJochHDF5Reader reader;
|
||||
REQUIRE_NOTHROW(reader.ReadFile("source_vds_mapping_master.h5"));
|
||||
|
||||
// Range crosses file boundary:
|
||||
// global images 1,2,3 map to:
|
||||
// data_000001 image 1
|
||||
// data_000002 images 0,1
|
||||
auto source = reader.GetHDF5DataSource(1, 3);
|
||||
REQUIRE(source.size() == 2);
|
||||
|
||||
CHECK(source[0].filename == "source_vds_mapping_data_000001.h5");
|
||||
CHECK(source[0].dataset == "/entry/data/data");
|
||||
CHECK(source[0].source_first_image == 1);
|
||||
CHECK(source[0].virtual_first_image == 0);
|
||||
CHECK(source[0].image_count == 1);
|
||||
|
||||
CHECK(source[1].filename == "source_vds_mapping_data_000002.h5");
|
||||
CHECK(source[1].dataset == "/entry/data/data");
|
||||
CHECK(source[1].source_first_image == 0);
|
||||
CHECK(source[1].virtual_first_image == 1);
|
||||
CHECK(source[1].image_count == 2);
|
||||
}
|
||||
|
||||
remove("source_vds_mapping_master.h5");
|
||||
remove("source_vds_mapping_data_000001.h5");
|
||||
remove("source_vds_mapping_data_000002.h5");
|
||||
remove("source_vds_mapping_data_000003.h5");
|
||||
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
|
||||
}
|
||||
|
||||
TEST_CASE("JFJochReader_HDF5DataSource_Legacy", "[HDF5][Full]") {
|
||||
DiffractionExperiment x(DetJF(1));
|
||||
|
||||
x.FilePrefix("source_legacy_mapping").ImagesPerTrigger(5).ImagesPerFile(2).OverwriteExistingFiles(true);
|
||||
x.BitDepthImage(16).SetFileWriterFormat(FileWriterFormat::NXmxLegacy).PixelSigned(true);
|
||||
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
|
||||
|
||||
std::vector<int16_t> image(x.GetPixelsNum(), 31);
|
||||
|
||||
RegisterHDF5Filter();
|
||||
{
|
||||
StartMessage start_message;
|
||||
x.FillMessage(start_message);
|
||||
FileWriter writer(start_message);
|
||||
|
||||
for (int i = 0; i < x.GetImageNum(); i++) {
|
||||
image[5678] = static_cast<int16_t>(300 + i);
|
||||
|
||||
DataMessage message{};
|
||||
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.number = i;
|
||||
REQUIRE_NOTHROW(writer.WriteHDF5(message));
|
||||
}
|
||||
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = x.GetImageNum();
|
||||
writer.WriteHDF5(end_message);
|
||||
writer.Finalize();
|
||||
}
|
||||
|
||||
{
|
||||
JFJochHDF5Reader reader;
|
||||
REQUIRE_NOTHROW(reader.ReadFile("source_legacy_mapping_master.h5"));
|
||||
|
||||
auto source = reader.GetHDF5DataSource(1, 3);
|
||||
REQUIRE(source.size() == 2);
|
||||
|
||||
CHECK(source[0].filename == "source_legacy_mapping_data_000001.h5");
|
||||
CHECK(source[0].dataset == "/entry/data/data");
|
||||
CHECK(source[0].source_first_image == 1);
|
||||
CHECK(source[0].virtual_first_image == 0);
|
||||
CHECK(source[0].image_count == 1);
|
||||
|
||||
CHECK(source[1].filename == "source_legacy_mapping_data_000002.h5");
|
||||
CHECK(source[1].dataset == "/entry/data/data");
|
||||
CHECK(source[1].source_first_image == 0);
|
||||
CHECK(source[1].virtual_first_image == 1);
|
||||
CHECK(source[1].image_count == 2);
|
||||
}
|
||||
|
||||
remove("source_legacy_mapping_master.h5");
|
||||
remove("source_legacy_mapping_data_000001.h5");
|
||||
remove("source_legacy_mapping_data_000002.h5");
|
||||
remove("source_legacy_mapping_data_000003.h5");
|
||||
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
|
||||
}
|
||||
|
||||
TEST_CASE("JFJochReader_ProcessingHDF5_FromVDS_MapsToDataFiles", "[HDF5][Full]") {
|
||||
DiffractionExperiment x(DetJF(1));
|
||||
|
||||
x.FilePrefix("proc_source_vds").ImagesPerTrigger(5).ImagesPerFile(2).OverwriteExistingFiles(true);
|
||||
x.BitDepthImage(16).SetFileWriterFormat(FileWriterFormat::NXmxVDS).PixelSigned(true);
|
||||
x.Compression(CompressionAlgorithm::NO_COMPRESSION);
|
||||
|
||||
std::vector<int16_t> image(x.GetPixelsNum(), 51);
|
||||
|
||||
RegisterHDF5Filter();
|
||||
{
|
||||
StartMessage start_message;
|
||||
x.FillMessage(start_message);
|
||||
FileWriter writer(start_message);
|
||||
|
||||
for (int i = 0; i < x.GetImageNum(); i++) {
|
||||
image[5678] = static_cast<int16_t>(500 + i);
|
||||
|
||||
DataMessage message{};
|
||||
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.number = i;
|
||||
REQUIRE_NOTHROW(writer.WriteHDF5(message));
|
||||
}
|
||||
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = x.GetImageNum();
|
||||
writer.WriteHDF5(end_message);
|
||||
writer.Finalize();
|
||||
}
|
||||
|
||||
std::vector<HDF5DataSourceMessage> source_data;
|
||||
{
|
||||
JFJochHDF5Reader reader;
|
||||
REQUIRE_NOTHROW(reader.ReadFile("proc_source_vds_master.h5"));
|
||||
source_data = reader.GetHDF5DataSource(1, 3);
|
||||
|
||||
REQUIRE(source_data.size() == 2);
|
||||
CHECK(source_data[0].filename == "proc_source_vds_data_000001.h5");
|
||||
CHECK(source_data[1].filename == "proc_source_vds_data_000002.h5");
|
||||
}
|
||||
|
||||
{
|
||||
DiffractionExperiment proc_x = x;
|
||||
proc_x.FilePrefix("proc_from_vds")
|
||||
.ImagesPerTrigger(3)
|
||||
.SetFileWriterFormat(FileWriterFormat::NXmxIntegrated)
|
||||
.OverwriteExistingFiles(true);
|
||||
|
||||
StartMessage start_message;
|
||||
proc_x.FillMessage(start_message);
|
||||
start_message.number_of_images = 3;
|
||||
start_message.images_per_file = 3;
|
||||
start_message.write_images = false;
|
||||
start_message.write_master_file = true;
|
||||
start_message.hdf5_source_data = source_data;
|
||||
|
||||
FileWriter writer(start_message);
|
||||
|
||||
for (int i = 0; i < 3; i++) {
|
||||
DataMessage message{};
|
||||
message.number = i;
|
||||
message.original_number = i + 1;
|
||||
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.spot_count = 200 + i;
|
||||
REQUIRE_NOTHROW(writer.WriteHDF5(message));
|
||||
}
|
||||
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = 3;
|
||||
writer.WriteHDF5(end_message);
|
||||
writer.Finalize();
|
||||
}
|
||||
|
||||
{
|
||||
HDF5ReadOnlyFile file("proc_from_vds_master.h5");
|
||||
HDF5DataSet data(file, "/entry/data/data");
|
||||
HDF5Dcpl dcpl(data);
|
||||
|
||||
REQUIRE(dcpl.GetLayout() == HDF5DataSetLayout::VIRTUAL);
|
||||
|
||||
auto mappings = dcpl.GetVirtualMappings();
|
||||
REQUIRE(mappings.size() == 2);
|
||||
|
||||
CHECK(mappings[0].filename == "proc_source_vds_data_000001.h5");
|
||||
CHECK(mappings[1].filename == "proc_source_vds_data_000002.h5");
|
||||
}
|
||||
|
||||
{
|
||||
JFJochHDF5Reader reader;
|
||||
REQUIRE_NOTHROW(reader.ReadFile("proc_from_vds_master.h5"));
|
||||
|
||||
auto img0 = reader.LoadImage(0);
|
||||
REQUIRE(img0);
|
||||
CHECK(img0->Image()[5678] == 501);
|
||||
|
||||
auto img2 = reader.LoadImage(2);
|
||||
REQUIRE(img2);
|
||||
CHECK(img2->Image()[5678] == 503);
|
||||
}
|
||||
|
||||
remove("proc_source_vds_master.h5");
|
||||
remove("proc_source_vds_data_000001.h5");
|
||||
remove("proc_source_vds_data_000002.h5");
|
||||
remove("proc_source_vds_data_000003.h5");
|
||||
remove("proc_from_vds_master.h5");
|
||||
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
|
||||
}
|
||||
|
||||
TEST_CASE("JFJochReader_ReadReflections_VDS", "[HDF5][Full]") {
|
||||
DiffractionExperiment x(DetJF(1));
|
||||
|
||||
x.FilePrefix("read_reflections_vds")
|
||||
.ImagesPerTrigger(4)
|
||||
.ImagesPerFile(1)
|
||||
.OverwriteExistingFiles(true)
|
||||
.BitDepthImage(16)
|
||||
.PixelSigned(true)
|
||||
.SetFileWriterFormat(FileWriterFormat::NXmxVDS)
|
||||
.IndexingAlgorithm(IndexingAlgorithmEnum::FFT)
|
||||
.Compression(CompressionAlgorithm::NO_COMPRESSION);
|
||||
|
||||
std::vector<int16_t> image(x.GetPixelsNum(), 0);
|
||||
|
||||
RegisterHDF5Filter();
|
||||
|
||||
{
|
||||
StartMessage start_message;
|
||||
x.FillMessage(start_message);
|
||||
|
||||
FileWriter writer(start_message);
|
||||
ScanResultGenerator scan_result(x);
|
||||
|
||||
for (int i = 0; i < x.GetImageNum(); i++) {
|
||||
DataMessage message{};
|
||||
message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum());
|
||||
message.number = i;
|
||||
|
||||
if (i == 1 || i == 3) {
|
||||
message.integrated_reflections = 2;
|
||||
message.reflections = {
|
||||
Reflection{
|
||||
.h = static_cast<int32_t>(10 + i),
|
||||
.k = 20,
|
||||
.l = 30,
|
||||
.delta_phi_deg = 0.1f,
|
||||
.predicted_x = 100.0f + static_cast<float>(i),
|
||||
.predicted_y = 200.0f + static_cast<float>(i),
|
||||
.d = 1.5f,
|
||||
.I = 1000.0f + static_cast<float>(i),
|
||||
.bkg = 10.0f,
|
||||
.sigma = 2.0f,
|
||||
.rlp = 1.0f,
|
||||
.partiality = 0.9f,
|
||||
.zeta = 0.01f
|
||||
},
|
||||
Reflection{
|
||||
.h = static_cast<int32_t>(40 + i),
|
||||
.k = 50,
|
||||
.l = 60,
|
||||
.delta_phi_deg = 0.2f,
|
||||
.predicted_x = 300.0f + static_cast<float>(i),
|
||||
.predicted_y = 400.0f + static_cast<float>(i),
|
||||
.d = 2.5f,
|
||||
.I = 2000.0f + static_cast<float>(i),
|
||||
.bkg = 20.0f,
|
||||
.sigma = 3.0f,
|
||||
.rlp = 1.0f,
|
||||
.partiality = 0.8f,
|
||||
.zeta = 0.02f
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
REQUIRE_NOTHROW(writer.WriteHDF5(message));
|
||||
scan_result.Add(message);
|
||||
}
|
||||
|
||||
EndMessage end_message;
|
||||
end_message.max_image_number = x.GetImageNum();
|
||||
scan_result.FillEndMessage(end_message);
|
||||
|
||||
writer.WriteHDF5(end_message);
|
||||
writer.Finalize();
|
||||
}
|
||||
|
||||
{
|
||||
JFJochHDF5Reader reader;
|
||||
REQUIRE_NOTHROW(reader.ReadFile("read_reflections_vds_master.h5"));
|
||||
|
||||
auto reflections = reader.ReadReflections();
|
||||
|
||||
REQUIRE(reflections.size() == 4);
|
||||
|
||||
CHECK(reflections[0].empty());
|
||||
|
||||
REQUIRE(reflections[1].size() == 2);
|
||||
CHECK(reflections[1][0].h == 11);
|
||||
CHECK(reflections[1][0].k == 20);
|
||||
CHECK(reflections[1][0].l == 30);
|
||||
CHECK(reflections[1][0].I == Catch::Approx(1001.0f));
|
||||
CHECK(reflections[1][0].predicted_x == Catch::Approx(101.0f));
|
||||
CHECK(reflections[1][0].predicted_y == Catch::Approx(201.0f));
|
||||
|
||||
CHECK(reflections[2].empty());
|
||||
|
||||
REQUIRE(reflections[3].size() == 2);
|
||||
CHECK(reflections[3][0].h == 13);
|
||||
CHECK(reflections[3][0].I == Catch::Approx(1003.0f));
|
||||
CHECK(reflections[3][1].h == 43);
|
||||
CHECK(reflections[3][1].I == Catch::Approx(2003.0f));
|
||||
}
|
||||
|
||||
{
|
||||
JFJochHDF5Reader reader;
|
||||
REQUIRE_NOTHROW(reader.ReadFile("read_reflections_vds_master.h5"));
|
||||
|
||||
auto reflections = reader.ReadReflections(1, 4);
|
||||
|
||||
REQUIRE(reflections.size() == 3);
|
||||
|
||||
REQUIRE(reflections[0].size() == 2); // original image 1
|
||||
CHECK(reflections[0][0].h == 11);
|
||||
|
||||
CHECK(reflections[1].empty()); // original image 2
|
||||
|
||||
REQUIRE(reflections[2].size() == 2); // original image 3
|
||||
CHECK(reflections[2][0].h == 13);
|
||||
}
|
||||
|
||||
remove("read_reflections_vds_master.h5");
|
||||
remove("read_reflections_vds_data_000001.h5");
|
||||
remove("read_reflections_vds_data_000002.h5");
|
||||
remove("read_reflections_vds_data_000003.h5");
|
||||
remove("read_reflections_vds_data_000004.h5");
|
||||
|
||||
REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0);
|
||||
}
|
||||
@@ -2,6 +2,7 @@
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include <random>
|
||||
#include <future>
|
||||
#include <catch2/catch_all.hpp>
|
||||
|
||||
#include "../image_pusher/TCPStreamPusher.h"
|
||||
@@ -425,7 +426,10 @@ TEST_CASE("TCPImageCommTest_AutoPort_StarBind", "[TCP]") {
|
||||
TCPStreamPusher pusher("tcp://127.0.0.1:*", 1);
|
||||
TCPImagePuller puller(pusher.GetAddress()[0], 64 * 1024 * 1024);
|
||||
|
||||
std::thread receiver([&] {
|
||||
std::this_thread::sleep_for(std::chrono::seconds(2));
|
||||
REQUIRE(pusher.GetConnectedWriters() == 1);
|
||||
|
||||
std::future<void> receiver = std::async(std::launch::async, [&] {
|
||||
bool seen_end = false;
|
||||
uint64_t processed = 0;
|
||||
|
||||
@@ -474,7 +478,7 @@ TEST_CASE("TCPImageCommTest_AutoPort_StarBind", "[TCP]") {
|
||||
}
|
||||
|
||||
REQUIRE(pusher.EndDataCollection(end));
|
||||
receiver.join();
|
||||
REQUIRE_NOTHROW(receiver.get());
|
||||
puller.Disconnect();
|
||||
}
|
||||
|
||||
|
||||
@@ -300,13 +300,13 @@ TEST_CASE("XtalOptimizer_hexagonal") {
|
||||
xtal_opt.geom.BeamX_pxl(1007).BeamY_pxl(990).DetectorDistance_mm(200)
|
||||
.PoniRot1_rad(0.01).PoniRot2_rad(0.02);
|
||||
xtal_opt.crystal_system = gemmi::CrystalSystem::Hexagonal;
|
||||
xtal_opt.max_time = 30.0;
|
||||
xtal_opt.max_time = 60.0;
|
||||
auto start = std::chrono::high_resolution_clock::now();
|
||||
REQUIRE(XtalOptimizer(xtal_opt, spots));
|
||||
bool ret = XtalOptimizer(xtal_opt, spots);
|
||||
auto end = std::chrono::high_resolution_clock::now();
|
||||
std::cout << "XtalOptimizer took " << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
|
||||
<< " microseconds" << std::endl;
|
||||
|
||||
REQUIRE(ret);
|
||||
auto uc_i = latt_i.GetUnitCell();
|
||||
auto uc_o = xtal_opt.latt.GetUnitCell();
|
||||
|
||||
|
||||
+42
-31
@@ -37,21 +37,30 @@ void print_usage(Logger &logger) {
|
||||
logger.Info(" -s<num> Start image number (default: 0)");
|
||||
logger.Info(" -e<num> End image number (default: all)");
|
||||
logger.Info(" -v Verbose output");
|
||||
logger.Info(" -R[num] Rotation indexing (optional: min angular range deg)");
|
||||
logger.Info(" -F Use FFT indexing algorithm (shortcut for -XFFT)");
|
||||
logger.Info(" -X<txt> Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)");
|
||||
logger.Info(" -x No least-square beam center refinement");
|
||||
logger.Info(" -W Include images in the written HDF5 file (otherwise only analysis results are saved)");
|
||||
logger.Info(" -U Unmerged intensities are written to a text file");
|
||||
logger.Info("");
|
||||
|
||||
logger.Info(" Spot finding");
|
||||
logger.Info(" -T<num> Noise sigma level for spot finding (default: 3.0)");
|
||||
logger.Info(" -t<num> Photon count threshold for spot finding (default: 10)");
|
||||
logger.Info(" -d<num> High resolution limit for spot finding (default: 1.5)");
|
||||
logger.Info(" -c<num> Max spot count (default: 250)");
|
||||
logger.Info("");
|
||||
logger.Info(" Indexing");
|
||||
logger.Info(" -R[num] Rotation indexing (optional: min angular range deg)");
|
||||
logger.Info(" -X<txt> Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)");
|
||||
logger.Info(" -F Use FFT indexing algorithm (shortcut for -XFFT)");
|
||||
logger.Info(" -S<num> Space group number - used for both indexing and scaling");
|
||||
logger.Info(" -C<cell> Fix reference unit cell: -C\"a,b,c,alpha,beta,gamma\" (comma-separated, no spaces; quotes optional)");
|
||||
logger.Info(" -x No least-square beam center refinement");
|
||||
|
||||
logger.Info("");
|
||||
logger.Info(" Scaling and merging");
|
||||
logger.Info(" -D<num> High resolution limit for scaling/merging (default: 0.0; no limit)");
|
||||
logger.Info(" -S<num> Space group number");
|
||||
logger.Info(" -M Scale and merge (refine mosaicity) and write scaled.hkl + image.dat");
|
||||
logger.Info(" -P<txt> Partiality refinement fixed|rot|unity (default: fixed)");
|
||||
logger.Info(" -A Anomalous mode (don't merge Friedel pairs)");
|
||||
logger.Info(" -C<cell> Fix reference unit cell: -C\"a,b,c,alpha,beta,gamma\" (comma-separated, no spaces; quotes optional)");
|
||||
logger.Info(" -c<num> Max spot count (default: 250)");
|
||||
logger.Info(" -W HDF5 file with analysis results is written");
|
||||
logger.Info(" -T<num> Noise sigma level for spot finding (default: 3.0)");
|
||||
logger.Info(" -t<num> Photon count threshold for spot finding (default: 10)");
|
||||
}
|
||||
|
||||
void trim_in_place(std::string& t) {
|
||||
@@ -364,15 +373,9 @@ int main(int argc, char **argv) {
|
||||
|
||||
start_message.pixel_mask["default"] = pixel_mask.GetMask(experiment);
|
||||
start_message.max_spot_count = experiment.GetMaxSpotCount();
|
||||
|
||||
std::unique_ptr<FileWriter> writer;
|
||||
try {
|
||||
if (write_output)
|
||||
writer = std::make_unique<FileWriter>(start_message);
|
||||
} catch (const std::exception &e) {
|
||||
logger.Error("Failed to initialize file writer: {}", e.what());
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
start_message.write_images = write_output;
|
||||
start_message.file_format = FileWriterFormat::NXmxIntegrated;
|
||||
start_message.master_suffix = "process";
|
||||
|
||||
// 4. Processing Setup
|
||||
int total_images_in_file = reader.GetNumberOfImages();
|
||||
@@ -392,7 +395,22 @@ int main(int argc, char **argv) {
|
||||
|
||||
std::atomic<int> processed_count = 0;
|
||||
std::atomic<uint64_t> total_uncompressed_bytes = 0;
|
||||
std::atomic<uint64_t> max_image_number_sent = 0;
|
||||
|
||||
start_message.file_format = FileWriterFormat::NXmxIntegrated;
|
||||
start_message.write_master_file = true;
|
||||
start_message.write_images = false;
|
||||
start_message.number_of_images = images_to_process;
|
||||
start_message.images_per_file = images_to_process;
|
||||
start_message.hdf5_source_data = reader.GetHDF5DataSource(start_image, images_to_process);
|
||||
|
||||
std::unique_ptr<FileWriter> writer;
|
||||
try {
|
||||
if (!output_prefix.empty())
|
||||
writer = std::make_unique<FileWriter>(start_message);
|
||||
} catch (const std::exception &e) {
|
||||
logger.Error("Failed to initialize file writer: {}", e.what());
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
// Mimic JFJochReceiver lattice handling (IndexAndRefine handles the logic per thread,
|
||||
// but we need a central accumulator or use the pool's functionality if IndexAndRefine wraps it)
|
||||
@@ -430,7 +448,8 @@ int main(int argc, char **argv) {
|
||||
DataMessage msg{};
|
||||
|
||||
msg.image = img->image;
|
||||
msg.number = image_idx;
|
||||
msg.number = current_idx_offset;
|
||||
msg.original_number = image_idx;
|
||||
msg.image_collection_efficiency = dataset->efficiency[image_idx];
|
||||
|
||||
total_uncompressed_bytes += msg.image.GetUncompressedSize();
|
||||
@@ -458,13 +477,6 @@ int main(int argc, char **argv) {
|
||||
if (writer)
|
||||
writer->Write(msg);
|
||||
|
||||
// Update max sent tracking
|
||||
uint64_t current_max = max_image_number_sent.load();
|
||||
while (static_cast<uint64_t>(msg.number) > current_max) {
|
||||
if (max_image_number_sent.compare_exchange_weak(current_max, static_cast<uint64_t>(msg.number)))
|
||||
break;
|
||||
}
|
||||
|
||||
finished_count.fetch_add(1);
|
||||
|
||||
// Progress log
|
||||
@@ -510,7 +522,7 @@ int main(int argc, char **argv) {
|
||||
|
||||
// 5. Finalize Statistics and Write EndMessage
|
||||
EndMessage end_msg;
|
||||
end_msg.max_image_number = max_image_number_sent;
|
||||
end_msg.max_image_number = images_to_process;
|
||||
end_msg.images_collected_count = images_to_process;
|
||||
end_msg.images_sent_to_write_count = images_to_process;
|
||||
end_msg.end_date = time_UTC(std::chrono::system_clock::now());
|
||||
@@ -536,7 +548,6 @@ int main(int argc, char **argv) {
|
||||
logger.Info("Rotation Indexing found lattice");
|
||||
}
|
||||
|
||||
// --- Optional: run scaling (mosaicity refinement) on accumulated reflections ---
|
||||
// --- Optional: run scaling (mosaicity refinement) on accumulated reflections ---
|
||||
if (run_scaling) {
|
||||
logger.Info("Running scaling (mosaicity refinement) ...");
|
||||
@@ -606,7 +617,7 @@ int main(int argc, char **argv) {
|
||||
}
|
||||
|
||||
if (scale_result) {
|
||||
end_msg.scale_factor = scale_result->image_scale_g;
|
||||
end_msg.image_scale_factor = scale_result->image_scale_g;
|
||||
|
||||
logger.Info("Scaling completed in {:.2f} s ({} unique reflections)",
|
||||
scale_time, scale_result->merged.size());
|
||||
|
||||
@@ -109,7 +109,7 @@ void JFJochViewerMenu::openSelected() {
|
||||
this,
|
||||
"Open File", // Dialog title
|
||||
"", // Default folder
|
||||
"HDF5 Master Files (*_master.h5);; HDF5 Files (*.h5);;All Files (*)" // Filter for .h5 files
|
||||
"HDF5 Master Files (*_master.h5 *_process.h5);; HDF5 Files (*.h5);;All Files (*)" // Filter for .h5 files
|
||||
);
|
||||
|
||||
if (!fileName.isEmpty())
|
||||
|
||||
+36
-19
@@ -71,26 +71,34 @@ void FileWriter::WriteHDF5(const DataMessage& msg) {
|
||||
if (msg.number < 0)
|
||||
throw JFJochException(JFJochExceptionCategory::ArrayOutOfBounds, "No support for negative images");
|
||||
|
||||
const uint64_t file_number = (start_message.images_per_file == 0) ? 0 : msg.number / start_message.images_per_file;
|
||||
const uint64_t image_number = (start_message.images_per_file == 0) ? msg.number : msg.number % start_message.images_per_file;
|
||||
|
||||
if (closed_files.contains(file_number))
|
||||
return;
|
||||
|
||||
if (files.size() <= file_number)
|
||||
files.resize(file_number + 1);
|
||||
|
||||
if (!files[file_number]) {
|
||||
files[file_number] = std::make_unique<HDF5DataFile>(start_message, file_number);
|
||||
if (format == FileWriterFormat::NXmxIntegrated && master_file)
|
||||
files[file_number]->CreateFile(msg, master_file->GetFile());
|
||||
}
|
||||
files[file_number]->Write(msg, image_number);
|
||||
|
||||
if (files[file_number]->GetNumImages() == start_message.images_per_file) {
|
||||
CloseFile(file_number);
|
||||
if (format == FileWriterFormat::NXmxIntegrated && master_file) {
|
||||
if (files.empty() )
|
||||
files.resize(1);
|
||||
if (!files[0]) {
|
||||
files[0] = std::make_unique<HDF5DataFile>(start_message, 0, HDF5Metadata::MasterFileName(start_message));
|
||||
files[0]->CreateFile(msg, master_file->GetFile());
|
||||
}
|
||||
files[0]->Write(msg, msg.number);
|
||||
} else {
|
||||
CloseOldFiles(static_cast<uint64_t>(msg.number));
|
||||
const uint64_t file_number = (start_message.images_per_file == 0) ? 0 : msg.number / start_message.images_per_file;
|
||||
const uint64_t image_number = (start_message.images_per_file == 0) ? msg.number : msg.number % start_message.images_per_file;
|
||||
|
||||
if (closed_files.contains(file_number))
|
||||
return;
|
||||
|
||||
if (files.size() <= file_number)
|
||||
files.resize(file_number + 1);
|
||||
|
||||
if (!files[file_number])
|
||||
files[file_number] = std::make_unique<HDF5DataFile>(start_message, file_number,
|
||||
HDF5Metadata::DataFileName(start_message, file_number));
|
||||
files[file_number]->Write(msg, image_number);
|
||||
|
||||
if (files[file_number]->GetNumImages() == start_message.images_per_file) {
|
||||
CloseFile(file_number);
|
||||
} else {
|
||||
CloseOldFiles(static_cast<uint64_t>(msg.number));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -230,6 +238,15 @@ void FileWriter::WriteHDF5(const CompressedImage &msg) {
|
||||
void FileWriter::WriteHDF5(const EndMessage &msg) {
|
||||
if (master_file) {
|
||||
std::lock_guard<std::mutex> lock(hdf5_mutex);
|
||||
|
||||
if (format == FileWriterFormat::NXmxIntegrated) {
|
||||
try {
|
||||
CloseFile(0);
|
||||
} catch (...) {
|
||||
throw;
|
||||
}
|
||||
}
|
||||
|
||||
master_file->Finalize(msg);
|
||||
}
|
||||
}
|
||||
|
||||
+48
-40
@@ -20,9 +20,10 @@
|
||||
#include "HDF5NXmx.h"
|
||||
#include "../common/time_utc.h"
|
||||
|
||||
HDF5DataFile::HDF5DataFile(const StartMessage &msg, uint64_t in_file_number) {
|
||||
file_number = in_file_number;
|
||||
|
||||
HDF5DataFile::HDF5DataFile(const StartMessage &msg, uint64_t file_number, const std::string &filename) :
|
||||
filename(filename),
|
||||
file_number(file_number),
|
||||
write_images(msg.write_images.value_or(true)) {
|
||||
if (msg.overwrite.has_value())
|
||||
overwrite = msg.overwrite.value();
|
||||
|
||||
@@ -30,9 +31,14 @@ HDF5DataFile::HDF5DataFile(const StartMessage &msg, uint64_t in_file_number) {
|
||||
ypixel = 0;
|
||||
max_image_number = 0;
|
||||
nimages = 0;
|
||||
filename = HDF5Metadata::DataFileName(msg, file_number);
|
||||
image_low = file_number * msg.images_per_file;
|
||||
images_per_file = msg.images_per_file;
|
||||
|
||||
if (msg.file_format == FileWriterFormat::NXmxIntegrated) {
|
||||
image_low = 0;
|
||||
images_per_file = msg.number_of_images;
|
||||
} else {
|
||||
image_low = file_number * msg.images_per_file;
|
||||
images_per_file = msg.images_per_file;
|
||||
}
|
||||
|
||||
timestamp.reserve(images_per_file);
|
||||
exptime.reserve(images_per_file);
|
||||
@@ -110,7 +116,6 @@ HDF5DataFile::~HDF5DataFile() {
|
||||
if (data_file) {
|
||||
try {
|
||||
data_set.reset();
|
||||
data_set_image_number.reset();
|
||||
data_file.reset();
|
||||
if (manage_file) {
|
||||
std::error_code ec;
|
||||
@@ -124,42 +129,44 @@ HDF5DataFile::~HDF5DataFile() {
|
||||
}
|
||||
}
|
||||
|
||||
void HDF5DataFile::CreateFile(const DataMessage& msg, std::shared_ptr<HDF5File> in_data_file, bool integrated) {
|
||||
HDF5Dcpl dcpl;
|
||||
|
||||
HDF5DataType data_type(msg.image.GetMode());
|
||||
|
||||
xpixel = msg.image.GetWidth();
|
||||
ypixel = msg.image.GetHeight();
|
||||
|
||||
dcpl.SetCompression(msg.image.GetCompressionAlgorithm(), JFJochBitShuffleCompressor::DefaultBlockSize);
|
||||
dcpl.SetChunking( {1, ypixel, xpixel});
|
||||
|
||||
H5Pset_fill_time(dcpl.GetID(), H5D_FILL_TIME_NEVER);
|
||||
H5Pset_alloc_time(dcpl.GetID(), H5D_ALLOC_TIME_INCR);
|
||||
|
||||
switch (msg.image.GetMode()) {
|
||||
case CompressedImageMode::Int8:
|
||||
dcpl.SetFillValue8(INT8_MIN);
|
||||
break;
|
||||
case CompressedImageMode::Int16:
|
||||
dcpl.SetFillValue16(INT16_MIN);
|
||||
break;
|
||||
case CompressedImageMode::Int32:
|
||||
dcpl.SetFillValue32(INT32_MIN);
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
|
||||
void HDF5DataFile::CreateFile(const DataMessage& msg, std::shared_ptr<HDF5File> in_data_file) {
|
||||
data_file = in_data_file;
|
||||
|
||||
HDF5Group(*data_file, "/entry").NXClass("NXentry");
|
||||
HDF5Group(*data_file, "/entry/data").NXClass("NXdata");
|
||||
if (write_images) {
|
||||
HDF5Dcpl dcpl;
|
||||
|
||||
HDF5DataType data_type(msg.image.GetMode());
|
||||
|
||||
xpixel = msg.image.GetWidth();
|
||||
ypixel = msg.image.GetHeight();
|
||||
|
||||
dcpl.SetCompression(msg.image.GetCompressionAlgorithm(), JFJochBitShuffleCompressor::DefaultBlockSize);
|
||||
dcpl.SetChunking( {1, ypixel, xpixel});
|
||||
|
||||
H5Pset_fill_time(dcpl.GetID(), H5D_FILL_TIME_NEVER);
|
||||
H5Pset_alloc_time(dcpl.GetID(), H5D_ALLOC_TIME_INCR);
|
||||
|
||||
switch (msg.image.GetMode()) {
|
||||
case CompressedImageMode::Int8:
|
||||
dcpl.SetFillValue8(INT8_MIN);
|
||||
break;
|
||||
case CompressedImageMode::Int16:
|
||||
dcpl.SetFillValue16(INT16_MIN);
|
||||
break;
|
||||
case CompressedImageMode::Int32:
|
||||
dcpl.SetFillValue32(INT32_MIN);
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
|
||||
HDF5Group(*data_file, "/entry/data").NXClass("NXdata");
|
||||
HDF5DataSpace data_space({1, ypixel, xpixel}, {H5S_UNLIMITED, ypixel, xpixel});
|
||||
data_set = std::make_unique<HDF5DataSet>(*data_file, "/entry/data/data", data_type, data_space, dcpl);
|
||||
data_set->SetExtent({images_per_file, ypixel, xpixel});
|
||||
}
|
||||
|
||||
HDF5DataSpace data_space({1, ypixel, xpixel}, {H5S_UNLIMITED, ypixel, xpixel});
|
||||
data_set = std::make_unique<HDF5DataSet>(*data_file, "/entry/data/data", data_type, data_space, dcpl);
|
||||
data_set->SetExtent({images_per_file, ypixel, xpixel});
|
||||
for (auto &p: plugins)
|
||||
p->OpenFile(*data_file, msg, images_per_file);
|
||||
}
|
||||
@@ -186,7 +193,8 @@ void HDF5DataFile::Write(const DataMessage &msg, uint64_t image_number) {
|
||||
}
|
||||
|
||||
nimages++;
|
||||
data_set->WriteDirectChunk(msg.image.GetCompressed(), msg.image.GetCompressedSize(), {image_number, 0, 0});
|
||||
if (data_set)
|
||||
data_set->WriteDirectChunk(msg.image.GetCompressed(), msg.image.GetCompressedSize(), {image_number, 0, 0});
|
||||
|
||||
for (auto &p: plugins)
|
||||
p->Write(msg, image_number);
|
||||
|
||||
@@ -22,12 +22,14 @@ struct HDF5DataFileStatistics {
|
||||
};
|
||||
|
||||
class HDF5DataFile {
|
||||
std::string filename;
|
||||
const std::string filename;
|
||||
const uint64_t file_number;
|
||||
const bool write_images;
|
||||
|
||||
std::string tmp_filename;
|
||||
|
||||
std::shared_ptr<HDF5File> data_file = nullptr;
|
||||
std::unique_ptr<HDF5DataSet> data_set = nullptr;
|
||||
std::unique_ptr<HDF5DataSet> data_set_image_number = nullptr;
|
||||
std::vector<std::unique_ptr<HDF5DataFilePlugin>> plugins;
|
||||
size_t images_per_file;
|
||||
size_t xpixel;
|
||||
@@ -46,17 +48,16 @@ class HDF5DataFile {
|
||||
bool closed = false;
|
||||
|
||||
bool overwrite = false;
|
||||
int64_t file_number;
|
||||
bool new_file = true;
|
||||
bool manage_file = false;
|
||||
public:
|
||||
HDF5DataFile(const StartMessage &msg, uint64_t file_number);
|
||||
HDF5DataFile(const StartMessage &msg, uint64_t file_number, const std::string &filename);
|
||||
~HDF5DataFile();
|
||||
std::optional<HDF5DataFileStatistics> Close();
|
||||
void Write(const DataMessage& msg, uint64_t image_number);
|
||||
size_t GetNumImages() const;
|
||||
|
||||
void CreateFile(const DataMessage& msg, std::shared_ptr<HDF5File> data_file, bool integrated = false);
|
||||
void CreateFile(const DataMessage& msg, std::shared_ptr<HDF5File> data_file);
|
||||
};
|
||||
|
||||
#endif //HDF5DATAFILE_H
|
||||
|
||||
@@ -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());
|
||||
}
|
||||
|
||||
@@ -16,8 +16,6 @@ class HDF5DataFilePluginDetector : public HDF5DataFilePlugin {
|
||||
AutoIncrVector<float> efficiency;
|
||||
AutoIncrVector<uint64_t> packets_received;
|
||||
AutoIncrVector<uint64_t> packets_expected;
|
||||
AutoIncrVector<int64_t> pixel_sum;
|
||||
AutoIncrVector<float> processing_time{NAN};
|
||||
public:
|
||||
HDF5DataFilePluginDetector(const StartMessage& msg);
|
||||
void OpenFile(HDF5File &data_file, const DataMessage& msg, size_t images_per_file) override;
|
||||
|
||||
@@ -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());
|
||||
}
|
||||
|
||||
@@ -13,6 +13,7 @@ class HDF5DataFilePluginImageStats : public HDF5DataFilePlugin {
|
||||
AutoIncrVector<int64_t> min_value;
|
||||
AutoIncrVector<int64_t> error_pixels;
|
||||
AutoIncrVector<int64_t> saturated_pixels;
|
||||
AutoIncrVector<int64_t> 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;
|
||||
|
||||
@@ -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())
|
||||
|
||||
+151
-80
@@ -11,19 +11,19 @@
|
||||
#include "../common/time_utc.h"
|
||||
#include "gemmi/symmetry.hpp"
|
||||
|
||||
namespace {
|
||||
std::string GenFilename(const StartMessage &start) {
|
||||
return fmt::format("{:s}_master.h5", start.file_prefix);
|
||||
}
|
||||
std::string HDF5Metadata::MasterFileName(const StartMessage &start) {
|
||||
if (start.master_suffix.has_value())
|
||||
return fmt::format("{:s}_{:s}.h5", start.file_prefix, start.master_suffix.value());
|
||||
return fmt::format("{:s}_master.h5", start.file_prefix);
|
||||
}
|
||||
|
||||
NXmx::NXmx(const StartMessage &start)
|
||||
: start_message(start),
|
||||
filename(GenFilename(start)) {
|
||||
filename(HDF5Metadata::MasterFileName(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();
|
||||
}
|
||||
@@ -35,7 +35,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<HDF5File>(tmp_filename, v1_10);
|
||||
hdf5_file->Attr("file_name", filename);
|
||||
@@ -121,54 +122,12 @@ 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));
|
||||
VDS(start, "/entry/MX/nPeaks", {total_images}, HDF5DataType((uint32_t) 0));
|
||||
}
|
||||
|
||||
if (start.indexing_algorithm != IndexingAlgorithmEnum::None) {
|
||||
@@ -177,23 +136,9 @@ 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) {
|
||||
VDS(start, "/entry/detector/beam_center_x",
|
||||
"/entry/instrument/detector/refined_beam_center_x",
|
||||
{total_images},
|
||||
HDF5DataType(0.0f));
|
||||
VDS(start, "/entry/detector/beam_center_y",
|
||||
"/entry/instrument/detector/refined_beam_center_y",
|
||||
{total_images},
|
||||
HDF5DataType(0.0f));
|
||||
}
|
||||
if (!start.az_int_bin_to_q.empty()) {
|
||||
size_t azimuthal_bins = start.az_int_phi_bin_count.value_or(1);
|
||||
size_t q_bins = start.az_int_q_bin_count.value_or(1);
|
||||
@@ -217,20 +162,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<hsize_t>(mapping.virtual_first_image), 0, 0},
|
||||
{static_cast<hsize_t>(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<hsize_t>(mapping.source_first_image), 0, 0},
|
||||
{static_cast<hsize_t>(mapping.image_count), height, width}
|
||||
);
|
||||
|
||||
dcpl.SetVirtual(mapping.filename,
|
||||
source_dataset,
|
||||
source_data_space,
|
||||
virtual_data_space);
|
||||
}
|
||||
|
||||
auto data_dataset = std::make_unique<HDF5DataSet>(
|
||||
*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<int32_t>(1))
|
||||
.Attr("image_nr_high", static_cast<int32_t>(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<uint64_t>(start.images_per_file);
|
||||
const uint64_t image_in_file = image % static_cast<uint64_t>(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);
|
||||
}
|
||||
}
|
||||
|
||||
@@ -698,6 +712,19 @@ void NXmx::ADUHistogram(const EndMessage &end) {
|
||||
}
|
||||
}
|
||||
|
||||
template<class T>
|
||||
void SaveVectorIfMissing(HDF5Object &object,
|
||||
const std::string &path,
|
||||
const std::vector<T> &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)
|
||||
@@ -717,6 +744,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:
|
||||
@@ -726,6 +754,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;
|
||||
}
|
||||
@@ -744,9 +775,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) {
|
||||
@@ -786,3 +814,46 @@ void NXmx::UserData(const StartMessage &start) {
|
||||
std::shared_ptr<HDF5File> 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/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);
|
||||
}
|
||||
|
||||
+7
-4
@@ -1,14 +1,14 @@
|
||||
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#ifndef JUNGFRAUJOCH_HDF5NXMX_H
|
||||
#define JUNGFRAUJOCH_HDF5NXMX_H
|
||||
#pragma once
|
||||
|
||||
#include "../common/JFJochMessages.h"
|
||||
|
||||
#include "HDF5Objects.h"
|
||||
|
||||
namespace HDF5Metadata {
|
||||
std::string MasterFileName(const StartMessage &msg);
|
||||
std::string DataFileName(const StartMessage &msg, int64_t file_number);
|
||||
}
|
||||
|
||||
@@ -22,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<HDF5DataSet> VDS(const StartMessage &start,
|
||||
const std::string& name,
|
||||
const std::vector<hsize_t> &dim,
|
||||
@@ -48,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();
|
||||
@@ -61,5 +66,3 @@ public:
|
||||
|
||||
std::shared_ptr<HDF5File> GetFile();
|
||||
};
|
||||
|
||||
#endif //JUNGFRAUJOCH_HDF5NXMX_H
|
||||
|
||||
Reference in New Issue
Block a user