[bug] correct yaakub pCT conversion

jkosciessa · Aug 8, 2024 · 39e453e · 39e453e
1 parent f8d7286
commit 39e453e
Showing 1 changed file with 10 additions and 9 deletions.
diff --git 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