BMND: Direct Poisson Denoising by N-Dimensional Block Matching and Collaborative Filtering

Motivation

• Classical local and fixed-basis denoisers can smooth fine structure and often require variance-stabilizing transforms for Poisson data.The paper notes that such transforms struggle in low-count regimes.Source: §1 · p. 1

Contributions

• BMND extends block matching and collaborative filtering beyond the 2D and volumetric settings of BM3D and BM4D.The paper identifies exponential growth in block-matching complexity and memory footprint as dimensionality increases.Source: §1 · pp. 2–3
• BMND directly handles Poisson data without a variance-stabilizing transform, including at low photon counts.It also supports arbitrary-dimensional data and offers optional aggregation-aware mass conservation.Source: §1 · pp. 2–3
• BMND replaces dimension-specific reference-patch lookup tables with a traversal scheme for arrays with arbitrary numbers of axes.The implementation is released as an open-source library.Source: §1 · p. 3

Method Overview

• BMND runs hard-thresholding and Wiener stages, each combining block matching, collaborative transforms, shrinkage, inverse transforms, and aggregation.The first stage creates the pilot estimate used to guide second-stage matching and Wiener coefficients.Source: §2 · p. 3
• A common covariance framework links Gaussian and Poisson observation models to transform-domain variances used for thresholds and shrinkage.The selected patch group is treated as fixed when computing local covariance, although matching is data-dependent.Source: §2 · pp. 5, 8
• For arbitrary-dimensional data, BMND replaces dimension-specific reference-patch lookup tables with a traversal schedule balancing coverage, boundary validity, and computational cost.Matching considers SSD, Poisson deviance, Pearson discrepancy, and Anscombe SSD; low-count acceptance rules use moment calibration.Source: §2 · p. 6
• Propagated noise covariance supplies coefficient-specific variances for thresholding, including signal-dependent Poisson variance and repeated voxels in overlapping patches.Partial covariance propagation trades model accuracy against computational cost.Source: §2 · p. 6
• The Wiener stage corrects pilot-based signal-power estimates for residual noise and uses coefficient-dependent variances to filter the original noisy coefficients.Patch-level aggregation weights use propagated variance and predicted marginal risk to reflect reliability differences within a group.Source: §2 · p. 6
• Optional aggregation-aware mass conservation makes each filtered group reproduce its normalized overlap-add contribution, preserving observed total intensity after aggregation when weights are fixed.The correction accounts for aggregation weights, window, and overlap normalization.Source: §2 · p. 7

3 Observation Models

• Gaussian and Poisson models both define observation-domain covariance, which is propagated through the filtering pipeline.The Gaussian model is inherited from BM3D/BM4D, while scaled-Poisson modeling captures photon-count data stored in intensity units.Source: §3 · p. 9

3.1 Stationary Gaussian Noise

• For stationary Gaussian noise, covariance depends only on spatial displacement; white noise with an orthonormal transform gives every transformed coefficient variance σ2.Colored noise can instead produce off-diagonal patch covariance and coefficient-dependent variances.Source: §3.1 · p. 9

3.2 Scaled-Poisson Noise

• For scaled-Poisson noise, BMND uses a variance map based on a nonnegative surrogate of the clean stored-intensity mean rather than a power spectral density.The Poisson scale changes units: an intensity X represents X/sP expected events, with stored-intensity variance sP X.Source: §3.2 · p. 10

4 Block Matching

• BMND extends block matching with arbitrary-dimensional reference schedules and noise-aware matching, including finite-count calibration and correction for expected noise contributions.Local-window search, SSD comparison, and best-match retention follow BM3D/BM4D, while the cited extensions are specific to BMND.Source: §4 · p. 11

4.1 Reference-Patch Schedule

• Shifted schedules sample different patch alignments at a cost proportional to schedule passes, while BMND generates shifts from patch sizes, step sizes, and density parameters for arbitrary dimensions.This addresses the coverage limits of a fixed coarse lattice without refining every axis simultaneously.Source: §4.1 · p. 11

4.2 Metrics

