From 4c45d1c835b03bb517c49617d06c3d54eb2090c6 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Fri, 5 Jun 2026 16:22:35 +0200 Subject: [PATCH] Multi lattice search: WIP - Integration added by Claude Opus (doesn't work :) ) --- common/JFJochMessages.h | 2 +- docs/CBOR.md | 207 ++++++++++---------- frame_serialize/CBORStream2Deserializer.cpp | 26 ++- frame_serialize/CBORStream2Serializer.cpp | 24 ++- image_analysis/IndexAndRefine.cpp | 46 ++++- image_analysis/IndexAndRefine.h | 1 + image_analysis/indexing/AnalyzeIndexing.cpp | 59 +++++- image_analysis/indexing/AnalyzeIndexing.h | 3 +- receiver/JFJochReceiver.cpp | 1 + writer/HDF5DataFilePluginMX.cpp | 32 ++- writer/HDF5DataFilePluginMX.h | 2 + 11 files changed, 286 insertions(+), 117 deletions(-) diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index 94d3fa83..ccc45015 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -223,7 +223,7 @@ struct StartMessage { std::optional unit_cell; // user data std::optional space_group_number; // user data uint64_t max_spot_count; // user data - uint64_t max_extra_lattices; + uint64_t max_extra_lattices = 0; std::optional storage_cell_number; uint64_t storage_cell_delay_ns; diff --git a/docs/CBOR.md b/docs/CBOR.md index 58fb8bb9..2c4c8869 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -64,6 +64,7 @@ There are minor differences at the moment: | 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 | | +| max_extra_lattices | uint64 | Maximum number of extra lattices | | | 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\] | | @@ -118,108 +119,109 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s ## Image message -| Field name | Type | Description | Present in DECTRIS format | Optional | -|-----------------------------|----------------------|-----------------------------------------------------------------------------------------------------------------------------------|:-------------------------:|:--------:| -| type | String | value "image" | 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 | | -| image_id | uint64 | Number of image within the series; for MX lossy compression this is sequential excluding removed frames | X | | -| original_image_id | uint64 | Number of image within the series; for MX lossy compression this includes removed frames in the count | | | -| real_time | Rational | Exposure time | X | | -| start_time | Rational | Exposure start time (highly approximate) | X | | -| end_time | Rational | Exposure end time (highly approximate) | X | | -| spots | Array(object) | Spots: | | | -| - x | float | observed position in x (pixels) | | | -| - y | float | observed position in y (pixels) | | | -| - I | float | intensity (photons) | | | -| - maxc | int64 | max count (photons) | | | -| - ice_ring | bool | spot in resolution range for ice rings | | | -| - indexed | bool | indexed solution | | | -| - latt | int64 | Lattice to which the peak belongs (negative number = not indexed) | | | -| reflections | Array(object) | Reflections: | | | -| - h | int64 | Miller index | | | -| - k | int64 | Miller index | | | -| - l | int64 | Miller index | | | -| - x | float | prediced position in x (pixels) | | | -| - y | float | predicted position in y (pixels) | | | -| - obs_x | float | observed position in x (pixels) | | | -| - obs_y | float | observed position in y (pixels) | | | -| - d | float | resolution \[Angstrom\] | | | -| - I | float | integrated intensity (photons) | | | -| - bkg | float | mean background value (photons) | | | -| - sigma | float | standard deviation, estimated from counting statistics (photons) | | | -| - image | float | image number (present for each spot) | | | -| - dist_ewald | float | distance to Ewald sphere (present only for indexed spots) | | | -| - rlp | float | Reciprocal Lorentz and polarization corrections | | | -| - partiality | float | Partiality of the reflection | | | -| spot_count | uint64 | Spot count | | | -| spot_count_ice_rings | uint64 | Number of spots within identified rings (experimental) | | | -| spot_count_low_res | uint64 | Number of spots in low resolution (prior to filtering) | | | -| spot_count_indexed | uint64 | Number of spots which fit indexing solution within a given tolerance | | | -| az_int_profile | Array(float) | Azimuthal integration results, use az_int_bin_to_q from start message for legend | | | -| | | NaN is used for empty bins and has to be taken care by the receiver | | | -| az_int_std | Array(float) | Standard deviation for azimuthal integration. (NaN for less than 2 samples) | | | -| az_int_count | Array(uint64) | Number of pixels contributing to azimuthal bin | | | -| indexing_result | bool | Indexing successful | | | -| indexing_lattice_count | uint64 | Number of lattices found in indexing | | X | -| indexing_lattice | Array(9 * float) | Indexing result real lattice; present only if indexed | | X | -| indexing_unit_cell | object | Indexing result unit cell: a, b, c \[angstrom\] and alpha, beta, gamma \[degree\]; present only if indexed | | X | -| | | Unit cell is redundant to lattice - yet to simplify downstream programs to analyze results, both are provided | | | -| profile_radius | float | Profile radius of the image - describes distance of observed reflections from the Ewald sphere \[Angstrom^-1\] | | | -| integrated_reflections | int64 | Count of integrated reflections | | | -| mosaicity | float | Angular range of spots in image from a rotation scan \[degree\] | | | -| b_factor | float | Estimated B-factor (Angstrom^2) | | | -| compression_time | float | Time spent on compression/decompressing image \[s\] | | | -| preprocessing_time | float | Time spent on preparing the image for analysis \[s\] | | | -| azint_time | float | Time spent on azimuthal integration \[s\] | | | -| spot_finding_time | float | Time spent on spot finding \[s\] | | | -| indexing_time | float | Time spent on indexing \[s\] | | | -| refinement_time | float | Time spent on refinement of indexing solution and experimental geometry \[s\] | | | -| index_analysis_time | float | Time spent on analyzing idnexing solution, calculating profile radius and mosaicity \[s\] | | | -| bragg_prediction_time | float | Time spent on predicting Bragg spots \[s\] | | | -| integration_time | float | Time spent on Bragg integration \[s\] | | | -| image_scale_time | float | Time spent on on-the-fly scaling \[s\] | | | -| processing_time | float | Total processing time \[s\] | | | -| xfel_pulse_id | uint64 | Bunch ID (for pulsed source, e.g., SwissFEL) | | X | -| xfel_event_code | uint64 | Event code (for pulsed source, e.g., SwissFEL) | | X | -| lattice_type | object | Bravais lattice classification of the indexing result (present only if available) | | X | -| - 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 | | | -| jf_info | uint64 | Detector info field | | | -| receiver_aq_dev_delay | uint64 | Receiver internal delay | | | -| receiver_free_send_buf | uint64 | Receiver internal number of available buffer locations | | | -| receiver_buf_in_sending | uint64 | Receiver internal number of buffer locations currently in sending/writing | | | -| receiver_buf_in_preparation | uint64 | Receiver internal number of buffer locations currently in processing | | | -| storage_cell | uint64 | Storage cell number | | | -| saturated_pixel_count | uint64 | Saturated pixel count | | | -| pixel_sum | uint64 | Sum of all pixels, excl. error and saturation | | | -| error_pixel_count | uint64 | Error pixel count | | | -| strong_pixel_count | uint64 | Strong pixel count (first stage of spot finding) | | | -| min_viable_pixel_value | int64 | Minimal pixel value, excl. error and saturation | | | -| max_viable_pixel_value | int64 | Maximal pixel value, excl. error and saturation | | | -| resolution_estimate | float | Diffraction resolution estimation \[Angstrom\] | | X | -| data_collection_efficiency | float | Image collection efficiency \[\] | | | -| packets_expected | uint64 | Number of packets expected per image (in units of 2 kB) | | | -| packets_received | uint64 | Number of packets received per image (in units of 2 kB) | | | -| bkg_estimate | float | Mean value for pixels in resolution range from 3.0 to 5.0 A \[photons\] | | | -| beam_corr_x | float | Beam center correction X applied during processing \[pixel\] | | X | -| beam_corr_y | float | Beam center correction Y applied during processing \[pixel\] | | X | -| image_scale_factor | float | Scaling result: Image scale factor (g) | | X | -| image_scale_mosaicity | float | Scaling result: Image scale mosaicity \[deg\] | | X | -| image_scale_b_factor | float | Scaling result: Image scale B factor \[Angstrom^2\] | | X | -| image_scale_cc | float | Scaling result: Image scale CC | | X | -| adu_histogram | Array(uint64) | ADU histogram | | | -| roi_integrals | object | Results of ROI calculation | | X | -| - sum | int64 | Sum of pixels in ROI area \[photons\] | | | -| - sum_square | int64 | Sum of squares of pixels in ROI area \[photons\] | | | -| - pixels | uint64 | Valid pixels in ROI area | | | -| - max_count | int64 | Highest count in ROI area \[photons\] | | | -| - x_weighted_sum | int64 | ROI pixel X position multiplied by photon count \[photons * pixels\] | | | -| - y_weighted_sum | int64 | ROI pixel Y position multiplied by photon count \[photons * pixels\] | | | -| user_data | string | Optional user defined text information - this is image_appendix serialized to JSON format | X | | -| data | Map(string -> Image) | Image | X | | +| Field name | Type | Description | Present in DECTRIS format | Optional | +|-----------------------------|-----------------------|-----------------------------------------------------------------------------------------------------------------------------------|:-------------------------:|:--------:| +| type | String | value "image" | 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 | | +| image_id | uint64 | Number of image within the series; for MX lossy compression this is sequential excluding removed frames | X | | +| original_image_id | uint64 | Number of image within the series; for MX lossy compression this includes removed frames in the count | | | +| real_time | Rational | Exposure time | X | | +| start_time | Rational | Exposure start time (highly approximate) | X | | +| end_time | Rational | Exposure end time (highly approximate) | X | | +| spots | Array(object) | Spots: | | | +| - x | float | observed position in x (pixels) | | | +| - y | float | observed position in y (pixels) | | | +| - I | float | intensity (photons) | | | +| - maxc | int64 | max count (photons) | | | +| - ice_ring | bool | spot in resolution range for ice rings | | | +| - indexed | bool | indexed solution | | | +| - latt | int64 | Lattice to which the peak belongs (negative number = not indexed) | | | +| reflections | Array(object) | Reflections: | | | +| - h | int64 | Miller index | | | +| - k | int64 | Miller index | | | +| - l | int64 | Miller index | | | +| - x | float | prediced position in x (pixels) | | | +| - y | float | predicted position in y (pixels) | | | +| - obs_x | float | observed position in x (pixels) | | | +| - obs_y | float | observed position in y (pixels) | | | +| - d | float | resolution \[Angstrom\] | | | +| - I | float | integrated intensity (photons) | | | +| - bkg | float | mean background value (photons) | | | +| - sigma | float | standard deviation, estimated from counting statistics (photons) | | | +| - image | float | image number (present for each spot) | | | +| - dist_ewald | float | distance to Ewald sphere (present only for indexed spots) | | | +| - rlp | float | Reciprocal Lorentz and polarization corrections | | | +| - partiality | float | Partiality of the reflection | | | +| spot_count | uint64 | Spot count | | | +| spot_count_ice_rings | uint64 | Number of spots within identified rings (experimental) | | | +| spot_count_low_res | uint64 | Number of spots in low resolution (prior to filtering) | | | +| spot_count_indexed | uint64 | Number of spots which fit indexing solution within a given tolerance | | | +| az_int_profile | Array(float) | Azimuthal integration results, use az_int_bin_to_q from start message for legend | | | +| | | NaN is used for empty bins and has to be taken care by the receiver | | | +| az_int_std | Array(float) | Standard deviation for azimuthal integration. (NaN for less than 2 samples) | | | +| az_int_count | Array(uint64) | Number of pixels contributing to azimuthal bin | | | +| indexing_result | bool | Indexing successful | | | +| indexing_lattice_count | uint64 | Number of lattices found in indexing | | X | +| indexing_lattice | Array(9 * float) | Indexing result real lattice; present only if indexed | | X | +| extra_lattices | Array(Array(9*float)) | Additional indexed lattices (orientation variants); present only if found | | | +| indexing_unit_cell | object | Indexing result unit cell: a, b, c \[angstrom\] and alpha, beta, gamma \[degree\]; present only if indexed | | X | +| | | Unit cell is redundant to lattice - yet to simplify downstream programs to analyze results, both are provided | | | +| profile_radius | float | Profile radius of the image - describes distance of observed reflections from the Ewald sphere \[Angstrom^-1\] | | | +| integrated_reflections | int64 | Count of integrated reflections | | | +| mosaicity | float | Angular range of spots in image from a rotation scan \[degree\] | | | +| b_factor | float | Estimated B-factor (Angstrom^2) | | | +| compression_time | float | Time spent on compression/decompressing image \[s\] | | | +| preprocessing_time | float | Time spent on preparing the image for analysis \[s\] | | | +| azint_time | float | Time spent on azimuthal integration \[s\] | | | +| spot_finding_time | float | Time spent on spot finding \[s\] | | | +| indexing_time | float | Time spent on indexing \[s\] | | | +| refinement_time | float | Time spent on refinement of indexing solution and experimental geometry \[s\] | | | +| index_analysis_time | float | Time spent on analyzing idnexing solution, calculating profile radius and mosaicity \[s\] | | | +| bragg_prediction_time | float | Time spent on predicting Bragg spots \[s\] | | | +| integration_time | float | Time spent on Bragg integration \[s\] | | | +| image_scale_time | float | Time spent on on-the-fly scaling \[s\] | | | +| processing_time | float | Total processing time \[s\] | | | +| xfel_pulse_id | uint64 | Bunch ID (for pulsed source, e.g., SwissFEL) | | X | +| xfel_event_code | uint64 | Event code (for pulsed source, e.g., SwissFEL) | | X | +| lattice_type | object | Bravais lattice classification of the indexing result (present only if available) | | X | +| - 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 | | | +| jf_info | uint64 | Detector info field | | | +| receiver_aq_dev_delay | uint64 | Receiver internal delay | | | +| receiver_free_send_buf | uint64 | Receiver internal number of available buffer locations | | | +| receiver_buf_in_sending | uint64 | Receiver internal number of buffer locations currently in sending/writing | | | +| receiver_buf_in_preparation | uint64 | Receiver internal number of buffer locations currently in processing | | | +| storage_cell | uint64 | Storage cell number | | | +| saturated_pixel_count | uint64 | Saturated pixel count | | | +| pixel_sum | uint64 | Sum of all pixels, excl. error and saturation | | | +| error_pixel_count | uint64 | Error pixel count | | | +| strong_pixel_count | uint64 | Strong pixel count (first stage of spot finding) | | | +| min_viable_pixel_value | int64 | Minimal pixel value, excl. error and saturation | | | +| max_viable_pixel_value | int64 | Maximal pixel value, excl. error and saturation | | | +| resolution_estimate | float | Diffraction resolution estimation \[Angstrom\] | | X | +| data_collection_efficiency | float | Image collection efficiency \[\] | | | +| packets_expected | uint64 | Number of packets expected per image (in units of 2 kB) | | | +| packets_received | uint64 | Number of packets received per image (in units of 2 kB) | | | +| bkg_estimate | float | Mean value for pixels in resolution range from 3.0 to 5.0 A \[photons\] | | | +| beam_corr_x | float | Beam center correction X applied during processing \[pixel\] | | X | +| beam_corr_y | float | Beam center correction Y applied during processing \[pixel\] | | X | +| image_scale_factor | float | Scaling result: Image scale factor (g) | | X | +| image_scale_mosaicity | float | Scaling result: Image scale mosaicity \[deg\] | | X | +| image_scale_b_factor | float | Scaling result: Image scale B factor \[Angstrom^2\] | | X | +| image_scale_cc | float | Scaling result: Image scale CC | | X | +| adu_histogram | Array(uint64) | ADU histogram | | | +| roi_integrals | object | Results of ROI calculation | | X | +| - sum | int64 | Sum of pixels in ROI area \[photons\] | | | +| - sum_square | int64 | Sum of squares of pixels in ROI area \[photons\] | | | +| - pixels | uint64 | Valid pixels in ROI area | | | +| - max_count | int64 | Highest count in ROI area \[photons\] | | | +| - x_weighted_sum | int64 | ROI pixel X position multiplied by photon count \[photons * pixels\] | | | +| - y_weighted_sum | int64 | ROI pixel Y position multiplied by photon count \[photons * pixels\] | | | +| user_data | string | Optional user defined text information - this is image_appendix serialized to JSON format | X | | +| data | Map(string -> Image) | Image | X | | ## Metadata message @@ -286,6 +288,7 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | - 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 | | +| extra_lattices | Array(Array(9*float)) | Additional indexed lattices (orientation variants); present only if found | | | 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 | | diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index e566beb9..c3c5769f 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -725,6 +725,17 @@ namespace { std::vector tmp; GetCBORFloatArray(value, tmp); message.indexing_lattice = CrystalLattice(tmp); + } else if (key == "extra_lattices") { + size_t array_len = GetCBORArrayLen(value); + CborValue array_value; + cborErr(cbor_value_enter_container(&value, &array_value)); + for (size_t i = 0; i < array_len; i++) { + std::vector tmp; + GetCBORFloatArray(array_value, tmp); + if (tmp.size() == 9) + message.extra_lattices.emplace_back(tmp); + } + cborErr(cbor_value_leave_container(&value, &array_value)); } else if (key == "indexing_time") message.indexing_time_s = GetCBORFloat(value); else if (key == "spot_finding_time") @@ -1225,6 +1236,8 @@ namespace { message.unit_cell = ProcessUnitCellElement(value); else if (key == "max_spot_count") message.max_spot_count = GetCBORUInt(value); + else if (key == "max_extra_lattices") + message.max_extra_lattices = GetCBORUInt(value); else if (key == "threshold_energy") message.threshold_energy = ProcessThresholdEnergy(value); else if (key == "image_dtype") { @@ -1408,7 +1421,18 @@ namespace { message.rotation_lattice = CrystalLattice(tmp); } else if (key == "rotation_lattice_type") message.rotation_lattice_type = GetCBORLatticeMessage(value); - else + else if (key == "extra_lattices") { + size_t array_len = GetCBORArrayLen(value); + CborValue array_value; + cborErr(cbor_value_enter_container(&value, &array_value)); + for (size_t i = 0; i < array_len; i++) { + std::vector tmp; + GetCBORFloatArray(array_value, tmp); + if (tmp.size() == 9) + message.extra_lattices.emplace_back(tmp); + } + cborErr(cbor_value_leave_container(&value, &array_value)); + } else cbor_value_advance(&value); return true; } catch (const JFJochException &e) { diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index 0d2eae87..6c9bc346 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -153,6 +153,11 @@ inline void CBOR_ENC(CborEncoder &encoder, const char* key, const std::vector& v) { + cborErr(cbor_encode_tag(&encoder, TagFloatLE)); + 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& v) { cborErr(cbor_encode_text_stringz(&encoder, key)); cborErr(cbor_encode_tag(&encoder, TagUnsignedInt8Bit)); @@ -659,7 +664,7 @@ void CBORStream2Serializer::SerializeSequenceStart(const StartMessage& message) CBOR_ENC(mapEncoder, "channels", message.channels); CBOR_ENC(mapEncoder, "max_spot_count", message.max_spot_count); - + CBOR_ENC(mapEncoder, "max_extra_lattices", message.max_extra_lattices); CBOR_ENC(mapEncoder, "storage_cell_number", message.storage_cell_number); CBOR_ENC_RATIONAL(mapEncoder, "storage_cell_delay", message.storage_cell_delay_ns, 1000*1000*1000UL); CBOR_ENC(mapEncoder, "threshold_energy", message.threshold_energy); @@ -717,7 +722,14 @@ 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()); - + if (!message.extra_lattices.empty()) { + CborEncoder arrayEncoder; + cborErr(cbor_encode_text_stringz(&mapEncoder, "extra_lattices")); + cborErr(cbor_encoder_create_array(&mapEncoder, &arrayEncoder, message.extra_lattices.size())); + for (const auto &el : message.extra_lattices) + CBOR_ENC_FLOAT_ARRAY_NOKEY(arrayEncoder, el.GetVector()); + cborErr(cbor_encoder_close_container(&mapEncoder, &arrayEncoder)); + } 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); @@ -769,6 +781,14 @@ void CBORStream2Serializer::SerializeImageInternal(CborEncoder &mapEncoder, cons CBOR_ENC(mapEncoder, "indexing_lattice_count", message.indexing_lattice_count); if (message.indexing_lattice) CBOR_ENC(mapEncoder, "indexing_lattice", message.indexing_lattice->GetVector()); + if (!message.extra_lattices.empty()) { + CborEncoder arrayEncoder; + cborErr(cbor_encode_text_stringz(&mapEncoder, "extra_lattices")); + cborErr(cbor_encoder_create_array(&mapEncoder, &arrayEncoder, message.extra_lattices.size())); + for (const auto &el : message.extra_lattices) + CBOR_ENC_FLOAT_ARRAY_NOKEY(arrayEncoder, el.GetVector()); + cborErr(cbor_encoder_close_container(&mapEncoder, &arrayEncoder)); + } CBOR_ENC(mapEncoder, "profile_radius", message.profile_radius); CBOR_ENC(mapEncoder, "mosaicity", message.mosaicity_deg); CBOR_ENC(mapEncoder, "integrated_reflections", message.integrated_reflections); diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index b72712d4..82598783 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -8,6 +8,7 @@ #include "geom_refinement/XtalOptimizer.h" #include "indexing/AnalyzeIndexing.h" #include "indexing/FFTIndexer.h" +#include "indexing/MultiLatticeSearch.h" #include "lattice_search/LatticeSearch.h" #include "scale_merge/ScaleOnTheFly.h" @@ -56,7 +57,12 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetryRotat auto gon = result->axis; if (gon) { const float angle_deg = gon->GetAngle_deg(msg.number) + gon->GetWedge_deg() / 2.0f; - outcome.lattice_candidate = result->lattice.Multiply(gon->GetTransformationAngle(-angle_deg)); + const auto rot_to_image = gon->GetTransformationAngle(-angle_deg); + outcome.lattice_candidate = result->lattice.Multiply(rot_to_image); + + outcome.extra_lattice_candidates.reserve(result->extra_lattices.size()); + for (const auto &el : result->extra_lattices) + outcome.extra_lattice_candidates.push_back(el.Multiply(rot_to_image)); } outcome.experiment.BeamX_pxl(result->geom.GetBeamX_pxl()) @@ -113,6 +119,20 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data }; outcome.lattice_candidate = sym_result.conventional; } + + // Multi-lattice search for stills: keep extra lattices that share the + // reference unit cell but differ only in orientation. Quick refinement + // (orientation only) is done later in RefineGeometryIfNeeded. + if (indexer_result.lattice.size() > 1) { + auto ml_latt = MultiLatticeSearch(indexer_result.lattice); + for (auto &ml : ml_latt) { + if (outcome.extra_lattice_candidates.size() >= experiment.GetIndexingSettings().GetMaxExtraLattices()) + break; + RotMatrix rot(ml.rotation_vector.Length(), ml.rotation_vector.Normalize()); + outcome.extra_lattice_candidates.push_back(outcome.lattice_candidate->Multiply(rot)); + } + } + } } @@ -166,6 +186,28 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In if (outcome.symmetry.crystal_system == gemmi::CrystalSystem::Monoclinic) outcome.lattice_candidate->ReorderMonoclinic(); + // Quick orientation-only refinement of extra lattices (stills path). + // Cell, beam center, detector geometry are taken from the first lattice. + if (!experiment.IsRotationIndexing() && !outcome.extra_lattice_candidates.empty()) { + for (auto &el : outcome.extra_lattice_candidates) { + XtalOptimizerData data_extra{ + .geom = data.geom, + .latt = el, + .crystal_system = data.crystal_system, + .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), + .refine_beam_center = false, + .refine_distance_mm = false, + .refine_detector_angles = false, + .refine_unit_cell = false, + .refine_rotation_axis = false, + .index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(), + .max_time = 0.02 + }; + XtalOptimizerRotationOnly(data_extra, msg.spots, 0.1); + el = data_extra.latt; + } + } + if (outcome.beam_center_updated) { msg.beam_corr_x = data.beam_corr_x; msg.beam_corr_y = data.beam_corr_y; @@ -296,7 +338,7 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, if (!outcome.lattice_candidate.has_value()) return; - if (!AnalyzeIndexing(msg, outcome.experiment, *outcome.lattice_candidate)) + if (!AnalyzeIndexing(msg, outcome.experiment, *outcome.lattice_candidate, outcome.extra_lattice_candidates)) return; { diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index eab212a0..8df60d21 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -35,6 +35,7 @@ class IndexAndRefine { struct IndexingOutcome { std::optional lattice_candidate; + std::vector extra_lattice_candidates; DiffractionExperiment experiment; LatticeMessage symmetry{ .centering = 'P', diff --git a/image_analysis/indexing/AnalyzeIndexing.cpp b/image_analysis/indexing/AnalyzeIndexing.cpp index a8a682c8..a85bd286 100644 --- a/image_analysis/indexing/AnalyzeIndexing.cpp +++ b/image_analysis/indexing/AnalyzeIndexing.cpp @@ -295,7 +295,8 @@ namespace { bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experiment, - const CrystalLattice &latt) { + const CrystalLattice &latt, + const std::vector &extra_lattices) { auto start_time = std::chrono::steady_clock::now(); std::vector indexed_spots(message.spots.size()); @@ -356,13 +357,67 @@ bool AnalyzeIndexing(DataMessage &message, assert(indexed_spots.size() == message.spots.size()); for (int i = 0; i < message.spots.size(); i++) { message.spots[i].indexed = indexed_spots[i]; - message.spots[i].lattice = 0; + message.spots[i].lattice = indexed_spots[i] ? 0 : -1; } message.profile_radius = FitProfileRadius(message.spots); message.spot_count_indexed = nspots_indexed; message.indexing_lattice = latt; message.indexing_unit_cell = latt.GetUnitCell(); message.mosaicity_deg = CalcMosaicityXDS(experiment, message.spots, astar, bstar, cstar); + + // Assign remaining (unindexed) spots to extra lattices, in order. + // Spots already assigned to the main lattice (lattice == 0) are never + // overwritten. Each extra lattice gets index 1, 2, 3, ... + const size_t n_extra = std::min(extra_lattices.size(), experiment.GetIndexingSettings().GetMaxExtraLattices()); + message.extra_lattices.clear(); + message.extra_lattices.reserve(n_extra); + + for (size_t li = 0; li < n_extra; li++) { + const CrystalLattice &el = extra_lattices[li]; + + const Coord ea = el.Vec0(); + const Coord eb = el.Vec1(); + const Coord ec = el.Vec2(); + const Coord east = el.Astar(); + const Coord ebst = el.Bstar(); + const Coord ecst = el.Cstar(); + + const int64_t lattice_id = static_cast(li) + 1; + + for (int i = 0; i < message.spots.size(); i++) { + // Do not overwrite spots already assigned to a lattice + if (message.spots[i].lattice >= 0) + continue; + + auto recip = message.spots[i].ReciprocalCoord(geom); + + float h_fp = recip * ea; + float k_fp = recip * eb; + float l_fp = recip * ec; + + float h_frac = h_fp - std::round(h_fp); + float k_frac = k_fp - std::round(k_fp); + float l_frac = l_fp - std::round(l_fp); + + float norm_sq = h_frac * h_frac + k_frac * k_frac + l_frac * l_frac; + + if (norm_sq < indexing_tolerance_sq) { + Coord recip_pred = std::round(h_fp) * east + + std::round(k_fp) * ebst + + std::round(l_fp) * ecst; + message.spots[i].indexed = true; + message.spots[i].lattice = lattice_id; + message.spots[i].dist_ewald_sphere = geom.DistFromEwaldSphere(recip_pred); + message.spots[i].h = std::lround(h_fp); + message.spots[i].k = std::lround(k_fp); + message.spots[i].l = std::lround(l_fp); + } + } + + message.extra_lattices.push_back(el); + } + + message.indexing_lattice_count = static_cast(1 + message.extra_lattices.size()); outcome = true; } } diff --git a/image_analysis/indexing/AnalyzeIndexing.h b/image_analysis/indexing/AnalyzeIndexing.h index ae0213cc..a2b347be 100644 --- a/image_analysis/indexing/AnalyzeIndexing.h +++ b/image_analysis/indexing/AnalyzeIndexing.h @@ -11,6 +11,7 @@ constexpr static float min_percentage_spots = 0.20f; bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experiment, - const CrystalLattice &latt); + const CrystalLattice &latt, + const std::vector &extra_lattices = {}); diff --git a/receiver/JFJochReceiver.cpp b/receiver/JFJochReceiver.cpp index da9367c9..ebc92658 100644 --- a/receiver/JFJochReceiver.cpp +++ b/receiver/JFJochReceiver.cpp @@ -172,6 +172,7 @@ void JFJochReceiver::SendEndMessage() { .niggli_class = rotation_indexer_ret->search_result.niggli_class, .crystal_system = rotation_indexer_ret->search_result.system }; + message.extra_lattices = rotation_indexer_ret->extra_lattices; rotation_indexing_lattice = rotation_indexer_ret->lattice; } message.unit_cell = indexer.GetConsensusUnitCell(); diff --git a/writer/HDF5DataFilePluginMX.cpp b/writer/HDF5DataFilePluginMX.cpp index e49ac866..e0318c4f 100644 --- a/writer/HDF5DataFilePluginMX.cpp +++ b/writer/HDF5DataFilePluginMX.cpp @@ -46,7 +46,9 @@ inline std::string to_bravais_code(const std::optional &lm_opt) } HDF5DataFilePluginMX::HDF5DataFilePluginMX(const StartMessage &msg) - : max_spots(msg.max_spot_count), indexing(msg.indexing_algorithm != IndexingAlgorithmEnum::None) { + : max_spots(msg.max_spot_count), + max_extra_lattices(msg.max_extra_lattices), + indexing(msg.indexing_algorithm != IndexingAlgorithmEnum::None) { } void HDF5DataFilePluginMX::OpenFile(HDF5File &data_file, const DataMessage &msg, size_t images_per_file) { @@ -74,11 +76,11 @@ void HDF5DataFilePluginMX::OpenFile(HDF5File &data_file, const DataMessage &msg, spot_l.reserve(max_spots * images_per_file); spot_lattice.reserve(max_spots * images_per_file); spot_dist_ewald.reserve(max_spots * images_per_file); - if (indexing) + if (indexing) { spot_indexed.reserve(max_spots * images_per_file); - if (indexing) indexed_lattice.reserve(images_per_file * 9); - + extra_lattices.reserve(images_per_file * max_extra_lattices * 9); + } beam_corr_x.reserve(images_per_file); beam_corr_y.reserve(images_per_file); @@ -122,10 +124,11 @@ void HDF5DataFilePluginMX::Write(const DataMessage &msg, uint64_t image_number) spot_lattice.resize(max_spots * (max_image_number + 1), -1); spot_dist_ewald.resize(max_spots * (max_image_number + 1), NAN); - if (indexing) + if (indexing) { spot_indexed.resize(max_spots * (max_image_number + 1)); - if (indexing) indexed_lattice.resize((max_image_number + 1) * 9, NAN); + extra_lattices.resize((max_image_number + 1) * max_extra_lattices * 9, NAN); + } } uint32_t spot_cnt = std::min(msg.spots.size(), max_spots); @@ -172,6 +175,18 @@ void HDF5DataFilePluginMX::Write(const DataMessage &msg, uint64_t image_number) indexed_lattice[image_number * 9 + i] = NAN; } + for (size_t li = 0; li < max_extra_lattices; li++) { + const size_t base = (image_number * max_extra_lattices + li) * 9; + if (li < msg.extra_lattices.size()) { + auto tmp = msg.extra_lattices[li].GetVector(); + for (int i = 0; i < 9; i++) + extra_lattices[base + i] = tmp[i]; + } else { + for (int i = 0; i < 9; i++) + extra_lattices[base + i] = NAN; + } + } + if (msg.lattice_type) { bravais_lattice[image_number] = to_bravais_code(msg.lattice_type); niggli_class[image_number] = msg.lattice_type->niggli_class; @@ -228,6 +243,11 @@ void HDF5DataFilePluginMX::WriteFinal(HDF5File &data_file) { if (!indexed_lattice.empty()) data_file.SaveVector("/entry/MX/latticeIndexed", indexed_lattice, {(hsize_t) (max_image_number + 1), 9}) ->Units("Angstrom"); + if (!extra_lattices.empty()) + data_file.SaveVector("/entry/MX/extraLattices", extra_lattices, + {(hsize_t) (max_image_number + 1), (hsize_t) max_extra_lattices, 9}) + ->Units("Angstrom"); + if (!bkg_estimate.empty()) data_file.SaveVector("/entry/MX/bkgEstimate", bkg_estimate.vec()); if (!profile_radius.empty()) diff --git a/writer/HDF5DataFilePluginMX.h b/writer/HDF5DataFilePluginMX.h index c239de96..8a5dd6ef 100644 --- a/writer/HDF5DataFilePluginMX.h +++ b/writer/HDF5DataFilePluginMX.h @@ -8,6 +8,7 @@ class HDF5DataFilePluginMX : public HDF5DataFilePlugin { const size_t max_spots; + const size_t max_extra_lattices; const bool indexing; // spots std::vector spot_x; @@ -35,6 +36,7 @@ class HDF5DataFilePluginMX : public HDF5DataFilePlugin { // indexing AutoIncrVector indexed; std::vector indexed_lattice; + std::vector extra_lattices; // [nimages, max_extra_lattices, 9], NaN-filled // crystal AutoIncrVector profile_radius{NAN};