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

12 次查看(过去 30 天)
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')

回答(1 个)

Voss
Voss 2024-9-8
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.

类别

Help CenterFile Exchange 中查找有关 Create Block Masks 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by