Average B1+ curves across subjects, and scale along x-axis (#30)

* Cleanup imports * Average across subjects * One subplot per subject * Legend only on the top plot * Normalize x-tick to PAM50 vertebrae * Added smoothing * Fixed ylabel * Updated
shimming-toolbox · Jan 9, 2024 · cd5d649 · cd5d649
1 parent c5d4578
commit cd5d649
Show file tree

Hide file tree

Showing 2 changed files with 85 additions and 82 deletions.
diff --git a/.gitignore b/.gitignore
@@ -1,3 +1,9 @@
+# Dataset
+ds004906
+
+# Visual Studio
+.vscode
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]

diff --git a/data_processing.ipynb b/data_processing.ipynb
@@ -84,13 +84,12 @@
     "import json\n",
     "import subprocess\n",
     "import glob\n",
+    "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
-    "# import matplotlib.pyplot as plt\n",
-    "# from PIL import Image\n",
-    "# from IPython.display import display\n",
-    "# from tabulate import tabulate\n",
     "import nibabel as nib\n",
-    "# import pandas as pd"
+    "import pandas as pd\n",
+    "from scipy.interpolate import interp1d\n",
+    "from scipy.ndimage import uniform_filter1d"
    ]
   },
   {
@@ -104,8 +103,8 @@
     "\n",
     "!datalad install https://github.com/OpenNeuroDatasets/ds004906.git\n",
     "os.chdir(\"ds004906\")\n",
-    "# !datalad get . # uncomment for production\n",
-    "!datalad get sub-01/  # debugging\n",
+    "!datalad get . # uncomment for production\n",
+    "# !datalad get sub-01/  # debugging\n",
     "# Get derivatives containing manual labels\n",
     "!datalad get derivatives"
    ]
@@ -205,7 +204,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Extract the signal intensity on the GRE scan within the spinal cord\n",
+    "# Extract the signal intensity on the GRE scan within the spinal cord between levels C1 and T2 (included)\n",
     "\n",
     "for subject in subjects:\n",
     "    # TODO: loop across other shim methods    \n",
@@ -220,7 +219,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Register TFL B1maps to the GRE scan ⏳\n",
+    "# Register TFL flip angle maps to the GRE scan ⏳\n",
     "\n",
     "for subject in subjects:\n",
     "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
@@ -245,7 +244,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Warping spinal cord segmentation and vertebral level to each B1 map\n",
+    "# Warping spinal cord segmentation and vertebral level to each flip angle map\n",
     "\n",
     "for subject in subjects:\n",
     "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
@@ -254,52 +253,6 @@
     "        !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": "97f53bfb",
-   "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,
@@ -321,7 +274,7 @@
     "        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",
+    "            print(f\"ref_voltage [V]: {ref_voltage} ({subject}_acq-famp{shim_mode}_TB1TFL)\")\n",
     "\n",
     "        # Open flip angle map with nibabel\n",
     "        nii = nib.load(f\"{subject}_acq-famp{shim_mode}_TB1TFL.nii.gz\")\n",
@@ -351,7 +304,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Extract B1+ value along the spinal cord\n",
+    "# Extract B1+ value along the spinal cord between levels C1 and T2 (included)\n",
     "\n",
     "for subject in subjects:\n",
     "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
