Jungfraujoch/image_analysis/structure_refinement/DensityBlobs.h

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

#pragma once

#include <string>
#include <vector>

#include "gemmi/grid.hpp"
#include "gemmi/model.hpp"

#include "../../common/JFJochMessages.h"

// What a model leaves unexplained: the connected regions of positive mFo-DFc density above a cutoff ("blobs"),
// each placed beside the model, measured and given a likely reading - a water, an ion, something ligand-sized,
// or an error of the model itself (a misplaced or missing atom, an alternative conformation).
//
// The readings are rules on what is measured, meant to point a person at the density, not to identify a
// compound: a blob called LIGAND is density of ligand size and shape where the model has nothing, which a
// buffer molecule or an unmodelled side chain can be as well as the compound that was soaked.
struct DensityBlob {
    enum class Kind { Water, Ion, Ligand, ModelError, Unassigned };
    Kind kind = Kind::Unassigned;
    gemmi::Position centroid;        // the copy nearest the model
    double peak_sigma = 0;           // the highest point, in the map's rms
    double volume = 0;               // A^3 above the cutoff
    double score = 0;                // the density above the cutoff integrated, sigma A^3
    double nearest_atom_distance = 0;
    std::string nearest_atom;        // e.g. "OG1 THR A42"
    int polar_contacts = 0;          // N or O of the model at 2.4-3.4 A
    int metal_contacts = 0;          // N, O or S of the model at 1.8-2.6 A
    double anomalous_sigma = 0;      // the anomalous difference map at the peak, 0 where there is none
};

const char *DensityBlobKindName(DensityBlob::Kind kind);

// The blobs as the end message carries them to the writer, for <prefix>_blobs.json.
DensityBlobsMessage DensityBlobsToMessage(const std::vector<DensityBlob> &blobs);

// The blobs of `difference` (an mFo-DFc map of the whole cell, rms `rms`) above 3 sigma that peak above 4 sigma,
// strongest (by score) first. `model`: what the map was made from, in the map's cell. `anomalous`: the anomalous
// difference map on any grid of the same cell, or null. Following the flood fill of Coot's unmodelled blobs as
// gemmi implements it (find_blobs_by_flood_fill).
std::vector<DensityBlob> FindDensityBlobs(const gemmi::Grid<float> &difference, double rms, const gemmi::Structure &model,
                                          double d_min, const gemmi::Grid<float> *anomalous = nullptr,
                                          double anomalous_rms = 0);