Jungfraujoch/image_analysis/indexing/SpindleBlindFraction.cpp

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include <algorithm>
#include <cmath>

#include "SpindleBlindFraction.h"

namespace {
    // Only rows short enough to be a plausible symmetry axis count. The cut is relative to the
    // crystal's own shortest row, not an absolute length, so it works the same for a 40 A cell and
    // a 200 A one: measured over 107 solved cells, 2.5 x the shortest row covers 88% of all
    // symmetry axes and 98% of crystals' shortest one.
    constexpr float MAX_ROW_LENGTH_RATIO = 2.5f;
    // The reference length is taken over the strong rows only. A long-cell still makes the pass
    // invent short spurious rows - measured at 0.4-0.6 of the true row's peak - and taking the
    // reference over every row would let one of those shrink the window until the real aligned row
    // fell outside it, turning a severe orientation into a reported zero.
    constexpr float MIN_MAGNITUDE_RATIO = 0.5f;
    // Where the search grid can no longer resolve the crystal's rows, the pass stops returning them
    // and starts returning short spurious ones instead, and the shortlist becomes internally
    // inconsistent: its strong rows are many times longer than its shortest entry. Measured over 384
    // frames of 22 solved crystals the ratio never exceeded 2.34 and was 1.00 at the median; on a
    // synthetic still whose cell is past the grid's reach it runs 4-13. Past this the score would be
    // reporting a cone it cannot see into, so it reports nothing instead.
    constexpr float MAX_REFERENCE_LENGTH_RATIO = 3.0f;
}

float SpindleThetaMax_deg(float wavelength_A, float d_min_A) {
    if (wavelength_A <= 0 || d_min_A <= 0)
        return 0;
    const float sin_theta = std::min(1.0f, wavelength_A / (2.0f * d_min_A));
    return static_cast<float>(std::asin(sin_theta) * 180.0 / M_PI);
}

float BlindConeSelfOverlap(float x) {
    if (x >= 1.0f)
        return 0.0f;
    if (x <= 0.0f)
        return 1.0f;
    return static_cast<float>(2.0 / M_PI) *
           (std::acos(x) - x * std::sqrt(1.0f - x * x));
}

std::optional<SpindleSeverity> SpindleBlindFraction(const std::vector<Coord> &rows,
                                                    const std::vector<float> &magnitudes,
                                                    const Coord &spindle,
                                                    float theta_max_deg) {
    if (rows.empty() || rows.size() != magnitudes.size() || theta_max_deg <= 0)
        return {};

    const float axis_length = spindle.Length();
    if (axis_length < 1e-6f)
        return {};
    const Coord axis = spindle.Normalize();

    float max_magnitude = 0;
    for (const auto &m : magnitudes)
        max_magnitude = std::max(max_magnitude, m);

    float reference_length = 0;
    for (size_t i = 0; i < rows.size(); i++)
        if (magnitudes[i] >= MIN_MAGNITUDE_RATIO * max_magnitude) {
            const float l = rows[i].Length();
            if (reference_length == 0 || l < reference_length)
                reference_length = l;
        }
    if (reference_length == 0)
        return {};

    float shortest_length = reference_length;
    for (const auto &r : rows)
        shortest_length = std::min(shortest_length, r.Length());
    if (reference_length > MAX_REFERENCE_LENGTH_RATIO * shortest_length)
        return {};

    SpindleSeverity ret;
    for (size_t i = 0; i < rows.size(); i++) {
        const float length = rows[i].Length();
        if (length <= 0 || length > MAX_ROW_LENGTH_RATIO * reference_length)
            continue;
        if (magnitudes[i] < MIN_MAGNITUDE_RATIO * max_magnitude)
            continue;
        const float cos_beta = std::min(1.0f, std::fabs(rows[i] * axis) / length);
        const float beta_deg = static_cast<float>(std::acos(cos_beta) * 180.0 / M_PI);
        // A row PERPENDICULAR to the spindle is as damaging as one along it, and far more common: a
        // 2-fold about it carries the blind cone onto the cone's opposite lobe, which the same sweep
        // leaves equally unmeasured. That is Friedel's rescue, which is no rescue - the cone is
        // double-sided. Both ends of the range are the bad case and the safe zone lies between them,
        // so the miss-angle is folded about 45 deg. Checked against a Monte-Carlo of the true
        // spherical overlap the folded form is within 0.006 to theta_max = 20 deg and 0.024 to 45 deg;
        // unfolded it is wrong by a full 1.0 at beta = 90 deg, reporting the worst case as the best.
        const float fold_deg = std::min(beta_deg, 90.0f - beta_deg);
        const float score = BlindConeSelfOverlap(fold_deg / theta_max_deg);
        if (score > ret.score || ret.row_length_A == 0) {
            ret.score = score;
            ret.row_length_A = length;
            ret.miss_angle_deg = beta_deg;
        }
    }
    return ret;
}