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')
0 个评论
回答(1 个)
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.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Create Block Masks 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!