Hi,
i met a problem while simulating with k-wave.
The case i test is keeping the 3-D medium's properties all the same but alpha coef which is the attenuation coeff at the middle of the area with a cylinder shape. I set c0 = 1540m/s rho0 = 1000, alpha_power = 1.5 or 1.01, and a map of alpha coef, within the cylinder, it's 1 and the other area is 0.5. i'm imaging the area with single angle planewave.
Theoretically, i should get an only black backscattered bmode image becausethere's no impedence mismatch which related to rho and c,so there should not be a reflector in the backscattered image.But what i actually got is it seems like there are two reflectors at the top and bottom edge of the circle area. Can anyone help explain this?
And during the testing, i met another problem about steering angle of the transducer,i wanted to tranmit steered planewaves and compound them into one image, but it seems like the steer is not electronic steering, it's mechanically steer the transducer, so when i do delay and sum, i should not add the steered angle there. Am i understanding it correct?
Thank you! and the following is the code.
clearvars;close all;
addpath(genpath([pwd,'/k-wave-toolbox-version-1.3']))
% simulation settings
DATA_CAST = 'gpuArray-single'; % set to 'single' or 'gpuArray-single' to speed up computations
RUN_SIMULATION = true; % set to false to reload previous results instead of running simulation
% =========================================================================
% DEFINE THE K-WAVE GRID
% =========================================================================
% set the size of the perfectly matched layer (PML)
pml_x_size = 20; % [grid points]
pml_y_size = 20; % [grid points]
pml_z_size = 10; % [grid points]
% set total number of grid points not including the PML
%Nx = 1024+44 - 2 * pml_x_size; % [grid points]
Nx = 512+44 - 2 * pml_x_size; % [grid points]
%Ny = 1024+44 - 2 * pml_y_size; % [grid points]
Ny = 512+44 - 2 * pml_y_size; % [grid points]
Nz = 128 - 2 * pml_z_size; % [grid points]
% set desired grid size in the x-direction not including the PML
x = 50e-3; % [m]
% calculate the spacing between the grid points
dx = x / Nx; % [m]
dy = dx;%1.923076923076923e-04; % [m]
%dz = 2*dx; % [m]
dz = 2e-4;
% create the k-space grid
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
% =========================================================================
% DEFINE THE MEDIUM PARAMETERS
% =========================================================================
% define the properties of the propagation medium
c0 = 1500; %for water % [m/s]
rho0 = 1000; % [kg/m^3]
% medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)] % need to change this paramenters to map also
medium.alpha_power = 1.01;
%medium.BonA = 20;
%BonA non linear
% create the time array
t_end = (Nx * dx) * 2.2 / c0; % [s]
%t_end = (Nx * dx) * 2 / c0; % [s]
kgrid.makeTime(c0, [], t_end);
%%
% =========================================================================
% DEFINE THE INPUT SIGNAL
% =========================================================================
% define properties of the input signal
source_strength = 10e5; % [Pa]
tone_burst_freq = 2e6; % [Hz]
tone_burst_cycles = 4;
% create the input signal using toneBurst
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);
input_signal = (source_strength ./ (c0 * rho0)) .* input_signal;
% =========================================================================
% DEFINE THE ULTRASOUND TRANSDUCER
% =========================================================================
% physical properties of the transducer
transducer.number_elements = 128; % total number of transducer elements
transducer.element_width = 3; % width of each element [grid points]
transducer.element_length = 24; % length of each element [grid points]
transducer.element_spacing = 0; % spacing (kerf width) between the elements [grid points]
transducer.radius = inf; % radius of curvature of the transducer [m]
% calculate the width of the transducer in grid points
transducer_width = transducer.number_elements * transducer.element_width ...
+ (transducer.number_elements - 1) * transducer.element_spacing;
pitch = transducer.element_width*dx;
%transducer.transmit_apodization = 'Hanning';
%%
% use this to position the transducer in the middle of the computational grid
transducer.position = round([1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]);
% properties used to derive the beamforming delays
transducer.sound_speed = c0; % sound speed [m/s]
transducer.focus_distance = Inf; % focus distance [m]
transducer.elevation_focus_distance = 25e-3; % focus distance in the elevation plane [m]
%transducer.steering_angle = 0; % steering angle [degrees]
% apodization
transducer.transmit_apodization = 'Hanning';
transducer.receive_apodization = 'Rectangular';
% define the transducer elements that are currently active
transducer.active_elements = ones(transducer.number_elements, 1);
% append input signal used to drive the transducer
transducer.input_signal = input_signal;
%
% % create the transducer using the defined settings
% transducer = kWaveTransducer(kgrid, transducer);
% print out transducer properties
transducer = kWaveTransducer(kgrid, transducer);
transducer.properties;
% =========================================================================
% DEFINE THE MEDIUM PROPERTIES
% =========================================================================
% define a large image size to move across
Nx_tot = Nx;
Ny_tot = Ny;
Nz_tot = Nz;
alpha_coeff = 0.5;
alpha_coeff2 = 1;
scattering_map = randn([Nx_tot, Ny_tot, Nz_tot]);
% define a random distribution of scatterers for the medium
background_map_mean = 1;
background_map_std = 0;
background_map = background_map_mean + background_map_std * scattering_map;
phanton_map_std = 0;
% define properties
%background_speed = c0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map;
sound_speed_map = c0 .* background_map;
density_map = rho0 .* background_map;
alpha_coeff_map = alpha_coeff.* background_map;
radius = 4e-3; % [m]
x_pos = 25e-3; % [m]
y_pos = 25e-3; % [m]
for zi = 1:Nz_tot
alpha_coeff_slice = alpha_coeff_map(:,:,zi);
scattering_regionb = zeros(Nx_tot, Ny_tot);
scattering_regionb = makeDisc(Nx_tot, Ny_tot, round(x_pos/dx), round(y_pos/dy), round(radius/dx));
scattering_slice = scattering_map(:,:,zi);
scattering_c02 = alpha_coeff2.*(1 + phanton_map_std * scattering_slice);
alpha_coeff_slice(scattering_regionb == 1) = scattering_c02(scattering_regionb == 1);
alpha_coeff_map(:,:,zi) = alpha_coeff_slice;
end
alpha_coeff_slice_ori = alpha_coeff_slice;
% =========================================================================
% RUN THE SIMULATION
% =========================================================================
% set the input settings
input_args = {...
'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size, pml_z_size], ...
'DataCast', DATA_CAST, 'DataRecast', true, 'PlotSim', false};
i = 0;
medium.sound_speed = c0;%sound_speed_map;%(:, medium_position:medium_position + Ny - 1, :);
medium.density = rho0;%density_map;%(:, medium_position:medium_position + Ny - 1, :);
for degree= 0:0 % rotate
disp(['Computing scan angle ' num2str(degree)]);
rot_alpha_coeff_slice = rotate_kgrid_dbz(alpha_coeff_slice_ori, degree, alpha_coeff);
for steer_ang = -1:1:1
i = i+1;
transducer.steering_angle = steer_ang;
%transducer.input_signal = input_signal;
% create the transducer using the defined settings
disp(['Computing steer angle ' num2str(steer_ang)]);
for zi = 1:Nz_tot
alpha_coeff_map(:,:,zi) = rot_alpha_coeff_slice;
end
medium.alpha_coeff = alpha_coeff_map;%(:, medium_position:medium_position + Ny - 1, :);
sensor_data = kspaceFirstOrder3DG(kgrid, medium, transducer, transducer, input_args{:});
back_rf(:,:,i) = sensor_data;
%ps_data(:,:,i) = bf_planewave(sensor_data,1/kgrid.dt);
end
end
figure();imagesc(20*log10(abs(hilbert(back_rf(:,:,1)'))));title('raw rf data')