• For stationary Gaussian noise, SSD is the canonical matching metric; BMND also offers Poisson deviance, Pearson discrepancy, and Anscombe SSD for Poisson observations.The section establishes SSD's likelihood-based justification for Gaussian noise and describes the other metrics as alternatives for Poisson data.Source: §4.2 · pp. 13, 16
• For Poisson data, deviance is an exact likelihood-ratio metric, whereas Pearson is a lower-cost approximation that becomes exact at high counts and Anscombe SSD is a Gaussian approximation.The Pearson statistic is also the classical chi-square test for equality of two Poisson rates [36].Source: §4.2 · pp. 15–16
• At high counts, deviance, Pearson, and normalized Anscombe contributions share a χ2_1 limit, while the Anscombe approximation can be inaccurate for low-count data.Pearson’s leading term is exact, and the Pearson distance is a computationally cheaper approximation to the exact likelihood-based deviance in the high-count regime.Source: §4.2 · pp. 17, 18, 20

4.3 Exact Finite-Count Moment Calibration

• Conditioning on the pooled count yields exact finite-count moments for each metric that depend only on the metric and count, not on the unknown rate or image.Under the equal-rate null, conditioning replaces the unknown-rate problem with allocating the observed pooled events between the two measurements.Source: §4.3 · p. 21
• Entrywise moments remain exact for any patch size, but patch-level variance depends on conditional independence and is approximate when overlapping patches reuse voxels.The patch-level variance calculation assumes conditional independence across entry pairs.Source: §4.3 · p. 23
• The reference-only threshold uses plug-in moment calibration and optional heuristic structural allowance, so it is not an exact conditional quantile and candidates cannot change their acceptance threshold.Setting the structural allowance to zero leaves only finite-count moment calibration.Source: §4.3 · pp. 23–24
• Candidate standardization gives conditional mean zero and variance one under the count model and independence assumptions, but candidate-specific centering and scaling can change distance rankings.Thus, standardized and raw-distance rankings need not agree.Source: §4.3 · pp. 24–25

4.4 Removing the Noise Bias from SSD

• For stationary zero-mean Gaussian noise, expected noise-only SSD depends on patch displacement and covariance; independence between patch entries is not required.Source: §4.4 · p. 25
• For independent scaled-Poisson observations, expected noise-only SSD depends on the clean patch means, while the plug-in correction replaces unknown means with observed nonnegative entries.Source: §4.4 · p. 26
• Noise-bias-corrected SSD subtracts a model-specific expected or estimated noise contribution from raw SSD, producing a matching score rather than a metric.Unlike count-aware distances, the corrected score retains squared stored-intensity units and may be negative.Source: §4.4 · p. 26
• Setting correction strength γ to 0 disables subtraction, γ = 1 subtracts one expected or estimated contribution, and γ > 1 applies a stronger heuristic correction.Source: §4.4 · p. 26

5 Hard Thresholding

• BMND extends the BM3D/BM4D coefficient-thresholding framework with direct-Poisson variance maps, shared-source covariance, and partial covariance propagation.It also includes an optional soft-thresholding variant.Source: §5 · p. 27
• Coefficients are retained when their magnitude exceeds a threshold scaled by their estimated noise standard deviation, making the decision relative to uncertainty rather than stored units.The hard-thresholded group is inverse-transformed to form a spatial-domain pilot for the subsequent Wiener stage.Source: §5 · p. 27
• The diagonal Poisson variance approximation treats group entries as independent, but misses dependence when overlapping patches repeat the same source voxel.The shared-source model combines the weights on each underlying random variable before covariance is computed, preserving overlap-induced dependence.Source: §5 · pp. 27–28
• Partial covariance propagation evaluates shared-source dependence for only the first selected transform planes and uses the diagonal approximation for the rest.With no selected planes it reduces to the diagonal model; with all planes selected it matches the shared-source model for all coefficients.Source: §5 · pp. 28–29
• The number of retained coefficients after hard thresholding is used to define the group-dependent aggregation weight.This connects thresholding decisions to the later overlap-add aggregation stage.Source: §5 · p. 29

