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);