DEV Community: Mine Cuneyitoglu

A Bug That Doesn't Crash Is Worse Than One That Does — Building a QUS Python Refresher

Mine Cuneyitoglu — Mon, 24 Aug 2026 20:02:33 +0000

I recently rebuilt my quantitative ultrasound (QUS) signal processing knowledge from scratch, in Python, as a five-notebook refresher repo. Along the way I reproduced a bug worth writing up on its own, because it never raised an exception — it just quietly returned zero.

The setup

Quantitative ultrasound tries to turn raw RF echo data into numbers that describe tissue — not just an image a radiologist reads, but a spectral slope or a backscatter coefficient that correlates with scatterer size, concentration, or fibrosis. The foundational method here is the Lizzi-Feleppa framework from 1983: fit a line to the log-power spectrum within a transducer's usable bandwidth, per depth window, and read off the slope, intercept, and midband value.

To get this back into muscle memory, I built qus-python-refresher — five Jupyter notebooks that simulate an RF echo from point scatterers and walk through the full pipeline: pulse generation, envelope detection, moving-window spectral analysis, and Lizzi-Feleppa parameter estimation.

The bug

Point-scatterer echo simulation needs two separate time axes:

a short axis just long enough to hold the transmit pulse (t_pulse)
a long axis covering the full imaging depth (t), onto which delayed, scaled copies of the pulse get summed

Building the pulse on the wrong axis — reusing t instead of t_pulse — is an easy mistake. The Gaussian envelope still decays to numerical zero away from its center, so if you plot the "pulse" it looks completely normal. The array itself, though, is now the same length as the full echo instead of a few dozen samples.

The consequence shows up in the placement loop's bounds check:

if idx + len(pulse) < len(rf):
    rf[idx:idx+len(pulse)] += a * pulse

With a pulse the same length as rf, idx + len(pulse) >= len(rf) for every idx >= 0 — the condition can never be true. In a 200-scatterer test I ran to confirm this: 0 out of 200 scatterers get placed with the buggy pulse, versus 200 out of 200 with the correctly-scoped one. rf comes back as a normal, correctly-shaped, all-zero array. No exception. No warning. Just silence.

That's the failure mode worth internalizing: shape mismatches that happen to broadcast are more dangerous than the ones that crash loudly. A flat line in a B-mode image would be an easy bug to notice; buried a few layers deep in a larger pipeline, it isn't.

What's in the repo

01 — RF Pulse and Envelope. A single transmit pulse, Hilbert transform, envelope, phase, log compression. No scatterers yet, on purpose — isolates the RF-to-B-mode pipeline before the second time axis shows up.
02 — Scatterer Echo Simulation. The bug above, reproduced with real numbers. Also: diffuse (Rayleigh-regime) vs. specular scattering, and frequency-dependent attenuation.
03 — Moving-Window Spectral Analysis. scipy.signal.spectrogram's default window is Tukey(0.25), not rectangular — easy to assume otherwise if you never pass window= explicitly. This notebook compares window lengths (1/3/10 µs — a real axial-resolution-vs-frequency-resolution trade-off) and window types (rectangular vs. Hann, and the leakage difference between them).
04 — Lizzi-Feleppa Slope Estimation. Slope, intercept, midband fit, per depth window. The interesting result here: scaling echo amplitude by 5x shifts intercept and midband by very close to the theoretical 20·log10(5) = 13.98 dB — but only once you exclude depth windows that sit before the first scatterer, where both runs are noise-floor-dominated and the shift is 0 dB regardless of gain. Averaging over the full depth range without that exclusion gave a misleading ~9.9 dB, and figuring out why was a good reminder to sanity-check estimator behavior against a known-zero baseline before trusting a summary statistic. Slope itself, unlike intercept/midband, comes out essentially unchanged by the gain — which is exactly why reference phantom normalization exists.
05 — Attenuation and Spectral Downshift. Why an uncompensated spectral slope is depth-biased: attenuation removes high frequencies faster than low ones, so the local spectrum's shape — not just its amplitude — shifts with depth. Includes a spectral-centroid-vs-depth demo (the observable the spectral shift method is built on), without implementing the full reference-phantom-based compensation.

Shared numerics (pulse generation, scatterer field simulation, attenuation, spectral fitting) live in a small src/qus_refresher package rather than being copy-pasted across notebooks, so a fix in one place doesn't drift out of sync elsewhere. Each notebook is stored twice — a jupytext "percent format" .py as the reviewable source, and an executed .ipynb with real, embedded plot outputs.

There's also a standalone NOMENCLATURE.md — a working reference on RF/QUS terminology, formulas, and trade-offs, and a reading list starting from Oelze & Mamou's 2016 IEEE TUFFC review, through the original Lizzi et al. 1983 JASA paper, to Yao/Zagzebski's work on reference-phantom backscatter coefficient extraction. I also wrote it up as its own post — Ultrasound RF & QUS — Nomenclature and Concepts — if you'd rather read the reference on its own, independent of this repo's code.

Not in scope yet: reference phantom normalization / backscatter coefficient extraction — most of the background papers behind this repo are specifically about that method, and it's the natural next notebook rather than something I decided to skip.

On how this got built

I built this with Claude doing a lot of the actual scaffolding and notebook-writing in a pair-programming style — I described the physics and the pedagogical structure I wanted, reviewed and pushed back on the design (the shared src/ package split, the notebook boundaries, catching that misleading 9.9 dB average above), and every notebook was actually executed end-to-end rather than hand-assembled, so the numbers and plots in the repo are real outputs, not illustrative guesses. It's a genuinely useful way to work through this kind of refresher fast, as long as you're willing to interrogate the numbers it hands you back.

Try it

Repo: qus-python-refresher — public. environment.yml + pip install -e . gets you a working environment; each notebook runs standalone once the shared package is installed. Feedback, issues, and "actually your slope-invariance argument is wrong because..." are all welcome.

Ultrasound RF & QUS — Nomenclature and Concepts

Mine Cuneyitoglu — Mon, 24 Aug 2026 19:55:32 +0000

