clear all

close all

% simulation settings

data_cast = 'single'; % set to 'single' or 'gsingle' to speed up computations

run_simulation = true;

 

% =========================================================================

% DEFINE THE K-WAVE GRID

% =========================================================================

 

% set the size of the perfectly matched layer (PML)

PML_X_SIZE = 20; % [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 = 128 - 2*PML_X_SIZE;   % [grid points]

Ny = 64 - 2*PML_Y_SIZE;   % [grid points]

Nz = 64 - 2*PML_Z_SIZE;   % [grid points]

 

% set desired grid size in the x-direction not including the PML

x = 40e-3 

% calculate the spacing between the grid points

dx = x/(4*Nx); % [m]

dy = dx; % [m]

dz = dx; % [m]

 

% create the k-space grid

kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);

 

sensor.mask = zeros(Nx, Ny, Nz);    

x = [max(kgrid.x(:,1,1)) 0 -max(kgrid.x(:,1,1))/2 -max(kgrid.x(:,1,1))/4*3 -max(kgrid.x(:,1,1))+3*dx];

y = [-max(kgrid.y(1,:,1))+7*dy -max(kgrid.y(1,:,1))+7*dy -max(kgrid.y(1,:,1))+7*dy -max(kgrid.y(1,:,1))+7*dy -max(kgrid.y(1,:,1))+7*dy];

z = [0 0 0 0 0];

 

sensor.mask = [x; y; z];

voxelPlot(double(cart2grid(kgrid, sensor.mask)));

view([50, 20]);

 

% =========================================================================

% DEFINE THE MEDIUM PARAMETERS

% =========================================================================

 

% define the properties of the propagation medium

medium.sound_speed = 1529; % [m/s]

medium.density = 1056; % [kg/m^3]

medium.alpha_coeff = 0.1256;   %a [dB/(MHz^y cm)]

medium.alpha_power = 1.001; %b

medium.BonA = 8.1;

 c0=medium.sound_speed;

 rho0=medium.density;

 a0 = medium.alpha_coeff;

 b0 = medium.alpha_power;

t_end = (Nx*dx)*2.2/medium.sound_speed;

kgrid.t_array = makeTime(kgrid, medium.sound_speed, 0.1, t_end);

 

 

% =========================================================================

% DEFINE THE INPUT SIGNAL

% =========================================================================

% define a single source element

source.p_mask = zeros(Nx, Ny, Nz);

source.p_mask(1, Ny/4 , Nz/2) = 1;

 

source.p = zeros(size(kgrid.t_array));

source.p(100) =1;

figure

plot(kgrid.t_array,source.p)

figure

source_func_filtered = filterTimeSeries(kgrid, medium, source.p,'ZeroPhase', true);

plot(kgrid.t_array,source_func_filtered)

source.p =source_func_filtered;

 

% =========================================================================

% DEFINE THE MEDIUM PROPERTIES

% =========================================================================

 

% define a large image size to move across

Nx_tot = Nx;

Ny_tot = Ny; % + number_scan_lines*transducer.element_width;

Nz_tot = Nz;

 

%define a random distribution of scatterers for the medium

background_map_mean = 1;

background_map_std = 0.008;

background_map = background_map_mean + background_map_std*randn([Nx_tot, Ny_tot, Nz_tot]);

 

c1 = 1688;

rho1 = 1920;

sound_speed_map = c0*ones(Nx_tot, Ny_tot, Nz_tot).*background_map;

density_map = rho0*ones(Nx_tot, Ny_tot, Nz_tot).*background_map;

bone_region = zeros(Nx_tot, Ny_tot, Nz_tot);

scattering_rho1 = rho1*ones(Nx_tot, Ny_tot, Nz_tot);

scattering_c1 = c1*ones(Nx_tot, Ny_tot, Nz_tot);

bone_region(:, Ny/2 - 7:Ny/2 +7, Nz/2 - 7:Nz/2 + 7)=1;

density_map(bone_region == 1) = scattering_rho1(bone_region == 1); 

sound_speed_map(bone_region == 1) = scattering_c1(bone_region == 1); 

 

 

% =========================================================================

% RUN THE SIMULATION

% =========================================================================

 

 

% set the input settings

input_args = {'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

 

 

% load the current section of the medium

medium.sound_speed = sound_speed_map(:, 1:Ny, :);

medium.density = density_map(:, 1:Ny, :);

   

  % run the simulation  

sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});

 

 

 

end 

figure;

imagesc((0:Ny_tot-1)*dy*1e3, (0:Nx_tot-1)*dx*1e3, sound_speed_map(:, :, Nz/2));

axis image;

colormap(gray);

title('Scattering Phantom');

xlabel('Horizontal Position [mm]');

ylabel('Depth [mm]');

 

aa=zeros(44,44);

for i=1:44

for j=1:44

aa(i,j) = sound_speed_map(1, i, j);

end

end

figure;

imagesc((0:Ny_tot-1)*dy*1e3, (0:Nz_tot-1)*dz*1e3, aa);

axis image;

colormap(gray);

title('Scattering Phantom');

xlabel('Horizontal Position [mm]');

ylabel('Depth [mm]');