Convolve (Ellipsoid Fusion)

Ellipsoid convolution functions for combining geolocation estimates with outlier detection and multi-cluster support.

Overview

gri-convolve provides three functions of increasing sophistication for fusing collections of Ell (ellipsoid) objects into combined position estimates:

• convolve -- combine all input ellipsoids into a single fused result with no outlier rejection
• smart_convolve -- iteratively remove outliers by Mahalanobis distance before fusing
• cluster_convolve -- find multiple clusters within a dataset and fuse each independently

Each function operates on Ell objects from gri-ell, which pair a 3D position with a statistical covariance (or information matrix). The output is one or more fused Ell objects representing the combined position estimate and its uncertainty.

Recursive single-target tracking -- the IMM / SmartSegmentedIMM filters, the motion-model bank, EKF/UKF observable updates, and RTS smoothing -- is provided by the companion gri-kalman package. Convolution is the batch face and the filter is the recursive face of the same estimation problem: the convolver seeds and refines tracks, the filter does the online sequential update. See the gri-kalman README for the tracker API and EKF-vs-UKF guidance.

Requires Python 3.12+.

Mathematical Background

Information matrix fusion. Given N ellipsoids, each with position x_k and information matrix I_k (the inverse of the covariance matrix, in XYZ coordinates, 1/m^2, 1-sigma), the fused position and information matrix are:

S = sum(I_k)              (combined information matrix)
x = S^{-1} sum(I_k x_k)  (fused position)

This is the maximum-likelihood estimator under Gaussian assumptions.

Inflation methods. The raw fusion above underestimates uncertainty when inputs are inconsistent. Three modes control how the output covariance is inflated:

• "none" -- strict information matrix combination (no inflation); calibrated when the inputs' own covariances are
• "std" -- scale by the inputs' Mahalanobis scatter and add the sample covariance of their offsets
• "bart" (default) -- as "std", but each input's horizontal offset is projected onto that input's own minor axis first, so scatter along a major axis, where it is expected, adds nothing

Both inflating modes are conservative on consistent inputs (about 3.5x and 4.6x the calibrated variance for "bart" and "std" in a Monte Carlo; see the shipped overview.md).

Outlier detection. smart_convolve measures the fused point in each input's own metric:

d_k = sqrt((x_k - mu)^T I_k (x_k - mu) / 7.815)

where mu is the (uninflated) fused position, I_k is input k's own information matrix, and 7.815 is the 3-degree-of-freedom 95% chi-square, so d_k = 1 is the surface of input k's 95% ellipsoid. Inputs with d_k > max_norm are removed worst first, refitting each time.

Reference: Mahalanobis, P.C. (1936). "On the generalized distance in statistics."

Documentation

The wheel ships its documentation inside the package, in gri_convolve/docs/, so it is available wherever the package is installed:

• overview.md -- the three functions, the fusion and its inflation modes, rejection, clustering, the array fast path, and input checks
• api_summary.md -- every public function and signature (generated)

Print the directory with python -c "import gri_convolve, pathlib; print(pathlib.Path(gri_convolve.__file__).parent / 'docs')". Every example in those files is run by the test suite.

Installation

pip install gri-convolve

For development:

git clone https://gitlab.com/geosol-foss/python/gri-convolve.git
cd gri-convolve
uv sync

Quick Start

from gri_convolve import convolve, smart_convolve, cluster_convolve
from gri_ell import Ell
from gri_pos import Pos
import numpy as np

# Create some ellipsoids at nearby positions (alt_95_m defaults to 50 m 1-sigma)
e1 = Ell.from_2d(Pos.LLA(40.0, -105.0, 1600), 100, 50, 45, alt_95_m=200)
e2 = Ell.from_2d(Pos.LLA(40.001, -105.001, 1610), 120, 60, 30, alt_95_m=200)
e3 = Ell.from_2d(Pos.LLA(40.0005, -104.999, 1605), 90, 45, 50, alt_95_m=200)

# Simple fusion
fused = convolve([e1, e2, e3])
print(fused.lla)            # Fused position
print(fused.ellipse.sma_95) # Fused semi-major axis (95%, meters)

