Field II and beamforming

Wednesday, May 24, 15.00-17.00 on the ground floor in build 349, rooms 19, 25 and 34.

Purpose:

This exercise will teach how to simulate ultrasound using Field-II, as well as how to beamform the images using synthetic aperture methods.

Preparation:

Download Field-II, along with the user-guide. Study the section/slides on beamforming in the course notes.

Exercise:

Delay-and-sum beamforming

You should make an m-file for each point in the exercise, so can run it for each of the points. This will gradually make your beamformer, and you can follow each step in the process. This can late be used for other exercises and the assignments in the course.

1. Make a sketch of the two middle elements of the 1D linear transducer array with the dimensions below. Indicate the pitch, kerf, width, and height in your figure.

c = 1540;            % Speed of sound [m/s].
f0 = 9e6;            % Signal center frequency [Hz].
lambda = c / f0;     % Wavelength [m]

pitch = lambda;      % Distance between transducer elements. [m]
kerf = pitch * 0.1;  % Empty space between transducer elements [m]
width = pitch * 0.9; % Width of elements [m]
height = width;      % Height of elements [m]

2. Compute the lateral position of each transducer element and store this in element_x. These positions should be symmetric around x = 0, and monotonically increasing.
   
   Hint: Use the colon operator, or function.

N_elements = 128; % Number of transducer elements.
element_x = ...; % 1xN_elements array of element lateral position [m]

3. It is very important that the positions are correct.
   Run the following code to verify your element positions.

% The following assertions will give hints if there was a mistake.
assert(element_x(1) < 0, ['The first element lateral position "element_x(1)" should be negative. (it is currently: ' num2str(element_x(1)) ')']);
assert(all(diff(element_x) > 0), 'Element lateral positions should increase monotonically');
assert(element_x(end/2) < 0, ['The ' num2str(N_elements/2) 'th element lateral position should be negative (it is currently: ' num2str(element_x(end/2)) ')']);
assert(element_x(end/2+1) > 0, ['The ' num2str(N_elements/2+1) 'th element lateral position should be positive (it is currently: ' num2str(element_x(end/2+1)) ')']);
assert(all(abs(element_x(1:end/2) + element_x(end:-1:end/2+1)) <= eps('single')), 'The element positions should show symmetry around x = 0.');
assert(all(abs(diff(element_x) - pitch) <= eps('single')), 'The spacing between element positions should be equal to the pitch.');
disp('Element positions OK!');  % No assertion failed :)

4. The following code plots your beamforming geometry. In this case with one scatterer, and one "active transmitting element".
   We will use this geometry going forward.

% Active emitting element (choose one).
active_element_idx = N_elements/2;

% Scatterer position (choose where you would like it to be)
scatterer_pos_x = 0;
scatterer_pos_z = 20/1000;

active_element_x = element_x(active_element_idx);  % 1x1 array of active element lateral position [m].

figure(1);
plot(active_element_x*1e3, active_element_x*0, 'ok', 'DisplayName', 'Transmit element');
hold on;
plot(scatterer_pos_x*1e3, scatterer_pos_z*1e3, 'xb', 'DisplayName', 'Scatterer');
plot(element_x*1e3, element_x*0, '.r', 'DisplayName', 'Receive elements');
hold off;
axis equal;
title('Aperture');
xlabel('x [mm]');
ylabel('z [mm]');
grid on;
legend;

5. Use the following code that simulates an emission from the active element.
   
   How do the received signal traces change when you modify the active element, or the scatterer position?

% Setup Field-II simulation
if ~exist('field_initialized', 'var')
    field_init(0);
    fs = f0 * 30;  % Sampling rate [Hz].
    set_sampling(fs);
    field_initialized = true;
end

