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 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 $\times$ 32 (64 lines, 32 columns) interrogation window is used.

b) If PARAM.winsize = [64 64;32 32;16 16], a 64 $\times$ 64 interrogation window is first used. Then a 32 $\times$ 32 window and finally a 16 $\times$ 16 window are used.



The image #k is compared with image #(k+PARAM.iminc): I(:,:,k) is compared with I(:,:,k+PARAM.iminc)


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

  1. The displacement field is returned in PIXELS. Perform an appropriate calibration to get physical units.
  1. 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 ${\times}$ (5 MHz) = 6.66 MHz. The properties of the linear array were:

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

bmode, iq2doppler, dasmtx



References



About the author

Damien Garcia, Eng., Ph.D.
INSERM researcher
Creatis, University of Lyon, France

websites: BioméCardio, MUST



Date modified


Published with MATLAB® R2020a