Data that is already columns of latitude, longitude, altitude and 2D ellipse values goes straight in, with no Ells and no coordinate conversion:

from gri_convolve import smart_convolve_lla

# one value per input, or a scalar shared by all; semi-axes in meters at 95%
fused, used, discarded = smart_convolve_lla(lat, lon, alt, sma_95_m, smi_95_m, ori_deg)

# one row (lat, lon, alt, sma_95, smi_95, ori) per input instead:
fused, used, discarded = smart_convolve_lla(*np.asarray(rows).T)

Every function comes in three forms that give identical results:

convolve()

Fuses all input ellipsoids into a single result. No outlier detection.

fused = convolve(ells, inflation="bart")

Parameters:

• ells -- sequence or generator of Ell objects
• inflation -- "none", "std", or "bart" (default: "bart")

Returns: A single fused Ell.

smart_convolve()

Fuses with iterative outlier rejection. Computes the fused point, finds the input with the largest normalized Mahalanobis distance, and removes it if it exceeds max_norm. Repeats until all remaining inputs are within tolerance or fewer than min_pts remain.

result = smart_convolve(ells, max_norm=2.0, min_pts=3)
if result is not None:
    fused_ell, used_indices, discarded_indices = result

Parameters:

• ells -- sequence or generator of Ell objects
• max_norm -- maximum allowed normalized distance (default: 2.0)
• min_pts -- minimum inputs required for a valid result (default: 3)
• check_scale -- warn on a likely unit mismatch or near-total rejection (default: True; see Unit sanity checks)

Returns: (Ell, list[int], list[int]) or None if no valid cluster is found.

Pre-cluster your data before calling smart_convolve. Without pre-clustering, a large group of scattered noise points can cause valid clusters to be discarded first.

cluster_convolve()

Finds multiple clusters within a dataset by peeling them off densest first (batch=True, the default) or one at a time with smart_convolve (batch=False, the exact reference path). Leftover inputs go to the discard list; the used lists and the discard list together hold every input index exactly once.

locations, used_per_location, discarded = cluster_convolve(
    ells,
    max_norm=2.0,
    min_pts=3,
    max_pts=10,
    min_sma_m=50.0,
)

for loc, indices in zip(locations, used_per_location):
    print(f"Cluster at {loc.lla} using {len(indices)} inputs")

Parameters:

• ells -- sequence or generator of Ell objects
• max_norm -- maximum normalized Mahalanobis distance (default: 2.0)
• min_pts -- minimum inputs per cluster (default: 3)
• max_pts -- maximum inputs per cluster; splits larger groups (default: None)
• min_sma_m -- minimum semi-major axis for output ellipsoids in meters (default: 0)
• max_ori_spread -- sort by orientation before splitting for diversity (default: True)
• check_scale -- warn on a likely unit mismatch (default: True; see Unit sanity checks)
• min_cluster_size, min_cluster_frac, min_utilization -- size floor for returned clusters, and a coverage floor that keeps peeling (batch path only)
• alt_post_process -- callback for altitude correction (e.g., snap to terrain); gri_convolve.altitude.nearest is provided

Returns: (list[Ell], list[list[int]], list[int])

Unit sanity checks

The convolution math is scale-invariant -- a Mahalanobis distance is dimensionless, so entering every quantity in kilometers instead of meters changes nothing. Trouble appears only when the position scale and the covariance scale disagree. The common mistake is positions in meters (positions always resolve to ECEF meters) paired with covariance sigmas typed in kilometers: the information matrices come out about a million times too large, every residual distance is about a thousand times too large, so smart_convolve and cluster_convolve reject nearly every point one at a time -- slow, and with an almost-empty result.

To catch this, smart_convolve and cluster_convolve run two cheap checks on the full input (controlled by check_scale, default True) and emit a ConvolveScaleWarning when either trips:

