I've read Bradley Treeby's recent nonlinear ultrasound simulation paper (https://asa.scitation.org/doi/full/10.1121/10.0002177) and I'm interested in acquiring the model of fUS skull penetration in the paper.
I read in the summary that the MATLAB version is available as part of the k-wave toolbox but I struggled to find it. Where is it?
Thanks in advance!
Hi milksoup,
The models are included with k-Wave as kspaceFirstOrderAS. Or is it the geometry of the skull (etc) that you're after?
Brad.
Hi Bradley,
It's the geometry of the skull actually! Would be great to access that.
Thank you!
Hi,
Where can I find the skull geometry?
Great, thank you!
Hi Brad,
I've tried to utilise your skull geometry in my own simulation utilising elements from the transducer field and axisymmetric examples. I made some changes to it so that it would be quicker to run (computational parameters, grid parameters, size of grid).
I'm having trouble making the kind of visualisation I want; something that shows the arc source and its produced waves going through the skull.
Just to experiment, I tried applying the visualisation scripts from some example - the only one that worked was from the transducer field example - but the results I got were not clear and I don't know how to modify the code to get the desired visual.
Any guidance would be appreciated.
Here is the code (I included things I commented out too):
clearvars;
% --------------------
% PROPERTIES
% --------------------
% source parameters
source_f0 = 1.1e6; % source frequency [Hz]
source_roc = 100e-3; % bowl radius of curvature [m]
source_diameter = 64e-3; % bowl aperture diameter [m]
source_mag = 2e6; % source pressure [Pa]
source_cycles = 4;
% material geometry
skull_radius = 80e-3;
skull_thickness = 6.5e-3;
skin_thickness = 6.5e-3;
skin_tx_distance = 30e-3;
% material properties
water_temp = 22;
water_c0 = waterSoundSpeed(water_temp);
water_rho0 = waterDensity(water_temp);
water_BonA = waterNonlinearity(water_temp);
water_a0 = waterAbsorption(source_f0 * 1e-6, water_temp);
skin_c0 = 1624;
skin_rho0 = 1109;
skin_BonA = 6.7;
skin_a0 = 1.84 / (source_f0 * 1e-6).^2;
brain_c0 = 1546;
brain_rho0 = 1045;
brain_BonA = 6.7;
brain_a0 = (0.59 * 1.1^1.3) / (source_f0 * 1e-6).^2;
skull_c0 = 2820;
skull_rho0 = 1732;
skull_BonA = 374;
skull_a0 = 2.7 / (source_f0 * 1e-6).^2;
% grid parameters (reduce by x10)
axial_size = 12e-3; % total grid size in the axial dimension [m]
lateral_size = 8e-3/2; % total grid size in the lateral dimension [m]
% computational parameters (reduce by x10)
ppw = 5; % number of points per wavelength
source_x_offset = 4; % grid points to offset the source
%% --------------------
% GRID
% --------------------
% calculate the grid spacing based on the PPW and F0
dx = water_c0 / (ppw * source_f0); % [m]
dy = dx;
% compute the size of the grid
Nx = roundEven(axial_size / dx) + source_x_offset;
Ny = roundEven(lateral_size / dx);
% My addition: forcing to approximate values with smaller primes and
% reduced by x10
Nx = 448;
Ny = 147;
% this didn't work at first because the issue was with the "PMLInside" set
% to false, which makes it so that it is added on to the dimentions
% later in the function kspaceFirstOrder_expandGridMatrices
% what I will do next is change PMLInside to true - although
% I wonder if I can keep it out and set its value.
% create the computational grid
kgrid = kWaveGrid(Nx, dx, Ny, dy);
%% --------------------
% MEDIUM
% --------------------
% convert geometry to gp
skull_radius_gp = skull_radius / dx;
skull_thickness_gp = skull_thickness / dx;
skin_thickness_gp = skin_thickness / dx;
skin_transducer_distance_gp = skin_tx_distance / dx;
% calculate centers for respective
skin_center = source_x_offset + skin_transducer_distance_gp + skin_thickness_gp + skull_radius_gp;
% create a skin map
skin_mask = makeDisc(Nx * 2, Ny * 6, skin_center, Ny + 1, skull_radius_gp + skin_thickness_gp);
% create a skull map
skull_outer_mask = makeDisc(Nx * 2, Ny * 6, skin_center, Ny + 1, skull_radius_gp);
skull_inner_mask = makeDisc(Nx * 2, Ny * 6, skin_center, Ny + 1, skull_radius_gp - skull_thickness_gp);
% trim masks
skin_mask = skin_mask(1:Nx, Ny + 1:2*Ny);
skull_outer_mask = skull_outer_mask(1:Nx, Ny + 1:2*Ny);
skull_inner_mask = skull_inner_mask(1:Nx, Ny + 1:2*Ny);
% reference sound speed
medium.sound_speed_ref = brain_c0;
% sound speed
medium.sound_speed = water_c0 * ones(Nx, Ny);
medium.sound_speed(skin_mask == 1) = skin_c0;
medium.sound_speed(skull_outer_mask == 1) = skull_c0;
medium.sound_speed(skull_inner_mask == 1) = brain_c0;
% sound speed
medium.density = water_rho0 * ones(Nx, Ny);
medium.density(skin_mask == 1) = skin_rho0;
medium.density(skull_outer_mask == 1) = skull_rho0;
medium.density(skull_inner_mask == 1) = brain_rho0;
% nonlinearity
medium.BonA = water_BonA * ones(Nx, Ny);
medium.BonA(skin_mask == 1) = skin_BonA;
medium.BonA(skull_outer_mask == 1) = skull_BonA;
medium.BonA(skull_inner_mask == 1) = brain_BonA;
% absorption
medium.alpha_coeff = water_a0 * ones(Nx, Ny);
medium.alpha_coeff(skin_mask == 1) = skin_a0;
medium.alpha_coeff(skull_outer_mask == 1) = skull_a0;
medium.alpha_coeff(skull_inner_mask == 1) = brain_a0;
%
%
% *(the following is from transducer field patterns example)*
% create the time array
kgrid.makeTime(medium.sound_speed);
% filter the source to remove high frequencies not supported by the grid
% source.p = filterTimeSeries(kgrid, medium, source.p);
%
%
% *(the following is from transducer field patterns example)*
% define a curved transducer element
arc_pos = [20, 20]; % [grid points]
radius = 60; % [grid points]
diameter = 81; % [grid points]
focus_pos = [Nx/2, Nx/2]; % [grid points]
source.p_mask = makeArc([Nx, Ny], arc_pos, radius, diameter, focus_pos);
% define a time varying sinusoidal source
source.p = source_mag * sin(2 * pi * source_f0 * kgrid.t_array);
%% --------------------
% MY SENSOR (from trans field patt ex)
% --------------------
% create a sensor mask covering the entire computational domain using the
% opposing corners of a rectangle
%sensor.mask = [1, 1, Nx, Ny].';
% create a display mask to display the transducer
display_mask = source.p_mask;
% from axisymmetric example
% define a Cartesian sensor mask with points in the shape of a circle
sensor.mask = makeCartCircle(40 * dx, 50);
% remove points from sensor mask where y < 0
sensor.mask(:, sensor.mask(2, :) < 0) = [];
% from transducer field example
% set the record mode capture the final wave-field and the statistics at
% each sensor point
sensor.record = {'p_final', 'p_max', 'p_rms'};
% assign the input options
% varargin = {'DisplayMask', display_mask, 'PMLInside', true, 'PlotPML', false};
% actually I'll change this PML inside to true so I don't have to worry
% about it being ---- actually actually this isnt important - ill remove
% varargin from the vars in kspacefirstorderAS function below
% run the simulation
sensor_data = kspaceFirstOrderAS(kgrid, medium, source, sensor);
%% --------------------
% PLOT
% --------------------
%
% % from kwave_test2
% % visualisation
% figure(43);
% kwave_vis(Nx,Ny,dx,dy,sensor_data,medium,'5mm',0);
% disp(['5 mm bone:' num2str(sensor_data.p_max([focus_pos])) ' Pa']);
% % from axisymmetric
% % create plot axis
% x_vec = 1e3 * kgrid.x_vec;
% y_vec = 1e3 * (kgrid.y_vec - kgrid.y_vec(1));
%
% % create the simulation layout, removing the PML
% pml_size = 20;
% sim_layout = double(source.p | cart2grid(kgrid, sensor.mask, true));
% sim_layout = sim_layout(1 + pml_size:end - pml_size, 1:end - pml_size);
%
% % plot the simulation layout
% figure;
% imagesc(y_vec, x_vec, sim_layout, [0, 1]);
% axis image;
% colormap(flipud(gray));
% xlabel('y (radial) position [mm]');
% ylabel('x (axial) position [mm]');
%
% % plot the simulated sensor data
% figure;
% imagesc(sensor_data, [-1, 1]);
% colormap(getColorMap);
% ylabel('Sensor Position');
% xlabel('Time Step');
% colorbar;
%
% % plot a trace from the simulated sensor data
% figure;
% plot(sensor_data(5, :));
% xlabel('Time Step');
% ylabel('Pressure [Pa]');
% %from axisymmetric
% from transducer field example
% add the source mask onto the recorded wave-field
sensor_data.p_final(source.p_mask ~= 0) = 1;imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.p_rms, [-1 1]);
sensor_data.p_max(source.p_mask ~= 0) = 1;
sensor_data.p_rms(source.p_mask ~= 0) = 1;
% plot the final wave-field
figure;
subplot(1, 3, 1);
imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.p_final, [-1 1]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('Final Wave Field');
% plot the maximum recorded pressure
subplot(1, 3, 2);
imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.p_max, [-1 1]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('Maximum Pressure');
% plot the rms recorded pressure
subplot(1, 3, 3);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('RMS Pressure');
scaleFig(2, 1);
% plot the material property maps
figure;
subplot(1, 4, 1);
imagesc(kgrid.y_vec - kgrid.y_vec(1), kgrid.x_vec - kgrid.x_vec(1), medium.sound_speed);
axis image;
colorbar;
title('Sound Speed');
subplot(1, 4, 2);
imagesc(kgrid.y_vec - kgrid.y_vec(1), kgrid.x_vec - kgrid.x_vec(1), medium.density);
axis image;
colorbar;
title('Density');
subplot(1, 4, 3);
imagesc(kgrid.y_vec - kgrid.y_vec(1), kgrid.x_vec - kgrid.x_vec(1), medium.BonA);
axis image;
colorbar;
title('BonA');
subplot(1, 4, 4);
imagesc(kgrid.y_vec - kgrid.y_vec(1), kgrid.x_vec - kgrid.x_vec(1), medium.alpha_coeff);
axis image;
colorbar;
title('\alpha_0');
scaleFig(1.5, 1);I'm not exactly sure what you mean by the desired visualisation, but keep in mind these simulations are axisymmetric, thus your source definition should be axisymmetric (in other words, half an arc positioned along the edge of the domain). There are some axisymmetric examples in the new kWaveArray class that you could take a look at for reference.
You must log in to post.