k-Wave User Forum » Topic: Simulating bony skull attenuation- why is it amplifying instead?

k-Wave User Forum » Topic: Simulating bony skull attenuation- why is it amplifying instead? http://www.k-wave.org/forum/topic/simulating-bony-skull-attenuation-why-is-it-amplifying-instead Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:39:17 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Simulating bony skull attenuation- why is it amplifying instead?" http://www.k-wave.org/forum/topic/simulating-bony-skull-attenuation-why-is-it-amplifying-instead#post-6736 Tue, 29 Jan 2019 23:24:32 +0000 Bradley Treeby 6736@http://www.k-wave.org/forum/ Hi an803355, I can't immediately spot anything wrong in your code. Do you mean the focal pressure is higher when you include attenuation compared to without? Or that the focal pressure is higher when you include the skull compared to without? Or something else? Keep in mind the focal pressure will always be much higher than the surface pressure due to the focusing gain of the transducer. For single frequency simulations, I would suggest using <code>alpha_coeff = 2</code>, and then adjust the alpha_coeff values accordingly. One way is to use the new function <a href="http://www.k-wave.org/downloads/fitPowerLawParamsMulti.m"><code>fitPowerLawParamsMulti</code></a> (this will be included in a future release of k-Wave). The reason for using <code>alpha_coeff = 2</code> is explained in Sec. 2.2 in <a href="http://bug.medphys.ucl.ac.uk/papers/2018-Robertson-PMB.pdf">this paper</a>. If you're still stuck, can you post an example with an artificial skull that has the same problem, or email me the skull image? Brad an803355 on "Simulating bony skull attenuation- why is it amplifying instead?" http://www.k-wave.org/forum/topic/simulating-bony-skull-attenuation-why-is-it-amplifying-instead#post-6716 Thu, 17 Jan 2019 17:39:01 +0000 an803355 6716@http://www.k-wave.org/forum/ Bradley, here is my code. The image mask is not included, but essentially it is black and white skull curved strip. The masks apply well (I've checked the matrix). Not sure why the amplitude isn't attenuating at the focus though, and I'm getting really weird behaviour within the skull (seems to be amplifying everything) instead of attenuating Thanks for your help clearvars; % ========================================================================= % ACOUSTIC SIMULATION % ========================================================================= % define the PML size pml_size = 20; % [grid points] % define the grid parameters Nx = 256 - 2 * pml_size; % [grid points] Ny = 256 - 2 * pml_size; % [grid points] dx = 0.25e-3; % [m] dy = 0.25e-3; % [m] % create the computational grid kgrid = kWaveGrid(Nx, dx, Ny, dy); %define the properties of the propagation medium %% ALPHA COEFF medium.alpha_coeff = 0.05 * ones(Nx, Ny); % [dB/(MHz^y cm)] WATER %skull pix=loadimagecustom('skull216.png'); pix = pix .* (8/ max(pix(:))); %set all black pixels to 8, all white to 0. medium.alpha_coeff =medium.alpha_coeff + pix; % add 8 to the 0.05 water, making the skull layer medium.alpha_power = 1.43; %% SPEED medium.sound_speed = 1450 * ones(Nx, Ny); % WATER default fill! [m/s] % %skull pix=loadimagecustom('skull216.png'); pix = pix .* (700/ max(pix(:))); *add 700 to the water speed medium.sound_speed =medium.sound_speed +pix; % 2500 bone! [kg/m^3] %% DENSITY medium.density = 1000 * ones(Nx, Ny); % WATER default fill [kg/m^3] %skull pix=loadimagecustom('skull216.png'); pix = pix .* (732 / max(pix(:))); medium.density =medium.density +pix; % bone! [kg/m^3] %% define the source parameters diameter = 0.03; % [m] radius = 0.03; % [m] freq = 500000; % [Hz] amp = 0.53e6; % 0.53[MPa] % define a focused ultrasound transducer source.p_mask = makeArc([Nx, Ny], [1, Ny/2], round(radius / dx), round(diameter / dx) + 1, [Nx/2, Ny/2]); % create a display mask to display the transducer display_mask = source.p_mask; % calculate the time step using an integer number of points per period ppw = 12.08; % points per wavelength cfl = 0.3; % cfl number ppp = 41; % points per period T = 1 / freq; % period [s] dt = 4.878e-8; % time step [s] % calculate the number of time steps to reach steady state t_end = 5.0575e-5; 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); % set the sensor mask to cover the entire grid sensor.mask = ones(Nx, Ny); sensor.record = {'p', 'p_max_all'}; % 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,'PlotLayout', true, 'PlotPML', false, 'DisplayMask', ... display_mask, 'PlotScale', [-1, 1] * amp}; % run the acoustic simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % ========================================================================= % CALCULATE PRESSURE %=================== % extract the pressure amplitude at each position p = extractAmpPhase(sensor_data.p, 1/kgrid.dt, freq); % reshape the data, p = reshape(p, Nx, Ny); % ========================================================================= % VISUALISATION % ========================================================================= % plot the thermal dose and lesion map figure; % plot the acoustic pressure subplot(2, 3, 1); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, p * 1e-6); h = colorbar; xlabel(h, '[MPa]'); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Acoustic Pressure Amplitude'); set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]); % set colormap and enlarge figure window colormap(jet(256)); scaleFig(1.5, 1); an803355 on "Simulating bony skull attenuation- why is it amplifying instead?" http://www.k-wave.org/forum/topic/simulating-bony-skull-attenuation-why-is-it-amplifying-instead#post-6715 Thu, 17 Jan 2019 17:31:56 +0000 an803355 6715@http://www.k-wave.org/forum/ Hi @Brad. Yes that was a typo. I am setting my coefficient to about +8.00. Nevertheless, it is not adequately running. Question: Do I need to set a different "alpha power" for my skull mask? Or does the alpha power stay as, say 1.43, for the entire field? I'm a bit confused on this point. ] When I try to run " fitPowerLawParams(8, 1.43, 2850, 400000, 500000, true) ", the function gives me an error "not enough input arguments". A bit confused about the whole alpha power. Bradley Treeby on "Simulating bony skull attenuation- why is it amplifying instead?" http://www.k-wave.org/forum/topic/simulating-bony-skull-attenuation-why-is-it-amplifying-instead#post-6712 Thu, 17 Jan 2019 00:26:45 +0000 Bradley Treeby 6712@http://www.k-wave.org/forum/ Are you using a negative absorption coefficient, or is that a typo? If so, that would cause an increase in pressure. Brad. an803355 on "Simulating bony skull attenuation- why is it amplifying instead?" http://www.k-wave.org/forum/topic/simulating-bony-skull-attenuation-why-is-it-amplifying-instead#post-6711 Wed, 16 Jan 2019 21:30:35 +0000 an803355 6711@http://www.k-wave.org/forum/ Dear all, I am trying to simulate a focused 0.5 MHz transducer directed at the inside of a skull (initial pressure 0.5 MPa). Of course, the skull will attenuate the pressure significantly, but in my code, it seems to make it stronger at the focus! I am using parameters from this paper: <a href="http://bug.medphys.ucl.ac.uk/papers/2017-Robertson-PMB.pdf" rel="nofollow">http://bug.medphys.ucl.ac.uk/papers/2017-Robertson-PMB.pdf</a> Specifically, I am using a 2D image mask to define a skull bone layer and the following parameters: medium.alpha_coeff = 0.75 for WATER medium.alpha_coeff= -8.83 for SKULL (0.8 cm thick) medium.alpha_power = 1.43 I am using standard values for bone speed 2850 m/s and density was 1732 kg/m 3 For some reason, my simulation is generating a larger value at the focus (1 MPA), even though the wave had to pass through the skull layer! I am confused about what I am doing wrong. Any help would be appreciated. The mask appears to be applied properly. I suspect my coefficients are off or I used the wrong sensor mask method? PS: I used the "Heating By A Focused Ultrasound Transducer" as a template script, removed the heating parts and left the "pressure" simulation there, with the addition of a skull.