Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Documentation of pseudoCT processing #43

Open
jkosciessa opened this issue Jun 11, 2024 · 6 comments
Open

Documentation of pseudoCT processing #43

jkosciessa opened this issue Jun 11, 2024 · 6 comments
Labels
documentation Improvements or additions to documentation

Comments

@jkosciessa
Copy link
Collaborator

Now that we have pseudoCT pseudo-integrated, it would be worthwhile to document how we are processing the related tissue mask (e.g., smoothing choices) and the conversion into pseudo-Hounsfield, sound_speed and alpha units.

The relevant functions are:

@jkosciessa
Copy link
Collaborator Author

Setting up kwave medium...
Setting up kwave source...
Check stability...
Starting acoustic simulations...
Running k-Wave simulation...
  start time: 07-Aug-2024 12:13:11
  reference sound speed: 3482.4578m/s
  dt: 21.2766ns, t_end: 108.4043us, time steps: 5096
  input grid size: 486 by 432 by 384 grid points (243 by 216 by 192mm)
  maximum supported frequency: 168.33kHz

I am trying to debug a pseudoCT run and noticed that my max. supported freq is below the trageted 500 kHz. Anyone run into this before? The relevant computation is:

c_min = min(medium.sound_speed(:)) = 168.33
kgrid.k_max = 6.28e3
        disp(['  maximum supported frequency: ' scaleSI( kgrid.k_max .* c_min ./ (2*pi) ) 'Hz']);         

The reference sound speed in this case is above what we usually define given that it is computed from the pseudo-HU units.

This is the usual spec I get (which in this case also produces NaNs however):

reference sound speed: 2800m/s 
dt: 26.3158ns, t_end: 134.8158us, time steps: 5124
  input grid size: 486 by 432 by 384 grid points (243 by 216 by 192mm)
  maximum supported frequency: 1.5MHz

I presume it's because the minimum sound speed is unreasonable (below that of water), so a problem with the CT-to-tissueproperty transform.

@jkosciessa
Copy link
Collaborator Author

Update: determining the lower bound of soundspeed as the one of water indeed fixes at least the max. supported frequency issue. The relevant fix is bdf0ff9. Given that soundspeed is directly translated from the estimated density this may only be a superficial fix for a deeper problem.

@jkosciessa
Copy link
Collaborator Author

The pCT-to-tissue conversion algorithm documented in Yaakub et al. is now implemented in
f8d7286 and can be called by specifying parameters.pseudoCT_variant = 'yaakub'. Funcitonality and con-/divercence is not yet tested.

Ref: Yaakub, S. N. et al. Pseudo-CTs from T1-Weighted MRI for Planning of Low-Intensity Transcranial Focused Ultrasound Neuromodulation: An Open-Source Tool. Brain Stimulation. 16. 75–78 (2023).

@jkosciessa
Copy link
Collaborator Author

jkosciessa commented Aug 7, 2024

For parameters.pseudoCT_variant = 'carpino':

All following steps are applied to the skull mask only.

  • The k-Wave function hounsfield2density converts pseudo-HUs to density.
  • Original pseudoCT values are initially shifted by 1000 and thresholded at 0 to exclude air voxels.
  • Resulting density values are regularized to a minimum of the specified water density [originally: 1].
  • The sound speed c is calculated from density 𝜌 using the linear relationship: c = 1.33𝜌 + 167. This is identical to the k-Plan estimation.
  • The absorption coefficient 𝛼 is derived from HU values according to formula (3) in Yaakub et al; the offset of 1000 for pseudo-HUs is removed prior to entering the formula.

For both variants, 𝛼 is bounded based on estimates made at 500 kHz (i.e., 𝛼(f); see Aubry, J.-F., 2022 for prior benchmark simulations). We actually require alpha_0 in 𝛼(f) = alpha_0 x f[MHz] ^ y. We estimate alpha_0 = 𝛼(f)/0.5^y with y being the specified alpha_power_truefor the skull tissue.

Reference: Adapted from Carpino et al. (2024). Transcranial ultrasonic stimulation of the human amygdala to modulate threat learning. MSc thesis.

@jkosciessa
Copy link
Collaborator Author

jkosciessa commented Aug 8, 2024

parameters.pseudoCT_variant = 'carpino' still runs as expected and is the default.

For yaakub the original implementation of formula (1) assumes/hard-codes HU_min as 0. I am not sure this is intentional, given that HU_min will be a minimum of 300 due to the regularization.

I have set it to HU_min for now with commit 39e453e.

@jkosciessa
Copy link
Collaborator Author

jkosciessa commented Aug 8, 2024

The current yaakub settings would turn much of the tribecular skull bone into water. Overall, our pseudoCT images appear to have many values < 300 (and even below 0). This could be due to diverging mapping functions or UTE acquisition parameters. Other than that, yaakuub now also seems to run fine.

Screenshot 2024-08-08 at 11 27 54

@jkosciessa jkosciessa added the documentation Improvements or additions to documentation label Sep 30, 2024
@jkosciessa jkosciessa pinned this issue Sep 30, 2024
jkosciessa added a commit that referenced this issue Oct 1, 2024
jkosciessa added a commit that referenced this issue Oct 1, 2024
jkosciessa added a commit that referenced this issue Oct 3, 2024
jkosciessa added a commit that referenced this issue Oct 3, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
documentation Improvements or additions to documentation
Projects
None yet
Development

No branches or pull requests

1 participant