http://www.k-wave.org/forum/topic/treebys-recent-nonlinear-ultrasound-simulation-paper-and-model Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:33:30 +0000 http://bbpress.org/?v=1.0.2 q http://www.k-wave.org/forum/search.php http://www.k-wave.org/forum/topic/treebys-recent-nonlinear-ultrasound-simulation-paper-and-model#post-7974 Thu, 03 Dec 2020 10:04:44 +0000 Bradley Treeby 7974@http://www.k-wave.org/forum/ <p>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 <a href="http://www.k-wave.org/forum/topic/alpha-version-of-kwavearray-off-grid-sources"><code>kWaveArray</code></a> class that you could take a look at for reference. </p> http://www.k-wave.org/forum/topic/treebys-recent-nonlinear-ultrasound-simulation-paper-and-model#post-7967 Sun, 29 Nov 2020 21:32:36 +0000 milksoup 7967@http://www.k-wave.org/forum/ <p>Hi Brad, </p> <p>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). </p> <p>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. </p> <p>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. </p> <p>Any guidance would be appreciated.</p> <p>Here is the code (I included things I commented out too): </p> <pre><code>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&#39;t work at first because the issue was with the &quot;PMLInside&quot; 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].&#39;; % 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 &lt; 0 sensor.mask(:, sensor.mask(2, :) &lt; 0) = []; % from transducer field example % set the record mode capture the final wave-field and the statistics at % each sensor point sensor.record = {&#39;p_final&#39;, &#39;p_max&#39;, &#39;p_rms&#39;}; % assign the input options % varargin = {&#39;DisplayMask&#39;, display_mask, &#39;PMLInside&#39;, true, &#39;PlotPML&#39;, false}; % actually I&#39;ll change this PML inside to true so I don&#39;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,&#39;5mm&#39;,0); % disp([&#39;5 mm bone:&#39; num2str(sensor_data.p_max([focus_pos])) &#39; Pa&#39;]); % % 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(&#39;y (radial) position [mm]&#39;); % ylabel(&#39;x (axial) position [mm]&#39;); % % % plot the simulated sensor data % figure; % imagesc(sensor_data, [-1, 1]); % colormap(getColorMap); % ylabel(&#39;Sensor Position&#39;); % xlabel(&#39;Time Step&#39;); % colorbar; % % % plot a trace from the simulated sensor data % figure; % plot(sensor_data(5, :)); % xlabel(&#39;Time Step&#39;); % ylabel(&#39;Pressure [Pa]&#39;); % %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(&#39;x-position [mm]&#39;); xlabel(&#39;y-position [mm]&#39;); axis image; title(&#39;Final Wave Field&#39;); % 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(&#39;x-position [mm]&#39;); xlabel(&#39;y-position [mm]&#39;); axis image; title(&#39;Maximum Pressure&#39;); % plot the rms recorded pressure subplot(1, 3, 3); colormap(getColorMap); ylabel(&#39;x-position [mm]&#39;); xlabel(&#39;y-position [mm]&#39;); axis image; title(&#39;RMS Pressure&#39;); 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(&#39;Sound Speed&#39;); 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(&#39;Density&#39;); 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(&#39;BonA&#39;); 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(&#39;\alpha_0&#39;); scaleFig(1.5, 1);</code></pre> http://www.k-wave.org/forum/topic/treebys-recent-nonlinear-ultrasound-simulation-paper-and-model#post-7956 Thu, 26 Nov 2020 18:41:31 +0000 milksoup 7956@http://www.k-wave.org/forum/ <p>Great, thank you! </p> http://www.k-wave.org/forum/topic/treebys-recent-nonlinear-ultrasound-simulation-paper-and-model#post-7950 Thu, 26 Nov 2020 13:02:11 +0000 Bradley Treeby 7950@http://www.k-wave.org/forum/ <p>I've uploaded the code used to create the geometry <a href="http://www.k-wave.org/downloads/forum_axisym_skull_geom.m">here</a>.</p> <p>Brad </p> http://www.k-wave.org/forum/topic/treebys-recent-nonlinear-ultrasound-simulation-paper-and-model#post-7928 Sat, 21 Nov 2020 21:56:44 +0000 milksoup 7928@http://www.k-wave.org/forum/ <p>Hi,</p> <p>Where can I find the skull geometry? </p> http://www.k-wave.org/forum/topic/treebys-recent-nonlinear-ultrasound-simulation-paper-and-model#post-7913 Tue, 03 Nov 2020 04:19:28 +0000 milksoup 7913@http://www.k-wave.org/forum/ <p>Hi Bradley,</p> <p>It's the geometry of the skull actually! Would be great to access that. </p> <p>Thank you! </p> http://www.k-wave.org/forum/topic/treebys-recent-nonlinear-ultrasound-simulation-paper-and-model#post-7905 Sun, 01 Nov 2020 22:42:58 +0000 Bradley Treeby 7905@http://www.k-wave.org/forum/ <p>Hi milksoup,</p> <p>The models are included with k-Wave as <code>kspaceFirstOrderAS</code>. Or is it the geometry of the skull (etc) that you're after?</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/treebys-recent-nonlinear-ultrasound-simulation-paper-and-model#post-7893 Thu, 29 Oct 2020 17:46:35 +0000 milksoup 7893@http://www.k-wave.org/forum/ <p>I've read Bradley Treeby's recent nonlinear ultrasound simulation paper (<a href="https://asa.scitation.org/doi/full/10.1121/10.0002177" rel="nofollow">https://asa.scitation.org/doi/full/10.1121/10.0002177</a>) and I'm interested in acquiring the model of fUS skull penetration in the paper. </p> <p>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?</p> <p>Thanks in advance! </p>