`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 interro- gation 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: 6.5 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