Jungfraujoch/image_analysis/structure_refinement/DensityBlobs.cpp

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

#include "DensityBlobs.h"

#include <algorithm>
#include <cmath>

#include "gemmi/blob.hpp"
#include "gemmi/neighbor.hpp"

namespace {

constexpr double CUTOFF_SIGMA = 3.0;
constexpr double PEAK_SIGMA = 4.0;
// Closer than this to a model atom, the density is the model's own atom out of place.
constexpr double MODEL_ERROR_DISTANCE = 1.6;
// Waters and ions sit at least this far from the model's atoms.
constexpr double CONTACT_MIN = 2.2;
// An ion's peak, in sigma, where the anomalous map gives no help (with two metal-like contacts); an anomalous
// peak this high says heavier than water on its own.
constexpr double ION_PEAK_SIGMA = 8.0;
constexpr double ION_ANOMALOUS_SIGMA = 4.0;
// The smallest blob called ligand-like, A^3, and at low resolution one and a half of the resolution's spheres: at
// high resolution a strong single atom's blob already spans several of them.
constexpr double LIGAND_MIN_VOLUME = 8.0;

bool Polar(const gemmi::Element &el) {
    return el == gemmi::El::N || el == gemmi::El::O;
}

}  // namespace

const char *DensityBlobKindName(DensityBlob::Kind kind) {
    switch (kind) {
        case DensityBlob::Kind::Water: return "WATER";
        case DensityBlob::Kind::Ion: return "ION";
        case DensityBlob::Kind::Ligand: return "LIGAND";
        case DensityBlob::Kind::ModelError: return "MODEL_ERROR";
        case DensityBlob::Kind::Unassigned: break;
    }
    return "UNASSIGNED";
}

// Following Coot's unmodelled blobs, Emsley et al. (2010) Acta Cryst. D66, 486-501
std::vector<DensityBlob> FindDensityBlobs(const gemmi::Grid<float> &difference, double rms, const gemmi::Structure &model,
                                          double d_min, const gemmi::Grid<float> *anomalous, double anomalous_rms) {
    std::vector<DensityBlob> out;
    if (!(rms > 0) || model.models.empty())
        return out;
    gemmi::Grid<float> map = difference;
    for (float &v : map.data)
        v = static_cast<float>(v / rms);
    gemmi::BlobCriteria criteria;
    criteria.cutoff = CUTOFF_SIGMA;
    criteria.min_volume = 0.0;
    criteria.min_score = 0.0;
    criteria.min_peak = PEAK_SIGMA;
    const std::vector<gemmi::Blob> blobs = gemmi::find_blobs_by_flood_fill(map, criteria);

    gemmi::Model copy = model.models[0];
    gemmi::UnitCell cell = map.unit_cell;
    cell.set_cell_images_from_spacegroup(map.spacegroup);
    gemmi::NeighborSearch ns(copy, cell, 5.0);
    ns.populate(false);
    // A single atom's blob at this resolution: a sphere of half the resolution's radius.
    const double atom_volume = 4.0 / 3.0 * 3.14159265358979 * std::pow(0.5 * d_min, 3);

    for (const gemmi::Blob &b : blobs) {
        DensityBlob d;
        d.peak_sigma = b.peak_value;
        d.volume = b.volume;
        d.score = b.score;
        d.centroid = b.centroid;
        // the copy of the blob beside the atom nearest to it
        if (gemmi::NeighborSearch::Mark *m = ns.find_nearest_atom(b.centroid, 12.0)) {
            const gemmi::CRA cra = m->to_cra(copy);
            const gemmi::FTransform t = ns.get_image_transformation(m->image_idx);
            const gemmi::Fractional fc = cell.fractionalize(b.centroid);
            const gemmi::Fractional fi = t.apply(cell.fractionalize(cra.atom->pos));
            const gemmi::Fractional shifted(fc.x - std::round(fc.x - fi.x), fc.y - std::round(fc.y - fi.y),
                                            fc.z - std::round(fc.z - fi.z));
            d.centroid = cell.orthogonalize(gemmi::FTransform(t.inverse()).apply(shifted));
            d.nearest_atom_distance = d.centroid.dist(cra.atom->pos);
            d.nearest_atom = cra.atom->name + " " + cra.residue->name + " " + cra.chain->name +
                             cra.residue->seqid.str();
        } else
            d.nearest_atom_distance = INFINITY;
        ns.for_each(d.centroid, '\0', 3.4, [&](gemmi::NeighborSearch::Mark &m, double dist_sq) {
            const double r = std::sqrt(dist_sq);
            if (Polar(m.element) && r >= 2.4)
                d.polar_contacts++;
            if ((Polar(m.element) || m.element == gemmi::El::S) && r >= 1.8 && r <= 2.6)
                d.metal_contacts++;
        });
        if (anomalous != nullptr && anomalous_rms > 0)
            d.anomalous_sigma = anomalous->interpolate_value(b.peak_pos) / anomalous_rms;

        if (d.nearest_atom_distance < MODEL_ERROR_DISTANCE)
            d.kind = DensityBlob::Kind::ModelError;
        else if (d.nearest_atom_distance >= 1.8 &&
                 (d.anomalous_sigma >= ION_ANOMALOUS_SIGMA || (d.metal_contacts >= 2 && d.peak_sigma >= ION_PEAK_SIGMA)))
            d.kind = DensityBlob::Kind::Ion;
        else if (d.volume > std::max(1.5 * atom_volume, LIGAND_MIN_VOLUME) && d.nearest_atom_distance >= CONTACT_MIN)
            d.kind = DensityBlob::Kind::Ligand;
        else if (d.nearest_atom_distance >= CONTACT_MIN && d.polar_contacts >= 1)
            d.kind = DensityBlob::Kind::Water;
        out.push_back(d);
    }
    std::sort(out.begin(), out.end(), [](const DensityBlob &a, const DensityBlob &b) { return a.score > b.score; });
    return out;
}