OSCR

Brain MRI Dataset Featuring a Full Clinical Protocol With and Without Intentional Motion.

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
  1. [1] § Technical Validation › Registration and Analysis ↔ MotionCorrectedClinicalMRProtocol/analysis_img_quality.py, lines 15–74 · score 0.71 · FreeSurfer, uncorrected scan, ground truth, T1 MPR, motion correction, prospective
  2. [2] § Technical Validation › Observer scores ↔ MotionCorrectedClinicalMRProtocol/generate_plots_manuscript.py, lines 582–650 · score 0.63 · T1 STIR, T2 TSE, T2 FLAIR, 1–5, T1 MPR, motion
  3. [3] § Methods › Clinical MRI protocol ↔ MotionCorrectedClinicalMRProtocol/generate_plots_manuscript.py, lines 582–650 · score 0.61 · T1 STIR, T2 TSE, T2 FLAIR, T1 MPR, DWI, MR
  4. [4] § Data Records ↔ MotionCorrectedClinicalMRProtocol/recon_register.py, lines 8–88 · score 0.58 · T1 TIRM, T2 TSE, FreeSurfer, MPRAGE, FLAIR, motion correction
  5. [5] § Background & Summary ↔ MotionCorrectedClinicalMRProtocol/generate_plots_manuscript.py, lines 527–580 · score 0.55 · image quality metrics, observer scores, retrospective, PMC, prospective, motion correction

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,187 lines · 47 KB · MIT · 3 matches

  1. import numpy as np
  2. import matplotlib
  3. import matplotlib.pyplot as plt
  4. import glob
  5. import os
  6. import string
  7. import nibabel as nib
  8. from statsmodels.stats.multitest import multipletests
  9. from utils import Show_Stars, SortFiles, MakeBoxplot, DrawLines2, DrawLines
  10. from statistical_tests import PerformWilcoxonMotion, PerformWilcoxonAllImg
  11. plt.style.use('seaborn-whitegrid')
  12. matplotlib.rc('axes', edgecolor='black')
  13. from matplotlib.ticker import ScalarFormatter
  15. out_dir = '../Results/Plots/'
  16. in_dir_mot = '../Results/Motion_Estimates/'
  17. in_dir_met = '../Results/Metrics_Results/'
  18. in_dir_qs = '../ObserverQualityScores/'
  19. out_dir_metrics = '../Results/Metrics_Results/Comparison/'
  21. subdir = []
  22. for i in range(1,10):
  23. subdir.append('Subject_0'+str(i)+'/')
  24. for i in range(10,20):
  25. subdir.append('Subject_'+str(i)+'/')
  26. for i in range(20,23):
  27. subdir.append('Subject_'+str(i)+'/')
  28. sequs = ['T1_MPR', 'T2_FLAIR', 'T2_TSE', 'T1_TIRM', 'T2STAR', 'DIFF']
  30. save = '_2022_05_27'
  32. plot_motion = True
  33. plot_still = True
  34. plot_nod = True
  35. plot_DWI = True
  36. calc_ADC_hist = True
  38. ''' (1) Plot motion data: '''
  39. if plot_motion:
  41. # calculate p-values - first STILL, then NOD, then SHAKE, in each
  42. # cathegory: RMS, then median, then maximum
  43. p_values = []
  44. effect_size = []
  45. for mot in ['STILL', 'NOD', 'SHAKE']:
  46. RMS, median, maxim, descr = [], [], [], []
  47. for sequ in sequs:
  48. if mot == 'SHAKE' and sequ != 'T1_MPR':
  49. continue
  50. # find the most recent data:
  51. file = glob.glob(in_dir_mot+'MotionMetrics_'+mot+'/'+sequ+'*.txt')
  52. file = [f for f in file if 'mid' not in f] #sort out 'mid'
  53. if len(file)>0:
  54. file = SortFiles(file)
  56. metrics = np.loadtxt(file[0], unpack=True, skiprows=1)
  57. RMS.append(metrics[0])
  58. RMS.append(metrics[3])
  59. median.append(metrics[1])
  60. median.append(metrics[4])
  61. maxim.append(metrics[2])
  62. maxim.append(metrics[5])
  63. descr.append(sequ)
  64. descr.append(sequ)
  66. p_val, ind, alt, ES = PerformWilcoxonMotion(['RMS', 'Med', 'Max'], mot,
  67. [RMS, median, maxim],
  68. in_dir_mot, save,)
  69. p_values.append(p_val)
  70. effect_size.append(ES)
  72. #correct for multiple comparisons:
  73. p_values = np.concatenate(p_values)
  74. rej, p_values_cor, _, __ = multipletests(p_values, alpha=0.05,
  75. method='fdr_bh', is_sorted=False,
  76. returnsorted=False)
  78. p_values_cor_still = p_values_cor[0:18]
  79. p_values_cor_nod = p_values_cor[18:36]
  80. p_values_cor_shake = p_values_cor[36:]
  81. ind_p = np.array([[0,1], [2,3], [4,5], [6,7], [8,9], [10,11]])
  82. ind_sh = np.array([[0,1]])
  85. # plot:
  86. num = 0
  87. plt.figure(figsize=(13,10))
  88. for mot in ['STILL', 'NOD']:
  89. RMS, median, maxim, descr = [], [], [], []
  90. for sequ in sequs:
  91. # find the most recent data:
  92. file = glob.glob(in_dir_mot+'MotionMetrics_'+mot+'/'+sequ+'*.txt')
  93. file = [f for f in file if 'mid' not in f] # sort out mid
  94. if len(file)>0:
  95. file = SortFiles(file)
  97. # load the most recent file:
  98. metrics = np.loadtxt(file[0], unpack=True, skiprows=1)
  99. RMS.append(metrics[0])
  100. RMS.append(metrics[3])
  101. median.append(metrics[1])
  102. median.append(metrics[4])
  103. maxim.append(metrics[2])
  104. maxim.append(metrics[5])
  105. descr.append(sequ)
  106. descr.append(sequ)
  108. if mot == 'NOD':
  109. mot_s = 'SHAKE'
  110. RMS_s, median_s, maxim_s, descr_s = [], [], [], []
  111. sequ = 'T1_MPR'
  112. # find the most recent data:
  113. file = glob.glob(in_dir_mot+'MotionMetrics_'+mot_s+'/'+sequ+'*.txt')
  114. file = [f for f in file if 'mid' not in f] # sort out mid
  115. if len(file)>0:
  116. file = SortFiles(file)
  118. # load the most recent file
  119. metrics = np.loadtxt(file[0], unpack=True, skiprows=1)
  120. RMS_s.append(metrics[0])
  121. RMS_s.append(metrics[3])
  122. median_s.append(metrics[1])
  123. median_s.append(metrics[4])
  124. maxim_s.append(metrics[2])
  125. maxim_s.append(metrics[5])
  126. descr_s.append(sequ)
  127. descr_s.append(sequ)
  128. mean_RMS_s, mean_med_s, mean_max_s = [], [], []
  129. for i in range(len(RMS_s)):
  130. mean_RMS_s.append(np.mean(RMS_s[i]))
  131. mean_med_s.append(np.mean(median_s[i]))
  132. mean_max_s.append(np.mean(maxim_s[i]))
  134. # calculate mean values and define colors
  135. mean_RMS, mean_med, mean_max = [], [], []
  136. for i in range(len(RMS)):
  137. mean_RMS.append(np.mean(RMS[i]))
  138. mean_med.append(np.mean(median[i]))
  139. mean_max.append(np.mean(maxim[i]))
  140. colors = []
  141. labels = []
  142. for i in range(int(len(RMS)/2)):
  143. colors.append('tab:orange')
  144. colors.append('tab:blue')
  145. labels.append('without PMC')
  146. labels.append('with PMC')
  147. small = dict(markersize=3)
  149. # change 'TIRM' into STIR for descr:
  150. # and change 'TRACEW_B0' to 'DWI'
  151. for i in range(len(descr)):
  152. if descr[i]=='TRACEW_B0':
  153. descr[i]='DWI'
  154. if descr[i]=='T1_TIRM':
  155. descr[i]='T1_STIR'
  156. if descr[i]=='T2STAR':
  157. descr[i]='T2*'
  159. # plot the data:
  160. if mot=='STILL':
  161. ax0=plt.subplot2grid((4,5), (0,0), colspan=4)
  162. ax = ax0
  163. else:
  164. ax1=plt.subplot2grid((4,5), (2,0), colspan=4)
  165. ax = ax1
  166. MakeBoxplot(RMS, colors)
  167. for i in range(len(mean_RMS)):
  168. plt.plot(i+1, mean_RMS[i], '.', c=colors[i], ls='')
  169. if mot=='STILL':
  170. for y1, y2 in zip(RMS[4], RMS[5]):
  171. plt.plot([5,6], [y1, y2], 'darkslategray', lw=0.7)
  172. for y1, y2 in zip(RMS[8], RMS[9]):
  173. plt.plot([9,10], [y1, y2], 'darkslategray', lw=0.7)
  175. plt.title(mot+' scans', fontsize=17)
  176. plt.ylabel('RMS displ [mm]')
  177. plt.xticks(ticks=np.arange(1, len(RMS)+1), labels=descr)
  178. lim = plt.gca().get_ylim()
  179. plt.ylim(lim[0],(lim[1]-lim[0])*1.2+lim[0])
  180. ax.text(-0.1, .95, string.ascii_lowercase[num], transform=ax.transAxes,
  181. size=20, weight='bold')
  182. num += 1
  183. if mot=='STILL':
  184. use = 0
  185. low = -1
  186. if mot == 'NOD':
  187. use = 1
  188. low = -4
  190. max_RMS = [np.amax(RMS[i]) for i in range(len(RMS))]
  191. if mot=='STILL':
  192. Show_Stars(p_values_cor_still[0:6], ind_p, np.arange(1, len(RMS)+1),
  193. max_RMS, dh=0.05, col='black')
  195. ax01=plt.subplot2grid((4,5), (1,0), colspan=4)
  196. ax = ax01
  197. else:
  198. Show_Stars(p_values_cor_nod[0:6], ind_p, np.arange(1, len(RMS)+1),
  199. max_RMS, col='black')
  201. ax2=plt.subplot2grid((4,5), (3,0), colspan=4)
  202. ax = ax2
  203. MakeBoxplot(maxim, colors)
  204. for i in range(0,2):
  205. plt.plot(i+1, mean_max[i], '.', c=colors[i], ls='',
  206. label=labels[i])
  207. for i in range(2,len(mean_RMS)):
  208. plt.plot(i+1, mean_max[i], '.', c=colors[i], ls='')
  209. if mot=='STILL':
  210. for y1, y2 in zip(maxim[4], maxim[5]):
  211. plt.plot([5,6], [y1, y2], 'darkslategray', lw=0.7)
  212. for y1, y2 in zip(maxim[8], maxim[9]):
  213. plt.plot([9,10], [y1, y2], 'darkslategray', lw=0.7)
  214. plt.ylabel('Maximum displ [mm]')
  215. plt.xticks(ticks=np.arange(1, len(RMS)+1), labels=descr)
  216. lim = plt.gca().get_ylim()
  217. plt.ylim(lim[0],(lim[1]-lim[0])*1.1+lim[0])
  218. ax.text(-0.1, .95, string.ascii_lowercase[num], transform=ax.transAxes,
  219. size=20, weight='bold')
  220. num += 1
  221. if mot=='STILL':
  222. use = 0
  223. low = -3
  224. max_max = [np.amax(maxim[i]) for i in range(len(RMS))]
  225. if mot=='STILL':
  226. Show_Stars(p_values_cor_still[12:], ind_p, np.arange(1, len(RMS)+1),
  227. max_max, col='black')
  228. legend = plt.legend(loc='upper center', bbox_to_anchor=(1.16, 1.3),
  229. fontsize=13, frameon=True)
  230. else:
  231. Show_Stars(p_values_cor_nod[12:], ind_p, np.arange(1, len(RMS)+1),
  232. max_max, col='black')
  234. if mot == 'NOD':
  235. ax4=plt.subplot2grid((4,5), (2,4))
  236. ax = ax4
  237. ax1.get_shared_y_axes().join(ax1, ax4)
  238. MakeBoxplot(RMS_s, colors)
  239. for i in range(len(mean_RMS_s)):
  240. plt.plot(i+1, mean_RMS_s[i], '.', c=colors[i], ls='')
  241. plt.title(mot_s+' scans', fontsize=17)
  242. plt.xticks(ticks=np.arange(1, len(RMS_s)+1), labels=descr_s)
  243. lim = plt.gca().get_ylim()
  244. plt.ylim(lim[0],(lim[1]-lim[0])*1.1+lim[0])
  245. ax.text(-0.3, .95, string.ascii_lowercase[num],
  246. transform=ax.transAxes, size=20, weight='bold')
  247. num += 1
  248. max_RMS = [np.amax(RMS_s[i]) for i in range(len(RMS_s))]
  249. Show_Stars(np.array([p_values_cor_shake[0]]), ind_sh,
  250. np.arange(1, len(RMS)+1), max_RMS, col='black')
  252. ax5=plt.subplot2grid((4,5), (3,4))
  253. ax = ax5
  254. ax2.get_shared_y_axes().join(ax2, ax5)
  256. MakeBoxplot(maxim_s, colors)
  257. for i in range(0,2):
  258. plt.plot(i+1, mean_max_s[i], '.', c=colors[i], ls='',
  259. label=labels[i])
  260. for i in range(2,len(mean_RMS_s)):
  261. plt.plot(i+1, mean_max_s[i], '.', c=colors[i], ls='')
  262. plt.xticks(ticks=np.arange(1, len(RMS_s)+1), labels=descr_s)
  263. lim = plt.gca().get_ylim()
  264. plt.ylim(lim[0],(lim[1]-lim[0])*1.1+lim[0])
  265. ax.text(-0.3, .95, string.ascii_lowercase[num],
  266. transform=ax.transAxes, size=20, weight='bold')
  267. num += 1
  268. max_max = [np.amax(maxim[i]) for i in range(len(RMS_s))]
  269. Show_Stars(np.array([p_values_cor_shake[2]]), ind_sh,
  270. np.arange(1, len(RMS)+1), max_max, col='black')
  272. legend.get_frame().set_linewidth(2)
  273. plt.subplots_adjust(hspace=0.4, wspace=0.4)
  274. plt.savefig(out_dir+'Boxplot_motion_'+save+'.tiff', format='tiff',
  275. bbox_inches='tight', dpi=200)
  276. plt.savefig(out_dir+'Boxplot_motion_'+save+'.png', format='png',
  277. bbox_inches='tight', dpi=200)
  278. plt.show()
  282. ''' (2) Import image quality metrics: '''
  283. withRR = True
  284. onlyRR = False
  285. quality_scores = True
  286. show_stat_test = True
  288. sequs = ['T1_MPR', 'T2_FLAIR', 'T2_TSE', 'T1_TIRM', 'T2STAR', 'TRACEW_B1000']
  290. ssims, psnrs, tgs, qss, names_tes = [], [], [], [], []
  291. p_ssims, p_psnrs, p_tgs, p_qss, ind_ps = [], [], [], [], []
  292. for sequ in sequs:
  293. n = len(subdir)
  294. if sequ == 'T1_MPR':
  295. m = 10
  296. names = np.array(['MOCO_ON_STILL', 'MOCO_ON_SHAKE_RR', 'MOCO_ON_SHAKE',
  297. 'MOCO_ON_NOD_RR', 'MOCO_ON_NOD', 'MOCO_OFF_STILL',
  298. 'MOCO_OFF_SHAKE_RR', 'MOCO_OFF_SHAKE',
  299. 'MOCO_OFF_NOD_RR', 'MOCO_OFF_NOD'])
  300. else:
  301. m = 6
  302. names = np.array(['MOCO_ON_STILL', 'MOCO_ON_NOD_RR', 'MOCO_ON_NOD',
  303. 'MOCO_OFF_STILL', 'MOCO_OFF_NOD_RR', 'MOCO_OFF_NOD'])
  304. ssim, psnr, tg_, qs = np.zeros((n,m)), np.zeros((n,m)), np.zeros((n,m)), np.zeros((n,m))
  305. i=0
  306. print('Files used for calculation:')
  307. for sub in subdir:
  308. # skip all volunteers where the following sequences were not acquired
  309. if sequ in ['ADC', 'TRACEW_B0', 'TRACEW_B1000', 'T2_FLAIR', 'T2STAR']:
  310. tmp = os.listdir(in_dir_met+sub)
  311. tmp_ = ''
  312. if sequ not in tmp_.join(tmp):
  313. continue
  316. folder = in_dir_met+sub
  317. file_ = glob.glob(folder+'Values_*'+sequ+'*')
  318. file = [f for f in file_ if f.find('Retro')==-1]
  320. if len(file)>0:
  321. # get the most recent file:
  322. file = SortFiles(file, True)
  323. print(file[0])
  324. tmp1 = np.loadtxt(file[0], unpack=True, dtype=str, usecols=0,
  325. skiprows=1)
  326. tmp2, tmp3, tmp4 = np.loadtxt(file[0], unpack=True,
  327. usecols=(1,2,3), skiprows=1)
  329. if sequ == 'T2STAR' and np.size(tmp1) == 1:
  330. # A few subject only contain a MOCO OFF STILL T2STAR scan.
  331. # Those were only acquired for comparison with susceptibility
  332. # weighted scans and are not included in analysis of image
  333. # quality metrics.
  334. continue
  336. incl = []
  337. for na,s,p,t in zip(tmp1, tmp2, tmp3, tmp4):
  338. if sequ in ['T2STAR', 'ADC', 'TRACEW_B0', 'TRACEW_B1000']:
  339. if 'NOD' in na:
  340. na = na+'_RR'
  342. ind = np.where(names==na)[0][0]
  343. ssim[i,ind], psnr[i,ind], tg_[i,ind] = s, p, t
  344. incl.append(ind)
  345. # 0 entires represent non-existing scans
  347. # do the same for quality scores:
  348. if sequ not in ['T2STAR', 'ADC', 'TRACEW_B0', 'TRACEW_B1000']:
  349. file = glob.glob(in_dir_qs+sequ+' Score.txt')
  350. if len(file)>0:
  351. # get the most recent file:
  352. file = SortFiles(file)
  353. print(file[0])
  355. subj_names = np.loadtxt(file[0], unpack=True, dtype=str, usecols=0,
  356. skiprows=1)
  357. tmp1 = np.loadtxt(file[0], unpack=True, dtype=str, usecols=1,
  358. skiprows=1)
  359. tmp2, tmp3, tmp4 = np.loadtxt(file[0], unpack=True, usecols=(2,3,4),
  360. skiprows=1)
  362. tmp1 = tmp1[subj_names==sub[:-1]]
  363. tmp2 = tmp2[subj_names==sub[:-1]]
  364. tmp3 = tmp3[subj_names==sub[:-1]]
  365. tmp4 = tmp4[subj_names==sub[:-1]]
  368. incl = []
  369. for na,q1,q2,q3 in zip(tmp1, tmp2, tmp3, tmp4):
  370. if 'RETRO' not in na:
  371. ind = np.where(names==na[:-1])[0][0]
  372. # average the scores of the three raters with double weight for the radiologist
  373. qs[i,ind] = (q1+q2+2*q3)/4
  374. incl.append(ind)
  378. # for DWI quality ranks instead of quality scores:
  379. elif sequ == 'TRACEW_B1000':
  380. file = glob.glob(in_dir_qs+'DWI Rank.txt')
  381. if len(file)>0:
  382. # get the most recent file:
  383. file = SortFiles(file)
  384. print(file[0])
  386. subj_names = np.loadtxt(file[0], unpack=True, dtype=str, usecols=0,
  387. skiprows=1)
  388. tmp1 = np.loadtxt(file[0], unpack=True, dtype=str, usecols=1,
  389. skiprows=1)
  390. tmp2 = np.loadtxt(file[0], unpack=True, usecols=(2),
  391. skiprows=1)
  393. tmp1 = tmp1[subj_names==sub[:-1]]
  394. tmp2 = tmp2[subj_names==sub[:-1]]
  396. incl = []
  397. for na,q in zip(tmp1, tmp2):
  398. if 'NOD' in na:
  399. na = na+'_RR'
  400. if 'RETRO' not in na:
  401. ind = np.where(names==na)[0][0]
  402. qs[i,ind] = q
  403. incl.append(ind)
  404. else:
  405. qs = np.zeros_like(ssim)
  406. i+=1
  408. names = np.tile(names, (n,1))
  410. # check that names are the same for all volunteers:
  411. for i in range(1, len(subdir)):
  412. if np.all(names[0]==names[i], axis=0)==False:
  413. print('ERROR: Names of 0 and '+str(i)+' are not the same!')
  415. names = names[0]
  417. ''' sort out values == 0 (corresponding to missing scans) '''
  418. ssim[ssim==0] = np.nan
  419. psnr[psnr==0] = np.nan
  420. tg_[tg_==0] = np.nan
  421. qs[qs==0] = np.nan
  424. ''' perform parametric tests:'''
  425. # resort names and data (only for test), final resorting will be performed
  426. # after RR scans are potentially sorted out:
  427. names_ch = []
  428. for n in names:
  429. ind=n.find('_', 5)
  430. if 'OFF' in n:
  431. tmp = 'C'
  432. elif 'ON' in n:
  433. tmp = 'A'
  434. elif 'RETRO' in n:
  435. tmp = 'B'
  436. names_ch.append(n[ind+1:]+'_'+tmp)
  438. names_ch = np.array(names_ch)
  439. ind = np.argsort(names_ch)[::-1]
  441. ssim_p, psnr_p, tg_p, qs_p = ssim[:,ind], psnr[:,ind], tg_[:,ind], qs[:,ind]
  442. names_p = names_ch[ind]
  444. p_ssim, rej_ssim, ind_p, alt = PerformWilcoxonAllImg('SSIM', ssim_p, sequ,
  445. out_dir_metrics, save)
  446. p_psnr, rej_psnr, ind_p, alt = PerformWilcoxonAllImg('PSNR', psnr_p, sequ,
  447. out_dir_metrics, save)
  448. p_tg, rej_tg, ind_p, alt = PerformWilcoxonAllImg('TG', tg_p, sequ,
  449. out_dir_metrics, save)
  450. p_qs, rej_qs, ind_p, alt = PerformWilcoxonAllImg('QS', qs_p, sequ,
  451. out_dir_metrics, save)
  453. # sort out the relevant tests:
  454. if onlyRR == True:
  455. if sequ == 'T1_MPR':
  456. rel_tests = [0,1,4]
  457. p_ssim, p_psnr = p_ssim[rel_tests], p_psnr[rel_tests]
  458. p_tg, p_qs = p_tg[rel_tests], p_qs[rel_tests]
  459. ind_p = np.array([[0,1], [2,3], [4,5]])
  460. else:
  461. rel_tests = [0,1]
  462. p_ssim, p_psnr = p_ssim[rel_tests], p_psnr[rel_tests]
  463. p_tg, p_qs = p_tg[rel_tests], p_qs[rel_tests]
  464. ind_p = np.array([[0,1], [2,3]])
  466. if withRR == False:
  467. findRRs = np.char.find(names, 'RR')
  468. findRRs = (findRRs<0)
  469. names = names[findRRs]
  470. ssim, psnr, tg_ = ssim[findRRs], psnr[findRRs], tg_[findRRs]
  471. ssim = np.reshape(ssim, (len(subdir), int(len(ssim)/len(subdir))))
  472. psnr = np.reshape(psnr, (len(subdir), int(len(psnr)/len(subdir))))
  473. tg_ = np.reshape(tg_, (len(subdir), int(len(tg_)/len(subdir))))
  474. qs = np.reshape(qs, (len(subdir), int(len(qs)/len(subdir))))
  475. names = np.reshape(names, (len(subdir), int(len(names)/len(subdir))))
  477. if onlyRR == True:
  478. findRRs = np.char.find(names, 'RR')
  479. findstills = np.char.find(names, 'STILL')
  480. both = findRRs*findstills
  481. both = (both<0)
  482. names = names[both]
  483. ssim, psnr, tg_, qs = ssim[:,both], psnr[:,both], tg_[:,both], qs[:,both]
  486. std_ssim, std_psnr = np.nanstd(ssim, axis=0), np.nanstd(psnr, axis=0)
  487. std_tg_, std_qs = np.nanstd(tg_, axis=0), np.nanstd(qs, axis=0)
  488. mean_ssim, mean_psnr = np.nanmean(ssim, axis=0), np.nanmean(psnr, axis=0)
  489. mean_tg_, mean_qs = np.nanmean(tg_, axis=0), np.nanmean(qs, axis=0)
  491. # Final resorting:
  492. names_ch = []
  493. for n in names:
  494. ind=n.find('_', 5)
  495. if 'OFF' in n:
  496. tmp = 'C'
  497. elif 'ON' in n:
  498. tmp = 'A'
  499. elif 'RETRO' in n:
  500. tmp = 'B'
  501. names_ch.append(n[ind+1:]+'_'+tmp)
  503. names_ch = np.array(names_ch)
  504. ind = np.argsort(names_ch)[::-1]
  505. ssim, psnr, tg_, qs = ssim[:,ind], psnr[:,ind], tg_[:,ind], qs[:,ind]
  506. names = names[ind]
  507. mean_ssim, mean_psnr = mean_ssim[ind], mean_psnr[ind]
  508. mean_tg_, mean_qs = mean_tg_[ind], mean_qs[ind]
  510. names_te = []
  511. for i in range(int(len(names))):
  512. tmp = names[i]
  513. if 'RETRO' in tmp:
  514. app = 'retrospective MoCo'
  515. elif 'ON' in tmp:
  516. app = 'prospective MoCo'
  517. elif 'OFF' in tmp:
  518. app = 'no MoCo'
  519. if onlyRR == False:
  520. if 'RR' in tmp or 'STILL' in tmp:
  521. app = app + ' no REAC'
  522. else:
  523. app = app + ' with REAC'
  524. names_te.append(app)
  527. # sort out nans (whole rows for subjects wihtout that sequence)
  528. mask = np.where(np.isnan(ssim[:,0])==False)[0]
  529. ssim = ssim[mask]
  530. mask = np.where(np.isnan(psnr[:,0])==False)[0]
  531. psnr = psnr[mask]
  532. mask = np.where(np.isnan(tg_[:,0])==False)[0]
  533. tg_ = tg_[mask]
  534. mask = np.where(np.isnan(qs[:,0])==False)[0]
  535. qs = qs[mask]
  537. # normalise tg values by dividing with mean value of off still
  538. val = np.mean(tg_[:,0])
  539. tg_ = tg_/val
  542. # new color code:
  543. color_dict = {'retrospective MoCo': 'tab:green',
  544. 'prospective MoCo': 'tab:blue', 'no MoCo':'tab:orange'}
  545. if onlyRR == False:
  546. color_dict = {'retrospective MoCo no REAC': 'tab:green',
  547. 'prospective MoCo no REAC': 'tab:blue',
  548. 'no MoCo no REAC':'tab:orange',
  549. 'retrospective MoCo with REAC': 'tab:gray',
  550. 'prospective MoCo with REAC': 'tab:olive',
  551. 'no MoCo with REAC':'tab:cyan'}
  553. ssims.append(ssim)
  554. psnrs.append(psnr)
  555. tgs.append(tg_)
  556. qss.append(qs)
  557. names_tes.append(names_te)
  558. p_ssims.append(p_ssim)
  559. p_psnrs.append(p_psnr)
  560. p_tgs.append(p_tg)
  561. p_qss.append(p_qs)
  562. ind_ps.append(ind_p)
  566. ''' (3) Plot image quality metrics for still scans: '''
  567. if plot_still:
  568. metrics = [tgs, qss, np.array([qss[-1]])]
  569. p_values = [p_tgs, p_qss[0:4], np.array([p_qss[-1]])]
  570. labels = ['Tenengrad', 'Observer Scores', 'Quality Rank']
  571. colors = ['tab:orange', 'tab:blue', 'tab:orange', 'tab:blue', 'tab:orange',
  572. 'tab:blue', 'tab:orange', 'tab:blue', 'tab:orange', 'tab:blue',
  573. 'tab:orange', 'tab:blue']
  574. fig_labels = ['without PMC', 'with PMC']
  575. small = dict(markersize=3)
  576. x = np.arange(1,13)
  577. a = [0,2,4,6,8,10]
  578. b = [1,3,5,7,9,11]
  579. c = [1,3,5,7,9,11]
  580. d = [2,4,6,8,10,12]
  582. num = 0
  583. plt.figure(figsize=(10,9))
  584. for i, metric, p, lab in zip(range(0,3), metrics, p_values, labels):
  585. m_still = []
  586. means = []
  587. for m in metric:
  588. if len(m)>0:
  589. m_still.append(m[:,0])
  590. m_still.append(m[:,1])
  591. means.append(np.nanmean(m[:,0]))
  592. means.append(np.nanmean(m[:,1]))
  593. p_still = []
  594. for pval in p:
  595. p_still.append(pval[0])
  597. if i == 0:
  598. ax = plt.subplot2grid((2,6), (0,0), colspan=6)
  599. box1 = plt.boxplot(m_still, flierprops=small)
  600. for j in range(len(means)):
  601. plt.errorbar(x[j], means[j], yerr=None, color=colors[j], fmt='.', capsize=3)
  602. ticklabels = ['T1_MPR', 'T2_FLAIR', 'T2_TSE', 'T1_STIR', 'T2*', 'TRACEW']
  603. ticks = [1.5, 3.5, 5.5, 7.5, 9.5, 11.5]
  604. plt.xticks(labels=ticklabels, ticks=ticks, fontsize=14)
  606. elif i == 1:
  607. ax = plt.subplot2grid((2,6), (1,0), colspan=4)
  608. box1 = plt.boxplot(m_still[0:8], flierprops=small)
  609. for j in range(len(means)-2):
  610. plt.errorbar(x[j], means[j], yerr=None, color=colors[j], fmt='.', capsize=3)
  611. for j in range(0,2):
  612. plt.errorbar(x[j], means[j], yerr=None, color=colors[j], fmt='.',
  613. capsize=3, label=fig_labels[j])
  614. ticklabels = ['T1_MPR', 'T2_FLAIR', 'T2_TSE', 'T1_STIR']
  615. ticks = [1.5, 3.5, 5.5, 7.5]
  616. plt.xticks(labels=ticklabels, ticks=ticks, fontsize=14)
  618. else:
  619. ax = plt.subplot2grid((2,6), (1,5), colspan=1)
  620. box1 = plt.boxplot(m_still, flierprops=small)
  621. for j in range(len(means)):
  622. plt.errorbar(x[j], means[j], yerr=None, color=colors[j], fmt='.', capsize=3)
  623. ticklabels = ['DWI']
  624. ticks = [1.5]
  625. plt.xticks(labels=ticklabels, ticks=ticks, fontsize=14)
  626. plt.yticks([2,3,4])
  628. for patch, patch2, color in zip(box1['boxes'], box1['medians'], colors):
  629. patch.set(color=color, lw=1.7)
  630. patch2.set(color='k', lw=1.7)
  632. plt.ylabel(lab, fontsize=15)
  633. if show_stat_test == True:
  634. maxi = []
  635. for v in m_still:
  636. maxi.append(np.amax(v))
  637. if i == 0:
  638. indices = [[0,1], [2,3], [4,5], [6,7], [8,9], [10,11]]
  639. Show_Stars(np.array(p_still), indices[0:len(p_still)], x, maxi,
  640. col='black')
  641. elif i ==1:
  642. indices = [[0,1], [2,3], [4,5], [6,7]]
  643. Show_Stars(np.array(p_still), indices[0:len(p_still)], x[0:8], maxi[0:8],
  644. col='black')
  645. else:
  646. indices = [[0,1]]
  647. Show_Stars(np.array(p_still), indices[0:len(p_still)], x, maxi,
  648. col='black')
  649. lim = plt.gca().get_ylim()
  650. plt.ylim(lim[0],(lim[1]-lim[0])*1.05+lim[0])
  652. DrawLines2(a[0:len(p_still)],b[0:len(p_still)],c[0:len(p_still)],
  653. d[0:len(p_still)],m_still, lw=0.7, col='darkslategray')
  655. if i == 2:
  656. ax.text(-0.6, 0.95, string.ascii_lowercase[num],
  657. transform=ax.transAxes, size=24, weight='bold')
  658. else:
  659. ax.text(-0.1, 0.95, string.ascii_lowercase[num],
  660. transform=ax.transAxes, size=24, weight='bold')
  661. num += 1
  663. plt.yticks(fontsize=13)
  664. plt.tick_params('both', length=0)
  665. plt.gca().yaxis.set_major_formatter(ScalarFormatter(useOffset=False))
  666. if i == 0:
  667. xlim = plt.gca().get_xlim()
  668. if i == 1:
  669. legend = plt.legend( loc='lower left', ncol=2,
  670. bbox_to_anchor=(0.4, -0.3), fontsize=14,
  671. frameon=True)
  672. plt.yticks(ticks=[2.5, 3, 3.5, 4, 4.5, 5])
  674. legend.get_frame().set_linewidth(2)
  675. plt.subplots_adjust(hspace=0.2, wspace=0.3)
  676. plt.savefig(out_dir+'Metrics_still'+save+'.tiff', format='tiff',
  677. bbox_inches='tight', dpi=200)
  678. plt.savefig(out_dir+'Metrics_still'+save+'.png', format='png',
  679. bbox_inches='tight', dpi=200)
  680. plt.show()
  684. ''' (4) Plot image quality metrics for nodding scans: '''
  685. if plot_nod:
  686. # first MPRAGE and FLAIR
  687. metrics = [ssims, psnrs, tgs, qss]
  688. p_values = [p_ssims, p_psnrs, p_tgs, p_qss[0:4]]
  689. labels = ['SSIM', 'PSNR', 'Tenengrad', 'Observer Scores']
  690. color_dict = {'prospective MoCo no REAC': 'tab:blue',
  691. 'no MoCo no REAC':'tab:orange',
  692. 'prospective MoCo with REAC': 'tab:green',
  693. 'no MoCo with REAC':'tab:cyan'}
  694. name_dict = {'prospective MoCo no REAC': 'with PMC without REAC',
  695. 'no MoCo no REAC':'without PMC without REAC',
  696. 'prospective MoCo with REAC': 'with PMC with REAC',
  697. 'no MoCo with REAC':'without PMC with REAC'}
  699. small = dict(markersize=3)
  700. x = np.concatenate((np.arange(1,9), np.arange(10,14)), axis=None)
  701. a = [0,2,4,6]
  702. b = [1,3,5,7]
  703. c = [1,3,5,7]
  704. d = [2,4,6,8]
  705. a_ = [0,2]
  706. b_ = [1,3]
  707. c_ = [10,12]
  708. d_ = [11,13]
  710. num = 0
  711. plt.figure(figsize=(11,8))
  712. for i, metric, p, lab in zip(range(0,4), metrics, p_values, labels):
  713. means = []
  714. names = []
  715. # first MPR
  716. tmp = metric[0]
  717. m_mpr = tmp[:,2:]
  718. names = names_tes[0][2:]
  720. # then FLAIR
  721. tmp = metric[1]
  722. m_fl = tmp[:,2:]
  723. names = np.concatenate((names, names_tes[1][2:]), axis=None)
  725. p_mot = []
  726. indices = []
  727. for pval, ind in zip(p[0:2], ind_ps[0:2]):
  728. p_mot.append(pval[1:])
  729. indices.append(ind[1:])
  731. means = np.concatenate((np.nanmean(m_mpr, axis=0),
  732. np.nanmean(m_fl, axis=0)), axis=None)
  734. colors = []
  735. leg_labels = []
  736. for n in names:
  737. colors.append(color_dict[n])
  738. leg_labels.append(name_dict[n])
  740. ax = plt.subplot(2,2,i+1)
  741. box1 = plt.boxplot(m_mpr, flierprops=small, widths=0.5)
  742. for patch, patch2, color in zip(box1['boxes'], box1['medians'], colors):
  743. patch.set(color=color, lw=1.5)
  744. patch2.set(color='k', lw=1.5)
  745. box2 = plt.boxplot(m_fl, positions=range(10,14), flierprops=small,
  746. widths=0.5)
  747. for patch, patch2, color in zip(box2['boxes'], box2['medians'], colors):
  748. patch.set(color=color, lw=1.5)
  749. patch2.set(color='k', lw=1.5)
  750. for j in range(len(means)):
  751. plt.errorbar(x[j], means[j], yerr=None, color=colors[j], fmt='.',
  752. capsize=3)
  753. if i == 2:
  754. for j in range(0,4):
  755. plt.errorbar(x[j], means[j], yerr=None, color=colors[j],
  756. fmt='.', capsize=3, label=leg_labels[j])
  757. legend = plt.legend( loc='lower left', ncol=2,
  758. bbox_to_anchor=(0.25, -0.4), fontsize=12,
  759. frameon=True)
  760. plt.ylabel(lab, fontsize=15)
  761. if show_stat_test == True:
  762. maxi = []
  763. for v in m_mpr.T:
  764. maxi.append(np.amax(v))
  765. Show_Stars(np.array(p_mot[0]), indices[0]-2, x, maxi,
  766. arange_dh='PAPER', col='black')
  767. for v in m_fl.T:
  768. maxi.append(np.amax(v))
  769. Show_Stars(np.array(p_mot[1]), indices[1]-2, x[8:], maxi,
  770. arange_dh='PAPER', col='black')
  771. lim = plt.gca().get_ylim()
  772. plt.ylim(lim[0],(lim[1]-lim[0])*1.01+lim[0])
  773. DrawLines(a,b,c,d,m_mpr, col='darkslategray')
  774. DrawLines(a_,b_,c_,d_,m_fl, col='darkslategray')
  775. ticklabels = ['T1_MPR\nSHAKE', 'T1_MPR\nNOD', 'T2_FLAIR\nNOD']
  776. ticks = [2.5, 6.5, 11.5]
  777. plt.xticks(labels=ticklabels, ticks=ticks, fontsize=12)
  778. ax.text(-0.18, 0.9, string.ascii_lowercase[num], transform=ax.transAxes,
  779. size=21, weight='bold')
  780. num += 1
  781. plt.yticks(fontsize=13)
  782. plt.tick_params('both', length=0)
  783. plt.gca().yaxis.set_major_formatter(ScalarFormatter(useOffset=False))
  784. if i == 3:
  785. plt.yticks(ticks=[1, 2, 3, 4, 5])
  787. legend.get_frame().set_linewidth(2)
  788. plt.subplots_adjust(hspace=0.2, wspace=0.3)
  789. plt.savefig(out_dir+'Metrics_3D'+save+'.tiff', format='tiff',
  790. bbox_inches='tight', dpi=200)
  791. plt.savefig(out_dir+'Metrics_3D'+save+'.png', format='png',
  792. bbox_inches='tight', dpi=200)
  793. plt.show()
  797. # now TSE, STIR, T2* and TRACEW
  798. metrics = [ssims, psnrs, tgs, qss, np.array(qss[-1])]
  799. p_values = [p_ssims, p_psnrs, p_tgs, p_qss[0:4], np.array([p_qss[-1]])]
  800. labels = ['SSIM', 'PSNR', 'Tenengrad', 'Observer Scores', 'Quality Rank']
  801. color_dict = {'prospective MoCo no REAC': 'tab:blue',
  802. 'no MoCo no REAC':'tab:orange',
  803. 'prospective MoCo with REAC': 'tab:green',
  804. 'no MoCo with REAC':'tab:cyan'}
  805. name_dict = {'prospective MoCo no REAC': 'with PMC without REAC',
  806. 'no MoCo no REAC':'without PMC without REAC',
  807. 'prospective MoCo with REAC': 'with PMC with REAC',
  808. 'no MoCo with REAC':'without PMC with REAC'}
  809. small = dict(markersize=3)
  811. num = 0
  812. plt.figure(figsize=(13,8))
  813. for i, metric, p, lab in zip(range(0,5), metrics, p_values, labels):
  814. means = []
  815. names = []
  816. m_nod = []
  818. if i == 4:
  819. tmp = metric
  820. m_nod= [tmp[:,2:4]]
  821. means = [np.nanmean(tmp[:,2:4], axis=0)]
  822. names = [names_tes[j][2:4]]
  824. p_mot = [p[0,1]]
  825. indices = [ind_ps[-1][1]]
  827. elif i in [0,1,2,3]:
  828. for j in range(2, 6):
  829. end = 6
  830. if j in [4,5]:
  831. end = 4
  832. tmp = metric[j]
  833. m_nod.append(tmp[:,2:end])
  834. means.append(np.nanmean(tmp[:,2:end], axis=0))
  835. names.append(names_tes[j][2:end])
  837. p_mot = []
  838. indices = []
  839. for pval, ind in zip(p[2:6], ind_ps[2:6]):
  840. p_mot.append(pval[1:])
  841. indices.append(ind[1:])
  843. if i in [0,1,2]:
  844. ax = plt.subplot2grid((2,4), (i//2,i%2*2), colspan=2)
  845. if i == 3:
  846. ax = plt.subplot2grid((2,6), (1,3), colspan=2)
  847. if i == 4:
  848. ax = plt.subplot2grid((2,6), (1,5), colspan=1)
  850. N = 1
  851. lims = []
  853. for m, mean, p, index, name in zip(m_nod, means, p_mot, indices, names):
  854. colors = []
  855. leg_labels = []
  856. for n in name:
  857. colors.append(color_dict[n])
  858. leg_labels.append(name_dict[n])
  859. if i in [0,1,2]:
  860. positions = np.arange(N, N+len(m.T))
  861. box1 = plt.boxplot(m, positions=positions, flierprops=small,
  862. widths=0.5)
  863. for j in range(len(mean)):
  864. plt.errorbar(positions[j], mean[j], yerr=None, color=colors[j],
  865. fmt='.', capsize=3)
  866. if i == 2 and N == 1:
  867. for j in range(0,4):
  868. plt.errorbar(positions[j], mean[j], yerr=None,
  869. color=colors[j], fmt='.',
  870. capsize=3, label=leg_labels[j])
  871. legend = plt.legend( loc='lower left', ncol=2,
  872. bbox_to_anchor=(0.45, -0.4), fontsize=12,
  873. frameon=True)
  874. ticklabels = ['T2_TSE', 'T1_STIR', 'T2*', 'TRACEW']
  875. ticks = [2.5, 7.5, 11.5, 14.5]
  876. plt.xticks(labels=ticklabels, ticks=ticks, fontsize=12)
  878. if i == 3:
  879. positions = np.arange(N, N+len(m.T))
  880. box1 = plt.boxplot(m[0:8], positions=positions[0:8], flierprops=small,
  881. widths=0.5)
  882. for j in range(len(mean)):
  883. plt.errorbar(positions[j], mean[j], yerr=None, color=colors[j],
  884. fmt='.', capsize=3)
  885. ticklabels = ['T2_TSE', 'T1_STIR']
  886. ticks = [2.5, 7.5]
  887. plt.xticks(labels=ticklabels, ticks=ticks, fontsize=12)
  889. if i == 4:
  890. positions = np.arange(N, N+len(m.T))
  891. box1 = plt.boxplot(m, positions=positions, flierprops=small,
  892. widths=0.5)
  893. for j in range(len(mean)):
  894. plt.errorbar(positions[j], mean[j], yerr=None, color=colors[j],
  895. fmt='.', capsize=3)
  896. ticklabels = ['DWI']
  897. ticks = [1.5]
  898. plt.xticks(labels=ticklabels, ticks=ticks, fontsize=12)
  901. for patch, patch2, color in zip(box1['boxes'], box1['medians'], colors):
  902. patch.set(color=color, lw=1.5)
  903. patch2.set(color='k', lw=1.5)
  906. if show_stat_test == True:
  907. maxi = []
  908. for v in m.T:
  909. maxi.append(np.amax(v))
  910. if i in [0,1,2]:
  911. Show_Stars(np.array(p), index-2, positions, maxi,
  912. flexible_dh=True, col='black')
  913. if i == 3:
  914. Show_Stars(np.array(p), index-2, positions[0:8], maxi[0:8],
  915. flexible_dh=True, col='black')
  917. if i == 4:
  918. Show_Stars(np.array([p]), [index-2], positions, maxi,
  919. flexible_dh=True, col='black')
  921. lims.append(plt.gca().get_ylim())
  923. a = [0,2]
  924. b = [1,3]
  925. if N > 10 or i == 4:
  926. a = [0]
  927. b = [1]
  929. c = [positions[p] for p in a]
  930. d = [positions[p] for p in b]
  931. DrawLines(a,b,c,d,m, col='darkslategray')
  932. N = N+len(m.T)+1
  934. if show_stat_test:
  935. lims = np.array(lims)
  936. lim_min = np.amin(lims[:,0])
  937. lim_max = np.amax(lims[:,1])
  938. plt.ylim(lim_min,(lim_max-lim_min)*1.01+lim_min)
  940. plt.ylabel(lab, fontsize=15)
  941. plt.yticks(fontsize=13)
  942. plt.tick_params('both', length=0)
  943. plt.gca().yaxis.set_major_formatter(ScalarFormatter(useOffset=False))
  944. if i == 4:
  945. ax.text(-0.51, 0.9, string.ascii_lowercase[num], transform=ax.transAxes,
  946. size=21, weight='bold')
  947. else:
  948. ax.text(-0.18, 0.9, string.ascii_lowercase[num], transform=ax.transAxes,
  949. size=21, weight='bold')
  950. num += 1
  951. if i == 3:
  952. plt.yticks(ticks=[1, 2, 3, 4, 5])
  953. if i == 4:
  954. plt.yticks(ticks=[1, 2, 3, 4])
  956. legend.get_frame().set_linewidth(2)
  957. plt.subplots_adjust(hspace=0.2, wspace=0.6)
  958. plt.savefig(out_dir+'Metrics_2D'+save+'.tiff', format='tiff',
  959. bbox_inches='tight', dpi=200)
  960. plt.savefig(out_dir+'Metrics_2D'+save+'.png', format='png',
  961. bbox_inches='tight', dpi=200)
  962. plt.show()
  966. ''' (5) Plot DWI analysis: '''
  967. if plot_DWI:
  968. # calculate histograms of ACS difference images:
  969. sub_out_dir = '../Results/Metrics_Results/'
  970. im_dir = '../BIDSdata_defaced/'
  971. in_dir = '../RegistrationTransforms/'
  972. bm_dir = '../Brainmasks/'
  975. subdir = []
  976. for i in range(1,10):
  977. subdir.append('Subject_0'+str(i)+'/')
  978. for i in range(10,20):
  979. subdir.append('Subject_'+str(i)+'/')
  980. for i in range(20,23):
  981. subdir.append('Subject_'+str(i)+'/')
  984. if calc_ADC_hist:
  985. all_mean, all_std, all_sum = [], [], []
  986. for sub in subdir:
  987. #print('Subject of interest is '+sub)
  988. # sort out subjects for which no DWI scans available:
  989. tmp = os.listdir(sub_out_dir+sub)
  990. tmp_ = ''
  991. if 'ADC' not in tmp_.join(tmp):
  992. continue
  994. differences = []
  995. descr = []
  996. still_off = glob.glob(im_dir+sub+'/TCLMOCO_OFF_STILL_EP2D_DIFF_EXTTRACKING_ADC*.nii')[0]
  997. still = nib.load(still_off).get_fdata().astype(np.uint16)
  998. bm = glob.glob(bm_dir+sub+'bm_mov_*ADC*.nii')[0]
  999. #glob.glob(in_dir+sub+'bm_mov_*ADC*.nii')[0]
  1000. bm = nib.load(bm).get_fdata().astype(np.uint16)
  1002. still = still*bm
  1004. all_imgs = [x for x in glob.glob(im_dir+sub+'/*ADC*.nii') if x not in still_off]
  1005. for im in all_imgs:
  1006. img = nib.load(im).get_fdata().astype(np.uint16)
  1007. img = img*bm
  1009. diff = still.astype(np.float) - img.astype(np.float)
  1010. differences.append(diff)
  1011. descr.append(os.path.basename(im))
  1012. differences, descr = np.array(differences), np.array(descr)
  1013. ind = np.argsort(descr)
  1014. differences = differences[ind]
  1015. descr = descr[ind]
  1016. sums, means, stds = [], [], []
  1017. for i in range(0, len(descr)):
  1018. diff = differences[i][bm!=0] # only look at voxels inside the brain!
  1019. n, bins, tmp = plt.hist(diff, bins = 2000)
  1020. plt.title(descr[i][12:25])
  1021. sums.append(np.sum(np.abs(diff))/len(diff))
  1023. mids = 0.5*(bins[1:]+bins[:-1])
  1024. mean = np.average(mids, weights=n)
  1025. std = np.sqrt(np.average((mids-mean)**2, weights=n))
  1026. means.append(mean)
  1027. stds.append(std)
  1029. all_mean.append(means)
  1030. all_std.append(stds)
  1031. all_sum.append(sums)
  1033. # look at results for all subjects:
  1034. All_mean, All_std, All_sum = np.array(all_mean)[:,::-1], np.array(all_std)[:,::-1], np.array(all_sum)[:,::-1]
  1035. save_arr = np.array([All_mean, All_std, All_sum])
  1036. np.save(out_dir_metrics+'ADC_all_values_'+save, save_arr)
  1039. # load the values for the ADC histograms and plot them:
  1040. file = glob.glob(out_dir_metrics+'ADC_all_values_**.npy')
  1041. if len(file)>0:
  1042. # get the most recent file:
  1043. #file = SortFiles(file)
  1044. print(file[0])
  1046. All_mean, All_std, All_sum = np.load(file[0])
  1047. All_mean = np.array([All_mean[:,0], All_mean[:,2], All_mean[:,1]]).T
  1048. All_std = np.array([All_std[:,0], All_std[:,2], All_std[:,1]]).T
  1049. All_sum = np.array([All_sum[:,0], All_sum[:,2], All_sum[:,1]]).T
  1050. names_diff = ['prospective MoCo', 'no MoCo', 'prospective MoCo']
  1052. p_mean, rej_mean, ind, altern = PerformWilcoxonAllImg('Mean', All_mean, 'ADC',
  1053. out_dir_metrics,
  1054. save, option='diff')
  1055. p_std, rej_std, ind, altern = PerformWilcoxonAllImg('Std', All_std, 'ADC',
  1056. out_dir_metrics,
  1057. save, option='diff')
  1058. p_sum, rej_sum, ind, altern = PerformWilcoxonAllImg('Sum', All_sum, 'ADC',
  1059. out_dir_metrics,
  1060. save, option='diff')
  1062. mean_mean = np.mean(All_mean, axis=0)
  1063. mean_std = np.mean(All_std, axis=0)
  1064. mean_sum = np.mean(All_sum, axis=0)
  1066. color_dict = { 'prospective MoCo': 'tab:blue', 'no MoCo': 'tab:orange'}
  1067. names = ['with PMC', 'without PMC', 'with PMC']
  1069. colors_te = []
  1070. for n in names_diff:
  1071. colors_te.append(color_dict[n])
  1072. x = np.arange(1,len(mean_mean)+1)
  1075. plt.figure(figsize=(10,3))
  1076. ax=plt.subplot(1,2,1)
  1077. for i in range(len(x)):
  1078. plt.errorbar(x[i], mean_mean[i], yerr=None, color=colors_te[i], fmt='.',
  1079. capsize=3)
  1081. small = dict(markersize=3)
  1082. box1 = plt.boxplot(All_mean, flierprops=small)
  1083. for patch, patch2, color in zip(box1['boxes'], box1['medians'], colors_te):
  1084. patch.set(color=color, lw=1.7)
  1085. patch2.set(color='k', lw=1.7)
  1087. plt.xticks(labels=[], ticks=[])
  1088. plt.gca().yaxis.set_major_formatter(ScalarFormatter(useOffset=False))
  1089. for y1, y2 in zip(All_mean[:,1], All_mean[:,2]):
  1090. plt.plot([2,3], [y1, y2], 'darkslategray', lw=1)
  1091. plt.annotate('Still', xy=(0.17, -0.1), xytext=(0.17, -0.18),
  1092. xycoords='axes fraction', fontsize=12, ha='center', va='bottom',
  1093. bbox=dict(boxstyle='square', fc='white', ec='grey', lw=1.5),
  1094. arrowprops=dict(arrowstyle='-[, widthB=1.4, lengthB=0.7',
  1095. lw=1.5, color='grey'))
  1096. plt.annotate('Nod', xy=(0.67, -0.1), xytext=(0.67, -0.18),
  1097. xycoords='axes fraction', fontsize=12, ha='center', va='bottom',
  1098. bbox=dict(boxstyle='square', fc='white', ec='grey', lw=1.5),
  1099. arrowprops=dict(arrowstyle='-[, widthB=3.0, lengthB=0.7',
  1100. lw=1.5, color='grey'))
  1101. Show_Stars(p_mean, ind, x, np.amax(All_mean, axis=0), arange_dh='diff',
  1102. col='black')
  1103. lim = plt.gca().get_ylim()
  1104. plt.ylim(lim[0],(lim[1]-lim[0])*1.1+lim[0])
  1105. plt.ylabel('Mean of histogram [$\\frac{\mu m^2}{s}$]')
  1106. ax.text(-0.22, 0.9, string.ascii_lowercase[0], transform=ax.transAxes,
  1107. size=21, weight='bold')
  1109. ax=plt.subplot(1,2,2)
  1110. for i in range(1):
  1111. plt.errorbar(x[i], mean_std[i], yerr=None, color=colors_te[i], fmt='.',
  1112. capsize=3)
  1113. for i in range(1,len(x)):
  1114. plt.errorbar(x[i], mean_std[i], yerr=None, label=names[i],
  1115. color=colors_te[i], fmt='.', capsize=3)
  1117. small = dict(markersize=3)
  1118. box1 = plt.boxplot(All_std, flierprops=small)
  1119. for patch, patch2, color in zip(box1['boxes'], box1['medians'], colors_te):
  1120. patch.set(color=color, lw=1.7)
  1121. patch2.set(color='k', lw=1.7)
  1123. plt.xticks(labels=[], ticks=[])
  1124. plt.gca().yaxis.set_major_formatter(ScalarFormatter(useOffset=False))
  1125. for y1, y2 in zip(All_std[:,1], All_std[:,2]):
  1126. plt.plot([2,3], [y1, y2], 'darkslategray', lw=1)
  1127. plt.annotate('Still', xy=(0.17, -0.1), xytext=(0.17, -0.18),
  1128. xycoords='axes fraction', fontsize=12, ha='center', va='bottom',
  1129. bbox=dict(boxstyle='square', fc='white', ec='grey', lw=1.5),
  1130. arrowprops=dict(arrowstyle='-[, widthB=1.4, lengthB=0.7',
  1131. lw=1.5, color='grey'))
  1132. plt.annotate('Nod', xy=(0.67, -0.1), xytext=(0.67, -0.18),
  1133. xycoords='axes fraction', fontsize=12, ha='center', va='bottom',
  1134. bbox=dict(boxstyle='square', fc='white', ec='grey', lw=1.5),
  1135. arrowprops=dict(arrowstyle='-[, widthB=3.0, lengthB=0.7',
  1136. lw=1.5, color='grey'))
  1137. Show_Stars(p_std, ind, x, np.amax(All_std, axis=0), arange_dh='diff',
  1138. col='black')
  1139. lim = plt.gca().get_ylim()
  1140. plt.ylim(lim[0],(lim[1]-lim[0])*1.1+lim[0])
  1141. plt.ylabel('Std of histogram [$\\frac{\mu m^2}{s}$]')
  1142. legend = plt.legend(loc='upper center', ncol = 2,
  1143. bbox_to_anchor=(-0.2, -0.25), frameon=True)
  1144. ax.text(-0.22, 0.9, string.ascii_lowercase[1], transform=ax.transAxes,
  1145. size=21, weight='bold')
  1147. legend.get_frame().set_linewidth(2)
  1148. plt.subplots_adjust(hspace=0.2, wspace=0.3)
  1149. plt.savefig(out_dir+'ADC'+save+'.tiff', format='tiff', bbox_inches='tight',
  1150. dpi=200)
  1151. plt.savefig(out_dir+'ADC'+save+'.png', format='png', bbox_inches='tight',
  1152. dpi=200)
  1153. plt.show()
  1156. # calculate mean value of ground truth scans:
  1157. means = []
  1158. for sub in subdir:
  1159. if sub in ['HC_0'+str(i)+'/' for i in range(1,10)]:
  1160. continue
  1161. if sub in ['HC_'+str(i)+'/' for i in range(10,13)]:
  1162. continue
  1164. # sort out subjects for which no DWI scans available:
  1165. tmp = os.listdir(sub_out_dir+sub)
  1166. tmp_ = ''
  1167. if 'ADC' not in tmp_.join(tmp):
  1168. continue
  1170. file = glob.glob(im_dir+sub+'/TCLMOCO_OFF_STILL_EP2D_DIFF_EXTTRACKING_ADC*.nii')[0]
  1171. #glob.glob(im_dir+sub+'/*ADC*.nii')[0]
  1172. img = nib.load(file).get_fdata().astype(np.uint16)
  1173. bm_file = glob.glob(bm_dir+sub+'bm_mov_*ADC*.nii')[0]
  1174. bm = nib.load(bm_file).get_fdata().astype(np.uint16)
  1176. bm_fl = bm.flatten()
  1177. dat = img.flatten()
  1178. img_fl = dat[bm_fl>0]
  1180. means.append(np.mean(img_fl))
  1182. print('Mean Values of ground truth scans:', means)
  1183. print('Overall mean value:', np.mean(means))

generate_plots_manuscript.py at commit ffd7086, under MIT · at the source

Overview

Authors: Kathrine Skak Madsen1, Tim Ruschke2,3, Hannah Eichhorn4,5, Puk Rising2, Melanie Ganz2,3
ORCID iDs: Melanie Ganz
  1. Danish Research Centre for Magnetic Resonance, Department of Radiology and Nuclear Medicine, Copenhagen University Hospital - Amager and Hvidovre,Hvidovre, Denmark
  2. Neurobiology Research Unit, Copenhagen University Hospital - Rigshospitalet,Copenhagen, Denmark
  3. Department of Computer Science, University of Copenhagen,Copenhagen, Denmark
  4. Institute of Machine Learning in Biomedical Imaging, Helmholtz Munich, Neuherberg, Germany
  5. School of Computation, Information and Technology, Technical University of Munich,Munich, Germany
Institutions: Copenhagen University Hospital (Denmark); University of Copenhagen (Denmark); Technical University of Munich (Germany)
Journal: Scientific data, volume 13, issue 1, article 943
Dates: received 16 July 2025; accepted 26 March 2026; published online 23 April 2026
Type: Data paper · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41597-026-07144-z · PMID 42026101 · PMCID PMC13315939 · OpenAlex W7155394641
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism), methods / tools (subfield)
Methods: fMRI & imaging
Keywords: Neurology, Mathematics and computing, Scientific data
MeSH: Brain*, Magnetic Resonance Imaging*, Motion*, Artifacts, Humans, Image Processing, Computer-Assisted (* major topic)
Journal subjects: Data Descriptor
Topic: Advanced MRI Techniques and Applications (Radiology, Nuclear Medicine and Imaging, Medicine), according to OpenAlex
Funding: Elsass Fonden (18-3-0147)
Citations: not cited yet (Europe PMC); 34 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

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

melanieganz/mocoproject

License: MIT
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: ffd708616a24b5abe0531605e0e7a30bdb2df6cc, 10 May 2026
Languages: Python (28), Jupyter (2), Shell (1)
Size: 53 files, 31 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (docker/dockerfile), 2 notebooks
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (30 files), SciPy (14 files), Matplotlib (11 files), NiBabel (9 files), scikit-image (8 files), FreeSurfer (7 files), pydicom (4 files), statsmodels (4 files)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
33 files
At the source: github.com/melanieganz/MoCoProject/tree/main/MotionCorrectedClinicalMRProtocol

OSF vzh4g

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: the link answers
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
At the source: osf.io/vzh4g

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41597-026-07144-z.

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:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 31 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 paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • no repository, dataset or request procedure was recognized in it

Read it in the paper: doi.org/10.1038/s41597-026-07144-z.

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, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 3 keywords, 6 MeSH terms, 1 funder, 31 references.

Cite

This paper

Skak Madsen, K., Ruschke, T., Eichhorn, H., Rising, P., & Ganz, M. (2026). Brain MRI Dataset Featuring a Full Clinical Protocol With and Without Intentional Motion. Scientific data, 13(1), 943. https://doi.org/10.1038/s41597-026-07144-z

BibTeX

@article{skakmadsen2026brain,
author = {Skak Madsen, Kathrine and Ruschke, Tim and Eichhorn, Hannah and Rising, Puk and Ganz, Melanie},
title = {{Brain MRI Dataset Featuring a Full Clinical Protocol With and Without Intentional Motion}},
journal = {Scientific data},
year = {2026},
month = apr,
volume = {13},
number = {1},
pages = {943},
publisher = {Nature Publishing Group},
issn = {2052-4463},
doi = {10.1038/s41597-026-07144-z},
url = {https://doi.org/10.1038/s41597-026-07144-z},
pmid = {42026101},
pmcid = {PMC13315939}
}

RIS

TY - JOUR
AU - Skak Madsen, Kathrine
AU - Ruschke, Tim
AU - Eichhorn, Hannah
AU - Rising, Puk
AU - Ganz, Melanie
TI - Brain MRI Dataset Featuring a Full Clinical Protocol With and Without Intentional Motion
T2 - Scientific data
J2 - Sci Data
PY - 2026
DA - 2026/04/23
VL - 13
IS - 1
SP - 943
SN - 2052-4463
PB - Nature Publishing Group
DO - 10.1038/s41597-026-07144-z
UR - https://doi.org/10.1038/s41597-026-07144-z
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41597-026-07144-z",
"type": "article-journal",
"title": "Brain MRI Dataset Featuring a Full Clinical Protocol With and Without Intentional Motion",
"container-title": "Scientific data",
"author": [
{
"family": "Skak Madsen",
"given": "Kathrine"
},
{
"family": "Ruschke",
"given": "Tim"
},
{
"family": "Eichhorn",
"given": "Hannah"
},
{
"family": "Rising",
"given": "Puk"
},
{
"family": "Ganz",
"given": "Melanie"
}
],
"container-title-short": "Sci Data",
"volume": "13",
"issue": "1",
"page": "943",
"DOI": "10.1038/s41597-026-07144-z",
"PMID": "42026101",
"PMCID": "PMC13315939",
"ISSN": "2052-4463",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41597-026-07144-z",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
23
]
]
}
}

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.70368 [code]
Intra-MRI Head Motion Tracking and Correction: A Quantitative In Vivo Evaluation Framework.
Journal: NMR in biomedicine
In common: scikit-image, NiBabel, SciPy, 2 other tools, methods / tools, structural MRI / diffusion, 7 references
[2] doi:10.1002/nbm.70286 [code]
DEEP-DISORDER: Motion Correction in 3D MRI via Segment Reconstruction and Registration.
Journal: NMR in biomedicine
In common: structural MRI / diffusion, 8 references
[3] doi:10.1038/s41597-026-07377-y [code]
An open-access multi-site fMRI dataset for investigating conscious visual perception.
Journal: Scientific data
In common: FreeSurfer, scikit-image, NiBabel, 4 other tools, 2 references
[4] doi:10.1162/imag.a.1366 [code]
MICAFlow: Fast and robust MRI preprocessing bridging research neuroimaging and clinical practice.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: FreeSurfer, NiBabel, SciPy, 2 other tools, methods / tools, structural MRI / diffusion, 3 references
[5] doi:10.1162/imag.a.1262 [code]
Frame-wise multi-echo distortion correction for superior functional MRI.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: pydicom, FreeSurfer, scikit-image, 4 other tools, 1 reference
[6] doi:10.1002/alz.71649 [code]
Postmortem brain MRI reveals differential associations of subcortical and limbic volumes with cortical thinning and neurodegenerative pathologies.
Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association
In common: FreeSurfer, scikit-image, NiBabel, 4 other tools, structural MRI / diffusion, 1 reference
[7] doi:10.7554/elife.108408 [code]
Frequency and laminar profile of feature-specific visual activity revealed by interleaved EEG-fMRI.
Journal: eLife
In common: FreeSurfer, scikit-image, NiBabel, 4 other tools, 1 reference
[8] doi:10.7554/elife.107933 [code]
Modality-agnostic decoding of vision and language from fMRI.
Journal: eLife
In common: FreeSurfer, scikit-image, NiBabel, 4 other tools, 1 reference
[9] doi:10.1038/s41597-026-07350-9 [code]
An open multi-center MEG-EEG dataset for studying conscious visual perception.
Journal: Scientific data
In common: FreeSurfer, scikit-image, NiBabel, 4 other tools, structural MRI / diffusion
[10] doi:10.1038/s41597-026-06869-1 [code]
Individual Brain Charting: fifth release of high-resolution fMRI data for cognitive mapping.
Journal: Scientific data
In common: FreeSurfer, scikit-image, NiBabel, 3 other tools, methods / 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.