From 39e453eaf6e33599e2b3a25a378dc77ae05ec83c Mon Sep 17 00:00:00 2001 From: Julian Kosciessa Date: Thu, 8 Aug 2024 11:00:33 +0200 Subject: [PATCH] [bug] correct yaakub pCT conversion --- functions/setup_medium.m | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/functions/setup_medium.m b/functions/setup_medium.m index 4c3074c..0cdff83 100755 --- a/functions/setup_medium.m +++ b/functions/setup_medium.m @@ -104,17 +104,17 @@ case 'yaakub' % see https://github.com/sitiny/BRIC_TUS_Simulation_Tools/blob/main/tussim_skull_3D.m - c_min = 1500; % sound speed [m/s] - c_max = 3100; % max. speed of sound in skull (F. A. Duck, 2013.) [m/s] - rho_min = 1000; % density [kg/m^3] - rho_max = 1900; % max. skull density [kg/m3] + c_water = 1500; % sound speed [m/s] + c_skull = 3100; % max. speed of sound in skull (F. A. Duck, 2013.) [m/s] + rho_water = 1000; % density [kg/m^3] + rho_bone = 1900; % max. skull density [kg/m3] alpha_min = 4; % Aubry, J.-F. et al., 2022 (cortical bone; α (dB/cm at 500 kHz)) alpha_max = 8.7; % Fry 1978 at 0.5MHz: 1 Np/cm (8.7 dB/cm) HU_min = 300; % minimum HU considered as skull HU_max = 2000; % maximum skull HU for regularization % if observed max HU is lower than threshold, set HU_max to actual max - HU_max = max(max(pseudoCT(skull_idx),[],'all'), HU_max) + HU_max = max(max(pseudoCT(skull_idx),[],'all'), HU_max); % truncate CT HU (see Marsac et al., 2017) % do not use pCT-based properties for presumed non-skull tissue @@ -123,12 +123,13 @@ pseudoCT(pseudoCT > HU_max) = HU_max; % estimate density from CT HU based on Marsac et al., 2017 & Bancel et al., 2021 - density(skull_idx) = rho_min + (rho_max - rho_min) * ... - (pseudoCT(skull_idx) - 0) / (HU_max - 0); + % note: the original code hard-codes HU_min as 0, which may have been an error + density(skull_idx) = rho_water + (rho_bone - rho_water) * ... + (pseudoCT(skull_idx) - HU_min) / (HU_max - HU_min); % estimate sound speed - sound_speed(skull_idx) = c_min + (c_max - c_min) * ... - (medium.density - rho_min) / (rho_max - rho_min); + sound_speed(skull_idx) = c_water + (c_skull - c_water) * ... + (density(skull_idx) - rho_water) / (rho_bone - rho_water); end % estimate attenuation coefficients