diff --git a/common/Reflection.h b/common/Reflection.h index f1665564..b4573c34 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -18,6 +18,8 @@ struct Reflection { float delta_phi_deg; // phi angle from XDS - difference from middle of current frame (NOT an absolute angle) float predicted_x; float predicted_y; + float observed_x; + float observed_y; float d; float I; float bkg; diff --git a/docs/CBOR.md b/docs/CBOR.md index b14d9d83..d9b0a1b1 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -139,8 +139,10 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | - h | int64 | Miller index | | | | - k | int64 | Miller index | | | | - l | int64 | Miller index | | | -| - x | float | position in x (pixels) | | | -| - y | float | position in y (pixels) | | | +| - 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) | | | diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index db1c1026..d5ff05ff 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -491,6 +491,10 @@ namespace { r.predicted_x = GetCBORFloat(map_value); else if (key == "y") r.predicted_y = GetCBORFloat(map_value); + else if (key == "obs_x") + r.observed_x = GetCBORFloat(map_value); + else if (key == "obs_y") + r.observed_y = GetCBORFloat(map_value); else if (key == "d") r.d = GetCBORFloat(map_value); else if (key == "I") diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index 8834a467..6473c0c4 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -255,6 +255,8 @@ inline void CBOR_ENC(CborEncoder &encoder, const Reflection& r) { CBOR_ENC(mapEncoder, "phi", r.delta_phi_deg); CBOR_ENC(mapEncoder, "x", r.predicted_x); CBOR_ENC(mapEncoder, "y", r.predicted_y); + CBOR_ENC(mapEncoder, "obs_x", r.observed_x); + CBOR_ENC(mapEncoder, "obs_y", r.observed_y); CBOR_ENC(mapEncoder, "d", r.d); CBOR_ENC(mapEncoder, "I", r.I); CBOR_ENC(mapEncoder, "bkg", r.bkg); diff --git a/image_analysis/bragg_integration/BraggIntegrate2D.cpp b/image_analysis/bragg_integration/BraggIntegrate2D.cpp index 44c15869..b616ecc1 100644 --- a/image_analysis/bragg_integration/BraggIntegrate2D.cpp +++ b/image_analysis/bragg_integration/BraggIntegrate2D.cpp @@ -92,6 +92,9 @@ void IntegrateReflection(Reflection &r, const T *image, const std::vector bkg_values; bkg_values.reserve(static_cast((x1 - x0 + 1) * (y1 - y0 + 1))); @@ -108,6 +111,8 @@ void IntegrateReflection(Reflection &r, const T *image, const std::vector= r_2_sq && dist_sq < r_3_sq) { if (reflection_mask[y * xpixel + x] != 0) @@ -120,13 +125,16 @@ void IntegrateReflection(Reflection &r, const T *image, const std::vector 5)) { r.bkg = Median(bkg_values); r.I = static_cast(I_sum) - static_cast(I_npixel_integrated) * r.bkg; + if (I_sum > 0) { + r.observed_x = static_cast(I_sum_x) / static_cast(I_sum); + r.observed_y = static_cast(I_sum_y) / static_cast(I_sum); + } // sigma is max of the: - // - 1 (for very small numbers) + // - 1 (for zero photons) // - Poisson noise (sqrt(I_sum)) (for in between) // - minimum_sigma_in_regards_to_i of Intensity (for very large numbers) - r.sigma = 1.0; - r.sigma = std::max(r.sigma, r.I * minimum_sigma_in_regards_to_i); + r.sigma = std::max(1.0f, r.I * minimum_sigma_in_regards_to_i); if (I_sum > 0) r.sigma = std::max(r.sigma, std::sqrt(static_cast(I_sum))); r.observed = true; diff --git a/image_analysis/bragg_prediction/BraggPrediction.cpp b/image_analysis/bragg_prediction/BraggPrediction.cpp index a05b9b39..962ae45d 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.cpp +++ b/image_analysis/bragg_prediction/BraggPrediction.cpp @@ -127,6 +127,8 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal .delta_phi_deg = delta_phi_deg, .predicted_x = x, .predicted_y = y, + .observed_x = NAN, + .observed_y = NAN, .d = d, .dist_ewald = dist_ewald_sphere, .rlp = 1.0, diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.cu b/image_analysis/bragg_prediction/BraggPredictionGPU.cu index 094912fc..b15c34bc 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.cu @@ -110,6 +110,8 @@ namespace { out.delta_phi_deg = angle_from_ewald_sphere_deg(C.S0, recip_x, recip_y, recip_z, recip_sq); out.predicted_x = x; out.predicted_y = y; + out.observed_x = NAN; + out.observed_y = NAN; out.d = 1.0f / sqrtf(recip_sq); out.dist_ewald = dist_ewald; out.rlp = 1.0f; diff --git a/image_analysis/bragg_prediction/BraggPredictionRot.cpp b/image_analysis/bragg_prediction/BraggPredictionRot.cpp index 409c2d84..c349fe3d 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRot.cpp +++ b/image_analysis/bragg_prediction/BraggPredictionRot.cpp @@ -137,6 +137,8 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys .delta_phi_deg = phi * 180.0f / static_cast(M_PI), .predicted_x = x, .predicted_y = y, + .observed_x = NAN, + .observed_y = NAN, .d = d, .dist_ewald = dist_ewald_sphere, .rlp = lorentz_reciprocal, diff --git a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu index 7c5e4699..f86061f6 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu @@ -153,6 +153,8 @@ namespace { out[count].delta_phi_deg = phi * 180.0 / M_PI; out[count].predicted_x = x; out[count].predicted_y = y; + out[count].observed_x = NAN; + out[count].observed_y = NAN; out[count].d = 1.0f / sqrtf(p0_sq); out[count].dist_ewald = dist_ewald; out[count].rlp = lorentz; diff --git a/reader/JFJochHDF5Reader.cpp b/reader/JFJochHDF5Reader.cpp index 4e1783a2..62a78c31 100644 --- a/reader/JFJochHDF5Reader.cpp +++ b/reader/JFJochHDF5Reader.cpp @@ -140,7 +140,8 @@ bool ReadReflectionsFromGroup(HDF5Object &file, auto l = file.ReadOptVector(image_group_name + "/l"); auto predicted_x = file.ReadOptVector(image_group_name + "/predicted_x"); auto predicted_y = file.ReadOptVector(image_group_name + "/predicted_y"); - + auto obs_x = file.ReadOptVector(image_group_name + "/observed_x"); + auto obs_y = file.ReadOptVector(image_group_name + "/observed_y"); auto d = file.ReadOptVector(image_group_name + "/d"); auto int_sum = file.ReadOptVector(image_group_name + "/int_sum"); auto int_err = file.ReadOptVector(image_group_name + "/int_err"); @@ -175,6 +176,14 @@ bool ReadReflectionsFromGroup(HDF5Object &file, if (image_scale_corr.size() > i) image_scale_corr_val = image_scale_corr[i]; + float obs_x_val = NAN; + float obs_y_val = NAN; + + if (obs_x.size() > i && obs_y.size() > i) { + obs_x_val = obs_x[i]; + obs_y_val = obs_y[i]; + } + Reflection r{ .h = h.at(i), .k = k.at(i), @@ -182,6 +191,8 @@ bool ReadReflectionsFromGroup(HDF5Object &file, .delta_phi_deg = delta_phi_val, .predicted_x = predicted_x.at(i), .predicted_y = predicted_y.at(i), + .observed_x = obs_x_val, + .observed_y = obs_y_val, .d = d.at(i), .I = int_sum.at(i), .bkg = bkg.at(i), diff --git a/writer/HDF5DataFilePluginReflection.cpp b/writer/HDF5DataFilePluginReflection.cpp index 362b3e6e..5a35c797 100644 --- a/writer/HDF5DataFilePluginReflection.cpp +++ b/writer/HDF5DataFilePluginReflection.cpp @@ -16,7 +16,7 @@ void HDF5DataFilePluginReflection::Write(const DataMessage &msg, uint64_t image_ std::vector h, k, l; std::vector I, sigma, d, lp; - std::vector image, phi, pred_x, pred_y, bkg, partiality, zeta, scale_factor; + std::vector image, phi, pred_x, pred_y, obs_x, obs_y, bkg, partiality, zeta, scale_factor; h.reserve(msg.reflections.size()); k.reserve(msg.reflections.size()); @@ -26,6 +26,8 @@ void HDF5DataFilePluginReflection::Write(const DataMessage &msg, uint64_t image_ d.reserve(msg.reflections.size()); pred_x.reserve(msg.reflections.size()); pred_y.reserve(msg.reflections.size()); + obs_x.reserve(msg.reflections.size()); + obs_y.reserve(msg.reflections.size()); bkg.reserve(msg.reflections.size()); lp.reserve(msg.reflections.size()); partiality.reserve(msg.reflections.size()); @@ -44,6 +46,8 @@ void HDF5DataFilePluginReflection::Write(const DataMessage &msg, uint64_t image_ d.emplace_back(refl.d); pred_x.emplace_back(refl.predicted_x); pred_y.emplace_back(refl.predicted_y); + obs_x.emplace_back(refl.observed_x); + obs_y.emplace_back(refl.observed_y); bkg.emplace_back(refl.bkg); lp.emplace_back(1.0/refl.rlp); partiality.emplace_back(refl.partiality); @@ -62,6 +66,8 @@ void HDF5DataFilePluginReflection::Write(const DataMessage &msg, uint64_t image_ image_group.SaveVector("delta_phi", phi); image_group.SaveVector("predicted_x", pred_x); image_group.SaveVector("predicted_y", pred_y); + image_group.SaveVector("observed_x", obs_x); + image_group.SaveVector("observed_y", obs_y); image_group.SaveVector("int_sum", I); image_group.SaveVector("int_err", sigma); image_group.SaveVector("background_mean", bkg);