diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index cb59be4c..0bbedfda 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -109,6 +109,7 @@ struct DataMessage { std::optional indexing_time_s; std::optional profile_radius; + std::optional mosaicity_deg; std::optional b_factor; diff --git a/docs/CBOR.md b/docs/CBOR.md index f2130797..68f86fab 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -140,12 +140,12 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | - l | int64 | Miller index | | | | - x | float | position in x (pixels) | | | | - y | float | position in y (pixels) | | | -| - d | float | resolution \[angstrom\] | | | +| - 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) | | | +| - dist_ewald | float | distance to Ewald sphere (present only for indexed spots) | | | | 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) | | | @@ -156,7 +156,8 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | 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) | | | +| profile_radius | float | Profile radius of the image - describes distance of observed reflections from the Ewald sphere \[Angstrom^-1\] | | | +| mosaicity | float | Angular range of spots in image from a rotation scan \[degree\] | | | | b_factor | float | Estimated B-factor (Angstrom^2) | | | | indexing_time | float | Time spent on indexing \[s\] | | | | processing_time | float | Total processing time \[s\] | | | diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index 450089f4..1419c302 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -673,6 +673,8 @@ namespace { message.processing_time_s = GetCBORFloat(value); else if (key == "profile_radius") message.profile_radius = GetCBORFloat(value); + else if (key == "mosaicity") + message.mosaicity_deg = GetCBORFloat(value); else if (key == "b_factor") message.b_factor = GetCBORFloat(value); else if (key == "indexing_unit_cell") diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index 25a53729..2a1012b9 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -694,6 +694,7 @@ void CBORStream2Serializer::SerializeImageInternal(CborEncoder &mapEncoder, cons if (message.indexing_lattice) CBOR_ENC(mapEncoder, "indexing_lattice", message.indexing_lattice->GetVector()); CBOR_ENC(mapEncoder, "profile_radius", message.profile_radius); + CBOR_ENC(mapEncoder, "mosaicity", message.mosaicity_deg); CBOR_ENC(mapEncoder, "b_factor", message.b_factor); CBOR_ENC(mapEncoder, "indexing_time", message.indexing_time_s); CBOR_ENC(mapEncoder, "processing_time", message.processing_time_s); diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 0bbbafde..99e96df7 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -129,7 +129,7 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, msg.beam_corr_y = data.beam_corr_y; } - if (AnalyzeIndexing(msg, experiment_copy, *lattice_candidate)) { + if (AnalyzeIndexing(msg, experiment_copy, *lattice_candidate, experiment_copy.GetGoniometer())) { msg.lattice_type = symmetry; float ewald_dist_cutoff = 0.001f; diff --git a/image_analysis/indexing/AnalyzeIndexing.cpp b/image_analysis/indexing/AnalyzeIndexing.cpp index 81209e26..e503c81f 100644 --- a/image_analysis/indexing/AnalyzeIndexing.cpp +++ b/image_analysis/indexing/AnalyzeIndexing.cpp @@ -7,23 +7,126 @@ #include "FitProfileRadius.h" -inline bool ok(float x) { - if (!std::isfinite(x)) - return false; - if (x < 0.0) - return false; - return true; -}; +namespace { + inline bool ok(float x) { + if (!std::isfinite(x)) + return false; + if (x < 0.0) + return false; + return true; + } -bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experiment, const CrystalLattice &latt) { + inline float deg_to_rad(float deg) { + return deg * (static_cast(M_PI) / 180.0f); + } + inline float rad_to_deg(float rad) { + return rad * (180.0f / static_cast(M_PI)); + } + + // Wrap to [-180, 180] (useful for residuals) + inline float wrap_deg_pm180(float deg) { + while (deg > 180.0f) deg -= 360.0f; + while (deg < -180.0f) deg += 360.0f; + return deg; + } + + // Solve A cos(phi) + B sin(phi) + D = 0, return solutions in [phi0, phi1] (radians) + inline int solve_trig(float A, float B, float D, + float phi0, float phi1, + float out_phi[2]) { + const float R = std::sqrt(A * A + B * B); + if (!(R > 0.0f)) + return 0; + + const float rhs = -D / R; + if (rhs < -1.0f || rhs > 1.0f) + return 0; + + const float phi_ref = std::atan2(B, A); + const float delta = std::acos(rhs); + + float s1 = phi_ref + delta; + float s2 = phi_ref - delta; + + const float two_pi = 2.0f * static_cast(M_PI); + auto shift_near = [&](float x) { + while (x < phi0 - two_pi) x += two_pi; + while (x > phi1 + two_pi) x -= two_pi; + return x; + }; + + s1 = shift_near(s1); + s2 = shift_near(s2); + + int n = 0; + if (s1 >= phi0 && s1 <= phi1) out_phi[n++] = s1; + if (s2 >= phi0 && s2 <= phi1) { + if (n == 0 || std::fabs(s2 - out_phi[0]) > 1e-6f) out_phi[n++] = s2; + } + return n; + } + + // Find predicted phi (deg) for given g0 around phi_obs (deg) within +/- half_window_deg. + // Returns nullopt if no solution in the local window. + inline std::optional predict_phi_deg_local(const Coord &g0, + const Coord &S0, + const Coord &w_unit, + float phi_obs_deg, + float half_window_deg) { + const float phi0 = deg_to_rad(phi_obs_deg - half_window_deg); + const float phi1 = deg_to_rad(phi_obs_deg + half_window_deg); + + // Decompose g0 into parallel/perp to w + const float g_par_s = g0 * w_unit; + const Coord g_par = w_unit * g_par_s; + const Coord g_perp = g0 - g_par; + + const float g_perp2 = g_perp * g_perp; + if (g_perp2 < 1e-12f) + return std::nullopt; + + const float k2 = (S0 * S0); // |S0|^2 = (1/lambda)^2 + + // Equation: |S0 + g(phi)|^2 = |S0|^2 + const Coord p = S0 + g_par; + const Coord w_x_gperp = w_unit % g_perp; + + const float A = 2.0f * (p * g_perp); + const float B = 2.0f * (p * w_x_gperp); + const float D = (p * p) + g_perp2 - k2; + + float sols[2]{}; + const int nsol = solve_trig(A, B, D, phi0, phi1, sols); + if (nsol == 0) + return std::nullopt; + + // Pick the solution closest to phi_obs + const float phi_obs = deg_to_rad(phi_obs_deg); + float best_phi = sols[0]; + float best_err = std::fabs(sols[0] - phi_obs); + + if (nsol == 2) { + const float err2 = std::fabs(sols[1] - phi_obs); + if (err2 < best_err) { + best_err = err2; + best_phi = sols[1]; + } + } + + return rad_to_deg(best_phi); + } +} // namespace + +bool AnalyzeIndexing(DataMessage &message, + const DiffractionExperiment &experiment, + const CrystalLattice &latt, + const std::optional &rotation_axis) { size_t nspots = message.spots.size(); uint64_t indexed_spot_count = 0; std::vector indexed_spots(nspots); - std::vector distance_ewald_sphere(nspots); - // Check spots const Coord a = latt.Vec0(); const Coord b = latt.Vec1(); @@ -36,7 +139,8 @@ bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experime const auto geom = experiment.GetDiffractionGeometry(); const auto indexing_tolerance = experiment.GetIndexingSettings().GetTolerance(); const auto viable_cell_min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(); - // identify indexed spots + + // identify indexed spots for (int i = 0; i < message.spots.size(); i++) { auto recip = message.spots[i].ReciprocalCoord(geom); @@ -79,8 +183,54 @@ bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experime message.spot_count_indexed = indexed_spot_count; message.indexing_lattice = latt; message.indexing_unit_cell = latt.GetUnitCell(); + + message.mosaicity_deg = std::nullopt; + if (rotation_axis.has_value()) { + const auto gon_opt = experiment.GetGoniometer(); + if (gon_opt.has_value()) { + const auto &gon = *gon_opt; + + const Coord w = rotation_axis->GetAxis().Normalize(); + const Coord S0 = geom.GetScatteringVector(); + const float wedge_deg = gon.GetWedge_deg(); + const float start_deg = gon.GetStart_deg(); + const float inc_deg = gon.GetIncrement_deg(); + + double sum_sq = 0.0; + int count = 0; + + for (const auto &s: message.spots) { + if (!s.indexed) + continue; + + // Observed angle: use frame center + const float image_center = static_cast(s.image) + 0.5f; + const float phi_obs_deg = start_deg + inc_deg * image_center; + + // g0 at phi=0 assumption + const Coord g0 = astar * static_cast(s.h) + + bstar * static_cast(s.k) + + cstar * static_cast(s.l); + + // Local solve window: +/- 1 wedge (easy/robust first try) + const auto phi_pred_deg_opt = predict_phi_deg_local(g0, S0, w, phi_obs_deg, wedge_deg); + if (!phi_pred_deg_opt.has_value()) + continue; + + float dphi = wrap_deg_pm180(phi_obs_deg - phi_pred_deg_opt.value()); + sum_sq += static_cast(dphi) * static_cast(dphi); + count++; + } + + if (count > 0) { + message.mosaicity_deg = static_cast(std::sqrt(sum_sq / static_cast(count))); + } + } + } + return true; } + message.indexing_result = false; return false; } diff --git a/image_analysis/indexing/AnalyzeIndexing.h b/image_analysis/indexing/AnalyzeIndexing.h index 3a199dc5..6757f3be 100644 --- a/image_analysis/indexing/AnalyzeIndexing.h +++ b/image_analysis/indexing/AnalyzeIndexing.h @@ -10,7 +10,10 @@ constexpr static float min_percentage_spots = 0.20f; -bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experiment, const CrystalLattice &latt); +bool AnalyzeIndexing(DataMessage &message, + const DiffractionExperiment &experiment, + const CrystalLattice &latt, + const std::optional &rotation_axis); #endif //JFJOCH_ANALYZEINDEXING_H \ No newline at end of file diff --git a/writer/HDF5DataFilePluginMX.cpp b/writer/HDF5DataFilePluginMX.cpp index ddb435f7..150c9147 100644 --- a/writer/HDF5DataFilePluginMX.cpp +++ b/writer/HDF5DataFilePluginMX.cpp @@ -73,6 +73,7 @@ void HDF5DataFilePluginMX::OpenFile(HDF5File &data_file, const DataMessage &msg) strong_pixel_count.reserve(RESERVE_IMAGES); indexed.reserve(RESERVE_IMAGES); profile_radius.reserve(RESERVE_IMAGES); + mosaicity_deg.reserve(RESERVE_IMAGES); b_factor.reserve(RESERVE_IMAGES); indexed_lattice.reserve(9 * RESERVE_IMAGES); resolution_estimate.reserve(RESERVE_IMAGES); @@ -133,6 +134,7 @@ void HDF5DataFilePluginMX::Write(const DataMessage &msg, uint64_t image_number) if (indexing) { indexed[image_number] = msg.indexing_result.value_or(0); profile_radius[image_number] = msg.profile_radius.value_or(NAN); + mosaicity_deg[image_number] = msg.mosaicity_deg.value_or(NAN); b_factor[image_number] = msg.b_factor.value_or(NAN); resolution_estimate[image_number] = msg.resolution_estimate.value_or(NAN); beam_corr_x[image_number] = msg.beam_corr_x.value_or(NAN); @@ -201,6 +203,8 @@ void HDF5DataFilePluginMX::WriteFinal(HDF5File &data_file) { data_file.SaveVector("/entry/MX/bkgEstimate", bkg_estimate.vec()); if (!profile_radius.empty()) data_file.SaveVector("/entry/MX/profileRadius", profile_radius.vec())->Units("Angstrom^-1"); + if (!mosaicity_deg.empty()) + data_file.SaveVector("/entry/MX/mosaicity", profile_radius.vec())->Units("deg"); if (!b_factor.empty()) data_file.SaveVector("/entry/MX/bFactor", b_factor.vec())->Units("Angstrom^2"); diff --git a/writer/HDF5DataFilePluginMX.h b/writer/HDF5DataFilePluginMX.h index bf9f2fea..8e1b3b3c 100644 --- a/writer/HDF5DataFilePluginMX.h +++ b/writer/HDF5DataFilePluginMX.h @@ -38,6 +38,7 @@ class HDF5DataFilePluginMX : public HDF5DataFilePlugin { // crystal AutoIncrVector profile_radius; + AutoIncrVector mosaicity_deg; AutoIncrVector b_factor; // bkg_estimate