I have been using k wave in MATLAB to record the intensity patterns of some ultrasound transductors, at different planes of the 3D simulation grid. I have been defining a sensor for each plane, and then running the simulation with the kspaceFirstOrder3D function. But this requires more time as I define more sensors.
In order to reduce computing times, I would like to combine both sensors into a single mask, and then separating each dataset to the corresponding sensor coordinates
This is an extract of interest of my code:
$ Sensor 1 mask
sensor1.mask = zeros(Nx, Ny, Nz);
sensor1.mask(:, :, Nz/2 + 1) = 1;
sensor1.record = {'p', 'p_max'};
sensor1_data = kspaceFirstOrder3D(kgrid, medium, source, sensor1, input_args{:});
$ Sensor 2 mask
sensor2.mask = zeros(Nx, Ny, Nz);
sensor2.mask(Nx/2 + 1, :, :) = 1;
sensor2.record = {'p', 'p_max'};
sensor2_data = kspaceFirstOrder3D(kgrid, medium, source, sensor2, input_args{:});
$ Displaying intensity pattern for sensor 1
sensor1_data.p_max = reshape(sensor1_data.p_max, [Ny, Nx]);
figure('Name', 'YX Plane');
imagesc(y_axis, x_axis, sensor1_data.p_max);
$ Displaying intensity pattern for sensor 2
sensor2_data.p_max = reshape(sensor2_data.p_max, [Nz, Ny]);
figure('Name', 'YZ Plane');
imagesc(y_axis, z_axis, sensor2_data.p_max);
Where the axis are defined as
x_axis = (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3;
y_axis = (kgrid.y_vec - min(kgrid.y_vec(:))) * 1e3;
z_axis = (kgrid.z_vec - min(kgrid.z_vec(:))) * 1e3;
This is how I combine both sensors into one. But I don’t then know how to separate the data to the original planes/points, or to be able to plot them directly from the combined dataset.
$ Combined sensors
sensor12.mask = zeros(Nx, Ny, Nz);
sensor12.mask = sensor1.mask + sensor2.mask;
sensor12.mask(sensor12.mask > 1) = 1;
$ Record the simulation
sensor12.record = {'p', 'p_max'};
sensor12_data = kspaceFirstOrder3D(kgrid, medium, source, sensor12, input_args{:});
I would appreciate some help.
Here is the full code (sorry if the comments are in spanish)
clearvars
% Tamaño (en grid points) de la perfectly matching layer (PML)
% Se recomienda de 10 puntos para la simulación en 3D
pml_size = 10;
% Medio (agua)
medium.sound_speed = 1480; % [m/s]
medium.density = 1000; % [kg/m^3]
% Parámetros computacionales
source_freq = 0.5e6; % Frecuencia de la fuente
lambda = medium.sound_speed/source_freq; % Longitud de onda
ppw = 2; % Puntos por longitud de onda
t_end = 80e-6; % Tiempo total computación
cfl = 0.5; % Número CFL
dx = medium.sound_speed/(ppw*source_freq); % Tamaño de muestreo espacial
ppp = round(ppw / cfl); % Puntos por periodo
% Malla computacional 3D
% Para una longitud de eje de 200 mm y dx = c/ppw/f0, se obtendría Nx=270,
% muy cercano a 256, que es un número potencia de 2, recomendable para este
% tipo de simulaciones basados en FFT
Nx = 128; % number of grid points in the x direction
Ny = Nx; % number of grid points in the y direction
Nz = Nx; % number of grid points in the z direction
dy = dx; % grid point spacing in the y direction [m]
dz = dx; % grid point spacing in the z direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
% Array temporal
dt = 1 / (ppp * source_freq); % Tamaño de muestreo temporal
Nt = round(t_end / dt); % Número de puntos temporales
kgrid.setTime(Nt, dt);
% Coordenadas de la malla computacional 3D
% Tenemos que tener en cuenta que el kgrid.x_vec va de -100 a 100
% originalmente, así que lo corregimos sumando el mínimo (-(-100) = +100)
eje_x = (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3;
eje_y = (kgrid.y_vec - min(kgrid.y_vec(:))) * 1e3;
eje_z = (kgrid.z_vec - min(kgrid.z_vec(:))) * 1e3;
%%%%%%%%%%%%%%%%%%%
%%% TRANSDUCTOR %%%
%%%%%%%%%%%%%%%%%%%
% Parámetros físicos
radius_mm = 100e-3; % radio del transductor [mm]
aperture_mm = 110e-3; % apertura del transductor [mm]
radius_points = ceil(radius_mm/dx); % radio del transductor [grid points]
aperture_points = ceil(aperture_mm/dx); % apertura del transductor [grid points]
% Parámetros de posición. Se encuentra centrado en el eje YZ, y sus
% coordenadas son (10, Ny/2, Nz/2)
% El transductor propaga la onda a lo largo del eje x. La posición en los
% ejes z e y se le suma 1 porque MATLAB empieza a contar en 1, no en 0
sphere_offset = 10; % [grid points]
bowl_pos = [sphere_offset + 1, Ny/2+1, Nz/2+1]; % [grid points]
focus_pos = [sphere_offset + radius_points + 1, Ny/2+1, Nz/2+1]; % [grid points]
% El diámetro ha de ser un número impar
if mod(aperture_points, 2) == 0
aperture_points = aperture_points + 1;
end
% Máscara de la fuente
source.p_mask = makeBowl([Nx, Ny, Nz], bowl_pos, radius_points, aperture_points, focus_pos);
% Parámetros de la señal fuente
amp = 1; % [Pa]
source.p = amp * sin(2 * pi * source_freq * kgrid.t_array);
% Se filtra la fuente para eliminar altas frecuencias no soportadas por
% la malla computacional
%source.p = filterTimeSeries(kgrid, medium, source.p);
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% SENSOR Y SIMULACIÓN %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creamos un sensor 1, longitudinal, que cubra el plano XY, con z = Nz/2
sensor1.mask = zeros(Nx, Ny, Nz);
sensor1.mask(:, :, Nz/2 + 1) = 1;
% Creamos otro sensor 2, transversal, que cubra el plano YZ, con x = 0
sensor2.mask = zeros(Nx, Ny, Nz);
% Posición: en el foco geométrico (a 100 mm del extremo del transductor)
sensor2_x = 1 + sphere_offset + radius_points;
sensor2.mask(sensor2_x, :, :) = 1;
% Combinamos ambos sensores, 1 y 2, para ahorrar tiempo de computación
sensor12.mask = zeros(Nx, Ny, Nz);
sensor12.mask = sensor1.mask + sensor2.mask;
sensor12.mask(sensor12.mask > 1) = 1;
%%% Plot 3D del transductor y sensores %%%
% Coordenadas de los puntos del transductor
[xt, yt, zt] = ind2sub(size(source.p_mask), find(source.p_mask));
xt = xt*dx*1e3;
yt = yt*dx*1e3;
zt = zt*dx*1e3;
% Coordenadas de los puntos del sensor 1 (longitudinal)
[x1, y1, z1] = ind2sub(size(sensor1.mask), find(sensor1.mask));
x1 = x1*dx*1e3;
y1 = y1*dx*1e3;
z1 = z1*dx*1e3;
% Coordenadas de los puntos del sensor 2 (transversal)
[x2, y2, z2] = ind2sub(size(sensor2.mask), find(sensor2.mask));
x2 = x2*dx*1e3;
y2 = y2*dx*1e3;
z2 = z2*dx*1e3;
% Coordenadas de los puntos del sensor combinado
[x_12, y_12, z_12] = ind2sub(size(sensor12.mask), find(sensor12.mask));
x_12 = x_12*dx*1e3;
y_12 = y_12*dx*1e3;
z_12 = z_12*dx*1e3;
% Se crea una nueva figura y utiliza plot3 para graficar las coordenadas
figure('Name', 'Configuración 3D');
plot3(x_12, y_12, z_12, '.', 'Color',"#D95319");
hold on
plot3(xt, yt, zt, '.', 'Color', "#0072BD");
xlabel('{it x} / mm');
ylabel('{it y} / mm');
zlabel('{it z} / mm');
title('Sensores y transductor')
legend('Sensor', 'Transductor');
grid on
hold off
%%
%%%%%%%%%%%%%%%%%%%%%
%%% VISUALIZACIÓN %%%
%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Sensor 1 (plano longitudinal) %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
input_args = {'DisplayMask', source.p_mask, 'PMLSize', pml_size, 'DataCast', 'single', 'PlotSim', true};
sensor12.record = {'p', 'p_max'};
sensor12_data = kspaceFirstOrder3D(kgrid, medium, source, sensor12, input_args{:});
% Extraer resultados de sensor1 usando su máscara
sensor1_data_p_max = sensor12_data.p_max .* sensor1.mask;
sensor1_data_p = sensor12_data.p .* sensor1.mask;
% Extraer resultados de sensor2 usando su máscara
sensor2_data_p_max = sensor12_data.p_max .* sensor2.mask;
sensor2_data_p = sensor12_data.p .* sensor2.mask;
%sensor12_data.p_max = reshape(sensor12_data.p_max, [Ny, Nx]);
%%
% Extraemos un corte en el tiempo de presión
t_corte = Nt-2;
p1_t = sensor1_data.p(:,t_corte);
p1_t = reshape(p1_t, [Ny, Nx]);
%p1_media = mean(sensor1_data.p, 3);
% Mapa de presión máxima del plano longitudinal
figure('Name', 'Plano YX');
imagesc(eje_x, eje_y, sensor1_data.p_max);
xlabel('y / mm');
ylabel('x / mm');
title('Presión máxima detectada en el plano z = 100');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Presión / Pa');
axis image;
% Mapa de presión del plano longitudinal
figure('Name', 'Plano XY');
imagesc(eje_x, eje_y, p1_t);
xlabel('y / mm');
ylabel('x / mm');
title('Presión detectada en el plano z = 100');
colormap(jet(256));
c = colorbar;
ylabel(c, 'Presión / Pa');
axis image;
%%% Plot sensor 1 (plano longitudinal) %%%
% Extraemos la amplitud del sensor a lo largo de la línea central en y = Ny/2
amp_max_sensor1 = sensor1_data.p_max(:, Ny/2 + 1);
amp_p1_t = p1_t(:, Ny/2 + 1);
figure('Name', 'Plot de presión máxima longitudinal')
plot(eje_x, amp_max_sensor1);
hold on
% Añadimos una línea vertical de referencia de donde está el transductor.
% Para ello convertimos la posición de referencia a mm, multiplicando *dx
ref_transductor = sphere_offset*dx*1e3;
xline(ref_transductor, '--m');
% Añadimos una línea vertical que indique la posición del foco geométrico
% (radio). Lo hacemos con las variables que sabemos, ref_transductor y radio
xline(ref_transductor + radius_mm*1e3, '--b');
% Añadimos una línea vertical que indique la posición del máximo de
% amplitud de focalización en el eje longitudinal. Para ello buscamos el
% índice que corresponde al máximo del plot, max_index. Tras ello buscamos
% el valor de kgrid.x_vec que corresponde a tal índice. Luego corregimos de
% nuevo ese valor para plotearlo, teniendo en cuenta que va de -100 a 100 mm.
[~, max_index] = max(amp_max_sensor1);
max_amplitud_long = (kgrid.x_vec(max_index) - min(kgrid.x_vec(:)))* 1e3;
xline(max_amplitud_long, '--r');
title('Presión máxima eje x')
xlabel('x / mm')
ylabel('Presión máxima / Pa')
legend('Simulación', 'Transductor', 'Foco geométrico (radio)', 'Máximo amplitud')
figure('Name', 'Plot de presión longitudinal')
plot(eje_x, amp_p1_t);
hold on
xline(ref_transductor, '--m');
xline(ref_transductor + radius_mm*1e3, '--b');
xline(max_amplitud_long, '--r');
title('Presión eje x en estado estacionario')
xlabel('x / mm')
ylabel('Presión máxima / Pa')
legend('Simulación', 'Transductor', 'Foco geométrico (radio)', 'Máximo amplitud')