SVS_processing.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/neuropoly/single_voxel_mrs_b0_shimming/blob/main/SVS_processing.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_hjCE6GLlTI7"
      },
      "source": [
        "# SVS processing - interactive notebook"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PnaT-VdNlTI-"
      },
      "source": [
        "This notebook is designed to reproduce the processing pipeline and visualize the results for the paper titled \"On the Impact of B0 Shimming Algorithms on Single Voxel Magnetic Resonance Spectroscopy,\" authored by Behrouz Vejdani Afkham and Eva Alonso-Ortiz.\n",
        "\n",
        "The data for this paper is available online on OSF and will be downloaded in the notebook to reproduce the results.\n",
        "\n",
        "**Note**: While installing new packages, the user may be prompted to restart the session to use them. Please restart the session before proceeding to run the next cells."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pfU21cYmlTI_"
      },
      "source": [
        "## Contents:\n",
        "\n",
        "1. [Install dependencies](#1.-Install-dependencies)\n",
        "2. [Import required libraries](#2.-Import-required-libraries)\n",
        "3. [Download the data](#3.-Download-the-data)\n",
        "4. [B0 Fieldmaps calculation](#4.-Compute-fieldmaps)\n",
        "5. [B0 shimming](#5.-Perform-B0-shimming-using-the-computed-filedmaps)\n",
        "6. [B0 Fieldmaps comparison](#6.-B0-Fieldmaps-comparison)\n",
        "7. [MRS analysis](#7.-MRS-analysis)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Phs6gCJ1lTI_"
      },
      "source": [
        "## 1. Install dependencies"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "gy_EGEgylTI_"
      },
      "outputs": [],
      "source": [
        "%cd /content\n",
        "!git clone --recurse-submodules https://git.fmrib.ox.ac.uk/fsl/fsl_mrs.git\n",
        "%cd /content/fsl_mrs\n",
        "!pip install ."
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "%cd /content\n",
        "!git clone --single-branch --branch gradient_optimizer https://github.com/shimming-toolbox/shimming-toolbox.git\n",
        "%cd /content/shimming-toolbox/\n",
        "!pip install ."
      ],
      "metadata": {
        "id": "V3IxN-qd-NAT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "O0IDwpEwlTJA"
      },
      "source": [
        "## 2. Import required libraries"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "LzNltv9jwoXz"
      },
      "outputs": [],
      "source": [
        "%cd /content\n",
        "import os\n",
        "os.environ[\"LD_PRELOAD\"] = \"\";\n",
        "os.environ[\"APPTAINER_BINDPATH\"] = \"/content\"\n",
        "os.environ[\"MPLCONFIGDIR\"] = \"/content/matplotlib-mpldir\"\n",
        "os.environ[\"LMOD_CMD\"] = \"/usr/share/lmod/lmod/libexec/lmod\"\n",
        "os.environ[\"MODULEPATH\"] = \"/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/molecular_biology:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/workflows:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/visualization:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/structural_imaging:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/statistics:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/spine:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/spectroscopy:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/shape_analysis:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/segmentation:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/rodent_imaging:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/quantitative_imaging:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/quality_control:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/programming:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/phase_processing:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/machine_learning:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/image_segmentation:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/image_registration:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/image_reconstruction:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/hippocampus:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/functional_imaging:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/electrophysiology:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/diffusion_imaging:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/data_organisation:/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/body\"\n",
        "\n",
        "!curl -J -O https://raw.githubusercontent.com/NeuroDesk/neurocommand/main/googlecolab_setup.sh\n",
        "!chmod +x googlecolab_setup.sh\n",
        "!./googlecolab_setup.sh\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "fArgkjVwnurU"
      },
      "outputs": [],
      "source": [
        "import lmod\n",
        "await lmod.load('fsl/6.0.4')\n",
        "await lmod.list()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "q08UP-P5lTJB"
      },
      "outputs": [],
      "source": [
        "import zipfile\n",
        "import json\n",
        "import os\n",
        "import pathlib\n",
        "import shutil\n",
        "import datetime\n",
        "import numpy as np\n",
        "import nibabel as nib\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy import stats\n",
        "import matplotlib.gridspec as gridspec\n",
        "from shimmingtoolbox.prepare_fieldmap import prepare_fieldmap\n",
        "from shimmingtoolbox.load_nifti import read_nii\n",
        "from shimmingtoolbox.cli import b0shim\n",
        "\n",
        "try:\n",
        "    import seaborn as sns\n",
        "\n",
        "except ImportError:\n",
        "    !pip install seaborn\n",
        "    import seaborn as sns\n",
        "\n",
        "!pip install statsmodels\n",
        "!pip install osfclient\n",
        "print('Necessary libraries are imported')\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NEuAWZjBlTJB"
      },
      "source": [
        "## 3. Download the data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "fWe7aUwCt4G3"
      },
      "outputs": [],
      "source": [
        "!osf -p s3gwv fetch /SVS_MRS_DATA.zip"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "id": "2VUS0SmWlTJC"
      },
      "outputs": [],
      "source": [
        "# Unzip the data\n",
        "zipped_data = \"/content/\" + \"SVS_MRS_DATA.zip\"\n",
        "with zipfile.ZipFile(zipped_data,\"r\") as zipfile:\n",
        "  zipfile.extractall('/content/')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7ovVOFv1ve2J"
      },
      "source": [
        "## 4. Compute fieldmaps"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "2U5h7uEJlTJD"
      },
      "outputs": [],
      "source": [
        "if os.path.basename(os.getcwd()) != 'SVS_MRS_DATA':\n",
        "    os.chdir(os.getcwd()+'/SVS_MRS_DATA/')\n",
        "\n",
        "current_path = os.getcwd()\n",
        "print(current_path)\n",
        "shim_methods=['siemens','pi' ,'quadprog' ,'lsq','grad']\n",
        "subjects_list = [name for name in os.listdir() if os.path.isdir(name) and name.startswith('sub')]\n",
        "subjects_list.sort(key=lambda x: int(x.split('-')[1]))\n",
        "print(subjects_list)\n",
        "\n",
        "# Iterate over each subject\n",
        "for subject in subjects_list:\n",
        "    derivative_path = current_path + f\"/derivatives/{subject}/\"\n",
        "    fmap_path = f\"{derivative_path}fmap/\"\n",
        "    test_path = f\"{fmap_path}test/\"\n",
        "    if not os.path.exists(fmap_path):\n",
        "        os.makedirs(fmap_path)\n",
        "    if not os.path.exists(test_path):\n",
        "        os.makedirs(test_path)\n",
        "\n",
        "    mask_mrs_nii = nib.load(derivative_path + 'mask_mrs.nii.gz')\n",
        "    mask_mrs = mask_mrs_nii.get_fdata()\n",
        "    mask_mrs[mask_mrs!=1]= np.nan\n",
        "    mask_anat_nii = nib.load (derivative_path + 'bet_anat_mask.nii.gz')\n",
        "\n",
        "    # Iterate over each shim method\n",
        "    for method in shim_methods:\n",
        "        gre_path =  current_path + f\"/{subject}/fmap/\"\n",
        "\n",
        "        if method==\"siemens\":\n",
        "\n",
        "            mag = read_nii(gre_path + f\"{subject}_acq-siemens-after-Shim_magnitude1.nii.gz\")\n",
        "            phase1 = read_nii(gre_path + f\"{subject}_acq-siemens-after-Shim_phase1.nii.gz\", auto_scale=True)\n",
        "            phase2 = read_nii(gre_path + f\"{subject}_acq-siemens-after-Shim_phase2.nii.gz\", auto_scale=True)\n",
        "\n",
        "            # Load JSON data from the magnitude nifti\n",
        "            with open(gre_path + f\"{subject}_acq-siemens-after-Shim_magnitude1.json\", 'r') as mag_json:\n",
        "                json_data = json.load(mag_json)\n",
        "            phase_list_data = [phase1[2], phase2[2]]\n",
        "\n",
        "            # scale the phase data to (-pi, pi) range\n",
        "            scaled_phase_list_data = [((phase_list_data[i] - np.min(phase_list_data[i]))/(np.max(phase_list_data[i])-np.min(phase_list_data[i])))\n",
        "                                 *2*np.pi-np.pi for i in range(len(phase_list_data))]\n",
        "            # convert to nifti\n",
        "            scaled_phase_list_nii = [nib.Nifti1Image(scaled_phase_list_data[i], phase1[0].affine, phase1[0].header)\n",
        "                                 for i in range(len(scaled_phase_list_data))]\n",
        "\n",
        "            EchoTime1 = phase1[1]['EchoTime']\n",
        "            EchoTime2 = phase2[1]['EchoTime']\n",
        "            echo_times = [EchoTime1, EchoTime2]\n",
        "\n",
        "            # compute fieldmap for each method\n",
        "            fieldmap = prepare_fieldmap(scaled_phase_list_nii, echo_times, mag[2], unwrapper='prelude', gaussian_filter= True, nii_mask= mask_anat_nii, sigma=0.5)\n",
        "\n",
        "            # store the fieldmap under the corresponding path\n",
        "            fieldmap_nii = nib.Nifti1Image(fieldmap[0]*mask_mrs, phase1[0].affine, phase1[0].header )\n",
        "            nib.save(fieldmap_nii, fmap_path + f\"{subject}_acq-siemens-after-Shim-fieldmap.nii.gz\")\n",
        "\n",
        "            # store the corresponding json file for each fieldmap\n",
        "            with open(fmap_path + f\"{subject}_acq-siemens-after-Shim-fieldmap.json\", 'w') as fmap_json:\n",
        "                json.dump(json_data, fmap_json)\n",
        "\n",
        "        else:\n",
        "            mag_before= read_nii(gre_path + f\"{subject}_acq-{method}-before-Shim_magnitude1.nii.gz\")\n",
        "            phase1_before = read_nii(gre_path + f\"{subject}_acq-{method}-before-Shim_phase1.nii.gz\", auto_scale=True)\n",
        "            phase2_before = read_nii(gre_path + f\"{subject}_acq-{method}-before-Shim_phase2.nii.gz\", auto_scale=True)\n",
        "\n",
        "            mag_after = read_nii(gre_path + f\"{subject}_acq-{method}-after-Shim_magnitude1.nii.gz\")\n",
        "            phase1_after = read_nii(gre_path + f\"{subject}_acq-{method}-after-Shim_phase1.nii.gz\", auto_scale=True)\n",
        "            phase2_after = read_nii(gre_path + f\"{subject}_acq-{method}-after-Shim_phase2.nii.gz\", auto_scale=True)\n",
        "\n",
        "            # Load JSON data from the magnitude nifti\n",
        "            with open(gre_path + f\"{subject}_acq-{method}-before-Shim_magnitude1.json\", 'r') as mag_json_before:\n",
        "                json_data_before = json.load(mag_json_before)\n",
        "            # Load JSON data from the magnitude nifti\n",
        "            with open(gre_path + f\"{subject}_acq-{method}-after-Shim_magnitude1.json\", 'r') as mag_json_after:\n",
        "                json_data_after = json.load(mag_json_after)\n",
        "\n",
        "            phase_list_data_before = [phase1_before[2], phase2_before[2]]\n",
        "            phase_list_data_after = [phase1_after[2], phase2_after[2]]\n",
        "\n",
        "            # scale the phase data to (-pi, pi) range\n",
        "            scaled_phase_list_data_before = [((phase_list_data_before[i] - np.min(phase_list_data_before[i]))/(np.max(phase_list_data_before[i])-np.min(phase_list_data_before[i])))\n",
        "                                 *2*np.pi-np.pi for i in range(len(phase_list_data_before))]\n",
        "\n",
        "            scaled_phase_list_data_after = [((phase_list_data_after[i] - np.min(phase_list_data_after[i]))/(np.max(phase_list_data_after[i])-np.min(phase_list_data_after[i])))\n",
        "                     *2*np.pi-np.pi for i in range(len(phase_list_data_after))]\n",
        "\n",
        "            # convert to nifti\n",
        "            scaled_phase_list_nii_before = [nib.Nifti1Image(scaled_phase_list_data_before[i], phase1_before[0].affine, phase1_before[0].header)\n",
        "                                 for i in range(len(scaled_phase_list_data_before))]\n",
        "\n",
        "            scaled_phase_list_nii_after = [nib.Nifti1Image(scaled_phase_list_data_after[i], phase1_after[0].affine, phase1_after[0].header)\n",
        "                                 for i in range(len(scaled_phase_list_data_after))]\n",
        "\n",
        "            EchoTime1 = phase1[1]['EchoTime']\n",
        "            EchoTime2 = phase2[1]['EchoTime']\n",
        "            echo_times = [EchoTime1, EchoTime2]\n",
        "\n",
        "            # compute fieldmap for each method\n",
        "            fieldmap_before = prepare_fieldmap(scaled_phase_list_nii_before, echo_times, mag_before[2], unwrapper='prelude', gaussian_filter= True, nii_mask= mask_anat_nii, sigma=0.5)\n",
        "            fieldmap_after = prepare_fieldmap(scaled_phase_list_nii_after, echo_times, mag_after[2], unwrapper='prelude', gaussian_filter= True, nii_mask= mask_anat_nii, sigma=0.5)\n",
        "\n",
        "            # store the fieldmap under the corresponding path\n",
        "            fieldmap_nii_after = nib.Nifti1Image(fieldmap_after[0]*mask_mrs, phase1_after[0].affine, phase1_after[0].header )\n",
        "            nib.save(fieldmap_nii_after, fmap_path + f\"{subject}_acq-{method}-after-Shim-fieldmap.nii.gz\")\n",
        "\n",
        "            fieldmap_nii_before = nib.Nifti1Image(fieldmap_before[0], phase1_before[0].affine, phase1_before[0].header )\n",
        "            nib.save(fieldmap_nii_before, fmap_path + f\"{subject}_acq-{method}-before-Shim-fieldmap.nii.gz\")\n",
        "\n",
        "            # store the corresponding json file for each fieldmap\n",
        "            with open(fmap_path + f\"{subject}_acq-{method}-before-Shim-fieldmap.json\", 'w') as fmap_json_before:\n",
        "                json.dump(json_data_before, fmap_json_before)\n",
        "\n",
        "            with open(fmap_path + f\"{subject}_acq-{method}-after-Shim-fieldmap.json\", 'w') as fmap_json_after:\n",
        "                json.dump(json_data_after, fmap_json_after)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CEYGsih_ZXKs"
      },
      "outputs": [],
      "source": [
        "### Display field maps for Siemens and Grad methods for a representative subject.\n",
        "\n",
        "fig = plt.figure(figsize=(10, 3), facecolor='white')\n",
        "\n",
        "mask_mrs_sub3 = nib.load(current_path + '/derivatives/sub-3/mask_mrs.nii.gz').get_fdata()\n",
        "mask_mrs_sub3[mask_mrs_sub3 != 1] = np.nan\n",
        "fieldmap_Grad_sub3 = nib.load(current_path+\"/derivatives/sub-3/fmap/sub-3_acq-grad-after-Shim-fieldmap.nii.gz\").get_fdata()*mask_mrs_sub3\n",
        "fieldmap_Siemens_sub3 = nib.load(current_path+\"/derivatives/sub-3/fmap/sub-3_acq-siemens-after-Shim-fieldmap.nii.gz\").get_fdata()*mask_mrs_sub3\n",
        "magnitude = nib.load(current_path+\"/sub-3/fmap/sub-3_acq-grad-before-Shim_magnitude1.nii.gz\").get_fdata()\n",
        "\n",
        "plt.suptitle(' Comparison of the measured fieldmaps [Hz] after shimming', fontsize=14, fontweight='bold')\n",
        "# Define gridspec\n",
        "gs = gridspec.GridSpec(1, 2, width_ratios=[1, 1])\n",
        "\n",
        "# Plot a sagittal cut for Grad method\n",
        "ax1 = plt.subplot(gs[0])\n",
        "ax1.imshow(np.rot90(magnitude[31,:,:]), cmap='gray')\n",
        "img1= ax1.imshow(np.rot90(fieldmap_Grad_sub3[31,:,:]), vmin=-20, vmax=20, cmap='jet')\n",
        "\n",
        "ax1.set_xticks([])  # Turn off x-axis ticks\n",
        "ax1.set_yticks([])  # Turn off y-axis ticks\n",
        "ax1.set_title('Grad', fontsize=12)\n",
        "cbar = plt.colorbar(img1, ax= ax1, fraction=0.025, pad=0.03)\n",
        "cbar.ax.tick_params(labelsize=12)  # Set font size for color bar\n",
        "\n",
        "# Plot a sagittal cut for the Siemens method\n",
        "ax2 = plt.subplot(gs[1])\n",
        "ax2.imshow(np.rot90(magnitude[31,:,:]), cmap='gray')\n",
        "img2 = ax2.imshow(np.rot90(fieldmap_Siemens_sub3[31,:,:]), vmin=-20, vmax=20, cmap='jet')\n",
        "\n",
        "ax2.set_xticks([])  # Turn off x-axis ticks\n",
        "ax2.set_yticks([])  # Turn off y-axis ticks\n",
        "ax2.set_title('Siemens', fontsize=12)\n",
        "\n",
        "# cbar = plt.colorbar(img, ax=ax2)  # Add color bar\n",
        "cbar = plt.colorbar(img2, ax=ax2, fraction=0.025, pad=0.03)\n",
        "cbar.ax.tick_params(labelsize=12)  # Set font size for color bar\n",
        "\n",
        "# Save the figure with 300 dpi\n",
        "plt.savefig('fieldmaps_after_shimming.png', dpi=300)\n",
        "plt.show()\n",
        "print(magnitude.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KyXMxQD_kuYY"
      },
      "source": [
        "## 5. Perform B0 shimming using the computed filedmaps and MRS voxel mask"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "7WSOyK-yZXKu"
      },
      "outputs": [],
      "source": [
        "%%bash\n",
        "# Generate subjects list as an array\n",
        "subjects_list=$(find .  -type d -maxdepth 1 -name \"sub-*\"| sort -V)\n",
        "IFS=$'\\n' read -d '' -r -a subjects_array <<< \"$subjects_list\"\n",
        "echo \"subjects list are: $subjects_list\"\n",
        "\n",
        "path_derivatives=$(find . -type d -name derivatives)\n",
        "shim_methods=(\"pi\" \"quadprog\" \"lsq\" \"grad\")\n",
        "\n",
        "for method in \"${shim_methods[@]}\"; do\n",
        "\n",
        "    for subject in \"${subjects_array[@]}\"; do\n",
        "        echo \"Subject: $subject\"\n",
        "        echo \"Method: $method\"\n",
        "        subject_name=$(basename \"$subject\")\n",
        "        fmap_path=${path_derivatives}/${subject_name}\n",
        "        mask=$(find ${fmap_path} -type f -name \"mask_mrs.nii.gz\")\n",
        "        output_grad=\"${fmap_path}/fmap/static_shim_grad/MRS_mask\"\n",
        "        output_lsq=\"${fmap_path}/fmap/static_shim_lsq/MRS_mask\"\n",
        "        output_pi=\"${fmap_path}/fmap/static_shim_pi/MRS_mask\"\n",
        "        output_quadprog=\"${fmap_path}/fmap/static_shim_quadprog/MRS_mask\"\n",
        "\n",
        "        fieldmap=$(find ${path_derivatives}/${subject_name}/fmap -type f -name \"${subject_name}_acq-${method}-before-Shim-fieldmap.nii.gz\")\n",
        "        echo $fieldmap\n",
        "\n",
        "\n",
        "        all_subjects_path=\"${path_derivatives}/all_subjects\"\n",
        "        if [ ! -d \"$all_subjects_path\" ]; then\n",
        "          mkdir \"$all_subjects_path\"\n",
        "        fi\n",
        "        # Create a file to store the time taken for each method\n",
        "        time_file=\"${all_subjects_path}/${method}_MRS_mask_time.txt\"\n",
        "\n",
        "        if [ \"${method}\" = \"lsq\" ]; then\n",
        "            { time -p st_b0shim dynamic --fmap \"${fieldmap}\" --anat \"${fieldmap}\" --mask \"${mask}\" --mask-dilation-kernel-size \"5\"  --optimizer-method \"least_squares\" --slices \"volume\" --output-file-format-scanner \"chronological-coil\" --scanner-coil-order \"0,1,2\" --output-value-format \"absolute\" --output \"${output_lsq}\"; } 2>&1 | grep real | cut -d' ' -f2 >> \"${time_file}\" || exit\n",
        "\n",
        "        elif [ \"${method}\" = \"grad\" ]; then\n",
        "            { time -p st_b0shim dynamic --fmap \"${fieldmap}\" --anat \"${fieldmap}\" --mask \"${mask}\" --mask-dilation-kernel-size \"5\"  --optimizer-method \"gradient\" --optimizer-criteria \"ps_huber\" --slices \"volume\" --output-file-format-scanner \"chronological-coil\" --scanner-coil-order \"0,1,2\" --output-value-format \"absolute\" --output \"${output_grad}\"; } 2>&1 | grep real | cut -d' ' -f2 >> \"${time_file}\" || exit\n",
        "\n",
        "        elif [ \"${method}\" = \"pi\" ]; then\n",
        "            { time -p st_b0shim dynamic --fmap \"${fieldmap}\" --anat \"${fieldmap}\" --mask \"${mask}\" --mask-dilation-kernel-size \"5\" --optimizer-method \"pseudo_inverse\" --slices \"volume\" --output-file-format-scanner \"chronological-coil\" --scanner-coil-order \"0,1,2\" --output-value-format \"absolute\" --output \"${output_pi}\"; } 2>&1 | grep real | cut -d' ' -f2 >> \"${time_file}\" || exit\n",
        "\n",
        "        elif [ \"${method}\" = \"quadprog\" ]; then\n",
        "            { time -p st_b0shim dynamic --fmap \"${fieldmap}\" --anat \"${fieldmap}\" --mask \"${mask}\" --mask-dilation-kernel-size \"5\" --optimizer-method \"quad_prog\" --slices \"volume\" --output-file-format-scanner \"chronological-coil\" --scanner-coil-order \"0,1,2\" --output-value-format \"absolute\" --output \"${output_quadprog}\"; } 2>&1 | grep real | cut -d' ' -f2 >> \"${time_file}\" || exit\n",
        "\n",
        "        fi\n",
        "    done\n",
        "done"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Bo8ZMIXxZXKu"
      },
      "outputs": [],
      "source": [
        "### Display the average runtime for each of the shimming methods\n",
        "\n",
        "def calculate_average(file_path):\n",
        "    with open(file_path, 'r') as file:\n",
        "        times = [float(line.strip()) for line in file]\n",
        "        return sum(times) / len(times)\n",
        "\n",
        "path = current_path+\"/derivatives/all_subjects\"\n",
        "shim_methods = [\"pi\", \"quadprog\", \"lsq\", \"grad\"]\n",
        "\n",
        "for method in shim_methods:\n",
        "    file_path = os.path.join(path, f\"{method}_MRS_mask_time.txt\")\n",
        "    if os.path.exists(file_path):\n",
        "        average_time = np.round(calculate_average(file_path), 1)\n",
        "        print(f\"Average computation time for {method} in MRS mask: {average_time} seconds\")\n",
        "    else:\n",
        "        print(f\"File not found for {method}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Q0qCBPnAlTJD"
      },
      "source": [
        "## 6. B0 Fieldmaps comparison"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HoEPEgcgkuYY"
      },
      "source": [
        "### 6.1 Function to analyze the fieldmaps"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "metadata": {
        "id": "VMkUSY6K0T7j"
      },
      "outputs": [],
      "source": [
        "def analyse_fieldmap(source_dir, labels, fieldmaps, output, mask, subject, roi_name, fig_size, colors, y_offset):\n",
        "    \"\"\"\n",
        "    Description: Generates a violin plot representing the distribution of the given fieldmaps\n",
        "    labels: List of strings including fieldmaps label\n",
        "    fieldmaps: List of numpy arrays of fieldmaps in the same order as labels\n",
        "    output: Output folder to save the violin plot\n",
        "    mask: Numpy array of a Mask corresponding to the MRS voxel volume\n",
        "    \"\"\"\n",
        "    plt.figure(figsize=fig_size)\n",
        "    sns.set(font_scale=2)\n",
        "\n",
        "    output_path = os.path.join(source_dir, output)\n",
        "    if not os.path.exists(output_path):\n",
        "        os.makedirs(output_path)\n",
        "\n",
        "    if np.shape(mask)==np.shape(fieldmaps[0]):\n",
        "        dict_map_avg = {label: fieldmap[mask == 1] for label, fieldmap in zip(labels, fieldmaps)}\n",
        "    else:\n",
        "        dict_map_avg = {label: fieldmap for label, fieldmap in zip(labels, fieldmaps)}\n",
        "\n",
        "    def calculate_rmse(data):\n",
        "        mse = np.mean(np.square(data))\n",
        "        rmse = np.sqrt(mse)\n",
        "        return rmse\n",
        "\n",
        "    dict_means = {key: np.round(np.mean(value), 1) for key, value in dict_map_avg.items()}\n",
        "    dict_stds = {key: np.round(np.std(value), 1) for key, value in dict_map_avg.items()}\n",
        "    dict_rmse = {key: calculate_rmse(value) for key, value in dict_map_avg.items()}\n",
        "\n",
        "    ax = sns.violinplot(data=list(dict_map_avg.values()), palette=colors)\n",
        "\n",
        "    # Get the limits of the y-axis to position the text\n",
        "    ymin, ymax = ax.get_ylim()\n",
        "    for i, (key, mean) in enumerate(dict_means.items()):\n",
        "        mean_str = \"{:.1f}\".format(mean)\n",
        "        std_str = \"{:.1f}\".format(dict_stds[key])\n",
        "        rmse_str = \"{:.1f}\".format(dict_rmse[key])\n",
        "        ax.text(i + 0.05, ymin + y_offset * (ymax - ymin), f\"SD: {std_str}\",\n",
        "                horizontalalignment='left', verticalalignment='top', fontsize=14)\n",
        "\n",
        "    ax.set_ylabel('Frequency [Hz]', fontsize=14, fontweight='bold')\n",
        "    ax.set_xticks(range(len(labels)))\n",
        "    ax.set_xticklabels(labels, fontsize=14 , fontweight='bold', rotation=25)\n",
        "    ax.set_title(f\"{subject}: $\\Delta$B0 distribution_within_{roi_name}\", fontsize=18, fontweight='bold')\n",
        "\n",
        "    ax.tick_params(axis='y', labelsize=14)\n",
        "    plt.savefig(os.path.join(output_path, f\"{subject}_B0_distribution_within_{roi_name}_mask.png\"), dpi=300)\n",
        "    plt.show()\n",
        "\n",
        "    return dict_stds\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QWQO-9gCkuYZ"
      },
      "source": [
        "### 6.2 Display the fieldmaps distribution"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "vUpZ33DfZXKw"
      },
      "outputs": [],
      "source": [
        "all_Grad_fmaps_mrs_mask_unshimmed=[]\n",
        "all_PI_fmaps_mrs_mask_unshimmed=[]\n",
        "all_LSq_fmaps_mrs_mask_unshimmed=[]\n",
        "all_QuadProg_fmaps_mrs_mask_unshimmed=[]\n",
        "\n",
        "all_Siemens_fmaps_mrs_mask_shimmed=[]\n",
        "all_Grad_fmaps_mrs_mask_shimmed=[]\n",
        "all_PI_fmaps_mrs_mask_shimmed=[]\n",
        "all_LSq_fmaps_mrs_mask_shimmed=[]\n",
        "all_QuadProg_fmaps_mrs_mask_shimmed=[]\n",
        "sd_field_maps_mrs_mask = {}\n",
        "methods = ['pi',  'quadprog','lsq', 'grad']\n",
        "\n",
        "# Iterate over each subject\n",
        "for subject in subjects_list:\n",
        "    sd_field_maps_mrs_mask[subject] = {}\n",
        "    derivative_path = current_path + f\"/derivatives/{subject}/\"\n",
        "    mask_mrs_nii = nib.load(derivative_path + 'mask_mrs.nii.gz')\n",
        "    mask_mrs = mask_mrs_nii.get_fdata()\n",
        "\n",
        "    unshimmed_fieldmap_grad = nib.load (derivative_path + f\"fmap/{subject}_acq-grad-before-Shim-fieldmap.nii.gz\").get_fdata()\n",
        "    unshimmed_fieldmap_pi = nib.load (derivative_path + f\"fmap/{subject}_acq-pi-before-Shim-fieldmap.nii.gz\").get_fdata()\n",
        "    unshimmed_fieldmap_lsq= nib.load (derivative_path + f\"fmap/{subject}_acq-lsq-before-Shim-fieldmap.nii.gz\").get_fdata()\n",
        "    unshimmed_fieldmap_quadprog = nib.load (derivative_path + f\"fmap/{subject}_acq-quadprog-before-Shim-fieldmap.nii.gz\").get_fdata()\n",
        "\n",
        "    shimmed_fieldmap_siemens = nib.load (derivative_path + f\"fmap/{subject}_acq-siemens-after-Shim-fieldmap.nii.gz\").get_fdata()\n",
        "    shimmed_fieldmap_grad = nib.load (derivative_path + f\"fmap/{subject}_acq-grad-after-Shim-fieldmap.nii.gz\").get_fdata()\n",
        "    shimmed_fieldmap_pi = nib.load (derivative_path + f\"fmap/{subject}_acq-pi-after-Shim-fieldmap.nii.gz\").get_fdata()\n",
        "    shimmed_fieldmap_lsq = nib.load (derivative_path + f\"fmap/{subject}_acq-lsq-after-Shim-fieldmap.nii.gz\").get_fdata()\n",
        "    shimmed_fieldmap_quadprog  = nib.load (derivative_path + f\"fmap/{subject}_acq-quadprog-after-Shim-fieldmap.nii.gz\").get_fdata()\n",
        "\n",
        "    for method in methods:\n",
        "        fieldmap_unshimmed = f\"unshimmed_fieldmap_{method}\"\n",
        "        fieldmap_shimmed = f\"shimmed_fieldmap_{method}\"\n",
        "        fieldmap1 = eval(fieldmap_unshimmed)\n",
        "        fieldmap2 = eval(fieldmap_shimmed)\n",
        "        mean1 = np.mean(fieldmap1[mask_mrs > 0])\n",
        "        mean2 = np.mean(fieldmap2[mask_mrs > 0])\n",
        "        if mean1 > 0:\n",
        "            fieldmap1[mask_mrs > 0] -= mean1\n",
        "        else:\n",
        "            fieldmap1[mask_mrs > 0] += abs(mean1)\n",
        "        if mean2 > 0:\n",
        "            fieldmap2[mask_mrs > 0] -= mean2\n",
        "        else:\n",
        "            fieldmap2[mask_mrs > 0] += abs(mean2)\n",
        "\n",
        "    all_Grad_fmaps_mrs_mask_unshimmed.append(unshimmed_fieldmap_grad[mask_mrs==1].flatten())\n",
        "    all_PI_fmaps_mrs_mask_unshimmed.append(unshimmed_fieldmap_pi[mask_mrs==1].flatten())\n",
        "    all_LSq_fmaps_mrs_mask_unshimmed.append(unshimmed_fieldmap_lsq[mask_mrs==1].flatten())\n",
        "    all_QuadProg_fmaps_mrs_mask_unshimmed.append(unshimmed_fieldmap_quadprog[mask_mrs==1].flatten())\n",
        "\n",
        "    all_Siemens_fmaps_mrs_mask_shimmed.append(shimmed_fieldmap_siemens[mask_mrs==1].flatten())\n",
        "    all_Grad_fmaps_mrs_mask_shimmed.append(shimmed_fieldmap_grad[mask_mrs==1].flatten())\n",
        "    all_PI_fmaps_mrs_mask_shimmed.append(shimmed_fieldmap_pi[mask_mrs==1].flatten())\n",
        "    all_LSq_fmaps_mrs_mask_shimmed.append(shimmed_fieldmap_lsq[mask_mrs==1].flatten())\n",
        "    all_QuadProg_fmaps_mrs_mask_shimmed.append(shimmed_fieldmap_quadprog[mask_mrs==1].flatten())\n",
        "\n",
        "    fmaps = [unshimmed_fieldmap_pi, shimmed_fieldmap_pi,\n",
        "             unshimmed_fieldmap_quadprog, shimmed_fieldmap_quadprog,\n",
        "             unshimmed_fieldmap_lsq, shimmed_fieldmap_lsq,\n",
        "             unshimmed_fieldmap_grad, shimmed_fieldmap_grad,\n",
        "             shimmed_fieldmap_siemens]\n",
        "\n",
        "    labels = [\"Unshimmed PI\", \"Shimmed PI\",\n",
        "              \"Unshimmed QuadProg\", \"Shimmed QuadProg\",\n",
        "              \"Unshimmed LSq\", \"Shimmed LSq\",\n",
        "              \"Unshimmed Grad\", \"Shimmed Grad\",\n",
        "              \"Shimmed Siemens\"]\n",
        "\n",
        "    colors = ['orange', 'orange', 'teal', 'teal',\n",
        "              'pink', 'pink', 'cyan', 'cyan', 'gold']\n",
        "    output = derivative_path\n",
        "    sd = analyse_fieldmap(current_path, labels, fmaps, output, mask_mrs, subject, roi_name='MRS voxel',fig_size=(16,12), colors=colors, y_offset=0.3)\n",
        "    sd_field_maps_mrs_mask[subject] =sd\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "X1atfozwZXKw"
      },
      "outputs": [],
      "source": [
        "# Calculate standard deviation for each method across all subjects\n",
        "std_dict = {}\n",
        "for key in sd_field_maps_mrs_mask['sub-2'].keys():\n",
        "    values = [sd_field_maps_mrs_mask[sub][key] for sub in sd_field_maps_mrs_mask]\n",
        "    std_dict[key] = np.round(np.std(values),1)\n",
        "\n",
        "print(\"Standard Deviation for Each method across all subjects:\")\n",
        "for key, value in std_dict.items():\n",
        "    print(f\"{key}: {value} Hz\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pZ-gGInpZXKw"
      },
      "outputs": [],
      "source": [
        "### Display the fieldmaps distributions for all subjects\n",
        "shim_methods=[ \"PI\", \"QuadProg\", \"LSq\", \"Grad\"]\n",
        "\n",
        "for method in shim_methods:\n",
        "    # Create the variable name dynamically\n",
        "    var_name_unshimmed = f\"all_{method}_data_mrs_mask_unshimmed\"\n",
        "    var_name_shimmed = f\"all_{method}_data_mrs_mask_shimmed\"\n",
        "    var_Siemens= f\"all_Siemens_data_mrs_mask_shimmed\"\n",
        "    # Extract the data for the current method\n",
        "    data_unshimmed = locals()[f\"all_{method}_fmaps_mrs_mask_unshimmed\"]\n",
        "    data_shimmed = locals()[f\"all_{method}_fmaps_mrs_mask_shimmed\"]\n",
        "    data_Siemens = locals()[f\"all_Siemens_fmaps_mrs_mask_shimmed\"]\n",
        "\n",
        "    concat_data_unshimmed = np.concatenate(data_unshimmed)\n",
        "    concat_data_shimmed = np.concatenate(data_shimmed)\n",
        "    concat_data_Siemens = np.concatenate(data_Siemens)\n",
        "\n",
        "    # Assign the concat_data to the dynamically created variable\n",
        "    exec(f\"{var_name_unshimmed} = concat_data_unshimmed\")\n",
        "    exec(f\"{var_name_shimmed} = concat_data_shimmed\")\n",
        "    exec(f\"{var_Siemens} = concat_data_Siemens\")\n",
        "\n",
        "fmaps = [all_PI_data_mrs_mask_unshimmed, all_PI_data_mrs_mask_shimmed, all_QuadProg_data_mrs_mask_unshimmed, all_QuadProg_data_mrs_mask_shimmed,\n",
        "         all_LSq_data_mrs_mask_unshimmed, all_LSq_data_mrs_mask_shimmed, all_Grad_data_mrs_mask_unshimmed, all_Grad_data_mrs_mask_shimmed, all_Siemens_data_mrs_mask_shimmed]\n",
        "labels = [\"Unshimmed PI\", \"Shimmed PI\", \"Unshimmed QuadProg\", \"Shimmed QuadProg\", \"Unshimmed LSq\", \"Shimmed LSq\", \"Unshimmed Grad\", \"Shimmed Grad\", \"Shimmed Siemens\"]\n",
        "\n",
        "source_dir = os.getcwd()\n",
        "output = 'derivatives/all_subjects'\n",
        "subject = 'all subjects'\n",
        "colors = ['orange', 'orange',  'teal', 'teal',\n",
        "          'pink', 'pink',  'cyan', 'cyan', 'gold']\n",
        "distribution = analyse_fieldmap(source_dir, labels, fmaps, output, mask_mrs, subject, roi_name='MRS',fig_size=(16,14), colors=colors, y_offset=0.3)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OWxV_fOYkuYZ"
      },
      "source": [
        "### 6.3 Perform statistical test on fieldmaps"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JHG45IZckuYa"
      },
      "source": [
        "By comparing the variance, the following cell assesses whether the distribution of each shimming method significantly differs from the Siemens shimming.\n",
        "The documentation for the following tests can be found [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kstest.html) and [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.levene.html).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "oEeGCp55ZXK1"
      },
      "outputs": [],
      "source": [
        "import statsmodels.stats.multitest as smm\n",
        "# A Kolmogorov-Smirnov test to check the normality of the reference distribution (Siemens shimming)\n",
        "statistic, p_value_normality = stats.kstest(all_Siemens_data_mrs_mask_shimmed, 'norm')\n",
        "# Set significant level\n",
        "alpha = 0.01\n",
        "\n",
        "# Interpret the results\n",
        "# H0: the sample has a Gaussian distribution.\n",
        "# H1: the sample does not have a Gaussian distribution.\n",
        "if p_value_normality < alpha:\n",
        "    print('Reject the null hypothesis. The provided data does not follow a normal distribution.')\n",
        "else:\n",
        "    print(\"Fail to reject the null hypothesis. The provided data follow a normal distribution.\")\n",
        "\n",
        "comparison_methods = ['PI','Grad','QuadProg','LSq']\n",
        "Uncor_pvals = []\n",
        "\n",
        "# levene test\n",
        "for method in comparison_methods:\n",
        "    data = eval(f\"all_{method}_data_mrs_mask_shimmed\")\n",
        "    if np.var(data)< np.var(all_Siemens_data_mrs_mask_shimmed):\n",
        "      statistic, p_value = stats.levene(all_Siemens_data_mrs_mask_shimmed, data)\n",
        "      Uncor_pvals.append(p_value)\n",
        "    else:\n",
        "      comparison_methods.remove(method)\n",
        "      print(f\"The variance of the {method} is higher than Siemens and excluded from further comparisons.\")\n",
        "\n",
        "print(\"Uncorrected p-values are: \", Uncor_pvals)\n",
        "# Perform Multiple comparison correction with step-down method using Bonferroni adjustments\n",
        "Corrected_pvals = smm.multipletests(Uncor_pvals, alpha=0.01, method='holm', is_sorted=False, returnsorted=False)\n",
        "print(\"Corrected p-values are: \" , Corrected_pvals[1])\n",
        "\n",
        "for pval, method in zip(Corrected_pvals[1], comparison_methods):\n",
        "  if pval < 0.01:\n",
        "    print(f\"Reject the null hypothesis. The variance between the Siemens and {method} methods are significantly different.\")\n",
        "  else:\n",
        "    print(f\"Fail to reject the null hypothesis. There is no significant difference in variances between the Siemens and {method} methods.\")\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "4LXgQCZDZXK2"
      },
      "outputs": [],
      "source": [
        "# Comparing Grad method to other methods\n",
        "\n",
        "comparison_methods= ['PI', 'QuadProg', 'LSq']\n",
        "# levene test\n",
        "for method in comparison_methods:\n",
        "    data = eval(f\"all_{method}_data_mrs_mask_shimmed\")\n",
        "    statistic, p_value = stats.levene(all_Grad_data_mrs_mask_shimmed, data)\n",
        "\n",
        "    # Set significancy level\n",
        "    alpha = 0.01\n",
        "\n",
        "    # Interpret the results\n",
        "    if p_value < alpha:\n",
        "        print(f\"Reject the null hypothesis. The variance between the Grad and {method} groups are significantly different.\")\n",
        "    else:\n",
        "        print(f\"Fail to reject the null hypothesis. There is no significant difference in variances between the Grad and {method} groups.\")\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Zgixo125ZXK2"
      },
      "outputs": [],
      "source": [
        "# Comparing PI method to other methods\n",
        "\n",
        "comparison_methods= ['QuadProg', 'LSq']\n",
        "# levene test\n",
        "for method in comparison_methods:\n",
        "    data = eval(f\"all_{method}_data_mrs_mask_shimmed\")\n",
        "    statistic, p_value = stats.levene(all_PI_data_mrs_mask_shimmed, data)\n",
        "\n",
        "    # Set significancy level\n",
        "    alpha = 0.01\n",
        "\n",
        "    # Interpret the results\n",
        "    if p_value < alpha:\n",
        "        print(f\"Reject the null hypothesis. The variance between the PI and {method} groups are significantly different.\")\n",
        "    else:\n",
        "        print(f\"Fail to reject the null hypothesis. There is no significant difference in variances between the PI and {method} groups.\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DyoZVUyblTJD"
      },
      "source": [
        "## 7. MRS analysis"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "edgU0UwLkuYa"
      },
      "source": [
        "###  7.1. Water removal and fitting"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "uEz2anTYkuYb"
      },
      "outputs": [],
      "source": [
        "%%bash\n",
        "\n",
        "subjects_list=$(find .  -type d -maxdepth 1 -name \"sub-*\")\n",
        "IFS=$'\\n' read -d '' -r -a subjects_array <<< \"$subjects_list\"\n",
        "echo \"subjects list are: $subjects_list\"\n",
        "\n",
        "path_derivatives=$(find . -type d -name derivatives)\n",
        "\n",
        "for subject in \"${subjects_array[@]}\"; do\n",
        "    echo \"Analyzing subject: ${subject}\"\n",
        "    path_mrs=$(find ${subject} -type d -name mrs)\n",
        "    echo \"mrs path is: ${path_mrs}\"\n",
        "    path_t1=$(find ${subject} -type f -name *T1w.nii.gz)\n",
        "\n",
        "    sub_basename=$(basename \"${subject}\")\n",
        "    subject_derivatives=${path_derivatives}/${sub_basename}/mrs\n",
        "    echo \"${sub_basename} derivatives path is: ${subject_derivatives}\"\n",
        "\n",
        "    shim_methods=(\"siemens\" \"pi\" \"quadprog\" \"lsq\" \"grad\")\n",
        "    for method in \"${shim_methods[@]}\"; do\n",
        "        fsl_mrs_proc remove --file ${path_mrs}/${subject}_acq-press-${method}-shim_nuc-H_echo-135_svs.nii.gz --output ${subject_derivatives}/${method}_water_removed -r  --ppm 3.5 9\n",
        "        fsl_mrs  --data ${subject_derivatives}/${method}_water_removed/${subject}_acq-press-${method}-shim_nuc-H_echo-135_svs.nii.gz --basis ${path_derivatives}/LCModel_Siemens_UnEdited_PRESS_135_ALL.BASIS --t1 ${path_t1} --output ${subject_derivatives}/${method}_fit_result --baseline_order 2   --internal_ref 'Cr' --ignore H2O PCr GPC  --report --overwrite --ppmlim 0 4 --TE 135 --TR 2\n",
        "    done\n",
        "done\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UF0j-YDBkuYb"
      },
      "source": [
        "### 7.2 Plot the fittings"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5_da0nkJkuYb"
      },
      "outputs": [],
      "source": [
        "current_path = os.getcwd()\n",
        "index = 1\n",
        "\n",
        "# Create a figure with a 5x5 grid of subplots\n",
        "fig, axes = plt.subplots(5, 5, figsize=(40, 40))\n",
        "\n",
        "# Reduce the gap between subplots\n",
        "plt.subplots_adjust(wspace=0, hspace=0.01)\n",
        "\n",
        "shim_methods = ['siemens', 'pi', 'quadprog', 'lsq', 'grad']\n",
        "\n",
        "for i, subject in enumerate(subjects_list):\n",
        "    for j, method in enumerate(shim_methods):\n",
        "        ax = axes[i, j]\n",
        "        fit_path = os.path.join(current_path, \"derivatives\", subject, \"mrs\", f\"{method}_fit_result\")\n",
        "        ax.imshow(plt.imread(os.path.join(fit_path, 'fit_summary.png')))\n",
        "        ax.set_xticks([])\n",
        "        ax.set_yticks([])\n",
        "\n",
        "        # Add titles for the first row\n",
        "        if i == 0:\n",
        "            ax.set_title(method.capitalize(), fontsize=32,fontweight='bold', pad=20)\n",
        "\n",
        "        # Add labels for the first column\n",
        "        if j == 0:\n",
        "            ax.set_ylabel(subject, fontsize=32, fontweight='bold',labelpad=20)\n",
        "\n",
        "# Save the plot with tight layout\n",
        "plt.savefig(os.path.join(current_path, 'derivatives', 'all_subjects', 'mrs_fittings.png'), dpi=300, bbox_inches='tight')\n",
        "\n",
        "# Show the plot\n",
        "plt.show()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qOsM9FkgkuYb"
      },
      "source": [
        "### 7.3. Extract the metabolites fit results from the CSV files"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "rbgG8qEXkuYb"
      },
      "outputs": [],
      "source": [
        "# path to the current directory\n",
        "current_path = os.getcwd()\n",
        "\n",
        "# Iterate over each subject and shim method to extract fit results\n",
        "crlb_NAA_all = {}\n",
        "fwhm_NAA_all = {}\n",
        "snr_NAA_all = {}\n",
        "\n",
        "index = 1\n",
        "for subject in subjects_list:\n",
        "    crlb_NAA_all[subject] = {}\n",
        "    fwhm_NAA_all[subject] = {}\n",
        "    snr_NAA_all[subject] = {}\n",
        "\n",
        "    for method in shim_methods:\n",
        "        fit_path = os.path.join(current_path, f\"derivatives/{subject}/mrs/{method}_fit_result/\")\n",
        "\n",
        "        # Read summary.csv file\n",
        "        summary_file_path = os.path.join(fit_path, 'summary.csv')\n",
        "        if os.path.exists(summary_file_path):\n",
        "            df = pd.read_csv(summary_file_path)\n",
        "\n",
        "            # Get the fit results for NAA\n",
        "            CRLB_NAA = df.at[11, '%CRLB']\n",
        "            fwhm_NAA = df.at[11, 'FWHM']\n",
        "            snr_NAA = df.at[11, 'SNR']\n",
        "            crlb_NAA_all[subject][method] = CRLB_NAA\n",
        "            fwhm_NAA_all[subject][method] = fwhm_NAA\n",
        "            snr_NAA_all[subject][method] = snr_NAA"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AVF25dC0lTJF"
      },
      "source": [
        "### 7.4. Comparing NAA's FWHM, SNR and CRLB across differnt shim methods"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ewauhiqUkuYc"
      },
      "outputs": [],
      "source": [
        "import warnings\n",
        "warnings.filterwarnings('ignore')\n",
        "import matplotlib.markers as mmarkers\n",
        "\n",
        "# Create figure and axes\n",
        "fig, axs = plt.subplots(1, 3, figsize=(18, 6))\n",
        "\n",
        "\n",
        "# Compute the mean values for FWHM, SNR, and CRLB for NAA\n",
        "mean_fwhm = {key: np.round(np.mean([fwhm_NAA_all[sub][key] for sub in fwhm_NAA_all]), decimals=1) for key in fwhm_NAA_all['sub-1'].keys()}\n",
        "mean_snr = {key: np.round(np.mean([snr_NAA_all[sub][key] for sub in snr_NAA_all]), decimals=1) for key in snr_NAA_all['sub-1'].keys()}\n",
        "mean_crlb = {key: np.round(np.mean([crlb_NAA_all[sub][key] for sub in crlb_NAA_all]), decimals=1) for key in crlb_NAA_all['sub-1'].keys()}\n",
        "\n",
        "# Define a horizontal line marker\n",
        "hline_marker = mmarkers.MarkerStyle('_')\n",
        "hline_marker._transform = hline_marker.get_transform().scale(5, 5)\n",
        "\n",
        "# Plot FWHM\n",
        "axs[0].set_title(\"FWHM (Hz) NAA\", fontsize=18)\n",
        "for sub in fwhm_NAA_all:\n",
        "    axs[0].scatter(list(fwhm_NAA_all[sub].keys()), list(fwhm_NAA_all[sub].values()), label=sub)\n",
        "axs[0].scatter(list(mean_fwhm.keys()), list(mean_fwhm.values()), color='black', marker=hline_marker, label='mean')\n",
        "axs[0].legend(prop={'size': 15}, loc='upper right', bbox_to_anchor=(0.4, 1))\n",
        "axs[0].tick_params(axis='x', rotation=45, labelsize=14)\n",
        "axs[0].yaxis.set_minor_locator(plt.MultipleLocator(1))  # Set minor ticks for y-axis\n",
        "axs[0].grid(which='major', axis='x')\n",
        "\n",
        "# Plot SNR\n",
        "axs[1].set_title(\"SNR (A.U.) NAA\", fontsize=18)\n",
        "for sub in snr_NAA_all:\n",
        "    axs[1].scatter(list(snr_NAA_all[sub].keys()), list(snr_NAA_all[sub].values()), label=sub)\n",
        "axs[1].scatter(list(mean_snr.keys()), list(mean_snr.values()), color='black', marker=hline_marker, label='mean')\n",
        "axs[1].legend(prop={'size': 15}, loc='upper right', bbox_to_anchor=(0.9, 1))\n",
        "axs[1].tick_params(axis='x', rotation=45, labelsize=14)\n",
        "axs[1].yaxis.set_minor_locator(plt.MultipleLocator(1))  # Set minor ticks for y-axis\n",
        "axs[1].grid(which='major', axis='x')\n",
        "\n",
        "# Plot CRLB\n",
        "axs[2].set_title(\"CRLB (%) NAA\", fontsize=18)\n",
        "for sub in crlb_NAA_all:\n",
        "    axs[2].scatter(list(crlb_NAA_all[sub].keys()), list(crlb_NAA_all[sub].values()), label=sub)\n",
        "axs[2].scatter(list(mean_crlb.keys()), list(mean_crlb.values()), color='black', marker=hline_marker, label='mean')\n",
        "axs[2].legend(prop={'size': 15}, loc='upper right', bbox_to_anchor=(0.5, 1))\n",
        "axs[2].tick_params(axis='x', rotation=45, labelsize=14)\n",
        "axs[2].yaxis.set_minor_locator(plt.MultipleLocator(1))  # Set minor ticks for y-axis\n",
        "axs[2].grid(which='major', axis='x')\n",
        "\n",
        "# Set x tick labels\n",
        "axs[0].set_xticklabels([\"Siemens\", \"PI\", \"QuadProg\", \"LSq\", \"Grad\"])\n",
        "axs[1].set_xticklabels([\"Siemens\", \"PI\", \"QuadProg\", \"LSq\", \"Grad\"])\n",
        "axs[2].set_xticklabels([\"Siemens\", \"PI\", \"QuadProg\", \"LSq\", \"Grad\"])\n",
        "\n",
        "# Adjust layout\n",
        "fig.tight_layout()\n",
        "plt.savefig(os.path.join(current_path, 'derivatives', 'all_subjects', 'FWHM_SNR_comparison.png'), dpi=300, bbox_inches='tight')\n",
        "\n",
        "# Show the plot\n",
        "plt.show()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mG21mJzQkuYc"
      },
      "source": [
        "### 7.5. Create a CSV output  "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "jVBs4Dz5kuYc"
      },
      "outputs": [],
      "source": [
        "\n",
        "# Initialize lists for each criterion\n",
        "all_data_dicts_fwhm = []\n",
        "all_data_dicts_snr = []\n",
        "all_data_dicts_crlb = []\n",
        "current_path = os.getcwd()\n",
        "criteria_list = ['fwhm', 'snr', 'crlb']\n",
        "\n",
        "# Create a dictionary for easy access to mean values\n",
        "mean_criteria_dict = {\n",
        "    'fwhm': mean_fwhm,\n",
        "    'snr': mean_snr,\n",
        "    'crlb': mean_crlb\n",
        "}\n",
        "\n",
        "for criteria in criteria_list:\n",
        "    mean_criteria = mean_criteria_dict[criteria]\n",
        "\n",
        "    # Calculate the percentage difference for each dictionary\n",
        "    percentage_diff_dict = {}\n",
        "    siemens_value = mean_criteria['siemens']\n",
        "    percentage_diff = {key: ((value - siemens_value) / siemens_value) * 100 for key, value in mean_criteria.items() if key != 'siemens'}\n",
        "\n",
        "    combined_dict = {'criteria': criteria.upper()}\n",
        "    for key, value in mean_criteria.items():\n",
        "        combined_dict[f'{key}'] = round(value, 1)  # Round to one decimal place\n",
        "        if key != 'siemens':\n",
        "            combined_dict[f'{key} change (%)'] = round(percentage_diff.get(key, None), 1)  # Round to one decimal place\n",
        "\n",
        "    if criteria == 'fwhm':\n",
        "        all_data_dicts_fwhm.append(combined_dict)\n",
        "    elif criteria == 'snr':\n",
        "        all_data_dicts_snr.append(combined_dict)\n",
        "    elif criteria == 'crlb':\n",
        "        all_data_dicts_crlb.append(combined_dict)\n",
        "\n",
        "# Create DataFrames for each criterion\n",
        "df_fwhm = pd.DataFrame(all_data_dicts_fwhm)\n",
        "df_snr = pd.DataFrame(all_data_dicts_snr)\n",
        "df_crlb = pd.DataFrame(all_data_dicts_crlb)\n",
        "\n",
        "# Display the DataFrames\n",
        "df_fwhm_formatted = df_fwhm.apply(lambda x: x.apply(lambda y: '{:.1f}'.format(y) if isinstance(y, float) else y))\n",
        "df_snr_formatted = df_snr.apply(lambda x: x.apply(lambda y: '{:.1f}'.format(y) if isinstance(y, float) else y))\n",
        "df_crlb_formatted = df_crlb.apply(lambda x: x.apply(lambda y: '{:.1f}'.format(y) if isinstance(y, float) else y))\n",
        "\n",
        "# Display the formatted DataFrames\n",
        "print(\"FWHM DataFrame:\")\n",
        "display(df_fwhm_formatted)\n",
        "\n",
        "print(\"SNR DataFrame:\")\n",
        "display(df_snr_formatted)\n",
        "\n",
        "print(\"CRLB DataFrame:\")\n",
        "display(df_crlb_formatted)\n",
        "\n",
        "# Create separator rows\n",
        "separator_row_snr = pd.DataFrame({'criteria': ['SNR']})\n",
        "separator_row_crlb = pd.DataFrame({'criteria': ['CRLB']})\n",
        "separator_row_columns = pd.DataFrame([df_fwhm_formatted.columns], columns=df_fwhm_formatted.columns)\n",
        "\n",
        "# Concatenate the DataFrames with the separator rows\n",
        "combined_df = pd.concat([df_fwhm_formatted, separator_row_snr, separator_row_columns, df_snr_formatted, separator_row_crlb, separator_row_columns, df_crlb_formatted], axis=0)\n",
        "\n",
        "# Save the combined DataFrame to a CSV file\n",
        "combined_df.to_csv(current_path + '/derivatives/all_subjects/Metabolite_fit.csv', index=False)"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [
        "Phs6gCJ1lTI_",
        "O0IDwpEwlTJA",
        "NEuAWZjBlTJB",
        "7ovVOFv1ve2J",
        "KyXMxQD_kuYY",
        "Q0qCBPnAlTJD",
        "HoEPEgcgkuYY",
        "QWQO-9gCkuYZ",
        "OWxV_fOYkuYZ",
        "DyoZVUyblTJD",
        "edgU0UwLkuYa",
        "UF0j-YDBkuYb",
        "qOsM9FkgkuYb",
        "AVF25dC0lTJF",
        "mG21mJzQkuYc"
      ],
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.18"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}