mx/Jungfraujoch
1
1
Fork 0
Code Issues Pull Requests 1 Actions Packages Projects Releases 53 Wiki Activity

Merge: Simplify (but not yet gaining performance)

Browse Source
This commit is contained in:
leonarski_f
2026-05-13 07:48:17 +02:00
parent 2fa9293fda
commit c55e91ccfa
1 changed files with 22 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/Merge.cpp
+22 -41
View File
@@ -71,46 +71,15 @@ namespace {
return out;
}
MergeResult InitResult(const std::vector<HKLKey> &slot_to_hkl,
const std::vector<Obs> &obs) {
MergeResult 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, MergeResult &out, const std::vector<Obs> &obs) {
MergeResult Merge(const std::vector<HKLKey> &slot_to_hkl, const std::vector<Obs> &obs) {
struct Accum {
double sum_wI = 0.0;
double sum_w = 0.0;
double sum_wsigma2 = 0.0;
std::vector<double> d_values;
};
std::vector<Accum> acc(nhkl);
std::vector<Accum> acc(slot_to_hkl.size());
for (const auto &o: obs) {
const double I_corr = static_cast<double>(o.r->I) * o.r->scaling_correction;
@@ -126,16 +95,31 @@ namespace {
a.sum_wI += w * I_corr;
a.sum_w += w;
a.sum_wsigma2 += w * w * sigma_corr * sigma_corr;
if (std::isfinite(o.r->d) && o.r->d > 0.0f)
a.d_values.push_back(o.r->d);
}
for (int h = 0; h < static_cast<int>(nhkl); ++h) {
const auto &a = acc[h];
MergeResult out;
out.merged.resize(slot_to_hkl.size());
for (int h = 0; h < static_cast<int>(slot_to_hkl.size()); ++h) {
auto &a = acc[h];
if (a.sum_w <= 0.0)
continue;
out.merged[h].h = slot_to_hkl[h].h;
out.merged[h].k = slot_to_hkl[h].k;
out.merged[h].l = slot_to_hkl[h].l;
out.merged[h].I = a.sum_wI / a.sum_w;
out.merged[h].sigma = std::sqrt(a.sum_wsigma2) / a.sum_w;
if (a.d_values.empty())
continue;
std::ranges::nth_element(a.d_values,
a.d_values.begin() + static_cast<long>(a.d_values.size() / 2));
out.merged[h].d = a.d_values[a.d_values.size() / 2];
}
return out;
}
void Stats(const DiffractionExperiment &x, MergeResult &out, const std::vector<Obs> &obs) {
@@ -281,10 +265,7 @@ MergeResult MergeReflections(const std::vector<std::vector<Reflection> > &observ
const DiffractionExperiment &x) {
std::vector<HKLKey> slot_to_hkl;
auto obs = BuildObservations(observations, x, slot_to_hkl);
auto out = InitResult(slot_to_hkl, obs);
Merge(slot_to_hkl.size(), out, obs);
auto out = Merge(slot_to_hkl, obs);
Stats(x, out, obs);
return out;