# Prepare Python environment
import numpy as np
import plotly.express as px
import os
import nibabel as nib
import numpy as np
from plotly.subplots import make_subplots
import plotly.graph_objects as go
import math
import json
import scipy.io as sio
from pathlib import Path
data_dir = Path("../../../data/05-B0/data/fmap")
PI_UNICODE = "\U0001D70B"
fname_mag_e1 = data_dir / "sub-fmap_magnitude1.nii.gz"
fname_phase_e1 = data_dir / "sub-fmap_phase1.nii.gz"
fname_phase_e1_json = data_dir / "sub-fmap_phase1.json"
fname_phase_e2 = data_dir / "sub-fmap_phase2.nii.gz"
fname_mask = data_dir / "mask.nii.gz"
fname_fmap = data_dir / "fmap.nii.gz"
nii_mag_e1 = nib.load(fname_mag_e1)
nii_phase_e1 = nib.load(fname_phase_e1)
nii_phase_e2 = nib.load(fname_phase_e2)
nii_mask = nib.load(fname_mask)
nii_fmap = nib.load(fname_fmap)
#| label: b0Fig8jn
#| label: b0Fig7jn
# Prepare Python environment
import numpy as np
import plotly.express as px
import os
import nibabel as nib
import numpy as np
from plotly.subplots import make_subplots
import plotly.graph_objects as go
import math
import json
import scipy.io as sio
from pathlib import Path
data_dir = Path("../../../data/05-B0/data/fmap")
PI_UNICODE = "\U0001D70B"
fname_mag_e1 = data_dir / "sub-fmap_magnitude1.nii.gz"
fname_phase_e1 = data_dir / "sub-fmap_phase1.nii.gz"
fname_phase_e1_json = data_dir / "sub-fmap_phase1.json"
fname_phase_e2 = data_dir / "sub-fmap_phase2.nii.gz"
fname_mask = data_dir / "mask.nii.gz"
fname_fmap = data_dir / "fmap.nii.gz"
nii_mag_e1 = nib.load(fname_mag_e1)
nii_phase_e1 = nib.load(fname_phase_e1)
nii_phase_e2 = nib.load(fname_phase_e2)
nii_mask = nib.load(fname_mask)
nii_fmap = nib.load(fname_fmap)
import numpy as np
import plotly.express as px
import os
import nibabel as nib
import numpy as np
from plotly.subplots import make_subplots
import plotly.graph_objects as go
import math
import json
PI_UNICODE = "\U0001D70B"
def complex_difference(phase1, phase2):
""" Calculates the complex difference between 2 phase arrays (phase2 - phase1)
Args:
phase1 (numpy.ndarray): Array containing phase data in radians
phase2 (numpy.ndarray): Array containing phase data in radians. Must be the same shape as phase1.
Returns:
numpy.ndarray: The difference in phase between each voxels of phase2 and phase1 (phase2 - phase1)
"""
# Calculate phasediff using complex difference
comp_0 = np.exp(-1j * phase1)
comp_1 = np.exp(1j * phase2)
return np.angle(comp_0 * comp_1)
n=2
def plot_2_echo_fmap(phase1, phase2, echotime1, echotime2):
phase_diff = complex_difference(phase1, phase2)
fmap = phase_diff / (echotime2 - echotime1) / 2 / math.pi
# Attempt at subplots
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"}]], )
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,
colorbar=dict(title="Rad",
titleside="top",
tickmode="array",
tickvals=[-math.pi, 0, math.pi-0.01],
ticktext = [f"-{PI_UNICODE}", 0, f'{PI_UNICODE}'])), 1, 1)
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,
colorbar=dict(title="Rad",
titleside="top",
tickmode="array",
tickvals=[-math.pi, 0, math.pi-0.01],
ticktext = [f"-{PI_UNICODE}", 0, f'{PI_UNICODE}'])), 1, 2)
fig.update_xaxes(showticklabels=False)
fig.update_yaxes(showticklabels=False)
fig.update_layout({"height": 450, "width": 750})
fig.show()
def plot_2_echo_fmap_bottom(phase1, phase2, echotime1, echotime2):
phase_diff = complex_difference(phase1, phase2)
fmap = phase_diff / (echotime2 - echotime1) / 2 / math.pi
# Attempt at subplots
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"}]], )
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,
colorbar=dict(title="Rad",
titleside="top",
tickmode="array",
tickvals=[-math.pi, 0, math.pi-0.01],
ticktext = [f"-{PI_UNICODE}", 0, f'{PI_UNICODE}'])), 1, 1)
fig.add_trace(go.Heatmap(z=np.rot90(fmap, k=-1), colorscale='gray',colorbar_x=2/n - 0.02,
colorbar=dict(title="Hz",
titleside="top")), 1, 2)
fig.update_xaxes(showticklabels=False)
fig.update_yaxes(showticklabels=False)
fig.update_layout({"height": 450, "width": 750})
fig.show()
def get_circle(x, y, r):
if x < 1 or y < 1 or r < 1:
raise ValueError("Input parameters are too small")
my_array = np.zeros([x,y])
for i in range(x):
for j in range(y):
squared = (i-(x/2))**2 + (j-(y/2))**2
h = np.sqrt(squared)
if h < r:
my_array[i,j] = 1
return my_array
echo1 = get_circle(100, 100, 30) * -1
echo2 = get_circle(100, 100, 30) * 2
echo_time1 = 0.005
echo_time2 = 0.01
mask = nii_mask.get_fdata()[30:-30,8:105,30]
phase1 = (nii_phase_e1.get_fdata()[30:-30,8:105,30] / 4095 * 2 * math.pi - math.pi) * mask
phase2 = (nii_phase_e2.get_fdata()[30:-30,8:105,30] / 4095 * 2 * math.pi - math.pi) * mask
echo_time1 = 0.00338
echo_time2 = 0.00558
plot_2_echo_fmap(phase1, phase2, echo_time1, echo_time2)
Loading...
plot_2_echo_fmap_bottom(phase1, phase2, echo_time1, echo_time2)Loading...