[RFC] Spatial feature detection for imaging MS (with optional ion mobility)

Summary

Add a feature detector for mass-spectrometry imaging (MSI) data that picks features across the spatial (x, y) pixel grid, with an optional ion-mobility (IM) dimension (timsTOF MALDI). The design generalizes the two existing 1‑D hill detectors in OpenMS — Biosaur2Algorithm (links peaks along the RT axis) and PASEFHillCentroider (links along the IM-scan axis within a frame) — to 2‑D space.

Proposed class: PeakPickerMSI (a DefaultParamHandler, MSImagingExperiment → FeatureMap).

Core concept

Biosaur2Algorithm builds "hills" by linking same‑m/z peaks across consecutive RT scans; PASEFHillCentroider does the same across IM scans inside one frame. Both link along a 1‑D, ordered separation axis, where an ion forms a roughly unimodal peak — which is why both split hills at intensity valleys to deconvolve co‑eluting/co‑drifting species.

The spatial axis pair (x, y) is not a separation axis — it is a map of where the molecule physically is. An ion image can be multimodal, fragmented, or diffuse; nothing forces a peak shape. Two consequences:

1. Sequential tail‑linking → order‑independent union‑find connected components. In 2‑D a node can be adjacent to up to 8 neighbours and two already‑grown blobs can turn out to touch, which requires merging two existing groups — exactly what union‑find's union does and what Biosaur2's single rolling prev_peak_to_hill pointer cannot.
2. No intensity‑valley split in (x, y). The chromatography‑style "split co‑occurring isobars" operation belongs only on a real separation axis — i.e. the IM axis (Stage A), where PASEFHillCentroider already does it. A spatial watershed would be image segmentation (anatomy), not chemical deconvolution, so it is deliberately not part of this design.

A feature is therefore one connected component of the 8‑adjacency graph. The same molecule present in two non‑touching tissue regions yields two features (merging same‑m/z features across regions, if ever wanted, is a separate downstream consensus step).

1‑D concept (Biosaur2Algorithm / PASEFHillCentroider)	Spatial picker	Status
RT / IM‑scan axis (ordered)	2‑D pixel grid from MSImagingGeometry	changed
Hill = path down the scan axis	Ion image = connected component (8‑adjacency)	changed
Adjacent‑scan link (prev_fast_dict)	8‑neighbour pixel adjacency	changed
Sequential greedy + banned_prev_idx_set	union‑find, union every qualifier	dropped (1‑D artifacts)
mz_step = ppm·1e-6·max_mz, bins {fm-1,fm,fm+1}	identical	faithful port
IM bin floor(im/im_tol), gate `	Δbin	≤1`
min_hill_length (scans)	min_area (connected pixels)	faithful analog
splitHills_ 1‑D valley split	none (space is not a separation axis)	dropped
detectIsotopePatterns_ RT‑profile cosine	ion‑image cosine over shared pixels	extension / optional

Pipeline

A  per-pixel IM hill centroiding   (PASEFHillCentroider::centroidFrame; skipped in no-IM mode)
B  spatial union-find link         (8-neighbour; gate m/z + IM only) -> connected components = features
   B-guard  m/z coherence guard    (bound transitive m/z drift)
   filter   min_area               (connected-pixel noise filter)
D  collapse component -> Feature
E  (optional) charge/isotopes by spatial co-localization

Data model

• MSImagingExperiment = MSExperiment + MSImagingGeometry. One MSSpectrum per pixel (Bruker TIMS imaging: 1 frame == 1 pixel == 1 spectrum).
• MSImagingGeometry: zero‑based (x,y) → spectrum_index via getSpectrumIndex/hasPixel/getPixels/getWidth/getHeight; Pixel{UInt x, UInt y, Size spectrum_index}.
• Node = (pixel × centroid). After Stage A a pixel owns a list of centroids (one per m/z species), so one pixel contributes many nodes. Fields mirror PASEFHillCentroider::Centroid plus pixel coords:

struct Node {
  UInt   x, y;              // zero-based pixel coords (MSImagingGeometry::Pixel)
  Size   spectrum_index;
  double mz;                // Centroid.mz       (intensity-weighted)
  double intensity;         // Centroid.intensity (summed)
  double im_apex;           // Centroid.im_apex   (-1 in no-IM mode)
  double im_lower, im_upper;// Centroid IM range  (IM interval-overlap test)
  double mz_lower, mz_upper;// Centroid m/z range
};

Stage A — per-pixel IM hill centroiding (faithful port; skipped in no-IM mode)