A reference document for RF-domain ultrasound signal processing, beamforming, and quantitative ultrasound (QUS). Written as a compact, working manual for research engineers moving between the envelope and RF signal domains.

Companion document to qus-python-refresher, a five-notebook Python refresher that implements the RF pulse → echo → moving-window spectral analysis → Lizzi-Feleppa pipeline this document describes. I wrote up the repo itself in a separate post — A Bug That Doesn't Crash Is Worse Than One That Does — if you want the worked code and a genuinely silent bug alongside this reference.

1. Fundamental principles and trade-offs

1.1 Time-frequency uncertainty (Heisenberg-Gabor)

Δ t \cdot Δ f \geq \frac{1}{4 π}

Time and frequency are dual domains; sharpness in one implies uncertainty in the other. This physical constraint underlies most ultrasound design trade-offs:

Short pulse ↔ wide bandwidth. Broadband transducers require short pulses at the cost of transmitted energy and penetration.
Long spectral window ↔ coarse axial localization. FFT bin width is fs / N_window. Shorter windows yield finer axial resolution but coarser spectral resolution.
High center frequency ↔ shallow penetration. Frequency-dependent attenuation limits depth reach at higher f₀.

1.2 Nyquist and aliasing

Sampling rate must exceed 2 × f_max. RF signals typically sampled at 40-60 MHz; IQ (baseband) sampling requires only > BW due to complex sampling. Aliasing manifests as spectral folding and is the most common silent error in prototyping — the code runs, the result is wrong.

1.3 Born approximation (single scattering)

Most QUS models assume single scattering. Multiple scattering regimes (fatty liver, tissues with many interfaces) violate this assumption and complicate clinical interpretation.

1.4 Ergodicity assumption

Regional statistics computed from a single image assume time-averaging equals ensemble-averaging. Rarely validated in practice.

2. Spatial resolution

2.1 Axial resolution

Resolution along the beam axis, set by pulse duration:

R_{a x ia l} = \frac{c \cdot τ}{2}

where τ is pulse duration. Short pulses improve axial resolution. Typical values: ~0.6 mm at 5 MHz with a 2-cycle pulse.

2.2 Lateral resolution

Resolution perpendicular to the beam axis, set by beam width:

R_{l a t} \approx \frac{λ \cdot z}{D} = λ \cdot F_{n u mb er}

where z is focal depth and D is aperture width. F-number is defined as z/D. Small F-number implies wide aperture and improved lateral resolution.

Figure 1: Lateral beam profile narrowing toward the focus and widening beyond it, shown at 10/20/30/40 mm depth, with the −6 dB level marked.

2.3 Elevational resolution

Slice thickness of the imaging plane. In 1D arrays, determined by physical lens design (typically 3-5 mm). 2D matrix probes improve this through electronic focusing.

2.4 Resolution cell (resolution volume)

Intersection of axial, lateral, and elevational resolutions. The number of scatterers within the resolution cell determines the local speckle statistics regime:

Many scatterers (>10) → fully developed speckle, Rayleigh (Nakagami m ≈ 1)
Few scatterers + specular components → sub-Rayleigh, Nakagami m < 1

Figure 2: Ultrasound image axes, resolution cell, and full-width-half-maximum (FWHM). Previously published in Cüneyitoğlu Özkul & Mumcuoğlu, 2018.

2.5 Full Width Half Maximum (FWHM)

Standard measure of beam profile or PSF width, taken at half of maximum amplitude (approximately −6 dB). Used for cross-comparison of resolution across systems.

2.6 Depth of field (DOF)

Length of the focal zone. Scales with F-number squared:

D OF \approx 7.1 λ F_{n u mb er 2}

This approximation is for a linear (1D) array; curvilinear and phased arrays have different beam geometry (diverging aperture, sector scanning) and DOF does not follow the same F-number scaling without modification.

Small F-number improves lateral resolution but reduces DOF — a fundamental trade-off in fixed-focus systems and one motivation for plane-wave imaging.

3. Transmit chain

3.1 Center frequency (f₀)

Peak frequency of the transducer response. Typical clinical range 2-15 MHz. Higher f₀ yields better resolution but reduced penetration; lower f₀ is used for deep applications (cardiac, obstetric).

3.2 Fractional bandwidth

Spectral width around f₀ expressed as a fraction. BW = 0.6 means −6 dB points span 60% of f₀. Broadband transducers (BW > 0.7) benefit two independent things, not one dependent on the other: harmonic imaging (needs enough sensitivity at both f₀ and 2f₀ to transmit and receive on the same element) and reliable QUS spectral fitting (needs enough usable bandwidth for a stable slope/intercept/midband estimate — see section 7.1-7.3). QUS itself is ordinarily performed entirely at the fundamental frequency; it does not require harmonic content.

3.3 Pulse duration

Clinical pulses typically span 1-3 cycles; τ ≈ n_cycles / f₀. Short duration yields wide bandwidth and improved axial resolution.

3.4 Transmit apodization

Weighting of array elements during transmission (Hamming, Hann, Tukey windows) to shape the beam pattern and suppress side lobes.

3.5 Transmit focus

DAS (Delay-and-Sum): Focused to a specific depth via TX delay pattern — see section 4.1 for the full beamforming description.
Plane wave imaging: No transmit focus; the entire field of view is insonified by a plane (or diverging) wavefront. Focusing is deferred to receive-side software beamforming.

4. Beamforming approaches

4.1 Delay-and-Sum (DAS)

Classical approach. Applies element-wise delays and sums the received signals per image line. Focused TX and RX give line-by-line images. Frame rate is limited by the number of scan lines and pulse repetition time.

4.2 Plane wave / Synthetic aperture

TX is unfocused (plane or diverging wave); the entire volume is reconstructed simultaneously via receive-side software beamforming. Frame rates exceed 5000 Hz. Enables ultrafast imaging modalities like shear wave elastography.

Coherent compounding. Combining RF responses from multiple plane-wave angles by complex (pre-envelope) summation, preserving phase. Introduced to clinical practice by SuperSonic Imagine.

4.3 Adaptive beamforming

