SPTRACK Speckle tracking
SPTRACK returns the motion field by speckle tracking.
Contents
Syntax
[Di,Dj] = SPTRACK(I,PARAM) returns the motion field [Di,Dj] that occurs from frame #k I(:,:,k) to frame #(k+1) I(:,:,k+1).
I must be a 3-D array, with I(:,:,k) corresponding to image #k. I can contain more than two images (i.e. size(I,3)>2). In such a case, an ensemble correlation is used.
TRY IT! Enter sptrack in the command window for an example.
Di,Dj are the displacements (unit = pix) in the IMAGE coordinate system (i.e. the "matrix" axes mode). The i-axis is vertical, with values increasing from top to bottom. The j-axis is horizontal with values increasing from left to right. The coordinate (1,1) corresponds to the center of the upper left pixel. To display the displacement field, you may use: quiver(Dj,Di), axis ij
PARAM is a structure that contains the parameter values required for speckle tracking (see below for details).
[Di,Dj,id,jd] = SPTRACK(...) also returns the coordinates of the points where the components of the displacement field are estimated.
PARAM is a structure that must contain the following fields:
- PARAM.winsize: size of the interrogation windows (required)
PARAM.winsize must be a 2-COLUMN array. If PARAM.winsize contains several rows, a multi-grid, multiple-pass interrogation process is used.
Examples:
a) If PARAM.winsize = [64 32], then a 64 32 (64 lines, 32 columns) interrogation window is used.
b) If PARAM.winsize = [64 64;32 32;16 16], a 64 64 interrogation window is first used. Then a 32
32 window and finally a 16
16 window are used.
- PARAM.overlap: overlap between the interrogation windows (in %, default = 50)
- PARAM.iminc: image increment (for ensemble correlation, default = 1)
The image #k is compared with image #(k+PARAM.iminc): I(:,:,k) is compared with I(:,:,k+PARAM.iminc)
- PARAM.ROI: 2-D region of interest (default = the whole image)
PARAM.ROI must be a logical 2-D array with a size of [size(I,1),size(I,2)]. The default is all(isfinite(I),3).
Notes
- The displacement field is returned in PIXELS. Perform an appropriate calibration to get physical units.
- SPTRACK is based on a multi-step cross-correlation method. The SMOOTHN function (see Reference below) is used at each iterative step for the validation and post-processing.
Example #1: Speckle tracking of a rotating image
This example shows how to obtain the motion field of a rotating image
Create a 500-by-500 image.
I1 = conv2(rand(500,500),ones(10,10),'same');
Rotate it in the clockwise direction.
I2 = imrotate(I1,-3,'bicubic','crop'); subplot(121), imagesc(I1), axis square off title('1^{st} frame') subplot(122), imagesc(I2), axis square off title('2^{nd} frame') colormap gray

Recover the motion field with SPTRACK.
param = []; param.winsize = [64 64;32 32]; [Di,Dj,id,jd] = sptrack(cat(3,I1,I2),param);
Display the motion field.
figure plot(jd(1:2:end,1:2:end),id(1:2:end,1:2:end),'ko',... 'MarkerFaceColor','k','MarkerSize',4) hold on h = quiver(jd(1:2:end,1:2:end),id(1:2:end,1:2:end),... Dj(1:2:end,1:2:end),Di(1:2:end,1:2:end),2); hold off set(h,'LineWidth',1.5) axis equal ij tight title('This is the motion field')

Example #2: Speckle tracking with plane wave imaging
This example shows how to obtain the motion field of a rotating disk insonified with plane waves.
A rotating disk (diameter of 2 cm) was insonified by a series of 32 unsteered plane waves with a Verasonics scanner, and a linear transducer, at a PRF (pulse repetition frequency) of 10 kHz. The RF signals were downsampled at 4/3 (5 MHz) = 6.66 MHz. The properties of the linear array were:
- 128 elements
- center frequency = 5 MHz
- pitch = 0.298 mm
Download the experimental RF data. The 3-D array RF contains 128 columns (as the transducer contained 128 elements), and its length is 32 in the third dimension (as 32 plane waves were transmitted).
load('PWI_disk.mat')
The structure param contains some experimental properties that are required for the following steps.
disp('''param'' is a structure whose fields are:')
disp(param)
'param' is a structure whose fields are: fs: 6.6667e+06 pitch: 2.9800e-04 fc: 5000000 c: 1480 t0: 9.9500e-06 PRF: 10000 width: 2.6200e-04 Nelements: 128 TXdelay: [1×128 double] bandwidth: 15
Demodulate the RF signals with RF2IQ.
IQ = rf2iq(RF,param);
Create a 2.5-cm-by-2.5-cm image grid.
dx = 1e-4; % grid x-step (in m) dz = 1e-4; % grid z-step (in m) [x,z] = meshgrid(-1.25e-2:dx:1.25e-2,1e-2:dz:3.5e-2);
Create a Delay-And-Sum DAS matrix with DASMTX.
param.fnumber = []; % an f-number will be determined by DASMTX M = dasmtx(1i*size(IQ),x,z,param,'nearest'); spy(M) axis square title('The sparse DAS matrix')

Beamform the I/Q signals.
The 32 I/Q series can be beamformed simultaneously with the DAS matrix.
tic IQb = M*reshape(IQ,[],32); IQb = reshape(IQb,[size(x) 32]); disp(['Beamforming: time per frame: ' num2str(toc/32*1e3,'%.1f'), ' ms'])
Beamforming: time per frame: 29.8 ms
Create the B-mode images with BMODE and display the first ultrasound image.
I = bmode(IQb,30); % B-mode images image(x(1,:)*100,z(:,1)*100,I(:,:,1)) c = colorbar; c.YTick = [0 255]; c.YTickLabel = {'-30 dB','0 dB'}; colormap gray title('The 1^{st} B-mode image') ylabel('[cm]') axis equal tight ij set(gca,'XColor','none','box','off')

Create an ROI.
param.ROI = median(I,3)>64;
Track the speckles with SPTRACK.
param.winsize = [32 32; 24 24; 16 16]; % size of the subwindows param.iminc = 4; % image increment [Di,Dj,id,jd] = sptrack(I,param);
Display the motion field.
image(I(:,:,1)) colormap gray hold on h = quiver(jd,id,Dj,Di,3,'r'); set(h,'LineWidth',1) hold off title('Motion field (in pix) by speckle tracking') axis equal off ij

See also
References
- Garcia D, Lantelme P, Saloux É. Introduction to speckle tracking in cardiac ultrasound imaging. Handbook of speckle filtering and tracking in cardiovascular ultrasound imaging and video. Institution of Engineering and Technology, 2018. (PDF)
- Perrot V, Garcia D. Back to basics in ultrasound velocimetry tracking speckles by using a standard PIV algorithm. IEEE International Ultrasonics Symposium (IUS), 2018. (PDF)
- Garcia D. A fast all-in-one method for automated post-processing of PIV data. Experiments in Fluids, 2011. (PDF)
About the author
Damien Garcia, Eng., Ph.D. INSERM researcher Creatis, University of Lyon, France
websites: BioméCardio, MUST
Date modified