3D UltraSound Delay and Sum (DAS) and image reconstruction

Hi nitai,

Here's a simple 2D example demonstrating beamforming, in case it's helpful. The principles are identical in 3D.

Best wishes,
Ben

%--------------------------------------------------------------------------
% Simulate time series from a linear detection array
%--------------------------------------------------------------------------

% create the computational grid
PML_size = 20;          % size of the PML in grid points
Nx = 128 - 2*PML_size;  % number of grid points in the x (row) direction
Ny = 128 - 2*PML_size;  % number of grid points in the y (column) direction
dx = 0.05e-3;            % grid point spacing in the x direction [m]
dy = 0.05e-3;            % grid point spacing in the y direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500;           % [m/s]

% create the time array
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);

% define the sensor array
N_elements = 16;          % number of elements in the detector array
element_spacing = 3;      % spacing between elements [pixels]
sensor.mask = zeros(Nx, Ny);
sensor.mask((Nx/2 + 1 - N_elements*element_spacing/2) : element_spacing :(Nx/2 + N_elements*element_spacing/2), 1) = 1;

% create initial pressure distribution consisting of a few points
source_spacing_x = 20;    % x spacing between source points [pixels]
source_spacing_y = 30;    % y spacing between source points [pixels]
source.p0 = zeros(Nx,Ny);
source.p0(Nx/2+1 , Ny/2+1 - source_spacing_y) = 1;
source.p0(Nx/2+1 , Ny/2+1) = 1;
source.p0(Nx/2+1 , Ny/2+1 + source_spacing_y) = 1;
source.p0(Nx/2+1 - source_spacing_x, Ny/2+1) = 1;
source.p0(Nx/2+1 + source_spacing_x, Ny/2+1 - source_spacing_y) = 1;

% smooth the initial pressure distribution and restore the magnitude
source.p0 = smooth(kgrid, source.p0, true);

% set the input arguments: force the PML to be outside the computational
% grid; switch off p0 smoothing within kspaceFirstOrder2D
input_args = {'PMLInside', false, 'PMLSize', PML_size, 'Smooth', false, ...
    'PlotPML', false, 'PlotLayout', true, 'PlotScale', [-0.25 0.25]};

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

%--------------------------------------------------------------------------
% Form the A-lines
%--------------------------------------------------------------------------

% Reconstruction grid. Define the points for the A-line
% (ie. the distance from the central element to each pixel on the A-line)
pixel_distances = (0:0.25:Nx-1)*dy;      % [m]

% find the distances of the elements from the element at x=0
element_distances = kgrid.x_vec(sensor.mask(:,1)==1);   % [m]

% distance of single focus from array
focal_distance = (Ny/2+1)*dy;      % [m]

% Calculate which time points need to be summed together to synthesise
% focussing. For each pixel in the A-line, we need to choose one time point
% for each element. Uses Pythagoras' theorem.
A_line_onefocus = zeros(size(pixel_distances));
A_line_focussed = zeros(size(pixel_distances));
for loop = 1:N_elements

    % time points which give just one focal point at focal_distance
    time_points1 = ( sqrt( focal_distance.^2 + element_distances(loop).^2 ) - focal_distance + pixel_distances ) / medium.sound_speed;  

    % time points which focus all along the line x = 0
    time_points2 =  sqrt( pixel_distances.^2 + element_distances(loop).^2 ) / medium.sound_speed; 

    % convert time shifts to time indices (+1 as matlab indexing starts at 1 not 0)
    time_indices1 = round( time_points1 / dt ) + 1;
    time_indices2 = round( time_points2 / dt ) + 1;

    % sum the time series to make the A-lines
    A_line_onefocus = sensor_data(loop,time_indices1) + A_line_onefocus;
    A_line_focussed = sensor_data(loop,time_indices2) + A_line_focussed;

end

% SIMPLE SUM, A-line formed as sum, no time delays
A_line_sum = sum(sensor_data);

% NO FOCUSSING, A-line from just the sensor at x=0
A_line_unfocussed = sensor_data(N_elements/2,:);

figure
subplot(3,1,1)
plot(kgrid.y_vec,source.p0(Nx/2+1,:))
set(gca,'XLim',[min(kgrid.y_vec) max(kgrid.y_vec)])
title('sources on the x=0 axis')
grid on

subplot(3,1,2)
plot(medium.sound_speed*kgrid.t_array + min(kgrid.y_vec), A_line_unfocussed)
set(gca,'XLim',[min(kgrid.y_vec) max(kgrid.y_vec)])
title('A-line from single unfocussed element')
grid on

subplot(3,1,3)
plot(medium.sound_speed*kgrid.t_array + min(kgrid.y_vec), A_line_sum,'k')
hold on
plot(pixel_distances + min(kgrid.y_vec), A_line_onefocus,'r-')
plot(pixel_distances + min(kgrid.y_vec), A_line_focussed,'b--')
set(gca,'XLim',[min(kgrid.y_vec) max(kgrid.y_vec)])
title('A-lines with focussing')
legend('simple sum','one focus','multi-focus')
xlabel('x [m]')
grid on

k-Wave

A MATLAB toolbox for the time-domain
simulation of acoustic wave fields

3D UltraSound Delay and Sum (DAS) and image reconstruction

Reply

k-Wave

A MATLAB toolbox for the time-domain simulation of acoustic wave fields

3D UltraSound Delay and Sum (DAS) and image reconstruction

Reply

A MATLAB toolbox for the time-domain
simulation of acoustic wave fields