Hello to everyone,
I am a new user of k-Wave.
I am trying to simulate a image acquisition of a muscle + pleura + lung phantom in a water medium. Lung ultrasonography is based in reverberations artifacts. Multiple horizontal lines can be seen after the pleural line as result of the big impedance difference between tissue and air inside the lung. Ultrasound signal reflects multiple times between pleura and transducer, and this artifact helps to do lung diagnosis.
I tried to simulate that effect in a 2-D model. As lung air is much more abundant than tissue, I used air velocity and density as parameters of the lung. Although, result image is all black with only the pleura line visible.
Is it possible to reproduce this reverberation with K-Wave?
My code was that:
% % =========================================================================
% % 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]
% set total number of grid points not including the PML
Nx = 500+40-2*pml_x_size; % [grid points]
Ny = 240+40-2*pml_y_size; % [grid points]
% set desired grid size in the x-direction not including the PML
x = 1e-2; % [m]
y = 48e-4;
% calculate the spacing between the grid points
dx = x/Nx; % [m]
dy = dx; % [m]
% create the k-space grid
kgrid = kWaveGrid(Nx, dx, Ny, dy);
% define a large image size to move across
number_scan_lines = 200;
Nx_tot = Nx;
Ny_tot = Ny + number_scan_lines;
% scatter size
number_scatters1 = 32;
number_scatters5 = 512;
area_total=x*y*1000000;
n_total1=area_total*number_scatters1;
n_total5=area_total*number_scatters5;
matrix1 = gerar_matrix(Nx_tot,Ny_tot,n_total1);
matrix5 = gerar_matrix(Nx_tot,Ny_tot,n_total5);
% =========================================================================
% DEFINE THE MEDIUM PARAMETERS
% =========================================================================
clearvars -except T pulse_freq pml_x_size pml_y_size Nx Ny ...
x y dx dy kgrid resolution_x resolution_y area_total n_total ...
matrix1 matrix2 matrix3 matrix4 matrix5 matrix6 number_scan_lines Nx_tot Ny_tot number_scatters...
x_pos mapUp mapDown c0 rho0 medium.alpha_coeff medium.alpha_power...
medium.BonA t_end kgrid.t_array
T=37;
pulse_freq = 5e6;
vMuscle = 1471.4+3.9979*T-0.03513*T.^2
rhoMuscle = 1065
vPleura = 1613;
rhoPleura = 1120;
vAir = 331.45*((T+273.15)/273.15)^(1/2);
rhoAir = 101325/(287.058*(T+273.15));
vWater = speedSoundWater(T);
rhoWater = 1000;
% define the properties of the propagation medium1
c0 = vWater; % [m/s]
rho0 = rhoWater; % [kg/m^3]
medium.alpha_coeff = 0.002; % [dB/(MHz^y cm)]
medium.alpha_power = 2;
medium.BonA = 5.2;
% create the time array
t_end = (Nx*dx)*2.2/c0; % [s]
CFL = 0.05;
kgrid.t_array = makeTime(kgrid, c0, CFL, t_end);
% =========================================================================
% DEFINE THE INPUT SIGNAL
% =========================================================================
% define properties of the input signal
source_strength = 5e6;%*blackman(128); % [Pa]
tone_burst_freq = pulse_freq; % [Hz]
tone_burst_cycles = 4;
% lambda
lambda=c0/tone_burst_freq;
% create the input signal using toneBurst
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);
% scale the source magnitude by the source_strength divided by the
% impedance (the source is assigned to the particle velocity)
input_signal = source_strength/(c0*rho0)*input_signal;
% =========================================================================
% DEFINE THE ULTRASOUND SOURCE
% =========================================================================
source_width = 128; % [grid points]
source.p_mask = zeros(Nx, Ny);
source.p_mask(1, Ny/2 - source_width/2+1:Ny/2 + source_width/2) = 1;
sensor_focus = [45e-4,0]; %[m]
input_signal = focus(kgrid,input_signal,source.p_mask,sensor_focus,vWater);
source.p=input_signal;
% =========================================================================
% DEFINE THE ULTRASOUND SENSOR
% =========================================================================
sensor_offset=1;
sensor.mask = zeros(Nx, Ny);
sensor.mask(sensor_offset, Ny/2 - source_width/2+1: Ny/2 + source_width/2) = 1;
sensor.record={'p'};
% =========================================================================
% DEFINE THE MEDIUM PROPERTIES
% =========================================================================
density_map=ones(size(matrix1)).*rhoWater;
sound_speed_map=ones(size(matrix1)).*vWater;
sound_speed_map(matrix1 == 0) = 0;
density_map(matrix1 == 0) = 0;
% define a sphere for a highly scattering region
x_pos = 25e-5; % [m]
y_pos = 44e-4; % [m]
%%%%%%%%%%%%%%%
Linha1=makeLine(Nx_tot, Ny_tot,[(round(x_pos/dx)),round((y_pos-18e-4)/dy)],[(round(x_pos/dx)),round((y_pos+18e-4)/dy)]);
for i=1:125
Linha1=Linha1+makeLine(Nx_tot, Ny_tot,[(round(x_pos/dx)+i),round((y_pos-18e-4)/dy)],[(round(x_pos/dx)+i),round((y_pos+18e-4)/dy)]);
end
%%%%%%%%%%%%%%%%%%%%
Linha2=makeLine(Nx_tot, Ny_tot,[(round(x_pos/dx)+126),round((y_pos-18e-4)/dy)],[(round(x_pos/dx)+126),round((y_pos+18e-4)/dy)]);
for i=127:150
Linha2=Linha2+makeLine(Nx_tot, Ny_tot,[(round(x_pos/dx)+i),round((y_pos-18e-4)/dy)],[(round(x_pos/dx)+i),round((y_pos+18e-4)/dy)]);
end
%%%%%%%%%%%%%%%%%
Linha3=makeLine(Nx_tot, Ny_tot,[(round(x_pos/dx)+151),round((y_pos-18e-4)/dy)],[(round(x_pos/dx)+151),round((y_pos+18e-4)/dy)]);
for i=152:450
Linha3=Linha3+makeLine(Nx_tot, Ny_tot,[(round(x_pos/dx)+i),round((y_pos-18e-4)/dy)],[(round(x_pos/dx)+i),round((y_pos+18e-4)/dy)]);
end
%%%%%%%%%%%%%%%%%
% assign region
sound_speed_map(Linha1 == 1 ) = matrix5(Linha1 == 1)*vMuscle;
density_map(Linha1 == 1) = matrix5(Linha1 == 1)*rhoMuscle;
sound_speed_map(Linha2 == 1 ) = matrix5(Linha2 == 1)*vPleura;
density_map(Linha2 == 1) = matrix5(Linha2 == 1)*rhoPleura;
sound_speed_map(Linha3 == 1 ) = matrix5(Linha3 == 1)*vAir;
density_map(Linha3 == 1) = matrix5(Linha3 == 1)*rhoAir;
sound_speed_map(sound_speed_map==0)=vWater;
density_map(density_map==0)=rhoWater;
% plot the medium, truncated to the field of view
figure;
imagesc((0:number_scan_lines-1)*dy*1e3, (0:Nx_tot-1)*dx*1e3, density_map(:, 1 + Ny/2:end - Ny/2));
axis image;
colormap(gray);
set(gca, 'YLim', [1, 10]);
title('Scattering Phantom');
xlabel('Horizontal Position [mm]');
ylabel('Depth [mm]');
% ========================================================================
% RUN THE SIMULATION
% ========================================================================
% preallocate the storage
scan_lines = zeros(number_scan_lines, kgrid.Nt);
% set the input settings
DATA_CAST = 'gpuArray-single'; % set to 'single' or 'gpuArray-single' to speed up computations
input_args = {...
'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size], ...
'DataCast', DATA_CAST, 'DataRecast', true, 'PlotSim', false};
% set medium position
medium_position = 1;
for scan_line_index = 1:number_scan_lines
% update the command line status
disp(['|| =============== T = ' num2str(T) ' =============== ||']);
disp('');
disp(['|| =============== Freq = ' num2str(pulse_freq) ' =============== ||']);
disp('');
disp(['Computing scan line ' num2str(scan_line_index) ' of ' num2str(number_scan_lines)]);
% load the current section of the medium
medium.sound_speed = sound_speed_map(:, medium_position:medium_position + Ny - 1);
medium.density = density_map(:, medium_position:medium_position + Ny - 1);
% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
% extract the scan line from the sensor data
scan_lines(scan_line_index, :) = sum(bsxfun(@times, hanning(128),sensor_data.p)); %bsxfun(@times, tgc, scan_lines);
% update medium position
medium_position = medium_position + 1;
% save(['Sim_2_resultsWindow',num2str(number_scatters)])
end
save(['sim_Rectangle_lus_CFL005'])
% end
% end