US Speckle Models — Part 4: Correlated Speckle and Why iid Breaks

Part 3 ended on a specific claim: local Nakagami m estimation turns speckle statistics into a tissue characterization map. The math is elegant, the interpretation is clinically motivated, and the parameter maps look convincing when you plot them.

There is a second problem, independent of the log-compression issue we discussed in Part 3, and it's the one that quietly determines how much you can trust the m maps you produce.

The problem: the standard estimator treats pixels as independent samples. They are not. The point spread function of the imaging system couples every pixel to its neighbors, and the correlation length of that coupling is often comparable to the analysis window you're estimating m over. When that happens, the estimator you're using doesn't measure what you think it measures.

This is the part of the story my PhD thesis actually addressed. Everything up to this point was setup. This is the payload.

The independence assumption, stated cleanly

Any local statistics estimator — method of moments, maximum likelihood, or anything derived from them — carries an implicit assumption that the samples inside the analysis window are drawn independently from the same distribution. The Nakagami moment estimator from Part 3:

$$\hat{m} \;=\; \frac{\big(\mathbb{E}[X^2]\big)^2}{\mathrm{Var}(X^2)}, \qquad \hat{\omega} \;=\; \mathbb{E}[X^2]$$

works exactly when the samples

$$X_1, X_2, \ldots, X_N$$

inside the window are iid. Under iid, the sample moments converge to the population moments at the standard

$$1/\sqrt{N}$$

rate, and the estimator inherits that convergence.

Ultrasound violates the iid assumption in a structural way. The received signal at each pixel is the convolution of the underlying tissue reflectivity with the imaging system's point spread function. Neighboring pixels are not independent draws — they are overlapping weighted sums of the same underlying scatterer contributions.

Effective sample size drops. The estimator sees far fewer independent observations than the window contains. The variance of

$$\hat{m}$$

is systematically underestimated, and the mean of

$$\hat{m}$$

carries a bias whose shape depends on the PSF.

What the correlation actually looks like

Take a real B-mode frame. Pick a window inside a visually homogeneous region — parenchyma, no boundaries, no obvious structural transitions. Compute the 2D autocorrelation function of the intensity, normalized to unity at zero lag. Look at the result.

Under an iid model, the autocorrelation would be a delta function: unity at zero lag, exactly zero everywhere else. What you actually see is a finite-width bump — anisotropic, tighter along the axial direction, wider along the lateral direction. That bump is the imaging system's PSF, imprinted directly into the statistics of the image.

The Wiener–Khinchin theorem makes this exact, under a wide-sense stationary assumption for the local patch. The autocorrelation of the observation is the inverse Fourier transform of its power spectral density, and if the observation is a convolution of the underlying reflectivity with a linear PSF plus additive noise, the observation autocorrelation factors into contributions from each:

$$R_y(\tau) \;=\; R_x(\tau) \,*\, R_h(\tau) \,+\, R_n(\tau)$$

where

$$R_y$$

is the observation autocorrelation,

$$R_x$$

is the underlying reflectivity autocorrelation,

$$R_h$$

is the PSF autocorrelation, and

$$R_n$$

is the additive noise autocorrelation (electronic noise, typically much smaller than the speckle-driven term). If the underlying reflectivity in a homogeneous region is approximately white and the additive noise is small,

$$R_x$$

is a delta, and

$$R_y \approx R_h$$

— the autocorrelation of the image is essentially the autocorrelation of the PSF.

Measure the half-width of the observation autocorrelation, and you have a lower bound on the PSF's spatial support. In the parenchyma patch here, that lower bound comes out to roughly 12 pixels laterally and 6 pixels axially — an anisotropic footprint on the order of a 70-pixel neighborhood of correlated area. If your Nakagami analysis window is 32×32 or smaller, the number of effectively independent samples inside it is very small.

What this costs a naive m estimator

The consequences are concrete, and they show up as three specific failure modes:

1. Underestimated uncertainty. Confidence intervals derived from the iid assumption are too narrow. A significance test on regional m differences will reject the null too often — false positives on tissue characterization, and downstream decisions built on those false positives.

2. Smoothed boundaries. The PSF integrates across boundaries between statistically distinct regions. Windows that straddle a boundary produce m estimates that are neither the m of one side nor the other — they are a PSF-weighted mixture. Boundaries in the m map blur exactly as a function of PSF width. The tissue characterization loses spatial resolution.

3. PSF-shaped bias. In regions with strong PSF-scale intensity variation (near boundaries, near strong reflectors), the correlation between neighboring pixels systematically distorts the moment ratio the estimator relies on. The bias isn't random — it's a shape imprinted by the PSF onto the m map, and it can be mistaken for real tissue variation if you don't know what to look for. (I made this mistake myself early on — for a while I was reading PSF footprints as tissue heterogeneity in my own parametric maps, and it took actually measuring the correlation to see it wasn't real.)

Most Nakagami-based QUS papers publish local m maps without measuring the correlation length of the data they estimated on. That's not a small omission. The difference between a tissue map and a PSF-shaped bias is exactly the thing you're supposed to be measuring.

The fix isn't post-hoc

The wrong move — and the tempting one — is to treat this as a denoising problem. Take the raw m map, apply some smoothing that assumes correlated observations, and call the result cleaner. This does not work, because the PSF footprint is already baked into the observation before you started estimating. Smoothing the estimate doesn't unsmooth what the PSF did to the samples.

The right move is to fold the PSF into the observation model before you estimate anything. Restoration is not a post-processing step applied to the parameter map — it is a re-formulation of what you are estimating in the first place.

Concretely: model the observation as