Data-dependent element weights instead of uniform apodization. Improves contrast and resolution at higher computational cost.

MVDR / Capon: Minimum variance criterion.
DAX: Dual apodization with cross-correlation.

Predominantly research-domain at present.

4.4 Receive apodization

Element weighting on reception to shape the acceptance pattern. Trade-off between resolution and sensitivity to off-axis signals.

4.5 F-number

F# = focal_depth / aperture_width. Standard descriptor of beam narrowness.

5. Signal chain (RF to B-mode)

5.1 RF (radiofrequency)

Raw ADC output. Real-valued, zero-mean, carries the carrier frequency and phase information. Typically sampled at 40-60 MHz. Phase carries diagnostic information that is lost in later stages of the chain.

5.2 IQ / baseband / analytic signal

Complex baseband representation:

Z(t) = I(t) + jQ(t)
Envelope: |Z|
Phase: angle(Z)

Obtained via:

Quadrature demodulation: Multiply by cos(2πf₀t) and sin(2πf₀t), then low-pass filter.
Hilbert transform: Z = RF + j·H{RF}, then frequency shift to baseband.

Complex sampling reduces required rate to > BW, decreasing data volume by a factor of 3-4.

5.3 Envelope detection

Magnitude of the analytic signal. Discards phase, retaining only amplitude modulation. Envelope statistics (Rayleigh, Rician, Nakagami) are defined at this stage.

5.4 Log compression

Logarithmic mapping of envelope for display:

bmode = 20·log₁₀(envelope / max_envelope), clipped to [-60, 0] dB

Compresses the 60+ dB dynamic range for visualization. Must generally be inverted for downstream QUS analysis.

5.5 Scan conversion

Coordinate transformation from acquisition geometry (polar for sector/curvilinear probes) to Cartesian display grid.

5.6 Time Gain Compensation (TGC)

Depth-dependent gain applied at reception to compensate for attenuation. Operator-adjustable in clinical systems. Must be characterized and removed for quantitative measurements.

6. Attenuation and tissue properties

6.1 Attenuation coefficient (α)

Units: dB/cm/MHz. Approximate values in soft tissue:

Fat: 0.6
Liver: 0.5-0.7
Muscle: 1.0-1.5
Bone: >10

Attenuation increases approximately linearly with frequency, producing a depth-dependent spectral downshift that must be compensated in QUS.

6.2 Backscatter coefficient (BSC)

Scattered energy per unit volume per unit solid angle, as a function of frequency. The central quantity in QUS, extracted via reference phantom normalization.

6.3 Speed of sound (c)

Clinical convention: 1540 m/s. Actual values vary:

Fat: ~1450 m/s
Soft tissue: 1540 m/s
Bone: ~4000 m/s

Deviations impact beamforming accuracy and elastography interpretation.

6.4 Scattering regimes (ka dependence)

Ratio of scatterer size (a) to wavelength (λ) governs the scattering behavior:

ka << 1 (Rayleigh): σ_bs ∝ f⁴, isotropic scattering.
ka ~ 1 (Mie / resonance): Oscillatory frequency dependence.
ka >> 1 (geometric): Specular reflection, tissue interfaces.

Most soft-tissue scatterers fall in the Rayleigh regime for clinical ultrasound frequencies.

6.5 Acoustic impedance mismatch

Z = ρ·c. Reflection coefficient at an interface:

R = \frac{Z _{2} - Z _{1}}{Z _{2} + Z _{1}}

Typical values: ~1% at soft tissue interfaces, ~50%+ at soft tissue/bone, ~100% at soft tissue/air.

7. QUS parameters (Lizzi-Feleppa framework)

7.1 Spectral slope (dB/MHz)

Slope of a linear fit to the calibrated log-power spectrum within the usable transducer bandwidth. Sensitive to scatterer size, attenuation, and depth.

7.2 Spectral intercept (dB)

Extrapolated fit value at 0 MHz. Sensitive to effective scatterer concentration.

7.3 Midband fit (dB)

Fit value at f₀. Reflects total backscattered energy in the band.

7.4 Effective scatterer diameter (ESD)

Scatterer size estimate from spectral shape, using Gaussian or spherical form-factor models. Associated with the Insana-Wagner framework.

7.5 Effective scatterer concentration (ESC)

Number of scatterers per unit volume, derived from spectral amplitude.

7.6 Reference phantom method

To isolate tissue backscatter from system response (transducer response, TGC, beamforming), the in vivo spectrum is normalized by that of a phantom with known backscatter properties, imaged with identical acquisition parameters. This is a required step in quantitative QUS.

Key references: Yang & Zagzebski 1990; Yao & Zagzebski 1990.

7.7 Attenuation compensation

Correction for depth-dependent spectral downshift. Common methods:

Spectral shift method
Spectral difference method
Hybrid method (Kim & Varghese 2007)

8. Envelope statistics

8.1 Rayleigh distribution

Envelope of a zero-mean complex Gaussian process. Corresponds to fully developed speckle from many independent scatterers:

p (x) = \frac{x}{σ ^{2}} exp (- \frac{x ^{2}}{2 σ ^{2}})

8.2 Rician (Rice) distribution

Rayleigh with an added coherent component. Characterized by the K-factor:

K = coherent_power / diffuse_power

8.3 K distribution

Envelope model for sub-Rayleigh regimes with few effective scatterers per resolution cell.

8.4 Nakagami distribution

Generalized envelope model. The shape parameter m spans:

m < 1: pre-Rayleigh clustering (sub-Rayleigh)
m = 1: Rayleigh (fully developed)
m > 1: post-Rayleigh (coherent structures)

8.5 Homodyned-K distribution

Two-parameter model capturing coherent, diffuse, and sub-Rayleigh regimes. Preferred in modern QUS envelope statistics.

8.6 Speckle

Coherent interference pattern from unresolved scatterers within the resolution cell. Not stochastic noise: repeat imaging of the same region reproduces the same pattern. Carries diagnostic information exploitable via statistical modeling.

9. Simulation tools

9.1 Field II (Jensen 1996)

Approach: Linear systems using spatial impulse response (Tupholme-Stepanishen).