6 Wiener Shrinkage

• BMND extends inherited BM3D/BM4D pilot-guided grouping and Wiener filtering with coefficient-dependent gains, bias-corrected pilot powers, and covariance-aware risk.The same ordered patch origins define the pilot and observation groups, which are transformed together before filtering.Source: §6 · p. 29

6.1 Scalar Wiener Risk and Gain

• The scalar Wiener gain balances signal attenuation against noise transmission, while the practical gain uses estimated signal power and a variance scale.Its gain lies in [0, 1], and the oracle is recovered when the variance scale is one and pilot power equals signal power.Source: §6.1 · p. 30

6.2 Estimating Signal Power from the Pilot

• Pilot-power correction subtracts an estimated noise floor, but the true pilot-error variance is unavailable and positive truncation changes the estimator’s bias.Corrected modes create a dead zone: weak true coefficients can fall inside it when the pilot is inaccurate, trading residual-noise transmission against pilot-induced attenuation.Source: §6.2 · pp. 31–32

7 Aggregation Weights

• Overlapping estimates receive reliability weights at either group level or patch level, but this inverse-reliability rule is only a local principle, not exact global optimality.Different group supports, shared noisy voxels, and shrinkage bias remain dependencies even with fuller risk surrogates.Source: §7 · pp. 32–33

7.1 Hard-Thresholding Weighting Models

• Hard-thresholding weights use coefficient masks or optional soft attenuation, with variance models propagating coefficient noise through inverse nonlocal synthesis.BMND omits a hard-thresholding risk model because no sufficiently reliable independent pilot is available before the first stage to estimate rejected signal power.Source: §7.1 · p. 34

7.2 Wiener Weighting Models

• Wiener aggregation extends classic gain-energy weights with coefficient-dependent noise variance and, in the full-risk form, rejected signal energy.For the diagonal coefficient-domain filter, mutual correlations among signal or noise coefficients do not change the risk.Source: §7.2 · pp. 34–35

7.3 Windowed Synthesized-Domain Weighting Models

• Windowed synthesized-domain weights evaluate risk after inverse synthesis and aggregation windowing, with Wiener risk using selected noise and pilot-based signal-power models.The diagonal signal model is an approximation: transform orthogonality alone does not eliminate cross-coefficient moments.Source: §7.3 · pp. 36–38

8 Aggregation-Aware Mass Conservation

• Reliability-weighted overlap-add alone does not preserve total intensity, so BMND adds a constant-synthesis correction matching each filtered group’s normalized contribution to its noisy source group.The correction depends on the group’s reliability weights, aggregation window, and pointwise overlap denominator.Source: §8 · p. 38
• With every location covered, fixed weights and denominator, and each contributing group corrected, weighted overlap-add preserves the sum of the observed realization.The theorem does not preserve unknown clean mass or impose separate constraints on arbitrary subregions.Source: §8 · p. 40
• The correction uses two sweeps: first compute filtered coefficients, reliability weights, and the denominator; then correct and aggregate using those same weights and denominator.Recomputing weights after correction could change the aggregation functional and invalidate the constraints.Source: §8 · p. 40

9.1 Datasets

• The evaluation includes generated 2D and 3D diagnostic structures, natural-image benchmarks, fluorescence microscopy, XCAT phantoms, and ECG recordings.The ECG assessment uses channel zero from all MIT–BIH records and examines reconstruction quality and beat morphology under added Gaussian noise.Source: §9.1 · p. 41

9.3 Shift Schedule Analysis

• Controlled schedule experiments compare reference-patch workload and coverage, then assess denoising quality under Gaussian noise at σ255 values 15, 25, and 50.Reference-patch count is used as a proxy for processing workload, not execution time.Source: §9.3 · p. 42, Figure 3 · p. 44
• The generated (2, 2) schedule uses approximately one third fewer reference patches than the official schedule while keeping nearly the same coverage variation.Balanced (2, 2) has the lowest coverage CV among configurations using fewer reference patches than the official schedule.Source: Figure 2 · p. 43, §9.3 · pp. 43–44
• Schedules using approximately 60–70% of the official reference count preserve or slightly improve mean PSNR, while all reported differences remain below 0.05 dB.Balanced (4, 2) has the highest PSNR in this reference-count range, but its advantage over balanced (2, 2) is at most 0.003 dB.Source: §9.3 · p. 43, Figure 3 · p. 44