$$\vec{y} \;=\; A\,\vec{x} \,+\, \vec{n}, \qquad A = DB$$

where

$$\vec{y}$$

is the observed image (in vectorized form),

$$\vec{x}$$

is the underlying high-resolution reflectivity,

$$B$$

is the blur matrix representing the PSF,

$$D$$

is a decimation matrix, and

$$\vec{n}$$

is additive noise with a known correlation structure. On log-compressed B-mode data the additive Gaussian assumption for

$$\vec{n}$$

is a workable simplification — not exact, but appropriate enough that the algebra stays tractable and the recovered images look right on tissue.

The estimator becomes the maximum-a-posteriori solution:

$$\hat{\vec{x}} \;=\; \arg\max_{\vec{x}} \; p(\vec{y} \mid \vec{x}) \, p(\vec{x})$$

The likelihood carries the PSF and the noise correlation structure. Under a correlated Gaussian noise model with covariance matrix

$$C$$

, the likelihood is a multivariate Gaussian, and the exponent of the negative log-posterior comes out to:

$$(\vec{y} - A\vec{x})^T \, C^{-1} \, (\vec{y} - A\vec{x}) \;+\; \lambda \, \Phi(\vec{x})$$

where

$$\Phi(\vec{x})$$

is a prior penalizing local dissimilarities in the reflectivity, and

$$\lambda$$

controls the tradeoff.

The PSF is not a nuisance in this formulation — it is inside

$$A$$

. The pixel correlations are not ignored — they are inside

$$C^{-1}$$

. Every pixel in the analysis window contributes to the estimate without pretending it is independent of its neighbors.

The practical challenge is inverting

$$C$$

for a full image. It is a large matrix, and a naive inversion is not tractable. But

$$C$$

is banded and sparse — the correlation dies off quickly with distance — and the operation you actually need in gradient descent is

$$C^{-1}\vec{r}$$

for the current residual, which can be computed by an inner iterative solve. On a GPU this becomes feasible.

What my thesis actually built

The full formulation is in Ultrasonic Imaging, 2019. The short summary, in one paragraph:

I built a Bayesian restoration framework for freehand B-mode ultrasound that folds three ingredients into a single observation model: a measured PSF (not one assumed from theoretical beam geometry), a spatially correlated speckle noise model (with correlation structure estimated from data), and an image prior appropriate for tissue-like reflectivity. The framework produces two things: a single-image restoration (BR-CG) that recovers the underlying reflectivity from one B-mode acquisition given the measured system model, and a multi-image super-resolution extension (BSRR-CG) that combines several freehand frames with estimated sub-pixel displacements to produce a higher-resolution reconstruction. Both were validated against classical methods (Wiener, anisotropic diffusion, bilateral filtering) on tissue-mimicking phantoms and on clinical images evaluated by five radiologists in a visual grading study.

The framework was validated on log-compressed data from a commercial hospital scanner, not on lab RF or simulated envelope data. That constraint shaped every choice in the formulation, and it is what makes the method transfer to real clinical acquisitions — the domain the published QUS literature quietly avoids.

Three things about that formulation are worth flagging, because they are not obvious from the abstract:

The PSF was measured, not assumed. A tissue-mimicking phantom with a single high-echogenicity nylon wire target — approximately a point reflector at known depth — was imaged with the same scanner and settings used for the clinical acquisitions. The PSF was extracted by dividing the Fourier transform of the observed small blurred image by the Fourier transform of the expected ideal (a one-pixel dot), and transforming back. Nine axial by seventeen lateral pixels of extent, measured, not modelled.

The correlated noise model was estimated from data, using the same frequency-domain logic. Large uniform regions of the tissue-mimicking phantom were treated as approximately white reflectivity, so their observed autocorrelation is dominated by the PSF autocorrelation — exactly the Wiener–Khinchin relationship laid out above. The autocorrelation was measured over forty-nine non-overlapping patches for statistical stability. Correlation lengths came out to three axial and eight lateral pixels at 10% of peak. In the same phantom the PSF measurement and the noise correlation measurement are two views of the same physical system — one direct (the impulse response), one energetic (the autocorrelation). Folding both into the observation model closes the loop.

The restoration was single-image before it was multi-image. Multi-frame super-resolution is a natural extension once you have a single-frame observation model, but the single-frame case is where the correlated-noise likelihood actually does its work. The multi-frame case adds a sub-pixel registration problem on top; it doesn't remove the need for the correlated observation model. In practice, BSRR-CG — the multi-image correlated-Gaussian formulation — was the only method in the study that consistently improved contrast-to-noise ratio across all target sizes on the tissue-mimicking phantom, and it received the highest scores from the radiologists in the visual grading study.

Next: how do you measure a PSF you can't calculate from first principles

The observation-model reformulation is only as good as the PSF you fold into it. If you assume a PSF and estimate parameters, you inherit every error in your PSF assumption as a bias in your parameter map. If you measure a PSF, you inherit the measurement error — but the measurement error is usually much smaller than a first-principles guess, especially for commercial scanners whose beam-forming details are not published.

The next mini-series is about measuring a PSF you can't calculate. A phantom, a point target, an FFT, and a series of decisions about what to do when the measurement doesn't match the beautiful separable Gaussian your theory wanted. It turns out that's its own story.

---

References

• The full formulation and validation: Ultrasonic Imaging, 2019
• ORCID (all publications): 0000-0003-2633-9371

---

I take on selective remote consulting on medical imaging and adjacent systems engineering problems — particularly on data from commercial hospital scanners rather than research-grade RF, and on system-level integration where imaging meets motion control. Reach out via LinkedIn if that's the kind of problem you're working on.