Skip to content

Commit

Permalink
Fix missing figs
Browse files Browse the repository at this point in the history
  • Loading branch information
@mathieuboudreau
mathieuboudreau committed Oct 9, 2024
1 parent 22174a1 commit 550010d
Show file tree
Hide file tree
Showing 4 changed files with 327 additions and 299 deletions.
4 changes: 4 additions & 0 deletions 5 B0 Mapping/2 Dual echo B0 mapping/03-Single Frequency Population.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,10 @@ numbering:

To build intuition about field maps, let us imagine a sample at a constant offset frequency from {math}`f_{0}` . Note that this simplistic representation of the field typically does not occur due to how the susceptibilities of the neighboring regions interact with one another to create the _B_{sub}`0` field offset (see the [_B_{sub}`0` inhomogeneity section](#b0InhomoIntro), but is shown as such for learning purposes. The sample is shown as a circle in [](#b0Plot8). As the frequency is not at the [Larmor frequency](https://en.wikipedia.org/wiki/Larmor_precession), phase accumulation is observed at the different echo times and a phase difference map can be computed. The _B_{sub}`0` field map is then calculated using the echo times. Note that if {math}`\Delta`TE is too long, the phase could make more than a half revolution between the two echo times resulting in an erroneous _B_{sub}`0` field estimation. This is because phase is defined over {math}`-\pi` to {math}`+\pi` and the sampled points should respect the [Nyquist criteria](https://en.wikipedia.org/wiki/Nyquist_frequency). In practice, this example field (constant offset) could easily be corrected by adjusting the transmit and receive frequency of the scanner.

:::{figure} #fig5p8cella
:label: b0Plot8a
:::

:::{figure} #fig5p8cell
:label: b0Plot8
:enumerator: 5.8
Expand Down
4 changes: 4 additions & 0 deletions 5 B0 Mapping/2 Dual echo B0 mapping/04-Multi-Frequency Population.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,10 @@ numbering:

A brain dataset is used to show a concrete example of a field map that could be acquired in practice. [](#b0Plot9) shows two phase images where phase accumulation is shown due to frequency offsets that vary spatially. As mentioned previously, phase wraps are visible where phase transitions from {math}`-\pi` to {math}`+\pi` and will be discussed in more detail in the next chapter. The phase difference and _B_{sub}`0` field maps are also shown. Note that taking the phase difference eliminates the wraps in this example, however, there could be residual wraps when the field is more inhomogeneous.

:::{figure} #fig5p9cella
:label: b0Plot9a
:::

:::{figure} #fig5p9cell
:label: b0Plot9
:enumerator: 5.9
Expand Down
270 changes: 140 additions & 130 deletions Figures/B0/fig5p8.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21451,7 +21451,145 @@
},
"metadata": {},
"output_type": "display_data"
},
}
],
"source": [
"#| label: fig5p8cella\n",
"\n",
"# Prepare Python environment\n",
"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",
"\n",
"import scipy.io as sio\n",
"from pathlib import Path\n",
"\n",
"data_dir = Path(\"../../data/05-B0/data/fmap\")\n",
"\n",
"PI_UNICODE = \"\\U0001D70B\"\n",
"fname_mag_e1 = data_dir / \"sub-fmap_magnitude1.nii.gz\"\n",
"fname_phase_e1 = data_dir / \"sub-fmap_phase1.nii.gz\"\n",
"fname_phase_e1_json = data_dir / \"sub-fmap_phase1.json\"\n",
"fname_phase_e2 = data_dir / \"sub-fmap_phase2.nii.gz\"\n",
"fname_mask = data_dir / \"mask.nii.gz\"\n",
"fname_fmap = data_dir / \"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)\n",
"\n",
"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\"\n",
"\n",
"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)\n",
"n=2\n",
"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",
" # Attempt at subplots\n",
" fig = make_subplots(rows=1, cols=2, shared_xaxes=False, horizontal_spacing=0.1, subplot_titles=(\"Phase 1\", \"Phase 2\"), specs=[[{\"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",
" \n",
" fig.update_xaxes(showticklabels=False)\n",
" fig.update_yaxes(showticklabels=False)\n",
" fig.update_layout({\"height\": 450, \"width\": 750})\n",
" \n",
" fig.show()\n",
"\n",
"\n",
"def plot_2_echo_fmap_bottom(phase1, phase2, echotime1, echotime2):\n",
" phase_diff = complex_difference(phase1, phase2)\n",
" fmap = phase_diff / (echotime2 - echotime1) / 2 / math.pi\n",
"\n",
" # Attempt at subplots\n",
" fig = make_subplots(rows=1, cols=2, shared_xaxes=False, horizontal_spacing=0.1, subplot_titles=(\"Phase difference\", \"B0 field map\"), specs=[[{\"type\": \"Heatmap\"}, {\"type\": \"Heatmap\"}]], )\n",
" \n",
" fig.add_trace(go.Heatmap(z=np.rot90(phase_diff, 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(fmap, k=-1), colorscale='gray',colorbar_x=2/n - 0.02,\n",
" colorbar=dict(title=\"Hz\",\n",
" titleside=\"top\")), 1, 2)\n",
" \n",
" fig.update_xaxes(showticklabels=False)\n",
" fig.update_yaxes(showticklabels=False)\n",
" fig.update_layout({\"height\": 450, \"width\": 750})\n",
" \n",
" fig.show()\n",
"\n",
"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)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.plotly.v1+json": {
Expand Down Expand Up @@ -42889,135 +43027,7 @@
],
"source": [
"#| label: fig5p8cell\n",
"\n",
"# Prepare Python environment\n",
"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",
"\n",
"import scipy.io as sio\n",
"from pathlib import Path\n",
"\n",
"data_dir = Path(\"../../data/05-B0/data/fmap\")\n",
"\n",
"PI_UNICODE = \"\\U0001D70B\"\n",
"fname_mag_e1 = data_dir / \"sub-fmap_magnitude1.nii.gz\"\n",
"fname_phase_e1 = data_dir / \"sub-fmap_phase1.nii.gz\"\n",
"fname_phase_e1_json = data_dir / \"sub-fmap_phase1.json\"\n",
"fname_phase_e2 = data_dir / \"sub-fmap_phase2.nii.gz\"\n",
"fname_mask = data_dir / \"mask.nii.gz\"\n",
"fname_fmap = data_dir / \"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)\n",
"\n",
"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\"\n",
"\n",
"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)\n",
"n=2\n",
"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",
" # Attempt at subplots\n",
" fig = make_subplots(rows=1, cols=2, shared_xaxes=False, horizontal_spacing=0.1, subplot_titles=(\"Phase 1\", \"Phase 2\"), specs=[[{\"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",
" \n",
" fig.update_xaxes(showticklabels=False)\n",
" fig.update_yaxes(showticklabels=False)\n",
" fig.update_layout({\"height\": 450, \"width\": 750})\n",
" \n",
" fig.show()\n",
"\n",
"\n",
"def plot_2_echo_fmap_bottom(phase1, phase2, echotime1, echotime2):\n",
" phase_diff = complex_difference(phase1, phase2)\n",
" fmap = phase_diff / (echotime2 - echotime1) / 2 / math.pi\n",
"\n",
" # Attempt at subplots\n",
" fig = make_subplots(rows=1, cols=2, shared_xaxes=False, horizontal_spacing=0.1, subplot_titles=(\"Phase difference\", \"B0 field map\"), specs=[[{\"type\": \"Heatmap\"}, {\"type\": \"Heatmap\"}]], )\n",
" \n",
" fig.add_trace(go.Heatmap(z=np.rot90(phase_diff, 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(fmap, k=-1), colorscale='gray',colorbar_x=2/n - 0.02,\n",
" colorbar=dict(title=\"Hz\",\n",
" titleside=\"top\")), 1, 2)\n",
" \n",
" fig.update_xaxes(showticklabels=False)\n",
" fig.update_yaxes(showticklabels=False)\n",
" fig.update_layout({\"height\": 450, \"width\": 750})\n",
" \n",
" fig.show()\n",
"\n",
"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)\n",
"plot_2_echo_fmap_bottom(echo1, echo2, echo_time1, echo_time2)\n",
"\n"
"plot_2_echo_fmap_bottom(echo1, echo2, echo_time1, echo_time2)\n"
]
}
],
Expand Down
Loading

0 comments on commit 550010d

Please sign in to comment.