diff --git a/.DS_Store b/.DS_Store
index 7239436..a1b93e9 100644
Binary files a/.DS_Store and b/.DS_Store differ
diff --git a/05-B0/README.md b/05-B0/README.md
new file mode 100644
index 0000000..31287b9
--- /dev/null
+++ b/05-B0/README.md
@@ -0,0 +1,2 @@
+# mooc_b0_mapping
+Notebooks to create Figures for the B0 mapping chapter of the MOOC
diff --git a/05-B0/binder/postBuild b/05-B0/binder/postBuild
new file mode 100644
index 0000000..e69de29
diff --git a/05-B0/binder/requirements.txt b/05-B0/binder/requirements.txt
new file mode 100644
index 0000000..e13b1a0
--- /dev/null
+++ b/05-B0/binder/requirements.txt
@@ -0,0 +1,7 @@
+jupyter-book==0.13.0
+lxml_html_clean
+nibabel
+pandas
+plotly
+scipy
+tabulate
diff --git a/05-B0/binder/runtime.txt b/05-B0/binder/runtime.txt
new file mode 100644
index 0000000..3cd3102
--- /dev/null
+++ b/05-B0/binder/runtime.txt
@@ -0,0 +1 @@
+python-3.8
\ No newline at end of file
diff --git a/05-B0/content/_config.yml b/05-B0/content/_config.yml
new file mode 100644
index 0000000..b156fa5
--- /dev/null
+++ b/05-B0/content/_config.yml
@@ -0,0 +1,37 @@
+# Book settings
+# Learn more at https://jupyterbook.org/customize/config.html
+
+title: B0 mapping
+logo: logo.png
+
+# Force re-execution of notebooks on each build.
+# See https://jupyterbook.org/content/execute.html
+execute:
+ execute_notebooks: force
+
+launch_buttons:
+ notebook_interface: "jupyterlab"
+
+# Information about where the book exists on the web
+repository:
+ url: https://github.com/qMRLab/qMRI-MOOC_B0-module # Online location of your book
+ path_to_book: content # Optional path to your book, relative to the repository root
+ branch: main # Which branch of the repository should be used when creating links (optional)
+
+only_build_toc_files: true
+
+execute:
+ timeout: 600
+
+sphinx:
+ local_extensions : # A list of local extensions to load by sphinx specified by "name: path" items
+ recursive_update : false # A boolean indicating whether to overwrite the Sphinx config (true) or recursively update
+ config : # key-value pairs to directly over-ride the Sphinx configuration
+ bibtex_reference_style: author_year
+ bibtex_default_style: plain
+ bibtex_tooltips: true
+ html_static_path: ['_static']
+ html_css_files: ['custom.css']
+
+bibtex_bibfiles:
+ - ../paper.bib
diff --git a/05-B0/content/_static/custom.css b/05-B0/content/_static/custom.css
new file mode 100644
index 0000000..905f567
--- /dev/null
+++ b/05-B0/content/_static/custom.css
@@ -0,0 +1,52 @@
+p {
+ text-align: justify;
+ }
+ a {
+ color: #870000!important;
+ }
+ a:hover {
+ color: brick!important;
+ }
+ a:visited {
+ color: #6f42c1!important;
+ }
+ body, h1, h2, h3, h4, h5 {
+ font-family: 'STIX Two Text';
+ }
+ h1 {
+ color: #342727!important;
+ }
+ .caption {
+ text-align:justify;
+ line-height:1.25;
+ font-size:80%
+ }
+
+ #site-navigation {
+ display: none;
+ }
+
+ #main-content{
+ margin: 0 auto;
+ }
+
+
+ .fa, .far, .fas {
+ color: red;
+ }
+
+ .header-article::after{
+ content: url('https://neurolibre.org/assets/joss2-b7ea8dd1f24ce2b0825db17e3495576253ba4c64df3384a66354713d76c40a75.svg');
+ width: 150px;
+ height: 80px;
+ position: relative;
+ right: 10px;
+ }
+
+ .header-article::before{
+ content: url('https://raw.githubusercontent.com/neurolibre/brand/b5e10e051a0059824ca2f4204e1288927e551d49/svg/neurolibre_logo.svg');
+ width: 33px;
+ height: 33px;
+ position: absolute;
+ left: 10px;
+ }
\ No newline at end of file
diff --git a/05-B0/content/_toc.yml b/05-B0/content/_toc.yml
new file mode 100644
index 0000000..bb3e2d9
--- /dev/null
+++ b/05-B0/content/_toc.yml
@@ -0,0 +1,5 @@
+# Table of contents
+# Learn more at https://jupyterbook.org/customize/toc.html
+
+format: jb-article
+root: paper
diff --git a/05-B0/content/assets/before_after.css b/05-B0/content/assets/before_after.css
new file mode 100644
index 0000000..a764c37
--- /dev/null
+++ b/05-B0/content/assets/before_after.css
@@ -0,0 +1,50 @@
+.beforeAfter {
+ border-radius: 5px;
+ /* border: 2px solid silver; */
+ position: relative;
+ margin: 0;
+ padding: 0;
+}
+
+.after {
+ height:100%;
+ width:100%;
+ z-index: 1;
+ position: absolute;
+ }
+
+.before {
+ height: 100%;
+ width: 50%;
+ overflow: hidden;
+ border-right: 4px solid silver;
+ z-index: 40;
+ position: absolute;
+}
+
+.middle {
+ top: 45%;
+ left: 48%;
+ z-index: 50;
+ draggable: true;
+ position: absolute;
+}
+
+.middle:before {
+ content: '\f1de';
+ font-family: "Font Awesome 6 Free";
+ font-size: 30pt;
+ font-weight: 900;
+ background-color: silver;
+ border: 1pt solid black;
+ cursor: grab;
+ border-radius: 5pt;
+}
+
+.middle:hover:before {
+ filter: brightness(1.15)
+}
+
+.middle:active:before {
+ cursor: grabbing;
+}
\ No newline at end of file
diff --git a/05-B0/content/assets/before_after.js b/05-B0/content/assets/before_after.js
new file mode 100644
index 0000000..961f5b3
--- /dev/null
+++ b/05-B0/content/assets/before_after.js
@@ -0,0 +1,57 @@
+function beforeAfter() {
+ $(".middle").on("mousedown", function () {
+ $(document).on("mouseup",function() {
+ $(document).unbind("mousemove")
+ $(document).unbind("mouseup")
+ oldX = undefined
+ oldY = undefined
+ })
+ $(document).on("mousemove", function () {
+ adjustSize()
+ })
+ })
+
+ $(".before .modebar").css({'left':'0px', 'position':'absolute'})
+}
+
+window.fetch = new Proxy(window.fetch, {
+ apply(fetch, that, args) {
+ // Forward function call to the original fetch
+ const result = fetch.apply(that, args);
+
+ // Do whatever you want with the resulting Promise
+ result.then((response) => {
+ if (args[0] == '/_dash-update-component') {
+ setTimeout(function() {beforeAfter()}, 500)
+ }})
+ return result
+ }
+ }
+ )
+
+var oldX;
+var oldY;
+function adjustSize() {
+ if (event.screenX != 0) {
+ var clientX = event.clientX
+ var scrollLeft = $(document).scrollLeft()
+
+ if (oldX && !((parseFloat($('.before').css("width"))+((clientX-scrollLeft)-oldX)) > $('.beforeAfter').width())) {
+ var width = parseFloat($('.before').css("width"))+((clientX-scrollLeft)-oldX)
+ $('.before').css("width", width)
+ $('.middle').css("left", ($('.before').width()/$('.beforeAfter').width()*100)-2 + '%')
+
+ }
+
+ if (parseFloat($('.middle').css('left'))> $('.beforeAfter').width()) {
+ $('.middle').css("left", '99.5%')
+ $('.before').css("width", $('.beforeAfter').width())
+ }
+ }
+ oldX = clientX-scrollLeft
+ event.target.style.cursor = 'grab';
+}
+
+$(document).ready(function() {
+ setTimeout(function() {beforeAfter()}, 1000)
+})
\ No newline at end of file
diff --git a/05-B0/content/chapter4_1.ipynb b/05-B0/content/chapter4_1.ipynb
new file mode 100644
index 0000000..4fc12d0
--- /dev/null
+++ b/05-B0/content/chapter4_1.ipynb
@@ -0,0 +1,516 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1fbfb83d-7790-47df-a66e-e92e64013190",
+ "metadata": {
+ "editable": true,
+ "slideshow": {
+ "slide_type": ""
+ },
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "import math\n",
+ "import json\n",
+ "import nibabel as nib\n",
+ "import numpy as np\n",
+ "from numpy.fft import ifftn, fftn, ifft, fftshift, ifftshift\n",
+ "import os\n",
+ "import plotly.express as px\n",
+ "import plotly.graph_objects as go\n",
+ "from plotly.subplots import make_subplots\n",
+ "from scipy.signal import butter, lfilter, freqz, filtfilt\n",
+ "from scipy.io import loadmat\n",
+ "import warnings\n",
+ "PI_UNICODE = \"\\U0001D70B\"\n",
+ "CHI_UNICODE = \"\\U0001D712\"\n",
+ "MICRO_UNICODE = \"\\u00B5\"\n",
+ "GYRO_BAR_RATIO_H = 42.6e6 # [Hz/T]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6aab6266-6443-4118-b61b-2dc6b796856b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fname_mag_e1 = os.path.join(\"fmap\", \"sub-fmap_magnitude1.nii.gz\")\n",
+ "fname_mag_e1_json = os.path.join(\"fmap\", \"sub-fmap_magnitude1.json\")\n",
+ "fname_phase_e1 = os.path.join(\"fmap\", \"sub-fmap_phase1.nii.gz\")\n",
+ "fname_phase_e2 = os.path.join(\"fmap\", \"sub-fmap_phase2.nii.gz\")\n",
+ "fname_mask = os.path.join(\"fmap\", \"mask.nii.gz\")\n",
+ "fname_fmap = os.path.join(\"fmap\", \"fmap.nii.gz\")\n",
+ "mag_e1 = nib.load(fname_mag_e1).get_fdata()[30:-30,8:108,30]\n",
+ "phase_e1 = nib.load(fname_phase_e1).get_fdata()[30:-30,8:108,30]\n",
+ "phase_e2 = nib.load(fname_phase_e2).get_fdata()[30:-30,8:108,30]\n",
+ "fmap_hz = nib.load(fname_fmap).get_fdata()[30:-30,8:108,30]\n",
+ "# mag_e1 = nib.load(fname_mag_e1).get_fdata()\n",
+ "# phase_e1 = nib.load(fname_phase_e1).get_fdata()\n",
+ "# phase_e2 = nib.load(fname_phase_e2).get_fdata()\n",
+ "# mask = nib.load(fname_mask).get_fdata()\n",
+ "# fmap_hz = nib.load(fname_fmap).get_fdata()\n",
+ "\n",
+ "with open(fname_mag_e1_json, 'r') as json_data:\n",
+ " data = json.load(json_data)\n",
+ "\n",
+ "freq = data[\"ImagingFrequency\"]\n",
+ "fmap_t = fmap_hz / GYRO_BAR_RATIO_H * 1e6\n",
+ "fmap_ppm = fmap_hz / freq\n",
+ "\n",
+ "fig = go.Figure()\n",
+ "fig.add_trace(go.Heatmap(z=np.rot90(mag_e1, k=-1),\n",
+ " colorscale='gray',\n",
+ " colorbar=dict(\n",
+ " title=\"a.u.\",\n",
+ " titleside=\"top\",\n",
+ " tickmode=\"array\"\n",
+ " ))\n",
+ " )\n",
+ "fig.add_trace(go.Heatmap(z=np.rot90(phase_e2 / 4095 * 2 * math.pi - math.pi, k=-1),\n",
+ " colorscale='gray',\n",
+ " colorbar=dict(\n",
+ " title=\"Rad\",\n",
+ " titleside=\"top\",\n",
+ " tickmode=\"array\",\n",
+ " tickvals=[-math.pi, 0, math.pi-0.01],\n",
+ " ticktext = [f\"-{PI_UNICODE}\", 0, f'{PI_UNICODE}']\n",
+ " ),\n",
+ " visible=False))\n",
+ "fig.add_trace(go.Heatmap(z=np.rot90(fmap_hz, k=-1),\n",
+ " colorscale='gray',\n",
+ " colorbar=dict(\n",
+ " title=\"Hz\",\n",
+ " titleside=\"top\",\n",
+ " tickmode=\"array\"\n",
+ " ),\n",
+ " visible=False))\n",
+ "fig.add_trace(go.Heatmap(z=np.rot90(fmap_t, k=-1),\n",
+ " colorscale='gray',\n",
+ " colorbar=dict(\n",
+ " title=f\"{MICRO_UNICODE}T\",\n",
+ " titleside=\"top\",\n",
+ " tickmode=\"array\"\n",
+ " ),\n",
+ " visible=False))\n",
+ "fig.add_trace(go.Heatmap(z=np.rot90(fmap_ppm, k=-1),\n",
+ " colorscale='gray',\n",
+ " colorbar=dict(\n",
+ " title=\"ppm\",\n",
+ " titleside=\"top\",\n",
+ " tickmode=\"array\"\n",
+ " ),\n",
+ " visible=False))\n",
+ "\n",
+ "\n",
+ "\n",
+ "x0=0\n",
+ "y0=89\n",
+ "x1=10\n",
+ "y1=99\n",
+ "h = 2\n",
+ "rounded_bottom_left = f' M {x0+h}, {y0} Q {x0}, {y0} {x0}, {y0+h}'#\n",
+ "rounded_top_left = f' L {x0}, {y1-h} Q {x0}, {y1} {x0+h}, {y1}'\n",
+ "rounded_top_right = f' L {x1-h}, {y1} Q {x1}, {y1} {x1}, {y1-h}'\n",
+ "rounded_bottom_right = f' L {x1}, {y0+h} Q {x1}, {y0} {x1-h}, {y0}Z'\n",
+ "path = rounded_bottom_left + rounded_top_left+\\\n",
+ " rounded_top_right+rounded_bottom_right\n",
+ "\n",
+ "annotations = ['A', 'B', 'C', 'D', 'E']\n",
+ "shapes = []\n",
+ "for i_shape, annotation in enumerate(annotations):\n",
+ " shapes.append(dict(type='path',\n",
+ " path=path,\n",
+ " fillcolor='white',\n",
+ " layer='above',\n",
+ " line=dict(width=1),\n",
+ " label=dict(text=f\"{annotation}\")\n",
+ " )\n",
+ " )\n",
+ "\n",
+ "fig.add_shape(shapes[0])\n",
+ "# Add dropdown\n",
+ "fig.update_layout(\n",
+ " title_text=\"Magnitude\",\n",
+ " title_x=0.5,\n",
+ " height=500,\n",
+ " width=600,\n",
+ " updatemenus=[\n",
+ " dict(\n",
+ " buttons=list([\n",
+ " dict(\n",
+ " method=\"update\",\n",
+ " args=[{\"visible\": [True, False, False, False, False]},\n",
+ " {'shapes': [shapes[0]], \"title\": \"Magnitude\"}],\n",
+ " label=\"Magnitude\",\n",
+ " ),\n",
+ " dict(\n",
+ " method=\"update\",\n",
+ " args=[{\"visible\": [False, True, False, False, False]},\n",
+ " {'shapes': [shapes[1]], \"title\": \"Phase\"}],\n",
+ " label=\"Phase\",\n",
+ " \n",
+ " ),\n",
+ " dict(\n",
+ " method=\"update\",\n",
+ " args=[{\"visible\": [False, False, True, False, False]},\n",
+ " {'shapes': [shapes[2]], \"title\": \"B0 Fieldmap (Hz)\"}],\n",
+ " label=\"B0 field map (Hz)\",\n",
+ " ),\n",
+ " dict(\n",
+ " method=\"update\",\n",
+ " args=[{\"visible\": [False, False, False, True, False]},\n",
+ " {'shapes': [shapes[3]], \"title\": f\"B0 Fieldmap ({MICRO_UNICODE}Tesla)\"}],\n",
+ " label=f\"B0 field map ({MICRO_UNICODE}Tesla)\",\n",
+ " ),\n",
+ " dict(\n",
+ " method=\"update\",\n",
+ " args=[{\"visible\": [False, False, False, False, True]},\n",
+ " {'shapes': [shapes[4]], \"title\": \"B0 Fieldmap (ppm)\"}],\n",
+ " label=\"B0 field map (ppm)\",\n",
+ " )\n",
+ " ]),\n",
+ " direction=\"down\",\n",
+ " showactive=True,\n",
+ " )\n",
+ " ]\n",
+ ")\n",
+ "fig.update_xaxes(showticklabels=False)\n",
+ "fig.update_yaxes(showticklabels=False)\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "321d843b-3df3-47c1-b084-3e895f23ac0d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def dipole_kernel(b0_dir, voxel_size, n_voxels):\n",
+ " \"\"\" Create a dipole kernel\n",
+ " dipole kernel: (3*cos(theta)**2 - 1) / (4*pi*r**3)\n",
+ " => (3*r**2*cos(theta)**2 - r**2) / (4*pi*r**5)\n",
+ " => (3*b0_dir**2 - r**2) / (4*pi*r**2**2.5)\n",
+ "\n",
+ " Function inspired and derived from: https://onlinelibrary.wiley.com/doi/10.1002/mrm.28716\n",
+ " \"\"\"\n",
+ " eps = 0.00001\n",
+ " x, y, z = np.meshgrid(range(round(-n_voxels[0]/2+0.5), round(n_voxels[0]/2+0.5)), range(round(-n_voxels[1]/2+0.5), round(n_voxels[1]/2+0.5)), range(round(-n_voxels[2]/2+0.5), round(n_voxels[2]/2+0.5)), indexing='ij')\n",
+ "\n",
+ " x = x * voxel_size[0] + eps\n",
+ " y = y * voxel_size[1] + eps\n",
+ " z = z * voxel_size[2] + eps\n",
+ "\n",
+ " r2 = (x**2 + y**2 + z**2)\n",
+ "\n",
+ " d = np.prod(voxel_size) * ( 3 * ((x*b0_dir[0] + y*b0_dir[1] + z*b0_dir[2])**2) - r2 ) / (4 * math.pi * r2**2.5)\n",
+ "\n",
+ " d[np.isnan(d)] = eps\n",
+ " D = np.real(fftshift(fftn(ifftshift(d))))\n",
+ "\n",
+ " mid_voxel = n_voxels[0]//2\n",
+ " return d[n_voxels[1]//2], D[n_voxels[1]//2]\n",
+ "\n",
+ "b0_dir = (0, 0, 1)\n",
+ "voxel_size = np.array((1, 1, 1)) * 1e-3\n",
+ "n_voxels = (201,201,201)\n",
+ "d, D = dipole_kernel(b0_dir, voxel_size, n_voxels)\n",
+ "\n",
+ "fig = go.Figure()\n",
+ "fig = make_subplots(rows=1, cols=2, shared_xaxes=False, horizontal_spacing=0.13, vertical_spacing = 0.12, subplot_titles=(\"Dipole Kernel (d)\", \"Dipole Kernel (D)\"), specs=[[{\"type\": \"Heatmap\"}, {\"type\": \"Heatmap\"}]])\n",
+ "fig.add_trace(go.Heatmap(z=d, colorscale='gray', showscale=False, zmin=-1e-6, zmax=1e-6))\n",
+ "fig.add_trace(go.Heatmap(z=D, colorscale='gray', showscale=False), 1, 2)\n",
+ "fig.update_xaxes(showticklabels=False)\n",
+ "fig.update_yaxes(showticklabels=False)\n",
+ "fig.update_layout(\n",
+ " height=500,\n",
+ " width=900)\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "fcc2c8b3-3e90-4c98-b8f3-9ad58138ba0a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Load cylinder (Y 64)\n",
+ "susc = nib.load(os.path.join(\"field_simulations\", \"cylinder\", \"Chi.nii.gz\")).get_fdata()\n",
+ "fmap_hz_all = nib.load(os.path.join(\"field_simulations\", \"cylinder\", \"fmap_hz.nii.gz\")).get_fdata()\n",
+ "local_field_cyl = nib.load(os.path.join(\"field_simulations\", \"cylinder\", \"local_field.nii.gz\")).get_fdata()\n",
+ "\n",
+ "# Load brain (Z 210)\n",
+ "susc_brain = nib.load(os.path.join(\"field_simulations\", \"brain\", \"chi_masked.nii.gz\")).get_fdata()\n",
+ "fmap_hz_brain_all = nib.load(os.path.join(\"field_simulations\", \"brain\", \"fmap_masked.nii.gz\")).get_fdata()\n",
+ "local_field_brain = nib.load(os.path.join(\"field_simulations\", \"brain\", \"local_field.nii.gz\")).get_fdata()\n",
+ "\n",
+ "fig = make_subplots(rows=1, cols=3, shared_xaxes=False, horizontal_spacing=0.13, vertical_spacing = 0.12, subplot_titles=(f\"Susceptibility distribution ({CHI_UNICODE})\", \"Simulated B0 map\", \"Simulated B0 map
no background field\"), specs=[[{\"type\": \"Heatmap\"}, {\"type\": \"Heatmap\"}, {\"type\": \"Heatmap\"}]])\n",
+ "fig.add_trace(go.Heatmap(z=susc, colorscale='gray', colorbar_x=1/3 - 0.09, colorbar=dict(title=\"ppm\", titleside=\"top\")), 1, 1)\n",
+ "fig.add_trace(go.Heatmap(z=fmap_hz_all, colorscale='gray', colorbar_x=2/3 - 0.045, colorbar=dict(title=\"Hz\", titleside=\"top\")), 1, 2)\n",
+ "fig.add_trace(go.Heatmap(z=local_field_cyl, colorscale='gray', colorbar_x=1-0.004, colorbar=dict(title=\"Hz\", titleside=\"top\")), 1, 3)\n",
+ "fig.add_trace(go.Heatmap(z=np.rot90(susc_brain, k=-1), colorscale='gray', colorbar_x=1/3 - 0.09, zmin=-0.5, zmax=0.5, colorbar=dict(title=\"ppm\", titleside=\"top\"), visible=False), 1, 1)\n",
+ "fig.add_trace(go.Heatmap(z=np.rot90(fmap_hz_brain_all, k=-1), colorscale='gray', colorbar_x=2/3 - 0.045, zmin=1100, zmax=2300, colorbar=dict(title=\"Hz\", titleside=\"top\"), visible=False), 1, 2)\n",
+ "fig.add_trace(go.Heatmap(z=np.rot90(local_field_brain, k=-1), colorscale='gray', zmin=-4, zmax=4, colorbar_x=1-0.004, colorbar=dict(title=\"Hz\", titleside=\"top\"), visible=False), 1, 3)\n",
+ "\n",
+ "### Create buttons for drop down menu\n",
+ "labels = [\"Cylinders\", \"Brain\"]\n",
+ "buttons = []\n",
+ "for i, label in enumerate(labels):\n",
+ " if label == \"Cylinders\":\n",
+ " visibility = [True, True, True, False, False, False]\n",
+ " else:\n",
+ " visibility = [False, False, False, True, True, True]\n",
+ " button = dict(\n",
+ " label = label,\n",
+ " method = 'update',\n",
+ " args = [{'visible': visibility},\n",
+ " {'title': label}])\n",
+ " buttons.append(button)\n",
+ "\n",
+ "updatemenus = list([\n",
+ " dict(active=0,\n",
+ " x=-0.15,\n",
+ " buttons=buttons,\n",
+ " showactive=True,\n",
+ " )\n",
+ "])\n",
+ "\n",
+ "fig.update_xaxes(showticklabels=False)\n",
+ "fig.update_yaxes(showticklabels=False)\n",
+ "fig.update_layout({\"height\": 380, \"width\": 900},\n",
+ " title_text=\"B0 maps from susceptibility simulations\",\n",
+ " title_x=0.5,\n",
+ " updatemenus=updatemenus,\n",
+ " showlegend=False\n",
+ " )\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1aaf54a3-7715-4379-9ff2-874bb6e2d6b2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Note: Field was reduced a lot to be able to show the sinusoid\n",
+ "# Note: *2 after lowpass filter is because this is a single coil (sin instead of e^(-ix)) and demodulating by multiplying a sinusoid creates a 1/2 difference. In practice, since we have both x and y components, we can recover the full signal instead of doing X2.\n",
+ "GYRO_BAR_RATIO_H = 42.6e6 # [Hz/T]\n",
+ "b0 = 0.000002 # [T]\n",
+ "T2 = 0.3 # s\n",
+ "y_0_cst = 100\n",
+ "fs = 10000\n",
+ "\n",
+ "f_larmor = b0 * GYRO_BAR_RATIO_H\n",
+ "t = np.linspace(0, 1, fs + 1) # 1 second\n",
+ "\n",
+ "def butter_lowpass(cutoff, fs, order=5):\n",
+ " return butter(order, cutoff, fs=fs, btype='low', analog=False)\n",
+ "\n",
+ "def butter_lowpass_filter(data, cutoff, fs, order=5):\n",
+ " b, a = butter_lowpass(cutoff, fs, order=order)\n",
+ " y = filtfilt(b, a, data, method='gust')\n",
+ " return y\n",
+ "\n",
+ "# Lab frame\n",
+ "y_0 = y_0_cst * np.sin(2 * math.pi * f_larmor * t)\n",
+ "exp = np.exp(-t/T2)\n",
+ "y = y_0 * exp / y_0_cst\n",
+ "temp = y * (y_0 / y_0_cst)\n",
+ "y_demod = butter_lowpass_filter(temp, f_larmor, fs, order=5) * 2\n",
+ "\n",
+ "fig = go.Figure()\n",
+ "fig.add_scatter(x=t, y=y, name=\"FID\")\n",
+ "fig.add_scatter(x=t, y=y_demod, name=\"Demodulated FID\")\n",
+ "fig.add_scatter(x=t, y=exp, name=\"T2 decay\")\n",
+ "\n",
+ "# 2 isochromats\n",
+ "y_1amp = y_0_cst / 10\n",
+ "y_1 = y_1amp * np.sin(2 * math.pi * (f_larmor + 10) * t)\n",
+ "y = (y_0 + y_1) * exp / (y_0_cst + y_1amp)\n",
+ "temp = y * (y_0 / y_0_cst)\n",
+ "y_demod = butter_lowpass_filter(temp, f_larmor, fs, order=5) * 2\n",
+ "\n",
+ "fig.add_scatter(x=t, y=y, name=\"FID\", visible=False)\n",
+ "fig.add_scatter(x=t, y=y_demod, name=\"Demodulated FID\", visible=False)\n",
+ "fig.add_scatter(x=t, y=exp, name=\"T2 decay\", visible=False)\n",
+ "\n",
+ "# Multiple isochromats\n",
+ "n_freqs = 100\n",
+ "fid = y_0 * exp\n",
+ "y_sum_demod = butter_lowpass_filter(fid * (y_0 / y_0_cst), f_larmor, fs, order=5) * 2\n",
+ "for i in range(n_freqs):\n",
+ " amp = 1\n",
+ " freq_offset = 10\n",
+ " scale = freq_offset/n_freqs\n",
+ " mid = n_freqs // 2\n",
+ " y_1 = amp * np.sin(2 * math.pi * (f_larmor + scale*(mid - i)) * t) * exp\n",
+ " fid += y_1\n",
+ " y_demod = butter_lowpass_filter(y_1 * (y_0 / y_0_cst), f_larmor, fs, order=5) * 2\n",
+ " y_sum_demod += y_demod\n",
+ "\n",
+ "y_demod_scaled = y_sum_demod / (y_0_cst + (n_freqs * amp))\n",
+ "fid_scaled = fid / (y_0_cst + (n_freqs * amp))\n",
+ "\n",
+ "fig.add_scatter(x=t, y=fid_scaled, name=\"FID\", visible=False)\n",
+ "fig.add_scatter(x=t, y=y_demod_scaled, name=\"Demodulated FID\", visible=False)\n",
+ "fig.add_scatter(x=t, y=exp, name=\"T2 decay\", visible=False)\n",
+ "fig.update_traces(marker=dict(size=3))\n",
+ "fig.update_layout(\n",
+ " title=\"Single species\",\n",
+ " title_x=0.5,\n",
+ " updatemenus=[\n",
+ " dict(\n",
+ " buttons=list([\n",
+ " dict(\n",
+ " args=[{\"visible\": [True, True, True, False, False, False, False, False, False]},\n",
+ " {\"title\": \"Single species\"}],\n",
+ " label=\"Single species\",\n",
+ " method=\"update\"\n",
+ " ),\n",
+ " dict(\n",
+ " args=[{\"visible\": [False, False, False, True, True, True, False, False, False]},\n",
+ " {\"title\": \"Two species\"}],\n",
+ " label=\"Two species\",\n",
+ " method=\"update\"\n",
+ " ),\n",
+ " dict(\n",
+ " args=[{\"visible\": [False, False, False, False, False, False, True, True, True]},\n",
+ " {\"title\": \"Multiple Species\"}],\n",
+ " label=\"Multiple Species\",\n",
+ " method=\"update\"\n",
+ " )\n",
+ " ]),\n",
+ " direction=\"down\",\n",
+ " showactive=True,\n",
+ "\n",
+ " ),\n",
+ " ])\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "720b531c-3095-4ded-beb4-ea3ec8ac78cd",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def calc_dk(gx, gy, dt):\n",
+ " dkx = GYRO_BAR_RATIO_H * gx * dt\n",
+ " dky = GYRO_BAR_RATIO_H * gy * dt\n",
+ " return (dkx, dky)\n",
+ "\n",
+ "gy_bad = 100e-6\n",
+ "\n",
+ "end_time = 0.0912\n",
+ "n_times = 913\n",
+ "t = np.linspace(0, end_time, n_times)\n",
+ "dt = end_time / n_times\n",
+ "nx = 64\n",
+ "k = np.zeros([n_times, 2])\n",
+ "k_distorted = np.zeros([n_times, 2])\n",
+ "for it in range(1, n_times):\n",
+ " \n",
+ " if it <= 20:\n",
+ " gx = -40e-3\n",
+ " gy = -40e-3\n",
+ " else:\n",
+ " n_steps = 138\n",
+ " i = (it - 20) % n_steps\n",
+ " if i <= 0:\n",
+ " gx = 0\n",
+ " gy = 25e-3\n",
+ " elif i <= nx:\n",
+ " gx = 25e-3\n",
+ " gy = 0\n",
+ " elif i <= nx + 5:\n",
+ " gx = 0\n",
+ " gy = 25e-3\n",
+ " elif i <= (2*nx) + 5:\n",
+ " gx = -25e-3\n",
+ " gy = 0\n",
+ " elif i <= n_steps:\n",
+ " gx = 0\n",
+ " gy = 25e-3\n",
+ "\n",
+ " dkx, dky = calc_dk(gx, gy, dt)\n",
+ " kx = k[it - 1, 0] + dkx\n",
+ " ky = k[it - 1, 1] + dky\n",
+ " k[it, :] = [kx, ky]\n",
+ " dkx, dky = calc_dk(gx, gy + gy_bad, dt)\n",
+ " kx = k_distorted[it - 1, 0] + dkx\n",
+ " ky = k_distorted[it - 1, 1] + dky\n",
+ " k_distorted[it, :] = [kx, ky]\n",
+ "\n",
+ "\n",
+ "fig = go.Figure()\n",
+ "fig.add_trace(go.Scatter(x=k[:, 0], y=k[:, 1],\n",
+ " mode='lines',\n",
+ " line=dict(color='#636EFA'),\n",
+ " name='Theoretical trajectory'))\n",
+ "fig.add_trace(go.Scatter(x=k_distorted[:, 0], y=k_distorted[:, 1],\n",
+ " mode='lines',\n",
+ " line=dict(color='#fa6363'),\n",
+ " name='Inhomogeneous trajectory'))\n",
+ "frames = [dict(\n",
+ " data=[go.Scatter(x=k[:2*i, 0], y=k[:2*i, 1],\n",
+ " mode='lines',\n",
+ " line=dict(color='#636EFA'),\n",
+ " name='Theoretical trajectory'),\n",
+ " go.Scatter(x=k_distorted[:2*i, 0], y=k_distorted[:2*i, 1],\n",
+ " mode='lines',\n",
+ " line=dict(color='#fa6363'),\n",
+ " name='Inhomogeneous trajectory')],\n",
+ " name=str(i),\n",
+ " traces=[0,1]) for i in range(int(n_times/2))]\n",
+ "fig.frames = frames\n",
+ "\n",
+ "fig.update_xaxes(range=[-3500, 3500])\n",
+ "fig.update_yaxes(range=[-3700, 3700])\n",
+ "fig.update_xaxes(showticklabels=False)\n",
+ "fig.update_yaxes(showticklabels=False)\n",
+ "fig.update_layout(title=\"K-space Trajectory\",\n",
+ " title_x=0.5,\n",
+ " height=600,\n",
+ " width=700,\n",
+ " updatemenus=[dict(\n",
+ " type='buttons',\n",
+ " buttons=[dict(label='Play',\n",
+ " method='animate',\n",
+ " args=[None, dict(frame=dict(duration=15, redraw=False), transition=dict(duration=15), fromcurrent=True, mode='immediate')]),\n",
+ " dict(label='Pause',\n",
+ " method='animate',\n",
+ " args=[[None], dict(frame=dict(duration=0, redraw=False), mode='immediate')])\n",
+ " ])])\n",
+ "fig.show()"
+ ]
+ }
+ ],
+ "metadata": {
+ "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.0"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/05-B0/content/chapter4_2.ipynb b/05-B0/content/chapter4_2.ipynb
new file mode 100644
index 0000000..a67b3d0
--- /dev/null
+++ b/05-B0/content/chapter4_2.ipynb
@@ -0,0 +1,473 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "63c39eb3-2ac0-4456-a905-d6981285a9b6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import plotly.express as px\n",
+ "import os\n",
+ "import nibabel as nib\n",
+ "import numpy as np\n",
+ "from plotly.subplots import make_subplots\n",
+ "import plotly.graph_objects as go\n",
+ "import math\n",
+ "import json\n",
+ "PI_UNICODE = \"\\U0001D70B\""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "21fa46b1-3008-42fd-a05e-a09e0a5fb5be",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Parameters\n",
+ "num_frames = 200\n",
+ "omega_0 = 1 # Larmor frequency\n",
+ "omega_1 = 0.9 # Inhomogeneous spin\n",
+ "time_max = 5 # [s]\n",
+ "\n",
+ "# Initial phase of the spin\n",
+ "initial_phase = 0.5\n",
+ "\n",
+ "# Time array\n",
+ "time = np.linspace(0, time_max, num_frames)\n",
+ "\n",
+ "# Generate data for spins\n",
+ "x = np.cos(omega_0 * (2*math.pi) * time + initial_phase)\n",
+ "y = np.sin(omega_0 * (2*math.pi) * time + initial_phase)\n",
+ "x1 = np.cos(omega_1 * (2*math.pi) * time + initial_phase)\n",
+ "y1 = np.sin(omega_1 * (2*math.pi) * time + initial_phase)\n",
+ "\n",
+ "# Generate data for spins in rotating frame of reference\n",
+ "x_rot = np.cos((omega_0-omega_0) * (2*math.pi) * time + initial_phase)\n",
+ "y_rot = np.sin((omega_0-omega_0) * (2*math.pi) * time + initial_phase)\n",
+ "x1_rot = np.cos((omega_1-omega_0) * (2*math.pi) * time + initial_phase)\n",
+ "y1_rot = np.sin((omega_1-omega_0) * (2*math.pi) * time + initial_phase)\n",
+ "\n",
+ "# Calculate angles\n",
+ "angles = (np.arctan2(y,x))\n",
+ "angles1 = (np.arctan2(y1,x1))\n",
+ "angles_rot = (np.arctan2(y_rot,x_rot))\n",
+ "angles1_rot = (np.arctan2(y1_rot,x1_rot))\n",
+ "\n",
+ "# Create figure\n",
+ "fig = make_subplots(rows=1, cols=2, shared_xaxes=False, horizontal_spacing=0.1, subplot_titles=(\"Spin Rotating Through Time\", \"Phase Evolution of the Signal (rad)\"))\n",
+ "\n",
+ "# Add spin as an arrow\n",
+ "fig.add_trace(go.Scatter(\n",
+ " x=[0, x[0]],\n",
+ " y=[0, y[0]],\n",
+ " mode='lines+markers',\n",
+ " marker=dict(size=5),\n",
+ " line=dict(color='blue', width=3),\n",
+ " name='Spin at f0'),\n",
+ " row=1, col=1)\n",
+ "fig.add_trace(go.Scatter(\n",
+ " x=[0, x1[0]],\n",
+ " y=[0, y1[0]],\n",
+ " mode='lines+markers',\n",
+ " marker=dict(size=5),\n",
+ " line=dict(color='red', width=3),\n",
+ " name='Inhomogeneous Spin'),\n",
+ " row=1, col=1)\n",
+ "\n",
+ "# Add phase of the signal\n",
+ "fig.add_trace(go.Scatter(\n",
+ " x=[time[0]],\n",
+ " y=[angles[0]],\n",
+ " mode='markers',\n",
+ " marker=dict(color='blue', size=5),\n",
+ " name='Phase'),\n",
+ " row=1, col=2\n",
+ ")\n",
+ "fig.add_trace(go.Scatter(\n",
+ " x=[time],\n",
+ " y=[angles1],\n",
+ " mode='markers',\n",
+ " marker=dict(color='red', size=5),\n",
+ " name='Inhomogeneous Phase'),\n",
+ " row=1, col=2\n",
+ ")\n",
+ "# Rotating frame\n",
+ "fig.add_trace(go.Scatter(\n",
+ " x=[0, x_rot[0]],\n",
+ " y=[0, y_rot[0]],\n",
+ " mode='lines+markers',\n",
+ " marker=dict(size=5),\n",
+ " line=dict(color='blue', width=3),\n",
+ " name='Spin at f0', visible=False),\n",
+ " row=1, col=1)\n",
+ "fig.add_trace(go.Scatter(\n",
+ " x=[0, x1_rot[0]],\n",
+ " y=[0, y1_rot[0]],\n",
+ " mode='lines+markers',\n",
+ " marker=dict(size=5),\n",
+ " line=dict(color='red', width=3),\n",
+ " name='Inhomogeneous Spin', visible=False),\n",
+ " row=1, col=1)\n",
+ "\n",
+ "# Add phase of the signal\n",
+ "fig.add_trace(go.Scatter(\n",
+ " x=[time[0]],\n",
+ " y=[angles_rot[0]],\n",
+ " mode='markers',\n",
+ " marker=dict(color='blue', size=5),\n",
+ " name='Phase', visible=False),\n",
+ " row=1, col=2\n",
+ ")\n",
+ "fig.add_trace(go.Scatter(\n",
+ " x=[time],\n",
+ " y=[angles1_rot],\n",
+ " mode='markers',\n",
+ " marker=dict(color='red', size=5),\n",
+ " name='Inhomogeneous Phase', visible=False),\n",
+ " row=1, col=2\n",
+ ")\n",
+ "\n",
+ "fig.update_xaxes(range=[-1.1, 1.1], row=1, col=1)\n",
+ "fig.update_yaxes(range=[-1.1, 1.1], row=1, col=1)\n",
+ "fig.update_xaxes(range=[np.min(time), np.max(time)], row=1, col=2)\n",
+ "fig.update_yaxes(range=[np.min(angles) + 0.1*np.min(angles), np.max(angles) + 0.1*np.max(angles)], row=1, col=2)\n",
+ "\n",
+ "# Add frames\n",
+ "frames = [dict(\n",
+ " data=[go.Scatter(x=[0, x[i]],\n",
+ " y=[0, y[i]],\n",
+ " mode='lines+markers',\n",
+ " marker=dict(size=5),\n",
+ " line=dict(color='#636EFA', width=3),\n",
+ " name='Spin at f0'),\n",
+ " go.Scatter(x=[0, x1[i]],\n",
+ " y=[0, y1[i]],\n",
+ " mode='lines+markers',\n",
+ " marker=dict(size=5),\n",
+ " line=dict(color='#fa6363', width=3),\n",
+ " name='Inhomogeneous Spin'),\n",
+ " go.Scatter(x=time[:i],\n",
+ " y=angles[:i],\n",
+ " mode='lines',\n",
+ " marker=dict(size=5),\n",
+ " line=dict(color='blue', width=3),\n",
+ " name='Phase'),\n",
+ " go.Scatter(x=time[:i],\n",
+ " y=angles1[:i],\n",
+ " mode='lines',\n",
+ " marker=dict(size=5),\n",
+ " line=dict(color='red', width=3),\n",
+ " name='Inhomogeneous Phase'),\n",
+ " go.Scatter(x=[0, x_rot[i]],\n",
+ " y=[0, y_rot[i]],\n",
+ " mode='lines+markers',\n",
+ " marker=dict(size=5),\n",
+ " line=dict(color='#636EFA', width=3),\n",
+ " name='Spin at f0'),\n",
+ " go.Scatter(x=[0, x1_rot[i]],\n",
+ " y=[0, y1_rot[i]],\n",
+ " mode='lines+markers',\n",
+ " marker=dict(size=5),\n",
+ " line=dict(color='#fa6363', width=3),\n",
+ " name='Inhomogeneous Spin'),\n",
+ " go.Scatter(x=time[:i],\n",
+ " y=angles_rot[:i],\n",
+ " mode='lines',\n",
+ " marker=dict(size=5),\n",
+ " line=dict(color='blue', width=3),\n",
+ " name='Phase'),\n",
+ " go.Scatter(x=time[:i],\n",
+ " y=angles1_rot[:i],\n",
+ " mode='lines',\n",
+ " marker=dict(size=5),\n",
+ " line=dict(color='red', width=3),\n",
+ " name='Inhomogeneous Phase')\n",
+ " \n",
+ " \n",
+ " ],\n",
+ " name=str(i),\n",
+ " traces=[0,1,2,3,4,5,6,7]) for i in range(num_frames)]\n",
+ "\n",
+ "fig.frames = frames\n",
+ "\n",
+ "# Determine the maximum absolute value of coordinates\n",
+ "max_coord = max(abs(x.max()), abs(y.max()))\n",
+ "\n",
+ "fig.update_xaxes(title_text=\"x\", row=1, col=1)\n",
+ "fig.update_xaxes(title_text=\"time\", row=1, col=2)\n",
+ "fig.update_yaxes(title_text=\"y\", row=1, col=1)\n",
+ "fig.update_yaxes(title_text=\"rad\", tickmode = 'array',\n",
+ " tickvals = [-math.pi, 0, math.pi],\n",
+ " ticktext = [f\"-{PI_UNICODE}\", 0, f'{PI_UNICODE}'], row=1, col=2)\n",
+ "\n",
+ "\n",
+ "# Update layout\n",
+ "fig.update_layout(\n",
+ " height=600,\n",
+ " title=\"Spins rotating\",\n",
+ " xaxis=dict(autorange=False),\n",
+ " yaxis=dict(autorange=False),\n",
+ " updatemenus=\n",
+ " [dict(\n",
+ " type='buttons',\n",
+ " buttons=[dict(label='Play',\n",
+ " method='animate',\n",
+ " args=[None, dict(frame=dict(duration=50, redraw=False), fromcurrent=True, mode='immediate')]),\n",
+ " dict(label='Pause',\n",
+ " method='animate',\n",
+ " args=[[None], dict(frame=dict(duration=0, redraw=False), mode='immediate')])\n",
+ " ],\n",
+ " ),\n",
+ " dict(\n",
+ " buttons=[dict(\n",
+ " args=[{\"visible\": [True, True, True, True, False, False, False, False]}],\n",
+ " label=\"Laboratory Frame\",\n",
+ " method=\"update\"),\n",
+ " dict(\n",
+ " args=[{\"visible\": [False, False, False, False, True, True, True, True]}],\n",
+ " label=\"Rotating Frame\",\n",
+ " method=\"update\"\n",
+ " )],\n",
+ " direction=\"down\",\n",
+ " pad={\"b\": 70},\n",
+ " showactive=True,\n",
+ " x=-0.13,\n",
+ " xanchor=\"left\",\n",
+ " y=-0.15,\n",
+ " yanchor=\"top\"\n",
+ " )\n",
+ " ]\n",
+ ")\n",
+ "\n",
+ "# Show figure\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "867c0c25-97f5-4236-95e3-17982e3e5be8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fname_mag_e1 = os.path.join(\"fmap\", \"sub-fmap_magnitude1.nii.gz\")\n",
+ "fname_phase_e1 = os.path.join(\"fmap\", \"sub-fmap_phase1.nii.gz\")\n",
+ "fname_phase_e1_json = os.path.join(\"fmap\", \"sub-fmap_phase1.json\")\n",
+ "fname_phase_e2 = os.path.join(\"fmap\", \"sub-fmap_phase2.nii.gz\")\n",
+ "fname_mask = os.path.join(\"fmap\", \"mask.nii.gz\")\n",
+ "fname_fmap = os.path.join(\"fmap\", \"fmap.nii.gz\")\n",
+ "\n",
+ "nii_mag_e1 = nib.load(fname_mag_e1)\n",
+ "nii_phase_e1 = nib.load(fname_phase_e1)\n",
+ "nii_phase_e2 = nib.load(fname_phase_e2)\n",
+ "nii_mask = nib.load(fname_mask)\n",
+ "nii_fmap = nib.load(fname_fmap)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3827145b-5911-480c-ab37-1acdb99bb918",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Phase evolution though different echo times\n",
+ "mask = nii_mask.get_fdata()[30:-30,8:105,30]\n",
+ "fmap = nii_fmap.get_fdata()[30:-30,8:105,30] * mask # [Hz]\n",
+ "phase1 = (nii_phase_e1.get_fdata()[30:-30,8:105,30] / 4095 * 2 * math.pi - math.pi) * mask\n",
+ "\n",
+ "with open(fname_phase_e1_json, 'r') as json_data1:\n",
+ " data1 = json.load(json_data1)\n",
+ " \n",
+ "echo_time1 = data1['EchoTime']\n",
+ "phase0 = phase1 - (echo_time1 * (fmap * 2 * math.pi))\n",
+ "zmin = -math.pi\n",
+ "zmax = math.pi\n",
+ "\n",
+ "steps = 31\n",
+ "last_echo_time = 0.03\n",
+ "echo_times = np.linspace(0.0, last_echo_time, steps)\n",
+ "fig = go.Figure()\n",
+ "for i_echo, echo_time in enumerate(echo_times):\n",
+ " phase = phase0 + (fmap * echo_time * 2 * math.pi)\n",
+ " phase = np.angle(np.exp(1j*phase))\n",
+ " if i_echo >= len(echo_times) - 1:\n",
+ " fig.add_trace(go.Heatmap(z=np.rot90(phase, k=-1), visible=True, coloraxis = \"coloraxis\"))\n",
+ " else:\n",
+ " fig.add_trace(go.Heatmap(z=np.rot90(phase, k=-1), visible=False, coloraxis = \"coloraxis\"))\n",
+ "\n",
+ "fig.update_layout(coloraxis = {'colorscale':'gray'},\n",
+ " coloraxis_cmin=zmin, coloraxis_cmax=zmax)\n",
+ "fig.update_coloraxes(\n",
+ " colorbar=dict(title=\"Rad\",\n",
+ " titleside=\"top\",\n",
+ " tickmode=\"array\",\n",
+ " tickvals=[-math.pi, 0, math.pi-0.01],\n",
+ " ticktext = [f\"-{PI_UNICODE}\", 0, f'{PI_UNICODE}']))\n",
+ "\n",
+ "echo_times_str = [f\"{time:.2}\" for time in echo_times]\n",
+ "steps = []\n",
+ "for i in range(len(fig.data)):\n",
+ " step = dict(\n",
+ " method=\"update\",\n",
+ " label=echo_times_str[i],\n",
+ " args=[{\"visible\": [False] * len(fig.data)}], # layout attribute\n",
+ " )\n",
+ " step[\"args\"][0][\"visible\"][i] = True # Toggle i'th trace to \"visible\"\n",
+ " steps.append(step)\n",
+ "\n",
+ "sliders = [dict(\n",
+ " active=30,\n",
+ " currentvalue={\"prefix\": \"Echo Time: \"},\n",
+ " steps=steps\n",
+ ")]\n",
+ "\n",
+ "fig.update_layout(\n",
+ " sliders=sliders\n",
+ ")\n",
+ "\n",
+ "fig.update_layout(\n",
+ " title=dict(text=\"Phase at different echo times\", x=0.5)\n",
+ ")\n",
+ "fig.update_xaxes(showticklabels=False)\n",
+ "fig.update_yaxes(showticklabels=False)\n",
+ "fig.update_layout({\"height\": 550, \"width\": 500})\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5e51f4c7-1f99-4610-b610-c2eafe51d311",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def complex_difference(phase1, phase2):\n",
+ " \"\"\" Calculates the complex difference between 2 phase arrays (phase2 - phase1)\n",
+ "\n",
+ " Args:\n",
+ " phase1 (numpy.ndarray): Array containing phase data in radians\n",
+ " phase2 (numpy.ndarray): Array containing phase data in radians. Must be the same shape as phase1.\n",
+ "\n",
+ " Returns:\n",
+ " numpy.ndarray: The difference in phase between each voxels of phase2 and phase1 (phase2 - phase1)\n",
+ " \"\"\"\n",
+ "\n",
+ " # Calculate phasediff using complex difference\n",
+ " comp_0 = np.exp(-1j * phase1)\n",
+ " comp_1 = np.exp(1j * phase2)\n",
+ " return np.angle(comp_0 * comp_1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7d308a9a-1257-4535-b536-88646075b068",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def plot_2_echo_fmap(phase1, phase2, echotime1, echotime2):\n",
+ " phase_diff = complex_difference(phase1, phase2)\n",
+ " fmap = phase_diff / (echotime2 - echotime1) / 2 / math.pi\n",
+ " n=4\n",
+ " # Attempt at subplots\n",
+ " fig = make_subplots(rows=1, cols=n, shared_xaxes=False, horizontal_spacing=0.1, subplot_titles=(\"Phase 1\", \"Phase 2\", \"Phase difference\", \"B0 field map\"), specs=[[{\"type\": \"Heatmap\"}, {\"type\": \"Heatmap\"}, {\"type\": \"Heatmap\"}, {\"type\": \"Heatmap\"}]], )\n",
+ " \n",
+ " fig.add_trace(go.Heatmap(z=np.rot90(phase1, k=-1), colorscale='gray', colorbar_x=1/n - 0.05, zmin=-math.pi, zmax=math.pi,\n",
+ " colorbar=dict(title=\"Rad\",\n",
+ " titleside=\"top\",\n",
+ " tickmode=\"array\",\n",
+ " tickvals=[-math.pi, 0, math.pi-0.01],\n",
+ " ticktext = [f\"-{PI_UNICODE}\", 0, f'{PI_UNICODE}'])), 1, 1)\n",
+ " fig.add_trace(go.Heatmap(z=np.rot90(phase2, k=-1), colorscale='gray', colorbar_x=2/n - 0.02, zmin=-math.pi, zmax=math.pi,\n",
+ " colorbar=dict(title=\"Rad\",\n",
+ " titleside=\"top\",\n",
+ " tickmode=\"array\",\n",
+ " tickvals=[-math.pi, 0, math.pi-0.01],\n",
+ " ticktext = [f\"-{PI_UNICODE}\", 0, f'{PI_UNICODE}'])), 1, 2)\n",
+ " fig.add_trace(go.Heatmap(z=np.rot90(phase_diff, k=-1), colorscale='gray', colorbar_x=3/n - 0.02, zmin=-math.pi, zmax=math.pi,\n",
+ " colorbar=dict(title=\"Rad\",\n",
+ " titleside=\"top\",\n",
+ " tickmode=\"array\",\n",
+ " tickvals=[-math.pi, 0, math.pi-0.01],\n",
+ " ticktext = [f\"-{PI_UNICODE}\", 0, f'{PI_UNICODE}'])), 1, 3)\n",
+ " fig.add_trace(go.Heatmap(z=np.rot90(fmap, k=-1), colorscale='gray',\n",
+ " colorbar=dict(title=\"Hz\",\n",
+ " titleside=\"top\")), 1, 4)\n",
+ " fig.update_xaxes(showticklabels=False)\n",
+ " fig.update_yaxes(showticklabels=False)\n",
+ " fig.update_layout({\"height\": 450, \"width\": 1700})\n",
+ " \n",
+ " fig.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "909bbd00-9f7e-41f0-b59d-6ce94f02834c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def get_circle(x, y, r):\n",
+ " if x < 1 or y < 1 or r < 1:\n",
+ " raise ValueError(\"Input parameters are too small\")\n",
+ " \n",
+ " my_array = np.zeros([x,y])\n",
+ " for i in range(x):\n",
+ " for j in range(y):\n",
+ " squared = (i-(x/2))**2 + (j-(y/2))**2\n",
+ " h = np.sqrt(squared)\n",
+ " if h < r:\n",
+ " my_array[i,j] = 1\n",
+ " return my_array\n",
+ "\n",
+ "echo1 = get_circle(100, 100, 30) * -1\n",
+ "echo2 = get_circle(100, 100, 30) * 2\n",
+ "echo_time1 = 0.005\n",
+ "echo_time2 = 0.01\n",
+ "\n",
+ "plot_2_echo_fmap(echo1, echo2, echo_time1, echo_time2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "34fffc05-5890-48dc-b13f-3fa00a3e5c49",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "mask = nii_mask.get_fdata()[30:-30,8:105,30]\n",
+ "phase1 = (nii_phase_e1.get_fdata()[30:-30,8:105,30] / 4095 * 2 * math.pi - math.pi) * mask\n",
+ "phase2 = (nii_phase_e2.get_fdata()[30:-30,8:105,30] / 4095 * 2 * math.pi - math.pi) * mask\n",
+ "echo_time1 = 0.00338\n",
+ "echo_time2 = 0.00558\n",
+ "\n",
+ "plot_2_echo_fmap(phase1, phase2, echo_time1, echo_time2)"
+ ]
+ }
+ ],
+ "metadata": {
+ "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.0"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/05-B0/content/chapter4_3.ipynb b/05-B0/content/chapter4_3.ipynb
new file mode 100644
index 0000000..a1c6d5e
--- /dev/null
+++ b/05-B0/content/chapter4_3.ipynb
@@ -0,0 +1,291 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7ee2da8d-6336-4054-99c7-bc374f0a0e7b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import plotly.express as px\n",
+ "import os\n",
+ "import nibabel as nib\n",
+ "import numpy as np\n",
+ "from plotly.subplots import make_subplots\n",
+ "import plotly.graph_objects as go\n",
+ "import math\n",
+ "from sklearn.linear_model import LinearRegression\n",
+ "from dash import Dash, dcc, html\n",
+ "import pandas as pd\n",
+ "PI_UNICODE = \"\\U0001D70B\"\n",
+ "GYRO_BAR_RATIO_H = 42.6e6 # [Hz/T]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "35eff464-d7ad-4124-b892-a57b6c4d410a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Taken from real brain data\n",
+ "phase1, phase2, phase3, phase4 = (2.9298516102709202, -0.21864564255753116, -2.2884910587688285, 2.392827225041896)\n",
+ "beg = 0\n",
+ "end = 0.015\n",
+ "\n",
+ "t = np.array([0.00263, 0.00526, 0.009, 0.013])\n",
+ "n_echoes = len(t)\n",
+ "y = np.array([phase1, phase2, phase3, phase4])\n",
+ "y_unwrapped = np.array([phase1, phase2, phase3, phase4 - (2 * math.pi)])\n",
+ "\n",
+ "reg = LinearRegression().fit(t.reshape(-1, 1), y_unwrapped.reshape(-1,1))\n",
+ "# Slope of linear regression reshaped into the shape of original 3D phase.\n",
+ "fieldmap_rad = reg.coef_[0] # [rad / s]\n",
+ "fieldmap_intercept = reg.intercept_[0] # [rad / s]\n",
+ "\n",
+ "t_predict = np.array([beg, end])\n",
+ "y_predict = reg.predict(t_predict.reshape(-1,1))[:,0]\n",
+ "\n",
+ "\n",
+ "fig = go.Figure()\n",
+ "fig.add_trace(go.Scatter(x=t, y=y, mode='markers', marker=dict(color='blue'), name='Original Phase'))\n",
+ "fig.add_trace(go.Scatter(x=[t[3]], y=[phase4- (2*math.pi)], mode='markers', marker=dict(color='red'), name='Unwrapped Phase'))\n",
+ "fig.add_trace(go.Scatter(x=t_predict, y=y_predict, mode='lines', marker=dict(color='green'), name='Fit'))\n",
+ "fig.add_trace(go.Scatter(x=[beg, end], y=[math.pi, math.pi], mode='lines', line=dict(color='gray'), showlegend=False))\n",
+ "fig.add_trace(go.Scatter(x=[beg, end], y=[-math.pi, -math.pi], mode='lines', line=dict(color='gray'), showlegend=False))\n",
+ "\n",
+ "fig.update_xaxes(title_text=\"Time (ms)\", range=[beg, end])\n",
+ "fig.update_yaxes(title_text=\"Phase (rad)\", tickmode = 'array', range=[-7,7],\n",
+ " tickvals = [-2*math.pi, -math.pi, 0, math.pi, 2*math.pi],\n",
+ " ticktext = [f'-2{PI_UNICODE}', f'-{PI_UNICODE}', '0', f'{PI_UNICODE}', f'2{PI_UNICODE}'])\n",
+ "fig.update_layout({\"width\": 800})\n",
+ "\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "12c5e585-9ce1-4e46-afea-6f5960f09d31",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "t = np.linspace(0, 10, 1001)\n",
+ "a = 2\n",
+ "y = a * t - 10\n",
+ "\n",
+ "# \n",
+ "phase = np.angle(np.exp(1j*y))\n",
+ "\n",
+ "fig = go.Figure()\n",
+ "fig.add_trace(go.Scatter(x=t, y=phase, mode='lines', name='Wrapped phase'))\n",
+ "fig.add_scatter(x=t, y=y, mode='lines', name='True phase')\n",
+ "fig.update_traces(marker=dict(size=3))\n",
+ "fig.update_xaxes(title_text=\"x\")\n",
+ "fig.update_yaxes(title_text=\"rad\", tickmode = 'array',\n",
+ " tickvals = [-3*math.pi, -2*math.pi, -math.pi, 0, math.pi, 2*math.pi, 3*math.pi],\n",
+ " ticktext = [f'-3{PI_UNICODE}', f'-2{PI_UNICODE}', f'-{PI_UNICODE}', '0', f'{PI_UNICODE}', f'2{PI_UNICODE}', f'3{PI_UNICODE}'])\n",
+ "fig.update_layout({\"width\": 800})\n",
+ "fig.show()\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a3750f9e-0258-4118-810c-fbb96d376b4d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import dash\n",
+ "import dash_bootstrap_components as dbc\n",
+ "from dash import Dash, html, dcc\n",
+ "import plotly.express as px\n",
+ "\n",
+ "t = np.linspace(0, 3, 1001)\n",
+ "y_unwrapped = t + t**2 - 3 * t**3 + t**4\n",
+ "y_wrapped = np.mod(y_unwrapped, 2 * math.pi);\n",
+ "\n",
+ "y_delta = y_unwrapped.max() - y_unwrapped.min()\n",
+ "min_y = y_unwrapped.min() - 2 * math.pi - y_delta * 0.05\n",
+ "max_y = y_unwrapped.max() + 2 * math.pi + y_delta * 0.05\n",
+ "\n",
+ "fig1 = go.Figure()\n",
+ "fig1.add_trace(go.Scatter(x=t, y=y_wrapped,\n",
+ " mode='lines',\n",
+ " name='Wrapped '))\n",
+ "fig1.update_layout(title_text=\"Wrapped\", title_x=0.25, showlegend=True,\n",
+ " legend={\"x\": 0.03, \"y\": 0.95})\n",
+ "fig1.update_yaxes(range=[min_y, max_y], title_text=\"rad\", tickmode = 'array',\n",
+ " tickvals = [-2*math.pi, -math.pi, 0, math.pi, 2*math.pi, 3*math.pi, 4*math.pi, 5*math.pi],\n",
+ " ticktext = [f'-2{PI_UNICODE}', f'-{PI_UNICODE}', '0', f'{PI_UNICODE}', f'2{PI_UNICODE}', f'3{PI_UNICODE}', f'4{PI_UNICODE}', f'5{PI_UNICODE}'])\n",
+ "fig1.update_xaxes(fixedrange=True)\n",
+ "fig1.update_yaxes(fixedrange=True)\n",
+ "\n",
+ "fig2 = go.Figure()\n",
+ "fig2.add_trace(go.Scatter(x=t, y=y_unwrapped,\n",
+ " mode='lines',\n",
+ " name='Solution 1', showlegend=True))\n",
+ "fig2.add_trace(go.Scatter(x=t, y=y_unwrapped+2*math.pi,\n",
+ " mode='lines',\n",
+ " name='Solution 2', showlegend=True))\n",
+ "fig2.add_trace(go.Scatter(x=t, y=y_unwrapped-2*math.pi,\n",
+ " mode='lines',\n",
+ " name='Solution 3', showlegend=True))\n",
+ "fig2.update_layout(title_text=\"Unwrapped\", title_x=0.75,\n",
+ " legend={\"x\": 0.85, \"y\": 0.95})\n",
+ "fig2.update_yaxes(range=[min_y, max_y], title_text=\"rad\", tickmode = 'array',\n",
+ " tickvals = [-2*math.pi, -math.pi, 0, math.pi, 2*math.pi, 3*math.pi, 4*math.pi, 5*math.pi],\n",
+ " ticktext = [f'-2{PI_UNICODE}', f'-{PI_UNICODE}', '0', f'{PI_UNICODE}', f'2{PI_UNICODE}', f'3{PI_UNICODE}', f'4{PI_UNICODE}', f'5{PI_UNICODE}'])\n",
+ "fig2.update_xaxes(fixedrange=True)\n",
+ "fig2.update_yaxes(fixedrange=True)\n",
+ "app = Dash(\n",
+ " __name__,\n",
+ " external_stylesheets=[dbc.themes.BOOTSTRAP, dbc.icons.FONT_AWESOME],\n",
+ " external_scripts=[{'src':\"https://ajax.googleapis.com/ajax/libs/jquery/3.6.0/jquery.min.js\"}]\n",
+ ")\n",
+ "\n",
+ "def beforeAfterSlide(fig1, fig2, style=None):\n",
+ " bfA = []\n",
+ " if not style:\n",
+ " style = {'width':'100vw', 'height':'100vh'}\n",
+ " for key in style:\n",
+ " if '%' in style[key]:\n",
+ " if key in ['width', 'left']:\n",
+ " style[key] = style[key].replace('%','vw')\n",
+ " if key in ['top', 'height']:\n",
+ " style[key] = style[key].replace('%','vh')\n",
+ " bfA.append(html.Div(dcc.Graph(figure=fig2, style=style), className='after'))\n",
+ " bfA.append(html.Div(className='middle'))\n",
+ " bfA.append(html.Div(dcc.Graph(figure=fig1, style=style), className='before'))\n",
+ " return html.Div(bfA, className='beforeAfter', style=style)\n",
+ "\n",
+ "app.layout = html.Div(beforeAfterSlide(fig1, fig2, {'height':'75%', 'width':'75%', 'top':'10%', 'left':'6%'}))\n",
+ "\n",
+ "if __name__ == \"__main__\":\n",
+ " app.run_server(debug=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6e5ff4f8-d6aa-46f0-beba-ac9e58e225ed",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "seed = 2\n",
+ "np.random.seed(seed)\n",
+ "n_points = 20\n",
+ "len_x = 500\n",
+ "len_y = 500\n",
+ "map = np.zeros([len_x, len_y]) - math.pi\n",
+ "x, y = np.meshgrid(range(len_x), range(len_y))\n",
+ "\n",
+ "def gauss(a, std_x, std_y, center):\n",
+ " x0 = int(center[0] * len_x)\n",
+ " y0 = int(center[1] * len_y)\n",
+ " return a * np.exp(-((((x-x0)**2)/(std_x**2)) + (((y-y0)**2)/(std_y**2))))\n",
+ "\n",
+ "\n",
+ "for i_point in range(n_points):\n",
+ " point = (np.random.uniform(), np.random.uniform())\n",
+ " amp = np.random.randint(4, 9)\n",
+ " std_x = np.random.randint(50, 100)\n",
+ " std_y = np.random.randint(50, 100)\n",
+ " map += gauss(amp, std_x, std_y, point)\n",
+ "\n",
+ "wrapped = np.angle(np.exp(1j*map))\n",
+ "\n",
+ "fig = make_subplots(rows=1, cols=2, shared_xaxes=False, horizontal_spacing=0.1, subplot_titles=(\"Wrapped\", \"Unwrapped\"), specs=[[{\"type\": \"Heatmap\"}, {\"type\": \"Heatmap\"}]])\n",
+ "fig.add_trace(go.Heatmap(z=wrapped, colorscale='gray', colorbar_x=0.47, colorbar=dict(title=\"Rad\", titleside=\"top\", tickmode=\"array\", tickvals=[-math.pi, 0, math.pi], ticktext = [f'-{PI_UNICODE}', 0, f'{PI_UNICODE}'])), 1, 1)\n",
+ "fig.add_trace(go.Heatmap(z=map, colorscale='gray', colorbar=dict(title=\"Rad\", titleside=\"top\", tickmode=\"array\",\n",
+ " tickvals=[-2*math.pi, 0, 2*math.pi, 4*math.pi, 6*math.pi, 8*math.pi, 10*math.pi, 12*math.pi, 14*math.pi],\n",
+ " ticktext = [f'-2{PI_UNICODE}', 0, f'2{PI_UNICODE}', f'4{PI_UNICODE}', f'6{PI_UNICODE}', f'8{PI_UNICODE}', f'10{PI_UNICODE}', f'12{PI_UNICODE}', f'14{PI_UNICODE}'])), 1, 2)\n",
+ "fig.update_layout({\"height\": 500, \"width\": 1000})\n",
+ "fig.update_xaxes(showticklabels=False)\n",
+ "fig.update_yaxes(showticklabels=False)\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "277d6bb9-e3b3-448d-bc50-99f834fb83c9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "len_x = 500\n",
+ "len_y = 500\n",
+ "x, y = np.meshgrid(np.linspace(-len_x, len_x/2, len_x//2), np.linspace(-len_y, len_y, len_y), indexing='ij')\n",
+ "x1 = np.array(x)\n",
+ "y1 = np.array(y)\n",
+ "x, y = np.meshgrid(np.linspace(-len_x/2, len_x, len_x//2), np.linspace(-len_y, len_y, len_y), indexing='ij')\n",
+ "x2 = np.array(x)\n",
+ "y2 = np.array(y)\n",
+ "x = np.zeros([len_x, len_y])\n",
+ "y = np.zeros([len_x, len_y])\n",
+ "x[:len_x//2] = x1\n",
+ "y[:len_x//2,:] = -y1\n",
+ "x[len_x//2:,:] = -x2\n",
+ "y[len_x//2:,:] = -y2\n",
+ "\n",
+ "\n",
+ "phase = np.rot90(np.arctan2(y.astype(float),x.astype(float)),k=-1)\n",
+ "fig = make_subplots(rows=1, cols=1, subplot_titles=[(\"Phase Singularities\")], specs=[[{\"type\": \"Heatmap\"}]])\n",
+ "fig.add_trace(go.Heatmap(z=phase, colorscale='gray', colorbar=dict(title=\"Rad\", titleside=\"top\", tickmode=\"array\", tickvals=[-math.pi+0.01, 0, math.pi-0.01], ticktext = [f'-{PI_UNICODE}', 0, f'{PI_UNICODE}'])))\n",
+ "fig.add_trace(go.Scatter(x=[370, 370], y=[90, 410], mode='markers', marker=dict(size=10), line=dict(color='red')))\n",
+ "\n",
+ "fig.add_shape(type=\"line\", x0=370, y0=90, x1=370, y1=410, line=dict(color=\"red\",width=3))\n",
+ "fig.add_shape(type=\"path\", path=\"M 370,90 Q 200,250 370,410\", line=dict(color=\"red\",width=3))\n",
+ "fig.add_annotation(x=390, y=90, text=\"A\", showarrow=False)\n",
+ "fig.add_annotation(x=390, y=410, text=\"B\", showarrow=False, font=dict(color=\"white\"))\n",
+ "fig.update_layout({\"width\": 500, \"height\": 400})\n",
+ "fig.update_xaxes(showticklabels=False)\n",
+ "fig.update_yaxes(showticklabels=False)\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "edb89952-5341-4c4a-b8ca-6f5844a18141",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "b0 = np.array((0.064, 0.1, 0.3, 1.5, 3, 7)) # [T]\n",
+ "f0 = b0 * GYRO_BAR_RATIO_H # [Hz]\n",
+ "df = 1e-6 * f0\n",
+ "te = 1e3/(2*df)\n",
+ "\n",
+ "tmp = {\n",
+ " 'B0': list(b0),\n",
+ " 'TE (ms)': list(te)\n",
+ "}\n",
+ "df = pd.DataFrame(tmp)\n",
+ "# df.to_markdown()\n",
+ "df = df.T.style.hide(axis='columns')"
+ ]
+ }
+ ],
+ "metadata": {
+ "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.0"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/05-B0/content/chapter4_4.ipynb b/05-B0/content/chapter4_4.ipynb
new file mode 100644
index 0000000..ce4f18c
--- /dev/null
+++ b/05-B0/content/chapter4_4.ipynb
@@ -0,0 +1,404 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "925fbe02-42e5-4449-9411-6c21926305ca",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import plotly.express as px\n",
+ "import os\n",
+ "import nibabel as nib\n",
+ "import numpy as np\n",
+ "from plotly.subplots import make_subplots\n",
+ "import plotly.graph_objects as go\n",
+ "import math\n",
+ "from sklearn.linear_model import LinearRegression\n",
+ "import copy\n",
+ "PI_UNICODE = \"\\U0001D70B\"\n",
+ "DELTA_UNICODE = \"\\u0394\"\n",
+ "GYRO_UNICODE = \"\\U0001D6FE\"\n",
+ "GREEK_DELTA_UNICODE = \"\\u03B4\""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e0f59d72-c91a-454d-b99b-7032358665e5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "phase1, phase2, phase3 = (-4, -2, 0)\n",
+ "beg = 0\n",
+ "end = 0.01\n",
+ "t = np.array([0.00263, 0.00526, 0.009])\n",
+ "n_echoes = len(t)\n",
+ "y1 = np.array([phase1, phase2, phase3])\n",
+ "y2 = y1 + 2 * math.pi\n",
+ "\n",
+ "# Linear fit of the first line\n",
+ "reg1 = LinearRegression().fit(t.reshape(-1, 1), y1.reshape(-1,1))\n",
+ "fieldmap_rad1 = reg1.coef_[0] # [rad / s]\n",
+ "fieldmap_intercept1 = reg1.intercept_[0] # [rad / s]\n",
+ "t_predict1 = np.array([beg, end])\n",
+ "y_predict1 = reg1.predict(t_predict1.reshape(-1,1))[:,0]\n",
+ "\n",
+ "# Linear fit of the second line\n",
+ "reg2 = LinearRegression().fit(t.reshape(-1, 1), y2.reshape(-1,1))\n",
+ "fieldmap_rad2 = reg2.coef_[0] # [rad / s]\n",
+ "fieldmap_intercept2 = reg2.intercept_[0] # [rad / s]\n",
+ "t_predict2 = np.array([beg, end])\n",
+ "y_predict2 = reg2.predict(t_predict2.reshape(-1,1))[:,0]\n",
+ "\n",
+ "# Plot\n",
+ "fig = go.Figure()\n",
+ "fig.add_trace(go.Scatter(x=t, y=y1, mode='markers', marker=dict(color='blue'), name='Unwrapped solution 1'))\n",
+ "fig.add_trace(go.Scatter(x=t_predict1, y=y_predict1, mode='lines', marker=dict(color='blue'), name='Fit'))\n",
+ "fig.add_trace(go.Scatter(x=t, y=y2, mode='markers', marker=dict(color='red'), name='Unwrapped solution 2'))\n",
+ "fig.add_trace(go.Scatter(x=t_predict2, y=y_predict2, mode='lines', marker=dict(color='red'), name='Fit'))\n",
+ "fig.add_annotation(x=0.004, y=-3, text=f\"Slope: {GYRO_UNICODE}{DELTA_UNICODE}B*t\")\n",
+ "fig.add_annotation(x=0.004, y=3.3, text=f\"Slope: {GYRO_UNICODE}{DELTA_UNICODE}B*t\")\n",
+ "fig.update_xaxes(title_text=\"Time (ms)\", range=[beg, end])\n",
+ "fig.update_yaxes(title_text=\"Phase (rad)\", tickmode = 'array', range=[-8,8],\n",
+ " tickvals = [-2*math.pi, -math.pi, 0, math.pi, 2*math.pi],\n",
+ " ticktext = [f'-2{PI_UNICODE}', f'-{PI_UNICODE}', '0', f'{PI_UNICODE}', f'2{PI_UNICODE}'])\n",
+ "fig.update_layout({\"width\": 800})\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "09b870d2-8be9-4c7c-8e26-64292afc7874",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "phase1, phase2, phase3 = (-1, 2, 14)\n",
+ "phase_wrapped1 = 14-4*math.pi\n",
+ "phase_unwrapped1 = 14-2*math.pi\n",
+ "beg = 0\n",
+ "end = 0.016\n",
+ "\n",
+ "t1 = np.array([0.00263, 0.00526])\n",
+ "t2 = np.array([0.00263, 0.00526, 0.015])\n",
+ "n_echoes1 = len(t1)\n",
+ "n_echoes2 = len(t2)\n",
+ "y1 = np.array([phase1, phase2])\n",
+ "y2 = np.array([phase1, phase2, phase3])\n",
+ "\n",
+ "reg1 = LinearRegression().fit(t1.reshape(-1, 1), y1.reshape(-1,1))\n",
+ "# Slope of linear regression reshaped into the shape of original 3D phase.\n",
+ "fieldmap_rad1 = reg1.coef_[0] # [rad / s]\n",
+ "fieldmap_intercept1 = reg1.intercept_[0] # [rad / s]\n",
+ "t_predict1 = np.array([beg, end])\n",
+ "y_predict1 = reg1.predict(t_predict1.reshape(-1,1))[:,0]\n",
+ "\n",
+ "reg2 = LinearRegression().fit(t2.reshape(-1, 1), y2.reshape(-1,1))\n",
+ "# Slope of linear regression reshaped into the shape of original 3D phase.\n",
+ "fieldmap_rad2 = reg2.coef_[0] # [rad / s]\n",
+ "fieldmap_intercept2 = reg2.intercept_[0] # [rad / s]\n",
+ "t_predict2 = np.array([beg, end])\n",
+ "y_predict2 = reg2.predict(t_predict2.reshape(-1,1))[:,0]\n",
+ "\n",
+ "fig = go.Figure()\n",
+ "fig.add_trace(go.Scatter(x=t2, y=y2, mode='markers', marker=dict(color='green'), name='Unwrapped phase'))\n",
+ "fig.add_trace(go.Scatter(x=[t2[2]], y=[phase_wrapped1], mode='markers', marker=dict(color='blue'), name='Acquired wrapped phase'))\n",
+ "fig.add_trace(go.Scatter(x=[t2[2]], y=[phase_unwrapped1], mode='markers', marker=dict(color='orange'), name='Wrong unwrapped phase'))\n",
+ "fig.add_trace(go.Scatter(x=t_predict1, y=y_predict1, mode='lines', marker=dict(color='blue'), name='Fit with first 2 echoes'))\n",
+ "fig.add_trace(go.Scatter(x=t_predict2, y=y_predict2, mode='lines', marker=dict(color='green'), name='Fit with 3 echoes'))\n",
+ "fig.add_trace(go.Scatter(x=[beg, end], y=[math.pi, math.pi], mode='lines', line=dict(color='gray'), showlegend=False))\n",
+ "fig.add_trace(go.Scatter(x=[beg, end], y=[-math.pi, -math.pi], mode='lines', line=dict(color='gray'), showlegend=False))\n",
+ "fig.update_xaxes(title_text=\"Time (ms)\", range=[beg, end])\n",
+ "fig.update_yaxes(title_text=\"Phase (rad)\", tickmode = 'array', range=[-5,15],\n",
+ " tickvals = [-2*math.pi, -math.pi, 0, math.pi, 2*math.pi, 3*math.pi, 4*math.pi, 5*math.pi],\n",
+ " ticktext = [f'-2{PI_UNICODE}', f'-{PI_UNICODE}', '0', f'{PI_UNICODE}', f'2{PI_UNICODE}', f'3{PI_UNICODE}', f'4{PI_UNICODE}', f'5{PI_UNICODE}'])\n",
+ "fig.update_layout({\"width\": 800}, title_text=f\"Using fast {DELTA_UNICODE}TE to help temporally unwrapped 3rd echo\", title_x=0.5)\n",
+ "\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6a3b1dbd-4e48-4d61-b528-ca612768cafc",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def complex_difference(phase1, phase2):\n",
+ " \"\"\" Calculates the complex difference between 2 phase arrays (phase2 - phase1)\n",
+ "\n",
+ " Args:\n",
+ " phase1 (numpy.ndarray): Array containing phase data in radians\n",
+ " phase2 (numpy.ndarray): Array containing phase data in radians. Must be the same shape as phase1.\n",
+ "\n",
+ " Returns:\n",
+ " numpy.ndarray: The difference in phase between each voxels of phase2 and phase1 (phase2 - phase1)\n",
+ " \"\"\"\n",
+ "\n",
+ " # Calculate phasediff using complex difference\n",
+ " comp_0 = np.exp(-1j * phase1)\n",
+ " comp_1 = np.exp(1j * phase2)\n",
+ " return np.angle(comp_0 * comp_1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e138a1ad-3945-4f6e-a51c-dd13c9bf6309",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def umpire_3echoes(phases, times):\n",
+ " \"\"\"\n",
+ " This function performs unwrapping using the UMPIRE algorithm with 3 echoes. UMPIRE requires echo times that are unevenly spaced.\n",
+ " \"\"\"\n",
+ " \n",
+ " # Complex difference\n",
+ " dpTE2 = complex_difference(phases[1], phases[2])\n",
+ " dpTE1 = complex_difference(phases[0], phases[1])\n",
+ " dpd = complex_difference(dpTE1, dpTE2)\n",
+ " # print(\"Diff in phase diff:\" , dpd)\n",
+ " dTEs = np.array([times[1]-times[0], times[2]-times[1]])\n",
+ " dt_dpd = dTEs[1] - dTEs[0]\n",
+ " \n",
+ " # Slope\n",
+ " slope = dpd / dt_dpd\n",
+ " \n",
+ " # n wraps in differences\n",
+ " n_wraps_dp = np.round((dTEs - dTEs*slope) / (2*math.pi))\n",
+ " \n",
+ " # Remove wraps in differences\n",
+ " dpTE1_prime = dpTE1 - (2*n_wraps_dp[0]*math.pi)\n",
+ " dpTE2_prime = dpTE2 - (2*n_wraps_dp[1]*math.pi)\n",
+ " \n",
+ " # Calculate better slope\n",
+ " slope_prime1 = dpTE1_prime / dTEs[0]\n",
+ " slope_prime2 = dpTE2_prime / dTEs[1]\n",
+ " slope_avg = (slope_prime1 + slope_prime2) / 2\n",
+ " \n",
+ " # Calculate wraps in original phase\n",
+ " n_wraps = np.round((phases - t*slope_avg) / (2*math.pi))\n",
+ " \n",
+ " # Remove wraps\n",
+ " unwrapped_with_phase_offset = phases - 2*math.pi*n_wraps\n",
+ " \n",
+ " # # Calculate receiver offset\n",
+ " # r = (t[0] * unwrapped_with_phase_offset[1] - t[1] * unwrapped_with_phase_offset[0]) / dTEs[0]\n",
+ "\n",
+ " # # Remove receiver phase offset\n",
+ " # phase_no_offset = complex_difference(r, unwrapped_with_phase_offset)\n",
+ " # # Unwrap one last time\n",
+ " # ns = np.round((phase_no_offset - t*slope_avg) / (2*math.pi))\n",
+ " # unwrapped_umpire = phase_no_offset - 2*math.pi*ns\n",
+ " \n",
+ " return unwrapped_with_phase_offset"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "23f65935-9811-445a-a2b7-f00fca06ed02",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "t = np.array([0.003, 0.011, 0.020])\n",
+ "y_unwrapped = np.array([1.0, 9.05, 17.75])\n",
+ "wrapped = copy.deepcopy(y_unwrapped)\n",
+ "wrapped[0] = np.angle(np.exp(1j*wrapped[0]))\n",
+ "wrapped[1] = np.angle(np.exp(1j*wrapped[1]))\n",
+ "wrapped[2] = np.angle(np.exp(1j*wrapped[2]))\n",
+ "beg = 0.0\n",
+ "end = 0.021\n",
+ "\n",
+ "# Fit original data\n",
+ "reg1 = LinearRegression().fit(t.reshape(-1, 1), y_unwrapped.reshape(-1,1))\n",
+ "fieldmap_rad1 = reg1.coef_[0] # [rad / s]\n",
+ "fieldmap_intercept1 = reg1.intercept_[0] # [rad / s]\n",
+ "t_predict1 = np.array([beg, end])\n",
+ "y_predict1 = reg1.predict(t_predict1.reshape(-1,1))[:,0]\n",
+ "\n",
+ "# Unwrap with UMPIRE\n",
+ "unwrapped_umpire = umpire_3echoes(wrapped, t)\n",
+ "\n",
+ "# Fit unwrapped data of UMPIRE\n",
+ "reg2 = LinearRegression().fit(t.reshape(-1, 1), unwrapped_umpire.reshape(-1,1))\n",
+ "# Slope of linear regression reshaped into the shape of original 3D phase.\n",
+ "fieldmap_rad2 = reg2.coef_[0] # [rad / s]\n",
+ "fieldmap_intercept2 = reg2.intercept_[0] # [rad / s]\n",
+ "t_predict2 = np.array([beg, end])\n",
+ "y_predict2 = reg2.predict(t_predict2.reshape(-1,1))[:,0]\n",
+ "\n",
+ "# Plot\n",
+ "height_annotations = 7*math.pi + 0.5\n",
+ "fig = go.Figure()\n",
+ "fig.add_trace(go.Scatter(x=t, y=wrapped, mode='markers', marker=dict(color='blue'), name='Wrapped'))\n",
+ "fig.add_trace(go.Scatter(x=t, y=y_unwrapped, mode='markers', marker=dict(color='red'), name='True phase'))\n",
+ "fig.add_trace(go.Scatter(x=t, y=unwrapped_umpire, mode='markers', marker=dict(color='green'), name='Umpire', visible='legendonly'))\n",
+ "fig.add_trace(go.Scatter(x=t_predict2, y=y_predict2, mode='lines', marker=dict(color='green'), name='Umpire fit', visible='legendonly'))\n",
+ "fig.add_trace(go.Scatter(x=[beg, end], y=[math.pi, math.pi], mode='lines', line=dict(color='gray'), showlegend=False))\n",
+ "fig.add_trace(go.Scatter(x=[beg, end], y=[-math.pi, -math.pi], mode='lines', line=dict(color='gray'), showlegend=False))\n",
+ "fig.add_trace(go.Scatter(x=[t[0], t[1]], y=[height_annotations, height_annotations], mode='lines+markers', \n",
+ " marker=dict(symbol=\"line-ns-open\", color=\"black\",size=10),\n",
+ " line=dict(color='black'), showlegend=False))\n",
+ "fig.add_trace(go.Scatter(x=[t[1], t[2]], y=[height_annotations, height_annotations], mode='lines+markers', \n",
+ " marker=dict(symbol=\"line-ns-open\", color=\"black\",size=10),\n",
+ " line=dict(color='black'), showlegend=False))\n",
+ "fig.add_trace(go.Scatter(x=[t[1], t[2]], y=[height_annotations-3, height_annotations-3], mode='lines+markers', \n",
+ " marker=dict(symbol=\"line-ns-open\", color=\"black\",size=10),\n",
+ " line=dict(color='black'), showlegend=False))\n",
+ "fig.add_trace(go.Scatter(x=[0.002, 0.002], y=[y_unwrapped[0], y_unwrapped[1]], mode='lines+markers', \n",
+ " marker=dict(symbol=\"line-ew-open\", color=\"black\",size=10),\n",
+ " line=dict(color='black'), showlegend=False))\n",
+ "fig.add_annotation(x=(t[1]-t[0])/2 + t[0], y=height_annotations+1.15, text=f\"{DELTA_UNICODE}TE1\", showarrow=False)\n",
+ "fig.add_annotation(x=(t[2]-t[1])/2 + t[1], y=height_annotations+1.15, text=f\"{DELTA_UNICODE}TE2\", showarrow=False)\n",
+ "fig.add_annotation(x=(t[2]-t[1])/2 + t[1], y=height_annotations+1.15-3, text=f\"{DELTA_UNICODE}TE1 + {GREEK_DELTA_UNICODE}TE\", showarrow=False)\n",
+ "fig.add_annotation(x=0.0015, y=(y_unwrapped[1] + y_unwrapped[0]) / 2, text=f\">2{PI_UNICODE}\", showarrow=False)\n",
+ "fig.update_xaxes(title_text=\"Time (ms)\", range=[beg, end])\n",
+ "fig.update_yaxes(title_text=\"Phase (rad)\", tickmode = 'array', range=[-4,25],\n",
+ " tickvals = [-2*math.pi, 0, 2*math.pi, 4*math.pi, 6*math.pi],\n",
+ " ticktext = [f'-2{PI_UNICODE}', '0', f'2{PI_UNICODE}', f'4{PI_UNICODE}', f'6{PI_UNICODE}'])\n",
+ "fig.update_layout({\"width\": 800}, title_text=\"Phase unwrapping using UMPIRE algorithm\", title_x=0.5)\n",
+ "fig.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "122984e1-72ea-43f1-af05-d3a3f965179d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "t = np.array([0.003, 0.011, 0.020])\n",
+ "y_unwrapped = np.array([1.0, 9.05, 17.75])\n",
+ "wrapped = copy.deepcopy(y_unwrapped)\n",
+ "wrapped[0] = np.angle(np.exp(1j*wrapped[0]))\n",
+ "wrapped[1] = np.angle(np.exp(1j*wrapped[1]))\n",
+ "wrapped[2] = np.angle(np.exp(1j*wrapped[2]))\n",
+ "beg = 0.0\n",
+ "end = 0.021\n",
+ "\n",
+ "fig = go.Figure()\n",
+ "noises = np.arange(-0.5, 0.51, 0.01)\n",
+ "# Add traces, one for each slider step\n",
+ "for noise in noises:\n",
+ " # Noisy unwrapped data\n",
+ " unwrapped_noisy = copy.deepcopy(y_unwrapped)\n",
+ " unwrapped_noisy[1] += noise\n",
+ " fig.add_trace(\n",
+ " go.Scatter(\n",
+ " visible=False,\n",
+ " mode='markers',\n",
+ " marker=dict(color='red'),\n",
+ " name=\"True Phase\",\n",
+ " x=t,\n",
+ " y=unwrapped_noisy))\n",
+ "\n",
+ " # Fit of noisy unwrapped data\n",
+ " reg1 = LinearRegression().fit(t.reshape(-1, 1), unwrapped_noisy.reshape(-1,1))\n",
+ " fieldmap_rad1 = reg1.coef_[0] # [rad / s]\n",
+ " fieldmap_intercept1 = reg1.intercept_[0] # [rad / s]\n",
+ " t_predict1 = np.array([beg, end])\n",
+ " y_predict1 = reg1.predict(t_predict1.reshape(-1,1))[:,0]\n",
+ " fig.add_trace(go.Scatter(visible=False, x=t_predict1, y=y_predict1, mode='lines', marker=dict(color='red'), name='True linear fit'))\n",
+ "\n",
+ " # Noisy wrapped data\n",
+ " wrapped_noisy = copy.deepcopy(unwrapped_noisy)\n",
+ " wrapped_noisy[0] = np.angle(np.exp(1j*wrapped_noisy[0]))\n",
+ " wrapped_noisy[1] = np.angle(np.exp(1j*wrapped_noisy[1]))\n",
+ " wrapped_noisy[2] = np.angle(np.exp(1j*wrapped_noisy[2]))\n",
+ " fig.add_trace(\n",
+ " go.Scatter(\n",
+ " visible=False,\n",
+ " mode='markers',\n",
+ " marker=dict(color='blue'),\n",
+ " name=\"Wrapped Phase\",\n",
+ " x=t,\n",
+ " y=wrapped_noisy))\n",
+ "\n",
+ " # UMPIRE\n",
+ " unwrapped_umpire = umpire_3echoes(wrapped_noisy, t)\n",
+ " fig.add_trace(\n",
+ " go.Scatter(\n",
+ " visible=False,\n",
+ " x=t, y=unwrapped_umpire, mode='markers', marker=dict(color='green'), name='Umpire'))\n",
+ "\n",
+ " # Fit unwrapped data of UMPIRE\n",
+ " reg2 = LinearRegression().fit(t.reshape(-1, 1), unwrapped_umpire.reshape(-1,1))\n",
+ " # Slope of linear regression reshaped into the shape of original 3D phase.\n",
+ " fieldmap_rad2 = reg2.coef_[0] # [rad / s]\n",
+ " fieldmap_intercept2 = reg2.intercept_[0] # [rad / s]\n",
+ " t_predict2 = np.array([beg, end])\n",
+ " y_predict2 = reg2.predict(t_predict2.reshape(-1,1))[:,0]\n",
+ " \n",
+ " fig.add_trace(go.Scatter(visible=False, x=t_predict2, y=y_predict2, mode='lines', marker=dict(color='green'), name='Umpire fit'))\n",
+ " \n",
+ "active = len(noises)//2\n",
+ "fig.data[active].visible = True\n",
+ "fig.data[active+1].visible = True\n",
+ "fig.data[active+2].visible = True\n",
+ "fig.data[active+3].visible = True\n",
+ "fig.data[active+4].visible = True\n",
+ "\n",
+ "# Static plot\n",
+ "fig.add_trace(go.Scatter(x=[beg, end], y=[math.pi, math.pi], mode='lines', line=dict(color='gray'), showlegend=False))\n",
+ "fig.add_trace(go.Scatter(x=[beg, end], y=[-math.pi, -math.pi], mode='lines', line=dict(color='gray'), showlegend=False))\n",
+ "fig.update_xaxes(title_text=\"Time (ms)\", range=[beg, end])\n",
+ "fig.update_yaxes(title_text=\"Phase (rad)\", tickmode = 'array', range=[-4,25],\n",
+ " tickvals = [-2*math.pi, 0, 2*math.pi, 4*math.pi, 6*math.pi],\n",
+ " ticktext = [f'-2{PI_UNICODE}', '0', f'2{PI_UNICODE}', f'4{PI_UNICODE}', f'6{PI_UNICODE}'])\n",
+ "\n",
+ "# Create and add slider\n",
+ "phase_offsets = [f\"{i:.2}\" for i in noises]\n",
+ "steps = []\n",
+ "for i in range(len(noises)):\n",
+ " step = dict(\n",
+ " method=\"update\",\n",
+ " label=phase_offsets[i],\n",
+ " args=[{\"visible\": [False] * 5*len(noises) + [True] * (len(fig.data) - 5*len(noises))}], # layout attribute\n",
+ " )\n",
+ " step[\"args\"][0][\"visible\"][5*i] = True\n",
+ " step[\"args\"][0][\"visible\"][5*i + 1] = True\n",
+ " step[\"args\"][0][\"visible\"][5*i + 2] = True\n",
+ " step[\"args\"][0][\"visible\"][5*i + 3] = True\n",
+ " step[\"args\"][0][\"visible\"][5*i + 4] = True\n",
+ " steps.append(step)\n",
+ "\n",
+ "sliders = [dict(\n",
+ " active=active,\n",
+ " currentvalue={\"prefix\": \"2nd echo phase offset (rad): \"},\n",
+ " pad={\"t\": 50},\n",
+ " steps=steps\n",
+ ")]\n",
+ "\n",
+ "fig.update_layout({\"width\": 800},\n",
+ " sliders=sliders,\n",
+ " title_text=\"Effect of noise using UMPIRE phase unwrapping\", title_x=0.5)\n",
+ "\n",
+ "fig.show()"
+ ]
+ }
+ ],
+ "metadata": {
+ "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.0"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/05-B0/content/data.zip b/05-B0/content/data.zip
new file mode 100644
index 0000000..990ffe2
Binary files /dev/null and b/05-B0/content/data.zip differ
diff --git a/05-B0/content/figure.html b/05-B0/content/figure.html
new file mode 100644
index 0000000..7b237e8
--- /dev/null
+++ b/05-B0/content/figure.html
@@ -0,0 +1,14 @@
+
+
+
+
+
+
\ No newline at end of file
diff --git a/05-B0/content/index.ipynb b/05-B0/content/index.ipynb
new file mode 100644
index 0000000..fcefa5b
--- /dev/null
+++ b/05-B0/content/index.ipynb
@@ -0,0 +1,187 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# MRathon 2024 - NeuroLibre preprint project\n",
+ "\n",
+ "Hello world!\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Figure 1"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "tags": [
+ "hide_input",
+ "remove_output"
+ ]
+ },
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "from scipy.fft import fftn, ifftn, fftshift\n",
+ "import plotly.express as px\n",
+ "import plotly.graph_objects as go\n",
+ "from plotly.subplots import make_subplots\n",
+ "from IPython.display import display, HTML\n",
+ "from plotly.offline import init_notebook_mode, iplot, plot\n",
+ "\n",
+ "init_notebook_mode(connected=True)\n",
+ "vol_size = (128, 128, 128)\n",
+ "vol_centre = (vol_size[1] // 2, vol_size[2] // 2)\n",
+ "offset = (10, 20)\n",
+ "radius_inner = 10\n",
+ "radius_outer = 40\n",
+ "cylinder_base = np.zeros(vol_size)\n",
+ "cylinder_internal = np.zeros(vol_size)\n",
+ "\n",
+ "x = np.arange(-vol_size[1]//2,vol_size[1]//2,dtype=float)\n",
+ "y = x\n",
+ "z = x\n",
+ "X,Y,Z = np.meshgrid(x,y,z)\n",
+ "\n",
+ "cylinder_base[X**2+Y**2<=radius_outer**2] = 1\n",
+ "\n",
+ "cylinder_internal[(X-offset[0])**2+(Y-offset[1])**2<=radius_inner**2] += 0.1\n",
+ "\n",
+ "kx = X\n",
+ "ky = Y\n",
+ "kz = Z\n",
+ "KK = (kx**2 + ky**2 + kz**2)\n",
+ "KK[KK==0] = np.nan\n",
+ "dipole_kernel = 2./3 - (kx**2)/(KK)\n",
+ "dipole_kernel[(kx**2 + ky**2 + kz**2) == 0] = 0\n",
+ "\n",
+ "fieldmap = np.real(ifftn(fftshift(dipole_kernel) * (fftn(cylinder_internal))))\n",
+ "\n",
+ "idx = 64\n",
+ "\n",
+ "fig = make_subplots(rows=1, cols=3, shared_xaxes=False, horizontal_spacing=0.1,\n",
+ " subplot_titles=(\"Geometry\",\n",
+ " \"Dipole kernel\",\n",
+ " \"Fieldmap\"),\n",
+ " specs=[[{\"type\": \"Heatmap\"}, {\"type\": \"Heatmap\"}, {\"type\": \"Heatmap\"}]])\n",
+ "fig.add_trace(go.Heatmap(z=cylinder_base[:,:,idx]+cylinder_internal[:,:,64], colorscale='gray'), 1, 1)\n",
+ "fig.add_trace(go.Heatmap(z=dipole_kernel[:,:,idx], colorscale='gray'), 1, 2)\n",
+ "fig.add_trace(go.Heatmap(z=fieldmap[:,:,idx], colorscale='gray'), 1, 3)\n",
+ "\n",
+ "plot(fig, filename = 'geometry.html')\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "tags": [
+ "remove_input",
+ "report_output"
+ ]
+ },
+ "outputs": [],
+ "source": [
+ "display(HTML('geometry.html'))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Figure 2"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "tags": [
+ "hide_input",
+ "remove_output"
+ ]
+ },
+ "outputs": [],
+ "source": [
+ "fig = go.Figure()\n",
+ "\n",
+ "echo_times = [10, 30, 50]\n",
+ "for step in echo_times:\n",
+ " fig.add_trace(\n",
+ " go.Heatmap(\n",
+ " visible=False,\n",
+ " z=(fieldmap[:,:,idx]*step) % np.pi,\n",
+ " colorscale='gray'))\n",
+ "\n",
+ "fig.data[0].visible = True\n",
+ "\n",
+ "# Create and add slider\n",
+ "steps = []\n",
+ "for i in range(len(fig.data)):\n",
+ " step = dict(\n",
+ " method=\"update\",\n",
+ " args=[{\"visible\": [False] * len(fig.data)},\n",
+ " {\"title\": \"TE=\" + str(echo_times[i]) + \" ms\"}], # layout attribute\n",
+ " )\n",
+ " step[\"args\"][0][\"visible\"][i] = True # Toggle i'th trace to \"visible\"\n",
+ " steps.append(step)\n",
+ "\n",
+ "sliders = [dict(\n",
+ " active=10,\n",
+ " currentvalue={\"prefix\": \"Frequency: \"},\n",
+ " pad={\"t\": 50},\n",
+ " steps=steps\n",
+ ")]\n",
+ "\n",
+ "fig.update_layout(\n",
+ " sliders=sliders\n",
+ ")\n",
+ "\n",
+ "plot(fig, filename = 'cylinder.html')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "tags": [
+ "remove_input",
+ "report_output"
+ ]
+ },
+ "outputs": [],
+ "source": [
+ "display(HTML('cylinder.html'))"
+ ]
+ }
+ ],
+ "metadata": {
+ "celltoolbar": "Tags",
+ "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.0"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/05-B0/content/logo.png b/05-B0/content/logo.png
new file mode 100644
index 0000000..06d56f4
Binary files /dev/null and b/05-B0/content/logo.png differ
diff --git a/05-B0/paper.bib b/05-B0/paper.bib
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/05-B0/paper.bib
@@ -0,0 +1 @@
+
diff --git a/05-B0/paper.md b/05-B0/paper.md
new file mode 100644
index 0000000..50350e9
--- /dev/null
+++ b/05-B0/paper.md
@@ -0,0 +1,29 @@
+---
+title: B0 mapping
+
+tags:
+ - NeuroLibre
+authors:
+ - name: MRathon Group
+ orcid:
+ affiliation: 1
+
+affiliations:
+ - name: Singapore
+ index: 1
+
+
+date: 3 May 2024
+bibliography: paper.bib
+
+---
+
+# Summary
+
+# Statement of need
+
+# Figures
+
+# Acknowledgements
+
+## References