OSCR

A Realistic In Silico Brain Phantom for Quantifying Susceptibility Anisotropy-Induced Error in Susceptibility Separation.

Code ↔ Paper

5 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 5 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 › Susceptibility Separation Phantom Validation Using In Vivo Data ↔ Manuscript_Figures.ipynb, lines 904–969 · score 0.67 · linear regression, simulated local field, local field map, slope, transformation, vivo
  2. [2] § Results › Comparison Between Simulated Versus In Vivo Susceptibility Maps ↔ Manuscript_Figures.ipynb, lines 904–969 · score 0.61 · vivo local field, linear regression, simulated local field, scatter, maps
  3. [3] § Methods › Phantom Creation › Susceptibility Maps ↔ PhantomCreation.m, lines 71–187 · score 0.56 · standard deviation, Gaussian noise, weighted, phantom, map, susceptibility
  4. [4] § Methods › Phantom Creation › Susceptibility Maps ↔ func/GRESimulation.m, the whole file · a weak match · score 0.56 · standard deviation, Gaussian noise, field, map
  5. [5] § Methods › Data Simulation › 3 T Simulations ↔ func/GRESimulation.m, the whole file · a weak match · score 0.54 · adding Gaussian noise, GRE, signal, TR, TE, phase

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

Jupyter notebook · 1,458 lines · 56 KB · MIT · 2 matches

  1. # %% [markdown]
  2. # # Install dependencies
  4. # %%
  5. import os
  6. import numpy as np
  7. import nibabel as nib
  8. import matplotlib.pyplot as plt
  9. from scipy.stats import linregress
  10. from pathlib import Path
  11. from scipy import stats
  12. from sklearn.linear_model import LinearRegression
  13. from scipy.stats import gaussian_kde
  14. from matplotlib.lines import Line2D
  15. from matplotlib.patches import Patch
  16. from pathlib import Path
  19. # %% [markdown]
  20. # # Download the data
  22. # %%
  23. base_dir = base_dir = Path.cwd()
  25. # Verify that the parent of the entered path exists
  26. if not base_dir.parent.exists():
  27. raise FileNotFoundError(f"The specified path {base_dir.parent} does not exist.")
  28. else:
  29. # Create the last folder if it doesn't exist
  30. base_dir.mkdir(parents=True, exist_ok=True)
  31. print(f"Using data directory: {base_dir}")
  32. !pip install osfclient
  33. # Clone the OSF project into the specified directory
  34. !osf -p 9xwhz clone "{base_dir}"
  37. # %% [markdown]
  38. # # Impact of susceptibility anisotropy
  40. # %%
  42. # ---------------------------
  43. # Parameters and File Paths
  44. # ---------------------------
  45. num_algorithms = 4
  46. algorithm_names = [
  47. '$\\chi$-separation',
  48. 'R2*-QSM',
  49. 'APART-QSM',
  50. 'DECOMPOSE-QSM'
  51. ]
  53. # Load segmentation mask (assumed to contain region labels 1 to 11)
  54. segmentation2_path = base_dir / 'osfstorage' / 'Masks' / 'white_matter_mask.nii.gz'
  55. segmentation2 = nib.load(segmentation2_path).get_fdata()
  56. region_labels = list(range(1, 12)) # regions 1 to 11
  58. # Load simulated maps (common for all algorithms)
  59. simulated_with_path = base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'Chi_negative_with_anisotropy.nii.gz'
  60. simulated_without_path = base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'Chi_negative.nii.gz'
  61. simulated_with = nib.load(simulated_with_path).get_fdata()
  62. simulated_without = nib.load(simulated_without_path).get_fdata()
  64. # Define measured maps for each algorithm
  65. measured_maps = {
  66. 0: { # Algorithm 1: χ-separation
  67. 'x_minus_with_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'X-separation' / 'With_anisotropy' / 'ChiNegMap.nii',
  68. 'x_minus_without_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'X-separation' / 'Without_anisotropy' / 'ChiNegMap.nii',
  69. },
  70. 1: { # Algorithm 2: R2*-QSM
  71. 'x_minus_with_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'R2star-QSM' / 'With_anisotropy' / 'ChiNegMap.nii',
  72. 'x_minus_without_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'R2star-QSM' / 'Without_anisotropy' / 'ChiNegMap.nii',
  73. },
  74. 2: { # Algorithm 3: APART-QSM
  75. 'x_minus_with_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'APART-QSM' / 'With_anisotropy' / 'X_dia_abs.nii',
  76. 'x_minus_without_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'APART-QSM' / 'Without_anisotropy' / 'X_dia_abs.nii',
  77. },
  78. 3: { # Algorithm 4: DECOMPOSE-QSM
  79. 'x_minus_with_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'DECOMPOSE-QSM' / 'With_anisotropy' / 'Results_DCS_abs.nii',
  80. 'x_minus_without_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'DECOMPOSE-QSM' / 'Without_anisotropy' / 'results_DCS_abs.nii',
  81. },
  82. }
  84. fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12, 12), dpi=300)
  85. t_test_results = {}
  86. for idx in range(num_algorithms):
  87. ax = axes[idx // 2, idx % 2]
  88. algorithm_name = algorithm_names[idx]
  90. measured_with = nib.load(measured_maps[idx]['x_minus_with_anisotropy']).get_fdata() * -1
  91. measured_without = nib.load(measured_maps[idx]['x_minus_without_anisotropy']).get_fdata() * -1
  93. errors_with_regions = []
  94. errors_without_regions = []
  96. for region in region_labels:
  97. region_mask = segmentation2 == region
  98. if np.any(region_mask):
  99. # For "With Anisotropy"
  100. sim_val_with = np.mean(simulated_with[region_mask])
  101. meas_val_with = np.mean(measured_with[region_mask])
  102. if sim_val_with != 0:
  103. error_with = ((meas_val_with - sim_val_with) / sim_val_with) ** 2 * 100
  104. else:
  105. error_with = np.nan
  107. # For "Without Anisotropy"
  108. sim_val_without = np.mean(simulated_without[region_mask])
  109. meas_val_without = np.mean(measured_without[region_mask])
  110. if sim_val_without != 0:
  111. error_without = ((meas_val_without - sim_val_without) / sim_val_without) ** 2 * 100
  112. else:
  113. error_without = np.nan
  115. errors_with_regions.append(error_with)
  116. errors_without_regions.append(error_without)
  118. errors_with_regions = np.array(errors_with_regions)
  119. errors_without_regions = np.array(errors_without_regions)
  120. errors_with_regions = errors_with_regions[~np.isnan(errors_with_regions)]
  121. errors_without_regions = errors_without_regions[~np.isnan(errors_without_regions)]
  123. x_vals = np.linspace(0, 70, 200)
  124. max_height = 0.4
  126. # Plot KDE for "With Anisotropy"
  127. if errors_with_regions.size > 1:
  128. kde_with = gaussian_kde(errors_with_regions)
  129. density_with = kde_with(x_vals)
  130. scaling_with = max_height / np.max(density_with) if np.max(density_with) > 0 else 1
  131. density_with_scaled = density_with * scaling_with
  132. ax.fill_between(x_vals, 0, density_with_scaled, color='tab:blue', alpha=0.4)
  133. ax.plot(x_vals, density_with_scaled, color='tab:blue', alpha=0.7)
  134. mean_mspe_with = np.mean(errors_with_regions)
  135. elif errors_with_regions.size == 1:
  136. mean_mspe_with = errors_with_regions[0]
  137. ax.plot([mean_mspe_with], [max_height/2], marker='o', color='tab:blue')
  138. else:
  139. mean_mspe_with = None # no data available
  141. # Plot KDE for "Without Anisotropy"
  142. if errors_without_regions.size > 1:
  143. kde_without = gaussian_kde(errors_without_regions)
  144. density_without = kde_without(x_vals)
  145. scaling_without = max_height / np.max(density_without) if np.max(density_without) > 0 else 1
  146. density_without_scaled = density_without * scaling_without
  147. ax.fill_between(x_vals, 0, density_without_scaled, color='tab:orange', alpha=0.4)
  148. ax.plot(x_vals, density_without_scaled, color='tab:orange', alpha=0.7)
  149. mean_mspe_without = np.mean(errors_without_regions)
  150. elif errors_without_regions.size == 1:
  151. mean_mspe_without = errors_without_regions[0]
  152. ax.plot([mean_mspe_without], [max_height/2], marker='o', color='tab:orange')
  153. else:
  154. mean_mspe_without = None
  156. if mean_mspe_with is not None:
  157. ax.text(0.95, 0.8, f"Mean MSPE: {mean_mspe_with:.1f}%", transform=ax.transAxes,
  158. color='tab:blue', ha='right', va='center', fontsize=12)
  159. if mean_mspe_without is not None:
  160. ax.text(0.95, 0.7, f"Mean MSPE: {mean_mspe_without:.1f}%", transform=ax.transAxes,
  161. color='tab:orange', ha='right', va='center', fontsize=12)
  163. # ------------------------------------------------
  164. # Create and add a legend for the subplot
  165. # ------------------------------------------------
  166. legend_elements = [
  167. Line2D([0], [0], color='tab:blue', lw=2, label='With Anisotropy'),
  168. Line2D([0], [0], color='tab:orange', lw=2, label='Without Anisotropy')
  169. ]
  170. ax.legend(handles=legend_elements, loc='upper right')
  172. # Set title and axis labels.
  173. ax.set_title(algorithm_name, fontsize=14)
  174. ax.set_xlabel("MSPE of $\\chi^-$ (%)", fontsize=12)
  175. ax.set_ylabel("Density (a.u.)", fontsize=12)
  176. ax.set_xlim(0, 70)
  177. ax.set_ylim(0, max_height * 1.2)
  178. ax.grid(alpha=0.3)
  180. # ---------------------------
  181. # Final figure adjustments
  182. # ---------------------------
  183. plt.tight_layout(rect=[0, 0, 1, 0.96])
  184. plt.show()
  187. # %% [markdown]
  188. # # Impact of noise
  190. # %%
  191. # -------------------------------------------
  192. # Define algorithms, SNR levels, and errors
  193. # -------------------------------------------
  194. algorithms = [
  195. '$\\chi$-separation',
  196. 'R2*-QSM',
  197. 'APART-QSM',
  198. 'DECOMPOSE-QSM'
  199. ]
  200. snr_levels = [50, 100, 200, 300]
  202. # Initialize dictionaries
  203. x_positive_errors = {alg: [] for alg in algorithms}
  204. x_negative_errors = {alg: [] for alg in algorithms}
  206. def construct_path(*args):
  207. """Construct a path with 'osfstorage' as a subdirectory of base_dir."""
  208. return base_dir / 'osfstorage' / Path(*args)
  210. # -------------------------------------------
  211. # Load Simulated Maps
  212. # -------------------------------------------
  213. simulated_positive_path = construct_path('Susceptibility_Separation_Results', 'Chi_positive.nii.gz')
  214. simulated_negative_path = construct_path('Susceptibility_Separation_Results', 'Chi_negative_with_anisotropy.nii.gz')
  215. simulated_positive = nib.load(simulated_positive_path).get_fdata()
  216. simulated_negative = nib.load(simulated_negative_path).get_fdata()
  218. # -------------------------------------------
  219. # Load Segmentation Masks
  220. # -------------------------------------------
  221. segmentation_positive_path = construct_path('Masks', 'SegmentedModel.nii.gz')
  222. segmentation_negative_path = construct_path('Masks', 'white_matter_mask.nii.gz')
  223. segmentation_positive = nib.load(segmentation_positive_path).get_fdata()
  224. segmentation_negative = nib.load(segmentation_negative_path).get_fdata()
  226. # -------------------------------------------
  227. # Measured Maps: Positive & Negative
  228. # -------------------------------------------
  229. measured_maps_positive = {
  230. '$\\chi$-separation': {
  231. 300: construct_path('Noise', 'X-separation', 'SNR_300', 'ChiPosMap.nii.gz'),
  232. 200: construct_path('Noise', 'X-separation', 'SNR_200', 'ChiPosMap.nii.gz'),
  233. 100: construct_path('Noise', 'X-separation', 'SNR_100', 'ChiPosMap.nii.gz'),
  234. 50: construct_path('Noise', 'X-separation', 'SNR_50', 'ChiPosMap.nii.gz'),
  235. },
  236. 'R2*-QSM': {
  237. 300: construct_path('Noise', 'R2star-QSM', 'SNR_300', 'ChiPosMap.nii.gz'),
  238. 200: construct_path('Noise', 'R2star-QSM', 'SNR_200', 'ChiPosMap.nii.gz'),
  239. 100: construct_path('Noise', 'R2star-QSM', 'SNR_100', 'ChiPosMap.nii.gz'),
  240. 50: construct_path('Noise', 'R2star-QSM', 'SNR_50', 'ChiPosMap.nii.gz'),
  241. },
  242. 'APART-QSM': {
  243. 300: construct_path('Noise', 'APART-QSM', 'SNR_300', 'X_para.nii.gz'),
  244. 200: construct_path('Noise', 'APART-QSM', 'SNR_200', 'X_para.nii.gz'),
  245. 100: construct_path('Noise', 'APART-QSM', 'SNR_100', 'X_para.nii.gz'),
  246. 50: construct_path('Noise', 'APART-QSM', 'SNR_50', 'X_para.nii.gz'),
  247. },
  248. 'DECOMPOSE-QSM': {
  249. 300: construct_path('Noise', 'DECOMPOSE-QSM', 'SNR_300', 'results_PCS.nii.gz'),
  250. 200: construct_path('Noise', 'DECOMPOSE-QSM', 'SNR_200', 'results_PCS.nii.gz'),
  251. 100: construct_path('Noise', 'DECOMPOSE-QSM', 'SNR_100', 'results_PCS.nii.gz'),
  252. 50: construct_path('Noise', 'DECOMPOSE-QSM', 'SNR_50', 'results_PCS.nii.gz'),
  253. },
  254. }
  256. measured_maps_negative = {
  257. '$\\chi$-separation': {
  258. 300: construct_path('Noise', 'X-separation', 'SNR_300', 'ChiNegMap.nii.gz'),
  259. 200: construct_path('Noise', 'X-separation', 'SNR_200', 'ChiNegMap.nii.gz'),
  260. 100: construct_path('Noise', 'X-separation', 'SNR_100', 'ChiNegMap.nii.gz'),
  261. 50: construct_path('Noise', 'X-separation', 'SNR_50', 'ChiNegMap.nii.gz'),
  262. },
  263. 'R2*-QSM': {
  264. 300: construct_path('Noise', 'R2star-QSM', 'SNR_300', 'ChiNegMap.nii.gz'),
  265. 200: construct_path('Noise', 'R2star-QSM', 'SNR_200', 'ChiNegMap.nii.gz'),
  266. 100: construct_path('Noise', 'R2star-QSM', 'SNR_100', 'ChiNegMap.nii.gz'),
  267. 50: construct_path('Noise', 'R2star-QSM', 'SNR_50', 'ChiNegMap.nii.gz'),
  268. },
  269. 'APART-QSM': {
  270. 300: construct_path('Noise', 'APART-QSM', 'SNR_300', 'X_dia_abs.nii.gz'),
  271. 200: construct_path('Noise', 'APART-QSM', 'SNR_200', 'X_dia_abs.nii.gz'),
  272. 100: construct_path('Noise', 'APART-QSM', 'SNR_100', 'X_dia_abs.nii.gz'),
  273. 50: construct_path('Noise', 'APART-QSM', 'SNR_50', 'X_dia_abs.nii.gz'),
  274. },
  275. 'DECOMPOSE-QSM': {
  276. 300: construct_path('Noise', 'DECOMPOSE-QSM', 'SNR_300', 'results_DCS_abs.nii.gz'),
  277. 200: construct_path('Noise', 'DECOMPOSE-QSM', 'SNR_200', 'results_DCS_abs.nii.gz'),
  278. 100: construct_path('Noise', 'DECOMPOSE-QSM', 'SNR_100', 'results_DCS_abs.nii.gz'),
  279. 50: construct_path('Noise', 'DECOMPOSE-QSM', 'SNR_50', 'results_DCS_abs.nii.gz'),
  280. },
  281. }
  283. # -------------------------------------------
  284. # Compute Errors
  285. # -------------------------------------------
  286. for alg in algorithms:
  287. for snr in snr_levels:
  288. # Load measured maps (X+ and X-)
  289. measured_positive_path = measured_maps_positive[alg][snr]
  290. measured_negative_path = measured_maps_negative[alg][snr]
  292. measured_positive = nib.load(measured_positive_path).get_fdata()
  293. measured_negative = nib.load(measured_negative_path).get_fdata()
  295. # Multiply X- by -1
  296. measured_negative = -1 * measured_negative
  298. # -----------------------------
  299. # Positive regions: 1 to 9
  300. # -----------------------------
  301. positive_errors = []
  302. for region in range(1, 10):
  303. mask = (segmentation_positive == region)
  304. measured_mean = np.mean(measured_positive[mask])
  305. simulated_mean = np.mean(simulated_positive[mask])
  307. error = ((measured_mean - simulated_mean) / simulated_mean)**2 * 100
  308. positive_errors.append(error)
  310. avg_positive_error = np.mean(positive_errors)
  311. x_positive_errors[alg].append(avg_positive_error)
  313. # -----------------------------
  314. # Negative regions
  315. # - 1 to 9 (segmentation_positive)
  316. # - 1 to 10 (segmentation_negative)
  317. # -----------------------------
  318. negative_errors = []
  320. # 1 to 9
  321. for region in range(1, 10):
  322. mask = (segmentation_positive == region)
  323. measured_mean = np.mean(measured_negative[mask])
  324. simulated_mean = np.mean(simulated_negative[mask])
  326. error = ((measured_mean - simulated_mean) / simulated_mean)**2 * 100
  327. negative_errors.append(error)
  329. # 1 to 10 (segmentation_negative)
  330. for region in range(1, 11):
  331. mask = (segmentation_negative == region)
  332. measured_mean = np.mean(measured_negative[mask])
  333. simulated_mean = np.mean(simulated_negative[mask])
  335. error = ((measured_mean - simulated_mean) / simulated_mean)**2 * 100
  336. negative_errors.append(error)
  338. avg_negative_error = np.mean(negative_errors)
  339. x_negative_errors[alg].append(avg_negative_error)
  341. # -------------------------------------------
  342. # Plotting
  343. # -------------------------------------------
  344. hatch_patterns = ['x', 'o', '*', '+']
  345. snr_colors = {
  346. 50: '#66c2a5',
  347. 100: '#8da0cb',
  348. 200: '#fc8d62',
  349. 300: '#e78ac3',
  350. }
  352. bar_width = 0.2
  353. snr_index = np.arange(len(snr_levels))
  355. fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6), dpi=600)
  357. tick_font_size = 14
  359. # --- Plot χ+ ---
  360. for i, (alg, hatch) in enumerate(zip(algorithms, hatch_patterns)):
  361. for j, snr in enumerate(snr_levels):
  362. x_coord = snr_index[j] + i * bar_width
  363. y_value = x_positive_errors[alg][j]
  366. ax1.bar(
  367. x_coord,
  368. y_value,
  369. bar_width,
  370. hatch=hatch,
  371. edgecolor='black',
  372. color=snr_colors[snr],
  373. )
  375. ax1.set_xlabel('SNR (a.u.)', fontsize=tick_font_size)
  376. ax1.set_ylabel('MSPE of $\\chi^+$ (%)', fontsize=tick_font_size)
  377. ax1.set_xticks(snr_index + bar_width * 1.5)
  378. ax1.set_xticklabels(snr_levels, fontsize=tick_font_size)
  379. ax1.grid(alpha=0.3)
  380. ax1.tick_params(axis='both', which='major', labelsize=tick_font_size)
  382. # --- Plot χ- ---
  383. for i, (alg, hatch) in enumerate(zip(algorithms, hatch_patterns)):
  384. for j, snr in enumerate(snr_levels):
  385. x_coord = snr_index[j] + i * bar_width
  386. y_value = x_negative_errors[alg][j]
  389. ax2.bar(
  390. x_coord,
  391. y_value,
  392. bar_width,
  393. hatch=hatch,
  394. edgecolor='black',
  395. color=snr_colors[snr],
  396. )
  398. ax2.set_xlabel('SNR (a.u.)', fontsize=tick_font_size)
  399. ax2.set_ylabel('MSPE of $\\chi^-$ (%)', fontsize=tick_font_size)
  400. ax2.set_xticks(snr_index + bar_width * 1.5)
  401. ax2.set_xticklabels(snr_levels, fontsize=tick_font_size)
  402. ax2.grid(alpha=0.3)
  403. ax2.tick_params(axis='both', which='major', labelsize=tick_font_size)
  406. # For χ+ axis:
  407. algorithm_patches_pos = []
  408. for alg, hatch in zip(algorithms, hatch_patterns):
  409. patch = Patch(
  410. facecolor='white',
  411. edgecolor='black',
  412. hatch=hatch,
  413. label=alg
  414. )
  415. algorithm_patches_pos.append(patch)
  417. ax1.legend(
  418. handles=algorithm_patches_pos,
  419. fontsize=14,
  420. facecolor='white',
  421. framealpha=1,
  422. loc='upper right',
  423. handlelength=1.5,
  424. handleheight=1.5
  425. )
  427. # For χ- axis:
  428. algorithm_patches_neg = []
  429. for alg, hatch in zip(algorithms, hatch_patterns):
  430. patch = Patch(
  431. facecolor='white',
  432. edgecolor='black',
  433. hatch=hatch,
  434. label=alg
  435. )
  436. algorithm_patches_neg.append(patch)
  438. ax2.legend(
  439. handles=algorithm_patches_neg,
  440. fontsize=14,
  441. facecolor='white',
  442. framealpha=1,
  443. loc='upper right',
  444. handlelength=1.5,
  445. handleheight=1.5
  446. )
  448. plt.tight_layout()
  449. plt.show()
  452. # %% [markdown]
  453. # # Impact of noise on susceptibility anisotorpy
  455. # %%
  456. # Define the algorithms and SNR levels
  457. algorithms = ['x-separation', 'R2*-QSM', 'APART-QSM', 'DECOMPOSE-QSM']
  458. snr_levels = ['SNR 300', 'SNR 200', 'SNR 100', 'SNR 50']
  461. # Initialize a dictionary to store the full paths for measured images
  462. measured_image_paths = {
  463. 'x-separation': {
  464. 'SNR 300': base_dir / 'osfstorage' / 'Noise' / 'X-separation' / 'SNR_300' / 'ChiNegMap.nii.gz',
  465. 'SNR 200': base_dir / 'osfstorage' / 'Noise' / 'X-separation' / 'SNR_200' / 'ChiNegMap.nii.gz',
  466. 'SNR 100': base_dir / 'osfstorage' / 'Noise' / 'X-separation' / 'SNR_100' / 'ChiNegMap.nii.gz',
  467. 'SNR 50': base_dir / 'osfstorage' / 'Noise' / 'X-separation' / 'SNR_50' / 'ChiNegMap.nii.gz',
  468. },
  469. 'R2*-QSM': {
  470. 'SNR 300': base_dir / 'osfstorage' / 'Noise' / 'R2star-QSM' / 'SNR_300' / 'ChiNegMap.nii.gz',
  471. 'SNR 200': base_dir / 'osfstorage' / 'Noise' / 'R2star-QSM' / 'SNR_200' / 'ChiNegMap.nii.gz',
  472. 'SNR 100': base_dir / 'osfstorage' / 'Noise' / 'R2star-QSM' / 'SNR_100' / 'ChiNegMap.nii.gz',
  473. 'SNR 50': base_dir / 'osfstorage' / 'Noise' / 'R2star-QSM' / 'SNR_50' / 'ChiNegMap.nii.gz',
  474. },
  475. 'APART-QSM': {
  476. 'SNR 300': base_dir / 'osfstorage' / 'Noise' / 'APART-QSM' / 'SNR_300' / 'X_dia_abs.nii.gz',
  477. 'SNR 200': base_dir / 'osfstorage' / 'Noise' / 'APART-QSM' / 'SNR_200' / 'X_dia_abs.nii.gz',
  478. 'SNR 100': base_dir / 'osfstorage' / 'Noise' / 'APART-QSM' / 'SNR_100' / 'X_dia_abs.nii.gz',
  479. 'SNR 50': base_dir / 'osfstorage' / 'Noise' / 'APART-QSM' / 'SNR_50' / 'X_dia_abs.nii.gz',
  480. },
  481. 'DECOMPOSE-QSM': {
  482. 'SNR 300': base_dir / 'osfstorage' / 'Noise' / 'DECOMPOSE-QSM' / 'SNR_300' / 'results_DCS_abs.nii.gz',
  483. 'SNR 200': base_dir / 'osfstorage' / 'Noise' / 'DECOMPOSE-QSM' / 'SNR_200' / 'results_DCS_abs.nii.gz',
  484. 'SNR 100': base_dir / 'osfstorage' / 'Noise' / 'DECOMPOSE-QSM' / 'SNR_100' / 'results_DCS_abs.nii.gz',
  485. 'SNR 50': base_dir / 'osfstorage' / 'Noise' / 'DECOMPOSE-QSM' / 'SNR_50' / 'results_DCS_abs.nii.gz',
  486. },
  487. }
  489. # ----------------------------
  490. # 2) LOAD NIFTI DATA
  491. # ----------------------------
  493. # Load simulated (ground truth) chi, theta, and segmentation images
  494. simulated_nii = nib.load(base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'Chi_negative_with_anisotropy.nii.gz')
  495. simulated_img_with_anisotropy = simulated_nii.get_fdata()
  497. theta_nii = nib.load(base_dir / 'osfstorage' / 'Noise' / 'theta.nii.gz')
  498. theta_img = theta_nii.get_fdata()
  499. theta_img[np.isnan(theta_img)] = 0
  501. segmentation_nii = nib.load(base_dir / 'osfstorage' / 'Masks' / 'WM_fibers_seg.nii.gz')
  502. segmentation_img = segmentation_nii.get_fdata()
  504. theta_shape = theta_nii.shape
  505. theta_affine = theta_nii.affine
  507. # Load measured images for each algorithm and SNR
  508. measured_images = {}
  509. for algorithm in algorithms:
  510. measured_images[algorithm] = {}
  511. for snr in snr_levels:
  512. path = measured_image_paths[algorithm][snr]
  513. measured_images[algorithm][snr] = nib.load(path).get_fdata()
  515. # ----------------------------
  516. # 3) GATHER VALID VOXELS & COMPUTE ERRORS
  517. # ----------------------------
  519. all_valid_coords = []
  520. all_theta_values = []
  521. all_sim_values = []
  523. errors = {
  524. algorithm: {
  525. snr: [] for snr in snr_levels
  526. }
  527. for algorithm in algorithms
  528. }
  530. # Loop over the segmentation regions
  531. for region in range(1, 28):
  532. # Region mask
  533. region_mask = (segmentation_img == region)
  534. region_valid_mask = region_mask & (theta_img != 0)
  536. # Get the 3D coordinates of those valid voxels
  537. valid_coords = np.argwhere(region_valid_mask)
  538. if valid_coords.size == 0:
  539. continue
  541. # Extract the actual voxel values
  542. region_theta_vals = theta_img[region_valid_mask]
  543. region_sim_vals = simulated_img_with_anisotropy[region_valid_mask] * -1
  545. # Append these to our global lists
  546. all_valid_coords.append(valid_coords)
  547. all_theta_values.append(region_theta_vals)
  548. all_sim_values.append(region_sim_vals)
  550. # For each algorithm & SNR, calculate MSPE in this region and append
  551. for algorithm in algorithms:
  552. for snr in snr_levels:
  553. measured_vals = measured_images[algorithm][snr][region_valid_mask]
  554. epsilon = 1e-10
  555. relative_squared_error = np.sqrt(
  556. np.square(region_sim_vals - measured_vals) /
  557. (np.square(region_sim_vals) + epsilon)
  558. ) * 100
  559. errors[algorithm][snr].append(relative_squared_error)
  561. all_valid_coords = np.concatenate(all_valid_coords, axis=0)
  562. all_theta_values = np.concatenate(all_theta_values, axis=0)
  563. all_sim_values = np.concatenate(all_sim_values, axis=0)
  565. for algorithm in algorithms:
  566. for snr in snr_levels:
  567. errors[algorithm][snr] = np.concatenate(errors[algorithm][snr], axis=0)
  569. # ----------------------------
  570. # 4) BIN THE THETA VALUES
  571. # ----------------------------
  573. bins = np.arange(0, 100, 10)
  575. bin_indices = np.digitize(all_theta_values, bins)
  576. num_bins = len(bins) - 1
  578. # Calculate the mean bin angle
  579. bin_means_theta = []
  580. for b in range(1, len(bins)):
  581. mask_b = (bin_indices == b)
  582. if np.any(mask_b):
  583. bin_means_theta.append(np.mean(all_theta_values[mask_b]))
  584. bin_means_theta = np.array(bin_means_theta)
  586. # ----------------------------
  587. # 5) CALCULATE BIN MEANS OF MSPE FOR EACH ALGORITHM AND SNR
  588. # ----------------------------
  590. bin_means = {
  591. algorithm: {
  592. snr: [] for snr in snr_levels
  593. } for algorithm in algorithms
  594. }
  596. for b in range(1, len(bins)):
  597. mask_b = (bin_indices == b)
  598. if np.any(mask_b):
  599. for algorithm in algorithms:
  600. for snr in snr_levels:
  601. mspe_values_bin = errors[algorithm][snr][mask_b]
  602. bin_mean_error = np.mean(mspe_values_bin)
  603. bin_means[algorithm][snr].append(bin_mean_error)
  604. else:
  606. for algorithm in algorithms:
  607. for snr in snr_levels:
  608. bin_means[algorithm][snr].append(np.nan)
  610. # Convert bin means to numpy arrays
  611. for algorithm in algorithms:
  612. for snr in snr_levels:
  613. bin_means[algorithm][snr] = np.array(bin_means[algorithm][snr])
  615. # ----------------------------
  616. # 6) SAVE A SINGLE 3D MASK WITH BIN LABELS
  617. # ----------------------------
  619. all_bins_3d = np.zeros(theta_shape, dtype=np.uint8)
  621. for b in range(1, len(bins)):
  622. # Voxels in bin b
  623. mask_b = (bin_indices == b)
  625. # Coordinates of voxels in bin b
  626. bin_coords = all_valid_coords[mask_b]
  627. all_bins_3d[bin_coords[:, 0], bin_coords[:, 1], bin_coords[:, 2]] = b
  629. # ----------------------------
  630. # 7) PLOT THE RESULTS
  631. # ----------------------------
  633. markers = {
  634. 'SNR 300': ('#e78ac3', 'o'),
  635. 'SNR 200': ('#fc8d62', 'o'),
  636. 'SNR 100': ('#8da0cb', 'o'),
  637. 'SNR 50': ('#66c2a5', 'o')
  638. }
  640. fig, axs = plt.subplots(2, 2, figsize=(14, 10))
  641. algorithm_subplot_indices = {
  642. 'x-separation': (0, 0),
  643. 'R2*-QSM': (0, 1),
  644. 'APART-QSM': (1, 0),
  645. 'DECOMPOSE-QSM': (1, 1)
  646. }
  648. for algorithm in algorithms:
  649. ax = axs[algorithm_subplot_indices[algorithm]]
  651. # Dictionary to store percentage changes for each SNR
  652. percentage_changes = {}
  654. for snr in snr_levels:
  655. color, marker = markers[snr]
  657. # x-values = bin_means_theta
  658. x_vals = bin_means_theta
  659. # y-values = bin_means[algorithm][snr]
  660. y_vals = bin_means[algorithm][snr]
  662. # Scatter plot
  663. ax.scatter(
  664. x_vals,
  665. y_vals,
  666. color=color,
  667. marker=marker,
  668. label=snr
  669. )
  671. # 1) Sort the points by x so we can connect them in ascending order
  672. sort_idx = np.argsort(x_vals)
  673. x_sorted = x_vals[sort_idx]
  674. y_sorted = y_vals[sort_idx]
  676. # Calculate percentage change (max vs. min) for this SNR
  677. valid_mspe = y_vals[~np.isnan(y_vals)]
  678. if len(valid_mspe) > 0:
  679. max_mspe = y_vals[0]
  680. min_mspe = y_vals[-1]
  681. if max_mspe != 0:
  682. percentage_change = (max_mspe - min_mspe) / max_mspe * 100
  683. else:
  684. percentage_change = 0
  685. else:
  686. percentage_change = 0
  688. percentage_changes[snr] = percentage_change
  690. # Build the "MEV=" Variation text with Matplotlib colors
  691. x_position = 0.47
  692. y_position = 0.05
  694. ax.text(
  695. x_position, y_position,
  696. "MEV=",
  697. transform=ax.transAxes, fontsize=12, color='black', ha='left', va='top'
  698. )
  699. x_position += 0.1
  701. text_colors = ['#e78ac3', '#fc8d62', '#8da0cb', '#66c2a5']
  702. for idx, snr in enumerate(snr_levels):
  703. color_ = text_colors[idx]
  704. ax.text(
  705. x_position, y_position,
  706. f"{percentage_changes[snr]:.1f}%",
  707. transform=ax.transAxes, fontsize=12, color=color_, ha='left', va='top'
  708. )
  709. x_position += 0.1
  711. # Set titles and labels
  712. if algorithm == 'x-separation':
  713. ax.set_title(r"$\chi$-separation")
  714. else:
  715. ax.set_title(algorithm)
  717. ax.set_xlabel('Angle of Orientation (degrees)', fontsize=16)
  718. ax.set_ylabel('MSPE of $\\chi^-$ (%)', fontsize=16)
  719. ax.grid(True)
  720. ax.legend()
  721. ax.set_ylim(0, 130)
  722. ax.tick_params(axis='both', which='major', labelsize=16)
  724. plt.tight_layout()
  725. plt.show()
  728. # %% [markdown]
  729. # # Comparison between simulated vs in-vivo susceptibility maps
  731. # %%
  732. # Define regions of interest
  733. regions_of_interest = [1, 2, 3, 7, 8, 9]
  735. # Function to extract mean values per region
  736. def extract_mean_per_region(data_map, segmentation_map, regions):
  737. means = []
  738. for region in regions:
  739. region_mask = segmentation_map == region
  740. if np.any(region_mask):
  741. mean_value = np.mean(data_map[region_mask])
  742. means.append(mean_value)
  743. else:
  744. means.append(np.nan)
  745. return means
  747. # ------------------------------------------------------------------------------
  748. # Load segmentation map (for simulated data only)
  749. # ------------------------------------------------------------------------------
  750. segmentation_simulated = nib.load(
  751. base_dir / 'osfstorage' / 'Masks' / 'SegmentedModel.nii.gz'
  752. ).get_fdata()
  754. # ------------------------------------------------------------------------------
  755. # Hard-coded in-vivo mean values (from your table) for each algorithm and region
  756. # ------------------------------------------------------------------------------
  757. # Chi-positive in-vivo data
  758. chi_positive_in_vivo_data = {
  759. 'chi-separation': [0.0658392,0.113222,0.0550474,0.0475523,0.0270498,0.0299241],
  760. 'R2*-QSM': [0.0357698,0.0610981,0.0347473,0.0265003,0.0159944,0.0140221],
  761. 'APART-QSM': [0.0482199,0.0813391,0.0435077,0.0346825,0.0205208,0.02208],
  762. 'DECOMPOSE-QSM': [0.0337807,0.0538997,0.0325998,0.0301356,0.0157292,0.0112171]
  763. }
  765. # Chi-negative in-vivo data (already positive)
  766. chi_negative_in_vivo_data = {
  767. 'chi-separation': [-0.0175603,-0.0332972,-0.0260955,-0.0270215,-0.0303985,-0.0277441],
  768. 'R2*-QSM': [-0.00613275,-0.00185815,-0.00892519,-0.0140371,-0.0187294,-0.0130369],
  769. 'APART-QSM': [-0.0105914,-0.0195933,-0.0182477,-0.0190805,-0.0222699,-2.15E-02],
  770. 'DECOMPOSE-QSM': [-0.00877939,-0.00756195,-0.00899253,-0.00902613,-0.0163517,-0.00915724]
  771. }
  773. # ------------------------------------------------------------------------------
  774. # Initialize dictionaries to store mean values
  775. # ------------------------------------------------------------------------------
  776. chi_positive_in_vivo_means = {}
  777. chi_positive_simulated_means = {}
  778. chi_negative_in_vivo_means = {}
  779. chi_negative_simulated_means = {}
  781. # ------------------------------------------------------------------------------
  782. # File paths for simulated data
  783. # ------------------------------------------------------------------------------
  784. chi_positive_simulated_files = {
  785. 'chi-separation': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'X-separation' / 'With_anisotropy' / 'ChiPosMap.nii',
  786. 'R2*-QSM': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'R2star-QSM' / 'With_anisotropy' / 'ChiPosMap.nii',
  787. 'APART-QSM': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'APART-QSM' / 'With_anisotropy' / 'X_para.nii',
  788. 'DECOMPOSE-QSM': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'DECOMPOSE-QSM'/ 'With_anisotropy' / 'Results_PCS.nii'
  789. }
  791. chi_negative_simulated_files = {
  792. 'chi-separation': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'X-separation' / 'With_anisotropy' / 'ChiNegMap.nii',
  793. 'R2*-QSM': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'R2star-QSM' / 'With_anisotropy' / 'ChiNegMap.nii',
  794. 'APART-QSM': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'APART-QSM' / 'With_anisotropy' / 'X_dia_abs.nii',
  795. 'DECOMPOSE-QSM': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'DECOMPOSE-QSM'/ 'With_anisotropy' / 'Results_DCS_abs.nii'
  796. }
  798. algorithms = list(chi_positive_simulated_files.keys())
  800. # ------------------------------------------------------------------------------
  801. # Process chi-positive data
  802. # ------------------------------------------------------------------------------
  803. for algo in algorithms:
  804. # In-vivo: use the hard-coded values
  805. in_vivo_means = chi_positive_in_vivo_data[algo]
  806. # Simulated: load from file and extract
  807. simulated_map = nib.load(chi_positive_simulated_files[algo]).get_fdata()
  808. simulated_means = extract_mean_per_region(simulated_map, segmentation_simulated, regions_of_interest)
  810. chi_positive_in_vivo_means[algo] = in_vivo_means
  811. chi_positive_simulated_means[algo] = simulated_means
  813. # ------------------------------------------------------------------------------
  814. # Process chi-negative data
  815. # ------------------------------------------------------------------------------
  816. for algo in algorithms:
  817. # In-vivo: use the hard-coded values
  818. in_vivo_means = chi_negative_in_vivo_data[algo]
  819. # Simulated: load from file, multiply by -1
  820. simulated_map = nib.load(chi_negative_simulated_files[algo]).get_fdata()
  821. simulated_map *= -1
  822. simulated_means = extract_mean_per_region(simulated_map, segmentation_simulated, regions_of_interest)
  824. chi_negative_in_vivo_means[algo] = in_vivo_means
  825. chi_negative_simulated_means[algo] = simulated_means
  827. # ------------------------------------------------------------------------------
  828. # Define display names
  829. # ------------------------------------------------------------------------------
  830. display_names = {
  831. 'chi-separation': r'$\chi$-separation',
  832. 'R2*-QSM': 'R2*-QSM',
  833. 'APART-QSM': 'APART-QSM',
  834. 'DECOMPOSE-QSM': 'DECOMPOSE-QSM'
  835. }
  837. # ------------------------------------------------------------------------------
  838. # Plotting function
  839. # ------------------------------------------------------------------------------
  840. def plot_data(in_vivo_means_dict, simulated_means_dict, chi_type, ax):
  841. marker_styles = ['x', 'o', '*', '+'] # for the 4 algorithms in the order they appear in 'algorithms'
  843. for idx, algo in enumerate(algorithms):
  844. x = np.array(in_vivo_means_dict[algo])
  845. y = np.array(simulated_means_dict[algo])
  847. mask = ~np.isnan(x) & ~np.isnan(y)
  848. x = x[mask]
  849. y = y[mask]
  851. if len(x) > 1:
  852. slope, intercept, r_value, p_value, std_err = linregress(x, y)
  853. r_label = f"(r={r_value:.2f})"
  854. sorted_indices = np.argsort(x)
  855. x_sorted = x[sorted_indices]
  856. y_sorted = intercept + slope * x_sorted
  857. ax.plot(x_sorted, y_sorted, color='black', linestyle='--', label=None)
  858. else:
  859. r_label = "(R²=NaN)"
  861. label = f"{display_names[algo]} {r_label}"
  862. marker = marker_styles[idx]
  864. if marker == 'x':
  865. ax.scatter(x, y, marker=marker, s=80, color='black',
  866. linewidth=1.5, label=label)
  867. elif marker == 'o':
  868. ax.scatter(x, y, marker=marker, s=80,
  869. edgecolors='black', facecolors='none',
  870. linewidth=1.5, label=label)
  871. else:
  872. ax.scatter(x, y, marker=marker, s=80,
  873. color='black', linewidth=1.5, label=label)
  875. x_limits = ax.get_xlim()
  876. y_limits = ax.get_ylim()
  877. min_val = min(x_limits[0], y_limits[0])
  878. max_val = max(x_limits[1], y_limits[1])
  879. ax.set_xlim([min_val, max_val])
  880. ax.set_ylim([min_val, max_val])
  882. ax.plot([min_val, max_val], [min_val, max_val], '-', color='grey', alpha=0.3)
  883. ax.set_xlabel(f'Measured {chi_type} in-vivo (ppm)', fontsize=18)
  884. ax.set_ylabel(f'Simulated {chi_type} (ppm)', fontsize=18)
  885. ax.legend(loc='upper left', fontsize=14)
  886. ax.grid(True)
  887. ax.tick_params(axis='both', which='major', labelsize=14)
  889. # ------------------------------------------------------------------------------
  890. # Create the figure and subplots
  891. # ------------------------------------------------------------------------------
  892. fig, axes = plt.subplots(1, 2, figsize=(20, 6))
  894. plot_data(chi_positive_in_vivo_means, chi_positive_simulated_means, r'$\chi^+$', axes[0])
  895. plot_data(chi_negative_in_vivo_means, chi_negative_simulated_means, r'$\chi^-$', axes[1])
  897. for ax in axes:
  898. ax.tick_params(axis='both', which='major', labelsize=16)
  900. plt.tight_layout()
  901. plt.show()
  904. # %% [markdown]
  905. # # Comparison between simulated vs in-vivo field maps
  907. # %%
  908. # Construct file paths relative to the identified base directory
  909. simulated_map_path = base_dir / "osfstorage/LocalField/Simulated_Field.nii.gz"
  910. segmentation_map_path = base_dir / "osfstorage/Masks/white_matter_mask.nii.gz"
  912. # Load the NIfTI files for simulated local field map and the segmentation map
  913. simulated_map_nii = nib.load(simulated_map_path)
  914. segmentation_map_nii = nib.load(segmentation_map_path)
  916. # Extract the data arrays
  917. simulated_map = simulated_map_nii.get_fdata()
  918. segmentation_map = segmentation_map_nii.get_fdata()
  920. # Extract simulated values for scatter plot (regions 1 to 10 in the segmentation map)
  921. simulated_values = []
  922. for region in range(1, 11):
  923. region_mask = segmentation_map == region
  924. simulated_region_values = simulated_map[region_mask]
  925. simulated_values.append(np.mean(simulated_region_values))
  927. # Hard-coded in-vivo values (for the same 10 regions)
  928. in_vivo_values = np.array([
  929. 1.74348,
  930. -1.72861,
  931. -1.59869,
  932. 1.46938,
  933. -2.61981,
  934. -0.557609,
  935. -0.714494,
  936. -1.31104,
  937. -2.75768,
  938. 0.288124
  939. ])
  941. # Scatter plot with linear regression
  942. slope, intercept, r_value, p_value, std_err = stats.linregress(in_vivo_values, simulated_values)
  943. regression_line = slope * in_vivo_values + intercept
  945. # Create the scatter plot with a smaller figure size
  946. fig, ax = plt.subplots(figsize=(6, 6))
  948. # Sort the x-values (in-vivo local field values) and calculate corresponding y-values
  949. sorted_indices = np.argsort(in_vivo_values)
  950. sorted_in_vivo_values = in_vivo_values[sorted_indices]
  951. sorted_regression_line = slope * sorted_in_vivo_values + intercept
  953. # Scatter plot
  954. ax.scatter(in_vivo_values, simulated_values, color='black', zorder=2)
  956. # Plot the regression line
  957. ax.plot(sorted_in_vivo_values, sorted_regression_line, linestyle='--',
  958. color='black', linewidth=1.2, zorder=1)
  960. # Add labels, grid, and R² text
  961. ax.set_xlabel('In-vivo local field (Hz)', fontsize=10)
  962. ax.set_ylabel('Simulated local field (Hz)', fontsize=10)
  963. ax.grid(True, alpha=0.3)
  964. ax.text(0.05, 0.95,
  965. f'Correlation = {r_value:.2f}\nSlope = {slope:.1f}',
  966. ha='left', va='top', transform=ax.transAxes, fontsize=14)
  968. plt.tight_layout()
  969. plt.show()
  971. # %% [markdown]
  972. # # Comparison between simulated vs in-vivo T2 maps
  974. # %%
  975. # File paths for simulated T2 map and segmentation map
  976. simulated_t2_path = base_dir / "osfstorage/T2/T2_simulated_resampled.nii.gz"
  977. simulated_segmentation_path = base_dir / "osfstorage/Masks/SegmentedModel_resampled.nii.gz"
  979. # Load NIfTI files for simulated data
  980. simulated_t2_map = nib.load(simulated_t2_path).get_fdata()
  981. simulated_segmentation = nib.load(simulated_segmentation_path).get_fdata()
  983. # Replace NaN and Inf values with 0 in the simulated T2 map
  984. simulated_t2_map = np.nan_to_num(simulated_t2_map, nan=0, posinf=0, neginf=0)
  986. # Hard-coded in-vivo T2 values for regions 1-3 and 7-9 (6 regions)
  987. T2_in_vivo = np.array([67.2624, 44.9122, 52.7201, 59.789, 54.5211, 104.057])
  989. # Extract simulated T2 values for regions 1-3 and 7-9
  990. T2_simulated = []
  991. for region in list(range(1, 4)) + list(range(7, 10)):
  992. simulated_mask = simulated_segmentation == region
  993. if np.any(simulated_mask):
  994. simulated_mean = np.mean(simulated_t2_map[simulated_mask])
  995. T2_simulated.append(simulated_mean)
  996. T2_simulated = np.array(T2_simulated)
  998. # Reshape data for linear regression
  999. T2_in_vivo_reshaped = T2_in_vivo.reshape(-1, 1)
  1000. T2_simulated_reshaped = T2_simulated.reshape(-1, 1)
  1002. # Perform linear regression
  1003. reg_model = LinearRegression()
  1004. reg_model.fit(T2_in_vivo_reshaped, T2_simulated_reshaped)
  1005. slope = reg_model.coef_[0][0]
  1006. intercept = reg_model.intercept_[0]
  1008. # Generate regression line values
  1009. x_fit = np.linspace(T2_in_vivo.min(), T2_in_vivo.max(), 100)
  1010. y_fit = slope * x_fit + intercept
  1012. # Plotting
  1013. plt.figure(figsize=(8, 6))
  1014. plt.scatter(T2_in_vivo, T2_simulated, color='black', label='Data Points')
  1015. plt.plot(x_fit, y_fit, '--', color='black', label='Regression Line')
  1016. plt.xlabel(r'$T_2^{\mathrm{in-vivo}}$ (ms)', fontsize=16)
  1017. plt.ylabel(r'$T_2^{\mathrm{simulated}}$ (ms)', fontsize=16)
  1019. correlation = np.corrcoef(T2_in_vivo, T2_simulated)[0, 1]
  1020. plt.text(45, 120, f'Correlation: {correlation:.2f}\nSlope: {slope:.2f}', fontsize=18)
  1022. plt.grid(True)
  1023. plt.xticks(fontsize=14)
  1024. plt.yticks(fontsize=14)
  1025. plt.tight_layout()
  1026. plt.show()
  1029. # %% [markdown]
  1030. # # Supplimantary Material
  1032. # %%
  1033. num_algorithms = 4
  1034. algorithm_names = [
  1035. '$\\chi$-separation',
  1036. 'R2*-QSM',
  1037. 'APART-QSM',
  1038. 'DECOMPOSE-QSM'
  1039. ]
  1041. segmentation2_path = base_dir / 'osfstorage' / 'Masks' / 'white_matter_mask.nii.gz'
  1042. segmentation2 = nib.load(segmentation2_path).get_fdata()
  1044. region_labels = list(range(1, 10))
  1046. # Load simulated maps
  1047. simulated_with_path = base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'Chi_positive.nii.gz'
  1048. simulated_without_path = base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'Chi_positive.nii.gz'
  1050. simulated_with = nib.load(simulated_with_path).get_fdata()
  1051. simulated_without = nib.load(simulated_without_path).get_fdata()
  1053. # Define measured maps for each algorithm
  1054. measured_maps = {
  1055. 0: { # Algorithm 1: χ-separation
  1056. 'x_minus_with_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'X-separation' / 'With_anisotropy' / 'ChiPosMap.nii',
  1057. 'x_minus_without_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'X-separation' / 'Without_anisotropy' / 'ChiPosMap.nii',
  1058. },
  1059. 1: { # Algorithm 2: R2*-QSM
  1060. 'x_minus_with_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'R2star-QSM' / 'With_anisotropy' / 'ChiPosMap.nii',
  1061. 'x_minus_without_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'R2star-QSM' / 'Without_anisotropy' / 'ChiPosMap.nii',
  1062. },
  1063. 2: { # Algorithm 3: APART-QSM
  1064. 'x_minus_with_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'APART-QSM' / 'With_anisotropy' / 'X_para.nii',
  1065. 'x_minus_without_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'APART-QSM' / 'Without_anisotropy' / 'X_para.nii',
  1066. },
  1067. 3: { # Algorithm 4: DECOMPOSE-QSM
  1068. 'x_minus_with_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'DECOMPOSE-QSM' / 'With_anisotropy' / 'Results_PCS.nii',
  1069. 'x_minus_without_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'DECOMPOSE-QSM' / 'Without_anisotropy' / 'results_PCS.nii',
  1070. },
  1071. }
  1073. fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12, 12), dpi=300)
  1075. for idx in range(num_algorithms):
  1076. ax = axes[idx // 2, idx % 2]
  1077. algorithm_name = algorithm_names[idx]
  1079. # Load the measured maps for the current algorithm
  1080. measured_with = nib.load(measured_maps[idx]['x_minus_with_anisotropy']).get_fdata()
  1081. measured_without = nib.load(measured_maps[idx]['x_minus_without_anisotropy']).get_fdata()
  1082. errors_with_regions = []
  1083. errors_without_regions = []
  1085. # --------------------------------------------------
  1086. # 1) Compute errors
  1087. # --------------------------------------------------
  1088. for region in region_labels:
  1089. region_mask = (segmentation2 == region)
  1090. if region == 8:
  1091. continue
  1093. if np.any(region_mask):
  1094. # "With Anisotropy"
  1095. sim_val_with = np.mean(simulated_with[region_mask])
  1097. if sim_val_with < 0:
  1098. sim_val_with *= -1
  1099. meas_val_with = np.mean(measured_with[region_mask])
  1100. if sim_val_with != 0:
  1101. error_with = ((meas_val_with - sim_val_with) / sim_val_with) ** 2 * 100
  1102. else:
  1103. error_with = np.nan
  1105. # "Without Anisotropy"
  1106. sim_val_without = np.mean(simulated_without[region_mask])
  1107. meas_val_without = np.mean(measured_without[region_mask])
  1108. if sim_val_without != 0:
  1109. error_without = ((meas_val_without - sim_val_without) / sim_val_without) ** 2 * 100
  1110. else:
  1111. error_without = np.nan
  1113. errors_with_regions.append(error_with)
  1114. errors_without_regions.append(error_without)
  1116. # Convert to numpy arrays and remove NaNs
  1117. errors_with_regions = np.array(errors_with_regions)
  1118. errors_without_regions = np.array(errors_without_regions)
  1119. errors_with_regions = errors_with_regions[~np.isnan(errors_with_regions)]
  1120. errors_without_regions = errors_without_regions[~np.isnan(errors_without_regions)]
  1122. # Prepare for KDE
  1123. x_vals = np.linspace(0, 30000, 200)
  1124. max_height = 0.4
  1126. # Plot KDE for "With Anisotropy"
  1127. if errors_with_regions.size > 1:
  1128. kde_with = gaussian_kde(errors_with_regions)
  1129. density_with = kde_with(x_vals)
  1130. scaling_with = max_height / np.max(density_with) if np.max(density_with) > 0 else 1
  1131. density_with_scaled = density_with * scaling_with
  1132. ax.fill_between(x_vals, 0, density_with_scaled, alpha=0.4, color='tab:blue')
  1133. ax.plot(x_vals, density_with_scaled, alpha=0.7, color='tab:blue')
  1134. mean_mspe_with = np.mean(errors_with_regions)
  1135. elif errors_with_regions.size == 1:
  1136. mean_mspe_with = errors_with_regions[0]
  1137. ax.plot([mean_mspe_with], [max_height/2], marker='o', color='tab:blue')
  1138. else:
  1139. mean_mspe_with = None
  1141. # Plot KDE for "Without Anisotropy"
  1142. if errors_without_regions.size > 1:
  1143. kde_without = gaussian_kde(errors_without_regions)
  1144. density_without = kde_without(x_vals)
  1145. scaling_without = max_height / np.max(density_without) if np.max(density_without) > 0 else 1
  1146. density_without_scaled = density_without * scaling_without
  1147. ax.fill_between(x_vals, 0, density_without_scaled, alpha=0.4, color='tab:orange')
  1148. ax.plot(x_vals, density_without_scaled, alpha=0.7, color='tab:orange')
  1149. mean_mspe_without = np.mean(errors_without_regions)
  1150. elif errors_without_regions.size == 1:
  1151. mean_mspe_without = errors_without_regions[0]
  1152. ax.plot([mean_mspe_without], [max_height/2], marker='o', color='tab:orange')
  1153. else:
  1154. mean_mspe_without = None
  1156. # Show mean MSPE for with/without anisotropy
  1157. if mean_mspe_with is not None:
  1158. ax.text(0.95, 0.8, f"Mean MSPE: {mean_mspe_with:.1f}%", transform=ax.transAxes,
  1159. color='tab:blue', ha='right', va='center', fontsize=12)
  1160. if mean_mspe_without is not None:
  1161. ax.text(0.95, 0.7, f"Mean MSPE: {mean_mspe_without:.1f}%", transform=ax.transAxes,
  1162. color='tab:orange', ha='right', va='center', fontsize=12)
  1164. legend_elements = [
  1165. Line2D([0], [0], color='tab:blue', lw=2, label='With Anisotropy'),
  1166. Line2D([0], [0], color='tab:orange', lw=2, label='Without Anisotropy')
  1167. ]
  1168. ax.legend(handles=legend_elements, loc='upper right')
  1170. ax.set_title(algorithm_name, fontsize=14)
  1171. ax.set_xlabel("MSPE of $\\chi^+$ (%)", fontsize=12)
  1172. ax.set_ylabel("Density (a.u.)", fontsize=12)
  1173. ax.set_ylim(0, max_height * 1.2)
  1174. ax.grid(alpha=0.3)
  1176. # ---------------------------
  1177. # Final figure adjustments
  1178. # ---------------------------
  1179. plt.tight_layout(rect=[0, 0, 1, 0.96])
  1180. plt.show()
  1183. # %%
  1184. num_algorithms = 4
  1185. algorithm_names = [
  1186. '$\\chi$-separation',
  1187. 'R2*-QSM',
  1188. 'APART-QSM',
  1189. 'DECOMPOSE-QSM'
  1190. ]
  1192. segmentation2_path = base_dir / 'osfstorage' / 'Masks' / 'SegmentedModel.nii.gz'
  1193. segmentation2 = nib.load(segmentation2_path).get_fdata()
  1195. region_labels = list(range(1, 10))
  1197. # Load simulated maps
  1198. simulated_with_path = base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'Chi_positive.nii.gz'
  1199. simulated_without_path = base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'Chi_positive.nii.gz'
  1201. simulated_with = nib.load(simulated_with_path).get_fdata()
  1202. simulated_without = nib.load(simulated_without_path).get_fdata()
  1204. # Define measured maps for each algorithm
  1205. measured_maps = {
  1206. 0: { # Algorithm 1: χ-separation
  1207. 'x_minus_with_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'X-separation' / 'With_anisotropy' / 'ChiPosMap.nii',
  1208. 'x_minus_without_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'X-separation' / 'Without_anisotropy' / 'ChiPosMap.nii',
  1209. },
  1210. 1: { # Algorithm 2: R2*-QSM
  1211. 'x_minus_with_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'R2star-QSM' / 'With_anisotropy' / 'ChiPosMap.nii',
  1212. 'x_minus_without_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'R2star-QSM' / 'Without_anisotropy' / 'ChiPosMap.nii',
  1213. },
  1214. 2: { # Algorithm 3: APART-QSM
  1215. 'x_minus_with_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'APART-QSM' / 'With_anisotropy' / 'X_para.nii',
  1216. 'x_minus_without_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'APART-QSM' / 'Without_anisotropy' / 'X_para.nii',
  1217. },
  1218. 3: { # Algorithm 4: DECOMPOSE-QSM
  1219. 'x_minus_with_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'DECOMPOSE-QSM' / 'With_anisotropy' / 'Results_PCS.nii',
  1220. 'x_minus_without_anisotropy': base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'DECOMPOSE-QSM' / 'Without_anisotropy' / 'results_PCS.nii',
  1221. },
  1222. }
  1224. fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12, 12), dpi=300)
  1226. for idx in range(num_algorithms):
  1227. ax = axes[idx // 2, idx % 2]
  1228. algorithm_name = algorithm_names[idx]
  1230. # Load the measured maps for the current algorithm
  1231. measured_with = nib.load(measured_maps[idx]['x_minus_with_anisotropy']).get_fdata()
  1232. measured_without = nib.load(measured_maps[idx]['x_minus_without_anisotropy']).get_fdata()
  1233. errors_with_regions = []
  1234. errors_without_regions = []
  1236. # --------------------------------------------------
  1237. # 1) Compute errors
  1238. # --------------------------------------------------
  1239. for region in region_labels:
  1240. region_mask = (segmentation2 == region)
  1242. if region == 8:
  1243. continue
  1244. if np.any(region_mask):
  1245. # "With Anisotropy"
  1246. sim_val_with = np.mean(simulated_with[region_mask])
  1248. if sim_val_with < 0:
  1249. sim_val_with *= -1
  1250. meas_val_with = np.mean(measured_with[region_mask])
  1251. if sim_val_with != 0:
  1252. error_with = ((meas_val_with - sim_val_with) / sim_val_with) ** 2 * 100
  1253. else:
  1254. error_with = np.nan
  1256. # "Without Anisotropy"
  1257. sim_val_without = np.mean(simulated_without[region_mask])
  1258. meas_val_without = np.mean(measured_without[region_mask])
  1259. if sim_val_without != 0:
  1260. error_without = ((meas_val_without - sim_val_without) / sim_val_without) ** 2 * 100
  1261. else:
  1262. error_without = np.nan
  1264. errors_with_regions.append(error_with)
  1265. errors_without_regions.append(error_without)
  1267. # Convert to numpy arrays and remove NaNs
  1268. errors_with_regions = np.array(errors_with_regions)
  1269. errors_without_regions = np.array(errors_without_regions)
  1270. errors_with_regions = errors_with_regions[~np.isnan(errors_with_regions)]
  1271. errors_without_regions = errors_without_regions[~np.isnan(errors_without_regions)]
  1273. # Prepare for KDE
  1274. x_vals = np.linspace(0, 70, 70)
  1275. max_height = 0.4
  1277. # Plot KDE for "With Anisotropy"
  1278. if errors_with_regions.size > 1:
  1279. kde_with = gaussian_kde(errors_with_regions)
  1280. density_with = kde_with(x_vals)
  1281. scaling_with = max_height / np.max(density_with) if np.max(density_with) > 0 else 1
  1282. density_with_scaled = density_with * scaling_with
  1283. ax.fill_between(x_vals, 0, density_with_scaled, alpha=0.4, color='tab:blue')
  1284. ax.plot(x_vals, density_with_scaled, alpha=0.7, color='tab:blue')
  1285. mean_mspe_with = np.mean(errors_with_regions)
  1286. elif errors_with_regions.size == 1:
  1287. mean_mspe_with = errors_with_regions[0]
  1288. ax.plot([mean_mspe_with], [max_height/2], marker='o', color='tab:blue')
  1289. else:
  1290. mean_mspe_with = None
  1292. # Plot KDE for "Without Anisotropy"
  1293. if errors_without_regions.size > 1:
  1294. kde_without = gaussian_kde(errors_without_regions)
  1295. density_without = kde_without(x_vals)
  1296. scaling_without = max_height / np.max(density_without) if np.max(density_without) > 0 else 1
  1297. density_without_scaled = density_without * scaling_without
  1298. ax.fill_between(x_vals, 0, density_without_scaled, alpha=0.4, color='tab:orange')
  1299. ax.plot(x_vals, density_without_scaled, alpha=0.7, color='tab:orange')
  1300. mean_mspe_without = np.mean(errors_without_regions)
  1301. elif errors_without_regions.size == 1:
  1302. mean_mspe_without = errors_without_regions[0]
  1303. ax.plot([mean_mspe_without], [max_height/2], marker='o', color='tab:orange')
  1304. else:
  1305. mean_mspe_without = None
  1307. # Show mean MSPE for with/without anisotropy
  1308. if mean_mspe_with is not None:
  1309. ax.text(0.95, 0.8, f"Mean MSPE: {mean_mspe_with:.1f}%", transform=ax.transAxes,
  1310. color='tab:blue', ha='right', va='center', fontsize=12)
  1311. if mean_mspe_without is not None:
  1312. ax.text(0.95, 0.7, f"Mean MSPE: {mean_mspe_without:.1f}%", transform=ax.transAxes,
  1313. color='tab:orange', ha='right', va='center', fontsize=12)
  1315. legend_elements = [
  1316. Line2D([0], [0], color='tab:blue', lw=2, label='With Anisotropy'),
  1317. Line2D([0], [0], color='tab:orange', lw=2, label='Without Anisotropy')
  1318. ]
  1319. ax.legend(handles=legend_elements, loc='upper right')
  1321. ax.set_title(algorithm_name, fontsize=14)
  1322. ax.set_xlabel("MSPE of $\\chi^+$ (%)", fontsize=12)
  1323. ax.set_ylabel("Density (a.u.)", fontsize=12)
  1324. ax.set_ylim(0, max_height * 1.2)
  1325. ax.grid(alpha=0.3)
  1327. # ---------------------------
  1328. # Final figure adjustments
  1329. # ---------------------------
  1330. plt.tight_layout(rect=[0, 0, 1, 0.96])
  1331. plt.show()
  1334. # %%
  1335. methods = {
  1336. "χ-separation": (
  1337. base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'X-separation' / 'With_anisotropy' / 'ChiNegMap.nii',
  1338. base_dir / 'osfstorage' / '7T' / 'X-separation' / 'ChiNegMap.nii'
  1339. ),
  1340. "R2*-QSM": (
  1341. base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'R2star-QSM' / 'With_anisotropy' / 'ChiNegMap.nii',
  1342. base_dir / 'osfstorage' / '7T' / 'R2star-QSM' / 'ChiNegMap.nii'
  1343. ),
  1344. "APART-QSM": (
  1345. base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'APART-QSM' / 'With_anisotropy' / 'X_dia_abs.nii',
  1346. base_dir / 'osfstorage' / '7T' / 'APART-QSM' / 'X_dia_abs.nii'
  1347. ),
  1348. "DECOMPOSE-QSM": (
  1349. base_dir / 'osfstorage' / 'Susceptibility_Separation_Results' / 'DECOMPOSE-QSM' / 'With_anisotropy' / 'Results_DCS_abs.nii',
  1350. base_dir / 'osfstorage' / '7T' / 'DECOMPOSE-QSM' / 'Results_DCS_abs.nii'
  1351. ),
  1352. }
  1354. #-----------------------------------------------------------------
  1355. # Load segmentation file
  1356. #-----------------------------------------------------------------
  1357. segmentation_file = base_dir / 'osfstorage' / 'Masks' / 'white_matter_mask.nii.gz'
  1358. seg_nii = nib.load(segmentation_file)
  1359. seg_data = seg_nii.get_fdata()
  1360. labels = range(1, 12)
  1362. method_points = {}
  1364. all_mean_values = []
  1365. all_difference_values = []
  1367. for method_name, (nii_3T_path, nii_7T_path) in methods.items():
  1368. # --------------------
  1369. # Load the 3T and 7T
  1370. # --------------------
  1371. nii_3T = nib.load(nii_3T_path)
  1372. nii_7T = nib.load(nii_7T_path)
  1374. data_3T = nii_3T.get_fdata()*-1
  1375. data_7T = nii_7T.get_fdata()*-1
  1377. mean_values = np.zeros(len(labels))
  1378. diff_values = np.zeros(len(labels))
  1380. for i, label_val in enumerate(labels):
  1381. # Mask for current region
  1382. mask = (seg_data == label_val)
  1384. # Extract the region's data for 3T and 7T
  1385. region_3T = data_3T[mask]
  1386. region_7T = data_7T[mask]
  1388. # Compute the mean within that region
  1389. mean_3T = np.mean(region_3T) if region_3T.size > 0 else np.nan
  1390. mean_7T = np.mean(region_7T) if region_7T.size > 0 else np.nan
  1392. # Bland-Altman points
  1393. mean_values[i] = (mean_3T + mean_7T) / 2
  1394. diff_values[i] = (mean_3T - mean_7T)
  1396. # Store in dictionary
  1397. method_points[method_name] = (mean_values, diff_values)
  1399. # Accumulate for global axis limits
  1400. all_mean_values.extend(mean_values)
  1401. all_difference_values.extend(diff_values)
  1403. # Convert to NumPy for min/max
  1404. all_mean_values = np.array(all_mean_values)
  1405. all_difference_values = np.array(all_difference_values)
  1407. # Determine global x-limits and y-limits
  1408. x_min, x_max = np.min(all_mean_values), np.max(all_mean_values)
  1409. x_margin = 0.1 * (x_max - x_min) # 10% margin
  1410. x_lim = (x_min - x_margin, x_max + x_margin)
  1412. # Fixed y-limits as in your code
  1413. y_lim = (-0.05, 0.05)
  1415. # Create subplots
  1416. fig, axes = plt.subplots(2, 2, figsize=(10, 8))
  1417. axes = axes.flatten()
  1419. for ax, (method_name, (mean_values, difference_values)) in zip(axes, method_points.items()):
  1420. # Mean & LoA
  1421. mean_diff = np.mean(difference_values)
  1422. std_diff = np.std(difference_values, ddof=1)
  1423. loa_upper = mean_diff + 1.96 * std_diff
  1424. loa_lower = mean_diff - 1.96 * std_diff
  1426. ax.scatter(mean_values, difference_values, marker="x", s=50)
  1428. # Plot lines
  1429. ax.axhline(mean_diff, color="red", linestyle="--")
  1430. ax.axhline(loa_upper, color="green", linestyle="--")
  1431. ax.axhline(loa_lower, color="green", linestyle="--")
  1433. x_text = x_lim[0] + 0.05*(x_lim[1] - x_lim[0])
  1436. y_offset = 0.005 * (y_lim[1] - y_lim[0])
  1438. # (3) For each line, add text "above" the line value
  1439. def place_text_above_line(line_value, label_color):
  1440. y_text = line_value
  1441. # Clamp y_text so it doesn't exceed the top/bottom
  1442. y_text = max(min(y_text, y_lim[1] - 0.001), y_lim[0] + 0.0)
  1443. ax.text(x_text, y_text, f"{line_value:.4f}", color=label_color,
  1444. va="bottom", ha="left",)
  1446. place_text_above_line(mean_diff, "red") # Mean difference
  1447. place_text_above_line(loa_upper, "green") # +1.96 SD
  1448. place_text_above_line(loa_lower, "green") # -1.96 SD
  1450. # Labels, limits, etc.
  1451. ax.set_xlabel(r"Mean of $\chi^{-}_{\text{3 T}}$ and $\chi^{-}_{\text{7 T}} \:(ppm)$")
  1452. ax.set_ylabel(r"Difference $\chi^{-}_{\text{3 T}} - \chi^{-}_{\text{7 T}} \:(ppm)$")
  1453. ax.set_title(method_name)
  1454. ax.set_xlim(x_lim)
  1455. ax.set_ylim(y_lim)
  1457. plt.tight_layout()
  1458. plt.show()

Manuscript_Figures.ipynb at commit d009ab6, under MIT · at the source

Overview

Authors: Daniel Ridani1, Benjamin De Leener1,2,3, Eva Alonso‐Ortiz1,2,4
ORCID iDs: Daniel Ridani, Eva Alonso‐Ortiz
  1. NeuroPoly Lab, Institute of Biomedical Engineering, Polytechnique Montreal, Montreal, Quebec, Canada
  2. CHU Sainte‐Justine Research Center, Montreal, Quebec, Canada
  3. Department of Computer Engineering and Software Engineering, Polytechnique Montreal, Montreal, Quebec, Canada
  4. Department of Electrical Engineering, Polytechnique Montreal, Montreal, Quebec, Canada
Institutions: Polytechnique Montréal (Canada); Centre Hospitalier Universitaire Sainte-Justine (Canada); CHU Sainte-Justine Research Center (Canada)
Journal: Magnetic resonance in medicine, volume 96, issue 5, pages 2412-2425
Dates: received 31 December 2025; accepted 1 June 2026; published online 9 July 2026; in print November 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1002/mrm.70468 · PMID 42426959 · PMCID PMC13527262 · OpenAlex W7167935606
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism), computational (subfield)
Methods: Connectivity, Statistics, fMRI & imaging
Keywords: brain phantom, MRI simulations, susceptibility anisotropy, susceptibility separation
MeSH: Brain*, Image Processing, Computer-Assisted*, Magnetic Resonance Imaging*, Phantoms, Imaging*, Algorithms, Anisotropy, Computer Simulation, Humans, Reproducibility of Results, Signal-To-Noise Ratio (* major topic)
Topic: Advanced MRI Techniques and Applications (Radiology, Nuclear Medicine and Imaging, Medicine), according to OpenAlex
Funding: Canadian Network for Research and Innovation in Machining Technology, Natural Sciences and Engineering Research Council of Canada (3101671); Fonds de recherche du Québec; Polytechnique Montréal; Canada First Research Excellence Fund; Polytechnique Montreal; TransMedTech Institute
Citations: not cited yet (Europe PMC); 57 references in the paper

