Part V · Module 31

ULM end to end (capstone)

Ultrasound localization microscopy turns a diffraction-limited scanner into a microscope by treating every isolated microbubble as a point source, localizing it far below the wavelength and accumulating hundreds of thousands of positions into a vascular map. Every stage is an estimator from an earlier module: the clutter filter is an SVD (Module 06), motion correction is delay estimation (Module 11), detection is matched filtering (Module 10), localization is bounded by the CRLB (Module 12), the bubbles are those of Module 28 and the beams those of Module 23. This capstone runs the whole chain on simulated data, quantifies what each stage loses, and ends with the metrics (saturation, FRC) and the frontiers (high concentration, 3D, transcranial, functional) where the field is moving.

You will be able to

Design a ULM acquisition from first numbers: frame rate, angles, displacement per frame, linking distance, block length for the SVD filter, concentration and acquisition time
Choose SVD cut-offs from the singular spectrum and the Doppler content of the singular vectors, and predict which slow bubbles the filter removes
Set a detection threshold from a false-alarm rate, localize with weighted average, Gaussian fit, radial symmetry or interpolation, and compare each with the CRLB \(\sigma_x = w_x/\sqrt{\mathrm{ENR}}\)
Link localizations into tracks with nearest neighbour, greedy or Hungarian assignment and a Kalman predictor, and predict the linking error from the density and displacement
Compute saturation curves and acquisition times per vessel size, measure resolution with Fourier ring correlation without fooling yourself, and separate localization precision from resolution
Place the frontiers on the same map: PALA benchmarking, learned localization at high concentration, 3D ULM with matrix and row–column arrays, transcranial ULM with aberration correction, and functional ULM

Assumes

Module 06 · Eigendecomposition & SVD: the SVD clutter filter; Module 11 · Correlation & time-delay estimation: subsample delay fits, speckle tracking, the Kasai autocorrelator
Module 12 · Estimation theory and the Cramér–Rao bound: Section 12.5 on ULM localization; Module 10 · Matched filtering & pulse compression: detection in white noise; Module 23 · Beam patterns, k-space and compounding: plane-wave compounding
Module 28 · Microbubble physics, contrast imaging and the ULM signal source: Section 28.8, the density trade-off; Module 29 · Doppler and flow imaging: SVD in ultrafast Doppler; Module 19 · Array architectures: matrix and row–column arrays; Module 22 · Acoustic output, bioeffects and regulatory measurement

Exact formulaclosed form plotted as is, or published numbers re-plottedSimulationsolvers that run here in JS: IQ synthesis, SVD, detection, localization, assignment, FFT-based FRC, Monte CarloApprox. modelsimplified physics; assumptions under the badgeSchematicdrawing for intuition, not to scaleDefaults: \(f_0\) = 5 MHz, \(c\) = 1540 m/s, \(\lambda\) = 308 µm; Gaussian PSF with lateral FWHM \(\lambda\) and axial FWHM \(0.75\lambda\); beamformed pixel \(\lambda/4\) = 77 µm; 1000 compounded frames per second.

31.1The pipeline as a chain of estimators, and the acquisition budget

A ULM image is the histogram of a very long list of estimated positions. Each entry of that list survived seven decisions: the bubble was in the image plane during an acquisition, its echo was separated from tissue, it was detected, it was isolated enough to be localized, its position was estimated, it was linked to its neighbours in time, and the resulting track was long enough to be kept. Each decision has a rate (what fraction passes) and an error (what the survivors carry with them), and the final image inherits both. Module 28 drew this chain as a schematic (Figure 28.E); here it runs.

The acquisition sets the numbers that every later stage lives with. Plane or diverging waves at \(N_a\) angles are compounded into one frame (Module 23); the pulse repetition frequency is capped by the round trip to the deepest echo; the frame rate fixes how far a bubble moves between frames, and that displacement, compared with the distance between bubbles, decides whether tracking works:

\[ f_{\mathrm{PRF}} \le \frac{c}{2z_{\max}}, \qquad \mathrm{FR} = \frac{f_{\mathrm{PRF}}}{N_a}, \qquad \delta = \frac{v}{\mathrm{FR}}, \qquad v_{\max} = d_{\mathrm{link}}\,\mathrm{FR}, \qquad \mathrm{ENR} \propto N_a \]

(31.1)

\(f_{\mathrm{PRF}}\): transmit (firing) rate, Hz
\(z_{\max}\): deepest echo that must return before the next firing, m
\(N_a\): number of steered transmits compounded into one frame
\(\mathrm{FR}\): compounded frame rate, frames/s
\(v,\ \delta\): bubble speed (m/s) and displacement between frames (m)
\(d_{\mathrm{link}}\): largest displacement the tracker accepts between frames, m
\(\mathrm{ENR}\): energy-to-noise ratio of one bubble in one frame (Module 12); coherent compounding adds signal in amplitude and noise in power
\(c\): sound speed, 1540 m/s
\(\propto\): proportional to

At 30 mm depth \(f_{\mathrm{PRF}} \le\) 25.7 kHz, so 5 angles allow at most 5.1 kHz; most systems run near 1 kHz, limited by data transfer and beamforming. At 1 kHz a 1 mm/s capillary bubble moves 1 µm per frame and a 30 mm/s arteriole bubble 30 µm (λ/10). PALA's default gate of 2 pixels of λ at 15.6 MHz (about 200 µm) tracks up to about 200 mm/s at 1 kHz. Five angles raise the ENR by 7 dB and shrink the localization error by \(\sqrt5\).

Two second-order effects matter for fast bubbles. A bubble moving axially at 30 mm/s during the 0.8 ms that 5 angles take at 5 kHz advances by 24 µm, a phase of \(4\pi\delta/\lambda\) = 0.98 rad across the compounding window; the coherent sum loses only 0.5 dB, but the loss grows quickly past π. And the same phase rotation from frame to frame is exactly the Doppler signal (Module 29) that lets the SVD separate bubbles from tissue in Section 31.2: speed is both what the clutter filter needs and what the tracker fears.

The table maps each stage onto the module that built its tools, the number that sets its error, and where this capstone treats it.

Stage	Estimator	Sets the error	Typical choice	Tools from	Here
1 Acquire	compounded plane or diverging waves	ENR per frame, displacement per frame	3–10 angles, 0.5–1 kHz, minutes	19, 22, 23, 28	31.1
2 Clutter filter	SVD of the Casorati matrix, rank cut-offs	slowest bubble kept, residual clutter	blocks of 100–800 frames; drop the first few components	06, 29	31.2
3 Motion correction	rigid or affine registration of tissue frames	residual motion added to every position	speckle correlation per block, applied to localizations	11	31.2
4 Detect	local maxima of a matched filter or NCC with the PSF	false alarms, missed and overlapping bubbles	threshold from a false-alarm rate, or the N brightest maxima	10	31.3
5 Localize	weighted average, Gaussian fit, radial symmetry, interpolation	\(w/\sqrt{\mathrm{ENR}}\) plus estimator bias	a small window around the peak, about one FWHM	11, 12	31.4
6 Track	nearest neighbour, Hungarian, Kalman	wrong links when bubbles are closer than their displacement	gate of 2–4 pixels, minimum track length 5–20 frames	new here	31.5
7 Render and measure	density, velocity and direction maps; saturation; FRC	sampling (how many bubbles visited each pixel)	λ/10 to λ/20 grid, track interpolation	02, 04	31.6

Demo 31.1

End-to-end ULM on a simulated microvessel tree: SVD → detect → localize → track → render

SimulationEverything runs in JS on synthetic beamformed IQ: 64 × 64 pixels of λ/4 (4.9 × 4.9 mm), 1000 frames/s, blocks of 100 frames. Bubbles flow along six vessels (120, 60, 60, 30, 30, 20 µm; centreline speeds 25 to 1.5 mm/s; parabolic profile; bubble count per vessel proportional to diameter × length, a compromise for the demo since blood volume would scale as diameter²), each imaged as a Gaussian PSF (FWHM λ × 0.75λ) with log-normal amplitude (4 dB spread), random phase and the Doppler phase \(-4\pi z/\lambda\). Clutter: three static speckle fields with slow complex gains (0, −12, −24 dB; the first carries a 6 µm axial pulsation at 3 Hz, the others rotate at about 2 and 5 Hz, so the clutter has rank of about 3–4), 30 dB above the median bubble. Noise: white complex Gaussian per pixel (not band-limited). SVD: top-24 subspace iteration with Rayleigh–Ritz on the 100 × 100 frame Gram matrix of each block, then projection; this is the truncated-SVD filter, computed only for the components it needs. Detection: matched filter with the PSF envelope, local maxima, threshold from the false-alarm rate with the noise power estimated from the median. Assignment: own nearest-neighbour and Hungarian code.

detection / localization / kept tracktrue bubble positionrejected short track, true vessel outlinesingular values removed

SVD cut-off (components removed) 3

False-alarm rate per pixel 10⁻⁵

Localizer

Tracker

Linking gate d_link 0.25 λ

Minimum track length 10 frames

Median bubble ENR (re-acquires) 24 dB

Bubbles in the field (re-acquires) 20

Frame shown (last block) 50

Render

Press a control to run.

Try: set the SVD cut-off to 0: the clutter is 30 dB above the bubbles, the median-based noise estimate becomes the clutter level and the detector finds nothing. At 3 the spectrum shows why: three tissue components near 0–5 Hz, then components whose mean Doppler is tens of hertz (bubbles). Push the cut-off to 12: recall halves, false alarms rise (removing bubble components leaves negative ghosts of slow bubbles) and the small branches vanish from the render. The localization error sits well above the CRLB of the median bubble: the 3-point fit on λ/4 pixels costs a factor 2–3 (Section 12.5, Demo 31.3), log-normal amplitudes make half the bubbles weaker than the median, and bubbles closer than a PSF along the same vessel bias each other. Raise the count to 60: recall falls and wrong links rise, but nearest neighbour and Hungarian agree here because each bubble moves at most λ/12 per frame; the wrong links come from overlapping bubbles, not from assignment (Demo 31.4 isolates the regime where assignment matters). Press Acquire: saturation grows, slowly for the small vessels.

31.2Clutter filtering and motion correction: what the SVD keeps, what motion adds

Tissue echoes are 30–60 dB stronger than a bubble and almost perfectly coherent from frame to frame; bubbles are weak, sparse and moving. Module 06 built the filter that exploits the difference: stack a block of \(N_t\) frames as the columns of the Casorati matrix \(\mathbf{S}\) (pixels × frames), take its SVD and discard the first \(r\) components, which carry the energetic, slowly varying tissue. Module 29 showed the same filter in ultrafast Doppler. In ULM it decides which bubbles can ever be localized:

\[ \mathbf{S} = \sum_{k=1}^{N_t} \sigma_k\,\mathbf{u}_k\mathbf{v}_k^{H}, \qquad \mathbf{S}_{\mathrm{b}} = \mathbf{S} - \sum_{k=1}^{r} \sigma_k\,\mathbf{u}_k\mathbf{v}_k^{H} = \mathbf{S}\left(\mathbf{I} - \mathbf{V}_r\mathbf{V}_r^{H}\right), \qquad \bar f_k = \frac{\mathrm{FR}}{2\pi}\,\arg\sum_{t} w_k[t+1]\,w_k^{*}[t],\ \ \mathbf{w}_k = \mathbf{v}_k^{*} \]

(31.2)

\(\mathbf{S}\): Casorati matrix, \(N_p\) pixels × \(N_t\) frames of beamformed IQ
\(\sigma_k,\ \mathbf{u}_k,\ \mathbf{v}_k\): singular values (descending), spatial and temporal singular vectors
\(r\): number of tissue components removed (lower cut-off)
\(\mathbf{V}_r\): the first \(r\) temporal singular vectors; the filter is a projection in slow time
\(\bar f_k\): mean Doppler frequency of component \(k\), from the lag-one autocorrelation of its temporal signature (Kasai, Module 11)
\(\mathbf{S}_{\mathrm{b}}\): blood (bubble) signal after removing the first \(r\) components
\(\mathbf{I}\): identity matrix
\((\cdot)^{H},\ (\cdot)^{*}\): conjugate transpose, complex conjugate
\(t\): frame index
\(\mathbf{w}_k\): temporal signature of component \(k\) (conjugated temporal singular vector)
\(\arg\): phase angle, rad
\(N_p,\ N_t\): number of pixels and of frames in the block

Only \(\mathbf{V}_r\) is needed, and it comes from the small \(N_t \times N_t\) Gram matrix \(\mathbf{S}^{H}\mathbf{S}\): that is how fast SVD filters run, and how Demo 31.1 runs in a browser. In m31_ulm_pipeline.m the first three singular values sit at 0, −20.7 and −31.7 dB with mean Doppler 0.2, 0.3 and 3.7 Hz; the fourth is 17 dB lower (−48.5 dB) and its mean Doppler is 57 Hz: bubbles.

Choosing r. Three signatures mark the tissue subspace: a knee in the singular spectrum, a low mean Doppler \(\bar f_k\) of the temporal vectors, and spatial singular vectors that look like the B-mode image (Baranger et al., IEEE Trans. Med. Imaging 2018, use the similarity of \(\lvert\mathbf{u}_k\rvert\) between components to find the boundary automatically). The cost of each extra component is asymmetric. The cut-off sweep of m31_ulm_pipeline.m on the same block:

components removed r	0	1	2	3	4	6	8	12
recall (detected / true)	0.00	0.00	0.10	0.45	0.43	0.39	0.38	0.33
false alarms per frame	0.64	0.01	0.00	0.26	2.45	6.01	5.42	5.22
median localization error, µm	–	110	25	26	27	32	36	49

Too few components and the residual clutter sets the noise estimate, so the detector sees nothing (r = 0–2). The right cut-off (r = 3, the clutter rank of this simulation) gives the best recall and almost no false alarms. Each component beyond it removes some bubble energy: slow bubbles lose their mean, which leaves a negative imprint of their path in every frame of the block; those imprints are detected as false alarms and bias the localizations of their neighbours. PALA's in vivo example script keeps components 5 onward (cut-off [5 N]) and states that it must be adapted to the clutter and SNR.

Which bubbles does the filter remove? A bubble is tissue to the SVD if, during one block, it neither rotates its phase nor leaves its PSF. Its temporal bandwidth is roughly the larger of its Doppler frequency and its transit rate through the PSF:

\[ B_{\mathrm{b}} \approx \max\!\left(\frac{2\lvert v_z\rvert}{\lambda},\ \frac{\lvert v_x\rvert}{\mathrm{FWHM}_x}\right) \gtrsim \frac{1}{T_b} \;\Longrightarrow\; \lvert v_z\rvert_{\min} \approx \frac{\lambda}{2T_b}, \qquad \lvert v_x\rvert_{\min} \approx \frac{\mathrm{FWHM}_x}{T_b} \]

(31.3)

\(B_{\mathrm{b}}\): slow-time bandwidth of the bubble signal, Hz
\(v_z,\ v_x\): axial and lateral bubble velocity, m/s
\(T_b\): block duration \(N_t/\mathrm{FR}\), s
\(\mathrm{FWHM}_x\): lateral width of the PSF, m
\(\lambda\): wavelength, m
\(\gtrsim\): greater than or about equal to

With 100-frame blocks at 1 kHz (\(T_b\) = 0.1 s) and λ = 308 µm: about 1.5 mm/s axially and 3.1 mm/s laterally, so capillary flow (0.5–2 mm/s) sits at the edge. Longer blocks lower both limits (and need more motion stability); the Doppler frequency of a 1 mm/s axial bubble is only 6.5 Hz (Module 28). The estimate is order-of-magnitude: tissue motion widens the tissue band and sets the real limit.

Motion. Everything downstream assumes the tissue holds still while thousands of frames accumulate. Breathing, cardiac pulsation and probe drift move the vessels by tens to hundreds of micrometres, ten times the localization error. Motion is estimated from the tissue signal that the SVD just removed (block-wise speckle correlation, Module 11, or rigid/affine registration of B-mode blocks) and subtracted from the localizations; resampling the IQ is unnecessary. Whatever is not corrected adds in quadrature to the localization error:

\[ \sigma_{\mathrm{eff}}^{2} = \sigma_{\mathrm{loc}}^{2} + \sigma_{\mathrm{m}}^{2} + \sigma_{\mathrm{reg}}^{2}, \qquad \sigma_{\mathrm{m}} = \frac{A}{\sqrt2}\ \ \text{(sinusoid of amplitude } A\text{)}, \qquad s_{\mathrm{vessel}}^{2} = \frac{R^{2}}{5} + \sigma_{\mathrm{eff}}^{2} \]

(31.4)

\(\sigma_{\mathrm{loc}}\): localization error of one bubble (Section 31.4), m
\(\sigma_{\mathrm{m}}\): rms of the uncorrected motion along the direction considered, m
\(\sigma_{\mathrm{reg}}\): error of the motion estimate when correction is applied, m
\(s_{\mathrm{vessel}}\): standard deviation of localizations across a vessel of radius \(R\) when each crossing bubble is counted once (track rendering) and the profile is a 2D parabola, \(\mathrm{E}[u^{2}] = 1/5\), the model of the simulations here; a cylindrical lumen gives \(R^{2}/6\) for track counts and \(R^{2}/4\) when every frame's localization is counted (slow wall bubbles then appear in more frames and the positions are uniform over the cross-section)
\(A\): amplitude of the sinusoidal motion, m
\(R\): vessel radius, m

A 40 µm vessel has a lumen spread of 8.9 µm (8.2 to 10 µm for a cylinder); with \(\sigma_{\mathrm{loc}}\) = 10 µm its localizations spread over 13.4 µm (FWHM 32 µm if Gaussian). An uncorrected 100 µm breathing excursion adds 71 µm rms and the vessel appears 170 µm wide. m31_ulm_pipeline.m checks the quadrature sum without motion: across its 60 µm vessel the localizations spread by 25.3 µm while the true positions of the same bubbles spread by 7.9 µm (a few bubbles only in 0.3 s), implying 24.0 µm of perpendicular localization error against 23.6 µm measured directly.

Demo 31.2

Motion: two vessels 80 µm apart, with and without registration

SimulationLocalizations (20,000 over 10 s) are drawn on two parallel 40 µm vessels 80 µm apart and one 20 µm vessel, flux-weighted across each lumen, with Gaussian localization error, then displaced by a breathing motion \(A\sin(2\pi f t)\) along the chosen direction. Registration: every 100 ms a synthetic tissue speckle frame (64 × 64 pixels of λ/4, PSF as in Demo 31.1) is shifted by the true motion in the Fourier domain, corrupted by additive noise at the chosen tissue SNR, and registered to the first frame by FFT cross-correlation with a 3-point parabolic peak fit (Module 11); the estimated shift is interpolated in time and subtracted from each localization. In-plane rigid motion only; no out-of-plane motion, no deformation.

no motionuncorrectedregistered and corrected

Motion amplitude A 150 µm

