OSCR

Frame-wise multi-echo distortion correction for superior functional MRI.

Code ↔ Paper

20 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 20 matches · 2 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › Processing pipeline ↔ me_pipeline/params.py, lines 406–535 · score 0.93 · nuisance regression, FD threshold, Frame censoring, Bias field, N4, FSL
  2. [2] § Methods › Multi-echo distortion correction (MEDIC) › Temporal phase correction ↔ warpkit/unwrap.py, lines 495–579 · score 0.66 · ensure phase unwrapping, magnitude image, minimizes, correlation, temporal, frames
  3. [3] § Results › MEDIC dynamic distortion correction reduces the impact of head motion on functional connectivity estimates ↔ medic_analysis/scripts/paper_figures.py, lines 464–523 · score 0.63 · dorsolateral prefrontal cortex, seed map, low motion, DLPFC, exemplar, functional connectivity
  4. [4] § Methods › Processing pipeline ↔ tools/install_4dfp.sh, lines 89–139 · score 0.61 · step resampling, FSL, tools, MNI152, dfp, threshold
  5. [5] § Methods › Data acquisition › ABCD dataset ↔ me_pipeline/params.py, lines 406–535 · score 0.60 · frame censoring, nuisance, FSL, variables, filtered, motion
  6. [6] § Results › MEDIC distortion correction is superior on local and global anatomical alignment metrics ↔ medic_analysis/scripts/paper_figures.py, lines 1280–1402 · score 0.60 · T2w NMI, global metrics, T1w, correlation, alignment, TOPUP
  7. [7] § Results › MEDIC dynamic distortion correction improves functional connectivity in pediatric populations ↔ medic_analysis/scripts/paper_figures.py, lines 772–846 · score 0.59 · occipital cortex, Seed maps, functional connectivity, ABCD, distortion correction, dynamic
  8. [8] § Results › MEDIC frame-wise distortion correction produces superior anatomical alignment ↔ medic_analysis/scripts/paper_figures.py, lines 1094–1178 · score 0.58 · UMinn, WashU, MEDIC TOPUP, arrows, Penn, field maps
  9. [9] § Methods › Multi-echo distortion correction (MEDIC) ↔ warpkit/distortion.py, the whole file · a weak match · score 0.57 · undistorted space, phase images, reconstruction, field map, echo, MEDIC
  10. [10] § Methods › Multi-echo distortion correction (MEDIC) › Phase offset correction and unwrapping ↔ warpkit/distortion.py, the whole file · a weak match · score 0.55 · ROMEO algorithm, Phase unwrapping, echoes, MEDIC
  11. [11] § Results › MEDIC distortion correction is superior on local and global anatomical alignment metrics ↔ medic_analysis/scripts/paper_figures.py, lines 1280–1402 · score 0.55 · Segmentation metrics, alignment metric, AUC, spotlight, T2w, T1w
  12. [12] § Results › MEDIC distortion correction is superior on local and global anatomical alignment metrics ↔ medic_analysis/scripts/alignment_metrics.py, lines 128–210 · score 0.54 · normalized mutual information, alignment metric, gradient, spotlight, correlation, MEDIC
  13. [13] § Methods › Multi-echo distortion correction (MEDIC) › Phase offset correction and unwrapping ↔ warpkit/unwrap.py, lines 82–135 · score 0.54 · unwrapped phase, phase offset, MCPC, echoes
  14. [14] § Results › MEDIC dynamic distortion correction improves functional connectivity in pediatric populations ↔ medic_analysis/scripts/paper_figures.py, lines 772–846 · score 0.52 · occipital cortex, Seed maps, ABCD, dynamic, correlations, TOPUP
  15. [15] § Results › MEDIC distortion correction is superior on local and global anatomical alignment metrics ↔ medic_analysis/scripts/alignment_metrics.py, lines 128–210 · score 0.52 · normalized mutual information, NMI, metrics, gradient, correlation, alignment
  16. [16] § Methods › Multi-echo distortion correction (MEDIC) › Weighted field map computation ↔ warpkit/unwrap.py, lines 495–579 · score 0.52 · weighted linear regression, model, magnitude, voxels, echo, map
  17. [17] § Results › MEDIC dynamic distortion correction reduces the impact of head motion on functional connectivity estimates ↔ medic_analysis/scripts/paper_figures.py, lines 464–523 · score 0.52 · prefrontal cortex, low motion, DLPFC, Functional connectivity, distortion corrected, seeds
  18. [18] § Results › MEDIC captures magnetic field changes due to head motion ↔ medic_analysis/scripts/paper_figures.py, lines 1581–1706 · score 0.51 · motion parameters, head position, neutral, Rotation, field map, dynamic
  19. [19] § Results › MEDIC frame-wise distortion correction produces superior anatomical alignment ↔ medic_analysis/scripts/paper_figures.py, lines 1094–1178 · score 0.51 · UMinn, WashU, Penn, TOPUP, MEDIC
  20. [20] § Methods › Multi-echo distortion correction (MEDIC) › Displacement field inversion ↔ include/itk/itkModifiedInvertDisplacementFieldImageFilter.h, lines 44–179 · score 0.51 · Displacement field, ITK, inverted, inverse, space

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,752 lines · 61 KB · MIT · 9 matches

  1. """Main script for generating paper figures.
  3. See `paper_figures --help` for more information.
  4. """
  5. import argparse
  6. import json
  7. import os
  8. from pathlib import Path
  10. import matplotlib as mpl
  11. import matplotlib.patches as patches
  12. import matplotlib.pyplot as plt
  13. import nibabel as nib
  14. import numpy as np
  15. import pandas as pd
  16. import seaborn as sns
  17. from matplotlib.gridspec import GridSpec, GridSpecFromSubplotSpec
  18. from nilearn.plotting.cm import _cmap_d as nilearn_cmaps
  19. from PIL import Image
  20. from scipy.stats import ttest_rel
  21. from skimage.exposure import equalize_hist
  23. from medic_analysis.common.figures import data_plotter, hz_limits_to_mm, render_dynamic_figure
  25. from . import FIGURES_DIR, MM_TO_INCHES
  27. # Set global seaborn figure settings
  28. GLOBAL_SETTINGS = {
  29. "font": "Satoshi",
  30. "font_scale": 1,
  31. "palette": "pastel",
  32. "style": "white",
  33. "rc": {
  34. "figure.dpi": 150,
  35. "figure.titlesize": 7,
  36. "font.size": 7,
  37. "axes.titlesize": 6,
  38. "axes.titlepad": 0,
  39. "axes.labelsize": 6,
  40. "axes.labelpad": 0,
  41. "axes.linewidth": 0.5,
  42. "legend.title_fontsize": 6,
  43. "legend.fontsize": 6,
  44. "xtick.labelsize": 6,
  45. "xtick.major.pad": 1,
  46. "xtick.major.size": 1,
  47. "ytick.labelsize": 6,
  48. "ytick.major.pad": 1,
  49. "ytick.major.size": 1,
  50. "xtick.major.width": 0.5,
  51. "ytick.major.width": 0.5,
  52. "xtick.minor.width": 0.5,
  53. "ytick.minor.width": 0.5,
  54. },
  55. }
  56. sns.set_theme(**GLOBAL_SETTINGS)
  57. LOWER_FONT_SIZE = 5
  59. # Default paths for data
  60. DATA_DIR = Path(__file__).resolve().parent.parent.parent / "data"
  61. # Head position figure
  62. FIGURE1_DATA = "/home/usr/vana/GMT2/Andrew/HEADPOSITIONSUSTEST"
  63. # Concatenated head position figure
  64. FIGURE2_DATA = "/home/usr/vana/GMT2/Andrew/HEADPOSITIONCAT"
  65. # Group Template Analysis
  66. FIGURE3_DATA = str(DATA_DIR)
  67. # Alignment and Field map Comparison
  68. FIGURE4_DATA = str(DATA_DIR)
  69. # Spotlight Analysis figure
  70. FIGURE5_DATA = str(DATA_DIR)
  71. # Alignment metrics
  72. FIGURE6_DATA = str(DATA_DIR / "alignment_metrics.csv")
  73. # Field map metrics
  74. FIGURE7_DATA = str(DATA_DIR)
  75. # tSNR figure
  76. FIGURE10_DATA = str(DATA_DIR / "tsnr.csv")
  77. # dynamic field map videos
  78. FIGURE100_DATA = "/home/usr/vana/GMT2/Andrew/HEADPOSITIONSUSTEST/derivatives"
  79. AA_DATA_DIR = Path("/data/Daenerys/ASD_ADHD/NP1173/derivatives/me_pipeline2")
  80. WASHU_DATA_DIR = Path("/net/10.20.145.34/DOSENBACH02/GMT2/Andrew/SLICETEST/derivatives/me_pipeline")
  81. PENN_DATA_DIR = Path("/net/10.20.145.34/DOSENBACH02/GMT2/Andrew/UPenn/derivatives/me_pipeline")
  82. MINN_DATA_DIR = Path("/net/10.20.145.34/DOSENBACH02/GMT2/Andrew/UMinn/derivatives")
  83. MINN_DATA_DIR2 = Path("/data/nil-bluearc/GMT/Laumann/Pilot_ME_res/BIO10001/bids/derivatives/me_pipeline")
  86. def plot_box_plot(data, variable, label, ax):
  87. p = sns.color_palette("pastel")
  88. subdata = (
  89. data[[f"{variable}_medic", f"{variable}_topup"]]
  90. .rename(columns={f"{variable}_medic": "MEDIC", f"{variable}_topup": "TOPUP"})
  91. .melt(var_name=label)
  92. )
  93. sb = sns.boxplot(
  94. data=subdata,
  95. x="value",
  96. y=label,
  97. order=["MEDIC", "TOPUP"],
  98. ax=ax,
  99. fliersize=1,
  100. linewidth=0.5,
  101. palette=[p[1], p[0]],
  102. )
  103. sb.set_xlabel("")
  104. sb.set_ylabel(label)
  105. ax.tick_params(axis="x")
  106. ax.tick_params(axis="y")
  107. return sb
  110. def draw_seed(ax, x, y, radius=30, fc="black", ec="white", linewidth=1, zorder=3):
  111. dcoords = ax.transAxes.transform((x, y))
  112. ncoords = ax.transData.inverted().transform(dcoords)
  113. ax.add_patch(patches.Circle(ncoords, radius, fc=fc, ec=ec, linewidth=linewidth, zorder=zorder))
  116. def data_to_ax(ax, corrds):
  117. return ax.transAxes.inverted().transform(ax.transData.transform(corrds))
  120. def draw_arrow(ax, loc1, loc2, color="red", linewidth=1, head_width=2, head_length=4):
  121. dloc1 = ax.transAxes.transform(loc1)
  122. nloc1 = ax.transData.inverted().transform(dloc1)
  123. dloc2 = ax.transAxes.transform(loc2)
  124. nloc2 = ax.transData.inverted().transform(dloc2)
  125. ax.add_patch(
  126. patches.FancyArrowPatch(
  127. nloc1,
  128. nloc2,
  129. arrowstyle="-|>,head_width={},head_length={}".format(head_width, head_length),
  130. color=color,
  131. linewidth=linewidth,
  132. )
  133. )
  136. # figure 1
  137. def head_position_fieldmap(data):
  138. mpl.rcParams["axes.titlesize"] = 7
  139. # Get the data
  140. output_dir = Path(data) / "derivatives"
  142. raw_func_path = (
  143. Path("/home/usr/vana/GMT2/Andrew/HEADPOSITIONCAT")
  144. / "sub-MSCHD02"
  145. / "ses-01"
  146. / "func"
  147. / "sub-MSCHD02_ses-01_task-rest_run-01_echo-1_part-mag_bold.nii.gz"
  148. )
  150. # create a list of expected labels for each run
  151. labels = [
  152. "Neutral",
  153. "+Z Rotation",
  154. "-Z Rotation",
  155. "+X Rotation",
  156. "-X Rotation",
  157. "+Y Rotation",
  158. "-Y Rotation",
  159. "Neutral to +Z Rotation",
  160. "Neutral to -Z Rotation",
  161. "Neutral to +X Rotation",
  162. "Neutral to -X Rotation",
  163. "Neutral to +Y Rotation",
  164. "Neutral to -Y Rotation",
  165. "Neutral to -Z Translation",
  166. "-Z Translation",
  167. ]
  169. # indices for run
  170. static_head_position_run_idx = [0, 1, 2, 3, 4, 5, 6, 14]
  172. # load raw_func data
  173. raw_func = nib.load(raw_func_path)
  175. # Figure 1 - Head Rotation Data
  176. # load field map files
  177. medic_fieldmaps = Path(output_dir) / "fieldmaps" / "medic_aligned"
  178. # load topup field map in neutral position as reference
  179. topup_fieldmap = nib.load(Path(output_dir) / "fieldmaps" / "topup" / "run01" / "fout.nii.gz").get_fdata()
  180. # load static field map runs
  181. static_fieldmaps = []
  182. for idx in static_head_position_run_idx:
  183. run = idx + 1
  184. static_fieldmaps.append(nib.load(medic_fieldmaps / f"run{run:02d}" / "fmap.nii.gz").dataobj)
  185. # load mask
  186. mask = nib.load(Path(output_dir) / "references" / "me_epi_ref_bet_mask.nii.gz").get_fdata()
  188. # plot range
  189. vlims = (-50, 50)
  191. f_topup = plt.figure(figsize=(90 * MM_TO_INCHES, 30 * MM_TO_INCHES), layout="constrained")
  192. gsm = GridSpec(
  193. 1,
  194. 3,
  195. left=0.025,
  196. right=0.975,
  197. bottom=0.025,
  198. top=0.975,
  199. hspace=0.03,
  200. width_ratios=[72, 110, 110],
  201. )
  202. topup_fieldmap = topup_fieldmap
  203. fmin = topup_fieldmap.min()
  204. fmax = topup_fieldmap.max()
  205. axes_list = []
  206. for j in range(3):
  207. axes_list.append(f_topup.add_subplot(gsm[j]))
  208. data_plotter(
  209. [topup_fieldmap],
  210. vmin=fmin,
  211. vmax=fmax,
  212. colormaps="gray",
  213. figure=f_topup,
  214. axes_list=axes_list,
  215. )
  216. sbs = axes_list
  217. sbs[0].set_title(r"TOPUP field map for High motion data", pad=4, weight="normal", loc="left")
  218. f_topup.savefig(FIGURES_DIR / "topup_fmap.png", dpi=300)
  220. # plot static field maps
  221. f0 = plt.figure(figsize=(180 * MM_TO_INCHES, 130 * MM_TO_INCHES), layout="constrained")
  223. # grid spec for head motion images
  224. gsm = GridSpec(
  225. 3,
  226. 7,
  227. left=0.025,
  228. right=0.975,
  229. bottom=0.625,
  230. top=0.96,
  231. wspace=0.04,
  232. hspace=0.03,
  233. height_ratios=[110, 72, 72],
  234. )
  236. # plot movement data
  237. # get min max
  238. func_min = raw_func.dataobj[..., 0].min()
  239. func_max = raw_func.dataobj[..., 0].max()
  240. # create subplots from gridspec
  241. axes_list = []
  242. for i in range(7):
  243. for j in range(3):
  244. axes_list.append(f0.add_subplot(gsm[j, i]))
  245. # plot data
  246. data_plotter( # f0.suptitle("Motion-dependent field map differences (Position - Neutral Position)")
  247. [
  248. raw_func.dataobj[..., 50],
  249. raw_func.dataobj[..., 150],
  250. raw_func.dataobj[..., 250],
  251. raw_func.dataobj[..., 350],
  252. raw_func.dataobj[..., 450],
  253. raw_func.dataobj[..., 550],
  254. raw_func.dataobj[..., 650],
  255. ],
  256. vmin=func_min,
  257. vmax=func_max,
  258. colormaps="gray",
  259. figure=f0,
  260. axes_list=axes_list,
  261. )
  262. sbs = axes_list
  263. sbs[2].set_xlabel("50", labelpad=2)
  264. sbs[20].set_xlabel("650", labelpad=2)
  265. sbs[11].set_xlabel("Frame", labelpad=2)
  266. sbs[0].set_title(r"$\bf{a}$ Functional MRI timeseries: High motion", pad=4, weight="normal", loc="left")
  268. # draw arrow line
  269. f0.canvas.draw()
  270. start = sbs[2].xaxis.label.get_window_extent()
  271. middle = sbs[11].xaxis.label.get_window_extent()
  272. end = sbs[20].xaxis.label.get_window_extent()
  273. # transform to figure coordinates
  274. start = start.transformed(f0.transFigure.inverted())
  275. middle = middle.transformed(f0.transFigure.inverted())
  276. end = end.transformed(f0.transFigure.inverted())
  277. # get the midpoints
  278. start = np.average(start.get_points(), axis=0)
  279. middle = np.average(middle.get_points(), axis=0)
  280. end = np.average(end.get_points(), axis=0)
  281. arrow1 = patches.FancyArrowPatch(
  282. start + np.array([0.01, 0]),
  283. middle - np.array([0.02, 0]),
  284. arrowstyle="-",
  285. color="black",
  286. linewidth=0.5,
  287. )
  288. arrow2 = patches.FancyArrowPatch(
  289. middle + np.array([0.02, 0]),
  290. end - np.array([0.01, 0]),
  291. arrowstyle="-|>,head_width=2,head_length=4",
  292. color="black",
  293. linewidth=0.5,
  294. )
  295. f0.add_artist(arrow1)
  296. f0.add_artist(arrow2)
  298. # get bounding box of last image
  299. bbox = sbs[-1].get_window_extent()
  300. # transform to figure coordinates
  301. bbox = bbox.transformed(f0.transFigure.inverted())
  303. # create a grid spec for the figure
  304. bottom = 0.025
  305. top = 0.575
  306. left_edge_1 = 0.09
  307. right_edge_2 = bbox.x1
  308. pad = 0.02
  309. width = (right_edge_2 - left_edge_1 - pad) / 2
  310. right_edge_1 = left_edge_1 + width
  311. left_edge_2 = right_edge_1 + pad
  312. gs0 = GridSpec(1, 1, left=0.03, right=0.09, bottom=bottom, top=top)
  313. gs_bar = GridSpecFromSubplotSpec(1, 3, wspace=0, hspace=0, width_ratios=[2, 1, 6], subplot_spec=gs0[:, :])
  314. gs1 = GridSpec(
  315. 3,
  316. 3,
  317. left=left_edge_1,
  318. right=right_edge_1,
  319. bottom=bottom,
  320. top=top,
  321. hspace=0.025,
  322. wspace=0.025,
  323. width_ratios=[72, 110, 110],
  324. )
  325. gs2 = GridSpec(
  326. 3,
  327. 3,
  328. left=left_edge_2,
  329. right=right_edge_2,
  330. bottom=bottom,
  331. top=top,
  332. hspace=0.025,
  333. wspace=0.025,
  334. width_ratios=[72, 110, 110],
  335. )
  337. # create subplots
  338. cbar_ax = f0.add_subplot(gs_bar[1])
  339. axes_list = []
  340. for i in range(3):
  341. for j in range(3):
  342. axes_list.append(f0.add_subplot(gs1[i, j]))
  343. for j in range(3):
  344. axes_list.append(f0.add_subplot(gs2[i, j]))
  346. # plot the data
  347. data_plotter(
  348. [
  349. (static_fieldmaps[1][..., 0] - static_fieldmaps[0][..., 0]) * mask,
  350. (static_fieldmaps[2][..., 0] - static_fieldmaps[0][..., 0]) * mask,
  351. (static_fieldmaps[3][..., 0] - static_fieldmaps[0][..., 0]) * mask,
  352. (static_fieldmaps[4][..., 0] - static_fieldmaps[0][..., 0]) * mask,
  353. (static_fieldmaps[5][..., 0] - static_fieldmaps[0][..., 0]) * mask,
  354. (static_fieldmaps[6][..., 0] - static_fieldmaps[0][..., 0]) * mask,
  355. ],
  356. colorbar=True,
  357. colorbar_alt_range=True,
  358. figure=f0,
  359. vmin=vlims[0],
  360. vmax=vlims[1],
  361. axes_list=axes_list,
  362. cbar_ax=cbar_ax,
  363. )
  364. sbs = axes_list
  365. sbs[0].set_title(r"$\bf{b}$ " + f"{labels[1]} (15.0 deg)", pad=4, weight="normal", loc="left")
  366. sbs[3].set_title(r"$\bf{c}$ " + f"{labels[2]} (9.8 deg)", pad=4, weight="normal", loc="left")
  367. sbs[6].set_title(r"$\bf{d}$ " + f"{labels[3]} (10.6 deg)", pad=4, weight="normal", loc="left")
  368. sbs[9].set_title(r"$\bf{e}$ " + f"{labels[4]} (13.7 deg)", pad=4, weight="normal", loc="left")
  369. sbs[12].set_title(r"$\bf{f}$ " + f"{labels[5]} (10.8 deg)", pad=4, weight="normal", loc="left")
  370. sbs[15].set_title(r"$\bf{g}$ " + f"{labels[6]} (8.6 deg)", pad=4, weight="normal", loc="left")
  372. f0.savefig(FIGURES_DIR / "fieldmap_differences.png", dpi=300)
  373. current_dir = os.getcwd()
  374. os.chdir(FIGURES_DIR)
  375. Path("figure1.png").unlink(missing_ok=True)
  376. Path("figure1.png").symlink_to("fieldmap_differences.png")
  377. os.chdir(current_dir)
  378. sns.set_theme(**GLOBAL_SETTINGS)
  381. # figure 2
  382. def head_concatenation(data):
  383. # get dataset
  384. dataset = Path(data)
  386. # load medic and topup workbench screenshots
  387. medic_scan_path = dataset / "medic_scan.png"
  388. topup_scan_path = dataset / "topup_scan.png"
  389. truth_scan_path = dataset / "truth_scan.png"
  390. medic_dlpfc_path = dataset / "medic_dlpfc.png"
  391. topup_dlpfc_path = dataset / "topup_dlpfc.png"
  392. truth_dlpfc_path = dataset / "truth_dlpfc.png"
  393. medic_occipital_path = dataset / "medic_occipital.png"
  394. topup_occipital_path = dataset / "topup_occipital.png"
  395. truth_occipital_path = dataset / "truth_occipital.png"
  397. # load data
  398. clip1 = 50
  399. clip2 = 60
  400. clipy1 = 150
  401. clipy2 = 200
  402. medic_dlpfc = np.array(Image.open(medic_dlpfc_path))
  403. medic_dlpfc_left = medic_dlpfc[clipy1:-clipy2, clip1 : medic_dlpfc.shape[1] // 2 - clip2]
  404. medic_dlpfc_right = medic_dlpfc[clipy1:-clipy2, clip2 + medic_dlpfc.shape[1] // 2 : -clip1]
  405. medic_dlpfc = np.concatenate([medic_dlpfc_left, medic_dlpfc_right], axis=1)
  406. topup_dlpfc = np.array(Image.open(topup_dlpfc_path))
  407. topup_dlpfc_left = topup_dlpfc[clipy1:-clipy2, clip1 : topup_dlpfc.shape[1] // 2 - clip2]
  408. topup_dlpfc_right = topup_dlpfc[clipy1:-clipy2, clip2 + topup_dlpfc.shape[1] // 2 : -clip1]
  409. topup_dlpfc = np.concatenate([topup_dlpfc_left, topup_dlpfc_right], axis=1)
  410. truth_dlpfc = np.array(Image.open(truth_dlpfc_path))
  411. truth_dlpfc_left = truth_dlpfc[clipy1:-clipy2, clip1 : truth_dlpfc.shape[1] // 2 - clip2]
  412. truth_dlpfc_right = truth_dlpfc[clipy1:-clipy2, clip2 + truth_dlpfc.shape[1] // 2 : -clip1]
  413. truth_dlpfc = np.concatenate([truth_dlpfc_left, truth_dlpfc_right], axis=1)
  414. medic_occipital = np.array(Image.open(medic_occipital_path))
  415. medic_occipital_left = medic_occipital[clipy1:-clipy2, clip1 : medic_occipital.shape[1] // 2 - clip2]
  416. medic_occipital_right = medic_occipital[clipy1:-clipy2, clip2 + medic_occipital.shape[1] // 2 : -clip1]
  417. medic_occipital = np.concatenate([medic_occipital_left, medic_occipital_right], axis=1)
  418. topup_occipital = np.array(Image.open(topup_occipital_path))
  419. topup_occipital_left = topup_occipital[clipy1:-clipy2, clip1 : topup_occipital.shape[1] // 2 - clip2]
  420. topup_occipital_right = topup_occipital[clipy1:-clipy2, clip2 + topup_occipital.shape[1] // 2 : -clip1]
  421. topup_occipital = np.concatenate([topup_occipital_left, topup_occipital_right], axis=1)
  422. truth_occipital = np.array(Image.open(truth_occipital_path))
  423. truth_occipital_left = truth_occipital[clipy1:-clipy2, clip1 : truth_occipital.shape[1] // 2 - clip2]
  424. truth_occipital_right = truth_occipital[clipy1:-clipy2, clip2 + truth_occipital.shape[1] // 2 : -clip1]
  425. truth_occipital = np.concatenate([truth_occipital_left, truth_occipital_right], axis=1)
  426. medic_scan = np.array(Image.open(medic_scan_path))
  427. medic_scan_left = medic_scan[clipy1:-clipy2, clip1 : medic_scan.shape[1] // 2 - clip2]
  428. medic_scan_right = medic_scan[clipy1:-clipy2, clip2 + medic_scan.shape[1] // 2 : -clip1]
  429. medic_scan = np.concatenate([medic_scan_left, medic_scan_right], axis=1)
  430. topup_scan = np.array(Image.open(topup_scan_path))
  431. topup_scan_left = topup_scan[clipy1:-clipy2, clip1 : topup_scan.shape[1] // 2 - clip2]
  432. topup_scan_right = topup_scan[clipy1:-clipy2, clip2 + topup_scan.shape[1] // 2 : -clip1]
  433. topup_scan = np.concatenate([topup_scan_left, topup_scan_right], axis=1)
  434. truth_scan = np.array(Image.open(truth_scan_path))
  435. truth_scan_left = truth_scan[clipy1:-clipy2, clip1 : truth_scan.shape[1] // 2 - clip2]
  436. truth_scan_right = truth_scan[clipy1:-clipy2, clip2 + truth_scan.shape[1] // 2 : -clip1]
  437. truth_scan = np.concatenate([truth_scan_left, truth_scan_right], axis=1)
  439. # create a figure
  440. f = plt.figure(figsize=(180 * MM_TO_INCHES, 100 * MM_TO_INCHES), layout="constrained")
  442. # create gridspec
  443. gs = GridSpec(
  444. 3,
  445. 4,
  446. left=0.005,
  447. right=0.995,
  448. bottom=0.005,
  449. top=0.975,
  450. wspace=0.15,
  451. hspace=0.01,
  452. width_ratios=[9, 9, 1, 9],
  453. )
  454. gs_cbar = GridSpecFromSubplotSpec(
  455. 3,
  456. 3,
  457. wspace=0,
  458. hspace=0,
  459. width_ratios=[3, 2, 7],
  460. height_ratios=[1, 20, 1],
  461. subplot_spec=gs[:, 2],
  462. )
  464. # plot images
  465. mpl.rcParams["axes.edgecolor"] = "white"
  466. ax_medic_dlpfc = f.add_subplot(gs[0, 0])
  467. ax_medic_dlpfc.imshow(medic_dlpfc)
  468. draw_seed(ax_medic_dlpfc, x=0.13, y=0.62)
  469. ax_medic_dlpfc.set_xticks([])
  470. ax_medic_dlpfc.set_yticks([])
  471. ax_medic_dlpfc.set_title("MEDIC: Dynamic distortion correction", pad=6, loc="center")
  472. medic_title_pos = ax_medic_dlpfc.transAxes.inverted().transform(ax_medic_dlpfc.title.get_window_extent())
  473. ax_medic_dlpfc.text(
  474. 0.5,
  475. 1,
  476. "Exemplar participant",
  477. ha="center",
  478. va="center",
  479. fontsize=5,
  480. transform=ax_medic_dlpfc.transAxes,
  481. )
  482. ax_medic_dlpfc.set_xlabel("Correlation to standard: r = 0.41", labelpad=2)
  483. ax_topup_dlpfc = f.add_subplot(gs[0, 1])
  484. ax_topup_dlpfc.imshow(topup_dlpfc)
  485. draw_seed(ax_topup_dlpfc, x=0.13, y=0.62)
  486. ax_topup_dlpfc.set_xticks([])
  487. ax_topup_dlpfc.set_yticks([])
  488. ax_topup_dlpfc.set_title("TOPUP: Static distortion correction", pad=6, loc="center")
  489. topup_title_pos = ax_topup_dlpfc.transAxes.inverted().transform(ax_topup_dlpfc.title.get_window_extent())
  490. ax_topup_dlpfc.text(
  491. 0.5,
  492. 1,
  493. "Exemplar participant",
  494. ha="center",
  495. va="center",
  496. fontsize=5,
  497. transform=ax_topup_dlpfc.transAxes,
  498. )
  499. ax_topup_dlpfc.set_xlabel("Correlation to standard: r = 0.18", labelpad=2)
  500. ax_truth_dlpfc = f.add_subplot(gs[0, 3])
  501. ax_truth_dlpfc.imshow(truth_dlpfc)
  502. draw_seed(ax_truth_dlpfc, x=0.13, y=0.62)
  503. ax_truth_dlpfc.set_xticks([])
  504. ax_truth_dlpfc.set_yticks([])
  505. ax_truth_dlpfc.set_title("Standard: Low motion (TOPUP: static)", pad=6, loc="center")
  506. truth_title_pos = ax_truth_dlpfc.transAxes.inverted().transform(ax_truth_dlpfc.title.get_window_extent())
  507. ax_truth_dlpfc.text(
  508. 0.5,
  509. 1,
  510. "Exemplar participant",
  511. ha="center",
  512. va="center",
  513. fontsize=5,
  514. transform=ax_truth_dlpfc.transAxes,
  515. )
  516. ax_pos = ax_medic_dlpfc.get_position()
  517. f.text(
  518. ax_pos.x0,
  519. ax_pos.y1 + 0.06,
  520. r"$\bf{a}$ Functional connectivity (FC) seed maps: Dorsolateral prefrontal cortex (DLPFC)",
  521. ha="left",
  522. va="center",
  523. )
  525. ax_medic_occiptal = f.add_subplot(gs[1, 0])
  526. ax_medic_occiptal.imshow(medic_occipital)
  527. draw_seed(ax_medic_occiptal, x=0.44, y=0.32)
  528. ax_medic_occiptal.set_xticks([])
  529. ax_medic_occiptal.set_yticks([])
  530. ax_medic_occiptal.set_xlabel("Correlation to standard: r = 0.53", labelpad=2)
  531. f.text(
  532. medic_title_pos[0, 0],
  533. medic_title_pos[0, 1],
  534. "MEDIC",
  535. ha="left",
  536. va="bottom",
  537. fontsize=6,
  538. transform=ax_medic_occiptal.transAxes,
  539. )
  540. ax_topup_occipital = f.add_subplot(gs[1, 1])
  541. ax_topup_occipital.imshow(topup_occipital)
  542. draw_seed(ax_topup_occipital, x=0.44, y=0.32)
  543. ax_topup_occipital.set_xticks([])
  544. ax_topup_occipital.set_yticks([])
  545. ax_topup_occipital.set_xlabel("Correlation to standard: r = 0.38", labelpad=2)
  546. f.text(
  547. topup_title_pos[0, 0],
  548. topup_title_pos[0, 1],
  549. "TOPUP",
  550. ha="left",
  551. va="bottom",
  552. fontsize=6,
  553. transform=ax_topup_occipital.transAxes,
  554. )
  555. ax_truth_occipital = f.add_subplot(gs[1, 3])
  556. ax_truth_occipital.imshow(truth_occipital)
  557. draw_seed(ax_truth_occipital, x=0.44, y=0.32)
  558. ax_truth_occipital.set_xticks([])
  559. ax_truth_occipital.set_yticks([])
  560. f.text(
  561. truth_title_pos[0, 0],
  562. truth_title_pos[0, 1],
  563. "Standard",
  564. ha="left",
  565. va="bottom",
  566. fontsize=6,
  567. transform=ax_truth_occipital.transAxes,
  568. )
  570. ax_pos = ax_medic_occiptal.get_position()
  571. f.text(
  572. ax_pos.x0,
  573. ax_pos.y1 + 0.06,
  574. r"$\bf{b}$ Functional connectivity (FC) seed maps: Occipital cortex (extrastriate visual)",
  575. ha="left",
  576. va="center",
  577. )
  579. ax_medic_scan = f.add_subplot(gs[2, 0])
  580. ax_medic_scan.imshow(medic_scan)
  581. draw_seed(ax_medic_scan, x=0.20, y=0.37)
  582. ax_medic_scan.set_xticks([])
  583. ax_medic_scan.set_yticks([])
  584. ax_medic_scan.set_xlabel("Correlation to standard: r = 0.23", labelpad=2)
  585. f.text(
  586. medic_title_pos[0, 0],
  587. medic_title_pos[0, 1],
  588. "MEDIC",
  589. ha="left",
  590. va="bottom",
  591. fontsize=6,
  592. transform=ax_medic_scan.transAxes,
  593. )
  594. ax_topup_scan = f.add_subplot(gs[2, 1])
  595. ax_topup_scan.imshow(topup_scan)
  596. draw_seed(ax_topup_scan, x=0.20, y=0.37)
  597. ax_topup_scan.set_xticks([])
  598. ax_topup_scan.set_yticks([])
  599. ax_topup_scan.set_xlabel("Correlation to standard: r = 0.18", labelpad=2)
  600. f.text(
  601. topup_title_pos[0, 0],
  602. topup_title_pos[0, 1],
  603. "TOPUP",
  604. ha="left",
  605. va="bottom",
  606. fontsize=6,
  607. transform=ax_topup_scan.transAxes,
  608. )
  609. ax_truth_scan = f.add_subplot(gs[2, 3])
  610. ax_truth_scan.imshow(truth_scan)
  611. draw_seed(ax_truth_scan, x=0.20, y=0.37)
  612. ax_truth_scan.set_xticks([])
  613. ax_truth_scan.set_yticks([])
  614. f.text(
  615. truth_title_pos[0, 0],
  616. truth_title_pos[0, 1],
  617. "Standard",
  618. ha="left",
  619. va="bottom",
  620. fontsize=6,
  621. transform=ax_truth_scan.transAxes,
  622. )
  623. ax_pos = ax_medic_scan.get_position()
  624. f.text(
  625. ax_pos.x0,
  626. ax_pos.y1 + 0.06,
  627. r"$\bf{c}$ Functional connectivity (FC) seed maps: Somato-cognitive action network (SCAN)",
  628. ha="left",
  629. va="center",
  630. )
  631. mpl.rcParams["axes.edgecolor"] = "black"
  633. # create colorbar
  634. cbar_ax = f.add_subplot(gs_cbar[1, 1])
  635. pl = cbar_ax.imshow(
  636. np.array([[-0.6, 0.6], [0.6, -0.6]]),
  637. vmin=-0.6,
  638. vmax=0.6,
  639. aspect="auto",
  640. cmap=nilearn_cmaps["roy_big_bl"],
  641. )
  642. cbar = f.colorbar(
  643. mappable=pl,
  644. cax=cbar_ax,
  645. location="left",
  646. orientation="vertical",
  647. ticks=[-0.6, -0.3, 0, 0.3, 0.6],
  648. )
  649. cbar.ax.yaxis.set_ticks_position("right")
  650. cbar.ax.invert_yaxis()
  651. # create axis for colorbar
  652. cbar.ax.set_ylabel("Functional Connectivity z(r)", labelpad=2)
  654. # # for computing correlations
  655. # # load dconn data for medic and topup
  656. # low_dconn = nib.load(
  657. # Path("/net/10.20.145.34/DOSENBACH02/GMT2/Andrew/HEADPOSITIONCAT/Pilot_ME_res/cifti_correlation_concat")
  658. # / "MSCHD02_10run_concat_ME_MNI152_T1_2mm_Swgt_norm_bpss_resid_LR_surf_subcort_32k_fsLR_brainstem_smooth1.7_corr.dconn.nii" # noqa
  659. # )
  660. # medic_dconn = nib.load(
  661. # dataset
  662. # / "derivatives"
  663. # / "me_pipeline"
  664. # / "sub-MSCHD02"
  665. # / "ses-01wNEWPROC"
  666. # / "cifti_correlation"
  667. # / "sub-MSCHD02_b1_MNI152_T1_2mm_Swgt_norm_bpss_resid_LR_surf_subcort_32k_fsLR_brainstem_surfsmooth1.7_subcortsmooth1.7.dconn.nii" # noqa
  668. # )
  669. # topup_dconn = nib.load(
  670. # dataset
  671. # / "derivatives"
  672. # / "me_pipeline"
  673. # / "sub-MSCHD02"
  674. # / "ses-01wTOPUP"
  675. # / "cifti_correlation"
  676. # / "sub-MSCHD02_b1_MNI152_T1_2mm_Swgt_norm_bpss_resid_LR_surf_subcort_32k_fsLR_brainstem_surfsmooth1.7_subcortsmooth1.7.dconn.nii" # noqa
  677. # )
  678. # # get the lower triangle of the dconn data
  679. # print("Loading dconn data...")
  680. # low_dconn_data = low_dconn.dataobj[:59412, :59412][np.tril_indices(59412)]
  681. # medic_dconn_data = medic_dconn.dataobj[:59412, :59412][np.tril_indices(59412)]
  682. # topup_dconn_data = topup_dconn.dataobj[:59412, :59412][np.tril_indices(59412)]
  683. # # remove nans
  684. # low_dconn_data[np.isnan(low_dconn_data)] = 0
  685. # medic_dconn_data[np.isnan(medic_dconn_data)] = 0
  686. # topup_dconn_data[np.isnan(topup_dconn_data)] = 0
  687. # print("Done.")
  688. # # compute correlations
  689. # print("Computing correlations...")
  690. # medic_corr = pearsonr(low_dconn_data, medic_dconn_data)
  691. # print(medic_corr)
  692. # topup_corr = pearsonr(low_dconn_data, topup_dconn_data)
  693. # print(topup_corr)
  694. # print("Done.")
  696. f.savefig(FIGURES_DIR / "head_position_concat.png", dpi=300)
  697. current_dir = os.getcwd()
  698. os.chdir(FIGURES_DIR)
  699. Path("figure2.png").unlink(missing_ok=True)
  700. Path("figure2.png").symlink_to("head_position_concat.png")
  701. os.chdir(current_dir)
  702. sns.set_theme(**GLOBAL_SETTINGS)
  705. # figure 3
  706. def group_template_comparison(data):
  707. aa_dir = Path(data)
  708. with open(aa_dir / "paircorr.json", "r") as f:
  709. data = json.load(f)
  711. # convert to dataframe
  712. df = pd.DataFrame(data)
  713. medic_similarities = df.MEDIC.to_numpy()
  714. topup_similarities = df.TOPUP.to_numpy()
  716. # get where medic is better and topup is better
  717. medic_better = medic_similarities > topup_similarities
  718. topup_better = medic_similarities < topup_similarities
  720. # get similarities where medic is better and topup is better
  721. medic_better_similarities_medic = medic_similarities[medic_better]
  722. topup_better_similarities_medic = topup_similarities[medic_better]
  723. medic_better_similarities_topup = medic_similarities[topup_better]
  724. topup_better_similarities_topup = topup_similarities[topup_better]
  726. # make figure
  727. f = plt.figure(figsize=(180 * MM_TO_INCHES, 128 * MM_TO_INCHES), layout="constrained")
  729. # create grid specs
  730. gs = GridSpec(
  731. 1,
  732. 5,
  733. left=0.005,
  734. right=0.995,
  735. bottom=0.505,
  736. top=0.995,
  737. wspace=0.075,
  738. hspace=0,
  739. width_ratios=[9, 1, 9, 1, 9],
  740. )
  741. gs_cbar = GridSpecFromSubplotSpec(
  742. 3,
  743. 3,
  744. wspace=0,
  745. hspace=0,
  746. width_ratios=[2, 1, 6],
  747. height_ratios=[2, 5, 2],
  748. subplot_spec=gs[3],
  749. )
  750. gs_bot = GridSpec(
  751. 1,
  752. 2,
  753. left=0.04,
  754. right=0.96,
  755. bottom=0.03,
  756. top=0.47,
  757. wspace=0.25,
  758. hspace=0.1,
  759. width_ratios=[1, 2],
  760. )
  761. gs_tstat = GridSpecFromSubplotSpec(1, 2, wspace=0.075, hspace=0, width_ratios=[9, 1], subplot_spec=gs_bot[1])
  762. gs_tstat_cbar = GridSpecFromSubplotSpec(
  763. 3,
  764. 3,
  765. wspace=0,
  766. hspace=0,
  767. width_ratios=[3, 1, 9],
  768. height_ratios=[1, 10, 1],
  769. subplot_spec=gs_tstat[1],
  770. )
  772. # plot surfaces
  773. clip_1 = 145
  774. clip_2 = 145
  775. medic_occipital_path = DATA_DIR / "medic_occipital_20008.png"
  776. medic_occipital = np.array(Image.open(medic_occipital_path))
  777. medic_occipital_left = medic_occipital[:, clip_1 : medic_occipital.shape[1] // 2 - clip_2]
  778. medic_occipital_right = medic_occipital[:, clip_2 + medic_occipital.shape[1] // 2 : -clip_1]
  779. medic_occipital = np.concatenate((medic_occipital_left, medic_occipital_right), axis=1)
  780. topup_occipital_path = DATA_DIR / "topup_occipital_20008.png"
  781. topup_occipital = np.array(Image.open(topup_occipital_path))
  782. topup_occipital_left = topup_occipital[:, clip_1 : topup_occipital.shape[1] // 2 - clip_2]
  783. topup_occipital_right = topup_occipital[:, clip_2 + topup_occipital.shape[1] // 2 : -clip_1]
  784. topup_occipital = np.concatenate((topup_occipital_left, topup_occipital_right), axis=1)
  785. group_template_abcd_path = DATA_DIR / "group_abcd_template_surface.png"
  786. group_template_abcd = np.array(Image.open(group_template_abcd_path))
  787. group_template_abcd_left = group_template_abcd[:, clip_1 : group_template_abcd.shape[1] // 2 - clip_2]
  788. group_template_abcd_right = group_template_abcd[:, clip_2 + group_template_abcd.shape[1] // 2 : -clip_1]
  789. group_template_abcd = np.concatenate((group_template_abcd_left, group_template_abcd_right), axis=1)
  790. mpl.rcParams["axes.edgecolor"] = "white"
  791. mpl.rcParams["font.size"] = 6
  792. axl_medic = f.add_subplot(gs[0])
  793. axl_topup = f.add_subplot(gs[2])
  794. axl_group = f.add_subplot(gs[4])
  795. axl_medic.imshow(medic_occipital)
  796. draw_seed(axl_medic, x=0.48, y=0.7)
  797. axl_medic.set_xticks([])
  798. axl_medic.set_yticks([])
  799. axl_medic.set_xlabel("Correlation to standard: r = 0.44", labelpad=6)
  800. axl_medic.set_title("MEDIC: Dynamic distortion correction", pad=6)
  801. axl_medic.text(0.5, 0.5, "Participant 1", ha="center", va="center", transform=axl_medic.transAxes)
  802. axl_topup.imshow(topup_occipital)
  803. draw_seed(axl_topup, x=0.48, y=0.7)
  804. axl_topup.set_xticks([])
  805. axl_topup.set_yticks([])
  806. axl_topup.set_xlabel("Correlation to standard: r = 0.04", labelpad=6)
  807. axl_topup.set_title("TOPUP: Static distortion correction", pad=6)
  808. axl_topup.text(0.5, 0.5, "Participant 1", ha="center", va="center", transform=axl_topup.transAxes)
  809. axl_group.imshow(group_template_abcd)
  810. draw_seed(axl_group, x=0.48, y=0.7)
  811. axl_group.set_xticks([])
  812. axl_group.set_yticks([])
  813. axl_group.set_title("Group-averaged standard (TOPUP: static)", pad=6)
  814. axl_group.text(0.5, 0.5, "ABCD (N = 3,928)", ha="center", va="center", transform=axl_group.transAxes)
  815. axl_pos = axl_medic.get_position()
  816. f.text(
  817. 0.01,
  818. axl_pos.y1 + 0.075,
  819. r"$\bf{a}$ Functional Connectivity (FC) seed maps: Occipital Cortex",
  820. ha="left",
  821. va="center",
  822. fontsize=7,
  823. )
  825. mpl.rcParams["axes.edgecolor"] = "black"
  826. mpl.rcParams["xtick.labelsize"] = 5
  827. mpl.rcParams["ytick.labelsize"] = 5
  828. cbar_ax = f.add_subplot(gs_cbar[1, 1])
  829. pl = cbar_ax.imshow(
  830. np.array([[-0.5, 0.5], [0.5, -0.5]]),
  831. vmin=-0.5,
  832. vmax=0.5,
  833. aspect="auto",
  834. cmap=nilearn_cmaps["roy_big_bl"],
  835. )
  836. cbar = f.colorbar(
  837. mappable=pl,
  838. cax=cbar_ax,
  839. location="left",
  840. orientation="vertical",
  841. ticks=[-0.5, -0.25, 0, 0.25, 0.5],
  842. )
  843. cbar.ax.yaxis.set_ticks_position("right")
  844. cbar.ax.invert_yaxis()
  845. # create axis for colorbar
  846. cbar.ax.set_ylabel("Functional Connectivity z(r)", labelpad=2)
  848. # plot group similarities
  849. settings = GLOBAL_SETTINGS.copy()
  850. settings["style"] = "darkgrid"
  851. sns.set_theme(**settings)
  852. mpl.rcParams["xtick.major.size"] = 0.5
  853. mpl.rcParams["xtick.major.pad"] = 2
  854. mpl.rcParams["ytick.major.size"] = 0.5
  855. mpl.rcParams["font.size"] = 6
  856. ax1 = f.add_subplot(gs_bot[0])
  857. clr_pal = sns.color_palette("pastel")
  858. sns.scatterplot(topup_better_similarities_topup, medic_better_similarities_topup, s=14, ax=ax1)
  859. sns.scatterplot(topup_better_similarities_medic, medic_better_similarities_medic, s=14, ax=ax1)
  860. sns.scatterplot([0.51], [0.33], s=14, linewidth=0.5, ax=ax1, color=clr_pal[0])
  861. sns.scatterplot([0.32], [0.52], s=14, linewidth=0.5, ax=ax1, color=clr_pal[1])
  862. ax1.text(0.515, 0.33, "Scan", ha="left", va="center", fontsize=5, transform=ax1.transData)
  863. ax1.text(0.325, 0.52, "Scan", ha="left", va="center", fontsize=5, transform=ax1.transData)
  864. ax1.axline((0, 0), slope=1, color="black", linestyle="--", linewidth=1)
  865. ax1.set_xlabel("Correlation (r)", labelpad=-6)
  866. ax1.set_ylabel("Correlation (r)", labelpad=-10)
  867. ax1.set_aspect("equal")
  868. vmax = 0.55
  869. vmin = 0.3
  870. ax1.set_xlim([vmin, vmax])
  871. ax1.set_ylim([vmin, vmax])
  872. ax1.set_xticks(np.arange(vmin, vmax, 0.05))
  873. ax1.set_yticks(np.arange(vmin, vmax, 0.05))
  874. ax1.set_xticklabels(
  875. [f"{x:.2f}" if np.isclose(x, vmin) or np.isclose(x, vmax) else "" for x in np.arange(vmin, vmax, 0.05)]
  876. )
  877. ax1.set_yticklabels([f"{x:.2f}" if np.isclose(x, vmax) else "" for x in np.arange(vmin, vmax, 0.05)])
  878. ax1.text(
  879. 0.075,
  880. 0.95,
  881. "MEDIC more similar to Group Average",
  882. ha="left",
  883. va="center",
  884. transform=ax1.transAxes,
  885. )
  886. ax1.text(
  887. 0.925,
  888. 0.05,
  889. "TOPUP more similar to Group Average",
  890. ha="right",
  891. va="center",
  892. transform=ax1.transAxes,
  893. )
  894. x = np.array([-0.1, 0.7])
  895. y = x
  896. y2 = np.ones(x.shape) * 0.6
  897. y3 = np.zeros(x.shape)
  898. colors = sns.color_palette("pastel")
  899. medic_color = colors[1]
  900. topup_color = colors[0]
  901. ax1.fill_between(x, y, y2, color=medic_color, alpha=0.2)
  902. ax1.fill_between(x, y3, y, color=topup_color, alpha=0.2)
  903. ax1_pos = ax1.get_position()
  904. f.text(
  905. 0.01,
  906. ax1_pos.y1 + 0.075,
  907. r"$\bf{b}$ Whole-brain FC similarity to group-averaged standard",
  908. ha="left",
  909. va="center",
  910. fontsize=7,
  911. )
  912. sns.set_theme(**GLOBAL_SETTINGS)
  914. # plot t-statistic surface
  915. tstat_surface_path = DATA_DIR / "group_tstat_surface.png"
  916. tstat_surface = np.array(Image.open(tstat_surface_path))
  917. tstat_surface_left = tstat_surface[:, clip_1 : tstat_surface.shape[1] // 2 - clip_2]
  918. tstat_surface_right = tstat_surface[:, clip_2 + tstat_surface.shape[1] // 2 : -clip_1]
  919. tstat_surface = np.concatenate((tstat_surface_left, tstat_surface_right), axis=1)
  920. mpl.rcParams["axes.edgecolor"] = "white"
  921. mpl.rcParams["xtick.labelsize"] = 5
  922. mpl.rcParams["ytick.labelsize"] = 5
  923. ax2 = f.add_subplot(gs_tstat[0])
  924. ax2.imshow(tstat_surface)
  925. ax2.set_xticks([])
  926. ax2.set_yticks([])
  927. ax2.text(0.5, 0.5, "N = 185 Scans", ha="center", va="center", transform=ax2.transAxes, fontsize=6)
  928. mpl.rcParams["axes.edgecolor"] = "black"
  929. cbar_ax = f.add_subplot(gs_tstat_cbar[1, 1])
  930. spectral_map = plt.cm.get_cmap("Spectral")
  931. spectral_rmap = spectral_map.reversed()
  932. pl = cbar_ax.imshow(np.array([[-6, 6], [6, -6]]), vmin=-6, vmax=6, aspect="auto", cmap=spectral_rmap)
  933. cbar = f.colorbar(
  934. mappable=pl,
  935. cax=cbar_ax,
  936. location="left",
  937. orientation="vertical",
  938. ticks=[-6, -3, 0, 3, 6],
  939. )
  940. cbar.ax.yaxis.set_ticks_position("right")
  941. cbar.ax.invert_yaxis()
  942. # create axis for colorbar
  943. cbar.ax.set_ylabel("t-statistic", labelpad=2)
  944. ax2_pos = ax2.get_position()
  945. f.text(
  946. ax2_pos.x0,
  947. ax1_pos.y1 + 0.075,
  948. r"$\bf{c}$ Whole-brain FC similarity to group-averaged standard",
  949. ha="left",
  950. va="center",
  951. )
  952. cbar.ax.text(
  953. 0.5,
  954. 1.075,
  955. "MEDIC > TOPUP",
  956. ha="center",
  957. va="center",
  958. transform=cbar.ax.transAxes,
  959. fontsize=6,
  960. )
  961. cbar.ax.text(
  962. 0.5,
  963. -0.075,
  964. "TOPUP > MEDIC",
  965. ha="center",
  966. va="center",
  967. transform=cbar.ax.transAxes,
  968. fontsize=6,
  969. )
  971. # # for computing correlations
  972. # # load dconn data for medic and topup
  973. # group_avg = nib.load("/data/nil-bluearc/GMT/Scott/ABCD/ABCD_4.5k_all.dconn.nii")
  974. # medic_dconn = nib.load(
  975. # AA_DATA_DIR
  976. # / "sub-20008"
  977. # / "ses-51692"
  978. # / "cifti_correlation"
  979. # / "sub-20008_b1_MNI152_T1_2mm_Swgt_norm_bpss_resid_LR_surf_subcort_32k_fsLR_brainstem_surfsmooth1.7_subcortsmooth1.7.dconn.nii" # noqa
  980. # )
  981. # topup_dconn = nib.load(
  982. # AA_DATA_DIR
  983. # / "sub-20008"
  984. # / "ses-51692wTOPUP"
  985. # / "cifti_correlation"
  986. # / "sub-20008_b1_MNI152_T1_2mm_Swgt_norm_bpss_resid_LR_surf_subcort_32k_fsLR_brainstem_surfsmooth1.7_subcortsmooth1.7.dconn.nii" # noqa
  987. # )
  988. # # get the lower triangle of the dconn data
  989. # print("Loading dconn data...")
  990. # group_avg_data = group_avg.dataobj[21891, :59412]
  991. # medic_dconn_data = medic_dconn.dataobj[21891, :59412]
  992. # topup_dconn_data = topup_dconn.dataobj[21891, :59412]
  993. # # remove nans
  994. # group_avg_data[np.isnan(group_avg_data)] = 0
  995. # medic_dconn_data[np.isnan(medic_dconn_data)] = 0
  996. # topup_dconn_data[np.isnan(topup_dconn_data)] = 0
  997. # print("Done.")
  998. # # compute correlations
  999. # print("Computing correlations...")
  1000. # medic_corr = pearsonr(group_avg_data, medic_dconn_data)
  1001. # print(medic_corr)
  1002. # topup_corr = pearsonr(group_avg_data, topup_dconn_data)
  1003. # print(topup_corr)
  1004. # print("Done.")
  1006. f.savefig(FIGURES_DIR / "group_template_compare.png", dpi=300)
  1007. current_dir = os.getcwd()
  1008. os.chdir(FIGURES_DIR)
  1009. Path("figure3.png").unlink(missing_ok=True)
  1010. Path("figure3.png").symlink_to("group_template_compare.png")
  1011. os.chdir(current_dir)
  1012. sns.set_theme(**GLOBAL_SETTINGS)
  1015. # figure 4
  1016. def fmap_comparison(data_dir):
  1017. mpl.rcParams["axes.titlesize"] = 7
  1018. mpl.rcParams["axes.labelsize"] = 7
  1019. # get data
  1020. data_path = Path(data_dir)
  1022. # load images
  1023. minn_example_medic_path = data_path / "UMinn_medic.png"
  1024. minn_example_medic = equalize_hist(np.array(Image.open(minn_example_medic_path)))
  1025. minn_example_topup_path = data_path / "UMinn_topup.png"
  1026. minn_example_topup = equalize_hist(np.array(Image.open(minn_example_topup_path)))
  1027. minn_example_fmap_path = data_path / "UMinn_fmap.png"
  1028. minn_example_fmap = np.array(Image.open(minn_example_fmap_path))
  1029. penn_example_medic_path = data_path / "Penn_medic.png"
  1030. penn_example_medic = equalize_hist(np.array(Image.open(penn_example_medic_path)))
  1031. penn_example_topup_path = data_path / "Penn_topup.png"
  1032. penn_example_topup = equalize_hist(np.array(Image.open(penn_example_topup_path)))
  1033. penn_example_fmap_path = data_path / "Penn_fmap.png"
  1034. penn_example_fmap = np.array(Image.open(penn_example_fmap_path))
  1035. washu_example_medic_path = data_path / "WashU_medic.png"
  1036. washu_example_medic = equalize_hist(np.array(Image.open(washu_example_medic_path)))
  1037. washu_example_topup_path = data_path / "WashU_topup.png"
  1038. washu_example_topup = equalize_hist(np.array(Image.open(washu_example_topup_path)))
  1039. washu_example_fmap_path = data_path / "WashU_fmap.png"
  1040. washu_example_fmap = np.array(Image.open(washu_example_fmap_path))
  1042. f = plt.figure(figsize=(180 * MM_TO_INCHES, 140 * MM_TO_INCHES), layout="constrained")
  1043. gs = GridSpec(
  1044. 4,
  1045. 3,
  1046. left=0.03,
  1047. right=0.97,
  1048. bottom=0.025,
  1049. top=0.96,
  1050. wspace=0.15,
  1051. hspace=0.125,
  1052. height_ratios=[7, 1, 7, 7],
  1053. )
  1054. gs_cbar = GridSpecFromSubplotSpec(
  1055. 3,
  1056. 3,
  1057. wspace=0,
  1058. hspace=0,
  1059. width_ratios=[1, 50, 1],
  1060. height_ratios=[11, 5, 11],
  1061. subplot_spec=gs[1, :],
  1062. )
  1063. ax_WashU_fmap = f.add_subplot(gs[0, 0])
  1064. ax_WashU_medic = f.add_subplot(gs[2, 0])
  1065. ax_WashU_topup = f.add_subplot(gs[3, 0])
  1066. ax_UMinn_fmap = f.add_subplot(gs[0, 1])
  1067. ax_UMinn_medic = f.add_subplot(gs[2, 1])
  1068. ax_UMinn_topup = f.add_subplot(gs[3, 1])
  1069. ax_Penn_fmap = f.add_subplot(gs[0, 2])
  1070. ax_Penn_medic = f.add_subplot(gs[2, 2])
  1071. ax_Penn_topup = f.add_subplot(gs[3, 2])
  1073. cbar_ax = f.add_subplot(gs_cbar[1, 1])
  1074. pl = cbar_ax.imshow(np.array([[-50, 50], [50, -50]]), vmin=-50, vmax=50, aspect="auto", cmap="icefire")
  1075. cbar = f.colorbar(
  1076. pl,
  1077. cax=cbar_ax,
  1078. location="top",
  1079. orientation="horizontal",
  1080. ticks=[-50, -25, 0, 25, 50],
  1081. )
  1082. # create axis for colorbar
  1083. cbar.ax.set_xlabel("Field map Difference (Hz)", labelpad=2)
  1084. alt_vmin, alt_vmax = hz_limits_to_mm(-50, 50)
  1085. cax = cbar.ax.twiny()
  1086. cbar.ax.xaxis.set_ticks_position("top")
  1087. cax.xaxis.set_ticks_position("bottom")
  1088. cax.set_xlim(alt_vmin, alt_vmax)
  1089. cbar.ax.invert_xaxis()
  1090. cax.invert_xaxis()
  1091. cax.xaxis.set_label_position("bottom")
  1092. cax.set_xlabel("Displacement difference (mm)", labelpad=1)
  1094. # plot images
  1095. ax_WashU_fmap.imshow(washu_example_fmap)
  1096. ax_WashU_fmap.set_title(r"$\bf{a}$ WashU Data", pad=6, loc="left")
  1097. ax_WashU_fmap.set_xticks([])
  1098. ax_WashU_fmap.set_yticks([])
  1099. ax_WashU_fmap.set_ylabel("Field Map Difference (MEDIC - TOPUP)", labelpad=4)
  1100. ax_WashU_medic.imshow(washu_example_medic)
  1101. ax_WashU_medic.set_xticks([])
  1102. ax_WashU_medic.set_yticks([])
  1103. ax_WashU_medic.set_ylabel("MEDIC", labelpad=4)
  1104. draw_arrow(
  1105. ax_WashU_medic,
  1106. data_to_ax(ax_WashU_medic, (1683, 1858)),
  1107. data_to_ax(ax_WashU_medic, (1457, 1691)),
  1108. )
  1109. draw_arrow(
  1110. ax_WashU_medic,
  1111. data_to_ax(ax_WashU_medic, (2045, 1620)),
  1112. data_to_ax(ax_WashU_medic, (1829, 1474)),
  1113. )
  1114. ax_WashU_topup.imshow(washu_example_topup)
  1115. ax_WashU_topup.set_xticks([])
  1116. ax_WashU_topup.set_yticks([])
  1117. ax_WashU_topup.set_ylabel("TOPUP", labelpad=4)
  1118. draw_arrow(
  1119. ax_WashU_topup,
  1120. data_to_ax(ax_WashU_topup, (1683, 1858)),
  1121. data_to_ax(ax_WashU_topup, (1457, 1691)),
  1122. )
  1123. draw_arrow(
  1124. ax_WashU_topup,
  1125. data_to_ax(ax_WashU_topup, (2045, 1620)),
  1126. data_to_ax(ax_WashU_topup, (1829, 1474)),
  1127. )
  1128. ax_UMinn_fmap.imshow(minn_example_fmap)
  1129. ax_UMinn_fmap.set_title(r"$\bf{b}$ UMinn Data", pad=6, loc="left")
  1130. ax_UMinn_fmap.set_xticks([])
  1131. ax_UMinn_fmap.set_yticks([])
  1132. ax_UMinn_medic.imshow(minn_example_medic)
  1133. draw_arrow(
  1134. ax_UMinn_medic,
  1135. data_to_ax(ax_UMinn_medic, (1129, 622)),
  1136. data_to_ax(ax_UMinn_medic, (898, 473)),
  1137. )
  1138. ax_UMinn_medic.set_xticks([])
  1139. ax_UMinn_medic.set_yticks([])
  1140. ax_UMinn_topup.imshow(minn_example_topup)
  1141. ax_UMinn_topup.set_xticks([])
  1142. ax_UMinn_topup.set_yticks([])
  1143. draw_arrow(
  1144. ax_UMinn_topup,
  1145. data_to_ax(ax_UMinn_topup, (1129, 622)),
  1146. data_to_ax(ax_UMinn_topup, (898, 473)),
  1147. )
  1148. ax_Penn_fmap.imshow(penn_example_fmap)
  1149. ax_Penn_fmap.set_title(r"$\bf{c}$ Penn Data", pad=6, loc="left")
  1150. ax_Penn_fmap.set_xticks([])
  1151. ax_Penn_fmap.set_yticks([])
  1152. ax_Penn_medic.imshow(penn_example_medic)
  1153. ax_Penn_medic.set_xticks([])
  1154. ax_Penn_medic.set_yticks([])
  1155. draw_arrow(ax_Penn_medic, data_to_ax(ax_Penn_medic, (817, 611)), data_to_ax(ax_Penn_medic, (940, 414)))
  1156. draw_arrow(
  1157. ax_Penn_medic,
  1158. data_to_ax(ax_Penn_medic, (1202, 390)),
  1159. data_to_ax(ax_Penn_medic, (1266, 164)),
  1160. )
  1161. ax_Penn_topup.imshow(penn_example_topup)
  1162. ax_Penn_topup.set_xticks([])
  1163. ax_Penn_topup.set_yticks([])
  1164. draw_arrow(ax_Penn_topup, data_to_ax(ax_Penn_topup, (817, 611)), data_to_ax(ax_Penn_topup, (940, 414)))
  1165. draw_arrow(
  1166. ax_Penn_topup,
  1167. data_to_ax(ax_Penn_topup, (1202, 390)),
  1168. data_to_ax(ax_Penn_topup, (1266, 164)),
  1169. )
  1171. # save figure
  1172. f.savefig(FIGURES_DIR / "fieldmap_comparison.png", dpi=300)
  1173. current_dir = os.getcwd()
  1174. os.chdir(FIGURES_DIR)
  1175. Path("figure4.png").unlink(missing_ok=True)
  1176. Path("figure4.png").symlink_to("fieldmap_comparison.png")
  1177. os.chdir(current_dir)
  1178. sns.set_theme(**GLOBAL_SETTINGS)
  1181. # figure 5
  1182. def spotlight_comparison(data):
  1183. mpl.rcParams["axes.titlesize"] = 7
  1184. mpl.rcParams["axes.labelsize"] = 7
  1185. # load t1 and t2 t stat maps
  1186. t1_tstat = nib.load(Path(data) / "local_corr_t1_tstat.nii.gz").get_fdata().squeeze()
  1187. t2_tstat = nib.load(Path(data) / "local_corr_t2_tstat.nii.gz").get_fdata().squeeze()
  1188. t1_atlas_exemplar = (
  1189. nib.load(AA_DATA_DIR / "sub-20008" / "T1" / "atlas" / "sub-20008_T1w_debias_avg_on_MNI152_T1_2mm.nii.gz")
  1190. .get_fdata()
  1191. .squeeze()
  1192. )
  1193. t2_atlas_exemplar = (
  1194. nib.load(AA_DATA_DIR / "sub-20008" / "T1" / "atlas" / "sub-20008_T2w_debias_avg_on_MNI152_T1_2mm.nii.gz")
  1195. .get_fdata()
  1196. .squeeze()
  1197. )
  1199. # choose slices to iterate over
  1200. slices = np.linspace(16, t1_tstat.shape[2] - 20, 9).astype(int)[::-1]
  1202. # create figure
  1203. f = plt.figure(figsize=(180 * MM_TO_INCHES, 100 * MM_TO_INCHES), layout="constrained")
  1205. # create gridspec
  1206. gs = GridSpec(3, 3, left=0, right=0.975, bottom=0.025, top=0.95, wspace=0.02, hspace=0.02)
  1208. # create subfigures
  1209. subfigs = f.subfigures(1, 3, width_ratios=[1, 5, 5], wspace=0.03)
  1211. cgs = GridSpec(1, 1, left=0.52, right=0.57, bottom=0.073, top=0.89)
  1212. cbar_ax = subfigs[0].add_subplot(cgs[:, :])
  1214. # create axes for each subfigure
  1215. axes_list1 = []
  1216. for i in range(3):
  1217. for j in range(3):
  1218. axes_list1.append(subfigs[1].add_subplot(gs[i, j]))
  1219. axes_list2 = []
  1220. for i in range(3):
  1221. for j in range(3):
  1222. axes_list2.append(subfigs[2].add_subplot(gs[i, j]))
  1224. # create new cmap
  1225. new_cmap = sns.diverging_palette(210, 30, l=70, center="dark", as_cmap=True)
  1227. # plot slices, iterate over list
  1228. for i, s in enumerate(slices):
  1229. ax1 = axes_list1[i]
  1230. ax1.imshow(t1_atlas_exemplar[..., s].T, cmap="gray", origin="lower")
  1231. a = ax1.imshow(t1_tstat[..., s].T, cmap=new_cmap, vmin=-10, vmax=10, origin="lower", alpha=0.75)
  1232. ax1.set_xticks([])
  1233. ax1.set_yticks([])
  1234. if i == 0:
  1235. source_plot = a
  1236. ax2 = axes_list2[i]
  1237. ax2.imshow(t2_atlas_exemplar[..., s].T, cmap="gray", origin="lower")
  1238. ax2.imshow(t2_tstat[..., s].T, cmap=new_cmap, vmin=-10, vmax=10, origin="lower", alpha=0.75)
  1239. ax2.set_xticks([])
  1240. ax2.set_yticks([])
  1241. axes_list1[0].set_title(r"$\bf{a}$ T1w alignment: Spotlight analysis", pad=4, loc="left")
  1242. axes_list2[0].set_title(r"$\bf{b}$ T2w alignment: Spotlight analysis", pad=4, loc="left")
  1243. # create axis for colorbar
  1244. cbar = subfigs[0].colorbar(
  1245. source_plot,
  1246. cax=cbar_ax,
  1247. location="left",
  1248. orientation="vertical",
  1249. )
  1250. cbar.ax.yaxis.set_label_position("right")
  1251. cbar.ax.set_ylabel("t-statistic", labelpad=4)
  1252. cbar.ax.text(
  1253. 0.5,
  1254. 1.05,
  1255. "MEDIC > TOPUP",
  1256. ha="center",
  1257. va="center",
  1258. fontsize=6,
  1259. transform=cbar.ax.transAxes,
  1260. )
  1261. cbar.ax.text(
  1262. 0.5,
  1263. -0.05,
  1264. "TOPUP > MEDIC",
  1265. ha="center",
  1266. va="center",
  1267. fontsize=6,
  1268. transform=cbar.ax.transAxes,
  1269. )
  1271. f.savefig(FIGURES_DIR / "spotlight_comparison.png", dpi=300)
  1272. current_dir = os.getcwd()
  1273. os.chdir(FIGURES_DIR)
  1274. Path("figure5.png").unlink(missing_ok=True)
  1275. Path("figure5.png").symlink_to("spotlight_comparison.png")
  1276. os.chdir(current_dir)
  1277. sns.set_theme(**GLOBAL_SETTINGS)
  1280. # figure 6
  1281. def alignment_metrics(data):
  1282. # plot stats
  1283. data = pd.read_csv(data)
  1285. # create figures
  1286. settings = GLOBAL_SETTINGS.copy()
  1287. settings["style"] = "darkgrid"
  1288. sns.set_theme(**settings)
  1289. mpl.rcParams["axes.labelsize"] = 6
  1290. mpl.rcParams["axes.labelpad"] = 2
  1291. mpl.rcParams["xtick.labelsize"] = 5
  1292. mpl.rcParams["xtick.major.pad"] = 2
  1293. mpl.rcParams["ytick.labelsize"] = 5
  1294. mpl.rcParams["ytick.major.pad"] = 2
  1295. fig = plt.figure(figsize=(180 * MM_TO_INCHES, 90 * MM_TO_INCHES), layout="constrained")
  1296. subfigs = fig.subfigures(1, 2, width_ratios=[2, 1], wspace=0.05)
  1297. subfigs2 = subfigs[0].subfigures(2, 1, hspace=0.05, height_ratios=[1, 3])
  1298. fig_local = subfigs2[0]
  1299. fig_local.suptitle(r"$\bf{a}$ Local Metrics", fontsize=7, ha="left", x=0.02, weight="normal")
  1300. axes_local = fig_local.subplots(1, 2)
  1301. fig_global = subfigs2[1]
  1302. fig_global.suptitle(r"$\bf{b}$ Global Metrics", fontsize=7, ha="left", x=0.02, weight="normal")
  1303. axes_global = fig_global.subplots(3, 2)
  1304. fig_roc = subfigs[1]
  1305. fig_roc.suptitle(r"$\bf{c}$ Segmentation Metrics", fontsize=7, ha="left", x=0.02, weight="normal")
  1306. axes_roc = fig_roc.subplots(4, 1)
  1308. # create list of metrics
  1309. metrics = [
  1310. "local_corr_mean_t1",
  1311. "local_corr_mean_t2",
  1312. "corr_t1",
  1313. "corr_t2",
  1314. "grad_corr_t1",
  1315. "grad_corr_t2",
  1316. "nmi_t1",
  1317. "nmi_t2",
  1318. "roc_gw",
  1319. "roc_ie",
  1320. "roc_vw",
  1321. "roc_cb_ie",
  1322. ]
  1324. # create list of titles
  1325. titles = [
  1326. "T1w R$^2$ Spotlight",
  1327. "T2w R$^2$ Spotlight",
  1328. "T1w R$^2$",
  1329. "T2w R$^2$",
  1330. "T1w Grad. Correlation",
  1331. "T2w Grad. Correlation",
  1332. "T1w NMI",
  1333. "T2w NMI",
  1334. "Gray/White Matter AUC",
  1335. "Brain/Exterior AUC",
  1336. "Ventricles/White Matter AUC",
  1337. "Cerebellum/Exterior AUC",
  1338. ]
  1339. formatted_titles = [
  1340. "T1w R$^2$\nSpotlight",
  1341. "T2w R$^2$\nSpotlight",
  1342. "T1w R$^2$",
  1343. "T2w R$^2$",
  1344. "T1w Grad.\nCorrelation",
  1345. "T2w Grad.\nCorrelation",
  1346. "T1w NMI",
  1347. "T2w NMI",
  1348. "Gray/White\nMatter AUC",
  1349. "Brain/Exterior\nAUC",
  1350. "Ventricles/White\nMatter AUC",
  1351. "Cerebellum/Exterior\nAUC",
  1352. ]
  1354. # create list of axes for plotting
  1355. axes_list = [
  1356. axes_local[0],
  1357. axes_local[1],
  1358. axes_global[0][0],
  1359. axes_global[0][1],
  1360. axes_global[1][0],
  1361. axes_global[1][1],
  1362. axes_global[2][0],
  1363. axes_global[2][1],
  1364. axes_roc[0],
  1365. axes_roc[1],
  1366. axes_roc[2],
  1367. axes_roc[3],
  1368. ]
  1370. # plot box plots
  1371. for m, t, a in zip(metrics, formatted_titles, axes_list):
  1372. plot_box_plot(data, m, t, a)
  1374. # print ttest results
  1375. table_str = "| Metric | MEDIC | TOPUP | t-statistic | p-value | df |\n"
  1376. table_str += "| ------ | ----- | ----- | ----------- | ------- | -- |\n"
  1377. for m, t in zip(metrics, titles):
  1378. medic_data = data[f"{m}_medic"]
  1379. topup_data = data[f"{m}_topup"]
  1380. res = ttest_rel(medic_data, topup_data)
  1381. # round stats
  1382. medic_mean = np.round(medic_data.mean(), 3)
  1383. medic_std = np.round(medic_data.std(), 3)
  1384. topup_mean = np.round(topup_data.mean(), 3)
  1385. topup_std = np.round(topup_data.std(), 3)
  1386. p_value = np.round(res.pvalue, 3)
  1387. p_value = p_value if p_value >= 0.001 else "<0.001"
  1388. t_stat = np.round(res.statistic, 3)
  1389. print(f"{t}:")
  1390. print(f"MEDIC={medic_mean} ({medic_std}); " f"TOPUP={topup_mean} ({topup_std}); " f"p={p_value}; t={t_stat}\n")
  1391. table_str += f"| {t} | {medic_mean} ({medic_std}) | {topup_mean} ({topup_std}) |"
  1392. table_str += f" {t_stat} | {p_value} | {data.shape[0] - 1} |\n"
  1393. print(table_str)
  1395. # save figure
  1396. fig.savefig(FIGURES_DIR / "alignment_metrics.png", dpi=300)
  1397. current_dir = os.getcwd()
  1398. os.chdir(FIGURES_DIR)
  1399. Path("figure6.png").unlink(missing_ok=True)
  1400. Path("figure6.png").symlink_to("alignment_metrics.png")
  1401. os.chdir(current_dir)
  1402. sns.set_theme(**GLOBAL_SETTINGS)
  1405. # figure 7
  1406. def resp_analysis(data):
  1407. settings = GLOBAL_SETTINGS.copy()
  1408. settings["style"] = "darkgrid"
  1409. sns.set_theme(**settings)
  1410. mpl.rcParams["axes.labelsize"] = 6
  1411. mpl.rcParams["axes.labelpad"] = 2
  1412. mpl.rcParams["xtick.labelsize"] = 6
  1413. mpl.rcParams["xtick.major.pad"] = 2
  1414. mpl.rcParams["ytick.labelsize"] = 6
  1415. mpl.rcParams["ytick.major.pad"] = 2
  1416. mpl.rcParams["legend.fontsize"] = 6
  1417. mpl.rcParams["legend.title_fontsize"] = 6
  1419. # load data
  1420. ps_csv = sorted(Path(data).glob("power_spectra_run_*.csv"))
  1421. r_csv = sorted(Path(data).glob("resp_data_run_*.csv"))
  1422. power_spectra = [pd.read_csv(f).set_index("Frequency (Hz)") for f in ps_csv]
  1423. resp_data = [pd.read_csv(f).set_index("VOLUME") for f in r_csv]
  1425. # create figure
  1426. f = plt.figure(figsize=(180 * MM_TO_INCHES, 140 * MM_TO_INCHES), layout="constrained")
  1427. n_rows = len(power_spectra)
  1428. gs = [
  1429. GridSpec(
  1430. 2,
  1431. 2,
  1432. left=0.075,
  1433. right=0.975,
  1434. bottom=0.667 + 0.05,
  1435. top=1.000 - 0.04,
  1436. wspace=0.2,
  1437. hspace=0.2,
  1438. ),
  1439. GridSpec(
  1440. 2,
  1441. 2,
  1442. left=0.075,
  1443. right=0.975,
  1444. bottom=0.333 + 0.05,
  1445. top=0.667 - 0.04,
  1446. wspace=0.2,
  1447. hspace=0.2,
  1448. ),
  1449. GridSpec(
  1450. 2,
  1451. 2,
  1452. left=0.075,
  1453. right=0.975,
  1454. bottom=0.000 + 0.05,
  1455. top=0.333 - 0.04,
  1456. wspace=0.2,
  1457. hspace=0.2,
  1458. ),
  1459. ]
  1460. pastel = sns.color_palette("pastel")
  1461. palette = [pastel[2], pastel[4]]
  1462. ypos = [1.000 - 0.015, 0.667 - 0.015, 0.333 - 0.015]
  1463. for r, l in zip(range(n_rows), ["a", "b", "c"]):
  1464. ps = power_spectra[r]
  1465. rd = resp_data[r]
  1466. if r == 0:
  1467. f.text(
  1468. 0.05,
  1469. ypos[r],
  1470. r"$\bf{a}$ " + f"Spectral Power Density",
  1471. ha="left",
  1472. va="center",
  1473. fontsize=7,
  1474. )
  1475. f.text(
  1476. 0.55,
  1477. ypos[r],
  1478. r"$\bf{b}$ " + f"Respiratory Signal",
  1479. ha="left",
  1480. va="center",
  1481. fontsize=7,
  1482. )
  1484. # plot power spectra
  1485. ax1 = f.add_subplot(gs[r][0, 0])
  1486. sns.lineplot(data=ps["resp_signal"], linewidth=0.6, color=pastel[3], ax=ax1)
  1487. ax1.set_ylim(0, 60)
  1488. ax1.set_xlabel("")
  1489. ax1.set_ylabel("Respiratory Belt\nPower Density", labelpad=2)
  1490. ax1.set_xticklabels([])
  1491. ax1.tick_params(axis="x")
  1492. ax1.tick_params(axis="y")
  1493. ax1.set_title(f"Run {r + 1}")
  1494. ax2 = f.add_subplot(gs[r][1, 0])
  1495. sns.lineplot(
  1496. data=ps[["fmap_signal", "fmap_signal_filtered"]],
  1497. dashes=False,
  1498. linewidth=0.6,
  1499. palette=palette,
  1500. ax=ax2,
  1501. )
  1502. ax2.set_ylim(0, 60)
  1503. ax2.legend(["Unfiltered", "Filtered"], loc="upper center")
  1504. ax2.set_ylabel("MEDIC Field Map\nPower Density", labelpad=2)
  1505. ax2.set_xlabel("Frequency (Hz)", labelpad=2)
  1506. ax2.tick_params(axis="x")
  1507. ax2.tick_params(axis="y")
  1509. # plot resp signals
  1510. tr = 1.761
  1511. rd.index = rd.index * tr
  1512. corr = np.corrcoef(rd["resp_signal"], rd["fmap_signal"])[0, 1]
  1513. ax3 = f.add_subplot(gs[r][0, 1])
  1514. sns.lineplot(data=rd["resp_signal"], linewidth=0.6, color=pastel[3], ax=ax3)
  1515. ax3.set_xlabel("")
  1516. ax3.set_ylabel("Signal from\nRespiratory Belt", labelpad=2)
  1517. ax3.set_xticklabels([])
  1518. ax3.tick_params(axis="x")
  1519. ax3.tick_params(axis="y")
  1520. ax3.set_title(f"R = {np.round(corr, 3)}")
  1521. ax3.set_ylim(-3, 3)
  1522. ax4 = f.add_subplot(gs[r][1, 1])
  1523. sns.lineplot(data=rd["fmap_signal"], linewidth=0.6, color=palette[1], ax=ax4)
  1524. ax4.set_xlabel("Time (seconds)", labelpad=2)
  1525. ax4.set_ylabel("Signal from\nMEDIC Field Map", labelpad=2)
  1526. ax4.tick_params(axis="x")
  1527. ax4.tick_params(axis="y")
  1528. ax4.set_ylim(-3, 3)
  1529. f.savefig(FIGURES_DIR / "resp_analysis.png", dpi=300)
  1530. current_dir = os.getcwd()
  1531. os.chdir(FIGURES_DIR)
  1532. Path("figure7.png").unlink(missing_ok=True)
  1533. Path("figure7.png").symlink_to("resp_analysis.png")
  1534. os.chdir(current_dir)
  1535. sns.set_theme(**GLOBAL_SETTINGS)
  1538. # figure 10
  1539. def tsnr_comparision(data):
  1540. # load tsnralignment_metrics
  1541. tsnr_table = pd.read_csv(data)
  1542. settings = GLOBAL_SETTINGS.copy()
  1543. settings["style"] = "darkgrid"
  1544. sns.set_theme(**settings)
  1545. mpl.rcParams["axes.labelsize"] = 7
  1546. mpl.rcParams["axes.labelpad"] = 2
  1547. mpl.rcParams["xtick.labelsize"] = 7
  1548. mpl.rcParams["xtick.major.pad"] = 2
  1549. mpl.rcParams["ytick.labelsize"] = 7
  1550. mpl.rcParams["ytick.major.pad"] = 2
  1551. f = plt.figure(figsize=(90 * MM_TO_INCHES, 45 * MM_TO_INCHES), layout="constrained")
  1552. ax = f.add_subplot(1, 1, 1)
  1553. mpl.rcParams["xtick.labelsize"] = 6
  1554. mpl.rcParams["ytick.labelsize"] = 6
  1555. plot_box_plot(tsnr_table, "mean_tsnr_masked", "tSNR", ax)
  1556. ax.set_xlim([0, 160])
  1557. m = "mean_tsnr_masked"
  1558. t = "tSNR"
  1559. medic_data = tsnr_table[f"{m}_medic"]
  1560. topup_data = tsnr_table[f"{m}_topup"]
  1561. res = ttest_rel(medic_data, topup_data)
  1562. # round stats
  1563. medic_mean = np.round(medic_data.mean(), 3)
  1564. medic_std = np.round(medic_data.std(), 3)
  1565. topup_mean = np.round(topup_data.mean(), 3)
  1566. topup_std = np.round(topup_data.std(), 3)
  1567. p_value = np.round(res.pvalue, 3)
  1568. p_value = p_value if p_value >= 0.001 else "<0.001"
  1569. t_stat = np.round(res.statistic, 3)
  1570. print(f"{t}:")
  1571. print(f"MEDIC={medic_mean} ({medic_std}); " f"TOPUP={topup_mean} ({topup_std}); " f"p={p_value}; t={t_stat}\n")
  1572. f.savefig(FIGURES_DIR / "tsnr.png", dpi=300)
  1573. current_dir = os.getcwd()
  1574. os.chdir(FIGURES_DIR)
  1575. Path("figure10.png").unlink(missing_ok=True)
  1576. Path("figure10.png").symlink_to("tsnr.png")
  1577. os.chdir(current_dir)
  1578. sns.set_theme(**GLOBAL_SETTINGS)
  1581. # figure 100
  1582. def head_position_videos(data):
  1583. settings = GLOBAL_SETTINGS.copy()
  1584. settings["rc"].update(
  1585. {
  1586. "axes.facecolor": "black",
  1587. "figure.facecolor": "black",
  1588. "axes.labelcolor": "white",
  1589. "axes.titlecolor": "white",
  1590. "text.color": "white",
  1591. "xtick.color": "white",
  1592. "ytick.color": "white",
  1593. }
  1594. )
  1595. sns.set_theme(**settings)
  1596. # load field map files
  1597. medic_fieldmaps = Path(data) / "fieldmaps" / "medic_aligned"
  1599. transient_head_position_run_idx = [7, 8, 9, 10, 11, 12, 13]
  1601. # load transient field map runs
  1602. transient_fieldmaps = []
  1603. for idx in transient_head_position_run_idx:
  1604. run = idx + 1
  1605. transient_fieldmaps.append(nib.load(medic_fieldmaps / f"run{run:02d}" / "fmap.nii.gz"))
  1607. labels = [
  1608. "Neutral",
  1609. "+Z Rotation",
  1610. "-Z Rotation",
  1611. "+X Rotation",
  1612. "-X Rotation",
  1613. "+Y Rotation",
  1614. "-Y Rotation",
  1615. "Neutral to +Z Rotation",
  1616. "Neutral to -Z Rotation",
  1617. "Neutral to +X Rotation",
  1618. "Neutral to -X Rotation",
  1619. "Neutral to +Y Rotation",
  1620. "Neutral to -Y Rotation",
  1621. "Neutral to -Z Translation",
  1622. "-Z Translation",
  1623. ]
  1625. # get labels
  1626. labels = [labels[i] for i in transient_head_position_run_idx]
  1627. # replace space with underscores
  1628. labels = [label.replace(" ", "_") for label in labels]
  1630. # make transients output directory
  1631. transients_out = Path(FIGURES_DIR) / "videos" / "transients"
  1632. transients_out.mkdir(parents=True, exist_ok=True)
  1634. # load motion parameters
  1635. motion_params = []
  1636. for idx in transient_head_position_run_idx:
  1637. run = idx + 1
  1638. motion_params.append(
  1639. np.loadtxt(Path(data) / "framewise_align" / "func" / f"run{run:02d}" / f"run{run:02d}.par")
  1640. )
  1641. motion_params[-1][:, :3] = np.rad2deg(motion_params[-1][:, :3])
  1643. # render transient field map videos
  1644. def set_moco_label(motion_params):
  1645. def set_figure_labels(fig, frame_num):
  1646. # set label on figure
  1647. fig.text(
  1648. 0.5,
  1649. 0.7,
  1650. f"Frame {frame_num}"
  1651. f"\nMotion Parameters:"
  1652. f"\nrot-x: {motion_params[frame_num, 0]:.2f} deg"
  1653. f"\nrot-y: {motion_params[frame_num, 1]:.2f} deg"
  1654. f"\nrot-z: {motion_params[frame_num, 2]:.2f} deg"
  1655. f"\ntx: {motion_params[frame_num, 3]:.2f} mm"
  1656. f"\nty: {motion_params[frame_num, 4]:.2f} mm"
  1657. f"\ntz: {motion_params[frame_num, 5]:.2f} mm",
  1658. ha="center",
  1659. )
  1661. # return figure
  1662. return fig
  1664. # return function
  1665. return set_figure_labels
  1667. gs0 = GridSpec(1, 1, left=0.05, right=0.1, bottom=0.05, top=0.95)
  1668. gs1 = GridSpec(
  1669. 1,
  1670. 3,
  1671. left=0.125,
  1672. right=0.95,
  1673. bottom=0.05,
  1674. top=0.95,
  1675. hspace=0.025,
  1676. wspace=0.025,
  1677. width_ratios=[72, 110, 110],
  1678. )
  1680. def fig_callback():
  1681. fig = plt.figure(figsize=(10, 6), layout="constrained")
  1682. axes_list = [fig.add_subplot(gs1[0, i]) for i in range(3)]
  1683. cbar_ax = fig.add_subplot(gs0[:, :])
  1684. cbar_ax.axis("off")
  1685. return fig, axes_list, cbar_ax
  1687. for fmap, moco, label in zip(transient_fieldmaps, motion_params, labels):
  1688. print(f"Processing {label}")
  1689. render_dynamic_figure(
  1690. str(transients_out / f"{label}.mp4"),
  1691. [fmap],
  1692. colorbar=True,
  1693. colorbar_aspect=60,
  1694. colorbar_pad=0,
  1695. colorbar_labelpad=-5,
  1696. colorbar_alt_range=True,
  1697. colorbar_alt_labelpad=0,
  1698. fraction=0.3,
  1699. vmin=-100,
  1700. vmax=150,
  1701. colormaps="icefire",
  1702. figure_fx=set_moco_label(moco),
  1703. fig_callback=fig_callback,
  1704. text_color="white",
  1705. )
  1706. sns.set_theme(**GLOBAL_SETTINGS)
  1709. def main():
  1710. parser = argparse.ArgumentParser(description="script for generating paper figures")
  1711. parser.add_argument("--figures", nargs="+", type=int, help="figures to generate, if not supplied will plot all")
  1712. parser.add_argument("--figure_1_data", default=FIGURE1_DATA, help="path to figure 1 data")
  1713. parser.add_argument("--figure_2_data", default=FIGURE2_DATA, help="path to figure 2 data")
  1714. parser.add_argument("--figure_3_data", default=FIGURE3_DATA, help="path to figure 3 data")
  1715. parser.add_argument("--figure_4_data", default=FIGURE4_DATA, help="path to figure 4 data")
  1716. parser.add_argument("--figure_5_data", default=FIGURE5_DATA, help="path to figure 5 data")
  1717. parser.add_argument("--figure_6_data", default=FIGURE6_DATA, help="path to figure 6 data")
  1718. parser.add_argument("--figure_7_data", default=FIGURE7_DATA, help="path to figure 7 data")
  1719. parser.add_argument("--figure_10_data", default=FIGURE10_DATA, help="path to figure 10 data")
  1720. parser.add_argument("--figure_100_data", default=FIGURE100_DATA, help="path to figure 100 data")
  1722. # get arguments
  1723. args = parser.parse_args()
  1725. if args.figures is None or 1 in args.figures:
  1726. head_position_fieldmap(args.figure_1_data)
  1728. if args.figures is None or 2 in args.figures:
  1729. head_concatenation(args.figure_2_data)
  1731. if args.figures is None or 3 in args.figures:
  1732. group_template_comparison(args.figure_3_data)
  1734. if args.figures is None or 4 in args.figures:
  1735. fmap_comparison(args.figure_4_data)
  1737. if args.figures is None or 5 in args.figures:
  1738. spotlight_comparison(args.figure_5_data)
  1740. if args.figures is None or 6 in args.figures:
  1741. alignment_metrics(args.figure_6_data)
  1743. if args.figures is None or 7 in args.figures:
  1744. resp_analysis(args.figure_7_data)
  1746. if args.figures is None or 10 in args.figures:
  1747. tsnr_comparision(args.figure_10_data)
  1749. if args.figures is not None and 100 in args.figures:
  1750. head_position_videos(args.figure_100_data)
  1752. plt.show()

paper_figures.py at commit cc2ea27, under MIT · at the source

Overview

Authors: Andrew N Van1,2, David F Montez2,3, Timothy O Laumann3, Philip N Cho2, Vahdeta Suljic2, Thomas Madison4,5, Noah J Baden2, Nadeshka Ramirez-Perez2, Kristen M Scheidter2, Julia S Monk2, Forrest I Whiting2, Babatunde Adeyemo2, Roselyne J Chauvin2, Samuel R Krimmel2, Athanasia Metoki2, Aishwarya Rajesh6, Jarod L Roland7, Taylor Salo8, Anxu Wang2,9, Kimberly B Weldon5
and 13 other authorsAristeidis Sotiras6,10, Joshua S Shimony6, Benjamin P Kay2, Steven M Nelson5,11, Brenden Tervo-Clemmens5,12, Scott A Marek6, Luca Vizioli13, Essa Yacoub13, Theodore D Satterthwaite8, Evan M Gordon6, Damien A Fair4,5,11, M Dylan Tisdall14, Nico UF Dosenbach1,2,6,15
ORCID iDs: Andrew N Van, Timothy O Laumann, Philip N Cho, Thomas Madison, Noah J Baden, Roselyne J Chauvin, Athanasia Metoki, Aishwarya Rajesh, Jarod L Roland, Taylor Salo, Kimberly B Weldon, Aristeidis Sotiras, Joshua S Shimony, Benjamin P Kay, Brenden Tervo-Clemmens, Scott A Marek, Theodore D Satterthwaite, Evan M Gordon, Damien A Fair, M Dylan Tisdall, Nico UF Dosenbach
15 affiliations
  1. Department of Biomedical Engineering, Washington University in St. Louis, St. Louis, MO, United States
  2. Department of Neurology, Washington University School of Medicine, St. Louis, MO, United States
  3. Department of Psychiatry, Washington University School of Medicine, St. Louis, MO, United States
  4. Institute of Child Development, University of Minnesota Medical School, Minneapolis, MN, United States
  5. Masonic Institute for the Developing Brain, University of Minnesota Medical School, Minneapolis, MN, United States
  6. Department of Radiology, Washington University School of Medicine, St. Louis, MO, United States
  7. Department of Neurosurgery, Washington University School of Medicine, St. Louis, MO, United States
  8. Lifespan Informatics and Neuroimaging Center (PennLINC), Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA, United States
  9. Division of Computation and Data Science, Washington University School of Medicine, St. Louis, MO, United States
  10. Institute for Informatics, Data Science & Biostatistics, Washington University School of Medicine, St. Louis, MO, United States
  11. Department of Pediatrics, University of Minnesota Medical School, Minneapolis, MN, United States
  12. Department of Psychiatry & Behavioral Sciences, University of Minnesota Medical School, Minneapolis, MN, United States
  13. Center for Magnetic Resonance Research, University of Minnesota Medical School, Minneapolis, MN, United States
  14. Department of Radiology, Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA, United States
  15. Department of Pediatrics, Washington University School of Medicine, St. Louis, MO, United States
Institutions: Washington University in St. Louis (United States); University of Minnesota (United States); University of Pennsylvania (United States)
Journal: Imaging neuroscience (Cambridge, Mass.), volume 4, article IMAG.a.1262
Dates: received 22 January 2024; accepted 22 April 2026; published online 29 May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1162/imag.a.1262 · PMID 42232073 · PMCID PMC13224312 · OpenAlex W4389154158
Open access: diamond, a free copy (OpenAlex)
Preprint: doi:10.1101/2023.11.28.568744
Status: code verified
Categories: fMRI (modality)
Methods: Connectivity, Statistics, fMRI & imaging
Keywords: distortion correction, fMRI, multi-echo
Topic: Advanced MRI Techniques and Applications (Radiology, Nuclear Medicine and Imaging, Medicine), according to OpenAlex
Funding: NIMH NIH HHS (R00 MH129616)
Citations: cited by 1 paper (Europe PMC); 62 references in the paper

Abstract

Functional MRI (fMRI) data are severely distorted by magnetic field (B0) inhomogeneities, which currently must be corrected using separately acquired field map data. However, changes in the head position of a participant across fMRI frames cause changes in the B0 field, preventing accurate correction of geometric distortions. Movement during field map acquisitions corrupts field maps, preventing distortion correction altogether. In this study, we use multi-echo (ME) fMRI data to dynamically sample and correct for magnetic field image distortions caused by head motion. Our distortion correction pipeline, MEDIC (Multi-Echo DIstortion Correction), leverages magnetic field inhomogeneity information found in the difference between echoes and uses it to correct for distortion on a frame-by-frame basis. Here, we demonstrate that MEDIC’s frame-wise distortion correction decreases the impact of head motion on resting-state functional connectivity (RSFC) maps and improves alignment to anatomy when compared with the prior gold standard approach (i.e., FSL TOPUP). Enhanced frame-wise distortion correction with MEDIC, without the requirement for field map collection, furthers the benefit of cutting-edge multi-echo fMRI imaging over single-echo fMRI.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 20 matches between paragraphs and lines of code.

vanandrew/warpkit

License: other
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 4ae7699d23d5d7404a9cc4386de2bd5d3f0335e7, 6 August 2026
Languages: Python (38), C/C++ (12), Shell (2), C++ (1)
Size: 103 files, 53 scripts
Software Heritage: archived
Found in: the text, “MEDIC is computationally efficient and open sour”
Holds: README, license file, CITATION.cff, environment (Dockerfile, pyproject.toml, uv.lock), tests, continuous integration
Not found: documentation
Tools: NiBabel (20 files), NumPy (15 files), SciPy (3 files), scikit-image (2 files)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
55 files
At the source: github.com/vanandrew/warpkit

DosenbachGreene/processing_pipeline

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 8be76287c68caf556ed96ba42d116c02a9b77c66, 20 October 2023
Languages: Python (21), Shell (10), C (2), C/C++ (1)
Size: 228 files, 34 scripts
Software Heritage: not archived
Found in: “Data and Code Availability”
Holds: README, license file, environment (Apptainer, compose.yml, Dockerfile, pyproject.toml, setup.cfg, setup.py), tests
Not found: CITATION.cff, continuous integration, documentation
Tools: FSL (4 files), NumPy (4 files), FreeSurfer (3 files), NiBabel (2 files), PyBIDS (2 files), Connectome Workbench (2 files), HeuDiConv (1 file), Matplotlib (1 file), pandas (1 file), SciPy (1 file), seaborn (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
36 files
At the source: github.com/DosenbachGreene/processing_pipeline

vanandrew/medic_analysis

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: cc2ea27a90085e50762a25820773851310da1e31, 21 January 2024
Languages: MATLAB (555), Python (15), C (8)
Size: 747 files, 578 scripts
Software Heritage: archived
Found in: “Data and Code Availability”
Holds: README, license file, environment (Dockerfile, pyproject.toml, setup.cfg, setup.py)
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: FieldTrip (175 files), Statistics and Machine Learning Toolbox (16 files), SPM (12 files), GIfTI library for MATLAB (10 files), NiBabel (9 files), NumPy (9 files), Connectome Workbench (8 files), FreeSurfer (5 files), SciPy (5 files), pandas (4 files), scikit-image (4 files), FSL (3 files), Matplotlib (3 files), Parallel Computing Toolbox (2 files), PyBIDS (2 files), seaborn (2 files), Image Processing Toolbox (1 file), Optimization Toolbox (1 file), Tools for NIfTI and ANALYZE image (MATLAB) (1 file), Nilearn (1 file), Pillow (1 file), pydicom (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
580 files
At the source: github.com/vanandrew/medic_analysis

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 3 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 665 scripts, each with its path and the digest of its content;
  • 20 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data and Code Availability

The implementation for MEDIC is found at https://github.com/vanandrew/warpkit. Code for the processing pipeline is found at https://github.com/DosenbachGreene/processing_pipeline. Code for data analysis and figure generation is found at https://github.com/vanandrew/medic_analysis. Data used in our study are found at https://dosenbachlab.wustl.edu/data.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 33 authors, 3 keywords, 1 funder, 62 references.

Cite

This paper

Van, A. N., Montez, D. F., Laumann, T. O., Cho, P. N., Suljic, V., Madison, T., Baden, N. J., Ramirez-Perez, N., Scheidter, K. M., Monk, J. S., Whiting, F. I., Adeyemo, B., Chauvin, R. J., Krimmel, S. R., Metoki, A., Rajesh, A., Roland, J. L., Salo, T., Wang, A., . . . Dosenbach, N. U. (2026). Frame-wise multi-echo distortion correction for superior functional MRI. Imaging neuroscience (Cambridge, Mass.), 4, IMAG.a.1262. https://doi.org/10.1162/imag.a.1262

BibTeX

@article{van2026frame,
author = {Van, Andrew N and Montez, David F and Laumann, Timothy O and Cho, Philip N and Suljic, Vahdeta and Madison, Thomas and Baden, Noah J and Ramirez-Perez, Nadeshka and Scheidter, Kristen M and Monk, Julia S and Whiting, Forrest I and Adeyemo, Babatunde and Chauvin, Roselyne J and Krimmel, Samuel R and Metoki, Athanasia and Rajesh, Aishwarya and Roland, Jarod L and Salo, Taylor and Wang, Anxu and Weldon, Kimberly B and Sotiras, Aristeidis and Shimony, Joshua S and Kay, Benjamin P and Nelson, Steven M and Tervo-Clemmens, Brenden and Marek, Scott A and Vizioli, Luca and Yacoub, Essa and Satterthwaite, Theodore D and Gordon, Evan M and Fair, Damien A and Tisdall, M Dylan and Dosenbach, Nico UF},
title = {{Frame-wise multi-echo distortion correction for superior functional MRI}},
journal = {Imaging neuroscience (Cambridge, Mass.)},
year = {2026},
month = may,
volume = {4},
pages = {IMAG.a.1262},
publisher = {MIT Press},
issn = {2837-6056},
doi = {10.1162/imag.a.1262},
url = {https://doi.org/10.1162/imag.a.1262},
pmid = {42232073},
pmcid = {PMC13224312}
}

RIS

TY - JOUR
AU - Van, Andrew N
AU - Montez, David F
AU - Laumann, Timothy O
AU - Cho, Philip N
AU - Suljic, Vahdeta
AU - Madison, Thomas
AU - Baden, Noah J
AU - Ramirez-Perez, Nadeshka
AU - Scheidter, Kristen M
AU - Monk, Julia S
AU - Whiting, Forrest I
AU - Adeyemo, Babatunde
AU - Chauvin, Roselyne J
AU - Krimmel, Samuel R
AU - Metoki, Athanasia
AU - Rajesh, Aishwarya
AU - Roland, Jarod L
AU - Salo, Taylor
AU - Wang, Anxu
AU - Weldon, Kimberly B
AU - Sotiras, Aristeidis
AU - Shimony, Joshua S
AU - Kay, Benjamin P
AU - Nelson, Steven M
AU - Tervo-Clemmens, Brenden
AU - Marek, Scott A
AU - Vizioli, Luca
AU - Yacoub, Essa
AU - Satterthwaite, Theodore D
AU - Gordon, Evan M
AU - Fair, Damien A
AU - Tisdall, M Dylan
AU - Dosenbach, Nico UF
TI - Frame-wise multi-echo distortion correction for superior functional MRI
T2 - Imaging neuroscience (Cambridge, Mass.)
J2 - Imaging Neurosci (Camb)
PY - 2026
DA - 2026/05/29
VL - 4
SP - IMAG.a.1262
SN - 2837-6056
PB - MIT Press
DO - 10.1162/imag.a.1262
UR - https://doi.org/10.1162/imag.a.1262
LA - en
ER -

CSL-JSON

{
"id": "10.1162/imag.a.1262",
"type": "article-journal",
"title": "Frame-wise multi-echo distortion correction for superior functional MRI",
"container-title": "Imaging neuroscience (Cambridge, Mass.)",
"author": [
{
"family": "Van",
"given": "Andrew N"
},
{
"family": "Montez",
"given": "David F"
},
{
"family": "Laumann",
"given": "Timothy O"
},
{
"family": "Cho",
"given": "Philip N"
},
{
"family": "Suljic",
"given": "Vahdeta"
},
{
"family": "Madison",
"given": "Thomas"
},
{
"family": "Baden",
"given": "Noah J"
},
{
"family": "Ramirez-Perez",
"given": "Nadeshka"
},
{
"family": "Scheidter",
"given": "Kristen M"
},
{
"family": "Monk",
"given": "Julia S"
},
{
"family": "Whiting",
"given": "Forrest I"
},
{
"family": "Adeyemo",
"given": "Babatunde"
},
{
"family": "Chauvin",
"given": "Roselyne J"
},
{
"family": "Krimmel",
"given": "Samuel R"
},
{
"family": "Metoki",
"given": "Athanasia"
},
{
"family": "Rajesh",
"given": "Aishwarya"
},
{
"family": "Roland",
"given": "Jarod L"
},
{
"family": "Salo",
"given": "Taylor"
},
{
"family": "Wang",
"given": "Anxu"
},
{
"family": "Weldon",
"given": "Kimberly B"
},
{
"family": "Sotiras",
"given": "Aristeidis"
},
{
"family": "Shimony",
"given": "Joshua S"
},
{
"family": "Kay",
"given": "Benjamin P"
},
{
"family": "Nelson",
"given": "Steven M"
},
{
"family": "Tervo-Clemmens",
"given": "Brenden"
},
{
"family": "Marek",
"given": "Scott A"
},
{
"family": "Vizioli",
"given": "Luca"
},
{
"family": "Yacoub",
"given": "Essa"
},
{
"family": "Satterthwaite",
"given": "Theodore D"
},
{
"family": "Gordon",
"given": "Evan M"
},
{
"family": "Fair",
"given": "Damien A"
},
{
"family": "Tisdall",
"given": "M Dylan"
},
{
"family": "Dosenbach",
"given": "Nico UF"
}
],
"container-title-short": "Imaging Neurosci (Camb)",
"volume": "4",
"page": "IMAG.a.1262",
"DOI": "10.1162/imag.a.1262",
"PMID": "42232073",
"PMCID": "PMC13224312",
"ISSN": "2837-6056",
"publisher": "MIT Press",
"URL": "https://doi.org/10.1162/imag.a.1262",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
29
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1016/j.neuron.2026.04.011 [code]
Precision fMRI reveals densely interdigitated network patches with conserved motifs in the lateral prefrontal cortex.
Journal: Neuron
In common: Connectome Workbench, GIfTI library for MATLAB, Tools for NIfTI and ANALYZE image (MATLAB), 12 other tools, fMRI, 7 references
[2] doi:10.1162/imag.a.1222 [code]
Network-based near-scalp personalized brain stimulation targets.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Connectome Workbench, GIfTI library for MATLAB, Tools for NIfTI and ANALYZE image (MATLAB), 11 other tools, 6 references
[3] doi:10.1002/hbm.70621 [code]
Optimising 7T-fMRI for Imaging Regions of Magnetic Susceptibility.
Journal: Human brain mapping
In common: HeuDiConv, Tools for NIfTI and ANALYZE image (MATLAB), FieldTrip, 4 other tools, fMRI, 10 references
[4] doi:10.1038/s41467-026-73668-y [code]
Convergent and divergent brain-cognition development in early adolescence.
Journal: Nature communications
In common: Connectome Workbench, GIfTI library for MATLAB, Tools for NIfTI and ANALYZE image (MATLAB), 13 other tools, fMRI, 3 references
[5] doi:10.1038/s41586-026-10631-3 [code]
A prognostic human brain network for diffuse midline glioma.
Journal: Nature
In common: pydicom, GIfTI library for MATLAB, Tools for NIfTI and ANALYZE image (MATLAB), 15 other tools, 1 reference
[6] doi:10.1038/s41467-026-71151-2 [code]
Common and distinct neural correlates of social interaction processing and theory of mind in narratives.
Journal: Nature communications
In common: Connectome Workbench, GIfTI library for MATLAB, Optimization Toolbox, 15 other tools, 1 reference
[7] doi:10.1162/imag.a.1198 [code]
MEPrep: A robust pipeline for multi-echo fMRI denoising and preprocessing.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: PyBIDS, FreeSurfer, FSL, 8 other tools, fMRI, 8 references
[8] doi:10.1111/ene.70678 [code]
Who Falls After a Stroke? Evidence From a Prospective Stroke Cohort.
Journal: European journal of neurology
In common: pydicom, GIfTI library for MATLAB, Tools for NIfTI and ANALYZE image (MATLAB), 13 other tools, 2 references
[9] doi:10.1162/netn.a.547 [code]
An evaluation of the efficacy of single-echo and multi-echo fMRI denoising strategies.
Journal: Network neuroscience (Cambridge, Mass.)
In common: FreeSurfer, FSL, Nilearn, 6 other tools, fMRI, 10 references
[10] doi:10.1016/j.celrep.2026.117404 [code]
Action and rest tremor map to distinct networks within the primary motor cortex.
Journal: Cell reports
In common: pydicom, GIfTI library for MATLAB, Tools for NIfTI and ANALYZE image (MATLAB), 13 other tools, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.