SpotAnalyze: Estimate resolution based on visible non-ice ring spots

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(),

image_analysis/spot_finding/SpotUtils.cpp

@@ -47,3 +47,22 @@ void FilterSpotsByCount(std::vector<SpotToSave> &input, int64_t count) {
                               });
     input.resize(output_size);
 }
+
+std::optional<float> GetResolution(const std::vector<SpotToSave> &spots) {
+    if (spots.size() < 6)
+        return std::nullopt;
+
+    std::vector<float> 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<size_t>(spots.size() * 0.95)];
+}

image_analysis/spot_finding/SpotUtils.h

@@ -19,4 +19,9 @@ void MarkIceRings(std::vector<SpotToSave> &spots, float tolerance_q_recipA);
 
 void FilterSpotsByCount(std::vector<SpotToSave> &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<float> GetResolution(const std::vector<SpotToSave> &spots);
+
 #endif //JFJOCH_SPOTANALYSIS_H