Motion direction (0° = across the vessels) 30°

Motion frequency 0.3 Hz

Localization error σ_loc 8 µm

Tissue SNR for registration 20 dB

Try: with the defaults (150 µm at 30°) the uncorrected profile fills the gap between the two 40 µm vessels (valley 73 % of the peaks) and the 20 µm vessel spreads to 94 µm rms, as (31.4) predicts (92 µm); after registration both return to the no-motion values. Set the direction to 90° (along the vessels): nothing is lost even without correction, because motion along a vessel only slides bubbles along it. Lower the tissue SNR to 0 dB: the registration error grows to about 7 µm and adds in quadrature (spread 10.1 instead of 9.1 µm). Raise σ_loc to 30 µm with no motion: the pair starts to merge, since 80 µm is then only 2.7 σ_loc.

31.3Detection: matched filter, threshold and the false-alarm budget

After filtering, each frame holds a few dozen bubbles in white-ish noise. The optimal detector for a known pulse shape in white Gaussian noise is the matched filter (Module 10): correlate the frame with the PSF and look for peaks. Normalized cross-correlation (NCC) with the PSF, used by several groups, is the same correlation divided by the local energy, which makes the threshold independent of bubble brightness at the price of detecting noise peaks wherever the energy is low. With an unknown phase the statistic is the magnitude, and its noise-only distribution sets the threshold:

\[ y(\mathbf{x}) = \sum_{\mathbf{p}} s_{\mathrm{b}}(\mathbf{p})\,g(\mathbf{p}-\mathbf{x}), \qquad P_{\mathrm{FA}} = e^{-\eta^{2}/\sigma_y^{2}}, \qquad P_{\mathrm{D}} = Q_1\!\left(\sqrt{2\,\mathrm{ENR}},\ \sqrt{2\ln(1/P_{\mathrm{FA}})}\right), \qquad \bar N_{\mathrm{FA}} \approx P_{\mathrm{FA}}\,\frac{A_{\mathrm{FOV}}}{\mathrm{FWHM}_x\,\mathrm{FWHM}_z} \]

(31.5)

\(y(\mathbf{x})\): matched-filter output at pixel \(\mathbf{x}\) (complex)
\(g\): PSF envelope used as template
\(\eta\): threshold on \(\lvert y\rvert\)
\(\sigma_y^{2}\): noise power at the filter output; \(\lvert y\rvert^{2}\) is exponential under noise only, so the median gives \(\sigma_y^{2} = \mathrm{median}\lvert y\rvert^{2}/\ln 2\)
\(Q_1\): Marcum Q-function: detection probability of a constant-amplitude, random-phase echo
\(\mathrm{ENR}\): output SNR of the matched filter = energy-to-noise ratio of the bubble
\(\bar N_{\mathrm{FA}}\): expected false alarms per frame, counting independent resolution cells in the field \(A_{\mathrm{FOV}}\)
\(\mathbf{x},\ \mathbf{p}\): output pixel and summation pixel, m
\(s_{\mathrm{b}}(\mathbf{p})\): beamformed bubble-filtered IQ signal at pixel \(\mathbf{p}\)
\(P_{\mathrm{FA}},\ P_{\mathrm{D}}\): false-alarm and detection probability per resolution cell
\(A_{\mathrm{FOV}}\): area of the field of view, m²
\(\mathrm{FWHM}_x,\ \mathrm{FWHM}_z\): lateral and axial PSF widths, m
\(\ln,\ e\): natural logarithm, Euler’s number

\(P_{\mathrm{FA}} = 10^{-5}\) needs \(\eta^{2}/\sigma_y^{2}\) = 11.5. Detection is then 50 % at an ENR of 10.4 dB, 82 % at 12 dB and 99.9 % at 14.8 dB: the transition is sharp, so the bubbles that matter are lost to merging and filtering, not to the threshold. The 4.9 × 4.9 mm field of Demo 31.1 holds 341 cells of λ × 0.75λ, so noise alone gives 0.003 false alarms per frame; the 0.26 measured in m31_ulm_pipeline.m come from residual clutter and the filter's ghosts.

Two practical detectors depart from the textbook. PALA's localization code keeps a fixed number of the brightest local maxima per frame (its numberOfParticles, 30–100) with a safeguard on the number of local maxima inside the FWHM mask: in effect a fixed-count detector whose false-alarm rate rises whenever the frame holds fewer bubbles than the count. Others threshold the NCC with the PSF. Either way the detector cannot separate two bubbles closer than about one FWHM: they produce one peak. In m31_ulm_pipeline.m only 35 % of the detected bubbles have no neighbour within 1.5 FWHM, and the missed ones are mostly pairs that produced one peak. This is the density trade-off of Module 28, \(P_{\mathrm{iso}} = e^{-nA}\): at most \(1/(eA)\) isolated bubbles per unit area per frame, 1.65 per mm² at 5 MHz.

31.4Localization: four estimators and one bound

Module 12 derived the bound for an isolated bubble with a Gaussian PSF of amplitude standard deviations \(w_x, w_z\) (FWHM/2.3548), unknown complex amplitude and white noise of density \(N_0\):

\[ \sigma_{x} \ge \frac{w_x}{\sqrt{\mathrm{ENR}}}, \qquad \sigma_{z} \ge \frac{w_z}{\sqrt{\mathrm{ENR}}}, \qquad \mathrm{ENR} = \frac{\lvert a\rvert^{2}\,\pi\,w_x w_z}{N_0}, \qquad \mathrm{SNR}_{\mathrm{peak}} = \mathrm{ENR}\,\frac{\Delta^{2}}{\pi w_x w_z} \]

(31.6)

\(w_x,\ w_z\): lateral and axial standard deviation of the amplitude PSF, m (FWHM/2.3548)
\(a\): complex bubble amplitude (unknown, nuisance)
\(N_0\): complex noise power per unit area (per-pixel variance × pixel area)
\(\Delta\): pixel size, m
\(\mathrm{SNR}_{\mathrm{peak}}\): SNR of the brightest pixel
\(\sigma_x,\ \sigma_z\): lateral and axial localization standard deviation, m (lower bound set by the PSF width and the ENR)
\(\lvert\cdot\rvert\): magnitude

At 5 MHz (\(w_x\) = 131 µm, \(w_z\) = 98 µm) an ENR of 20 dB gives 13.1 × 9.8 µm, 30 dB gives 4.1 × 3.1 µm; reaching 10 µm laterally takes 22.3 dB. On λ/4 pixels the peak pixel is 8.3 dB below the ENR. The bound is local: below about 15 dB noise peaks start to win (Section 12.4), and the bound no longer describes the error.

