diff --git a/common/IndexingSettings.cpp b/common/IndexingSettings.cpp index 401e9b050..ce5ad34d6 100644 --- a/common/IndexingSettings.cpp +++ b/common/IndexingSettings.cpp @@ -44,10 +44,16 @@ IndexingSettings &IndexingSettings::Algorithm(IndexingAlgorithmEnum input) { return *this; } +// The accepted range, not the default (which stays 500 A - see IndexingSettings.h). The FFT can +// only recover a basis vector up to this length, so the ceiling is exactly the longest cell the +// indexer can ever find; at 500 it excluded 1091 of the PDB's 206950 X-ray entries (0.53%) outright. +// 1200 A leaves 8. Nothing pays for the wider range: the histogram is sized from the VALUE in use, +// and only a caller that asks for more - a given cell that needs it, or the long-axis rescue - gets +// a longer transform. IndexingSettings &IndexingSettings::FFT_MaxUnitCell_A(float input) { check_finite("FFT indexing max unit cell (A)", input); check_min("FFT indexing max unit cell (A)", input, 50); - check_max("FFT indexing max unit cell (A)", input, 500); + check_max("FFT indexing max unit cell (A)", input, fft_max_unit_cell_limit_A); fft_max_unit_cell_A = input; return *this; } diff --git a/common/IndexingSettings.h b/common/IndexingSettings.h index b387c3106..34e2d864b 100644 --- a/common/IndexingSettings.h +++ b/common/IndexingSettings.h @@ -26,6 +26,14 @@ class IndexingSettings { float max_angle_from_ewald_deg = 2.0; float unit_cell_dist_tolerance_vs_reference = 0.05; // relative static constexpr float unit_cell_angle_tolerance_deg = 5.0; // degree +public: + // The longest cell the FFT search can be asked to reach, and so the longest it can ever find: + // FFTIndexer sizes its histogram from fft_max_unit_cell_A and the transform's last usable bin IS + // that length. Callers that widen the bound (a given cell, the long-axis rescue) must clamp to + // this rather than let the setter throw - a rescue that recovers an implausible axis must not + // take the whole run down with it. + static constexpr float fft_max_unit_cell_limit_A = 1200.0; +private: int64_t indexing_threads = 4; // Threads splitting the candidate-cell refinement WITHIN one indexer call. 1 (the default) is the // right answer whenever indexers already run one per image across all workers; it is raised only diff --git a/image_analysis/indexing/FFTIndexerCPU.cpp b/image_analysis/indexing/FFTIndexerCPU.cpp index 92b0bd55f..8061b62f5 100644 --- a/image_analysis/indexing/FFTIndexerCPU.cpp +++ b/image_analysis/indexing/FFTIndexerCPU.cpp @@ -110,8 +110,13 @@ void FFTIndexerCPU::ExecuteFFT(const std::vector &coord, size_t nspots) { double len = len_coeff * static_cast(j); if (len <= static_cast(min_length_A)) continue; - const int lo = std::max(0, j - bg_half); - const int hi = std::min(out_len, j + bg_half + 1); + // Slide the background window inward at the ends rather than truncating it. A peak within + // bg_half bins of either end - which is where the LONGEST cells sit, the last usable bin + // being max_length_A itself - otherwise gets its background from a one-sided window, and + // the prominence it is judged on is biased by however much the spectrum slopes there. + int lo = j - bg_half, hi = j + bg_half + 1; + if (lo < 0) { hi = std::min(out_len, hi - lo); lo = 0; } + if (hi > out_len) { lo = std::max(0, lo - (hi - out_len)); hi = out_len; } const double background = (pref[hi] - pref[lo]) / static_cast(hi - lo); const double prominence = mag[j] - background; diff --git a/image_analysis/indexing/FFTIndexerGPU.cu b/image_analysis/indexing/FFTIndexerGPU.cu index e6b011f2a..22b0b1400 100644 --- a/image_analysis/indexing/FFTIndexerGPU.cu +++ b/image_analysis/indexing/FFTIndexerGPU.cu @@ -40,9 +40,18 @@ __global__ void calculate_fft_result( FFTResult result{.magnitude = 0.0f, .direction = i, .length = -1}; for (int j = 0; j < out_len; ++j) { - const int want_hi = min(j + bg_half, out_len - 1); + // Constant-width window, slid inward at the ends instead of truncated (see + // FFTIndexerCPU): a peak within bg_half of either end - where the LONGEST cells sit - + // otherwise gets a one-sided background and a biased prominence. Both bounds stay + // monotonically non-decreasing in j, so the running sum below is still valid. + int want_lo = j - bg_half; + int want_hi = j + bg_half; + if (want_lo < 0) { want_hi = min(out_len - 1, want_hi - want_lo); want_lo = 0; } + if (want_hi > out_len - 1) { + want_lo = max(0, want_lo - (want_hi - (out_len - 1))); + want_hi = out_len - 1; + } while (whi < want_hi) { ++whi; winsum += complex_abs(d_output[offset + whi]); } - const int want_lo = max(0, j - bg_half); while (wlo < want_lo) { winsum -= complex_abs(d_output[offset + wlo]); ++wlo; } const float len = len_coeff * static_cast(j);