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

Improve handling/calculating merge mask

Browse Source
This commit is contained in:
leonarski_f
2026-05-17 14:44:33 +02:00
parent 1bf70dcae7
commit bb40c9518b
5 changed files with 42 additions and 31 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/IndexAndRefine.cpp
+7 -11
View File
@@ -304,15 +304,11 @@ ScalingResult IndexAndRefine::ScaleAllImages(const std::vector<MergedReflection>
return result;
}
MergeResult IndexAndRefine::Merge(bool calc_statistics, bool apply_cc_limit) const {
MergeResult out;
std::vector<uint8_t> merge_mask(reflections.size(), 1);
if (apply_cc_limit)
CalcMergeMask(experiment, scale_cc, merge_mask);
out.merged = MergeAll(experiment, reflections, merge_mask);
if (calc_statistics)
out.statistics = MergeStats(experiment, out.merged, reflections, merge_mask);
return out;
const std::vector<std::vector<Reflection> > &IndexAndRefine::GetReflections() const {
return reflections;
}
const std::vector<float> &IndexAndRefine::GetImageCC() const {
return scale_cc;
}
image_analysis/IndexAndRefine.h
+2 -1
View File
@@ -68,7 +68,8 @@ public:
IndexAndRefine& ReferenceIntensities(std::vector<MergedReflection> &reference);
ScalingResult ScaleAllImages(const std::vector<MergedReflection> &reference, size_t nthreads = 0);
MergeResult Merge(bool statistics, bool apply_cc_limit) const;
std::optional<RotationIndexerResult> Finalize();
const std::vector<std::vector<Reflection>> &GetReflections() const;
const std::vector<float> &GetImageCC() const;
};
image_analysis/scale_merge/Merge.cpp
+11 -6
View File
@@ -16,19 +16,24 @@
// Need to use consistent lattice to calc resolution (which would also allow to calculate completeness)
// But this is not possible for current still workflow (while fine for rotation workflow)
void CalcMergeMask(const DiffractionExperiment &x,
size_t CalcMergeMask(const DiffractionExperiment &x,
const std::vector<float> &scale_cc,
std::vector<uint8_t> &result_mask) {
if (scale_cc.size() != result_mask.size())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Scale CC size mismatch");
size_t n_rejected = 0;
auto min_cc = x.GetScalingSettings().GetMinCCForImage();
if (min_cc <= 0.0)
return;
for (int i = 0; i < result_mask.size(); ++i)
result_mask[i] = result_mask[i] && (scale_cc[i] < min_cc);
if (min_cc > 0.0) {
for (int i = 0; i < result_mask.size(); ++i) {
if (result_mask[i] && (scale_cc[i] < min_cc)) {
result_mask[i] = 0;
++n_rejected;
}
}
}
return n_rejected;
}
std::vector<MergedReflection> MergeAll(const DiffractionExperiment &x,
image_analysis/scale_merge/Merge.h
+1 -1
View File
@@ -32,7 +32,7 @@ struct MergeResult {
MergeStatistics statistics;
};
void CalcMergeMask(const DiffractionExperiment &x,
size_t CalcMergeMask(const DiffractionExperiment &x,
const std::vector<float> &scale_cc,
std::vector<uint8_t> &result_mask);
tools/jfjoch_process.cpp
+21 -12
View File
@@ -604,6 +604,8 @@ int main(int argc, char **argv) {
if (run_scaling || !reference_data.empty()) {
logger.Info("Running scaling (mosaicity refinement) ...");
std::vector<float> scale_cc;
if (reference_data.empty()) {
// If reference data are given, there is live scaling (no need to go again)
ScalingResult scale_result(0);
@@ -611,14 +613,15 @@ int main(int argc, char **argv) {
auto scale_start = std::chrono::steady_clock::now();
for (int i = 0; i < scaling_iter; i++) {
auto iter_start = std::chrono::steady_clock::now();
auto merge_result = indexer.Merge(false, false);
scale_result = indexer.ScaleAllImages(merge_result.merged);
auto merge_result = MergeAll(experiment, indexer.GetReflections());
scale_result = indexer.ScaleAllImages(merge_result);
scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat");
auto iter_end = std::chrono::steady_clock::now();
double iter_time = std::chrono::duration<double>(iter_end - iter_start).count();
logger.Info("Scaling iteration {} took {:.3f} seconds", i, iter_time);
}
scale_cc = scale_result.image_cc;
end_msg.image_scale_factor = scale_result.image_scale_g;
end_msg.image_scale_cc = scale_result.image_cc;
end_msg.image_scale_mosaicity = scale_result.mosaicity_deg;
@@ -627,15 +630,23 @@ int main(int argc, char **argv) {
auto scale_end = std::chrono::steady_clock::now();
double scale_time = std::chrono::duration<double>(scale_end - scale_start).count();
logger.Info("Scaling completed in {:.2f} s", scale_time);
}
} else
scale_cc = indexer.GetImageCC();
auto merge_start = std::chrono::steady_clock::now();
auto merge_result = indexer.Merge(true, true);
std::vector<uint8_t> merge_mask(images_to_process, 1);
auto rejected = CalcMergeMask(experiment, scale_cc, merge_mask);
logger.Info("Rejected {} images due to low CC with reference", rejected);
auto merged_reflections = MergeAll(experiment, indexer.GetReflections(), merge_mask);
auto merged_statistics = MergeStats(experiment, merged_reflections, indexer.GetReflections(), merge_mask);
auto merge_end = std::chrono::steady_clock::now();
double merge_time = std::chrono::duration<double>(merge_end - merge_start).count();
logger.Info("Merge completed in {:.2f} s ({} unique reflections)", merge_time,
merge_result.merged.size());
merged_reflections.size());
const bool fixed_space_group = space_group || experiment.GetGemmiSpaceGroup().has_value();
@@ -653,7 +664,7 @@ int main(int argc, char **argv) {
sg_opts.min_total_compared = 100;
sg_opts.test_systematic_absences = true;
const auto sg_search = SearchSpaceGroup(merge_result.merged, sg_opts);
const auto sg_search = SearchSpaceGroup(merged_reflections, sg_opts);
logger.Info("");
{
@@ -686,9 +697,7 @@ int main(int argc, char **argv) {
} */
}
// Print resolution-shell statistics table
const auto &stats = merge_result.statistics;
stats.Print(logger);
merged_statistics.Print(logger);
{
const std::string hkl_path = output_prefix + "_intensities.hkl";
@@ -696,13 +705,13 @@ int main(int argc, char **argv) {
if (!hkl_file) {
logger.Error("Cannot open {} for writing", hkl_path);
} else {
for (const auto &r: merge_result.merged) {
for (const auto &r: merged_reflections) {
hkl_file << r.h << " " << r.k << " " << r.l << " "
<< r.I << " " << r.sigma
<< "\n";
}
hkl_file.close();
logger.Info("Wrote {} reflections to {}", merge_result.merged.size(), hkl_path);
logger.Info("Wrote {} reflections to {}", merged_reflections.size(), hkl_path);
}
}
@@ -717,7 +726,7 @@ int main(int argc, char **argv) {
cif_meta.unit_cell = rotation_indexer_ret->lattice.GetUnitCell();
cif_meta.data_block_name = output_prefix;
WriteMmcifReflections(cif_path, merge_result.merged, cif_meta);
WriteMmcifReflections(cif_path, merged_reflections, cif_meta);
logger.Info("Wrote mmCIF reflections to {}", cif_path);
} catch (const std::exception &e) {
logger.Error("Failed to write mmCIF: {}", e.what());