9.4 Finite-count Structural Allowance Analysis

• Approximately 0.08 dB is the strongest mean PSNR improvement from positive βs, observed on generated 2D data; most other gains are below 0.04 dB, and XCAT is near zero.κs has smaller effects, generally below 0.01 dB; the experiments use βs = 0.1 and κs = 1.Source: §9.4 · p. 44, Figure 4 · p. 45

9.5 Mass Conservation

• At λpk = 0.1, whole-image bias is approximately −97% without conservation and remains close to zero with Wiener-stage conservation, while two-stage conservation reduces object-support bias to approximately −5%.Preserving the total does not ensure unbiased regional intensities.Source: §9.5 · pp. 46, 48, Figure 6 · p. 47
• At λpk = 0.1, Wiener-only conservation leaves whole-sinogram bias near zero but reduces reconstructed-image PSNR by approximately 4.4 dB; HT-only gives the largest reconstruction gains at low counts.Two-stage conservation preserves observed sinogram mass while improving reconstruction quality; benefits diminish and can become small PSNR losses at higher counts.Source: §9.5 · p. 48, Figure 7 · p. 49

9.6 Ablation Study

• All evaluated noise-aware matching strategies improve mean PSNR over SSD with fixed calibration; candidate standardization gives the largest improvements on generated, natural-image, and XCAT datasets.Poisson deviance, Pearson, and Anscombe-transformed SSD have similar marginal effects under the same calibration; fixed and reference-based finite-count calibration are nearly identical.Source: Figure 8 · p. 50, §9.6 · pp. 50–51
• Soft thresholding gives higher marginal mean PSNR than hard thresholding across dataset groups, while lowering λHT from 3 to 1.5 substantially decreases mean PSNR in every group.Increasing λHT to 4.5 produces smaller losses on generated, natural-image, and XCAT datasets, with little change on FMD.Source: Figure 9 · p. 51, §9.6 · pp. 51–52
• Both modified Wiener gains improve mean PSNR over classic gains in every dataset group, with the noise-floor formulation producing the larger improvement.Variance and risk weighting help most in the Wiener stage on natural images and XCAT; effects on generated 2D and 3D data are inconclusive, and risk weighting shows no clear added benefit over variance weighting.Source: Figure 10 · p. 52, §9.6 · p. 52, Figure 11 · p. 53
• Modified Wiener gains are strongest under the highest evaluated noise on natural images, XCAT, and FMD, while generated datasets retain substantial benefits across the count range.At low counts, candidate standardization improves over fixed and finite-count calibration on generated, natural-image, and XCAT data; direct Poisson BMND also has higher PSNR and SSIM than Anscombe BMND in each displayed example.Source: §9.6 · pp. 52, 55, Figure 12 · p. 54
• Noise-aware matching and modified Wiener gains have positive marginal mean PSNR effects across evaluated dataset groups; noise-aware Wiener aggregation gives smaller improvements, clearest on natural images and anatomical phantom data.The modified Wiener gains improve more under strong noise, while matching calibration matters most at low counts.Source: §9.6 · p. 55

9.7 FMD Denoising

• At every averaging level, modality-specific BMND exceeds published VST + BM3D and PURE-LET scores in both metrics, but remains below the best published learning-method results.The published methods were not reevaluated.Source: Table 2 · p. 57, §9.7 · p. 58
• Modality-specific profiles gain 0.55–0.76 dB PSNR over the global profile; their advantage over Anscombe BMND falls from 1.00 dB at single acquisitions to 0.13 dB at 16-image averaging.Source: §9.7 · p. 58
• For single acquisitions, direct Poisson BMND improves PSNR over noisy input by 6.82, 7.68, and 7.88 dB for confocal, two-photon, and widefield images, respectively, while Anscombe gains are similar only for the first two.The widefield advantage over Anscombe BMND decreases with averaging, from 3.91 dB to 0.30 dB PSNR and from 0.2173 to 0.0006 SSIM.Source: §9.7 · p. 58