Abstract

Purpose: To create a realistic in silico brain phantom for positive and negative magnetic susceptibility that incorporates susceptibility anisotropy, enabling the evaluation of how susceptibility anisotropy influences susceptibility separation algorithm performance.

Methods: We expanded an existing QSM validation phantom by creating separate maps for positive and negative susceptibility, with the option of modeling susceptibility anisotropy. Multi‐echo gradient echo data were simulated to evaluate four susceptibility separation techniques (χ‐separation, DECOMPOSE‐QSM, APART‐QSM, and R2*‐QSM). To assess the impact of noise, simulations were performed at different SNR levels (50, 100, 200, 300).

Results: Our findings showed that the error in negative susceptibility estimates increased by up to 53% when susceptibility anisotropy was present, compared to the case without susceptibility anisotropy, with χ‐separation being the algorithm that was most sensitive to anisotropy. Robustness to noise varied across the assessed algorithms, with APART‐QSM and χ‐separation having the highest and lowest sensitivity to noise, respectively.

Conclusion: The modified phantom is open‐source and can serve as a numerical ground truth for evaluating susceptibility separation methods. Our findings emphasize the importance of incorporating susceptibility anisotropy into susceptibility separation models to improve their accuracy.

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

Repository

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

neuropoly/Susceptibility-Separation-Phantom

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: d009ab67660fd6e5368d9f9811e04e19eca1636a, 17 April 2026
Languages: MATLAB (14), Jupyter (1)
Size: 20 files, 15 scripts
Software Heritage: not archived
Found in: “Data Availability Statement”
Holds: README, license file, 1 notebook
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Tools for NIfTI and ANALYZE image (MATLAB) (11 files), Image Processing Toolbox (1 file), Parallel Computing Toolbox (1 file), Matplotlib (1 file), NiBabel (1 file), NumPy (1 file), scikit-learn (1 file), SciPy (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
17 files
At the source: github.com/neuropoly/Susceptibility-Separation-Phantom

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:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 15 scripts, each with its path and the digest of its content;
  • 5 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

Datasets cited

Data Availability Statement

The code used to create susceptibility and relaxation maps, and simulate GRE signals is available at: https://github.com/neuropoly/Susceptibility‐Separation‐Phantom (https://github.com/neuropoly/Susceptibility-Separation-Phantom), while the data used to generate the figures described in this paper can be found here: https://osf.io/9xwhz/.

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 2, 28 September 2026

  • Publisher: n/a → Wiley

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 3 authors, 4 keywords, 10 MeSH terms, 6 funders, 55 references.

Cite

This paper

Ridani, D., De Leener, B., & Alonso‐Ortiz, E. (2026). A Realistic In Silico Brain Phantom for Quantifying Susceptibility Anisotropy-Induced Error in Susceptibility Separation. Magnetic resonance in medicine, 96(5), 2412-2425. https://doi.org/10.1002/mrm.70468

BibTeX

@article{ridani2026realistic,
author = {Ridani, Daniel and De Leener, Benjamin and Alonso‐Ortiz, Eva},
title = {{A Realistic In Silico Brain Phantom for Quantifying Susceptibility Anisotropy-Induced Error in Susceptibility Separation}},
journal = {Magnetic resonance in medicine},
year = {2026},
month = jul,
volume = {96},
number = {5},
pages = {2412--2425},
publisher = {Wiley},
issn = {0740-3194},
doi = {10.1002/mrm.70468},
url = {https://doi.org/10.1002/mrm.70468},
pmid = {42426959},
pmcid = {PMC13527262}
}

RIS

TY - JOUR
AU - Ridani, Daniel
AU - De Leener, Benjamin
AU - Alonso‐Ortiz, Eva
TI - A Realistic In Silico Brain Phantom for Quantifying Susceptibility Anisotropy-Induced Error in Susceptibility Separation
T2 - Magnetic resonance in medicine
J2 - Magn Reson Med
PY - 2026
DA - 2026/07/09
VL - 96
IS - 5
SP - 2412
EP - 2425
SN - 0740-3194
PB - Wiley
DO - 10.1002/mrm.70468
UR - https://doi.org/10.1002/mrm.70468
LA - en
ER -

CSL-JSON

{
"id": "10.1002/mrm.70468",
"type": "article-journal",
"title": "A Realistic In Silico Brain Phantom for Quantifying Susceptibility Anisotropy-Induced Error in Susceptibility Separation",
"container-title": "Magnetic resonance in medicine",
"author": [
{
"family": "Ridani",
"given": "Daniel"
},
{
"family": "De Leener",
"given": "Benjamin"
},
{
"family": "Alonso‐Ortiz",
"given": "Eva"
}
],
"container-title-short": "Magn Reson Med",
"volume": "96",
"issue": "5",
"page": "2412-2425",
"DOI": "10.1002/mrm.70468",
"PMID": "42426959",
"PMCID": "PMC13527262",
"ISSN": "0740-3194",
"publisher": "Wiley",
"URL": "https://doi.org/10.1002/mrm.70468",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
9
]
]
}
}

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.1002/nbm.70349 [code]
Susceptibility Source Separation Unveils Paramagnetic and Diamagnetic Trajectories in Healthy Brains From 5 to 90 Years.
Journal: NMR in biomedicine
In common: Tools for NIfTI and ANALYZE image (MATLAB), SciPy, Matplotlib, 1 other tool, structural MRI / diffusion, 9 references
[2] doi:10.1002/nbm.70301 [code]
Reducing Variability in Deep Gray Matter QSM Using Differential ROI Referencing: A Phantom and In Vivo Evaluation.
Journal: NMR in biomedicine
In common: structural MRI / diffusion, 7 references
[3] doi:10.1038/s41586-026-10631-3 [code]
A prognostic human brain network for diffuse midline glioma.
Journal: Nature
In common: Tools for NIfTI and ANALYZE image (MATLAB), Parallel Computing Toolbox, Image Processing Toolbox, 5 other tools, 1 reference
[4] doi:10.1002/nbm.70277 [code]
Hierarchical Bayesian Modelling Improves Microstructural Parameter Mapping in Diffusion and Exchange MRI Data.
Journal: NMR in biomedicine
In common: Tools for NIfTI and ANALYZE image (MATLAB), Parallel Computing Toolbox, NiBabel, 3 other tools, computational, structural MRI / diffusion, 1 reference
[5] 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: Tools for NIfTI and ANALYZE image (MATLAB), Parallel Computing Toolbox, Image Processing Toolbox, 5 other tools
[6] doi:10.1002/mrm.70405 [code]
DeepRelaxo: Fast Mono-Exponential Magnitude Brain R2* Mapping With Reduced Echoes Using Self-Supervised Deep Learning.
Journal: Magnetic resonance in medicine
In common: NiBabel, SciPy, Matplotlib, 1 other tool, computational, 3 references
[7] doi:10.1038/s41467-026-74215-5 [code]
Multi-metric evaluations of acute psychedelic effects on fMRI brain entropy.
Journal: Nature communications
In common: Tools for NIfTI and ANALYZE image (MATLAB), Parallel Computing Toolbox, Image Processing Toolbox, 4 other tools, computational
[8] doi:10.1002/hbm.70602 [code]
Neuroimaging Correlates of Post-Stroke Pain After Ischemic Stroke: Secondary Analysis of the INSPiRE-TMS Trial.
Journal: Human brain mapping
In common: Tools for NIfTI and ANALYZE image (MATLAB), Parallel Computing Toolbox, Image Processing Toolbox, 4 other tools
[9] doi:10.1111/ene.70678 [code]
Who Falls After a Stroke? Evidence From a Prospective Stroke Cohort.
Journal: European journal of neurology
In common: Tools for NIfTI and ANALYZE image (MATLAB), Parallel Computing Toolbox, Image Processing Toolbox, 4 other tools
[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: Tools for NIfTI and ANALYZE image (MATLAB), Parallel Computing Toolbox, Image Processing Toolbox, 4 other tools

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.