• Scale pre-flight (before the rejection loop): if the tightest-packed inputs still sit many 1-sigma widths from their nearest neighbor -- meaning no cluster can form -- it warns that the sigmas look too small for the point spacing, and flags a ratio near 1000 as a likely kilometer/meter mix-up. Because it measures nearest-neighbor spacing (a local quantity), legitimate multi-site data with many separate clusters does not trip it.
• Utilization backstop (after): if nearly all inputs land in the discard list, it warns that most of the data was rejected. (cluster_convolve suppresses this when a cluster size floor is set, since trimming small clusters discards points by design.)

The checks never change the result -- they only warn. Filter or catch them with the exported ConvolveScaleWarning category, or pass check_scale=False to silence them (e.g. for data that is intentionally spread far in sigma units):

import warnings
from gri_convolve import ConvolveScaleWarning, smart_convolve

with warnings.catch_warnings():
    warnings.simplefilter("error", ConvolveScaleWarning)  # promote to an exception
    result = smart_convolve(ells)

Degenerate inputs are rejected

Separate from the warnings above, and not optional. Every semi-axis of an input ellipsoid must be finite and greater than zero. A zero-width axis -- most often alt_95_m=0, from passing an altitude where its uncertainty belongs -- makes the covariance rank deficient, and there is no fused answer to give.

This used to fail in whichever of four ways the arithmetic happened to land in, decided by float noise in the ENU-to-ECEF rotation at that latitude and longitude: a bare LinAlgError: Singular matrix at build time, the same error much later from the inflation step, a silent None with every point discarded, or a silently wrong answer built from a subset. It now raises ConvolveInputError up front, naming the offending rows:

from gri_convolve import ConvolveInputError, smart_convolve

try:
    result = smart_convolve(ells)
except ConvolveInputError as exc:
    print(exc)
    # 1 of 5000 input covariances have a variance of zero or less on the
    # diagonal: rows [4999]. A zero-width axis (commonly alt_95_m=0) has no
    # invertible covariance; the convolve solves in 3D, so even a pinned
    # altitude needs a real uncertainty.

ConvolveInputError subclasses numpy.linalg.LinAlgError (itself a ValueError), so code already guarding a convolve call against a singular matrix keeps working and simply gets a message that says which input is at fault.

A zero altitude uncertainty is never replaced by the default: an input that states zero is refused. The default is for inputs that state nothing -- see below.

Altitude uncertainty

The solver works in 3D, but most inputs are 2D ellipses. An Ell built without alt_95_m gets DEFAULT_ALT_SIGMA_M, 50 m at 1 sigma (alt_95_m 97.998 m). Every convolve function taking Ells, and ells_to_arrays, accepts alt_sigma_m (1 sigma) or alt_95_m (95%), one value or one per input:

• Nothing passed: inputs with no altitude uncertainty get the default; the rest keep their own.
• A value passed: it fills those, and also replaces the ones inputs carry, with a ConvolveAltitudeWarning naming the rows it changed. An input whose 3D covariance correlates horizontal and vertical keeps its own (and warns): replacing only its Up variance would leave cross terms that no longer match.

The altitude itself belongs in the position (Pos.LLA(lat, lon, alt)), never in alt_95_m.

Eccentricity is not degeneracy: a semi-major/semi-minor ratio of 1e7 is accepted. params_to_arrays builds the information matrix directly from the ellipse parameters instead of inverting a covariance, so it keeps full float64 precision at any ratio, where inverting loses accuracy as the square of the ratio.

Units and Conventions

• Positions are in ECEF XYZ (meters) internally
• Information matrices are in XYZ, 1/m^2, 1-sigma
• Covariance matrices are in ENU, m^2, 1-sigma
• Output ellipse parameters (SMA, SMI, orientation) are at 95% confidence
• Mahalanobis distances are normalized to 95% scale for max_norm comparisons

Dependencies

• gri-ell: Ellipsoid objects with position and covariance
• gri-pos: Position objects (XYZ, LLA coordinates)
• gri-utils: Coordinate conversions and constants
• numpy: Array operations
• scipy: Nearest-neighbor search (KDTree) for the unit sanity check

Other Projects

Current list of other GRI FOSS Projects we are building and maintaining.

License

MIT License. See LICENSE for details.