9.8 1-Dimensional ECG Data Denoising

• BMND gives the largest mean SNR improvement at every input noise level, decreasing from 11.17 dB at 0 dB input SNR to 6.79 dB at 18 dB.Source: Figure 15 · p. 60, §9.8 · p. 60
• Across confirmation subjects, BMND outperforms every baseline after averaging over input SNR levels, with NLM and SSCF the closest competitors.The evaluation uses controlled added Gaussian noise and held-out subjects.Source: §9.8 · pp. 58, 60
• BMND has the lowest QRS waveform RMSE and relative amplitude error at every input SNR level, while signed amplitude bias remains between −1.04% and 0.16%.The table reports QRS waveform RMSE, relative amplitude error, and signed relative amplitude bias.Source: §9.8 · pp. 58, 60

• BMND formulates block-matching collaborative filtering for Gaussian and Poisson noise across one-dimensional signals, two-dimensional images, and three-dimensional volumes.Source: §10 · p. 60
• Noise-aware matching and Wiener gains improve marginal mean reconstruction quality across evaluated dataset groups, while exact covariance treatment is inconsistent and bias-aware weighting adds no clear benefit over variance weighting.Source: §10 · p. 63
• In the Shepp–Logan experiment, hard-thresholding-only conservation gives the largest low-count reconstruction PSNR gains; two-stage conservation preserves observed total intensity while improving reconstruction quality.Preserving the observed total does not guarantee unbiased regional intensities or improved PSNR.Source: §10 · p. 63
• On FMD, modality-specific BMND exceeds evaluated Anscombe and published classical scores at every averaging level but remains below the best published learning-method scores; its advantage over Anscombe is largest for single-acquisition widefield images.A common formulation still requires application-specific parameter selection.Source: §10 · p. 63
• Under controlled Gaussian noise, one-dimensional BMND improves mean SNR and QRS waveform and amplitude errors over evaluated baselines, but the ECG experiment does not assess physiological recording artifacts.Source: §10 · p. 63
• Future work proposes evaluating higher-dimensional PET projection data and explicit mixed Poisson–Gaussian modeling, which would extend beyond the present scaled-Poisson formulation and fluorescence moment-matching approximation.Source: §10 · p. 63

• The experiment parameters are taken from modern BM3D/BM4D implementations when not otherwise specified.Source: A Standard Configuration · p. 68
• Table 4 lists standard Gaussian BM3D and BM4D profile values used by BMND.Source: Table 4 · p. 68

• The ablation uses a Cartesian profile space, with level count defined as the number of independently crossed options.Unlisted profile fields retain the dimension-specific standard configuration.Source: Table 5 · p. 69, B Ablation Study · p. 69
• Figure 5 reports marginal PSNR effects of covariance-plane count k relative to k = 0, with source-unit bootstrap intervals.Source: Figure 18 · p. 70
• Figures 19–22 report marginal effects on PSNR for hard-thresholding and Wiener-stage aggregation-weight domains versus coefficient-domain weighting, and scopes versus group-level weighting.Error bars denote 95% source-unit bootstrap intervals.Source: Figure 19 · p. 70, Figure 20 · p. 71, Figure 21 · p. 71
• Figure 23 reports marginal effects of group mass conservation on PSNR relative to denoising without mass conservation, comparing conservation in the hard-thresholding stage, Wiener stage, or both.Error bars denote 95% source-unit bootstrap intervals.Source: Figure 23 · p. 72

C ECG Experiment

• Table 6 describes the parameter spaces evaluated in the controlled ECG experiment, with candidate counts limited to combinations varied during development.Fixed settings are reported for reproducibility.Source: Table 6 · p. 73
• Curvature was mapped to integer polynomial orders from 1 through the selected maximum.Source: C ECG Experiment · p. 74
• The ECG experiment used 558 total development candidates.Source: C ECG Experiment · p. 74