Covers: Transducer geometry, focused/phased-array beamforming, point scatterer response, B-mode imaging simulation.

Does not cover: Nonlinearity, refraction, reflection at heterogeneous boundaries, wave-medium interaction beyond linear scattering.

Typical runtime: Seconds to minutes.

Use for: B-mode phantom simulation, beamformer development, QUS pipeline testing, teaching.

9.2 k-Wave (Treeby & Cox 2010)

Approach: Full-wave PDE solver on a k-space pseudospectral grid.

Covers: Nonlinearity (Westervelt equation), heterogeneous attenuation, refraction, reflection, complex boundary geometries, photoacoustic sources.

Trade-off: Computationally intensive; overkill for simple point-scatterer B-mode simulation.

Typical runtime: Minutes to hours (accelerated with GPU).

Use for: Photoacoustic reconstruction, HIFU treatment planning, aberration studies, nonlinear harmonic imaging, situations requiring full wave physics.

9.3 Selection guide

Task	Preferred tool
B-mode phantom simulation	Field II
QUS pipeline testing	Field II
Beamformer development	Field II
HIFU / therapeutic ultrasound	k-Wave
Photoacoustic imaging	k-Wave
Aberration correction	k-Wave
Nonlinear harmonic imaging	k-Wave
Rapid prototyping	Custom NumPy/SciPy pipeline

9.4 Python ecosystem

k-wave-python — Official Python port of k-Wave (waltsims/k-wave-python). Provides both a pure NumPy/CuPy solver and an interface to the compiled C++/CUDA k-Wave binaries.
jwave — JAX-based, differentiable acoustic simulation. Suited to machine-learning integrated workflows.
stride — Medical ultrasound tomography framework.
USBMD — PyTorch-based differentiable beamforming library.
PICMUS dataset — Plane-wave imaging benchmark with reference beamformer implementations.

10. Frequently overlooked considerations

10.1 Log-domain estimation bias

Applying log(S + noise) is nonlinear; bias grows as SNR decreases. The order of averaging and log transformation affects the result. Reference: Alam & Lizzi, UMB 2006.

10.2 Focal zone dependence

Beam width varies with depth. Outside the focal zone, lateral resolution degrades, and QUS parameters become depth-biased — see Figure 1 (section 2.2) for the beam-width-vs-depth behavior this is describing.

10.3 Scattering anisotropy

Most QUS models assume isotropic scatterers. Real tissue (muscle fibers, collagen) is anisotropic, so imaging angle influences spectral slope estimates.

10.4 Noise distribution across the chain

RF domain: approximately Gaussian
Envelope: Rayleigh (envelope of Gaussian)
Log domain: approximately Gumbel

Estimator design should reflect the distribution at the analysis stage.

10.5 Window boundary effects

Rectangular windowing produces spectral leakage. Hann or Hamming windows are generally preferred for QUS spectral estimation.

10.6 Log-domain windowing bias

Shorter windows reduce local SNR, amplifying log-domain bias in slope estimation. Both mean and variance of the estimator can shift with window size.

10.7 Silent IQ demodulation errors

Incorrect LPF cutoff causing spectral leakage
Baseband sampling below Nyquist causing aliasing
Hardware amplitude/phase imbalance between I and Q channels

11. Reading list (recommended order)

Oelze & Mamou (2016). "Review of Quantitative Ultrasound: Envelope Statistics and Backscatter Coefficient Imaging and Contributions to Diagnostic Ultrasound." IEEE TUFFC. Bridges envelope statistics and RF spectral analysis; most efficient starting point.
Lizzi, Greenebaum, Feleppa, Elbaum, Coleman (1983). "Theoretical framework for spectrum analysis in ultrasonic tissue characterization." JASA 73:1366-1373. Foundational reference for the Lizzi-Feleppa method.
Yao, Zagzebski, Madsen (1990). "Backscatter coefficient measurements using a reference phantom to extract depth-dependent instrumentation factors." Ultrasonic Imaging. Practical protocol for reference phantom normalization.
Alam & Lizzi (2006). "On the statistics of ultrasonic spectral parameters." UMB. Log-domain bias analysis; relevant to Bayesian estimation workflows.
Bercoff, Tanter, Fink (2004). "Supersonic shear imaging: a new technique for soft tissue elasticity mapping." IEEE TUFFC. Foundational for plane wave imaging with coherent compounding.
Mamou & Oelze (eds.) (2013). Quantitative Ultrasound in Soft Tissues. Springer. Comprehensive treatment; chapters 1-2 are the most efficient overview.

This started as a personal working reference, updated as new concepts and experimental observations get incorporated. The runnable side of it — pulse simulation, the moving-window spectral analysis, the Lizzi-Feleppa fits — lives in qus-python-refresher, written up separately here.

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

Mine Cuneyitoglu — Thu, 30 Jul 2026 16:44:09 +0000

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:

m^= \frac{( E [ X ^{2} ] ) ^{2}}{Var ( X ^{2} )}, ω^= E [X^{2}]

works exactly when the samples

X_{1}, X_{2}, \dots, X_{N}

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

1/ 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

m^

is systematically underestimated, and the mean of

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} (τ) = R_{x} (τ) * R_{h} (τ) + R_{n} (τ)

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

y = A x + n, A = D B

where

y

is the observed image (in vectorized form),

x

is the underlying high-resolution reflectivity,

B

is the blur matrix representing the PSF,

D

is a decimation matrix, and

n

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

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:

x^= ar g x max p (y ∣ x) p (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:

(y - A x)^{T} C^{- 1} (y - A x) + λ Φ (x)

where

Φ (x)

is a prior penalizing local dissimilarities in the reflectivity, and

λ

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} 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.

US Speckle Models — Part 3: m as a Tissue Map

Mine Cuneyitoglu — Sat, 25 Jul 2026 06:46:01 +0000

Two posts ago I argued that ultrasound speckle isn't random noise — it's the interference pattern of many sub-resolution scatterers, and its statistics carry information about the tissue that produced it. Rayleigh distribution was the classical model: mathematically elegant, physically motivated, but built on assumptions that break in real tissue.

