k-Wave User Forum » Topic: HIFU Wave through a heterogenous medium

k-Wave User Forum » Topic: HIFU Wave through a heterogenous medium http://www.k-wave.org/forum/topic/hifu-wave-through-a-heterogenous-medium Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:32:55 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php scmmw on "HIFU Wave through a heterogenous medium" http://www.k-wave.org/forum/topic/hifu-wave-through-a-heterogenous-medium#post-9242 Tue, 13 Jan 2026 16:35:52 +0000 scmmw 9242@http://www.k-wave.org/forum/ Hi, I'm trying to simulate a HIFU wave propagating through a heterogeneous medium (water, soft tissue, muscle, thyroid, and nodule). My simulation works as expected at 2 MHz, producing a lesion at the focus point. However, at 3 MHz, side lesions appear instead of a lesion at the focus. I'm unsure how to resolve this issue. Any help would be greatly appreciated! Here is my script (I have a lot of visualisations, this is just to help me understand the process) clearvars; % ========================================================================= % ACOUSTIC SIMULATION % ========================================================================= %% define the PML size pml_size = 20; % [grid points] % define the grid parameters Nx = 524 - 2 * pml_size; % [grid points] Ny = 427 - 2 * pml_size; % [grid points] dx = 206e-6; % [m] dy = 206e-6; % [m] % create the computational grid kgrid = kWaveGrid(Nx, dx, Ny, dy); % define the properties of the propagation medium medium.sound_speed = zeros(Nx,Ny); % [m/s] medium.density = zeros(Nx,Ny); % [kg/m^3] medium.alpha_coeff = zeros(Nx,Ny); % [dB/(MHz^y cm)] water.sound_speed = 1480; water.density = 1000; % [kg/m^3] water.alpha_coeff = 0.0022; % [dB/(MHz^y cm)] water.alpha_power = 1; tissue.sound_speed = 1540; tissue.density = 1020; % [kg/m^3] tissue.alpha_coeff = 0.75; % [dB/(MHz^y cm)] tissue.alpha_power = 1.5; muscle.sound_speed = 1588; muscle.density = 1090; % [kg/m^3] muscle.alpha_coeff = 0.63; % [dB/(MHz^y cm)] muscle.alpha_power = 1.3; thyroid.sound_speed = 1585; thyroid.density = 1050; % [kg/m^3] thyroid.alpha_coeff = 1.34; % [dB/(MHz^y cm)] thyroid.alpha_power = 1.5; nodule.sound_speed = 1583; nodule.density = 1080; % [kg/m^3] nodule.alpha_coeff = 0.31; % [dB/(MHz^y cm)] nodule.alpha_power = 1.1; % now break down the grid % --- Water region ------------------------------------------------------- medium.sound_speed(1:156, :) = water.sound_speed; medium.density(1:156, :) = water.density; medium.alpha_coeff(1:156, :) = water.alpha_coeff; % --- Tissue region ------------------------------------------------------- medium.sound_speed(157:341, :) = tissue.sound_speed; medium.density(157:341, :) = tissue.density; medium.alpha_coeff(157:341, :) = tissue.alpha_coeff; %medium.alpha_power(1:185, :) = tissue.alpha_power; % --- Muscle region ------------------------------------------------------- medium.sound_speed(342:416, :) = muscle.sound_speed; medium.density(342:416, :) = muscle.density; medium.alpha_coeff(342:416, :) = muscle.alpha_coeff; %medium.alpha_power(186:260, :) = muscle.alpha_power; % --- Thyroid region ------------------------------------------------------ medium.sound_speed(417:427, :) = thyroid.sound_speed; medium.density(417:427, :) = thyroid.density; medium.alpha_coeff(417:427, :) = thyroid.alpha_coeff; %medium.alpha_power(261:270, :) = thyroid.alpha_power; % --- Nodule region ------------------------------------------------------- medium.sound_speed(427:end, :) = nodule.sound_speed; medium.density(427:end, :) = nodule.density; medium.alpha_coeff(427:end, :) = nodule.alpha_coeff; %medium.alpha_power(271:end, :) = nodule.alpha_power; medium.sound_speed = smooth(medium.sound_speed, true); medium.density = smooth(medium.density, true); medium.alpha_coeff = smooth(medium.alpha_coeff, true); % define the source parameters diameter = 64e-3; % [m] radius = 63.2e-3; % [m] aperture = 20e-3; % [m] freq = 2e6; % [Hz] amp = 0.5e6; % [Pa] %% This may need running a few times to optimise alphafit to get best fit power value i.e. lines overlap as best as possible. alpha_coeff_vec = [water.alpha_coeff, tissue.alpha_coeff, muscle.alpha_coeff, thyroid.alpha_coeff, nodule.alpha_coeff]; alpha_power_vec = [water.alpha_power,tissue.alpha_power, muscle.alpha_power, thyroid.alpha_power, nodule.alpha_power]; c_vec = [water.sound_speed, tissue.sound_speed, muscle.sound_speed, thyroid.sound_speed, nodule.sound_speed]; alphaFit = 1.5; alpha_coeff_fit_vec = fitPowerLawParamsMulti(alpha_coeff_vec, alpha_power_vec, c_vec, freq, alphaFit, true); medium.alpha_power = alphaFit; %% % --- Tissue --- medium.alpha_coeff(1:156, :) = alpha_coeff_fit_vec(1); % --- Tissue --- medium.alpha_coeff(157:341, :) = alpha_coeff_fit_vec(2); % --- Muscle --- medium.alpha_coeff(342:416, :) = alpha_coeff_fit_vec(3); % --- Thyroid --- medium.alpha_coeff(417:427, :) = alpha_coeff_fit_vec(4); % --- Nodule --- medium.alpha_coeff(427:end, :) = alpha_coeff_fit_vec(5); %% close all source.p_mask = zeros(Nx,Ny); % define a focused ultrasound transducer null_mask = makeArc([Nx, Ny], [1, Ny/2], round(radius / dx), round(aperture / dx), [Nx/2, Ny/2]); % define aperture in transducer source.p_mask = source.p_mask + makeArc([Nx, Ny], [1, Ny/2], round(radius / dx), round(diameter / dx), [Nx/2, Ny/2]) - null_mask; % calculate the time step using an integer number of points per period ppw = mean(mean(medium.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 ) / mean(mean(medium.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); % 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, 'PlotPML', false, 'DisplayMask', ... 'off', 'PlotScale', [-1, 1] * amp}; % run the acoustic simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % ========================================================================= % CALCULATE HEATING % ========================================================================= % convert the absorption coefficient to nepers/m alpha_np = db2neper(medium.alpha_coeff, medium.alpha_power) * ... (2 * pi * freq).^medium.alpha_power; % extract the pressure amplitude at each position p = extractAmpPhase(sensor_data.p, 1/kgrid.dt, freq); % reshape the data, and calculate the volume rate of heat deposition p = reshape(p, Nx, Ny); Q = alpha_np .* p.^2 ./ (medium.density .* medium.sound_speed); % Find the location of maximum heat deposition [maxQ, idxMaxQ] = max(Q(:)); [rowMax, colMax] = ind2sub(size(Q), idxMaxQ); % Place the thermal sensor at that location sensor.mask = zeros(Nx, Ny); sensor.mask(rowMax, colMax) = 1; % Report coordinates fprintf('Thermal sensor placed at row=%d, col=%d (x=%.3f m, y=%.3f m)\n', ... rowMax, colMax, kgrid.x_vec(rowMax), kgrid.y_vec(colMax)); % ========================================================================= % THERMAL SIMULATION % ========================================================================= % clear the input structures %clear source sensor; %medium % set the background temperature and heating term source.Q = Q; source.T0 = 37; medium.thermal_conductivity = zeros(Nx,Ny); % [W/(m.K)] medium.specific_heat = zeros(Nx,Ny); % [J/(kg.K)] water.thermal_conductivity = 0.598; % [W/(m.K)] water.specific_heat = 4200; % [J/(kg.K)] tissue.thermal_conductivity = 0.5; % [W/(m.K)] tissue.specific_heat = 3600; % [J/(kg.K)] muscle.thermal_conductivity = 0.49; muscle.specific_heat = 3421; thyroid.thermal_conductivity = 0.5; thyroid.specific_heat = 3826; nodule.thermal_conductivity = 0.5; nodule.specific_heat = 3826; % --- Water region ------------------------------------------------------- medium.thermal_conductivity(1:156, :) = water.thermal_conductivity; medium.specific_heat(1:156, :) = water.specific_heat; % --- Tissue region ------------------------------------------------------- medium.thermal_conductivity(157:341, :) = tissue.thermal_conductivity; medium.specific_heat(157:341, :) = tissue.specific_heat; % --- Muscle region ------------------------------------------------------- medium.thermal_conductivity(342:416, :) = muscle.thermal_conductivity; medium.specific_heat(342:416, :) = muscle.specific_heat; % --- Thyroid region ------------------------------------------------------ medium.thermal_conductivity(417:427, :) = thyroid.thermal_conductivity; medium.specific_heat(417:427, :) = thyroid.specific_heat; % --- Nodule region ------------------------------------------------------- medium.thermal_conductivity(427:end, :) = nodule.thermal_conductivity; medium.specific_heat(427:end, :) = nodule.specific_heat; % create kWaveDiffusion object kdiff = kWaveDiffusion(kgrid, medium, source, sensor); % set source on time and off time on_time = 10; % [s] off_time = 200; % [s] % set time step size dt = 0.1; % Number of timesteps for the entire heating period num_steps = round(on_time / dt); % Initialize a container for CEM43 values at each timestep cem43_values = zeros(1, num_steps); % Initialize variable to store the time when CEM43 exceeds 240 time_cem43_exceeded = NaN; % Default to NaN (not exceeded yet) % Iterate over each timestep for step = 1:num_steps % Take one timestep kdiff.takeTimeStep(1, dt); % Record CEM43 value at the sensor location cem43_values(step) = max(kdiff.cem43(sensor.mask == 1)); % Assuming you want the max CEM43 at the sensor position % Display the CEM43 value and time %fprintf('Time: %.2f s, CEM43: %.2f\n', step * dt, cem43_values(step)); % Check if CEM43 exceeds 240 if cem43_values(step) >= 240 && isnan(time_cem43_exceeded) time_cem43_exceeded = step * dt; % Record the time when it exceeds 240 fprintf('CEM43 exceeded 240 at time %.2f s.\n', time_cem43_exceeded); end end % After the loop, you can use time_cem43_exceeded to check when the threshold was met if ~isnan(time_cem43_exceeded) fprintf('CEM43 exceeded 240 at %.2f seconds during the simulation.\n', time_cem43_exceeded); else disp('CEM43 did not exceed 240 during the simulation.'); end % store the current temperature field T1 = kdiff.T; if max(max(T1)) >= 100 disp('Boiling occured'); end % turn off heat source and take time steps kdiff.Q = 0; kdiff.takeTimeStep(round(off_time / dt), dt); % store the current temperature field T2 = kdiff.T; [max_dose, idx_dose] = max(kdiff.cem43(:)); % Find max thermal dose and index [row_dose, col_dose] = ind2sub(size(kdiff.cem43), idx_dose); % Get coordinates % Get the x and y coordinates of the maximum thermal dose center_x_dose = kgrid.x_vec(row_dose); center_y_dose = kgrid.y_vec(col_dose); % Display the result disp(['Center of Thermal Dose: x = ', num2str(center_x_dose), ' m, y = ', num2str(center_y_dose), ' m']); disp(['Center of Thermal Dose: x = ', num2str(row_dose), ' , y = ', num2str(col_dose), ' ']); %% break % ========================================================================= % VISUALISATION % ========================================================================= % Define the intensity threshold for ablated tissue threshold = 0.99; % Values close to 1 indicate ablated tissue % Create a binary mask of the ablated region binary_mask = kdiff.lesion_map >= threshold; % Find the rows and columns where the ablated tissue exists [row_idx, col_idx] = find(binary_mask); % Add half the grid spacing to account for pixel edges min_x = min(kgrid.x_vec(row_idx)) - abs(kgrid.x_vec(2) - kgrid.x_vec(1)) / 2; max_x = max(kgrid.x_vec(row_idx)) + abs(kgrid.x_vec(2) - kgrid.x_vec(1)) / 2; min_y = min(kgrid.y_vec(col_idx)) - abs(kgrid.y_vec(2) - kgrid.y_vec(1)) / 2; max_y = max(kgrid.y_vec(col_idx)) + abs(kgrid.y_vec(2) - kgrid.y_vec(1)) / 2; % Compute the adjusted dimensions length_ablated = (max_x - min_x) * 1e3; % in mm width_ablated = (max_y - min_y) * 1e3; % in mm % Display the results fprintf('Ablated Tissue Dimensions:\n'); fprintf('Length: %.2f mm\n', length_ablated); fprintf('Width: %.2f mm\n', width_ablated); % Optional: Visualize the lesion map with bounding box figure; imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, kdiff.lesion_map, [0, 1]); colorbar; xlabel('y-position [mm]'); ylabel('x-position [mm]'); axis image; title('Ablated Tissue'); % Use adjusted bounds in the rectangle function rectangle('Position', [min_y * 1e3, min_x * 1e3, ... width_ablated, length_ablated], ... 'EdgeColor', 'r', 'LineWidth', 1.5); % Create time axis for the recorded data t_axis = (0:kdiff.time_steps_taken - 1) * dt; % Plot temperature at every grid point (optional: you can visualize data at a specific point) figure; imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, kdiff.sensor_data(:,:,end)); % Plot temperature at last time step colorbar; xlabel('y-position [mm]'); ylabel('x-position [mm]'); title('Temperature at Last Time Step'); % Plot temperature time profile figure; plot(t_axis, kdiff.sensor_data, 'LineWidth', 1.5); % Enhanced line width for better visibility xlabel('Time [s]'); ylabel('Temperature [^\circC]'); % Add grid lines grid on; % Turns on the grid % Customize y-axis ticks ylim([floor(min(kdiff.sensor_data(:))), ceil(max(kdiff.sensor_data(:)))]); % Ensure y-axis limits encompass the data yticks(floor(min(kdiff.sensor_data(:))):1:ceil(max(kdiff.sensor_data(:)))); % Set ticks with a step of 1 % Optionally add a title title('Temperature Time Profile'); % 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'); % plot the volume rate of heat deposition subplot(2, 3, 2); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, Q * 1e-7); h = colorbar; xlabel(h, '[kW/cm^2]'); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Volume Rate Of Heat Deposition'); % plot the temperature after heating subplot(2, 3, 3); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, T1); h = colorbar; xlabel(h, '[degC]'); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Temperature After Heating'); % plot the temperature after cooling subplot(2, 3, 4); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, T2); h = colorbar; xlabel(h, '[degC]'); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Temperature After Cooling'); % plot thermal dose subplot(2, 3, 5); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, kdiff.cem43, [0, 1000]); h = colorbar; xlabel(h, '[CEM43]'); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Thermal Dose'); % plot lesion map subplot(2, 3, 6); imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, kdiff.lesion_map, [0, 1]); colorbar; ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('Ablated Tissue'); % set colormap and enlarge figure window colormap(jet(256)); scaleFig(1.5, 1);