# 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
from contextlib import contextmanager
import sys, os
from pathlib import Path
@contextmanager
def suppress_stdout():
with open(os.devnull, "w") as devnull:
old_stdout = sys.stdout
sys.stdout = devnull
try:
yield
finally:
sys.stdout = old_stdout
import os
from pathlib import Path
def find_myst_yml_directories(start_dir=None):
"""
Recursively search for directories containing myst.yml file.
Args:
start_dir (str or Path): Starting directory (defaults to current directory)
Returns:
list: List of full paths to directories containing myst.yml
"""
if start_dir is None:
start_dir = Path.cwd()
else:
start_dir = Path(start_dir)
myst_dirs = []
def _search_directory(current_dir):
# Check if myst.yml exists in current directory
myst_file = current_dir / "myst.yml"
if myst_file.exists():
myst_dirs.append(str(current_dir.resolve()))
# Don't search subdirectories if we found myst.yml here
return
# Recursively search all subdirectories
for item in current_dir.iterdir():
if item.is_dir():
try:
_search_directory(item)
except (PermissionError, OSError):
# Skip directories we can't access
continue
_search_directory(start_dir)
return myst_dirs
def find_myst_yml_directories_upwards(start_dir=None):
"""
Search for myst.yml in current directory, if not found go to parent and repeat.
Args:
start_dir (str or Path): Starting directory (defaults to current directory)
Returns:
str or None: Full path of directory containing myst.yml, or None if not found
"""
if start_dir is None:
current_dir = Path.cwd()
else:
current_dir = Path(start_dir)
# Keep going up until we reach the filesystem root
while current_dir != current_dir.parent: # Stop at root
myst_file = current_dir / "myst.yml"
if myst_file.exists():
return str(current_dir.resolve())
# Move to parent directory
current_dir = current_dir.parent
return None
with suppress_stdout():
repo_path = Path(find_myst_yml_directories_upwards())
print(repo_path)
data_req_path = repo_path / "binder" / "data_requirement.json"
data_path = repo_path / "data"
dataset_path = data_path / "qmrlab-mooc"
data_dir = dataset_path / "05-B0" / "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
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...