diff --git a/image_analysis/SpotAnalyze.cpp b/image_analysis/SpotAnalyze.cpp index 4766a8784..accf1fd86 100644 --- a/image_analysis/SpotAnalyze.cpp +++ b/image_analysis/SpotAnalyze.cpp @@ -31,6 +31,8 @@ void SpotAnalyze(const DiffractionExperiment &experiment, output.spots = spots_out; + output.resolution_estimate = GetResolution(output.spots); + if ((indexer != nullptr) && spot_finding_settings.indexing) { auto latt_f = indexer->Run(experiment, output); auto latt = latt_f.get(); @@ -83,10 +85,6 @@ void SpotAnalyze(const DiffractionExperiment &experiment, auto res = BraggIntegrate2D(experiment_copy, image, data.latt, ewald_dist_cutoff, output.number); - float res_estimate = sqrtf((res.b_factor.value_or(12.0f) - 12.0f) / 4.0f); - if (res_estimate > 1.0 && res_estimate < 4.0) - output.resolution_estimate = res_estimate; - constexpr size_t kMaxReflections = 10000; if (res.reflections.size() > kMaxReflections) { output.reflections.assign(res.reflections.begin(), diff --git a/image_analysis/spot_finding/SpotUtils.cpp b/image_analysis/spot_finding/SpotUtils.cpp index 6bcb809ee..f8e534ddf 100644 --- a/image_analysis/spot_finding/SpotUtils.cpp +++ b/image_analysis/spot_finding/SpotUtils.cpp @@ -47,3 +47,22 @@ void FilterSpotsByCount(std::vector &input, int64_t count) { }); input.resize(output_size); } + +std::optional GetResolution(const std::vector &spots) { + if (spots.size() < 6) + return std::nullopt; + + std::vector resolutions; + resolutions.reserve(spots.size()); + + for (const auto &spot: spots) { + if (!spot.ice_ring) + resolutions.push_back(spot.d_A); + } + + std::ranges::sort(resolutions); + + if (spots.size() < 20) + return resolutions[2]; + return resolutions[static_cast(spots.size() * 0.95)]; +} diff --git a/image_analysis/spot_finding/SpotUtils.h b/image_analysis/spot_finding/SpotUtils.h index af9f9a0f2..2395ece3f 100644 --- a/image_analysis/spot_finding/SpotUtils.h +++ b/image_analysis/spot_finding/SpotUtils.h @@ -19,4 +19,9 @@ void MarkIceRings(std::vector &spots, float tolerance_q_recipA); void FilterSpotsByCount(std::vector &input, int64_t count); +// For less than 6 spots - don't provide resolution estimation +// For less than 20 spots - take the third best high resolution +// For 20 spots or more - take the 95% percentile best high resolution spot +std::optional GetResolution(const std::vector &spots); + #endif //JFJOCH_SPOTANALYSIS_H