Nakagami generalized Rayleigh. Same physical intuition, one extra parameter, and the freedom to adapt to whatever scatterer density the tissue actually has.

Last post I laid out the Nakagami PDF and what each term does. This post is where it stops being a math exercise and starts being clinically useful.

From distribution to map

The Nakagami PDF is:

p (x; m, ω) = \frac{2 m ^{m}}{Γ ( m ) ω ^{m}} x^{2 m - 1} e^{- m x^{2} / ω}, x \geq 0

Two parameters carry the tissue information:

m — shape parameter. Low m spreads the distribution out (few scatterers, high uncertainty in the envelope amplitude). High m sharpens the peak (many scatterers, predictable behavior).
ω — scale parameter, equal to $E [x^{2}]$ . Sets the local mean intensity, normalizes for depth attenuation and gain.

At m = 1, the whole expression collapses to Rayleigh. Nakagami includes Rayleigh as a single point in a continuum, then extends both directions — sub-Rayleigh (fewer scatterers, edges and boundaries) and super-Rayleigh (denser, more organized tissue).

Here's the pivot that makes it clinically useful: m can be estimated locally. Slide a small window across the image, compute the method-of-moments estimator inside it:

m^= \frac{( E [ X ^{2} ] ) ^{2}}{Var ( X ^{2} )}, ω^= E [X^{2}]

The output isn't a denoised image. It's a map — one m value per window, spatially varying, tracking whatever tissue characteristic changes as you move across the scan.

Canonical envelope-domain values

Published literature converges on a rough three-regime interpretation, all measured on raw envelope amplitude before any display processing:

m ≈ 0.5 — lesion boundaries, cyst walls, tissue interfaces. Few scatterers per resolution cell, high variance in the envelope.
m ≈ 1.0 — fully developed speckle. What Rayleigh assumed everywhere. Homogeneous soft tissue at sufficient depth.
m > 1.5 — dense parenchyma, fibrous tissue, structured collagen. Many scatterers, predictable statistics.

Same image, three regimes, one model. That's the promise.

From the thesis: real B-mode, real fit

Here's what happens when you actually run this on a real hospital scan.

Three 64×64 windows on a B-mode frame from my PhD data:

A — inside a dense, near-uniform lesion. High m: histogram is narrow, PDF nearly Gaussian-looking.
B — a homogeneous parenchyma region. Mid-range m, classic developed-speckle shape.
C — a heterogeneous zone with varying intensity. Low m, right-skewed distribution.

The relative ordering matches the canonical interpretation: dense → high m, heterogeneous → low m. But the absolute m values are much higher than the 0.5–2 range the literature reports. The reason is what the image is: this is not envelope amplitude. This is the display frame — log-compressed, gain-adjusted, gamma-corrected, 8-bit quantized.

Log compression squashes the dynamic range. A distribution that looked broad and skewed in the envelope domain becomes narrower and more symmetric after compression. Method-of-moments picks that up as a higher m estimate. It's not wrong — it's just measuring m on a different variable than the textbook does.

For most practical work, the relative m ordering carries the tissue information. Absolute values shift with the specific compression curve of the scanner, but the ranking of regions doesn't.

I'll come back to why this matters at the end.

For contrast: what the textbook world looks like

If I simulate the same three regimes in the envelope domain — Nakagami-sampled directly, no compression, on a synthetic sector scan geometry — this is what I get:

Two things worth noting from the contrast:

First, the m values here (0.5, 1.0, 1.7) sit exactly where the canonical literature puts them. That's because this is synthetic envelope amplitude, drawn directly from Nakagami distributions with those m values. No log compression, no display processing in between. The values are the ground truth, not an estimate contaminated by a display pipeline.

Second, the PDF shapes are visibly cleaner — sharper mode separation, less overlap between regimes. Real clinical data doesn't look like this. Once you go through display processing, distributions compress toward each other, m estimates shift, and the model has to accommodate what the scanner actually outputs, not what a textbook assumes.

Where a lot of QUS quietly breaks

This is where a lot of published quantitative ultrasound work quietly breaks when it hits real hospital data.

Most method development happens on simulated envelope data or on lab-grade research scanners that expose raw RF or beamformed envelope. That's the world where the canonical m values live. It's a clean world, and the math works there.

Real clinics don't run research scanners. They run commercial systems that output DICOM images with log compression, dynamic range mapping, and a proprietary post-processing pipeline you don't have access to. If you take a method validated on envelope-domain simulation and drop it on a hospital scanner's video output, three things happen:

Absolute m values shift out of the calibration range the method assumes
Distributions overlap more than the training data, so region classification degrades
The relationship between m and underlying tissue reflectivity is no longer direct — there's a monotonic-but-nonlinear display curve sitting in between

The math is right. The domain assumption isn't.

You can still do useful QUS on log-compressed data. But you have to design for it — pick estimators that are compression-tolerant, work with relative rather than absolute m, and calibrate on the specific scanner's output curve. That's a different research program than "assume envelope, apply Nakagami."

Next: the other thing Nakagami by itself misses

There's a second problem, independent of log compression. Nakagami models the amplitude distribution at each pixel. It says nothing about the relationship between pixels. And in ultrasound, that relationship is never zero: the PSF spreads every scatterer's signal across a neighborhood, and neighboring pixels are always weighted combinations of overlapping scatterer sets.

Standard local m estimation treats those pixels as if they were independent samples. They aren't. Next post: why the iid assumption breaks, what it costs you, and what you have to fold into the observation model to fix it.

References

My paper on Bayesian restoration for B-mode ultrasound, accounting for correlated speckle noise: Ultrasonic Imaging, 2019
ORCID (all publications): 0000-0003-2633-9371

I take on selective remote consulting on ultrasound image quality, quantitative ultrasound, and restoration problems — particularly the messy real-clinical-data variety this post is about. Reach out via LinkedIn if that's the kind of problem you're working on.

US Speckle Models — Part 2: The Nakagami Distribution

Mine Cuneyitoglu — Thu, 16 Jul 2026 17:53:20 +0000

