FLUST example: Simulation of a spinning disk phantom

In this example, FLUST is used to simulate IQ Doppler signals from a spinning disk phantom, using unfocused ultrasound transmissions. This phantom represents a special case where the start and end of the flowlines are both at x = 0 in the lower part of the phantom, creating a discontinuity at this border.


% Updated 27/02/2023,
% Joergen Avdal (jorgen.avdal@ntnu.no),
% Ingvild Kinn Ekroll (ingvild.k.ekroll@ntnu.no)

% If using FLUST for scientific publications, please cite the original paper
% Avdal et al: Fast Flow-Line-Based Analysis of Ultrasound Spectral and
% Vector Velocity Estimators, tUFFC, 2019. DOI: 10.1109/TUFFC.2018.2887398

% FLUST is a simulation tool based on flowlines, useful for producing many
% realizations of the same flowfield. The motivation for making FLUST was
% to faciliate accurate assessment of statistical properties of velocity
% estimators like bias and variance.

% Using FLUST requires the user to provide/select two functions, one to create
% the flowlines, and one to calculate the PSFs from a vector of spatial
% positions.

% After running this script, variable realTab contains realizations, and
% variables PSFs and PSFstruct contains point spread functions of the last flow line.

% FLUST accounts for interframe motion in plane wave sequences, assuming
% uniform firing rate, but not in scanned sequences.

% How to use FLUST:
% 1) Provide/select function to calculate PSFs from a vector of spatial positions.
% 2) Run FLUST on phantom of interest.
% 3) Apply your favorite velocity estimator to realizations.
% 4) Assess statistical properties of estimator, optimize estimator.
% 5) Publish results, report statistical properties, make results
%    reproducible, cite original paper DOI: 10.1109/TUFFC.2018.2887398

clear all;
close all;



First, we define the firing rate of the simulated signal, the number of realizations and the ensemble size, and whether or not this is a simulation of contrast bubbles or not:

s = struct();

s.firing_rate = 12000; % firing rate of output signal, (Doppler PRF) = (firing rate)/(nr of firings)
s.nrReps = 10;         % nr of realizations
s.nrSamps = 40;       % nr of slow time samples in each realization (Ensemble size)

s.contrastMode = 0;      % is set to 1, will simulate contrast scatterers propagating in flow field
s.contrastDensity = 0.1; % if using contrastMode, determines the density of scatterers, typically < 0.2


Then, a recommended distance for PSF function seedpoints is calculated based on a user defined accepted value for interpolation error in percent:

s.interpErrorLimit = 4; % FLUST will set s.dr to attain interpolation error smaller than this value in percent


Choose whether or not to utilize your GPU, and how many scanlines/channels to process simultaneously when generating realizations from calculated PSFs. Optimal values depend primarily on GPU specifications and available memory.

s.chunksize = 16; % chunking on scanlines, adjust according to available memory.
s.useGPU = 1;


Set desired PSF function and acquisition parameters. If not set, predefined values are used (see individual PSF function definitions). Chunksize here refers to the number of simultaneously processed PSFs.

s.PSF_function = @PSFfunc_LinearProbe_PlaneWaveImaging_txApod;

% Transducer and acquisition parameters. Print s.PSF_params after running simulation to see which parameters can be set.
s.PSF_params = [];

% Transducer params
s.PSF_params.trans.f0 = 12.0e6;
s.PSF_params.trans.pulse_duration = 1.5;
s.PSF_params.trans.element_height = 3e-3;
s.PSF_params.trans.lens_el = 0.01;
s.PSF_params.trans.pitch = 100e-6;
s.PSF_params.trans.kerf = 10e-6;
s.PSF_params.trans.N = 256;

% Acquisition params
s.PSF_params.acq.alphaTx = [0]*pi/180;
s.PSF_params.acq.alphaRx = [0]*pi/180;
s.PSF_params.acq.F_number = 1.0;
s.PSF_params.acq.tx_apod = 'none';

% Image/scan region params
s.PSF_params.scan.rx_apod = 'blackman';
s.PSF_params.scan.xStart = -5e-3;
s.PSF_params.scan.xEnd = 5e-3;
s.PSF_params.scan.Nx = 128;
s.PSF_params.scan.zStart = 5e-3;
s.PSF_params.scan.zEnd = 15e-3;
s.PSF_params.scan.Nz = 128;

% Runtime params
s.PSF_params.run.chunkSize = 125;


Set desired phantom and phantom parameters. If not set, predefined values are used (see individual phantom function definitions).

s.phantom_function = @Phantom_spinningDisk;

% Phantom parameters. Print s.phantom_params after running simulation to see which parameters can be set.
s.phantom_params = [];
s.phantom_params.maxVel = 0.7;
s.phantom_params.minVel = 0.2;
s.phantom_params.diameter = 0.007; %m
s.phantom_params.depth = 0.01;

% To output true velocities in phantom, define grid
myX = linspace(s.PSF_params.scan.xStart,s.PSF_params.scan.xEnd,s.PSF_params.scan.Nx);
myZ = linspace(s.PSF_params.scan.zStart,s.PSF_params.scan.zEnd,s.PSF_params.scan.Nz);
myY = 0;
[X,Z,Y] =  meshgrid(myX,myZ,myY);

% make phantom, get true velocity and phantom parameters
[flowField, s.phantom_params, GT] = s.phantom_function(s.phantom_params,X,Y,Z); % flowField should have timetab and postab fields

FLUST main loop

Run your FLUST simulation!

Suggested upper limit for s.dr for axial interpolation 3e-05
Suggested upper limit for s.dr for lateral interpolation 9e-05
s.dr changed to 3e-05

VISUALIZE FIRST REALIZATION using the built-in beamformed data object

firstRealization = realTab(:,:,:,1,1);

b_data = uff.beamformed_data();
b_data.scan = PSFstruct.scan;
b_data.data = reshape(firstRealization,size(firstRealization,1)*size(firstRealization,2),1,1,size(firstRealization,3));
b_data.frame_rate = 20;
b_data.plot([],'Flow from FLUST', 20);

True velocities

The true phantom velocities are available for visualization. Future versions of this example will provide a demonstration of how to use the built in velocity estimator performance assessment tools in FLUST.

 % GT_rsh = reshape( GT, [s.PSF_params.scan.Nz s.PSF_params.scan.Nx 1 3] ); 
 % figure(2); imagesc(X(:), Z(:), GT_rsh(:,:,1,1)); axis equal; 
 % title('Vx') % Example, looking at x component of velocity field 
 % figure(3); imagesc(X(:), Z(:), GT_rsh(:,:,1,3)); axis equal; 
 % title('Vz') % Example, looking at z component of velocity field 

SAVE realizations and metadata

 % tag = 'myTag'; 
 % saveBinary( tag);