@@ -368,38 +321,82 @@
    "source": [
     "# Make figure of B1+ values along the spinal cord across shim methods\n",
     "\n",
-    "import pandas as pd\n",
-    "import matplotlib.pyplot as plt\n",
-    "import glob\n",
+    "def smooth_data(data, window_size=20):\n",
+    "    \"\"\" Apply a simple moving average to smooth the data. \"\"\"\n",
+    "    return uniform_filter1d(data, size=window_size, mode='nearest')\n",
+    "\n",
+    "# Fixed grid for x-axis\n",
+    "x_grid = np.linspace(0, 1, num_points)\n",
     "\n",
-    "subject='sub-01'\n",
-    "os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
+    "# z-slices corresponding to levels C1 to T2 on the PAM50 template. These will be used to scale the x-label of each subject.\n",
+    "original_vector = np.array([984, 938, 907, 870, 833, 800, 769, 735, 692, 646])\n",
     "\n",
-    "# Pattern to match all files starting with 'TB1TFL_' and ending with '.csv'\n",
-    "file_pattern = f\"TB1map_*.csv\"\n",
+    "# Normalize the PAM50 z-slice numbers to the 1-0 range (to show inferior-superior instead of superior-inferior)\n",
+    "min_val = original_vector.min()\n",
+    "max_val = original_vector.max()\n",
+    "normalized_vector = 1 - ((original_vector - min_val) / (max_val - min_val))\n",
     "\n",
-    "# Initialize a matplotlib figure\n",
-    "plt.figure(figsize=(10, 6))\n",
+    "# Use this normalized vector as x-ticks\n",
+    "custom_xticks = normalized_vector\n",
     "\n",
-    "# Iterate over all files that match the pattern\n",
-    "for file_path in glob.glob(file_pattern):\n",
-    "    # Read the CSV file\n",
-    "    df = pd.read_csv(file_path)\n",
+    "# Number of subjects determines the number of rows in the subplot\n",
+    "n_rows = len(subjects)\n",
     "\n",
-    "    # Assuming 'WA()' is the column name, extract data from this column\n",
-    "    wa_data = df['WA()']\n",
+    "# Create a figure with multiple subplots\n",
+    "fig, axes = plt.subplots(n_rows, 1, figsize=(10, 6 * n_rows))\n",
+    "\n",
+    "# Check if axes is an array or a single object\n",
+    "if n_rows == 1:\n",
+    "    axes = [axes]\n",
+    "\n",
+    "# Iterate over each subject and create a subplot\n",
+    "for i, subject in enumerate(subjects):\n",
+    "    ax = axes[i]\n",
     "    \n",
-    "    # Plot the data\n",
-    "    plt.plot(wa_data, label=file_path.split('/')[-1])  # Using file name as label\n",
+    "    os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
+    "\n",
+    "    for shim_method in shim_modes:\n",
+    "        # Initialize list to collect data for this shim method\n",
+    "        method_data = []\n",
+    "\n",
+    "        for file_path in glob.glob(f\"TB1map_*{shim_method}*.csv\"):\n",
+    "            df = pd.read_csv(file_path)\n",
+    "            wa_data = df['WA()']\n",
+    "\n",
+    "            # Normalize the x-axis to a 0-1 scale for each subject\n",
+    "            x_subject = np.linspace(0, 1, len(wa_data))\n",
+    "\n",
+    "            # Interpolate to the fixed grid\n",
+    "            interp_func = interp1d(x_subject, wa_data, kind='linear', bounds_error=False, fill_value='extrapolate')\n",
+    "            resampled_data = interp_func(x_grid)\n",
+    "\n",
+    "            # Apply smoothing\n",
+    "            smoothed_data = smooth_data(resampled_data)\n",
+    "\n",
+    "            method_data.append(smoothed_data)\n",
+    "\n",
+    "        # If there's data for this shim method, plot it\n",
+    "        if method_data:\n",
+    "            # Plotting each file's data separately\n",
+    "            for resampled_data in method_data:\n",
+    "                ax.plot(x_grid, resampled_data, label=f\"{shim_method} - {file_path.split('/')[-1]}\")\n",
+    "\n",
+    "    # Set custom x-ticks\n",
+    "    ax.set_xticks(custom_xticks)\n",
+    "    ax.set_xticklabels([''] * len(original_vector))\n",
+    "\n",
+    "    ax.set_title(f'{subject}')\n",
+    "#     ax.set_xlabel('Normalized Index')\n",
+    "    ax.set_ylabel('B1+ efficiency [nT/V]')\n",
+    "\n",
+    "    # Add legend only to the first subplot\n",
+    "    if i == 0:\n",
+    "        ax.legend()\n",
     "\n",
-    "# Adding title, labels, and legend\n",
-    "plt.title('Data from Column WA() Across Multiple Files')\n",
-    "plt.xlabel('Index')\n",
-    "plt.ylabel('Values')\n",
-    "plt.legend()\n",
-    "plt.grid(True)\n",
+    "    ax.grid(True)\n",
     "\n",
-    "# Display the plot\n",
+    "# Adjust the layout so labels and titles do not overlap\n",
+    "plt.tight_layout()\n",
     "plt.show()"
    ]
   }