The estimators in common use, all applied to a small window around the brightest pixel of \(\lvert s_{\mathrm{b}}\rvert\) (PALA's code benchmarks seven: no shift, weighted average, cubic, Lanczos and spline interpolation, Gaussian fit and radial symmetry):

\[ \begin{gathered} \hat x_{\mathrm{WA}} = \frac{\sum_i A_i x_i}{\sum_i A_i}, \qquad \hat\delta_{\mathrm{G}} = \frac{\ln A_{-1} - \ln A_{+1}}{2\left(\ln A_{-1} - 2\ln A_{0} + \ln A_{+1}\right)}, \\ \hat{\mathbf{x}}_{\mathrm{RS}} = \arg\min_{\mathbf{c}} \sum_k w_k\, d^{2}(\mathbf{c}, \ell_k), \qquad \hat{\mathbf{x}}_{\mathrm{ML}} = \arg\max_{\mathbf{x}} \frac{\bigl\lvert \langle \mathbf{y}, \mathbf{g}_{\mathbf{x}} \rangle \bigr\rvert^{2}}{\lVert \mathbf{g}_{\mathbf{x}} \rVert^{2}} \end{gathered} \]

(31.7)

\(A_i,\ x_i\): envelope and position of pixel \(i\) in the window
\(\hat\delta_{\mathrm{G}}\): sub-pixel offset of the parabola through \(\ln A\) at the peak and its two neighbours (exact for a noise-free Gaussian)
\(\ell_k,\ w_k\): line through the \(k\)-th pixel corner along the local intensity gradient, and its weight (gradient² over distance); radial symmetry (Parthasarathy, Nature Methods 2012) finds the point closest to all lines
\(\mathbf{g}_{\mathbf{x}}\): PSF template centred at \(\mathbf{x}\); the ML estimate with unknown complex amplitude maximizes the normalized matched-filter output
\(\hat x_{\mathrm{WA}}\): weighted-average (centroid) position estimate, m
\(A_{-1},\ A_0,\ A_{+1}\): envelope at the peak pixel and its two neighbours along one axis
\(\ln\): natural logarithm
\(\hat{\mathbf{x}}_{\mathrm{RS}},\ \hat{\mathbf{x}}_{\mathrm{ML}}\): radial-symmetry and maximum-likelihood position estimates
\(\mathbf{c},\ d(\mathbf{c}, \ell_k)\): candidate centre and its distance to line \(\ell_k\)
\(\mathbf{y}\): local IQ image patch
\(\langle\cdot,\cdot\rangle,\ \lVert\cdot\rVert\): inner product and norm
\(\arg\min,\ \arg\max\): argument that minimizes or maximizes the cost

m31_localization_methods.m runs all of them (plus cubic interpolation) on one bubble and compares with the bound. On λ/4 pixels with white pixel noise, at 25 dB: weighted average 2.86×, Gaussian 3-point 3.39×, cubic interpolation 4.14×, radial symmetry 1.76×, ML 0.98× the CRLB. On λ/2 pixels the Gaussian fit improves to 1.34× and radial symmetry to 1.19×, because three pixels then span most of the PSF. Two lessons follow. Fine pixels do not help a 3-point estimator: it uses three samples, and their noise is independent while the PSF is wide. And the weighted average floors at 11–15 µm at high ENR because it is biased toward the window centre (pixel locking: the sub-pixel part of its estimates has a standard deviation of 0.19 pixel instead of the 0.29 of a uniform distribution). The ML estimate reaches the bound from about 20 dB at every pixel size but costs a search per bubble; that cost, at millions of bubbles, is why simpler estimators persist. Note that beamformed noise is band-limited, not white; Module 12 showed that this narrows the gap between estimators.

Overlap. The bound assumes one bubble. Two equal bubbles at 30 dB each merge into a single peak in 63 % of cases at 0.5 FWHM separation and 40 % at 0.75 FWHM; at 1 FWHM they are detected separately but each estimate is pulled 44 µm toward its partner, more than ten times the bound; beyond 2 FWHM the bias vanishes (m31_localization_methods.m, part 2). This is the regime where learned and sparsity-based localizers earn their place (Section 31.7).

Demo 31.3

Localization precision versus ENR: weighted average, Gaussian fit, radial symmetry and ML against the CRLB

SimulationMonte Carlo in JS: one bubble per trial with uniform random sub-pixel position and random phase, Gaussian PSF (FWHM λ × 0.75λ at 5 MHz), white circular complex noise per pixel scaled to the ENR of (31.6); peak search within ±3 pixels of the true position. WA: box of ±FWHM/2; Gaussian: 3-point log-parabola; RS: radial symmetry on a ±0.75λ patch; ML: separable template search over ±1.5 pixels or ±0.45λ, whichever is larger (31 coarse steps, then 21 fine steps), on a ±λ window. Same estimators as m31_localization_methods.m.

CRLBweighted averageGaussian 3-pointradial symmetrymaximum likelihood

Pixel size

Trials per ENR 200

ENR for histogram and readout 30 dB

Try: at λ/4 the ML curve lies on the bound from 20 dB while the Gaussian fit runs about 3× above it; switch to λ/2 and the Gaussian fit closes most of the gap; switch to λ/8 and it gets worse, not better. Look at the histogram at 42 dB: the weighted average piles its estimates toward the pixel centre (pixel locking), the other estimators spread evenly. Below 14 dB every curve leaves the bound: noise peaks are now picked instead of the bubble (the threshold effect of Module 12).

31.5Tracking: nearest neighbour, optimal assignment and a Kalman predictor

Linking localizations into tracks does three jobs at once. It rejects false alarms (noise does not move along vessels for 15 frames), it measures velocity, and it lets the renderer draw continuous paths instead of dots. Between two frames the problem is a bipartite matching: localizations \(\mathbf{x}_i\) of frame \(t\) against \(\mathbf{y}_j\) of frame \(t+1\). Nearest neighbour takes the closest candidate for each bubble independently; the Hungarian algorithm (Kuhn–Munkres) finds the one-to-one matching of minimal total cost; a predictor moves each \(\mathbf{x}_i\) to where its velocity says it should be before matching. For Poisson-distributed bubbles the crowding is captured by one number:

\[ \rho = n\,d^{2}, \qquad P_{\mathrm{c,NN}} = e^{-\pi\rho}, \qquad \min_{\pi} \sum_i C_{i\pi(i)}, \quad C_{ij} = \bigl\lVert \mathbf{y}_j - (\mathbf{x}_i + \hat{\mathbf{v}}_i\,\Delta t) \bigr\rVert^{2} \ \text{(gated)}, \qquad P_{\mathrm{c,pred}} \approx \frac{1}{1 + \pi\rho\, d_{\mathrm{eff}}^{2}/d^{2}} \]

(31.8)

\(n\): areal density of localizations per frame, m⁻²
\(d\): displacement between frames, \(v/\mathrm{FR}\), m
\(\rho\): bubbles per displacement squared (dimensionless crowding)
\(P_{\mathrm{c,NN}}\): probability that no other bubble is closer than the true partner: nearest neighbour is right
\(\pi\): a one-to-one assignment of frame-\(t\) to frame-\((t+1)\) localizations; the Hungarian algorithm finds the minimum in \(O(N^{3})\)
\(\hat{\mathbf{v}}_i\): predicted velocity (zero for plain matching)
\(d_{\mathrm{eff}}\): rms prediction residual: velocity change plus localization noise; a Gaussian residual turns the exponential into \(1/(1+\pi\rho\,d_{\mathrm{eff}}^{2}/d^{2})\)
\(\Delta t\): time between frames, \(1/\mathrm{FR}\), s
\(\mathbf{x}_i,\ \mathbf{y}_j\): position of localization \(i\) in frame \(t\) and of candidate \(j\) in frame \(t+1\), m
\(C_{ij}\): assignment cost (squared distance to the predicted position), m²
\(P_{\mathrm{c,pred}}\): probability of a correct link with prediction
\(\lVert\cdot\rVert\): Euclidean norm
\(\min_{\pi}\): minimum over all assignments \(\pi\)

m31_tracking_hungarian.m (60 bubbles on a torus, 40 trials, noise 0.1 d, velocity change 0.2 d): at ρ = 0.1 nearest neighbour is right 72.9 % of the time against 73.0 % predicted; greedy global assignment only 62.7 %, Hungarian 85.6 %, Hungarian with prediction 97.5 %. At ρ = 1 they reach 5.6, 6.0, 12.0 and 80.4 %. With prediction \(d_{\mathrm{eff}}\) = 0.40 d.

Where ULM sits on this curve. At 1 kHz, 5 localizations per mm² and 10 mm/s, ρ = 5 × 0.01² = 5×10⁻⁴: every method is right more than 99.8 % of the time. Assignment matters when \(d\) grows: low frame rates (clinical systems, 3D with fewer firings per volume), fast arterial flow, or high concentrations. In ordinary 2D preclinical ULM, wrong links come from the detector instead: two bubbles that merge for a few frames, a bubble missed in one frame, a false alarm next to a real track. That is what Demo 31.1 shows at high concentration, where nearest neighbour and Hungarian give the same result. The remedies are gap closing (allowing a link across a missing frame), a minimum track length, and a motion model.

Note the counter-intuitive middle row: greedy global assignment, which takes the globally shortest pair first, does worse than plain nearest neighbour up to \(\rho\) = 0.5 (at \(\rho\) = 1 all three plain rules are below 13 % and greedy edges ahead, 6.0 against 5.6 %). It locks in short wrong pairs early and forces their true partners into worse choices. Optimality of the total cost is what makes Hungarian better, and a prediction is what makes the cost meaningful.

\[ \mathbf{s}_{k+1} = \begin{bmatrix} \mathbf{I} & \Delta t\,\mathbf{I} \\ \mathbf{0} & \mathbf{I} \end{bmatrix}\mathbf{s}_k + \mathbf{w}_k, \qquad \mathbf{y}_k = \begin{bmatrix} \mathbf{I} & \mathbf{0} \end{bmatrix}\mathbf{s}_k + \mathbf{e}_k, \qquad \sigma_{\hat v} = \sigma_{\mathrm{loc}}\,\mathrm{FR}\,\sqrt{\frac{12}{L\,(L^{2}-1)}} \]

(31.9)

\(\mathbf{s}_k\): state of one bubble: position and velocity (4 numbers in 2D)
\(\mathbf{w}_k,\ \mathbf{e}_k\): process noise (acceleration, covariance \(\mathbf{Q}\)) and localization error (covariance \(\mathbf{R} = \sigma_{\mathrm{loc}}^{2}\mathbf{I}\))
\(\sigma_{\hat v}\): error of the least-squares velocity of a straight track of \(L\) equally spaced localizations
\(\Delta t\): time between frames, \(1/\mathrm{FR}\), s
\(\mathbf{I},\ \mathbf{0}\): identity and zero matrices (2×2 in 2D)
\(\mathbf{y}_k\): measured (localized) position at frame \(k\), m
\(\sigma_{\mathrm{loc}},\ \mathrm{FR}\): localization error (m) and frame rate (frames/s), as in (31.1), (31.4)
\(L\): number of localizations in the track

The Kalman filter predicts \(\mathbf{s}\) and its covariance one frame ahead; the predicted position defines the cost in (31.8) and the innovation covariance defines a statistical gate (a Mahalanobis distance) instead of a fixed radius. The velocity formula is the slope error of a straight-line fit: with \(\sigma_{\mathrm{loc}}\) = 10 µm, \(L\) = 15 frames and 1 kHz, \(\sigma_{\hat v}\) = 0.6 mm/s, which is why ULM velocity maps resolve capillary speeds only after averaging many tracks. PALA's tracker links with Tinevez's simpletracker (Hungarian, via Munkres) with a gate of 2 pixels and keeps tracks of at least 15 frames before interpolating and smoothing them for rendering.

Demo 31.4

Linking accuracy versus crowding: nearest neighbour, greedy, Hungarian and Hungarian with prediction

SimulationSame model as m31_tracking_hungarian.m, run in JS: N bubbles on a unit torus (Poisson positions, no edges), each moving by a fixed displacement d in a random direction; between frames the velocity turns by a random vector of rms σ_v; localizations carry Gaussian errors σ_loc per axis; no births or deaths. Own Hungarian code (shortest augmenting path with potentials, O(N³)); the prediction uses the velocity measured on the previous, correct link (the predict step of a constant-velocity Kalman filter without its covariance weighting). 12 trials per point.

e^−πρ (NN theory)nearest neighbourgreedyHungarianHungarian + prediction

Localization noise σ_loc/d 0.10

Velocity change σ_v/d per frame 0.20

Bubbles per frame 50

Crowding ρ shown in the scene 0.2

Scene method

Try: nearest neighbour follows e^−πρ closely; Hungarian lifts the curve, and the prediction moves it right by about a factor \(d^{2}/d_{\mathrm{eff}}^{2}\). Set the velocity change to 1.5 d: prediction loses its advantage (the bubble turns unpredictably). Set the noise to 0.5 d: every method degrades, prediction most, since its residual contains the noise of three localizations. In the scene, wrong links (red) appear first where two bubbles happen to cross.

31.6Rendering, saturation and resolution: when is the image finished, and how sharp is it?

Rendering. Kept tracks are drawn on a grid of λ/10 to λ/20 (PALA renders at λ/10 by default; Hingot et al. 2019 used 5 µm pixels). Tracks are interpolated so that consecutive points are less than a render pixel apart, and each track adds one count per pixel it crosses. Three maps come out of the same tracks: density (tracks per pixel, displayed with a compressive scale because big vessels carry orders of magnitude more bubbles), velocity (mean track speed per pixel), and direction (the sign of the axial velocity, which separates arterioles from venules running side by side in the cortex). The density map is a histogram, and histograms have sampling problems.

Saturation. A pixel shows a vessel only after a bubble has crossed it. If bubbles cross a vessel of diameter \(d\) at the flux \(\eta\,C\,Q\) and each one paints one of the \(\max(d,l)/l\) pixel columns across it, every column is visited as a Poisson process. The time constant is the acquisition time of Hingot et al. (Sci. Rep. 2019, their Eq. 4: \(T_{\mathrm{acq}} = d/(C\,l\,Q)\)), generalized here with the efficiency \(\eta\) and with \(\max(d,l)\) for vessels thinner than a pixel:

\[ T(d) = \frac{\max(d,l)}{\eta\,C\,Q(d)\,l}, \qquad Q = \frac{\pi d^{2}}{4}\,v, \qquad S(t) = 1 - e^{-t/T}, \qquad T_{3\mathrm{D}}(d) = \frac{\max(d,l)^{2}}{\eta\,C\,Q(d)\,l^{2}} \]

(31.10)

\(T(d)\): saturation time constant of a vessel of diameter \(d\) imaged with pixels of size \(l\), s
\(C\): bubble concentration in blood, µL⁻¹
\(Q,\ v\): volumetric flow (µL/s) and mean speed in the vessel
\(\eta\): fraction of crossing bubbles that end up in kept tracks (in plane, isolated, detected, linked)
\(S(t)\): fraction of the vessel's pixels visited after acquiring for \(t\); 90 % at \(t = T\ln 10 = 2.30\,T\)
\(T_{3\mathrm{D}}\): the same for voxels of size \(l^{3}\) (Hingot et al. Eq. 4′): every bubble now paints one of \((d/l)^{2}\) columns
\(d,\ l\): vessel diameter and pixel size, µm
\(t\): acquisition time, s
\(\max\): larger of the two arguments
\(e\): Euler’s number

With \(C\) = 200 µL⁻¹ (2×10⁵ per mL, the concentration Hingot et al. report), \(v\) = 0.2 mm/s per µm of diameter and \(\eta\) = 1 (m31_saturation_frc.m): \(T\) = 0.6 s for a 100 µm vessel, 15.9 s for 20 µm with 5 µm pixels, and 255 s for a 5 µm capillary whatever the pixel size, since one bubble every 255 s crosses it. Real \(\eta\) is far below 1, which is why Hingot et al. found that vessels below 25 µm may need more than 10 minutes. \(T\) scales as \(d^{-2}\) above the pixel size and \(d^{-3}\) below it (with \(v \propto d\)); in 3D the extra factor \(d/l\) makes small-voxel volumes much slower.

A Monte Carlo of the 20 µm vessel (four columns of 5 µm) reproduces \(S(T)\) = 0.634 against 0.632. When bubbles follow a parabolic flow profile, the wall columns are visited less often and the vessel fills more slowly (0.862 instead of 0.900 at 2.3 T): images look narrower than the lumen until they saturate, which biases diameter measurements early in an acquisition.

Demo 31.5

Saturation curves and acquisition time per vessel size

Approx. modelPoisson column-filling model (31.10): speed proportional to diameter, uniform lateral bubble positions, constant concentration and efficiency, no vessel overlap in the image plane. Markers in the right panel: global 90 % saturation times reported in the text of Hingot et al. (2019) for their rat-brain data (about 10 s at 100 µm pixels, 60 s at 20 µm, close to 4 min at 5 µm), shown for scale; they describe a whole vascular network at their own \(C\) and efficiency, not one vessel.

5 µm10 µm20 µm50 µm100 µmHingot et al. 2019 (global, reported)

Concentration C 200 /µL

Efficiency η (crossing bubbles kept) 0.10

Speed per diameter 0.20 mm/s per µm

Pixel size l (left panel) 10 µm

3D voxels (T × max(d,l)/l)

Try: with the defaults a 50 µm vessel reaches 90 % in about 30 s and a 10 µm one in about 12 minutes. Halve the pixel size: big vessels slow down in proportion, capillaries (already one pixel wide) do not. Raise the concentration tenfold: every curve moves left tenfold, but Section 31.3 says the number of isolated bubbles per frame peaks at \(n = 1/A\), so in practice η falls. Tick 3D: the 50 µm vessel at 10 µm voxels slows by a factor 5.

Resolution is not precision. "λ/10 resolution" is usually a statement about localization precision. Resolution, the smallest separation at which two vessels appear as two, also depends on how densely the localizations sample each vessel. Fourier ring correlation (FRC), borrowed from optical localization microscopy and cryo-electron microscopy, measures it from the data: split the localizations into two independent halves, render both, and correlate their spectra ring by ring:

\[ \mathrm{FRC}(q) = \frac{\sum_{\lvert\mathbf{k}\rvert \in q} F_1(\mathbf{k})\,F_2^{*}(\mathbf{k})}{\sqrt{\sum_{\lvert\mathbf{k}\rvert \in q} \lvert F_1(\mathbf{k})\rvert^{2}\,\sum_{\lvert\mathbf{k}\rvert \in q} \lvert F_2(\mathbf{k})\rvert^{2}}}, \qquad T_{1/2\,\mathrm{bit}}(n_q) = \frac{0.2071 + 1.9102/\sqrt{n_q}}{1.2071 + 0.9102/\sqrt{n_q}}, \qquad T_{2\sigma} = \frac{2}{\sqrt{n_q/2}}, \qquad T_{\mathrm{fix}} = \frac17 \]

(31.11)

\(F_1,\ F_2\): 2D Fourier transforms of the two half-images
\(q\): ring of spatial frequency (cycles/m); resolution is \(1/q\) at the first crossing of a threshold
\(n_q\): number of Fourier pixels in the ring
\(T_{1/2\,\mathrm{bit}}\): half-bit information threshold of van Heel and Schatz (J. Struct. Biol. 2005); tends to 0.172 for large rings
\(T_{2\sigma}\): two standard deviations of the correlation of pure noise (Hermitian symmetry halves the independent pixels)
\(T_{\mathrm{fix}}\): the fixed 1/7 ≈ 0.143 threshold used in optical localization microscopy (for example in the FRC implementation of Nieuwenhuizen et al., Nature Methods 2013)
\(\mathrm{FRC}(q)\): Fourier ring correlation at ring \(q\), dimensionless
\(\mathbf{k}\): spatial-frequency vector, cycles/m
\((\cdot)^{*}\): complex conjugate

Hingot et al. (IEEE Trans. Med. Imaging 2021) adapted FRC to ULM: they split by odd and even tracks, read resolution at the 2σ and half-bit curves (they note that fixed thresholds of 0.5 and 1/7 have been used), and report 11–34 µm resolution for localization precisions of 11–15 µm across five datasets. They model resolution as \(\sqrt{\delta_{\mathrm{loc}}^{2} + \delta_{\mathrm{Nyq}}^{2}}\), with a sampling term that shrinks as localizations accumulate; reaching a 10 µm sampling term took more than a minute in every one of their datasets.

How the halves are split matters. Localizations of one track lie on one streamline, a sub-resolution structure that the other half does not share unless the split puts points of the same track on both sides. m31_saturation_frc.m (40 vessels of 10–80 µm, σ = 5 µm, 1/7 threshold): splitting by track gives 64 µm with 300 tracks, 49 µm with 1000 and 36 µm with 10,000; splitting by localization reports 28–31 µm from the start. The track split is honest about undersampling; the localization split measures the precision along streamlines and ignores the sampling term. At σ = 20 µm both converge (about 60 µm), because precision then dominates. At low counts the half-bit curve is the most conservative of the three thresholds (88 µm against 64 µm for 1/7 with 300 tracks).

Demo 31.6

Fourier ring correlation: track split versus localization split, three thresholds

SimulationJS version of part 2 of m31_saturation_frc.m on a smaller field: 1.28 × 1.28 mm rendered at 5 µm (256 × 256), 16 random smooth vessels of 10–80 µm, tracks allocated in proportion to flux (∝ d³), each track 0.2–0.6 mm long on one flux-weighted streamline, a localization every 10 µm with isotropic Gaussian error σ. Halves split by odd/even track or odd/even localization; FFT by Lab.fft; FRC smoothed over 3 rings; first crossing searched from 213 µm downward.

split by tracksplit by localizationthresholds: half-bit (dashed), 2σ (dotted), 1/7 (dash-dot)

Number of tracks 1000

Localization error σ 5 µm

Try: start at 100 tracks: the track-split curve falls early (the halves contain different streamlines and many vessels are not sampled at all) while the localization split already looks sharp. Increase to 10,000 tracks: the two converge as sampling stops being the limit. Then raise σ to 25 µm at 10,000 tracks: both curves now fall at the same frequency, set by precision. The thresholds disagree most when counts are low; report which one you used.

31.7Frontiers: validation, high concentration, 3D, through the skull, and function

Validation. Every number in this module was measured against a ground truth that only simulation provides. The field adopted the same discipline with PALA (Heiles et al., Nature Biomed. Eng. 2022): seven localization algorithms scored with 18 metrics (localization errors, true and false positives, Jaccard index, processing time, gridding index, saturation) combined into a global score, on two in silico datasets (isolated PSFs; a simulated microvascular flow) and an in vivo rat brain, with three supplementary in vivo datasets (rat brain bolus, rat kidney, mouse tumour), code and data released. The lesson for anyone proposing a new localizer, tracker or filter: report it on shared data, with the same metrics, against the same baselines, and separate detection (how many) from precision (how well).

High concentration and learned localization. The density optimum \(n^{*} = 1/A\) of Section 31.3 caps the useful localizations per frame, and with them the acquisition time of (31.10). Any method that localizes overlapping bubbles shrinks the effective exclusion area \(A\). Deep-ULM (van Sloun et al., IEEE Trans. Med. Imaging 2021) trains a convolutional network on simulated frames with densities up to 260 bubbles per cm² to output positions on an 8× finer grid, and reports in vivo vessel widths of about 20–30 µm in a rat spinal cord, with centroid-based ULM breaking down as density rises (its supplementary density sweep). Sparsity-based recovery (for example SUSHI, one of the comparisons in Deep-ULM's supplementary figures) does the same with a model instead of training. Both inherit a forward model: the network its simulator (PSF, noise, clutter residue), the sparse solver its dictionary. Domain gap is to learned ULM what model error is to sound-speed tomography (Module 32).

3D. A 2D image plane integrates bubbles over the elevation slice (0.5–1 mm), so vessels that cross the slice appear as dots and the in-plane density includes bubbles from out of plane. Volumetric ULM removes both problems at a cost: a 32 × 32 matrix array (1024 elements; Heiles et al., IEEE Trans. Med. Imaging 2019) needs 1024 channels or multiplexing, while a row–column array addresses an \(N \times N\) aperture with \(2N\) channels (Jensen et al., IEEE TUFFC 2020: 62 + 62 channels at 3 MHz, localization precision of about 20 µm laterally and 9 µm axially on a point phantom), with weaker focusing and edge waves (Module 19). Saturation is slower in 3D by the factor \(\max(d,l)/l\) of (31.10), but tracks no longer leave the plane, so η rises.

Through the skull. Transcranial ULM in patients (Demené et al., Nature Biomed. Eng. 2021) imaged deep brain vessels in adults at resolutions of up to 25 µm, using adaptive speckle tracking to correct both micrometric brain motion and the aberrations of transcranial propagation; it separated systolic from diastolic flow and showed vortices in a small deep aneurysm. The skull adds a phase screen (Module 33). Its defocusing part widens the PSF and costs precision as \(1/\sqrt{\mathrm{ENR}}\); its linear part steers the beam and moves every bubble behind it:

\[ \tau_{\mathrm{err}}(x_e) \approx \tau_0 + \alpha\,x_e + \beta\,x_e^{2} \;\Longrightarrow\; \Delta x \approx z\,c\,\alpha, \qquad \Delta z \approx \frac{c\,\tau_0}{2}\ \ \text{(two-way)}, \qquad \beta \ \to\ \text{defocus: } w_x \uparrow,\ \sigma_x \uparrow \]

(31.12)

\(\tau_{\mathrm{err}}(x_e)\): residual receive delay error across the aperture after the beamformer's own delays, s
\(\alpha\): its linear gradient across the aperture, s/m
\(\Delta x,\ \Delta z\): lateral and axial displacement of every localization behind that patch of skull, m
\(\beta\): quadratic part: defocus, which widens the PSF (Module 32, (32.1))
\(x_e\): element position across the aperture, m
\(\tau_0\): constant (mean) delay error, s
\(z,\ c\): depth (m) and sound speed (m/s)
\(w_x,\ \sigma_x\): PSF width and localization std, as in (31.6)

A 10 ns delay difference across a 10 mm aperture (α = 10⁻⁶ s/m) at 50 mm depth shifts the bubble laterally by 77 µm, far beyond the localization precision, and the shift is the same for every bubble in the isoplanatic patch: an aberration is a bias field in ULM, which distorts the map without blurring it, while the defocus term blurs. Averaging more bubbles does not remove a bias; only correction does (Modules 32 and 33).

Function. Functional ULM (Renaudin et al., Nature Methods 2022) turns the map into a movie: in rats it reached 6.5 µm spatial and 1 s temporal resolution by accumulating localizations in 5 s sliding windows (1 s step) over 20 repeated stimulations (100 s of data per image) at 1000 frames/s, and found activation down to subcortical nuclei, with the largest relative flow increase (+49 %) in intraparenchymal vessels of the barrel cortex. The trade-off is a counting problem: a pixel's bubble flux in a window is a Poisson count, so

\[ \mathrm{CV} = \frac{1}{\sqrt{\nu\,t_w\,N_{\mathrm{rep}}}}, \qquad \text{detectable relative change} \approx \frac{2\sqrt2}{\sqrt{\nu\,t_w\,N_{\mathrm{rep}}}} \ \ (\text{two conditions, } 2\sigma) \]

(31.13)

\(\nu\): rate of tracks crossing the pixel (or vessel segment), s⁻¹
\(t_w\): time window per image, s
\(N_{\mathrm{rep}}\): number of repeated stimulations averaged
\(\mathrm{CV}\): coefficient of variation of the count

A vessel segment crossed by 1 track/s, with 5 s windows and 20 repetitions, collects 100 tracks: CV = 10 %, and a 28 % change is detectable at 2σ; a +49 % change in a vessel crossed by 0.3 track/s (30 counts, CV 18 %) is near that limit. This is why fULM pools pixels into vascular compartments and repeats stimuli, and why concentration (Section 31.3) and efficiency η matter as much for function as for anatomy.

Frontier	Bottleneck it attacks	Physics from this module	What it trades
Learned and sparse localization	acquisition time (density cap)	\(n^{*} = 1/A\), overlap bias (31.4)	forward-model or training-domain error
3D matrix and row–column ULM	out-of-plane bubbles and vessels	slice overlap; \(T_{3\mathrm{D}}\) of (31.10)	channels (matrix) or focusing quality (RCA); slower saturation
Transcranial ULM	skull aberration, low frequency, motion	bias field (31.12); \(w/\sqrt{\mathrm{ENR}}\) at large λ	needs correction (32, 33) and long acquisitions
Functional ULM	temporal resolution of a histogram	Poisson counting (31.13)	repetitions, pooling of pixels
Better tracking	wrong links, velocity precision	ρ, Kalman, (31.8)–(31.9)	model assumptions about flow
Slow-flow filters (nonlinear pulses, adaptive SVD)	capillaries removed with tissue	(31.3), Module 28	extra firings per frame, motion residue

31.8Visual summary

Thirteen questions a ULM expert answers in five seconds. Each figure is computed from the formula named in its badge, or re-plotted from the numbers printed by this module's MATLAB scripts.

31.9How experts think about it

A ULM expert reads a super-resolved map as the output of a chain of estimators and asks three questions of every claim. Which stage set the number (precision is the localizer, resolution is precision plus sampling, coverage is flux times time)? Which bubbles never made it into the list (slow ones in the SVD, overlapping ones in the detector, short tracks in the filter)? And what is biased rather than noisy (pixel locking, overlap, motion, aberration), since averaging more bubbles cures noise and never cures bias?

A ULM image is a histogram of estimates: every stage has a pass rate and an error.

recall 0.45links 93 % correct64 % kept in tracks

m31_ulm_pipeline.m at the best SVD cut-off detects 45 % of the bubbles in the field (the rest merge with neighbours or are too slow), 96 % of detections are real, greedy linking is 93 % correct, and the 10-frame track filter keeps 64 % of the localizations. Coverage after 0.3 s: 45 % of the lumen pixels. Quote all of these, not only the precision.

SimulationPass rates printed by m31_ulm_pipeline.m (Octave), r = 3, ENR 24 dB, 20 bubbles.

Rates depend on concentration, flow speeds and clutter; the chain is the same everywhere, the numbers are not.

Displacement per frame, not frame rate, is the tracking number: δ = v/FR.

30 µm at 30 mm/s, 1 kHzgate 2 px

The tracker compares \(\delta\) with the gate and with the spacing between bubbles (31.1), (31.8). At 1 kHz the fastest arteriolar bubbles move λ/10 per frame, far inside a 2-pixel gate; halving the frame rate doubles δ and quadruples ρ. Compounding more angles costs frame rate linearly and buys ENR linearly.

Exact formula

Ignores motion within a compounded frame (0.5 dB at 30 mm/s with 5 angles at 5 kHz) and out-of-plane motion.

To the SVD, a bubble that neither rotates its phase nor leaves its PSF within one block is tissue.

v_z,min ≈ λ/(2T_b)v_x,min ≈ FWHM/T_b1.5 / 3.1 mm/s at 0.1 s

The filter projects out the slow-time subspace of the strongest components (31.2). A bubble whose signal is nearly constant over the block lives in that subspace (31.3). Longer blocks lower both limits, provided tissue motion stays coherent over the block; otherwise the tissue rank grows and the cut-off must rise.

Exact formula

Order of magnitude only; the real limit is the tissue bandwidth, which motion widens.

Pick the cut-off where the singular spectrum drops and the Doppler of v_k jumps; every extra component costs bubbles and adds ghosts.

σ₄ 17 dB below σ₃f̄ 3.7 → 57 Hzr = 4: FA × 9

In m31_ulm_pipeline.m the third component has mean Doppler 3.7 Hz and the fourth 57 Hz, 17 dB weaker. At r = 3: recall 0.45, 0.26 false alarms per frame. At r = 4: recall 0.43, 2.45 false alarms; at r = 12: 0.33 and 5.2. Removing bubble energy leaves negative imprints of slow bubbles that the detector finds.

SimulationCut-off sweep printed by m31_ulm_pipeline.m (block 1).

One clutter model; real tissue has a gradual spectrum, which is why spatial-similarity criteria and block-wise cut-offs exist.

Uncorrected motion adds in quadrature, and only the component across a vessel hurts it.

σ_m = A/√2100 µm → vessel 170 µm wide

\(\sigma_{\mathrm{eff}}^{2} = \sigma_{\mathrm{loc}}^{2} + \sigma_{\mathrm{m}}^{2} + \sigma_{\mathrm{reg}}^{2}\) (31.4). Demo 31.2: 150 µm at 30° spreads a 20 µm vessel to 94 µm rms (predicted 92 µm); speckle registration brings it back to 9.1 µm. Motion along a vessel slides bubbles along it and costs nothing.

Exact formula

Rigid in-plane motion; out-of-plane motion and deformation cannot be undone by a shift.

Set the threshold from a false-alarm rate; detection then switches on over about 4 dB.

η²/σ² = ln(1/P_FA) = 11.550 % → 99.9 %: 10.4 → 14.8 dB

The magnitude of the matched-filter output is Rayleigh under noise and Rician with a bubble (31.5). Because the transition is sharp, recall is lost to merging and filtering, not to the threshold, and the median of \(\lvert y\rvert^{2}\)/ln 2 is a robust noise estimate when bubbles are sparse.

Exact formula

White Gaussian noise after filtering; residual clutter and ghosts are neither, so measured false-alarm rates run higher.

Precision is PSF width over √ENR; resolution is precision plus sampling.

σ_x = w_x/√ENR10 µm needs 22.3 dB at 5 MHz

The CRLB (31.6) with \(w_x\) = 131 µm gives 13 µm at 20 dB and 4 µm at 30 dB. At 15 MHz \(w_x\) is three times smaller, so the same ENR gives three times better precision, if attenuation allows the ENR. Resolution measured by FRC also needs enough tracks per vessel (Section 31.6).

Exact formula

Isolated bubble, Gaussian PSF, white noise, above the threshold region (about 15 dB).

A 3-point fit wastes most of the PSF on fine pixels; ML or radial symmetry recover it.

λ/4: Gauss 3.4×, RS 1.8×, ML 0.98×λ/2: Gauss 1.3×

m31_localization_methods.m at 25 dB: the 3-point log-parabola uses three samples whose noise is independent, so finer pixels make it worse; radial symmetry uses the gradients of a whole patch; the ML estimate correlates the full template and reaches the bound. The pipeline's localization error (3× the bound) is mostly this, not overlap.

SimulationRMSE/CRLB at 25 dB printed by m31_localization_methods.m, white pixel noise.

White pixel noise exaggerates the gap; band-limited beamformed noise narrows it (Module 12).

Weighted averages lock to pixels: look at the histogram of sub-pixel positions before trusting a map.

std 0.19 px vs 0.29 uniformfloor 11–15 µm

A centroid in a finite box is pulled toward the box centre, so estimates cluster at pixel centres and the RMSE stops improving with ENR (11.3 µm at 40 dB on λ/4 pixels). The rendered map then shows a grid pattern at the beamforming pitch. PALA's in vivo figure includes a "gridding index" for exactly this reason.

SimulationStandard deviation of the sub-pixel part of the estimates at 30 dB, λ/4, printed by m31_localization_methods.m.

Background subtraction and smaller boxes reduce, but do not remove, the bias.

Two bubbles within one FWHM are one bubble to a single-source localizer.

63 % merged at 0.5 FWHM44 µm bias at 1 FWHM

m31_localization_methods.m, two equal bubbles at 30 dB: below 1 FWHM most pairs give one peak, located near the midpoint; at 1 FWHM both are detected but each estimate is pulled 44 µm toward the other, ten times the bound; beyond 2 FWHM the bias disappears. This sets the density cap and motivates multi-source estimators.

SimulationMerging rate and bias printed by m31_localization_methods.m (part 2).

Equal amplitudes, lateral separation; unequal pairs bias the weak bubble more and the strong one less.

Linking fails as e^−πnd²; in 2D ULM at 1 kHz that is 10⁻⁴, so blame the detector first.

ρ = n d²ρ = 0.1: NN 0.729 (theory 0.730)ULM: ρ ≈ 5×10⁻⁴

Nearest neighbour is right when no other bubble is closer than the true partner (31.8). With 5 localizations per mm² and 10 mm/s at 1 kHz, ρ = 5×10⁻⁴: 99.8 % correct. Wrong links in practice come from merged, missed and false detections, which is why Demo 31.1 shows no difference between NN and Hungarian.

SimulationCurve: e^−πρ; markers: nearest neighbour printed by m31_tracking_hungarian.m.

Poisson positions and constant displacement; flow along vessels makes bubbles queue, which changes the neighbour statistics.

Greedy global assignment is worse than nearest neighbour; optimal assignment plus a prediction is the fix.

ρ = 0.2: greedy 0.415, NN 0.521Hungarian 0.716, + prediction 0.957

Taking the globally shortest pair first locks in short wrong pairs and pushes their partners into worse ones. The Hungarian algorithm minimizes the total cost; a constant-velocity prediction (the Kalman predict step, (31.9)) shrinks the residual from d to 0.40 d, worth a factor 6 in tolerable crowding.

SimulationFractions at ρ = 0.2 printed by m31_tracking_hungarian.m.

Prediction helps only if velocity is coherent over frames; turbulent or branching flow reduces the gain.

Velocity from a track improves as L^−3/2: short tracks give speeds, long tracks give profiles.

σ_v = σ_locFR√(12/(L(L²−1)))0.6 mm/s at L = 15

The least-squares slope of \(L\) equally spaced positions (31.9). With 10 µm localizations at 1 kHz, a 15-frame track gives 0.6 mm/s, comparable to capillary speed itself; a 60-frame track gives 0.07 mm/s. Velocity maps are averages over many tracks per pixel for this reason.

Exact formula

Straight, constant-speed track and independent errors; curved tracks need a model, and smoothing correlates errors.

Capillaries need minutes because flux, not frame rate, paints the pixels.

T = max(d,l)/(ηCQl)5 µm: 255 s at η = 1T ∝ d⁻² … d⁻³

Hingot et al.'s model (31.10): each crossing bubble fills one pixel column. At 200 bubbles/µL and 1 mm/s, one bubble crosses a 5 µm capillary every 255 s; with realistic efficiency it is tens of minutes. Faster frame rates add samples of the same bubbles, not new bubbles.

Exact formula

Uniform lateral filling; parabolic flow fills wall pixels later (0.862 instead of 0.900 at 2.3 T in the Monte Carlo).

Split FRC by track, and say which threshold you used.

300 tracks: 64 µm (track) vs 31 µm (loc.)half-bit 88 µm

Points of one track share a streamline; splitting by localization puts it in both halves and the FRC measures precision along streamlines, not sampling of the vessel (m31_saturation_frc.m, σ = 5 µm). Hingot et al. (2021) split by odd and even tracks and read the 2σ and half-bit curves. At low counts the three thresholds differ by tens of micrometres.

SimulationFRC resolution (1/7) printed by m31_saturation_frc.m, σ = 5 µm.

FRC is a global number; local artefacts and motion bias it, and it says nothing about missing vessels.

Through the skull, an aberration is a bias field: averaging more bubbles does not remove it.

Δx ≈ z c α10 ns over 10 mm at 50 mm: 77 µm

The linear part of the residual delay error steers the receive beam, so every bubble behind that patch of skull moves by the same amount (31.12); the quadratic part defocuses and only costs precision. Demené et al. (2021) corrected aberration and motion with adaptive speckle tracking before localizing. Correction methods are Modules 32 and 33.

Exact formula

Isoplanatic patch; beyond it the shift varies across the image and warps the map.

At most 1/(eA) isolated bubbles per area per frame: beating the density cap means beating overlap.

1.65 /mm² at 5 MHzk → k/2: × 4

\(n e^{-nA}\) peaks at \(n = 1/A\) (Module 28). Any localizer that separates bubbles at \(k\) FWHM instead of 1 shrinks \(A\) by \(k^{2}\) and raises the cap by \(1/k^{2}\), which divides acquisition time by the same factor (31.10). This is the quantitative case for learned and sparse localizers, and their burden of proof.

Approx. modelPoisson density, elliptical exclusion area π k² F_xF_z at 5 MHz.

Real vessels concentrate bubbles far above the image average; local density, not mean density, sets overlap.

Functional ULM is Poisson counting: CV = 1/√(ν t_w N_rep).

100 tracks → CV 10 %2σ change 28 %

A pixel's flux in a time window is a count (31.13); comparing two conditions doubles the variance. Renaudin et al. (2022) used 5 s windows and 20 repeated stimuli (100 s per image) and pooled vessels into compartments: the counting statistics, not the localization precision, set the functional sensitivity.

Poisson arrivals; physiological variability between repetitions adds variance on top.

31.10Cheat sheet

Concept	Formula	Ultrasound use
Frame-rate budget	\(\mathrm{FR} = f_{\mathrm{PRF}}/N_a\), \(f_{\mathrm{PRF}} \le c/(2z_{\max})\)	25.7 kHz at 30 mm; 5 angles → ≤ 5.1 kHz
Displacement per frame	\(\delta = v/\mathrm{FR}\)	30 µm at 30 mm/s, 1 kHz
Compounding gain	\(\mathrm{ENR} \propto N_a\), \(\sigma \propto N_a^{-1/2}\)	9 angles: 9.5 dB, error ÷ 3
SVD clutter filter	\(\mathbf{S}(\mathbf{I} - \mathbf{V}_r\mathbf{V}_r^{H})\)	cut where σ drops and \(\bar f_k\) jumps
Slowest surviving bubble	\(v_z \approx \lambda/(2T_b)\), \(v_x \approx \mathrm{FWHM}/T_b\)	1.5 / 3.1 mm/s for 0.1 s blocks
Motion in the error budget	\(\sigma_{\mathrm{eff}}^{2} = \sigma_{\mathrm{loc}}^{2} + \sigma_{\mathrm{m}}^{2} + \sigma_{\mathrm{reg}}^{2}\), \(\sigma_{\mathrm{m}} = A/\sqrt2\)	100 µm breathing → 71 µm rms
Spread across a vessel	\(s^{2} = R^{2}/5 + \sigma_{\mathrm{eff}}^{2}\) (2D parabola, one count per bubble; cylinder: \(R^{2}/6\) or \(R^{2}/4\))	infers σ from a vessel of known size
Threshold	\(\eta^{2}/\sigma_y^{2} = \ln(1/P_{\mathrm{FA}})\)	11.5 for \(10^{-5}\); \(\sigma_y^{2}\) = median/ln 2
Detection probability	\(Q_1(\sqrt{2\,\mathrm{ENR}}, \sqrt{2\ln(1/P_{\mathrm{FA}})})\)	50 % at 10.4 dB, 99.9 % at 14.8 dB
Localization CRLB	\(\sigma_x = w_x/\sqrt{\mathrm{ENR}}\)	13 µm at 20 dB, 5 MHz
Gaussian 3-point fit	\(\hat\delta = \frac{\ln A_{-1} - \ln A_{+1}}{2(\ln A_{-1} - 2\ln A_0 + \ln A_{+1})}\)	3.4× CRLB at λ/4, 1.3× at λ/2
Density cap	\(n e^{-nA}\), max \(1/(eA)\) at \(n = 1/A\)	1.65 isolated /mm² per frame at 5 MHz
Linking	\(P_{\mathrm{c,NN}} = e^{-\pi n d^{2}}\)	ρ = 5×10⁻⁴ at 1 kHz, 5 /mm², 10 mm/s
Track velocity error	\(\sigma_v = \sigma_{\mathrm{loc}}\mathrm{FR}\sqrt{12/(L(L^{2}-1))}\)	0.6 mm/s for L = 15, 10 µm
Saturation	\(T = \max(d,l)/(\eta C Q l)\), \(S = 1 - e^{-t/T}\)	5 µm capillary: 255 s at η = 1
FRC thresholds	half-bit \(\frac{0.2071 + 1.9102/\sqrt n}{1.2071 + 0.9102/\sqrt n}\), \(2/\sqrt{n/2}\), 1/7	split by track
Aberration bias	\(\Delta x \approx z c \alpha\)	77 µm for 10 ns over 10 mm at 50 mm
Functional counting	\(\mathrm{CV} = 1/\sqrt{\nu t_w N_{\mathrm{rep}}}\)	100 tracks → 10 %

31.11Interview and qualifying-exam questions

I1Walk me through a ULM pipeline and name the number that limits each stage. · screen

Acquire compounded plane waves (ENR per frame and displacement per frame, \(\delta = v/\mathrm{FR}\)); SVD clutter filter (slowest bubble kept, about \(\lambda/(2T_b)\) axially); motion correction (residual adds in quadrature); detection (false-alarm rate and merging below one FWHM); localization (\(w/\sqrt{\mathrm{ENR}}\) plus estimator bias); tracking (\(e^{-\pi n d^{2}}\), gate, minimum length); rendering (saturation \(T = \max(d,l)/(\eta CQl)\), FRC). Follow-up: which stage would you improve first for capillary maps? Acquisition time, through η: more isolated bubbles per frame (overlap-tolerant localization) and fewer slow bubbles lost in the filter.

I2How do you choose the SVD cut-off, and what goes wrong on each side? · screen

Look for the drop in the singular spectrum, the jump in the mean Doppler of the temporal vectors (tissue at a few hertz, bubbles at tens), and spatial vectors that stop looking like the B-mode (Baranger et al. 2018 automate this). Too low: residual clutter dominates the noise estimate and the detector finds nothing or fires on tissue. Too high: slow bubbles are removed, and their removed mean leaves negative imprints that become false alarms and bias neighbours (m31_ulm_pipeline.m: r = 3 gives 0.26 false alarms per frame, r = 4 gives 2.45). Follow-up: what if tissue moves a lot? The tissue rank grows; use shorter or block-wise SVD, or register before filtering.

I3Derive the localization CRLB for a Gaussian PSF and give numbers at 5 MHz. · onsite

For \(s = a\,g(x - x_0)\) in complex white noise of density \(N_0\), the Fisher information on \(x_0\) is \((2/N_0)\lvert a\rvert^{2}\int\!\!\int (\partial_x g)^{2}\); with \(g = e^{-x^{2}/2w_x^{2} - z^{2}/2w_z^{2}}\) that is \(\mathrm{ENR}/w_x^{2}\) with \(\mathrm{ENR} = \lvert a\rvert^{2}\pi w_x w_z/N_0\), and the unknown phase costs nothing because \(\partial_{x_0}s\) is orthogonal to \(s\). So \(\sigma_x \ge w_x/\sqrt{\mathrm{ENR}}\): 13 µm at 20 dB, 4 µm at 30 dB (\(w_x\) = 131 µm). Follow-up: why do real localizers miss it? Few-sample estimators, pixel locking, overlap, and the threshold region below about 15 dB.

I4Your Gaussian-fit localizer runs at three times the CRLB. Is the CRLB wrong? · onsite

Probably not: a 3-point fit on λ/4 pixels uses three samples of a PSF four pixels wide, with independent noise in each. m31_localization_methods.m gives 3.4× at 25 dB for white pixel noise, 1.3× at λ/2, while the ML estimate (normalized matched-filter maximum) reaches 0.98× and radial symmetry 1.8×. Check that it is not overlap (split errors by nearest-neighbour distance) and not bias (sub-pixel histogram). Follow-up: is ML affordable? Separable templates and a two-stage grid make it a few thousand multiply-adds per bubble.

I5Nearest neighbour or Hungarian for ULM tracking? · onsite

Depends on \(\rho = n d^{2}\). NN is right with probability \(e^{-\pi\rho}\); at 1 kHz, 5 localizations per mm² and 10 mm/s, ρ = 5×10⁻⁴ and every method is above 99.8 %. Hungarian matters at low frame rates, fast flow or high density: at ρ = 0.2 NN gets 52 %, Hungarian 72 %, Hungarian with a constant-velocity prediction 96 % (m31_tracking_hungarian.m); greedy global assignment is worse than NN (42 %). Follow-up: what breaks links in practice then? Missed detections, merged bubbles and false alarms: gap closing and a Kalman gate help more than the assignment rule.

I6How long must I acquire to map 10 µm vessels with 5 µm pixels? · onsite

\(T = \max(d,l)/(\eta C Q l)\) (Hingot et al. 2019). With \(C\) = 200 /µL, \(v\) = 2 mm/s: \(Q\) = 1.57×10⁻⁴ µL/s, flux 0.031 bubbles/s, two columns, so \(T\) = 64 s at η = 1 and 90 % after 2.3 T = 2.4 min; with η = 0.1 about 25 minutes. Follow-up: does a faster frame rate help? No; flux sets new bubbles per second. Higher concentration helps until overlap caps η (Section 31.3).

I7A paper claims "λ/10 resolution." What do you ask? · research

Is it precision (from repeated localizations or track smoothness) or resolution (separation of two vessels, or FRC)? If FRC: split by tracks or by localizations, which threshold (half-bit, 2σ, 1/7), how many tracks, and was the image saturated? m31_saturation_frc.m shows the localization split reporting 28–31 µm when the track split gives 64 µm with 300 tracks. Then: motion correction residual, aberration (bias, not blur), and whether vessels below the saturation diameter are missing rather than resolved. Follow-up: one number you would add? The saturation curve.

I8Why does uncorrected motion matter more than localization precision in vivo? · research

Errors add in quadrature: \(\sigma_{\mathrm{eff}}^{2} = \sigma_{\mathrm{loc}}^{2} + \sigma_{\mathrm{m}}^{2} + \sigma_{\mathrm{reg}}^{2}\). A breathing excursion of 100 µm is 71 µm rms against 5–15 µm of localization error, so a vessel appears 170 µm wide. Correction by speckle registration of the tissue signal per block reduces it to the registration error (a few micrometres at good tissue SNR), but out-of-plane motion and deformation remain. Motion along a vessel is harmless; across it is fatal. Follow-up: how would you validate? Measure the spread across a vessel of known size: \(s^{2} = R^{2}/5 + \sigma_{\mathrm{eff}}^{2}\) for a 2D parabolic profile counted once per bubble (\(R^{2}/6\) to \(R^{2}/4\) for a cylinder, depending on how localizations are counted).

I9What changes in transcranial ULM, and why is aberration worse than blur? · research

Lower frequency (larger PSF, so worse \(w/\sqrt{\mathrm{ENR}}\)), skull attenuation (lower ENR), longer acquisitions, and a phase screen. The quadratic part of the screen defocuses (precision); the linear part steers the beam and displaces every bubble behind it by \(z c \alpha\): 77 µm for 10 ns across 10 mm at 50 mm. A bias does not average out, so the map is warped with full confidence. Demené et al. (2021) corrected aberration and brain motion with adaptive speckle tracking. Follow-up: how would you estimate the screen? From the bubbles themselves as guide stars, or from the distortion matrix (Module 33).

I10Make the case for, and against, deep-learning localization at high concentration. · research

For: the density cap \(1/(eA)\) limits isolated bubbles per frame; a localizer that separates bubbles at half an FWHM raises the cap fourfold and cuts acquisition time fourfold. Deep-ULM trained on densities up to 260 /cm² and reports 20–30 µm in vivo where centroid methods break down. Against: it inherits its simulator (PSF shape, clutter residue, noise), errors become domain gap, and it is hard to give a per-bubble uncertainty. I would ask for PALA-style benchmarks on shared data, recall and precision separately, and a test on a PSF it was not trained on. Follow-up: a physics-anchored alternative? Sparse recovery or unrolled networks with the beamformer as forward model.

I11Design a 3D ULM system: matrix or row–column array? · research

Matrix (e.g. 32 × 32 = 1024 elements, Heiles et al. 2019): good focusing in both directions, but 1024 channels or multiplexing, small aperture per channel budget, lower ENR per element. Row–column (Jensen et al. 2020: 62 + 62 channels): \(2N\) channels for an \(N \times N\) aperture, but one-way focusing in each direction, edge waves and lower contrast, so stronger clutter filtering and careful apodization. Either way saturation is slower by \(\max(d,l)/l\), volume rates are lower (larger ρ, so assignment and prediction matter), and the data rate is the system bottleneck. Follow-up: what does 3D fix? Out-of-plane bubbles and vessels, so tracks are longer and η rises.

31.12MATLAB lab

Base MATLAB only (no toolboxes); each script runs as-is in MATLAB R2019b+ and GNU Octave 8. Paste into the editor and press Run. The numbers below are the Octave output; MATLAB's random generator gives slightly different values.

Download the scripts: m31_ulm_pipeline.m · m31_localization_methods.m · m31_tracking_hungarian.m · m31_saturation_frc.m. Base MATLAB R2019b+ or GNU Octave 8, no toolboxes; all scripts in the repository.

m31_ulm_pipeline.mThe whole chain on simulated IQ: synthesis, SVD (spectrum and Doppler of the temporal vectors, cut-off sweep), matched-filter detection, Gaussian and weighted-average localization against truth, isolated versus crowded errors, spread across a vessel, greedy tracking, track filter, density map and saturation · change ENRdB, Nb, rList, Pfa, gate, Lmin · runtime ≈ 7 s

%% m31_ulm_pipeline.m
% Purpose : a complete, minimal ULM pipeline on simulated beamformed IQ data:
%           synthesis -> SVD clutter filter -> matched-filter detection ->
%           sub-pixel localization (weighted average, Gaussian 3-point) ->
%           greedy frame-to-frame assignment -> track filter -> density map.
%           Every stage is compared with the ground truth that only a
%           simulation provides (true bubble positions and identities).
% Model   : 64 x 64 pixels of lambda/4 (5 MHz, lambda = 308 um), 1000 frames/s,
%           3 blocks of 100 frames. Six vessels (120, 60, 60, 30, 30, 20 um;
%           centreline speeds 25 ... 1.5 mm/s, parabolic profile). Bubbles:
%           Gaussian PSF, FWHM lambda (lateral) x 0.75 lambda (axial),
%           log-normal amplitude (4 dB), random phase, Doppler phase -4*pi*z/lambda.
%           Clutter: 3 speckle fields with slow complex gains, 30 dB above the
%           median bubble. White complex noise set by the median-bubble ENR.
%           Same model as Demo 31.1 (different random numbers).
% Change  : ENRdB, Nb (bubbles in the field), rList (SVD cut-offs), Pfa,
%           gate (linking distance, lambda), Lmin (minimum track length).
% Expect  : three clutter components (0-4 Hz mean Doppler, then a drop of
%           17 dB to bubble components at 50-100 Hz); r = 3 is the best cut-off
%           (recall 0.45, 0.26 false alarms per frame); r = 1 or 2 detects
%           almost nothing, r >= 4 loses bubbles and adds false alarms.
%           Localization error about 3x the median-bubble CRLB, mostly the
%           3-point estimator on lambda/4 pixels with white pixel noise (it
%           uses 3 samples of a PSF 4 pixels wide), also for isolated bubbles.
%           Greedy links 93 % correct. Across the 60 um vessel the localizations
%           spread 25.3 um while the true positions of the same (few) bubbles
%           spread 7.9 um: the quadrature difference, 24.0 um, matches the
%           actual rms perpendicular error (23.6 um). Runtime about 10 s (Octave).
rng(1);

%% Parameters
lam_um = 308; pix = 0.25; Nx = 64; Nz = 64; Np = Nx*Nz;  % pixel = lambda/4
FR = 1000; Tb = 100; nBlk = 3;                            % frames/s, frames per block
wx = 1/2.3548/pix; wz = 0.75/2.3548/pix;                  % PSF amplitude std (pixels)
ENRdB = 24; Nb = 20;                                      % median-bubble ENR, bubbles
Pfa = 1e-5; gate = 0.25; Lmin = 10;                       % detection, linking (lambda), track length
rList = [0 1 2 3 4 6 8 12];                               % SVD cut-offs to compare (block 1)
rUse = 3;                                                 % cut-off for the full pipeline
sigp = sqrt(pi*wx*wz/10^(ENRdB/10));                      % per-pixel complex noise std

%% Vessels: centrelines resampled to 0.02 lambda, normals, lengths
P = {[0.5 2;4 3;7 5.2;10 7.8;12.5 10.2;15.5 12.6], [7 5.2;5.6 7.8;4.6 10.2;3.6 12.6;3 15.5], ...
     [10 7.8;12 6;13.6 4.2;15.5 2.2], [4.6 10.2;6.4 11.8;8.2 13.6;9.4 15.5], ...
     [13.6 4.2;12 2.6;10.6 1.6;9 0.6], [12.5 10.2;11.6 12.2;12.2 14;12 15.5]};
dV = [120 60 60 30 30 20]; vV = [25 10 8 3 2.5 1.5];      % diameter (um), centreline speed (mm/s)
nV = numel(P); VX = cell(1,nV); VZ = VX; NXv = VX; NZv = VX; Lv = zeros(1,nV);
for k = 1:nV
    Q = P{k}; t = [0; cumsum(sqrt(sum(diff(Q).^2, 2)))];
    tt = linspace(0, t(end), 2000)';
    x = interp1(t, Q(:,1), tt, 'pchip'); z = interp1(t, Q(:,2), tt, 'pchip');
    cum = [0; cumsum(sqrt(diff(x).^2 + diff(z).^2))];
    su = (0:0.02:cum(end))';
    VX{k} = interp1(cum, x, su); VZ{k} = interp1(cum, z, su); Lv(k) = su(end);
    gx = gradient(VX{k}); gz = gradient(VZ{k}); gl = sqrt(gx.^2 + gz.^2);
    NXv{k} = -gz./gl; NZv{k} = gx./gl;
end
rV = dV/lam_um/2;                                         % radius (lambda)
vf = vV/(lam_um/1000)/FR;                                 % speed (lambda per frame)
wV = dV.*Lv; cw = cumsum(wV)/sum(wV);                     % bubbles per vessel ~ d x L (demo compromise)

%% Clutter fields (speckle = white noise convolved with the PSF)
kx = exp(-((-ceil(3*wx):ceil(3*wx)).^2)/(2*wx^2)); kz = exp(-((-ceil(3*wz):ceil(3*wz))'.^2)/(2*wz^2));
ampC = 10.^([30 18 6]/20); T = cell(1,3);
for k = 1:3
    W = conv2(kz, kx, randn(Nz,Nx) + 1i*randn(Nz,Nx), 'same');
    T{k} = W/sqrt(mean(abs(W(:)).^2))*ampC(k);
end
phc = 2*pi*rand(1,3); fC = [2.1 -5.3];                    % clutter gain phases, rotation (Hz)

%% Bubbles: vessel, arclength, lateral position u, amplitude, phase, id
bk = zeros(Nb,1); bs = bk; bu = bk; ba = bk; bph = bk; bid = (1:Nb)'; nid = Nb;
for i = 1:Nb
    bk(i) = find(rand <= cw, 1); bs(i) = rand*Lv(bk(i));
    u = 2*rand - 1; while rand > 1 - u^2, u = 2*rand - 1; end; bu(i) = u;
    ba(i) = 10^(4*randn/20); bph(i) = 2*pi*rand;
end

%% 1) Synthesis
tic;
nF = nBlk*Tb; IQ = zeros(Nz, Nx, nF); truth = cell(1, nF);
xg = 0:Nx-1; zg = (0:Nz-1)';
for f = 1:nF
    t = (f-1)/FR;
    c0 = exp(-1i*4*pi*0.02*sin(2*pi*3*t + phc(1)));
    fr = T{1}*c0 + T{2}*exp(1i*(2*pi*fC(1)*t + phc(2))) + T{3}*exp(1i*(2*pi*fC(2)*t + phc(3)));
    tl = zeros(0,4);
    for i = 1:Nb
        k = bk(i); bs(i) = bs(i) + vf(k)*(1 - bu(i)^2);
        if bs(i) > Lv(k)                                   % respawn at an inlet
            k = find(rand <= cw, 1); bk(i) = k; bs(i) = 0;
            u = 2*rand - 1; while rand > 1 - u^2, u = 2*rand - 1; end; bu(i) = u;
            ba(i) = 10^(4*randn/20); bph(i) = 2*pi*rand; nid = nid + 1; bid(i) = nid;
        end
        j = min(numel(VX{k}), max(1, round(bs(i)/0.02) + 1));
        x = VX{k}(j) + bu(i)*rV(k)*NXv{k}(j); z = VZ{k}(j) + bu(i)*rV(k)*NZv{k}(j);
        X = x/pix - 0.5; Z = z/pix - 0.5;                  % pixel coordinates (0-based)
        a = ba(i)*(1 + 0.1*randn)*exp(1i*(bph(i) - 4*pi*z));
        fr = fr + a*exp(-(zg - Z).^2/(2*wz^2))*exp(-(xg - X).^2/(2*wx^2));
        tl(end+1,:) = [x z bid(i) k];                      %#ok<SAGROW>
    end
    IQ(:,:,f) = fr + sigp/sqrt(2)*(randn(Nz,Nx) + 1i*randn(Nz,Nx));
    truth{f} = tl;
end
fprintf('Synthesis: %d frames of %d x %d pixels (%.1f s)\n', nF, Nz, Nx, toc);

%% 2) SVD of each block; singular spectrum and Doppler of the temporal vectors
U = cell(1,nBlk); Sg = U; Vv = U;
for b = 1:nBlk
    X = reshape(IQ(:,:,(b-1)*Tb + (1:Tb)), Np, Tb);
    [U{b}, Sg{b}, Vv{b}] = svd(X, 'econ');
end
sv = diag(Sg{1}); w = conj(Vv{1});                         % temporal signature of component k
fd = abs(angle(sum(w(2:end,:).*conj(w(1:end-1,:)), 1)))*FR/(2*pi);
fprintf('\nBlock 1: singular values (dB re sigma_1) and |mean Doppler| of components 1-8\n');
fprintf('  k      : %s\n', sprintf('%7d', 1:8));
fprintf('  dB     : %s\n', sprintf('%7.1f', 20*log10(sv(1:8)/sv(1))));
fprintf('  |fD| Hz: %s\n', sprintf('%7.1f', fd(1:8)));

%% 3-4) Detection and localization (anonymous helpers, base MATLAB)
kmx = exp(-((-ceil(2.5*wx):ceil(2.5*wx)).^2)/(2*wx^2)); kmz = exp(-((-ceil(2.5*wz):ceil(2.5*wz))'.^2)/(2*wz^2));
thr = log(1/Pfa);                                          % eta^2/sigma^2 for Pfa = exp(-eta^2/sigma^2)
filt = @(b, r) reshape(reshape(IQ(:,:,(b-1)*Tb + (1:Tb)), Np, Tb) - U{b}(:,1:r)*Sg{b}(1:r,1:r)*Vv{b}(:,1:r)', Nz, Nx, Tb);
o3 = @(a, b, c) max(-0.5, min(0.5, 0.5*(log(a) - log(c))/min(log(a) - 2*log(b) + log(c), -eps)));
res = zeros(numel(rList), 4);                              % recall, FA/frame, RMSE (um), median error (um)
for pass = 1:2
    if pass == 1, rr = rList; blocks = 1; else, rr = rUse; blocks = 1:nBlk; end
    for ir = 1:numel(rr)
        L = zeros(0, 9);                                   % [x z frame id ex ez xWA zWA vessel]
        nTrue = 0;
        for b = blocks
            F = filt(b, rr(ir));
            for t = 1:Tb
                f = (b-1)*Tb + t; A = abs(F(:,:,t));
                M = abs(conv2(kmz, kmx, F(:,:,t), 'same')).^2;
                s2 = median(M(:))/log(2);
                C = M(4:end-3, 4:end-3); pk = C > s2*thr;
                for di = -1:1, for dj = -1:1
                    if di || dj, pk = pk & C >= M((4:end-3)+di, (4:end-3)+dj); end
                end, end
                [ii, jj] = find(pk); ii = ii + 3; jj = jj + 3;
                tl = truth{f}; inner = tl(:,1) > 0.75 & tl(:,1) < 15.25 & tl(:,2) > 0.75 & tl(:,2) < 15.25;
                nTrue = nTrue + sum(inner);
                for q = 1:numel(ii)
                    blk = A(ii(q)-1:ii(q)+1, jj(q)-1:jj(q)+1); [~, m] = max(blk(:)); [di, dj] = ind2sub([3 3], m);
                    i0 = ii(q) + di - 2; j0 = jj(q) + dj - 2;
                    i0 = min(Nz-3, max(3, i0)); j0 = min(Nx-3, max(3, j0));
                    gxo = o3(A(i0,j0-1), A(i0,j0), A(i0,j0+1)); gzo = o3(A(i0-1,j0), A(i0,j0), A(i0+1,j0));
                    Wb = A(i0-2:i0+2, j0-2:j0+2); [jw, iw] = meshgrid(-2:2, -2:2);
                    xw = (j0 - 1 + sum(Wb(:).*jw(:))/sum(Wb(:)) + 0.5)*pix; zw = (i0 - 1 + sum(Wb(:).*iw(:))/sum(Wb(:)) + 0.5)*pix;
                    x = (j0 - 1 + gxo + 0.5)*pix; z = (i0 - 1 + gzo + 0.5)*pix;
                    d = sqrt((tl(:,1) - x).^2 + (tl(:,2) - z).^2); [dm, im] = min(d);
                    if dm < 0.5, L(end+1,:) = [x z f tl(im,3) x-tl(im,1) z-tl(im,2) xw zw tl(im,4)]; %#ok<SAGROW>
                    else, L(end+1,:) = [x z f -1 0 0 xw zw 0]; end %#ok<SAGROW>
                end
            end
        end
        tp = L(:,4) >= 0; [~, uq] = unique(L(tp,3)*1e6 + L(tp,4)); E = L(tp,:); E = E(uq,:);
        er = sqrt(E(:,5).^2 + E(:,6).^2)*lam_um;
        if pass == 1
            if isempty(er), er = NaN; end
            res(ir,:) = [size(E,1)/max(1,nTrue), sum(~tp)/(numel(blocks)*Tb), sqrt(mean(er.^2)), median(er)];
        end
    end
end
fprintf('\nSVD cut-off sweep on block 1 (Pfa = %.0e, Gaussian 3-point localizer)\n', Pfa);
fprintf('   r   recall   FA/frame   RMSE(um)   median(um)\n');
for ir = 1:numel(rList), fprintf('  %2d   %5.2f    %7.2f    %7.1f    %7.1f\n', rList(ir), res(ir,:)); end

crx = wx*pix/sqrt(10^(ENRdB/10))*lam_um; crz = wz*pix/sqrt(10^(ENRdB/10))*lam_um;
tp = L(:,4) >= 0;
eGx = L(tp,5)*lam_um; eGz = L(tp,6)*lam_um;
idx = find(tp); tl0 = zeros(numel(idx),2);
for q = 1:numel(idx), tl = truth{L(idx(q),3)}; m = find(tl(:,3) == L(idx(q),4), 1); tl0(q,:) = tl(m,1:2); end
eWx = (L(tp,7) - tl0(:,1))*lam_um; eWz = (L(tp,8) - tl0(:,2))*lam_um;
fprintf('\nFull pipeline, r = %d, %d blocks (%.1f s of data): %d localizations, %d true positives\n', rUse, nBlk, nF/FR, size(L,1), sum(tp));
fprintf('  CRLB of the median bubble (ENR %d dB): sigma_x = %.1f um, sigma_z = %.1f um\n', ENRdB, crx, crz);
fprintf('  Gaussian 3-point: RMSE x = %.1f um, z = %.1f um; median |error| = %.1f um\n', sqrt(mean(eGx.^2)), sqrt(mean(eGz.^2)), median(sqrt(eGx.^2 + eGz.^2)));
fprintf('  Weighted average: RMSE x = %.1f um, z = %.1f um; median |error| = %.1f um\n', sqrt(mean(eWx.^2)), sqrt(mean(eWz.^2)), median(sqrt(eWx.^2 + eWz.^2)));

%% Isolated versus crowded bubbles (nearest other true bubble beyond 1.5 lateral FWHM)
nn = zeros(numel(idx),1); amp = nn;
for q = 1:numel(idx)
    tl = truth{L(idx(q),3)}; m = tl(:,3) == L(idx(q),4);
    d = sqrt((tl(~m,1) - tl0(q,1)).^2 + ((tl(~m,2) - tl0(q,2))/0.75).^2); if isempty(d), d = Inf; end
    nn(q) = min(d);
end
iso = nn > 1.5; eG = sqrt(eGx.^2 + eGz.^2);
fprintf('  isolated (%d %%): RMSE x = %.1f, z = %.1f um, median |error| %.1f um; crowded: median |error| %.1f um\n', ...
    round(100*mean(iso)), sqrt(mean(eGx(iso).^2)), sqrt(mean(eGz(iso).^2)), median(eG(iso)), median(eG(~iso)));

%% Spread across the 60 um vessel (vessel 2) versus the prediction
% Lateral positions u are drawn with p(u) ~ 1 - u^2 (flux weighting of a 2D parabolic profile), so the
% model spread of the lumen is R/sqrt(5). The true positions of the matched bubbles give the actual
% lumen spread of this sample (a long steady-state run would drift toward uniform u, R/sqrt(3), because
% slow wall bubbles stay in the field longer; a cylindrical vessel gives R/sqrt(6) and R/2 instead).
in2 = tp & L(:,9) == 2; q2 = find(in2); perp = zeros(numel(q2),1); perpT = perp;
for q = 1:numel(q2)
    x = L(q2(q),1); z = L(q2(q),2); dd = (VX{2} - x).^2 + (VZ{2} - z).^2; [~, j] = min(dd);
    perp(q) = ((x - VX{2}(j))*NXv{2}(j) + (z - VZ{2}(j))*NZv{2}(j))*lam_um;
    m = find(idx == q2(q), 1); xt = tl0(m,1); zt = tl0(m,2);
    perpT(q) = ((xt - VX{2}(j))*NXv{2}(j) + (zt - VZ{2}(j))*NZv{2}(j))*lam_um;
end
sModel = dV(2)/2/sqrt(5);                                  % std of u*R with p(u) ~ 1 - u^2
fprintf('\nVessel 2 (60 um): %d localizations; spread across the vessel = %.1f um\n', numel(perp), std(perp));
fprintf('  true positions of the same bubbles: %.1f um (model R/sqrt(5) = %.1f um)\n', std(perpT), sModel);
fprintf('  implied perpendicular localization error = %.1f um (actual rms perpendicular error %.1f um)\n', ...
    sqrt(max(0, std(perp)^2 - std(perpT)^2)), sqrt(mean((perp - perpT).^2)));

%% 5) Tracking: greedy global assignment within the gate, then track filter
frames = L(:,3); ntr = 0; tid = zeros(size(L,1),1); nl = 0; nlt = 0; nlc = 0;
for f = 1:nF
    cur = find(frames == f); if isempty(cur), continue; end
    prv = find(frames == f-1);
    if ~isempty(prv)
        D = sqrt((L(prv,1) - L(cur,1)').^2 + (L(prv,2) - L(cur,2)').^2);
        [dsort, ord] = sort(D(:)); usedA = false(numel(prv),1); usedB = false(numel(cur),1);
        for q = 1:numel(ord)
            if dsort(q) > gate, break; end
            [a, b2] = ind2sub(size(D), ord(q));
            if ~usedA(a) && ~usedB(b2)
                usedA(a) = true; usedB(b2) = true; tid(cur(b2)) = tid(prv(a)); nl = nl + 1;
                if L(prv(a),4) >= 0 && L(cur(b2),4) >= 0, nlt = nlt + 1; nlc = nlc + (L(prv(a),4) == L(cur(b2),4)); end
            end
        end
    end
    for q = cur', if tid(q) == 0, ntr = ntr + 1; tid(q) = ntr; end, end
end
len = accumarray(tid, 1); keep = find(len >= Lmin);
fprintf('\nTracking (greedy, gate %.2f lambda = %.0f um): %d links, %.1f %% of true-true links correct\n', gate, gate*lam_um, nl, 100*nlc/max(1,nlt));
fprintf('  %d tracks, %d kept (>= %d frames), %.0f %% of localizations kept\n', ntr, numel(keep), Lmin, 100*sum(len(keep))/size(L,1));

%% 6) Render (lambda/12 grid), each track counted once per pixel
RG = 192; rp = 16/RG; den = zeros(RG);
for k = keep'
    q = find(tid == k); [~, o] = sort(L(q,3)); q = q(o); st = false(RG);
    for m = 2:numel(q)
        n = max(1, ceil(hypot(L(q(m),1) - L(q(m-1),1), L(q(m),2) - L(q(m-1),2))/(0.3*rp)));
        xs = L(q(m-1),1) + (0:n)/n*(L(q(m),1) - L(q(m-1),1)); zs = L(q(m-1),2) + (0:n)/n*(L(q(m),2) - L(q(m-1),2));
        ix = floor(xs/rp) + 1; iz = floor(zs/rp) + 1; ok = ix >= 1 & ix <= RG & iz >= 1 & iz <= RG;
        st(sub2ind([RG RG], iz(ok), ix(ok))) = true;
    end
    den = den + st;
end
mask = false(RG); [gxr, gzr] = meshgrid(((1:RG) - 0.5)*rp, ((1:RG) - 0.5)*rp);
for k = 1:nV
    for j = 1:4:numel(VX{k})
        mask = mask | ((gxr - VX{k}(j)).^2 + (gzr - VZ{k}(j)).^2 <= (max(rV(k), 0.5*rp) + 0.5*rp)^2);
    end
end
fprintf('  saturation: %.1f %% of true-lumen pixels visited by at least one kept track\n', 100*mean(den(mask) > 0));
fprintf('Total runtime %.1f s\n', toc);

%% Figures
figure('Position', [50 50 1100 650]);
subplot(2,3,1); imagesc((0:Nx-1)*pix*lam_um/1000, (0:Nz-1)*pix*lam_um/1000, 20*log10(abs(IQ(:,:,50))/max(max(abs(IQ(:,:,50)))))); caxis([-40 0]);
colormap(gray); axis image; title('raw IQ frame 50 (dB)'); xlabel('x (mm)'); ylabel('z (mm)');
subplot(2,3,2); stem(1:20, 20*log10(sv(1:20)/sv(1)), 'filled'); hold on; plot([rUse rUse] + 0.5, [-70 0], 'r--');
xlabel('component k'); ylabel('\sigma_k (dB re \sigma_1)'); title('block 1 singular values'); ylim([-70 5]);
Ff = filt(1, rUse); subplot(2,3,3); imagesc((0:Nx-1)*pix*lam_um/1000, (0:Nz-1)*pix*lam_um/1000, 20*log10(abs(Ff(:,:,50))/max(max(abs(Ff(:,:,50)))))); caxis([-30 0]);
axis image; hold on; tl = truth{50}; plot(tl(:,1)*lam_um/1000, tl(:,2)*lam_um/1000, 'r+'); q = L(:,3) == 50; plot(L(q,1)*lam_um/1000, L(q,2)*lam_um/1000, 'co');
title('filtered frame 50, detections (o), truth (+)'); xlabel('x (mm)'); ylabel('z (mm)');
subplot(2,3,4); plot(L(:,1)*lam_um/1000, L(:,2)*lam_um/1000, 'b.', 'MarkerSize', 2); axis ij; axis image; axis([0 4.93 0 4.93]);
title('all localizations'); xlabel('x (mm)'); ylabel('z (mm)');
subplot(2,3,5); hold on; for k = keep', q = find(tid == k); plot(L(q,1)*lam_um/1000, L(q,2)*lam_um/1000, '-'); end
axis ij; axis image; axis([0 4.93 0 4.93]); box on; title(sprintf('%d kept tracks', numel(keep))); xlabel('x (mm)'); ylabel('z (mm)');
subplot(2,3,6); imagesc(((1:RG) - 0.5)*rp*lam_um/1000, ((1:RG) - 0.5)*rp*lam_um/1000, log(1 + den)); axis image; colormap(gca, hot);
title('density map (log), \lambda/12 grid'); xlabel('x (mm)'); ylabel('z (mm)');

What to look for: singular values 0, −20.7, −31.7, −48.5 dB with mean Doppler 0.2, 0.3, 3.7 and 56.6 Hz; the cut-off sweep (recall 0.00, 0.00, 0.10, 0.45, 0.43, 0.39, 0.38, 0.33 for r = 0, 1, 2, 3, 4, 6, 8, 12; false alarms 0.26 per frame at r = 3, 2.45 at r = 4). Full pipeline at r = 3: 2851 localizations, 2728 true positives; Gaussian fit RMSE 32.9 / 26.3 µm (median 27.2 µm), weighted average 33.2 / 26.4 µm, against a median-bubble CRLB of 8.3 / 6.2 µm; isolated bubbles (35 %) still have a median error of 25.0 µm, crowded ones 28.5 µm. Across the 60 µm vessel the localizations spread 25.3 µm against 7.9 µm for the true positions of the same bubbles: 24.0 µm of perpendicular error implied, 23.6 µm measured. Tracking: 92.6 % of true links correct, 53 of 519 tracks kept, 64 % of localizations; saturation 45.1 %.

m31_localization_methods.mWeighted average, Gaussian 3-point, radial symmetry, cubic interpolation and ML against the CRLB at λ/4 and λ/2 pixels; pixel locking; merging and bias of two overlapping bubbles · change pixList, enrList, nTrial, sepList · runtime ≈ 16 s

%% m31_localization_methods.m
% Purpose : compare the sub-pixel localizers used in ULM (the PALA family)
%           with the Cramer-Rao bound, then show what overlap does to them.
% Model   : one bubble = Gaussian PSF (amplitude std w = FWHM/2.3548,
%           FWHM lambda lateral, 0.75 lambda axial), unknown complex amplitude
%           (random phase), white circular complex noise per pixel.
%           ENR = |a|^2 * pi * wx * wz / (sigma_pix^2 * pix^2) (integrated SNR).
%           Bound (Module 12): sigma_x >= wx / sqrt(ENR), sigma_z >= wz / sqrt(ENR).
% Estimators, all started from the brightest pixel of |IQ|:
%   WA  : amplitude-weighted average of |IQ| in a 5 x 5 box (lambda/4 pixels)
%   GAU : separable 3-point Gaussian fit (parabola through log|IQ|)
%   RS  : radial symmetry centre (Parthasarathy 2012) on a 7 x 7 patch of |IQ|
%   INT : cubic interpolation of |IQ| on a 0.05-pixel grid, then the maximum
%   MLE : maximum of |<y, g(x0, z0)>|^2 over a 9 x 9 window (grid of 0.1 px
%         over +-1.5 px, then 0.01 px; separable templates), the ML estimate for unknown complex amplitude
% Part 2: two equal bubbles at lateral separation s (in FWHM), ENR 30 dB each,
%         random phases: bias and merging of a single-bubble localizer.
% Change  : pixList, enrList, nTrial, sepList.
% Expect  : MLE reaches the bound from 20 dB at both pixel sizes (0.98x and
%           1.04x at 25 dB); at lambda/4 with white pixel noise GAU is 3.4x,
%           WA 2.9x, INT 4.1x and RS 1.8x the bound at 25 dB; at lambda/2
%           GAU and RS come within 1.2-1.3x. WA floors at 11-15 um (pixel
%           locking: its sub-pixel histogram std is 0.19 px instead of 0.29).
%           Overlap: at 0.5-0.75 FWHM 40-63 % of pairs merge into one peak;
%           at 1 FWHM the estimate is pulled 44 um toward the partner;
%           beyond 2 FWHM the bias vanishes.
% Runtime : about 10-20 s in Octave.
rng(1);

%% Parameters (lengths in wavelengths)
lam_um = 308;
fwx = 1.0; fwz = 0.75; wxl = fwx/2.3548; wzl = fwz/2.3548;
pixList = [1/4 1/2];
enrList = [10 15 20 25 30 35 40];
nTrial = 400;
names = {'WA ', 'GAU', 'RS ', 'INT', 'MLE'}; nM = numel(names);
o3 = @(a, b, c) max(-0.5, min(0.5, 0.5*(log(a) - log(c))/min(log(a) - 2*log(b) + log(c), -eps)));
rmse = zeros(numel(pixList), numel(enrList), nM); frac = cell(1, nM);
tic;
for ip = 1:numel(pixList)
    pix = pixList(ip); wx = wxl/pix; wz = wzl/pix;            % PSF std in pixels
    N = 2*ceil(3*wx) + 9; c0 = (N + 1)/2; g = (1:N) - c0;     % patch, centre index c0
    h = 4;                                                    % MLE half-window (pixels)
    for ie = 1:numel(enrList)
        sigp = sqrt(pi*wx*wz/10^(enrList(ie)/10));
        err = zeros(nTrial, nM); fx = zeros(nTrial, nM);
        for k = 1:nTrial
            x0 = rand - 0.5; z0 = rand - 0.5;                 % true offset from the centre pixel
            Y = exp(1i*2*pi*rand)*exp(-(g' - z0).^2/(2*wz^2))*exp(-(g - x0).^2/(2*wx^2)) + sigp/sqrt(2)*(randn(N) + 1i*randn(N));
            A = abs(Y); A([1:c0-4, c0+4:N], :) = 0; A(:, [1:c0-4, c0+4:N]) = 0;   % search the peak within +-3 pixels
            [~, m] = max(A(:)); [i0, j0] = ind2sub([N N], m); A = abs(Y);
            est = zeros(nM, 1);
            % WA (5 x 5 at lambda/4, 3 x 3 at lambda/2)
            q = 2*(pix < 0.3) + 1*(pix >= 0.3); B = A(i0-q:i0+q, j0-q:j0+q); [jj, ii] = meshgrid(-q:q, -q:q);
            est(1) = j0 + sum(B(:).*jj(:))/sum(B(:)) - c0 + 1i*(i0 + sum(B(:).*ii(:))/sum(B(:)) - c0);
            % GAU
            est(2) = j0 + o3(A(i0,j0-1), A(i0,j0), A(i0,j0+1)) - c0 + 1i*(i0 + o3(A(i0-1,j0), A(i0,j0), A(i0+1,j0)) - c0);
            % RS: gradients at the midpoints of a 7 x 7 patch (Roberts cross), weighted LS intersection of the gradient lines
            Pt = A(i0-3:i0+3, j0-3:j0+3);
            du = Pt(1:6, 2:7) - Pt(2:7, 1:6); dv = Pt(1:6, 1:6) - Pt(2:7, 2:7);
            gx = (du - dv)/2; gy = -(du + dv)/2; g2 = gx.^2 + gy.^2;
            [xm, ym] = meshgrid((0:5) - 2.5, (0:5) - 2.5);
            xc0 = sum(g2(:).*xm(:))/sum(g2(:)); yc0 = sum(g2(:).*ym(:))/sum(g2(:));
            wgt = g2./max(0.3, sqrt((xm - xc0).^2 + (ym - yc0).^2));
            gn = sqrt(g2) + eps; nx = -gy./gn; ny = gx./gn; cc = nx.*xm + ny.*ym;
            M2 = [sum(wgt(:).*nx(:).^2) sum(wgt(:).*nx(:).*ny(:)); sum(wgt(:).*nx(:).*ny(:)) sum(wgt(:).*ny(:).^2)];
            v = M2 \ [sum(wgt(:).*nx(:).*cc(:)); sum(wgt(:).*ny(:).*cc(:))]; v = max(-1.5, min(1.5, v));
            est(3) = j0 + v(1) - c0 + 1i*(i0 + v(2) - c0);
            % INT: cubic interpolation on a 0.05-pixel grid within +-1 pixel
            [Xq, Zq] = meshgrid(-1:0.05:1, -1:0.05:1);
            Ai = interp2(-2:2, (-2:2)', A(i0-2:i0+2, j0-2:j0+2), Xq, Zq, 'cubic');
            [~, m] = max(Ai(:)); est(4) = j0 + Xq(m) - c0 + 1i*(i0 + Zq(m) - c0);
            % MLE: |<y,g>|^2 with separable templates, coarse then fine grid
            Yw = Y(i0-h:i0+h, j0-h:j0+h); gg = (-h:h)';
            cx = 0; cz = 0;
            for stage = 1:2
                if stage == 1, cand = (-15:15)*0.1; else, cand = (-10:10)*0.01; end
                Gx = exp(-(gg - (cx + cand)).^2/(2*wx^2)); Gz = exp(-(gg - (cz + cand)).^2/(2*wz^2));
                C = abs(Gz.'*Yw*Gx).^2 ./ (sum(Gz.^2, 1)'*sum(Gx.^2, 1));
                [~, m] = max(C(:)); [a1, a2] = ind2sub(size(C), m); cz = cz + cand(a1); cx = cx + cand(a2);
            end
            est(5) = j0 + cx - c0 + 1i*(i0 + cz - c0);
            err(k,:) = (real(est).' - x0)*pix*lam_um + 1i*(imag(est).' - z0)*pix*lam_um;
            fx(k,:) = real(est).' - round(real(est).');           % sub-pixel position of the estimate
        end
        rmse(ip, ie, :) = sqrt(mean(real(err).^2, 1));
        if ip == 1 && enrList(ie) == 30, for m = 1:nM, frac{m} = fx(:,m); end, end
    end
end
crlb = wxl*lam_um./sqrt(10.^(enrList/10));
fprintf('Lateral RMSE (um) versus ENR, %d trials per point; CRLB sigma_x = w_x/sqrt(ENR)\n', nTrial);
for ip = 1:numel(pixList)
    fprintf('\npixel = lambda/%d (%.0f um)\n   ENR(dB)  CRLB  %s\n', round(1/pixList(ip)), pixList(ip)*lam_um, sprintf('%8s', names{:}));
    for ie = 1:numel(enrList), fprintf('   %4d   %6.1f %s\n', enrList(ie), crlb(ie), sprintf('%8.1f', squeeze(rmse(ip,ie,:)))); end
    fprintf('   ratio to CRLB at 25 dB: %s\n', sprintf('%8.2f', squeeze(rmse(ip, enrList == 25, :))/crlb(enrList == 25)));
end
fprintf('\nPixel locking at lambda/4, 30 dB: std of the sub-pixel part of x-estimates (uniform = %.3f px)\n', sqrt(1/12));
for m = 1:nM, fprintf('   %s  %.3f px\n', names{m}, std(frac{m})); end

%% Part 2: two overlapping bubbles (lambda/4 pixels, ENR 30 dB each)
pix = 1/4; wx = wxl/pix; wz = wzl/pix; N = 41; c0 = 21; g = (1:N) - c0; h = 4;
sepList = [0.5 0.75 1 1.25 1.5 2 3];                       % lateral separation in FWHM
sigp = sqrt(pi*wx*wz/10^(30/10)); nT2 = 300;
fprintf('\nTwo equal bubbles, lateral separation s, ENR 30 dB each; single-bubble MLE started at the brightest pixel\n');
fprintf('   s (FWHM)  s (um)   merged %%   bias toward partner (um)   RMSE (um)\n');
bias2 = zeros(size(sepList)); rm2 = bias2; mg = bias2;
for is = 1:numel(sepList)
    s = sepList(is)*fwx/pix;                                % separation in pixels
    e = zeros(nT2,1); merged = 0;
    for k = 1:nT2
        xa = -s/2 + rand - 0.5; xb = xa + s; za = rand - 0.5; zb = za;
        Y = exp(1i*2*pi*rand)*exp(-(g' - za).^2/(2*wz^2))*exp(-(g - xa).^2/(2*wx^2)) + ...
            exp(1i*2*pi*rand)*exp(-(g' - zb).^2/(2*wz^2))*exp(-(g - xb).^2/(2*wx^2)) + sigp/sqrt(2)*(randn(N) + 1i*randn(N));
        A = abs(Y);
        % count local maxima above half the global maximum (merging test)
        Cm = A(2:end-1, 2:end-1); pk = Cm > 0.5*max(A(:));
        for di = -1:1, for dj = -1:1, if di || dj, pk = pk & Cm >= A((2:end-1)+di, (2:end-1)+dj); end, end, end
        merged = merged + (sum(pk(:)) < 2);
        % localize the peak nearest bubble a
        [ii, jj] = find(pk); ii = ii + 1; jj = jj + 1; [~, q] = min((jj - c0 - xa).^2 + (ii - c0 - za).^2);
        i0 = ii(q); j0 = jj(q); Yw = Y(i0-h:i0+h, j0-h:j0+h); gg = (-h:h)'; cx = 0; cz = 0;
        for stage = 1:2
            if stage == 1, cand = (-15:15)*0.1; else, cand = (-10:10)*0.01; end
            Gx = exp(-(gg - (cx + cand)).^2/(2*wx^2)); Gz = exp(-(gg - (cz + cand)).^2/(2*wz^2));
            C = abs(Gz.'*Yw*Gx).^2 ./ (sum(Gz.^2, 1)'*sum(Gx.^2, 1));
            [~, m] = max(C(:)); [a1, a2] = ind2sub(size(C), m); cz = cz + cand(a1); cx = cx + cand(a2);
        end
        e(k) = (j0 + cx - c0 - xa)*pix*lam_um;              % positive = toward the partner (which is at +x)
    end
    bias2(is) = mean(e); rm2(is) = sqrt(mean(e.^2)); mg(is) = 100*merged/nT2;
    fprintf('   %5.2f    %5.0f    %6.1f        %8.1f                %7.1f\n', sepList(is), sepList(is)*fwx*lam_um, mg(is), bias2(is), rm2(is));
end
fprintf('Runtime %.1f s\n', toc);

%% Figures
figure('Position', [50 50 1100 380]);
subplot(1,3,1); semilogy(enrList, crlb, 'k-', 'LineWidth', 2); hold on;
mk = {'o-', 's-', 'd-', '^-', 'v-'};
for m = 1:nM, semilogy(enrList, squeeze(rmse(1,:,m)), mk{m}); end
xlabel('ENR (dB)'); ylabel('lateral RMSE (\mum)'); title('\lambda/4 pixels'); legend(['CRLB', names], 'Location', 'southwest'); grid on;
subplot(1,3,2); hist([frac{1} frac{2} frac{5}], -0.5:0.05:0.5); legend('WA', 'GAU', 'MLE');
xlabel('sub-pixel part of x estimate (pixels)'); ylabel('count'); title('pixel locking, 30 dB, \lambda/4');
subplot(1,3,3); plot(sepList, bias2, 'o-', sepList, rm2, 's-'); hold on;
plot(sepList, mg, 'k^--'); xlabel('separation (FWHM)'); ylabel('\mum  /  % merged'); legend('bias toward partner', 'RMSE', '% merged'); grid on;
title('two bubbles, 30 dB each');

What to look for: ratio to the CRLB at 25 dB, λ/4 pixels: WA 2.86, GAU 3.39, RS 1.76, INT 4.14, MLE 0.98; λ/2 pixels: 2.51, 1.34, 1.19, 2.41, 1.04. The weighted average floors at 11.3 µm (λ/4) and 14.7 µm (λ/2) at 40 dB; its sub-pixel spread is 0.192 pixel against 0.289 for a uniform distribution. Two bubbles at 30 dB: 62.7 % merged at 0.5 FWHM, 39.7 % at 0.75; bias toward the partner 37.7, 52.3, 43.6, 9.5, 2.9, 0.3 µm at 0.5, 0.75, 1, 1.25, 1.5, 2 FWHM.

m31_tracking_hungarian.mOwn Hungarian implementation checked against brute force; nearest neighbour, greedy, Hungarian and Hungarian with constant-velocity prediction versus crowding ρ = n d², with e^−πρ and 1/(1 + πρ d_eff²) · change N, rhoList, sl, sv · runtime ≈ 10 s

%% m31_tracking_hungarian.m
% Purpose : frame-to-frame linking of ULM localizations. Compares nearest
%           neighbour (NN), greedy global assignment, the optimal (Hungarian)
%           assignment and Hungarian with a constant-velocity prediction (the
%           predict step of a Kalman filter), against the Poisson prediction
%           for NN: P(correct) = exp(-pi * n * d^2).
% Model   : N bubbles on a unit torus (no edges), Poisson positions, density
%           n = N. Each bubble moves by a displacement of fixed length d in a
%           random direction; between frames its velocity turns by a random
%           vector of rms sv*d; localizations carry Gaussian errors sl*d per
%           axis. No births or deaths (every bubble has a partner).
%           Dimensionless crowding: rho = n * d^2 (bubbles per displacement^2).
% Hungarian: e-maxx / Jonker-style shortest augmenting path with potentials,
%           O(N^3), written out below (vectorized over columns), checked
%           against brute force on small problems.
% Change  : N, rhoList, nTrial, sl, sv.
% Expect  : NN follows exp(-pi rho) (0.729 measured vs 0.730 at rho = 0.1);
%           greedy global assignment is WORSE than NN (0.627): it locks in
%           short wrong pairs early. Hungarian is best without prediction
%           (0.856 at 0.1, 0.716 at 0.2). Prediction shrinks the residual to
%           0.40 d and keeps 0.96 correct at rho = 0.2 and 0.80 at rho = 1.
%           At 1 kHz, rho in ULM is 1e-5 to 2e-2 (table at the end), so in
%           practice links fail from missed or overlapping detections, not from
%           assignment. Runtime about 10 s (Octave).
rng(1);

N = 60; rhoList = [0.01 0.02 0.05 0.1 0.2 0.5 1 2]; nTrial = 40;
sl = 0.1; sv = 0.2;                                       % noise and velocity change, in units of d
names = {'NN', 'greedy', 'Hungarian', 'Hung.+pred'}; nM = 4;

%% Check the Hungarian code against brute force (6 x 6)
pm = perms(1:6); worst = 0;
for trial = 1:200
    C = rand(6);
    % ---- Hungarian (rows -> columns), minimal total cost ----
    n = 6; m = 6; u = zeros(n,1); v = zeros(1,m+1); p = zeros(1,m+1); way = zeros(1,m+1);
    for i = 1:n
        p(1) = i; j0 = 1; minv = inf(1,m+1); used = false(1,m+1);
        while true
            used(j0) = true; i0 = p(j0); cols = find(~used); cols = cols(cols > 1);
            cur = C(i0, cols-1) - u(i0) - v(cols); upd = cur < minv(cols);
            minv(cols(upd)) = cur(upd); way(cols(upd)) = j0;
            [delta, k] = min(minv(cols)); j1 = cols(k);
            uj = find(used); u(p(uj)) = u(p(uj)) + delta; v(uj) = v(uj) - delta; minv(~used) = minv(~used) - delta;
            j0 = j1; if p(j0) == 0, break; end
        end
        while true, j1 = way(j0); p(j0) = p(j1); j0 = j1; if j0 == 1, break; end, end
    end
    asg = zeros(1,n); asg(p(2:end)) = 1:m;
    costH = sum(C(sub2ind([n n], 1:n, asg)));
    costB = min(sum(C(sub2ind([n n], repmat(1:n, size(pm,1), 1), pm)), 2));
    worst = max(worst, abs(costH - costB));
end
fprintf('Hungarian versus brute force on 200 random 6 x 6 problems: max cost difference = %.2e\n', worst);

%% Linking experiment
res = zeros(numel(rhoList), nM); tic;
for ir = 1:numel(rhoList)
    d = sqrt(rhoList(ir)/N);                               % displacement for this crowding
    ok = zeros(1, nM); tot = 0;
    for trial = 1:nTrial
        x0 = rand(N,2);                                    % positions at frame t-1
        th = 2*pi*rand(N,1); v1 = d*[cos(th) sin(th)];     % velocity t-1 -> t
        v2 = v1 + sv*d/sqrt(2)*randn(N,2); v2 = d*v2./sqrt(sum(v2.^2,2));   % turned velocity t -> t+1, length d
        x1 = mod(x0 + v1, 1); x2 = mod(x1 + v2, 1);
        y0 = mod(x0 + sl*d*randn(N,2), 1); y1 = mod(x1 + sl*d*randn(N,2), 1); y2 = mod(x2 + sl*d*randn(N,2), 1);
        perm = randperm(N); y2 = y2(perm,:);               % unknown order in frame t+1; true partner of i is find(perm == i)
        truthIdx = zeros(N,1); truthIdx(perm) = 1:N;
        wrapd = @(a) a - round(a);                         % torus difference
        DX = wrapd(y1(:,1) - y2(:,1)'); DZ = wrapd(y1(:,2) - y2(:,2)'); D2 = DX.^2 + DZ.^2;
        vp = [wrapd(y1(:,1) - y0(:,1)) wrapd(y1(:,2) - y0(:,2))];   % velocity from the previous (correct) link
        PX = wrapd(y1(:,1) + vp(:,1) - y2(:,1)'); PZ = wrapd(y1(:,2) + vp(:,2) - y2(:,2)'); P2 = PX.^2 + PZ.^2;
        % NN: each bubble takes its nearest candidate (no exclusivity)
        [~, nnIdx] = min(D2, [], 2); ok(1) = ok(1) + sum(nnIdx == truthIdx);
        % greedy global: sort all pairs, accept if both are free
        [~, ord] = sort(D2(:)); ua = false(N,1); ub = false(N,1); gIdx = zeros(N,1); cnt = 0;
        for q = 1:numel(ord)
            [a, b] = ind2sub([N N], ord(q));
            if ~ua(a) && ~ub(b), ua(a) = true; ub(b) = true; gIdx(a) = b; cnt = cnt + 1; if cnt == N, break; end, end
        end
        ok(2) = ok(2) + sum(gIdx == truthIdx);
        % Hungarian on squared distances, then on squared prediction residuals
        for cm = 1:2
            if cm == 1, C = D2; else, C = P2; end
            n = N; m = N; u = zeros(n,1); v = zeros(1,m+1); p = zeros(1,m+1); way = zeros(1,m+1);
            for i = 1:n
                p(1) = i; j0 = 1; minv = inf(1,m+1); used = false(1,m+1);
                while true
                    used(j0) = true; i0 = p(j0); cols = find(~used); cols = cols(cols > 1);
                    cur = C(i0, cols-1) - u(i0) - v(cols); upd = cur < minv(cols);
                    minv(cols(upd)) = cur(upd); way(cols(upd)) = j0;
                    [delta, k] = min(minv(cols)); j1 = cols(k);
                    uj = find(used); u(p(uj)) = u(p(uj)) + delta; v(uj) = v(uj) - delta; minv(~used) = minv(~used) - delta;
                    j0 = j1; if p(j0) == 0, break; end
                end
                while true, j1 = way(j0); p(j0) = p(j1); j0 = j1; if j0 == 1, break; end, end
            end
            asg = zeros(N,1); asg(p(2:end)) = 1:m;
            ok(2 + cm) = ok(2 + cm) + sum(asg == truthIdx);
        end
        tot = tot + N;
    end
    res(ir,:) = ok/tot;
end
deff = sqrt(sv^2 + 12*sl^2);   % rms 2D prediction residual / d: velocity change plus noise of y(t+1) - 2y(t) + y(t-1) (1+4+1 per axis)
fprintf('\nFraction of correct links, N = %d per frame, %d trials, noise %.2f d, velocity change %.2f d rms\n', N, nTrial, sl, sv);
fprintf('   rho=n d^2  exp(-pi rho)  %s   1/(1+pi rho deff^2)\n', sprintf('%11s', names{:}));
for ir = 1:numel(rhoList)
    fprintf('   %6.2f     %8.3f    %s     %8.3f\n', rhoList(ir), exp(-pi*rhoList(ir)), sprintf('%11.3f', res(ir,:)), 1/(1 + pi*rhoList(ir)*deff^2));
end
fprintf('  (prediction residual: Gaussian, rms deff = sqrt(sv^2 + 12 sl^2) d = %.2f d, so NN with prediction gives 1/(1 + pi rho deff^2))\n', deff);
fprintf('Runtime %.1f s\n', toc);

%% What rho means for ULM numbers
FR = 1000; vb = [2 10 30];                                 % bubble speed, mm/s
nmm = [1 5 20];                                            % localizations per mm^2 per frame
fprintf('\nrho = n (v/FR)^2 at %d frames/s:\n   n (/mm^2)  v = %s mm/s\n', FR, sprintf('%6.0f', vb));
for k = 1:numel(nmm), fprintf('   %6.0f    %s\n', nmm(k), sprintf('%9.2e', nmm(k)*(vb/FR).^2)); end

%% Figure
figure('Position', [50 50 800 380]);
semilogx(rhoList, exp(-pi*rhoList), 'k-', 'LineWidth', 2); hold on;
mk = {'o-', 's-', 'd-', '^-'}; for m = 1:nM, semilogx(rhoList, res(:,m), mk{m}); end
semilogx(rhoList, 1./(1 + pi*rhoList*deff^2), 'k--');
xlabel('\rho = n d^2 (bubbles per displacement squared)'); ylabel('fraction of correct links');
legend(['e^{-\pi\rho} (NN theory)', names, '1/(1+\pi\rho d_{eff}^2)'], 'Location', 'southwest'); grid on; ylim([0 1.02]);
title('Linking accuracy versus crowding');

What to look for: maximum cost difference against brute force 0. At ρ = 0.1: e^−πρ = 0.730, NN 0.729, greedy 0.627, Hungarian 0.856, with prediction 0.975; at ρ = 1: 0.043, 0.056, 0.060, 0.120, 0.804. The prediction residual is 0.40 d, and Hungarian with prediction beats the NN-with-prediction formula (0.804 against 0.665 at ρ = 1). The closing table gives ρ for ULM conditions: 4×10⁻⁶ to 1.8×10⁻² at 1 kHz.

m31_saturation_frc.mSaturation time constants per vessel and pixel size (Hingot model) with a Poisson Monte Carlo, then FRC of a simulated vessel network with track and localization splits and three thresholds · change C, kv, eta, nTrackList, sigList · runtime ≈ 20 s

%% m31_saturation_frc.m
% Purpose : two ways to ask "is my ULM image finished, and how sharp is it?"
%   Part 1 - saturation. Time to fill a vessel of diameter d with pixels of
%           size l (Hingot et al., Sci. Rep. 2019, Eq. 4 and 6):
%               T(d) = max(d, l) / (eta * C * Q(d) * l),  Q = pi d^2 v / 4,
%           i.e. each bubble paints one pixel column of the max(d,l)/l columns
%           across the vessel. Saturation of that vessel: S(t) = 1 - exp(-t/T).
%           Checked by a Poisson Monte Carlo with uniform and with flux-weighted
%           (parabolic) lateral bubble positions.
%   Part 2 - Fourier ring correlation (FRC). Localizations along a random
%           vessel network are split into two halves, rendered, and their
%           ring-wise spectral correlation is read against three thresholds:
%           fixed 1/7, half-bit (van Heel and Schatz 2005) and 2-sigma.
%           Split by track (odd/even tracks, as in Hingot et al. 2021) versus
%           split by localization (odd/even points of every track).
% Model   : v = kv * d (kv = 0.2 mm/s per um: 8 um -> 1.6 mm/s, 50 um -> 10 mm/s),
%           C = 200 bubbles per uL (2e5 per mL), detection efficiency eta.
%           FRC: 2.56 x 2.56 mm field, 5 um render pixels (512 x 512); 40
%           vessels of 10-80 um, tracks allocated in proportion to flux (d^3),
%           points every 10 um along a fixed streamline, Gaussian errors sigma.
% Change  : C, kv, eta, lpix; nTrackList, sigList.
% Expect  : T = 255 s for 5 um capillaries at any pixel size, 15.9 s for a
%           20 um vessel with 5 um pixels, 0.6 s for 100 um (eta = 1). The
%           Monte Carlo matches 1 - exp(-t/T) for uniform lateral positions
%           (0.634 vs 0.632 at t = T); a parabolic profile fills the wall
%           pixels later (0.862 instead of 0.900 at 2.3 T). FRC: track-split
%           resolution improves from 64 to 36 um (sigma 5 um, 1/7) as tracks
%           grow from 300 to 10,000: it measures sampling as well as precision.
%           Splitting by localization puts the same streamline in both halves
%           and reports 28-31 um from the start: optimistic when undersampled.
%           Half-bit is the most conservative threshold at low counts.
% Runtime : about 20 s in Octave.
rng(1);

%% Part 1: saturation times
C = 200;            % bubbles per microlitre (2e5 per mL)
kv = 0.2;           % mm/s per um of diameter
eta = 1;            % fraction of passing bubbles that end up in kept tracks
dList = [5 10 20 50 100];  lList = [5 10 20];
Tsat = @(d, l) max(d, l)./(eta*C*(pi*(d*1e-6).^2/4.*(kv*d*1e-3))*1e9 .* l);   % C (/uL) x Q (uL/s; 1 m^3 = 1e9 uL)
fprintf('Saturation time constant T (s) = max(d,l)/(eta C Q l), C = %d /uL, v = %.1f mm/s per um, eta = %.2f\n', C, kv, eta);
fprintf('   d (um)   v (mm/s)   Q (nL/s)   %s\n', sprintf('   l = %2d um', lList));
for d = dList
    Q = pi*(d*1e-6)^2/4*(kv*d*1e-3)*1e12;  % nL/s
    fprintf('   %5d    %6.1f    %8.4f   %s\n', d, kv*d, Q, sprintf('%10.1f  ', Tsat(d, lList)));
end
fprintf('   time to 90 %% = T ln 10 = 2.30 T; 3D (voxel l^3) multiplies T by max(d,l)/l (Hingot Eq. 4'')\n');

% Monte Carlo check for d = 20 um, l = 5 um (4 pixel columns)
d = 20; l = 5; T = Tsat(d, l); ncol = d/l; tgrid = linspace(0, 4*T, 41); nRep = 2000;
rate = eta*C*pi*(d*1e-6)^2/4*(kv*d*1e-3)*1e9;           % bubbles per second through the vessel
Su = zeros(size(tgrid)); Sp = Su;
for r = 1:nRep
    tb = cumsum(-log(rand(1, ceil(rate*tgrid(end)*3) + 10)))/rate; tb = tb(tb <= tgrid(end));
    cu = ceil(rand(size(tb))*ncol);                        % uniform across the lumen
    up = zeros(size(tb)); for k = 1:numel(tb), u = 2*rand - 1; while rand > 1 - u^2, u = 2*rand - 1; end; up(k) = u; end
    cp = min(ncol, floor((up + 1)/2*ncol) + 1);            % flux-weighted (parabolic profile)
    for it = 1:numel(tgrid)
        Su(it) = Su(it) + numel(unique(cu(tb <= tgrid(it))))/ncol;
        Sp(it) = Sp(it) + numel(unique(cp(tb <= tgrid(it))))/ncol;
    end
end
Su = Su/nRep; Sp = Sp/nRep;
i1 = find(tgrid >= T, 1); i23 = find(tgrid >= 2.3*T, 1);
fprintf('\nMonte Carlo, d = %d um, l = %d um (T = %.1f s): S(T) = %.3f uniform, %.3f flux-weighted (model %.3f)\n', d, l, T, Su(i1), Sp(i1), 1 - exp(-1));
fprintf('                                     S(2.3T) = %.3f uniform, %.3f flux-weighted (model %.3f)\n', Su(i23), Sp(i23), 0.9);

%% Part 2: FRC on a simulated vessel network
Np = 512; pr = 5; Lf = Np*pr;                              % render pixels (um), field (um)
nV = 40; dv = exp(log(10) + rand(1,nV)*log(8));            % diameters 10-80 um, log-uniform
VX = cell(1,nV); VZ = VX; VN = VX;
for k = 1:nV                                               % smooth random curves across the field
    th = rand*pi; c = rand(1,2)*Lf; s = linspace(-Lf, Lf, 4000)';
    a = 150*randn; f = 1/(400 + 800*rand);
    x = c(1) + s*cos(th) - a*sin(2*pi*f*s)*sin(th); z = c(2) + s*sin(th) + a*sin(2*pi*f*s)*cos(th);
    in = x > 0 & x < Lf & z > 0 & z < Lf; x = x(in); z = z(in);
    if numel(x) < 50, x = c(1) + [0;1]; z = c(2) + [0;1]; end
    cum = [0; cumsum(sqrt(diff(x).^2 + diff(z).^2))]; [cum, iu] = unique(cum); x = x(iu); z = z(iu);
    su = (0:2:cum(end))'; VX{k} = interp1(cum, x, su); VZ{k} = interp1(cum, z, su);
    gx = gradient(VX{k}); gz = gradient(VZ{k}); gl = sqrt(gx.^2 + gz.^2) + eps; VN{k} = [-gz./gl gx./gl];
end
flux = dv.^3; cwf = cumsum(flux)/sum(flux);
q = (0:Np-1); q(q >= Np/2) = q(q >= Np/2) - Np; [QX, QZ] = meshgrid(q, q); QR = sqrt(QX.^2 + QZ.^2);
ring = round(QR); nR = Np/2; nring = accumarray(ring(:) + 1, 1); nring = nring(1:nR);
thrHalf = (0.2071 + 1.9102./sqrt(nring))./(1.2071 + 0.9102./sqrt(nring));
thr2s = 2./sqrt(nring/2);
qf = (0:nR-1)'/(Np*pr);                                    % cycles per um
nTrackList = [300 1000 3000 10000]; sigList = [5 10 20];
fprintf('\nFRC resolution (um) = 1/q at the first crossing; render pixel %d um, Nyquist limit %d um\n', pr, 2*pr);
fprintf('   sigma  tracks   locs     split     1/7     half-bit   2-sigma\n');
keepCurves = {};
for sg = sigList
    for nT = nTrackList
        P1 = cell(1, nT); nloc = 0;                                                     % [iz ix track-parity point-parity]
        for t = 1:nT
            k = find(rand <= cwf, 1); n = numel(VX{k}); len = round((200 + 400*rand)/2);   % 0.2-0.6 mm of track, points every 2 um of curve
            if n < len + 2, len = n - 2; end
            j0 = 1 + floor(rand*(n - len)); jj = (j0:5:j0+len)';                        % a point every 10 um
            u = 2*rand - 1; while rand > 1 - u^2, u = 2*rand - 1; end                   % fixed streamline (flux-weighted)
            px = VX{k}(jj) + u*dv(k)/2*VN{k}(jj,1) + sg*randn(numel(jj),1);
            pz = VZ{k}(jj) + u*dv(k)/2*VN{k}(jj,2) + sg*randn(numel(jj),1);
            ix = floor(px/pr) + 1; iz = floor(pz/pr) + 1; ok = ix >= 1 & ix <= Np & iz >= 1 & iz <= Np;
            P1{t} = [iz(ok) ix(ok) mod(t,2)*ones(sum(ok),1) mod((1:sum(ok))',2)];
        end
        PP = cat(1, P1{:}); nloc = size(PP,1);
        im1 = accumarray(PP(PP(:,3)==1,1:2), 1, [Np Np]); im2 = accumarray(PP(PP(:,3)==0,1:2), 1, [Np Np]);     % split by track
        im1l = accumarray(PP(PP(:,4)==1,1:2), 1, [Np Np]); im2l = accumarray(PP(PP(:,4)==0,1:2), 1, [Np Np]);   % split by localization
        for sp = 1:2
            if sp == 1, F1 = fft2(im1); F2 = fft2(im2); lab = 'tracks'; else, F1 = fft2(im1l); F2 = fft2(im2l); lab = 'locs  '; end
            num = accumarray(ring(:) + 1, real(F1(:).*conj(F2(:)))); d1 = accumarray(ring(:) + 1, abs(F1(:)).^2); d2 = accumarray(ring(:) + 1, abs(F2(:)).^2);
            frc = num(1:nR)./sqrt(d1(1:nR).*d2(1:nR)); frc = conv(frc, ones(3,1)/3, 'same'); frc(1) = 1;
            r = zeros(1,3); th = {1/7*ones(nR,1), thrHalf, thr2s};
            for k = 1:3                                    % first crossing above 1/(512 um), i.e. from ring 6
                ic = find(frc(6:end) < th{k}(6:end), 1) + 5;
                if isempty(ic), r(k) = 2*pr; elseif ic == 6, r(k) = Inf; else, r(k) = 1/qf(ic); end
            end
            fprintf('   %4d   %6d  %7d    %s   %6.1f    %6.1f    %6.1f\n', sg, nT, nloc, lab, r);
            if sg == 10 && (nT == 1000 || nT == 10000), keepCurves{end+1} = {frc, sprintf('%d tracks, %s', nT, strtrim(lab))}; end %#ok<SAGROW>
        end
    end
end

%% Figures
figure('Position', [50 50 1100 380]);
subplot(1,3,1); semilogx(dList, Tsat(dList, 5), 'o-', dList, Tsat(dList, 10), 's-', dList, Tsat(dList, 20), 'd-'); grid on;
set(gca, 'YScale', 'log'); xlabel('vessel diameter d (\mum)'); ylabel('T (s)'); legend('l = 5 \mum', 'l = 10 \mum', 'l = 20 \mum'); title('saturation time constant');
subplot(1,3,2); plot(tgrid/T, Su, 'o', tgrid/T, Sp, 's', tgrid/T, 1 - exp(-tgrid/T), 'k-'); grid on;
xlabel('t / T'); ylabel('fraction of pixel columns visited'); legend('MC uniform', 'MC flux-weighted', '1 - e^{-t/T}', 'Location', 'southeast'); title('d = 20 \mum, l = 5 \mum');
subplot(1,3,3); hold on; for k = 1:numel(keepCurves), plot(qf*1000, keepCurves{k}{1}); end
plot(qf*1000, thrHalf, 'k--', qf*1000, thr2s, 'k:', qf*1000, ones(nR,1)/7, 'k-.'); grid on; ylim([-0.1 1]);
legend([cellfun(@(c) c{2}, keepCurves, 'UniformOutput', false), {'half-bit', '2\sigma', '1/7'}]);
xlabel('spatial frequency (cycles/mm)'); ylabel('FRC'); title('FRC, \sigma = 10 \mum');

What to look for: T = 254.6 s for 5 µm at every pixel size; 63.7 / 31.8 / 31.8 s for 10 µm and 15.9 / 8.0 / 4.0 s for 20 µm at l = 5 / 10 / 20 µm; 0.6 s for 100 µm at 5 µm. Monte Carlo at d = 20, l = 5 µm: S(T) = 0.634 (uniform) and 0.609 (flux-weighted) against 0.632; S(2.3T) = 0.899 and 0.862 against 0.900. FRC (1/7 threshold): σ = 5 µm, track split 64.0, 49.2, 47.4, 36.1 µm for 300, 1000, 3000, 10,000 tracks; localization split 31.2, 28.4, 29.1, 28.4 µm. σ = 20 µm, 10,000 tracks: 59.5 (track) and 61.0 µm (localization).

31.13Self-check

Q1Depth 40 mm, 7 angles, ideal firing rate. Maximum compounded frame rate, and the displacement of a 20 mm/s bubble per frame?

\(f_{\mathrm{PRF}} \le c/(2z)\) = 1540/0.08 = 19.25 kHz; FR ≤ 2.75 kHz; δ = 20 mm/s / 2750 = 7.3 µm.

Q2Blocks of 200 frames at 500 Hz, 5 MHz. Slowest axial and lateral speeds the SVD keeps (order of magnitude)?

\(T_b\) = 0.4 s: \(v_z \approx \lambda/(2T_b)\) = 308/0.8 µm/s = 0.39 mm/s; \(v_x \approx\) FWHM/\(T_b\) = 0.77 mm/s.

Q3What ENR gives 5 µm lateral precision at 15 MHz (FWHM = λ)?

λ = 102.7 µm, \(w_x\) = 43.6 µm; ENR = (43.6/5)² = 76 = 18.8 dB.

Q4Matched-filter threshold for a false-alarm rate of 10⁻⁶, in dB above the noise power, and the false alarms per frame in 2000 independent cells?

\(\eta^{2}/\sigma_y^{2} = \ln 10^{6}\) = 13.8 = 11.4 dB; 2000 × 10⁻⁶ = 0.002 per frame.

Q510 localizations per mm² per frame, 50 mm/s, 500 Hz. Crowding ρ and nearest-neighbour link accuracy?

d = 0.1 mm, ρ = 10 × 0.01 = 0.1, \(e^{-0.314}\) = 0.73: a quarter of the links would be wrong without assignment or prediction.

Q6A 20 µm vessel, 5 µm pixels, 400 bubbles/µL, 4 mm/s, η = 0.2. Time to 90 % saturation?

\(Q\) = π(20 µm)²/4 × 4 mm/s = 1.257×10⁻³ µL/s; flux \(\eta C Q\) = 0.1 bubbles/s; \(T\) = (20/5)/0.1 = 40 s; 90 % after 2.3 T = 92 s.

Q7A breathing motion of 60 µm amplitude at 45° to a vessel, σ_loc = 8 µm, no correction. Spread across a 30 µm vessel?

Component across: 60 cos 45° = 42.4 µm, rms 30 µm. \(s^{2} = R^{2}/5 + \sigma_{\mathrm{loc}}^{2} + \sigma_{\mathrm{m}}^{2}\) = 15²/5 + 64 + 900 = 1009, s = 31.8 µm (10.4 µm without motion).

31.14Figures from the literature

Measured data behind the models of this module. Figure titles and numbers were checked against the publishers' article pages, the authors' open preprints, or (for PALA) the figure labels in the authors' public scripts; links only, no images reproduced.

The first in vivo ultrafast ULM: rat brain through a thinned skull

C. Errico, J. Pierre, S. Pezet, Y. Desailly, Z. Lenkei, O. Couture and M. Tanter, "Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging," Nature (2015). doi:10.1038/nature16066

Look at Fig. 1 (principle of ultrafast ULM: frames, bubble detection and localization, accumulation), Fig. 2 (spatial resolution and quantification in the rat cortex through a thinned skull window) and Fig. 3 (thinned skull versus intact skull). Why: the abstract's numbers (more than 500 frames per second, 75,000 images, vessels below 10 µm, velocities from 1 mm/s to several cm/s, more than 10 mm deep) are the operating point that Sections 31.1–31.6 reconstruct from first principles.

Saturation and acquisition time are set by microvascular flow

V. Hingot, C. Errico, B. Heiles, L. Rahal, M. Tanter and O. Couture, "Microvascular flow dictates the compromise between spatial resolution and acquisition time in ultrasound localization microscopy," Sci. Rep. (2019), open access. doi:10.1038/s41598-018-38349-x

Look at Fig. 3 (global saturation curves for pixel sizes and acquisition times; time constants versus pixel size fitted by \(p_1/x\), \(p_1\) = 384 s·µm), Fig. 4 (bubble counts over time in 80, 20 and 5 µm vessels) and Fig. 5 (reconstruction times of individual vessels of 90, 30 and 5 µm, and saturation time versus vessel diameter for 150 vessels). Why: the measured version of (31.10) and Demo 31.5; 192,000 frames over 4 minutes, about 6 million localizations, and vessels below 25 µm that may need more than 10 minutes.

Benchmarking localization algorithms on shared data (PALA)

B. Heiles, A. Chavignon, V. Hingot, P. Lopez, E. Teston and O. Couture, "Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy," Nat. Biomed. Eng. (2022). doi:10.1038/s41551-021-00824-8; code and data: github.com/AChavignon/PALA

Look at Fig. 2 (global score and RMSE of the seven algorithms on simulated data), Fig. 3 (lateral and axial error distributions, true/false positives and Jaccard index), Fig. 5 (rendered simulated flow against the target) and Fig. 6 (in vivo rat brain rendered with radial symmetry, and the gridding index). Figure numbers as labelled in the authors' figure scripts (fig2_…, fig3_…, fig5_…, fig6_…). Why: Sections 31.3–31.4 and Demo 31.3 on real benchmarks; the gridding index is the pixel locking of Rule 9.

Measuring resolution in ULM with Fourier ring correlation

V. Hingot, A. Chavignon, B. Heiles and O. Couture, "Measuring image resolution in ultrasound localization microscopy," IEEE Trans. Med. Imaging (2021). doi:10.1109/TMI.2021.3097150

Look at Fig. 2 (the FRC steps: split tracks into two subsets, render, compare spectra, read the 2σ crossing), Fig. 3 (FRC curves with 2σ and half-bit thresholds for brain, tumour and kidney datasets), Fig. 5 (localization precision from deviations around smoothed tracks) and Fig. 7 (FRC resolution against the sampling, or Nyquist, dimension). Why: the measured version of Demo 31.6 and Rule 15: 11–34 µm resolution for 11–15 µm precision, and the split by odd and even tracks.

Transcranial ULM of brain vessels in patients

C. Demené, J. Robin, A. Dizeux, B. Heiles, M. Pernot, M. Tanter and F. Perren, "Transcranial ultrafast ultrasound localization microscopy of brain vasculature in patients," Nat. Biomed. Eng. (2021). doi:10.1038/s41551-021-00697-x

Look at Fig. 1 (transcranial ULM of deep brain vessels in patients), Fig. 2 (estimation of the ULM resolution in a standard imaging case), Fig. 3 (haemodynamics: diastolic versus systolic flow) and Fig. 4 (a deep-seated aneurysm). Why: Section 31.7 and Rule 16 in practice: resolution up to 25 µm after correcting brain motion and transcranial aberration with adaptive speckle tracking.

Functional ULM: brain-wide activation at a microscopic scale

N. Renaudin, C. Demené, A. Dizeux, N. Ialy-Radio, S. Pezet and M. Tanter, "Functional ultrasound localization microscopy reveals brain-wide neurovascular activity on a microscopic scale," Nat. Methods (2022). doi:10.1038/s41592-022-01549-5

Look at Fig. 1 (activation maps at 6.5 µm during whisker stimulation, cortical and subcortical), Fig. 2 (flux, speed and diameter changes per vascular compartment; +49 % in intraparenchymal vessels) and Fig. 4 (SVD of sparse dynamic ULM data isolating the activation). Why: the counting trade-off of (31.13): 5 s windows, 20 stimulations, 100 s of data per image.

Learned localization at high bubble density (Deep-ULM)

R. J. G. van Sloun, O. Solomon, M. Bruce, Z. Z. Khaing, H. Wijkstra, Y. C. Eldar and M. Mischi, "Super-resolution ultrasound localization microscopy through deep learning," IEEE Trans. Med. Imaging (2021). doi:10.1109/TMI.2020.3037790; preprint arXiv:1804.07661

In the preprint: Fig. 1 (recoveries on synthetic data at increasing density), Fig. 3 (rat spinal cord in vivo with vessel profiles), Supplementary Fig. 1 (detection rate and precision versus density for Deep-ULM, centroid ULM and sparse recovery). The published version may number them differently. Why: the density cap of Rule 17 attacked with a learned localizer, and the domain-gap caveat of Section 31.7.

31.15Go deeper

O. Couture et al., "Ultrasound localization microscopy and super-resolution: a state of the art," IEEE Trans. Ultrason. Ferroelectr. Freq. Control (2018). doi:10.1109/TUFFC.2018.2850811. Review of the field's first years.
K. Christensen-Jeffries et al., "Super-resolution ultrasound imaging," Ultrasound Med. Biol. (2020). Multi-group review of methods and terminology.
B. Heiles, M. Correia, V. Hingot, M. Pernot, J. Provost, M. Tanter and O. Couture, "Ultrafast 3D ultrasound localization microscopy using a 32 × 32 matrix array," IEEE Trans. Med. Imaging (2019). doi:10.1109/TMI.2018.2890358.
J. A. Jensen et al., "Three-dimensional super-resolution imaging using a row–column array," IEEE Trans. Ultrason. Ferroelectr. Freq. Control (2020). doi:10.1109/TUFFC.2019.2948563.
J. Baranger et al., "Adaptive spatiotemporal SVD clutter filtering for ultrafast Doppler imaging using similarity of spatial singular vectors," IEEE Trans. Med. Imaging (2018). Automatic cut-off selection.
R. Parthasarathy, "Rapid, accurate particle tracking by calculation of radial symmetry centers," Nature Methods (2012). The radial-symmetry localizer.
R. P. J. Nieuwenhuizen et al., "Measuring image resolution in optical nanoscopy," Nature Methods (2013). doi:10.1038/nmeth.2448. FRC for localization microscopy.
M. van Heel and M. Schatz, "Fourier shell correlation threshold criteria," J. Struct. Biol. (2005). doi:10.1016/j.jsb.2005.05.009. Half-bit and σ thresholds.
J.-Y. Tinevez, simpletracker (MATLAB File Exchange / GitHub): the Hungarian and nearest-neighbour linkers used by PALA.
Lab cross-links: SVD clutter filtering (Module 06); delay estimation and subsample fits (Module 11); CRLB and ULM localization (Module 12); matched filtering (Module 10); compounding (Module 23); 2D and row–column arrays (Module 19); acoustic output (Module 22); microbubbles and the density trade-off (Module 28); Doppler and SVD in ultrafast imaging (Module 29); sound speed (Module 32) and aberration correction (Module 33), which decide whether transcranial and deep ULM maps are where they appear.

← 30 Elastography 32 Sound-speed estimation →