7t_spine_coil_demo.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "7t-spine-coil-demo.ipynb",
      "private_outputs": true,
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "authorship_tag": "ABX9TyM7WIgVbG1xFiCfCDcel1qp",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/evaalonsoortiz/7t-spine-coil/blob/main/7t_spine_coil_demo.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iTKC2dBVcmLo"
      },
      "source": [
        "# **7T ACDC Spine Coil**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nbRwbnAHc6Wo"
      },
      "source": [
        "This Colab notebook is a demo for 7T ACDC Spine Coil Evaluation"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "l8PlGv4RdPNN"
      },
      "source": [
        "#Setup\n",
        "\n",
        "To setup, we will first clone the 7t-spine-coil and spinalcordtoolbox repos and then instal the spinalcordtoolbox.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5JQagzfBMLin"
      },
      "source": [
        "! git clone https://github.com/evaalonsoortiz/7t-spine-coil\n",
        "\n",
        "% cd 7t-spine-coil\n",
        "! pip install -r requirements.txt\n",
        "\n",
        "% cd ../\n",
        "\n",
        "!git clone --depth 1 --branch 5.3.0 https://github.com/spinalcordtoolbox/spinalcordtoolbox"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-NZsf_wm0LYP"
      },
      "source": [
        "# Go to directory and install SCT\n",
        "%cd spinalcordtoolbox/\n",
        "!yes | ./install_sct"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "f4v6DH9i0lXD"
      },
      "source": [
        "# Activate environment\n",
        "import os\n",
        "os.environ['PATH'] += ':/content/spinalcordtoolbox/bin'\n",
        "os.environ['SCT_DIR'] = '/content/spinalcordtoolbox'"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KsyLcaIpR6_A"
      },
      "source": [
        "%matplotlib inline\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import matplotlib.image as mpimg\n",
        "import plotly.express as px\n",
        "import nibabel as nib\n",
        "import sys\n",
        "import os\n",
        "import cv2\n",
        "import re\n",
        "import pandas as pd\n",
        "from os.path import join"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VRQSrBdld7m0"
      },
      "source": [
        "# Segment the spinal cord"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AtTR_dK51mBS"
      },
      "source": [
        "% cd ../7t-spine-coil\n",
        "! sct_deepseg_sc -i tfl_sag_2mmISO_384mm_REF720_0013/acdc_spine_7t_006_20210512_145613749_tfl_sag_2mmISO_384mm_REF720_20210512132653_13.nii -c t2s -qc tfl_sag_2mmISO_384mm_REF720_0013/qc\n",
        "seg = nib.load('acdc_spine_7t_006_20210512_145613749_tfl_sag_2mmISO_384mm_REF720_20210512132653_13_seg.nii')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hzZlOrWWrgPx"
      },
      "source": [
        "# Plot B1+ map with overlaid SC ROI"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "O4cCkHUi4BoH"
      },
      "source": [
        "img = nib.load('tfl_sag_2mmISO_384mm_REF720_0014/acdc_spine_7t_006_20210512_145613749_tfl_sag_2mmISO_384mm_REF720_20210512132653_14.nii')\n",
        "\n",
        "img_data = img.get_fdata()\n",
        "seg_data = seg.get_fdata()\n",
        "seg_data = np.ma.masked_where(seg_data < 0.3, seg_data)\n",
        "\n",
        "# Overlay the two images\n",
        "fig, ax = plt.subplots()\n",
        "plt.axis('off')\n",
        "img1 = ax.imshow(np.rot90(img_data[:,:,28]/10), cmap=plt.cm.gray)\n",
        "ax.imshow(np.rot90(seg_data[:,:,28]), cmap=plt.cm.jet, interpolation='none')\n",
        "plt.colorbar(img1)\n",
        "img1.set_clim(vmin=0, vmax=90)\n",
        "plt.title('B1+ [degrees]')\n",
        "plt.show()\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "T7RV31nxFAfF"
      },
      "source": [
        "# Extract B1+ Metrics"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ubR816PdFD2E"
      },
      "source": [
        "! sct_extract_metric -i tfl_sag_2mmISO_384mm_REF720_0014/acdc_spine_7t_006_20210512_145613749_tfl_sag_2mmISO_384mm_REF720_20210512132653_14.nii -f acdc_spine_7t_006_20210512_145613749_tfl_sag_2mmISO_384mm_REF720_20210512132653_13_seg.nii -method wa\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SL--zD3pFolN"
      },
      "source": [
        "## Plot Metrics"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XU2BypAfFrH1"
      },
      "source": [
        "# load metrics\n",
        "df = pd.read_csv(\"extract_metric.csv\", sep=\",\")\n",
        "\n",
        "df = df[(df != 0).all(1)]\n",
        "\n",
        "df['WA()'] = pd.to_numeric(df['WA()'], errors='coerce')\n",
        "df['WA()'] = df['WA()'].div(10)\n",
        "\n",
        "# plot average signal in degrees\n",
        "fig = px.line(df, y = 'WA()')\n",
        "\n",
        "fig.update_xaxes(title_text='slice [I->S]')\n",
        "fig.update_yaxes(title_text='B1+ [degrees]')\n",
        "\n",
        "fig.show()\n",
        "\n",
        "\n",
        "\n",
        "# plot average signal in [nT/V]\n",
        "\n",
        "# first compute voltage at socket\n",
        "# Open the file for reading\n",
        "with open('tfl_sag_2mmISO_384mm_REF720_0014/acdc_spine_7t_006_20210512_145613749_tfl_sag_2mmISO_384mm_REF720_20210512132653_14.json') as fd:\n",
        "\n",
        "    # Iterate over the lines\n",
        "    for line in fd:\n",
        "\n",
        "        # Capture one-or-more characters of non-whitespace after the initial match\n",
        "        match = re.search(r'\"TxRefAmp\":(.*?),', line)\n",
        "\n",
        "        # Did we find a match?\n",
        "        if match:\n",
        "            # Yes, process it\n",
        "            refV = float(match.group(1))\n",
        "\n",
        "SocketV = refV*10**-0.095  # [V]   \n",
        "gamma = 42.58e6 # [Hz/T]\n",
        "\n",
        "# plot average signal in nT/V\n",
        "# note: these results differ by a factor of 2 wrt to those that Kyle/Gaspard report. I am not sure where that factor of 2 comes from.\n",
        "fig = px.line(df['WA()'].mul((1/(360*gamma*0.001*SocketV))*10**9), y = 'WA()')\n",
        "\n",
        "fig.update_xaxes(title_text='slice [I->S]')\n",
        "fig.update_yaxes(title_text='B1+ [nT/V]')\n",
        "\n",
        "fig.show()\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6C7sG8LFXXdO"
      },
      "source": [
        "Investigating our missing factor 2 in $B_1$ scaling:\n",
        "$\\large FA (rad) = \\int\\limits_0^\\tau \\gamma B_1(t) t dt$\n",
        "\n",
        "\n",
        "From Kyle's $B_1$ mapping code:\n",
        "```\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",
        "```\n",
        "Assuming a 1ms $\\pi$-pulse at the acquisition voltage (with constant $B_1(t)$):\n",
        "$\\large B_1 = \\frac{\\pi}{\\gamma\\times 10^{-3}} = \\frac{\\pi}{2\\pi \\times 42.58 \\times 10^6 \\times 10^{-3}} = \\frac{1}{2 \\times 42.58 \\times 10^6 \\times 10^{-3}}$\n",
        "\n",
        "\n",
        "Scaling for any FA different from $\\pi$:\n",
        "$\\large B_1 =  \\frac{FA (rad)}{\\pi} \\times \\frac{1}{2 \\times 42.58 \\times 10^6 \\times 10^{-3}} = \\frac{FA (rad)}{2\\pi \\times 42.58 \\times 10^6 \\times 10^{-3}}$\n",
        "\n",
        "\n",
        "What's currently used in Kyle's code: \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",
        "RequestedFA = 90; % saturation flip angle -- hard-coded in sequence\n",
        "VoltageAtSocket = AcqVoltage .* 10^-0.095; % Voltage at the socket \n",
        "B1eff_mag = (AcquiredFA ./ RequestedFA) .* (pi ./ (GAMMA .* 1e-3 .* VoltageAtSocket)); % [T/V]\n",
        "B1eff_mag = B1eff_mag .* 1e9; % [T/V] to [nT/V]\n",
        "```\n",
        "(Ignoring voltage scaling and conversion to nT for clarity)\n",
        "$\\large B_1 = \\frac{FA(^{\\circ})}{Requested FA (90^{\\circ})} \\times \\frac{\\pi}{2\\pi \\times 42.58 \\times 10^6 \\times 10^{-3}}$\n",
        "\n",
        "Coversion to rad:\n",
        "\n",
        "$\\large B_1 = \\frac{FA(rad)}{\\frac{\\pi}{2}} \\times \\frac{\\pi}{2\\pi \\times 42.58 \\times 10^6 \\times 10^{-3}} = \\frac{FA(rad)}{\\pi \\times 42.58 \\times 10^6 \\times 10^{-3}}$\n",
        "\n",
        "A factor 1/2 is missing compared to our results."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xxULY9Bz_dsZ"
      },
      "source": [
        "# Plot SNR map with overlaid SC ROI"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iUFDAiBh_iSy"
      },
      "source": [
        "img = nib.load('coilQA_sag_FH_384mm_REF720_1p5x1p5x3mm/acdc_spine_7t_006_20210512_145613749_coilQA_sag_FH_384mm_REF720_1p5x1p5x3mm_20210512132653_24.nii')\n",
        "\n",
        "img_data = img.get_fdata()\n",
        "\n",
        "fig = px.imshow(np.rot90(img_data[:,:,3,0]),color_continuous_scale='gray',zmin=0, zmax=200)\n",
        "fig.update_layout(\n",
        "    height=800, width = 800,\n",
        "    title_text='SNR (1.5 x 1.5 x 3mm^3)'\n",
        ")\n",
        "fig.show()\n"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}