Skip to content

Commit

Permalink
Clarified code for computing B1+ maps
Browse files Browse the repository at this point in the history 
Fixes #22
  • Loading branch information
@jcohenadad
jcohenadad committed Jan 8, 2024
1 parent ec36eed commit 9f50c12
Showing 1 changed file with 69 additions and 16 deletions.
85 changes: 69 additions & 16 deletions data_processing.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -254,20 +254,65 @@
" !sct_apply_transfo -i ../anat/{subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -d {subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -w warp_{subject}_acq-CoV_T2starw_crop2{subject}_acq-anat{shim_mode}_TB1TFL.nii.gz -x nn -o {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "87a06e2d",
"metadata": {},
"outputs": [],
"source": [
"subject='sub-01'\n",
"os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
"shim_mode = 'CoV'\n",
"\n",
"# Convert the flip angle maps to B1+ efficiency maps [nT/V]\n",
"# % Calculate the B1eff using a 1ms, pi-pulse at the acquisition voltage,\n",
"# % then scale the efficiency by the ratio of the measured flip angle to the\n",
"# % requested flip angle in the pulse sequence.\n",
"GAMMA = 2.675e8; # [rad / (s T)]\n",
"requested_fa = 90 # saturation flip angle -- hard-coded in sequence\n",
"\n",
"\n",
"# Fetch the reference voltage from the JSON sidecar to the TFL B1map sequence\n",
"with open(f\"{subject}_acq-famp{shim_mode}_TB1TFL.json\", \"r\") as f:\n",
" metadata = json.load(f)\n",
" ref_voltage = metadata.get(\"TxRefAmp\", \"N/A\")\n",
" print(f\"ref_voltage [V]: {ref_voltage}\")\n",
"\n",
"# Open flip angle map with nibabel\n",
"nii = nib.load(f\"{subject}_acq-famp{shim_mode}_TB1TFL.nii.gz\")\n",
"acquired_fa = nii.get_fdata()\n",
"\n",
"# Siemens maps are in units of flip angle * 10 (in degrees)\n",
"acquired_fa = acquired_fa / 10\n",
"\n",
"# Account for the power loss between the amplifier and the socket. That number was given by Siemens.\n",
"voltage_at_socket = ref_voltage * 10 ** -0.095\n",
"\n",
"# Compute B1 map in [T/V]\n",
"b1_map = (acquired_fa / requested_fa) * (np.pi / (GAMMA * 1e-3 * voltage_at_socket))\n",
"\n",
"# Convert to [nT/V]\n",
"b1_map = b1_map * 1e9\n",
"\n",
"# Save as NIfTI file\n",
"nii_b1 = nib.Nifti1Image(b1_map, nii.affine, nii.header)\n",
"nib.save(nii_b1, f\"{subject}_acq-{shim_mode}_B1plusmap.nii.gz\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a835cdee",
"metadata": {},
"outputs": [],
"source": [
"# Convert the B1 map to nT/V units (by Kyle Gilbert)\n",
"# Convert the flip angle maps to B1+ efficiency maps [nT/V] (inspired by code from Kyle Gilbert)\n",
"# The approach consists in calculating the B1+ efficiency using a 1ms, pi-pulse at the acquisition voltage,\n",
"# then scale the efficiency by the ratio of the measured flip angle to the requested flip angle in the pulse sequence.\n",
"\n",
"# GAMMA = 2.675e8\n",
"# B1eff_mag = (AcquiredFA ./ RequestedFA) .* (pi ./ (GAMMA .* 1e-3 .* VoltageAtSocket)); % [T/V]\n",
"# B1eff_mag = B1eff_mag .* 1e9; % [T/V] to [nT/V]\n",
"# The constants sum up to 130.492, so to convert the B1map to nT/V, it has to be divided by 10 (to get it back into units of FA)\n",
"# then multiplied by 130.492 and divided by the VoltageAtSocket\n",
"GAMMA = 2.675e8; # [rad / (s T)]\n",
"requested_fa = 90 # saturation flip angle -- hard-coded in sequence\n",
"\n",
"for subject in subjects:\n",
" os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
Expand All @@ -276,19 +321,27 @@
" with open(f\"{subject}_acq-famp{shim_mode}_TB1TFL.json\", \"r\") as f:\n",
" metadata = json.load(f)\n",
" ref_voltage = metadata.get(\"TxRefAmp\", \"N/A\")\n",
" print(f\"ref_voltage: {ref_voltage}\")\n",
" print(f\"ref_voltage [V]: {ref_voltage}\")\n",
"\n",
" # Open flip angle map with nibabel\n",
" nii = nib.load(f\"{subject}_acq-famp{shim_mode}_TB1TFL.nii.gz\")\n",
" acquired_fa = nii.get_fdata()\n",
"\n",
" # Siemens maps are in units of flip angle * 10 (in degrees)\n",
" acquired_fa = acquired_fa / 10\n",
"\n",
" voltage_at_socket = ref_voltage * 10**-0.095 # TODO: justify this 0.095 number https://github.com/shimming-toolbox/rf-shimming-7t/issues/22\n",
" # VoltageAtSocket=np.around(VoltageAtSocket, decimals=2)\n",
" # Account for the power loss between the amplifier and the socket. That number was given by Siemens.\n",
" voltage_at_socket = ref_voltage * 10 ** -0.095\n",
"\n",
" # Divide the flip angle map by a factor 10 to get it in the proper unit (https://github.com/shimming-toolbox/rf-shimming-7t/issues/22)\n",
" !sct_maths -i {subject}_acq-famp{shim_mode}_TB1TFL.nii.gz -div 10 -o {subject}_acq-famp{shim_mode}_TB1TFL_nTpV.nii.gz\n",
" # Compute B1 map in [T/V]\n",
" b1_map = (acquired_fa / requested_fa) * (np.pi / (GAMMA * 1e-3 * voltage_at_socket))\n",
"\n",
" # Divide the flip angle map by the voltage at the socket\n",
" !sct_maths -i {subject}_acq-famp{shim_mode}_TB1TFL_nTpV.nii.gz -div {voltage_at_socket} -o {subject}_acq-famp{shim_mode}_TB1TFL_nTpV.nii.gz\n",
" # Convert to [nT/V]\n",
" b1_map = b1_map * 1e9\n",
"\n",
" # Multiply the output by 130.492, which corresponds to ...? (https://github.com/shimming-toolbox/rf-shimming-7t/issues/22)\n",
" !sct_maths -i {subject}_acq-famp{shim_mode}_TB1TFL_nTpV.nii.gz -mul 130.492 -o {subject}_acq-famp{shim_mode}_TB1TFL_nTpV.nii.gz"
" # Save as NIfTI file\n",
" nii_b1 = nib.Nifti1Image(b1_map, nii.affine, nii.header)\n",
" nib.save(nii_b1, f\"{subject}_acq-{shim_mode}_B1plusmap.nii.gz\")"
]
},
{
Expand All @@ -303,7 +356,7 @@
"for subject in subjects:\n",
" os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
" for shim_mode in shim_modes:\n",
" !sct_extract_metric -i {subject}_acq-famp{shim_mode}_TB1TFL_nTpV.nii.gz -f {subject}_acq-anat{shim_mode}_TB1TFL_seg.nii.gz -method wa -vert 1:9 -vertfile {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz -perslice 1 -o TB1TFL_{shim_mode}.csv"
" !sct_extract_metric -i {subject}_acq-{shim_mode}_b1plusmap.nii.gz -f {subject}_acq-anat{shim_mode}_TB1TFL_seg.nii.gz -method wa -vert 1:9 -vertfile {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz -perslice 1 -o TB1TFL_{shim_mode}.csv"
]
},
{
Expand Down

0 comments on commit 9f50c12

Please sign in to comment.