Dear Dr. Treeby and Dr. Cox,
Like, my previous question in:
http://www.k-wave.org/forum/topic/understanding-speed-of-sound-of-the-medium#post-6190
For the 2nd case, I want to simulate plane wave ultrasound (128 elements array) propagation through homogeneous medium (water)that I generate by integrating 2D images (see the 2nd section of my code and the link for the image source).
For this simulation, I got an error that said:
Error using matlab.graphics.axis.Axes/set
and a warning: MATLAB has disabled some advanced graphics rendering featuresby switching to software openGL.
By the way, I use MATLAB 2016a version.
Thank you.
Here is my code:
% =========================================================================
% ULTRASOUND SIMULATION
% Linear array transducer (128 element)
% Transmit: Planewave, 128 active elements
% Receive : 128 active elements
% Medium : water (256 x 512 x 256)
% =========================================================================
clear all;
%simulation setting
data_cast = 'single';
% run the simulation
% true : run the simulation
% false : load previous results from disk
run_simulation = true;
% =========================================================================
% DEFINE THE K-WAVE GRID
% =========================================================================
% set the size of the perfectly matched layer (PML)
PML_X_SIZE = 10; % [grid points]
PML_Y_SIZE = 10; % [grid points]
PML_Z_SIZE = 10; % [grid points]
% set total number of grid points not including the PML
Nx = 256 - 2*PML_X_SIZE; % [grid points]
Ny = 512 - 2*PML_Y_SIZE; % [grid points]
Nz = 256 - 2*PML_Z_SIZE; % [grid points]
% set desired grid size in the x-direction not including the PML
x = 38e-3; % [m]
y = 50e-3; % [m]
z = 38e-3; % [m]
% calculate the spacing between the grid points
dx = z/Nz; % [m]
dy = dx; % [m]
dz = dx; % [m]
% create the k-space grid
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
% =========================================================================
% DEFINE THE MEDIUM PARAMETERS
% =========================================================================
c0 = 1498; % c in water [m/s]
rho0 = 1000; % density of water [kg/m^3]
% create the time array
t_end = (Nx*dx)*2.2/c0; % [s]
kgrid.t_array = makeTime(kgrid, c0, [], t_end);
% =========================================================================
% DEFINE THE INPUT SIGNAL
% =========================================================================
% define properties of the input signal
source_strength = 1e6; % [Pa]
tone_burst_freq = 2e6; % [Hz]
tone_burst_cycles = 10;
% 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 TRANSDUCER [SOURCE AND SENSOR]
% =========================================================================
% physical properties of the transducer
transducer.number_elements = 128; % total number of transducer elements
transducer.element_width = 2; % width of each element [grid points]
transducer.element_length = 25; % 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;
% 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]);
transducer.sound_speed = c0; % sound speed [m/s]
transducer.elevation_focus_distance = 16e-3; % focus distance in the elevation plane [m]
% apodization
transducer.transmit_apodization = 'Hanning';
transducer.receive_apodization = 'Rectangular';
% define the transducer elements that are currently active
number_active_elements = 128;
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 = makeTransducer(kgrid, transducer);
% print out transducer properties
transducer.properties;
% =========================================================================
% DEFINE THE MEDIUM PROPERTIES
% =========================================================================
% define a random distribution of scatterers for the medium
background_map_mean = 1;
background_map_std = 0.008;
% % define the properties of the propagation medium
load c_water_171106;
medium.sound_speed = c_water_171106; % [m/s]
load rho_water_171106;
medium.density = rho_water_171106; % [kg/m^3]
% =========================================================================
% RUN THE SIMULATION
% =========================================================================
% set the input settings
% plot medium map
input_args = {...
'PlotLayout', true, 'PMLInside', false, 'PlotSim', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
'DataCast', data_cast};
% run the simulation if set to true, otherwise, load previous results from disk
if run_simulation
RF_data_PW_water_171107 = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:});
% save the scan lines to disk
save RF_data_PW_water_171107 RF_data_PW_water_171107;
else
% % load the scan lines from disk
load RF_data_PW_water_171107
end=============================================================
and this is the code for creating the medium:
clear all
close all
[data]=imread('water_2D.bmp');
data_bw=double(data(:,:,1));
[xx,yy]=size(data_bw);
clear data
for i=1:xx
for j=1:yy
if data_bw(i,j)==1
c(i,j)=1498; % average c of the medium (phantom and water) [m/s]
rho(i,j)=1000; % average density of the medium (phantom and water) [kg/m^3]
elseif data_bw(i,j)==0
c(i,j)=1498; % c in water [m/s]
rho(i,j)=1000; % density of water [kg/m^3]
end;
end;
end;
cit_3d=zeros(xx,yy,492);
for i=1:492
c_3d(:,:,i)=c;
rho_3d(:,:,i)=rho;
end;
n=zeros(236,492,236);
m=zeros(236,492,236);
[x,y,z]=size(n);
for i=1:x
for j=1:y
for k=1:z
c_water_171106(i,j,k)=c_3d(i,k,j);
rho_water_171106(i,j,k)=rho_3d(i,k,j);
end
end
end
save c_water_171106 c_water_171106
save rho_water_171106 rho_water_171106==========================================================
and here is the image:
https://drive.google.com/file/d/1may-MoszFDxJgkQpyW6mjoD57Rs9LJJ0/view?usp=sharing