Hi Treeby and colleagues,
I am trying to simulate ultrasound stimulation in the human brain.
I defined the Kgrid as the size of T1 weighted image, in this case, 192*256*256, with dx=dy=dz=0.8e-3.
I defined the medium using segmentation result which separate skull, soft tissue and scalp, and set each density, sound_speed and alpha_coeff according to literatures.
I set up the source using makeBowl function to mimic a ultrasound transducer.
I set up the sensor so that it covers the whole Kgrid.
I set up the input signal as what 'example_diff_focused_ultrasound_heating' did.
The simulation worked well at the first half of computing, while went very strange during the second half, as shown in the figure below.
The following two figures show how the simulation result changes over time. You can access to the figure by the URL.
https://ibb.co/d60B7ch
https://ibb.co/LRFQDS0
This only happens when I introduced the skull into medium.
Once I removed the skull from the medium, i.e., left the skull part as free water, the simulation became reasonable.
The simulation also works well if I select one X-Y plane and do a 2D simulation, with skull defined in medium.
Can you advise what maight account for the strange simulation result above, and how I can modify my code to fix things? I have attached my code for simulation below. Can you also have a look on it? Thanks in advance!
P.S.:If you need the segmentation file for replication, please tell me and I'll upload it to somewhere like the google drive.
%% Code for simulation
clear; close all; clc;
%% get segmentation
T1_nii = load_nii('final_tissues.nii.gz');
T1_recon = T1_nii.img;
%% specify K-grid
Nx = size(T1_recon,1); Ny = size(T1_recon,2); Nz = size(T1_recon,3); %get T1WI resolution
%x - 1:left to max:right, y - 1:back to max:front, z - 1:bottom to max:top
dx = 0.8e-3; dy = 0.8e-3; dz = 0.8e-3; %get T1WI step
if ~rem(Ny,2)
center_y = Ny/2;
else
center_y = (Ny+1)/2;
end
if ~rem(Nz,2)
center_z = Nz/2;
else
center_z = (Nz+1)/2;
end
% create the computational grid
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
%% specify medium, parameter from Fomenko et al., 2020, Elife
medium.density = 1000.*ones(Nx, Ny, Nz);
medium.sound_speed = 1482.*ones(Nx, Ny, Nz);
medium.alpha_coeff = 0.21.*ones(Nx, Ny, Nz);
medium.alpha_power = 1.1;
medium.alpha_mode = 'no_dispersion';
% brain tissues, gm&wm
medium.density(T1_recon==1|T1_recon==2) = 1040;
medium.sound_speed(T1_recon==1|T1_recon==2) = 1552;
medium.alpha_coeff(T1_recon==1|T1_recon==2) = 0.21;
% scalp
medium.density(T1_recon==5) = 1100;
medium.sound_speed(T1_recon==5) = 1732;
medium.alpha_coeff(T1_recon==5) = 0.21;
% bone
% medium.density(T1_recon==7|T1_recon==8) = 1732;
% medium.sound_speed(T1_recon==7|T1_recon==8) = 2850;
% medium.alpha_coeff(T1_recon==7|T1_recon==8) = 2.7;
%% define the transucer
aperture_grid = round(64e-3/dx)+1;
curvature_grid = round(64e-3/dx);
pos_transducer = [1, center_y ,center_z]; %where to put the transducer
source.p_mask = makeBowl([Nx, Ny, Nz], pos_transducer, ...
curvature_grid, aperture_grid, [curvature_grid, center_y, center_z], 'Plot', true);
%% define input signal
% define the acoustic parameters
freq = 5e5; % [Hz] 500KHz
amp = 0.5e6; % [Pa]
% calculate the time step using an integer number of points per period
sound_speed = 1482; % [m/s]
ppw = sound_speed / (freq * dx); % points per wavelength
cfl = 0.3; % cfl number
ppp = ceil(ppw / cfl); % points per period
T = 1 / freq; % period [s]
dt = T / ppp; % time step [s]
% calculate the number of time steps to reach steady state
t_end = sqrt( kgrid.x_size.^2 + kgrid.y_size.^2 + kgrid.z_size.^2 ) / sound_speed;
Nt = round(t_end / dt);
% create the time array
kgrid.setTime(Nt, dt);
% define the input signal
source.p = createCWSignals(kgrid.t_array, freq, amp, 0, 0); % ramp_length=0
% figure;
% plot(source.p);
%% define the sensor mask
% set the sensor mask to cover the entire grid
sensor.mask = ones(Nx, Ny, Nz);
% sensor.record = {'p', 'p_max_all', 'I', 'I_avg'};
sensor.record = {'p', 'p_max_all'};
%% run simulation
% record the last 3 cycles in steady state
num_periods = 3;
T_points = round(num_periods * T / kgrid.dt);
sensor.record_start_index = Nt - T_points + 1;
% set the input arguements
input_args = {'PMLInside', false, 'PlotPML', false, 'PMLAlpha', 1e6, 'DisplayMask', ...
'off', 'PlotScale', [-1, 1] * amp, 'Smooth', [false, true, true]};
% run the acoustic simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});