Part 1 argued that Rayleigh is the right model for one specific regime — many independent scatterers, homogeneous tissue, no interest in absolute quantitative recovery — and breaks in three important ways when tissue doesn't satisfy those conditions. This part is about the generalization that fixes the first two of those problems (scatterer count assumptions and the fixed-shape constraint) while explicitly not fixing the third (pixel correlation). That last one waits until Part 4.

The generalization is the Nakagami distribution.

Where it comes from

The Nakagami distribution wasn't invented for ultrasound. Minoru Nakagami developed it in the 1940s to model the amplitude of received radio signals in HF propagation, where the received envelope showed variability the Rayleigh model couldn't capture — sometimes more, sometimes less. He wanted one flexible family that could sit anywhere along the spread-vs-concentration axis, with a single shape parameter to move between regimes.

Medical ultrasound picked it up several decades later for the same reason: the envelope of the backscattered signal doesn't always look Rayleigh, and a single flexible family beats switching between separate models for boundaries, parenchyma, and fibrous tissue.

The PDF:

p (x; m, ω) = \frac{2 m ^{m}}{Γ ( m ) ω ^{m}} x^{2 m - 1} e^{- m x^{2} / ω}, x \geq 0

Two parameters, one normalizing constant, one exponential kernel. Each piece is doing specific work.

Reading the terms

x — envelope amplitude. The variable being modeled. This is what the scanner reconstructs from the beamformed complex signal and hands off as a pixel intensity, up to whatever post-processing sits between envelope detection and display. In the raw envelope domain, x is real and non-negative. In the display domain, x is the pixel value after log compression, gain, and gamma — related to the true envelope amplitude but not identical to it. Part 3 gets into what that difference costs.

m — shape parameter. This is the term that made Nakagami interesting for tissue characterization. It's dimensionless, non-negative, and controls how concentrated the distribution is around its mode:

Low m (m < 1) → distribution spreads out, heavy tail, mode close to zero. Physically: few scatterers per resolution cell, high uncertainty in the envelope amplitude. Boundaries, interfaces, sparsely-scattering regions.
m = 1 → Rayleigh case (details below). Fully developed speckle regime.
High m (m > 1) → sharp peak, tight distribution, mode moves away from zero. Physically: many scatterers, more organized structure, predictable statistics.

m is not fixed by physics — it's estimated from the data. Give it a homogeneous parenchyma window and it comes back near 1. Give it a boundary or cyst wall and it comes back below 1. Give it dense fibrous tissue and it comes back above 1. The distribution adapts to what's actually there.

ω — scale parameter. Equal to

E [x^{2}]

, the mean squared envelope amplitude. It sets the intensity level of the distribution — think of it as absorbing depth-dependent attenuation, gain settings, and any global scaling of the image. When you compare m between regions, ω is what you normalize by to strip out those global effects.

Γ(m) — the gamma function. The normalizing constant. It ensures the PDF integrates to 1 across all valid x. For integer m,

Γ (m) = (m - 1)!

; for non-integer m it interpolates smoothly. Its role here is structural — nothing physical, just what makes the math a proper probability distribution.

exp(−mx²/ω) — the exponential kernel. This is where the shape work happens. As m increases, this term penalizes large deviations from the mean more aggressively, which is what makes high-m distributions look tightly concentrated. As m decreases, the penalty softens, and the distribution spreads. All the intuition about "few scatterers → high uncertainty → wide distribution" and "many scatterers → predictable → tight distribution" is carried by this single term.

The prefactor

x^{2 m - 1}

shapes the low-x behavior — pulls the distribution away from zero as m grows, allows density at zero when m is small.

Rayleigh as a special case

The claim that Nakagami reduces to Rayleigh at m = 1 is often stated but rarely shown. It's worth doing explicitly, because the substitution reveals how

ω

and Rayleigh's

σ

relate.

Set m = 1 in the Nakagami PDF. Since

Γ (1) = 1

p (x; 1, ω) = \frac{2}{ω} x exp! (- \frac{x ^{2}}{ω})

Compare to the Rayleigh PDF from Part 1:

p (x; σ) = \frac{x}{σ ^{2}} exp! (- \frac{x ^{2}}{2 σ ^{2}})

These match if

ω = 2 σ^{2}

. So Nakagami's scale parameter is exactly twice the variance of the Gaussian I/Q components — the same physical quantity, written in a different convention.

The offset issue from Part 1 carries through: at m = 1, the Nakagami mean is still

σ π /2

in the underlying units, and the envelope-to-scatterer-strength mapping still has that

π /2

factor sitting in the middle. Nakagami doesn't fix that — it wasn't designed to. What Nakagami fixes is the shape rigidity.

What Nakagami buys you, and what it doesn't

The right way to describe Nakagami is that it has no built-in assumptions about scatterer count. Both parameters, m and ω, are free — estimated from the local data — and the model adapts to whatever regime the tissue is in. This is different from a model that has a hand-tuned hyperparameter you have to pick a priori.

The standard estimator is method of moments:

m^= \frac{( E [ x ^{2} ] ) ^{2}}{Var ( x ^{2} )}, ω^= E [x^{2}]

Both quantities come directly from the sample. There's a maximum-likelihood variant that's slightly more efficient in small windows, but for most practical purposes moment-based estimation is what the literature uses and what tools implement.

What Nakagami does not buy you:

Pixel independence. The estimator above assumes samples inside the window are i.i.d. In ultrasound they aren't — the PSF spreads every scatterer's contribution across a neighborhood, and pixels are correlated by construction. Local m estimates carry a PSF-shaped bias that grows with the ratio of PSF width to window size. Part 4 is about this.
A separation between distribution and display. Nakagami models the envelope amplitude. If your data has gone through log compression, gain adjustment, and gamma correction — which is what every commercial hospital scanner does before writing DICOM — the pixel values aren't envelope amplitudes anymore. You can still fit Nakagami, but the m you get is a pseudo-m on the display domain, not the true envelope-domain m. Part 3 is about this.
A guarantee that m corresponds to a specific tissue property. m tracks the scatterer density in a statistical sense, but the mapping isn't one-to-one with any single clinical variable. Interpreting m maps as diagnostic requires calibration against ground truth for the specific tissue and scanner combination.

