Separating and indexing (k-wave) datasets to different planes/coordinates

I have been using k wave in MATLAB to record intensity patterns emited by a source, at different planes of a 3D simulation grid. I have defined a sensor for each plane, and then ran the simulation with the kspaceFirstOrder3D function. But this requires more time as I define more sensors, since I have to run a simulation for each sensor.

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:

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;

$ 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);

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, if it helps (sorry 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};

% Grabamos los últimos periodos 
%num_periods = 4;
%T_points = round(num_periods * (1/source_freq) / kgrid.dt);
%sensor1.record_start_index = Nt-100;


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')

Answer 1

[Voss]

One way to separate the the data for the two sensors is using the following approach

% Extraer resultados de sensor1 y sensor2 usando sus máscaras

% x,y,z grid points over the entire grid

[x_idx,y_idx,z_idx] = ndgrid(1:Nx,1:Ny,1:Nz);

% x,y,z where sensor12.mask is 1

tmp = logical(sensor12.mask);

sensor12_xyz = [x_idx(tmp) y_idx(tmp) z_idx(tmp)];

% x,y,z where sensor1.mask is 1

tmp = logical(sensor1.mask);

sensor1_xyz = [x_idx(tmp) y_idx(tmp) z_idx(tmp)];

% x,y,z where sensor2.mask is 1

tmp = logical(sensor2.mask);

sensor2_xyz = [x_idx(tmp) y_idx(tmp) z_idx(tmp)];

% indices where sensor12_xyz is in sensor1_xyz

[~,idx] = ismember(sensor1_xyz,sensor12_xyz,'rows');

% store the results at those indices as sensor1_data

sensor1_data.p = sensor12_data.p(idx,:);

sensor1_data.p_max = reshape(sensor12_data.p_max(idx),Nx,Ny); % reshape p_max for use in imagesc

% indices where sensor12_xyz is in sensor2_xyz

[~,idx] = ismember(sensor2_xyz,sensor12_xyz,'rows');

% store the results at those indices as sensor2_data

sensor2_data.p = sensor12_data.p(idx,:);

sensor2_data.p_max = reshape(sensor12_data.p_max(idx),Ny,Nz); % reshape p_max for use in imagesc

which should be done immediately after getting the results from kspaceFirstOrder3D.

And a little further down, this

p1_t = reshape(p1_t, [Ny, Nx]);

should be this

p1_t = reshape(p1_t, [Nx, Ny]);

This is based on the fact that p1_t's second dimension is expected to correspond to the y spatial dimension, which is clear from this subsequent operation:

amp_p1_t = p1_t(:, Ny/2 + 1);

That change makes no difference when Nx and Ny equal.

The complete modified script is attached.