% Aperture geometry.
Tx = xdc_linear_array (N_elements, width, height, kerf, 1, 1, [0, 0, 0]); % Transmit aperture.
Rx = xdc_linear_array (N_elements, width, height, kerf, 1, 1, [0, 0, 0]); % Receive aperture.
xdc_focus_times(Tx, 0, zeros(1, N_elements));  % Clear focusing delays, as we implement them manually.
xdc_focus_times(Rx, 0, zeros(1, N_elements));

% Set active element using Tx apodization mask.
element_mask = zeros(1, N_elements);
element_mask(active_element_idx) = 1;
xdc_apodization(Tx, 0, element_mask);

% A 4-cycle sinusoid emission (you may vary the number of cycles).
cycles = 4;
excitation_osci = sin((0:1/fs:cycles/f0)'*f0*pi*2);  % Pulse oscillation.  
excitation_envelope = hann(numel(excitation_osci)); % Pulse envelope.
excitation = excitation_osci .* excitation_envelope;
excitation = excitation - mean(excitation);
xdc_excitation(Tx, excitation');

% Set scatterer positions.
scatter_pos = [scatterer_pos_x(:), scatterer_pos_x(:)*0, scatterer_pos_z(:)];
scatter_amplitude = ones(size(scatter_pos, 1), 1);

% Use Field-II to simulate the received signals from the emission!
[rf_data, t0] = calc_scat_multi(Tx, Rx, scatter_pos, scatter_amplitude);
rf_data = rf_data ./ max(abs(rf_data(:))); % Rescale simulation output.

% Time axis for excitation.
t_axis0 = (0:1/fs:(numel(excitation)-1)/fs)';

% Move t = 0 to the middle of the excitation.
t_offset = mean(t_axis0);
t0 = t0 - t_offset;
t_axis = t0+(0:size(rf_data, 1)-1)'/fs;
t_axis0 = t_axis0 - t_offset;

figure(2);
subplot(211);
plot(t_axis0*1e6, excitation, 'DisplayName', 'Excitation');
grid on;
xlabel('time [\mus]');
ylabel('Amplitude [V]');
title('Transmitted signal');
legend;
subplot(212);
imagesc(1:N_elements, t_axis * 1e6, rf_data, [-1, 1]);
% imagesc(1:N_elements, t_axis * 1e6, abs(hilbert(rf_data)));
hold on;
tof = plot(nan, nan, '.-r', 'DisplayName', 'Predicted time of flight');
hold off;
colorbar;
xlabel('Element number');
ylabel('time [\mus]');
title('Received signal from each element');

6. Calculate the distance from the active element to the scatterer.

dist_tx = ...; % scalar [m]
% Tip: Use the variables active_element_x, scatterer_pos_x, scatterer_pos_z

7. Calculate the distance from each receiving element to the scatterer.

dist_rx = ...; % 1xN_elements [m]
% Tip: Use the variables element_x, scatterer_pos_x, scatterer_pos_z

8. Time to verify your calculations are correct!
   Use the following code to plot the predicted time of flight on top of the received signal traces (previous figure).
   They should be aligned perfectly.

% Run the following assertion first. It gives a hint if there was a mistake.
assert(abs(dist_tx-dist_rx(active_element_idx)) <= eps('single'), ...
'"dist_tx" should be equal to "dist_rx(active_element_idx)"');

% Compute the delays.
distance = dist_tx+dist_rx;  % distance [m]
delays = distance/c;  % time of flight [s]

figure(2);
subplot(212);
tof.XData = 1:N_elements;
tof.YData = delays * 1e6;  % Plot predicted time of flight in microseconds.
legend;

9. Now to the delay-and-sum beamformation. Data is delayed using your computed delays, which brings each trace into alignment.

rf_data_delayed = zeros(numel(t_axis0), size(rf_data, 2));
% Delay and sum
for i = 1:numel(delays)
    % Perform the delay.
    rf_data_delayed(:, i) = interp1(t_axis, rf_data(:, i), delays(i)+t_axis0, 'spline', 0);
    % rf_data_delayed(:, i) = rf_data_delayed(:, i) .* dist_rx(i) .* dist_tx;
end
% Perform the sum.
DAS_result = sum(rf_data_delayed, 2);
DAS_result = DAS_result / max(abs(DAS_result)) * max(abs(excitation));  % Scale to match first plot.

rf_data_delayed = rf_data_delayed / max(abs(rf_data_delayed(:)));

figure(3);
subplot(211);

imagesc(1:N_elements, t_axis0*1e6, rf_data_delayed, [-1, 1]);
title('Received signals after applying delays.');
colorbar;
xlabel('Element number');
ylabel('time w.r.t. delay [\mus]');

hold on;
plot(1:N_elements, delays*0, '.-r', 'DisplayName', 'Predicted time of flight');
hold off;

subplot(212);
plot(t_axis0*1e6, excitation, '--', 'DisplayName', 'Excitation');
hold on;
plot(t_axis0*1e6, DAS_result, 'DisplayName', 'Delay-and-sum result');
hold off;
title('Beamforming (Sum of signals after applying delays)');
ylabel('Amplitude [$V/V$]', 'Interpreter', 'latex');
xlabel('time w.r.t. delay [\mus]');
grid on;
legend;

There may be small time-quantization artifacts, but these are usually suppressed by increasing the simulation sampling rate "fs".


10. Next step is envelope detection. Notice how the pulse amplitude in the center is 0, but this is actually the exact position of the scatterer!
    We would like to have a peak at the scatterer position to get a more useful image.
    The function "hilbert" computes the analytic (complex) representation of our data, after which the envelope can be found as the absolute value of this.

figure(4);
DAS_result_analytic = hilbert(DAS_result);

plot(t_axis0*1e6, real(DAS_result_analytic), 'Display', 'real(analytic(DAS)) = DAS');
hold on;
plot(t_axis0*1e6, imag(DAS_result_analytic), 'Display', 'imag(analytic(DAS))');
plot(t_axis0*1e6, abs(DAS_result_analytic), 'Display', 'abs(analytic(DAS)) = detected envelope');
plot(t_axis0*1e6, excitation_envelope, '--', 'Display', 'True envelope');
hold off;
ylabel('Amplitude [$V/V$]', 'Interpreter', 'latex');
xlabel('time w.r.t. delay [\mus]');
legend;
grid on;
title('Envelope detection');

Interestingly, the hilbert transformation can be applied before beamforming the RF data (as below), or it can be applied after (as in the example above).


11. Time for beamforming an image.
    With envelope detection, we only need to gather the central value (t=0 after applying delay), as that is a peak. But in practice we do not know where exactly the scatterer is. Therefore we need to calculate "dist_rx", "dist_tx" for any pixel as-if the scatterer is in this pixel. Insert your expressions for these distances in the code below.

apodization = ones(1, N_elements);

% Define image dimensions.
h = 40/1000;  % [m]
w = (element_x(end)-element_x(1))*1.5;
Nz = ceil(8*h/lambda);  % Sample at a proper rate.
Nx = ceil(4*w/lambda);
x = linspace(-w/2, w/2, Nx);  % Lateral positions.
z = linspace(0, h, Nz)';  % Axial positions.

% Prepare animated figure.
figure(6);
sc = imagesc(x*1e3, z*1e3, 0);
hold on;
plot(element_x(active_element_idx)*1e3, element_x(active_element_idx)*0, 'ok', 'DisplayName', 'Transmit element');
plot(scatterer_pos_x(:)*1e3, scatterer_pos_z(:)*1e3, 'xb', 'DisplayName', 'Scatterer');
plot(element_x(:)*1e3, element_x(:)*0, '.r', 'DisplayName', 'Receive elements');
hold off;
legend;
axis image;
colormap gray;
xlim(scatterer_pos_x*1e3 + [-2, 2]);  % Zoom on scatterer (disable if you want).
ylim(scatterer_pos_z*1e3 + [-2, 2]);
xlabel('x [mm]');
ylabel('z [mm]');

% Delay and sum beamformer.
rf_data_analytic = hilbert(rf_data);
beamformed_image = 0;
for i = 1:size(rf_data, 2)
    tx_pos_x = element_x(active_element_idx);
    rx_pos_x = element_x(i);
    
    % Compute the distance from tx_pos to (x, z).
    dist_tx = ...;
    % Compute the distance from rx_pos to (x, z).
    dist_rx = ...;
    % Hint:
    %  The two arrays should be in meters, and the size should be: [Nz, Nx].
    %  Use the same calculation as in parts 6 and 7 of the exercises, 
    %   but replace "scatterer_pos_x" with "x", and "active_element_x" with "tx_pos_x" and so on.
    %  Computations should be vectorized for efficiency (don't use loops!).
    %   This means you need to use element-wise operations, e.g. .^ instead of ^.
    delays = (dist_tx + dist_rx) * 1/c;
    
    backprojection = interp1(t_axis, rf_data_analytic(:, i), delays, 'spline', 0);
    beamformed_image = beamformed_image + apodization(i) .* backprojection;
    
    % Update plot
    sc.CData = abs(beamformed_image);
    title(['Delay and sum beamforming for element 1 to ' num2str(i)]);
    drawnow;
end
title(['Delay and sum beamforming result']);

12. The B-mode image is the log-compressed version of the envelope image. These are used because ultrasound has a high dynamic range.
    Compute the B-mode image.
    
    

B_mode = ...;
% Hint: Use abs(beamformed_image).
%       Log compression means calculating 20*log10 of the values, and subtracting the peak value.

13. Use the following code to display your B-mode image.

figure(7);
imagesc(x*1e3, z*1e3, B_mode, [-60, 0]);
hold on;
plot(element_x(active_element_idx)*1e3, element_x(active_element_idx)*0, 'ok', 'DisplayName', 'Transmit element');
plot(scatterer_pos_x(:)*1e3, scatterer_pos_z(:)*1e3, 'xb', 'DisplayName', 'Scatterer');
plot(element_x(:)*1e3, element_x(:)*0, '.r', 'DisplayName', 'Receive elements');
hold off;
legend;
axis image;
colormap gray;
xlim(scatterer_pos_x*1e3 + [-2, 2]);
ylim(scatterer_pos_z*1e3 + [-2, 2]);
xlabel('x [mm]');
ylabel('z [mm]');
title('B-mode image [dB]');
colorbar;

Now that the image is more detailed, you will notice there are large "wings" on your scatterer. These are called sidelobes, and they ruin your image contrast and resolution.


14. Apodization is used to suppress sidelobes. The reason for the high sidelobes is essentially that the receiving aperture suddenly stops, which generates a sinc-like beam shape (an approximation says the beam is the FFT of the aperture). The solution is to use [https://en.wikipedia.org/wiki/Window_function windowing].
    Modify the apodization in exercise 11 to beamform using your desired window.

% Some examples (choose one, and re-run exercise 11 and 12 with your choice):
apodization = hanning(N_elements);  % Hann window. Easy to compute. A very common choice.
apodization = hamming(N_elements);  % Similar to the hann window, but will suppress the first side-lobe better.
apodization = blackmanharris(N_elements); % Larger main-lobe, but much smaller side-lobes overall.

The windows can be analyzed with "wvtool(apodization)". See window functions [https://se.mathworks.com/help/signal/windows.html here]. There is a trade-off between side-lobe level (contrast), and main-lobe width (resolution).
Essentially, choosing a window is an art form, but almost any choice is better than no apodization. (Even though no apodization gives the smallest possible main lobe).


15. This is it. Now try beamforming a more complex image. Use the following to generate a point phantom, and the resulting RF data.

scatterer_pos_x = 0.5*linspace(-15, 15, 5)/1000;
scatterer_pos_z = 0.5*linspace(10, 40, 5)'/1000;

% Convert the two lines to their cartesian product grid.
scatterer_pos_x = scatterer_pos_x + 0 * scatterer_pos_z;
scatterer_pos_z = 0 * scatterer_pos_x + scatterer_pos_z;

% Set scatterer positions.
scatter_pos = [scatterer_pos_x(:), scatterer_pos_x(:)*0, scatterer_pos_z(:)];
scatter_amplitude = ones(size(scatter_pos, 1), 1);
[rf_data, t0] = calc_scat_multi(Tx, Rx, scatter_pos, scatter_amplitude);
rf_data = rf_data ./ max(abs(rf_data(:))); % Rescale simulation output.

% Time axis for excitation.
t_axis0 = (0:(numel(excitation)-1))'/fs;
% Move t = 0 to the middle of the excitation.
t0 = t0 - mean(t_axis0);
t_axis = t0+(0:size(rf_data, 1)-1)'/fs;
t_axis0 = t_axis0 - mean(t_axis0);

figure(8);
imagesc(1:N_elements, t_axis * 1e6, abs(hilbert(rf_data)));
colorbar;
xlabel('Element number');
ylabel('time [\mus]');
title('Received signal envelope from each element.');

16. Now try beamforming the image with these scatterers.

apodization = ...;  % <- Insert your apodization

% Delay and sum beamformer.
w = diff(scatterer_pos_x([1, end]))*1.2;
h = scatterer_pos_z(end)*1.2;

Nz = ceil(8*h/lambda);
Nx = ceil(4*w/lambda);
x = linspace(-w/2, w/2, Nx);  % Lateral positions.
z = linspace(0, h, Nz)';  % Axial positions.

rf_data_analytic = hilbert(rf_data);

figure(9);
sc = imagesc(x*1e3, z*1e3, 0);
hold on;
plot(element_x(active_element_idx)*1e3, element_x(active_element_idx)*0, 'ok', 'DisplayName', 'Transmit element');
plot(scatterer_pos_x(:)*1e3, scatterer_pos_z(:)*1e3, 'xb', 'DisplayName', 'Scatterer');
plot(element_x(:)*1e3, element_x(:)*0, '.r', 'DisplayName', 'Receive elements');
hold off;
xlabel('x [mm]');
ylabel('z [mm]');
axis image; axis xy;
colormap gray;
legend;
drawnow;

beamformed_image = 0;
for i = 1:size(rf_data, 2)
    tx_pos_x = element_x(active_element_idx);
    rx_pos_x = element_x(i);
    
    % Compute the distance from tx_pos to (x, z).
    dist_tx = ...;  % <- Insert your expression
    % Compute the distance from rx_pos to (x, z).
    dist_rx = ...;  % <- Insert your expression

    delays = (dist_tx + dist_rx) * 1/c;
    backprojection = interp1(t_axis, rf_data_analytic(:, i), delays, 'spline', 0);
    beamformed_image = beamformed_image + apodization(i) .* backprojection .* dist_rx .* dist_tx;
    if mod(i-1, 1) == 0
        sc.CData = abs(beamformed_image);
        drawnow;
    end
end
colorbar;
title('Press any key to continue to B-mode display.');
pause
B_mode = ...;
sc.CData = B_mode;
caxis([-60, 0]);
title('B-mode of wire phantom.');

Note we scaled the backprojection by the distance traveled by the transmitted wave, and that of the received (scattered) wave to keep all scatterers at the same amplitude
Try e.g. to plot the center line with plot(z*1e3, B_mode(:, (end-1)/2)); xlabel('z [mm]'); ylabel('Amplitude [dB]'); ylim([-60, 0]);.
You should also notice that scatterers on the side of the image appear diagonal/rotated. Finally, the size of the scatterer increases with depth. We usually solve these three issues by using dynamic apodization.


17. Dynamic apodization part 1:
    Move the aperture window with the point being beamformed - this removes the rotated point effect, since then the receiving aperture is always below the point being beamformed!
    Place the following code in your beamforming loop from exercise 16. What do you observe?

dynamic_apodization = interp1(element_x, apodization, rx_pos_x-x, 'linear', 0);
beamformed_image = beamformed_image + dynamic_apodization .* backprojection .* dist_rx .* dist_tx;

Of course, this dynamic aperture cannot move far outside the physical aperture, and still give meaningful results.
You can visualize the dynamic aperture in the main loop with sc.CData = dynamic_apodization;


18. Dynamic apodization part 2:
    Grow the aperture with depth - this gives us constant resolution, as the resolution depends on the F-number: aperture width / focal depth.
    It also gives us the "dist_rx" gain compensation term, because now we only include receivers directly below the pixel, within a certain angular span (more depth = more receivers within the angle = more signal, less noise). In other words, we get (approximately) uniform resolution and signal-to-noise ratio with dynamic apodization.
    Place the following code in your beamforming loop from exercise 16. What do you observe when you vary F_num_rx?

F_num_rx = 0.5;  % You get to choose the F-number, thus resolution!
% Since F-number is width divided by depth, simply grow width with depth to maintain your desired F#:
dynamic_apodization = interp1(linspace(-0.5, 0.5, numel(apodization)), apodization, F_num_rx*(rx_pos_x-x)./z, 'linear', 0);

beamformed_image = beamformed_image + dynamic_apodization .* backprojection .* dist_tx;

Of course the dynamic virtual aperture cannot grow larger than the physical aperture. The uniform resolution can only be maintained till a certain depth. The Tx apodization is used to reject grating lobes.


19. The final exercise is to use more than one emission. Play around with the apodization and virtual sources to see how good an image you can get.

% Adjust the beamforming parameters below how you like.
F_num_rx = 0.7;
F_num_tx = 0.7;
% The elements you want to emit with (change how you like)
active_element_idx = round(linspace(1, N_elements, 12));
% active_element_idx = [1, N_elements/4, N_elements/2, 3/4*N_elements, N_elements];
apodization_rx = apodization; % <- Your window of choice here.
apodization_tx = apodization;

% Delay and sum beamformer.
w = (max(scatterer_pos_x(:))-min(scatterer_pos_x(:)))*1.2;
h = max(scatterer_pos_z(:))*1.2;

S = 2;  % Sparse beamforming, to speed up calculations.
Nz = ceil(8/S*h/lambda);
Nx = ceil(4/S*w/lambda);
x = linspace(-w/2, w/2, Nx);  % Lateral positions.
z = linspace(0, h, Nz)';  % Axial positions.

rf_data_analytic = hilbert(rf_data);

figure(8);
sc = imagesc(x*1e3, z*1e3, 0);
hold on;
plot(element_x(active_element_idx)*1e3, element_x(active_element_idx)*0, 'ok', 'DisplayName', 'Transmit elements');
plot(scatterer_pos_x(:)*1e3, scatterer_pos_z(:)*1e3, 'xb', 'DisplayName', 'Scatterer');
plot(element_x(:)*1e3, element_x(:)*0, '.r', 'DisplayName', 'Receive elements');
hold off;
axis image; axis xy;
colormap gray;
legend;
drawnow;

beamformed_image = 0;

% Simulate and beamform multiple emissions.
for j = 1:numel(active_element_idx)
    tx_pos_x = element_x(active_element_idx(j));
    dynamic_apodization_tx = interp1(linspace(-0.5, 0.5, numel(apodization_tx)), apodization_tx, F_num_tx*(tx_pos_x-x)./z, 'linear', 0);
    
    % Set active element
    element_mask = zeros(1, N_elements);
    element_mask(active_element_idx(j)) = 1;
    xdc_apodization(Tx, 0, element_mask);

    % Simulate the emission.
    [rf_data, t0] = calc_scat_multi(Tx, Rx, scatter_pos, scatter_amplitude);
    t0 = t0 - (numel(excitation)-1)/fs/2;
    t_axis = t0+(0:size(rf_data, 1)-1)'/fs;
    rf_data_analytic = hilbert(rf_data);

    % Compute the distance from tx_pos to (x, z).
    dist_tx = ...; % <- Your code here.
    
    for i = 1:size(rf_data, 2)
        rx_pos_x = element_x(i);

        
        % Compute the distance from rx_pos to (x, z).
        dist_rx = ...; % <- Your code here.

        delays = (dist_tx + dist_rx) * 1/c;
        backprojection = interp1(t_axis, rf_data_analytic(:, i), delays, 'spline', 0);
        dynamic_apodization_rx = interp1(linspace(-0.5, 0.5, numel(apodization_rx)), apodization_rx, F_num_rx*(rx_pos_x-x)./z, 'linear', 0);
        beamformed_image = beamformed_image + dynamic_apodization_rx .* dynamic_apodization_tx .* backprojection;
        if mod(i-1, 10) == 0
            sc.CData = abs(beamformed_image);
            title(['Beamforming Tx: ' num2str(j) '/' num2str(numel(active_element_idx)) ', Rx: ' num2str(i) '/' num2str(N_elements)]);
            drawnow;
        end
    end
end
B_mode = ...; % <- Your code here.
sc.CData = B_mode;
caxis([-60, 0]);
colorbar;
title('B-mode synthetic aperture imaging');

With dynamic transmit and receive apodization we get approximately a shift-invariant aperture, meaning image properties such as resolution and SNR (signal-to-noise ratio) should be the same across the image. In practice, there is one additional gain compensation term to counteract the attenuation that is observed, but this depends on what is being imaged and therefore cannot be known in advance. Usually the probe operator sets a gain curve for this. In practice, single-element emissions give too weak an echo due to low energy. We usually use a subset (sub-aperture) of elements to synthesize a virtual audio source, which is used as-if it is an active element. Conversely, if the whole probe is used in emission (e.g. in high frame-rate plane-wave imaging), you easily run into issues with the probe heating up. You can use the following code to see the uniformity of the lateral resolution (FWHM, width at half of max = -6dB)

normalize = true;  % Normalize to show resolution, not ampltiude uniformity.
depths = unique(scatterer_pos_z);
laterals = unique(scatterer_pos_x);
xline_axis = linspace(-1, 1, 401)/1000;
for i = 1:numel(depths)
    x_line = B_mode(find(z>=depths(i), 1), :);

    for j = 1:numel(laterals)
        x_line_sub = interp1(x, x_line, laterals(j) + xline_axis, 'spline', 0);
        if normalize
            x_line_sub = x_line_sub - max(x_line_sub(:));
        end
        subplot(100+numel(laterals)*10+j);
        plot(xline_axis*1e3, x_line_sub, 'DisplayName', ['z = ' num2str(depths(i)*1e3) 'mm']);
        xlabel('x offset [mm]');
        title(['PSF at x = ' num2str(laterals(j)*1e3) 'mm']);
        ylim([20*log10(1/2), 0]);
        hold on;
    end
end
subplot(100+numel(laterals)*10+1);
ylabel('Amplitude [dB]');
for j = 1:numel(laterals)
    subplot(100+numel(laterals)*10+j); hold off;
end
legend('Location', 'South');

20. Bonus: Tissue consists of millions of scatters. Try and beamform using the phantom definition below for a more complicated image. It may take some minutes.

scatterer_pos_x = 10/1000 * (rand(20000, 1)-0.5);
scatterer_pos_z = (5+20 * rand(20000, 1))/1000;

scatter_pos = [scatterer_pos_x(:), scatterer_pos_x(:)*0, scatterer_pos_z(:)];
scatter_amplitude = ones(size(scatter_pos, 1), 1);
scatter_amplitude((scatterer_pos_z-15/1000).^2 + scatterer_pos_x.^2 <= (3/1000).^2) = 0;
scatter_amplitude(1) = 10;

% Remember to disable the scatter-plotting of scatterers!