Below are the visuals that accompanied the original discussion — Nakagami PDFs at several m values, and a formulation diagram tying the model to the physical scattering setup.

Nakagami probability density function shown for a range of shape parameter m values.

Physical formulation tying the Nakagami model to the underlying scattering setup.

Nakagami models the amplitude distribution at each pixel. What it doesn't model is the relationship between pixels — and that relationship is never zero in ultrasound, because the PSF couples them.

Before we get to correlation, though, there's a step in between: what happens when you estimate m locally across the image, sliding a window and producing a spatial map of the shape parameter. That's Part 3 — where Nakagami stops being a noise model and becomes a tissue characterization tool.

References

My paper on Bayesian restoration for B-mode ultrasound, using a system model that folds in both Nakagami-like envelope statistics and correlated noise: Ultrasonic Imaging, 2019
ORCID (all publications): 0000-0003-2633-9371
Previous in series: Part 1 — Rayleigh Distribution

I take on selective remote consulting on ultrasound image quality, quantitative ultrasound, and restoration problems — particularly on data from commercial hospital scanners rather than research-grade RF. Reach out via LinkedIn if that's the kind of problem you're working on.

US Speckle Models — Part 1: Rayleigh Distribution

Mine Cuneyitoglu — Fri, 10 Jul 2026 18:27:38 +0000

Ultrasound speckle isn't random noise. It's physics.

The bright-and-dark grainy texture you see in every B-mode image is not a sensor artifact and not electronic noise. It's the deterministic interference pattern of many sub-resolution acoustic scatterers, imaged through a finite-aperture point spread function. If you scanned the exact same tissue twice with the exact same probe position and settings, you'd get the exact same speckle. It's a signal, not a nuisance.

The reason speckle looks random is that the underlying scatterer configuration is unknown to us — we can't resolve individual scatterers, and small changes in probe position redistribute the interference pattern in ways we can't track. So we treat it statistically, and the first statistical model anyone reaches for is Rayleigh.

This is Part 1 of a series on statistical models for ultrasound speckle. Here I'll cover where the Rayleigh model comes from, why it's a natural starting point, and the three specific places it stops holding up in real tissue.

Where Rayleigh comes from

Consider a single resolution cell containing many independent scatterers. Each scatterer contributes a complex-valued backscatter with some amplitude and phase. The total received signal is the coherent sum of all these contributions.

If we split that sum into in-phase (I) and quadrature (Q) components:

S = I + j Q = k = 1 \sum N a_{k} e^{j ϕ_{k}}

where the scatterers have random amplitudes and random phases uniformly distributed on

[0, 2 π]

, and N is large.

By the Central Limit Theorem, both I and Q converge to zero-mean Gaussian distributions with equal variance

σ^{2}

. The envelope amplitude — what the B-mode reconstruction actually visualizes — is the magnitude:

A = I^{2} + Q^{2}

The magnitude of two independent zero-mean Gaussians with equal variance follows a Rayleigh distribution:

p (A; σ) = \frac{A}{σ ^{2}} exp! (- \frac{A ^{2}}{2 σ ^{2}}), A \geq 0

One parameter,

σ

, related to the underlying scatterer population. Mathematically elegant, physically motivated, and analytically tractable. This is why the model has stayed in textbooks for decades.

Where it breaks

Real tissue satisfies the CLT assumption in some places, and violates it in others. Three specific ways this shows up in practice:

1. Scatterer count isn't always "many"

The CLT argument assumes N is large enough for the Gaussian limit to hold — in practice, on the order of 10 or more scatterers per resolution cell. In homogeneous parenchyma at sufficient depth, this is a reasonable assumption. At lesion boundaries, cyst walls, tissue interfaces, and structured collagen, it isn't.

When N is small, the sum in the equation above doesn't converge to a Gaussian, and the envelope doesn't converge to Rayleigh. What you get instead is a sub-Rayleigh distribution — broader, with heavier tails and more variance in the low-amplitude region. The Rician distribution captures the case of few strong scatterers plus diffuse background; the Nakagami distribution generalizes further and covers both sub- and super-Rayleigh cases with a single shape parameter.

This isn't a small correction. Regions where Rayleigh breaks are exactly the regions clinicians care about — boundaries, lesions, structural changes — and forcing a Rayleigh fit there systematically misrepresents the tissue statistics.

2. Envelope detection has a built-in offset

The mean of a Rayleigh distribution is:

E [A] = σ \frac{π}{2} \approx 1.253 σ

This creates a subtle but important issue for quantitative work. If you want to recover the underlying scatterer strength parameter

σ

from an image, you can't just take the local pixel mean and use it directly — the envelope detection operation inflates the mean by a factor of

π /2

relative to the parameter that actually characterizes the tissue.

For qualitative B-mode viewing this doesn't matter — the image still looks like tissue. For quantitative analysis, it matters a lot. Any metric that treats pixel intensity as a direct proxy for scatterer density carries this systematic offset. It's not measurement error, it's not noise — it's a consequence of how the envelope operator maps the underlying complex signal to a positive real amplitude.

You can compensate for it if you know Rayleigh holds locally. You can't compensate for it cleanly when Rayleigh doesn't hold — because the relationship between pixel mean and scatterer strength changes with the underlying distribution.

3. Pixels aren't independent samples

Every scatterer's contribution is spread across a small spatial neighborhood by the point spread function. Neighboring pixels aren't independent samples of the same distribution — they're weighted combinations of overlapping scatterer sets. Even if the underlying distribution were perfectly Rayleigh at every location, the observed pixels would be correlated.

Standard estimators for

σ

— or for any distribution parameter, Rayleigh or otherwise — assume i.i.d. samples. When pixels are correlated, those estimators become biased, and their variance is systematically underestimated. Confidence intervals shrink, false-discovery rates rise, and the resulting parameter maps carry the PSF footprint mixed in with the tissue signal.

