From 6b788521432b988d8718e4d56fde2f5eea4e4641 Mon Sep 17 00:00:00 2001
From: jcohenadad <jcohen@polymtl.ca>
Date: Sat, 23 Dec 2023 14:39:41 -0500
Subject: [PATCH] Added nT/V conversion for TFL B1 maps

---
 data_processing.ipynb | 36 ++++++++++++++++++++++++++++++++----
 1 file changed, 32 insertions(+), 4 deletions(-)

diff --git a/data_processing.ipynb b/data_processing.ipynb
index af43a05..c597272 100644
--- a/data_processing.ipynb
+++ b/data_processing.ipynb
@@ -9,7 +9,7 @@
     "\n",
     "## Data\n",
     "\n",
-    "The data can be download at https://openneuro.org/datasets/ds004906\n",
+    "The data can be downloaded at https://openneuro.org/datasets/ds004906\n",
     "\n",
     "The structure of the input dataset is as follows (JSON sidecars are not listed for clarity):\n",
     "~~~\n",
@@ -64,6 +64,9 @@
     "- Convert the B1 map to nT/V units\n",
     "- Extract the B1 map value within the spinal cord\n",
     "\n",
+    "```{note}\n",
+    "Here is a note!\n",
+    "```\n",
     "Slow processes are indicated with the emoji ⏳"
    ]
   },
@@ -258,9 +261,34 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Convert the B1 map to nT/V units\n",
+    "# Convert the B1 map to nT/V units (by Kyle Gilbert)\n",
     "\n",
-    "# TODO"
+    "# 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",
+    "\n",
+    "for subject in subjects:\n",
+    "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
+    "    for shim_mode in shim_modes:\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: {ref_voltage}\")\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",
+    "\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",
+    "\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",
+    "\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"
    ]
   },
   {
@@ -275,7 +303,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.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 -append 1 -perslice 1 -o TB1TFL_{shim_mode}.csv"
+    "        !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"
    ]
   },
   {