use_im is decided once per dataset via IMTypes::determineIMFormat(spec) == IM_PEAK.

for pixel in geometry.getPixels():
  spec = experiment[pixel.spectrum_index]
  if not use_im:                       // NO-IM: raw spectrum centroids are the nodes; IM gate disabled
     nodes += {x,y,idx, peak.mz, peak.int, im=-1, ...} for each peak >= min_intensity
  else:                                // IM: collapse the (m/z x IM) cloud to IM-centroids
     build flat (mz[], int[], im[]) from spec + its IM FloatDataArray   // cf. centroidPASEFData_
     cents = PASEFHillCentroider::centroidFrame(mz, int, im, n, hill_params, /*scan_indices=*/nullptr)
     nodes += {x,y,idx, c.mz, c.intensity, c.im_apex, c.im_lower, c.im_upper, c.mz_lower, c.mz_upper} for c in cents

Notes: the imaging MSSpectrum exposes no per‑peak TIMS scan index, so scan_indices = nullptr (⇒ max_scan_gap = 0). If the loader already ran hill‑centroiding (BrukerTimsFile HILL_BASED MS1), the peaks are the centroids and the bounds live in named float arrays.

Stage B — spatial union-find linking (the heart)

Edge predicate gates on m/z (ppm) + IM only — never intensity (an ion image legitimately varies in brightness across space; co‑localization is reserved for Stage E).

mz_step = mz_tol_ppm * 1e-6 * global_max_mz          // == computeHillMzStep_
bin(mz) = floor(mz / mz_step);  im_bin(im) = floor(im / im_tol)
fast_dict: per-pixel  spectrum_index -> { mz_bin -> [node ids] }   // per pixel, not 1 rolling dict

spatialEdge(A, B):                                   // A, B in DIFFERENT adjacent pixels
  if |A.mz - B.mz| / min(A.mz,B.mz) * 1e6 > mz_tol_ppm: return false   // SYMMETRIC ppm
  if use_im:
     if |im_bin(A.im_apex) - im_bin(B.im_apex)| > 1: return false      // cheap prefilter (port)
     if A.im_upper < B.im_lower or B.im_upper < A.im_lower: return false// interval overlap (authoritative)
  return true                                        // NO intensity / cosine term

linkSpatial(nodes):
  uf = UnionFind(|nodes|)
  for A in nodes:
    for (dx,dy) in 8-neighbourhood:
      int qx = (int)A.x + dx, qy = (int)A.y + dy      // signed FIRST: UInt underflow at x==0
      if qx<0 or qy<0 or not geometry.hasPixel(qx,qy): continue
      sq = geometry.getSpectrumIndex(qx,qy)
      if sq == A.spectrum_index: continue             // never union within a pixel
      for fb in {bin(A.mz)-1, bin(A.mz), bin(A.mz)+1}:
        for B in fast_dict[sq][fb]:
          if spatialEdge(A,B): uf.union(A.id, B.id)    // union ALL qualifiers; order-free, idempotent
  return componentsOf(uf, nodes)

What changes vs. the 1‑D linker: drop descending‑intensity ordering, banned_prev_idx_set, and single‑winner tie‑breaks (sequential artifacts); symmetrize the ppm denominator; build one m/z‑bin index per pixel. What ports verbatim: the mz_step/binning, the 3‑bin probe, the IM bin gate.

B-guard — m/z coherence (mandatory under the component model)

A component is defined by transitive linking, so its m/z can drift across space (A~B~…~C with |A−C| ≫ tol). After union‑find, any component whose member‑m/z spread exceeds mz_cap_ppm (default 2·mz_tol_ppm) is re‑split by 1‑D single‑link clustering on member m/z. This is m/z‑identity control — distinct from the rejected intensity‑valley split. (Alternative: seed‑anchored linking that bounds drift at join time, at the cost of reintroducing a global intensity ordering.)

filter — min_area

Drop components with fewer than min_area connected pixels — the spatial analog of min_hill_length and the primary noise filter (real analyte = spatially coherent; noise = scattered).

Stage D — collapse component → Feature

Scalar summary in native Feature slots + scalar meta; footprint / ion image carried as parallel metavalue lists, not a ConvexHull2D (avoids overloading the conventional (RT, m/z) hull axes and the setHullPoints/addPoints recomputation footgun):

sort members row-major by (y, x)                      // deterministic, reproducible lists
IntList pixel_x, pixel_y;  DoubleList pixel_intensity // index-aligned == sparse COO ion image