Most published QUS work on Rayleigh-based analysis quietly assumes independence and hopes the correlation is small enough not to matter. Sometimes it is. Often it isn't. This is a topic in its own right, and I'll come back to it in a later post in the series.

Where the series goes from here

Rayleigh isn't wrong. It's the right model for a specific regime — many independent scatterers, homogeneous tissue, no interest in absolute quantitative recovery of scatterer strength. That regime exists, and Rayleigh describes it accurately. What it isn't is a general model for ultrasound speckle.

Later parts of the series work through what to reach for when the Rayleigh assumptions break:

Part 2 — Nakagami as a generalization, its shape parameter, and what each term in the PDF actually does
Part 3 — Local m estimation and the shift from noise modeling to tissue characterization
Part 4 — Correlated speckle, the PSF signature in the autocorrelation function, and what has to go into the observation model to handle it

Below are two visuals from the original discussion — Rayleigh compared to Nakagami PDFs at a few shape parameters, and a stylized picture of the "fully developed speckle" regime the CLT argument assumes.

References

My paper on Bayesian restoration for B-mode ultrasound, including the correlated speckle noise model previewed in point (3) above: Ultrasonic Imaging, 2019
ORCID (all publications): 0000-0003-2633-9371

US Speckle Reduction Part 2 — Classical vs DnCNN

Mine Cuneyitoglu — Fri, 26 Jun 2026 15:49:01 +0000

Speckle Reduction on Ultrasound Images — Classical vs DnCNN

Background

My PhD research focused on freehand B-mode ultrasound image quality —
speckle noise, spatial correlation modeling, and super-resolution restoration.
That was 8 years ago. This project is me rebuilding that work with modern tools.

Starting point: a clean benchmark on the
BUSI dataset
(780 breast ultrasound images — benign, malignant, normal).

The Problem with Standard Benchmarks

Most speckle reduction papers evaluate on the full image — background included.
But in clinical practice, what matters is image quality within the lesion.

So I added ROI-based evaluation using the segmentation masks provided in BUSI.
For images with multiple lesions, masks were combined via union.

def evaluate_with_roi(original, denoised, mask=None):
    global_psnr = psnr(original, denoised, data_range=1.0)
    global_ssim = ssim(original, denoised, data_range=1.0)

    if mask is not None and mask.max() > 0:
        # ROI bounding box
        rows = np.any(mask > 0, axis=1)
        cols = np.any(mask > 0, axis=0)
        rmin, rmax = np.where(rows)[0][[0, -1]]
        cmin, cmax = np.where(cols)[0][[0, -1]]

        if (rmax - rmin) >= 7 and (cmax - cmin) >= 7:
            roi_orig = original[rmin:rmax, cmin:cmax]
            roi_den  = denoised[rmin:rmax, cmin:cmax]
            roi_ssim = ssim(roi_orig, roi_den, data_range=1.0)

Methods

Noise model: Additive Gaussian (σ=0.05) as speckle proxy (Note: not physically accurate — Nakagami or correlated models are more realistic. PSF-aware correlated noise model is next on the list.)
Non-Local Means (NLM)
Median Filter
Gaussian Filter
DnCNN — 17-layer residual CNN, trained from scratch on BUSI

DnCNN — A Note on Training

Training was not straightforward.

Initial run with lr=1e-3 plateaued immediately after epoch 2:
Epoch [ 2/30] — Loss: 0.002472

Epoch [10/30] — Loss: 0.002408

Epoch [20/30] — Loss: 0.002383

Epoch [30/30] — Loss: 0.002363 ← barely moved

Dropping to lr=1e-4 with step decay fixed it:
Epoch [ 5/50] — Loss: 0.001487

Epoch [20/50] — Loss: 0.001023

Epoch [50/50] — Loss: 0.000899

This is worth noting: DnCNN is sensitive to initialization and learning rate.
Results varied across runs (PSNR range ~31.2–31.7 dB). Single-run DL benchmarks
should be interpreted with some caution.

Results (780 images)

Method	PSNR Global	SSIM Global	PSNR ROI	SSIM ROI
Noisy	26.25	0.649	26.34	0.681
Median Filter	29.71	0.810	29.32	0.811
Gaussian Filter	29.99	0.836	29.76	0.833
NLM	31.47	0.827	31.52	0.842
DnCNN	31.24	0.847	31.17	0.860

Key Findings

1. NLM vs DnCNN — it depends on the metric.
NLM wins on PSNR, DnCNN wins on SSIM. SSIM captures structural similarity
better than pixel-level fidelity — arguably more relevant for diagnostic imaging.

2. Global metrics can mislead.
Median filter looks decent globally (29.71 dB) but drops in ROI (29.32 dB) —
edge blurring at lesion boundaries. A paper reporting only global metrics
would miss this.

3. Dataset-scale evaluation matters.
On individual images, NLM sometimes outperforms DnCNN. The pattern only
stabilizes at dataset scale. Single-image comparisons in this field are unreliable.

What's Next

Super-resolution: Real-ESRGAN / SwinIR fine-tuned on ultrasound
Realistic noise model: PSF-aware correlated speckle from phantom data (I actually collected this data during my PhD — time to revisit it)
Thyroid nodule segmentation with SAM

Code on GitHub 👇

https://github.com/minebyte/us-speckle-reduction-benchmark

Classical Speckle Reduction Benchmark

Mine Cuneyitoglu — Fri, 26 Jun 2026 12:42:06 +0000

Classical speckle reduction benchmark on 780 breast ultrasound images (BUSI dataset).

Compared Non-Local Means, Median, and Gaussian filtering —
evaluated both globally and within lesion ROI using segmentation masks.

Key finding: global metrics can be misleading :) (why proper masking is important)

Median filter scores well globally (PSNR: 29.71 dB) but drops notably
within lesion boundaries (29.32 dB ROI) — clinically relevant edge blurring
that global evaluation misses.

NLM remains consistent across both evaluations — more reliable for diagnostic regions.

Next: DnCNN and U-Net comparison.

GitHub 👇
https://github.com/minebyte/us-speckle-reduction-benchmark

MedicalImaging #Ultrasound #ImageProcessing #ComputerVision #Python