Feature f
f.setMZ(weighted_mean_mz)            // native (dim 1)
f.setIntensity(apex_intensity)       // or intensity_sum (intensity_type param)
f.setRT(0.0)                         // MSI: RT slot meaningless
f.setCharge(0)                       // set by Stage E if it runs
f.setOverallQuality(n_pixels)        // area in the quality slot

f.setMetaValue("pixel_x", pixel_x)            // IntList
f.setMetaValue("pixel_y", pixel_y)            // IntList
f.setMetaValue("pixel_intensity", pixel_int)  // DoubleList (optional: omit and regenerate for huge footprints)
f.setMetaValue("n_pixels", n_pixels)
f.setMetaValue("intensity_sum", intensity_sum)
f.setMetaValue("centroid_x", cx);  f.setMetaValue("centroid_y", cy)
f.setMetaValue("apex_x", ax);      f.setMetaValue("apex_y", ay)
f.setMetaValue("mz_min", mz_lo);   f.setMetaValue("mz_max", mz_hi)
if use_im:
  f.setMetaValue(Constants::UserParam::IM, im_apex)
  f.setMetaValue("IM_mean", im_mean); f.setMetaValue("IM_min", im_lo); f.setMetaValue("IM_max", im_hi)

A dense IonImage (W×H grid + mask, already in OpenMS/IMAGING/IonImage.h) is reconstructed on demand from the three lists + geometry dims; it is not stored on the feature (keeps featureXML compact).

Stage E — optional charge / isotope assignment by spatial co-localization

The only place cosine correlation appears. Structurally mirrors detectIsotopePatterns_ (same bin lookup, IM gate, cos ≥ threshold), but RT‑profile overlap → pixel‑set overlap and RT‑profile cosine → ion‑image cosine over shared pixels.

ISO_DM = Constants::C13C12_MASSDIFF_U   // = 1.0033548378
for mono in features by descending intensity, for z in cmin..cmax:
  look for features near mono.mz + k*ISO_DM/z (k=1,2,...) within itol_ppm that
    (a) share IM bin, (b) spatially overlap, (c) ion-image cosine >= cos_thresh
  keep the longest consistent series -> assign charge z; isotopes -> subordinates / meta

Parameters

Name	Stage	Default	Meaning
mz_tol_ppm	A, B	8.0	m/z link tolerance; sets mz_step
im_tol	A, B	dataset	IM bin width / gate (picker-level; not a PASEFHillCentroider::Params field)
valley_factor	A	1.3	IM-axis hill valley factor (Stage A only)
min_hill_length	A	1	min IM scans per hill
min_intensity	A	0.0	input peak floor
min_area	filter	3–4	min connected pixels per feature
mz_cap_ppm	B-guard	2·mz_tol_ppm	max component m/z spread before re-split
itol_ppm, cmin, cmax, cos_thresh	E	—	isotope/charge search (optional)

Faithful port vs. extension

• Faithful ports: Stage A centroiding (centroidFrame + Centroid collapse); Stage B mz_step/binning/3‑bin probe and the IM bin gate; min_area ≈ min_hill_length; Stage E isotope skeleton + cosineCorrelation_ + ISO_DM.
• Adapted (1‑D → 2‑D): union‑find replaces sequential linking; symmetric ppm; per‑pixel m/z‑bin index; IM interval‑overlap authority.
• Extensions (no code basis): per‑pixel geometry loop + unified Node; the m/z‑coherence guard; spatial co‑localization in Stage E; the (x,y) metavalue footprint.
• Deliberately omitted: intensity‑valley / watershed spatial split.

Hard invariant: Stage B gates on m/z + IM only. Intensity / cosine appear nowhere in spatialEdge; cosine lives only in Stage E.

Open decisions

1. im_tol units — absolute IM units (matches linkScanToHills_ im/paseftol) vs. relative (max_im·tol, matches centroidPASEFData_). Pick and document.
2. pixel_intensity — store inline (full sparse image) vs. store footprint only and regenerate intensities from the experiment for very large/diffuse features.
3. m/z‑drift bounding — post‑hoc spread cap (default, fully order‑independent) vs. seed‑anchored linking.
4. Feature::setIntensity — apex (mirrors Biosaur2) vs. summed; expose intensity_type.

Implementation sketch / next steps

1. PeakPickerMSI (DefaultParamHandler) + a small reusable UnionFind.
2. Stage A delegates to existing PASEFHillCentroider::centroidFrame.
3. Stages B + D first (no‑IM and IM paths); Stage E behind a flag.
4. Synthetic fixture (a few pixels, two m/z species, one with two IM values) to test A+B+D before wiring E.

---

Design produced collaboratively; intended as an RFC for discussion before implementation.