OSCR

Automated Surface-Based Segmentation of Deep Gray Matter Regions Based on Diffusion Tensor Images Reveals Unique Age Trajectories Over the Healthy Lifespan.

Code ↔ Paper

6 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 6 matches
  1. [1] § Materials and Methods › Surface‐Based Deep GM Segmentation on DTI ↔ microbrain/subcort_segmentation/mbrain_segment.py, lines 898–962 · score 0.81 · Harvard Oxford, initial mesh, Initial voxel, probability map, DTI imaging, atlas
  2. [2] § Materials and Methods › DTI Preprocessing ↔ scripts/microbrain_run.py, lines 275–344 · score 0.79 · Gibbs ringing, tensor models, brain mask, MD map, BET, denoised
  3. [3] § Materials and Methods › Surface‐Based Deep GM Segmentation on DTI ↔ microbrain/surfing/mbrain_cortical_segmentation.py, lines 850–928 · score 0.72 · cortical segmentation, CSF probabilities, probability map, external, forces, internal
  4. [4] § Materials and Methods › Surface‐Based Deep GM Segmentation on DTI ↔ microbrain/subcort_segmentation/mbrain_segment.py, lines 898–962 · score 0.57 · Harvard Oxford, probability maps, deformed, deformation, native, mesh
  5. [5] § Materials and Methods › Surface‐Based Deep GM Segmentation on DTI ↔ microbrain/subcort_segmentation/mbrain_segment.py, lines 581–651 · score 0.55 · Initial segmentation, right thalamus, closest, intensity, force, atlas
  6. [6] § Materials and Methods › DTI Preprocessing ↔ scripts/microbrain_run.py, lines 275–344 · score 0.54 · tensor models, MD map, gradients, fit, b0, Preprocessing

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 · 962 lines · 48 KB · MIT · 3 matches

  1. #!/usr/local/bin/python
  2. import microbrain.utils.surf_util as sutil
  4. import vtk
  5. from shutil import which
  7. from os import system, environ, path
  8. import numpy as np
  9. import nibabel as nib
  10. import inspect
  12. from skimage.morphology import binary_erosion
  14. # FSL Harvard Atlas label indices
  15. LEFT_WHITE_IDX = 0
  16. LEFT_CORTEX_IDX = 1
  17. LEFT_VENT_IDX = 2
  18. LEFT_THAL_IDX = 3
  19. LEFT_CAUDATE_IDX = 4
  20. LEFT_PUT_IDX = 5
  21. LEFT_GLOB_IDX = 6
  22. BRAIN_STEM_IDX = 7
  23. LEFT_HIPPO_IDX = 8
  24. LEFT_AMYG_IDX = 9
  25. LEFT_ACCUM_IDX = 10
  27. RIGHT_WHITE_IDX = 11
  28. RIGHT_CORTEX_IDX = 12
  29. RIGHT_VENT_IDX = 13
  30. RIGHT_THAL_IDX = 14
  31. RIGHT_CAUDATE_IDX = 15
  32. RIGHT_PUT_IDX = 16
  33. RIGHT_GLOB_IDX = 17
  34. RIGHT_HIPPO_IDX = 18
  35. RIGHT_AMYG_IDX = 19
  36. RIGHT_ACCUM_IDX = 20
  38. # indices specific for microbrain ROIs
  39. LEFT_HIPPOAMYG_IDX = 88
  40. LEFT_STRIATUM_IDX = 66
  41. RIGHT_HIPPOAMYG_IDX = 188
  42. RIGHT_STRIATUM_IDX = 166
  43. VENT_IDX = 200
  46. def is_tool(name):
  47. """
  48. Checks to see if third party binary is installed and accessible by command line
  50. Parameters
  51. ----------
  52. name: string
  53. the name of the command to look for
  55. Returns
  56. -------
  57. Boolean: true if command exists, false if not in path
  58. """
  60. return which(name) is not None
  63. def fsl_ext():
  64. """
  65. Gets the preferred image extension used by FSL
  67. Parameters
  68. ----------
  69. none
  71. Returns
  72. -------
  73. extension: a string containing the preferred output
  74. """
  76. fsl_extension = ''
  77. if environ['FSLOUTPUTTYPE'] == 'NIFTI':
  78. fsl_extension = '.nii'
  79. elif environ['FSLOUTPUTTYPE'] == 'NIFTI_GZ':
  80. fsl_extension = '.nii.gz'
  81. return fsl_extension
  84. def get_fsl_standard_dir():
  85. """
  86. Return tissue directory in microbrain repository
  88. Returns
  89. -------
  90. tissue_dir: string
  91. tissue path
  92. """
  94. import microbrain # ToDo. Is this the only way?
  95. module_path = inspect.getfile(microbrain)
  97. fsl_standard_dir = path.dirname(module_path) + "/data/fsl/standard/"
  99. return fsl_standard_dir
  102. def get_fsl_atlas_dir():
  103. """
  104. Return tissue directory in microbrain repository
  106. Returns
  107. -------
  108. tissue_dir: string
  109. tissue path
  110. """
  112. import microbrain # ToDo. Is this the only way?
  113. module_path = inspect.getfile(microbrain)
  115. fsl_atlas_dir = path.dirname(module_path) + "/data/fsl/atlases/"
  117. return fsl_atlas_dir
  119. def register_probatlas_to_native(fsource, ftemplate, fatlas, regDir, cpu_num=0):
  120. """
  121. Given an image will register this image to the template (ANTs SyN)
  122. Then applies this transform to an atlas
  124. Parameters
  125. ----------
  126. fsource: string
  127. filename of image in native space
  128. ftemplate: string
  129. filename of image template in MNI space (or normal space)
  130. fatlas: string
  131. filename of 4D nifti file in MNI space to be registered
  132. regDir: string
  133. directory to store output
  135. Optional Parameters
  136. -------------------
  137. cpu_num: integer
  138. number of cpu threads to use for registration
  140. Returns
  141. -------
  142. fatlas_out: filename for atlas registered to native space
  143. """
  144. if not path.exists(regDir):
  145. system('mkdir ' + regDir)
  147. if cpu_num > 0:
  148. cpu_str = ' -n ' + str(cpu_num)
  149. else:
  150. cpu_str = ''
  152. method_suffix = '_ANTsReg_native'
  153. ftempbasename = path.basename(ftemplate)
  155. fatlas_basename = path.basename(fatlas)
  156. fatlas_out = regDir + \
  157. fatlas_basename.replace('.nii.gz', method_suffix + fsl_ext())
  159. ftransform_out = regDir + 'ants_native2mni_'
  161. print('Running Ants Registration')
  162. # Make template mask
  163. ftemplate_mask = regDir + \
  164. ftempbasename.replace('.nii.gz', '_mask' + fsl_ext())
  165. if not path.exists(ftemplate_mask):
  166. system('fslmaths ' + ftemplate + ' -bin ' + ftemplate_mask)
  168. if not path.exists(ftransform_out + '1InverseWarp.nii.gz'):
  169. system('antsRegistrationSyN.sh' +
  170. ' -d 3' +
  171. ' -f ' + ftemplate +
  172. ' -m ' + fsource +
  173. ' -o ' + ftransform_out +
  174. ' -x ' + ftemplate_mask +
  175. cpu_str)
  177. print('Warping probabilistic atlas to native space')
  178. if not path.exists(fatlas_out):
  179. system('antsApplyTransforms' +
  180. ' -t [' + ftransform_out + '0GenericAffine.mat,1] ' +
  181. ' -t ' + ftransform_out + '1InverseWarp.nii.gz ' +
  182. ' -r ' + fsource +
  183. ' -e 3 ' +
  184. ' -i ' + fatlas +
  185. ' -o ' + fatlas_out)
  186. else:
  187. print("ANTs nonlinear registeration of atlases already performed")
  189. return fatlas_out
  192. def initial_voxel_labels_from_harvard(fharvard, output_prefix, outDir, md_file=None, md_thresh=0.0015, fa_file=None, fa_thresh=0.5):
  193. """
  194. Generates 3D meshes from the FSL harvard brain atlas.
  196. Parameters
  197. ----------
  198. fharvard: string
  199. filename for the harvard probabilistic atlas to be used (usually registered to native space)
  200. output_prefix: string
  201. prefix used for filenames for microbrain subject
  202. outDir: string
  203. output directory to store all files
  205. Optional Parameters:
  206. --------------------
  207. md_file: string
  208. filename for mean diffusivity map used to remove voxels from labels (useful in case of poor registration to native space)
  209. md_thresh: float
  210. removes voxels in lables with MD > md_thresh
  212. fa_file: string
  213. filename for fractional anisotropy map used to remove voxels from labels (useful in case of poor registration to native space)
  214. fa_thresh: string
  215. removes voxels in lables with FA > fa_thresh
  217. Returns
  218. -------
  219. None
  220. """
  222. if md_file:
  223. md_data = nib.load(md_file).get_fdata()
  225. if fa_file:
  226. fa_data = nib.load(fa_file).get_fdata()
  228. finitlabels_prefix = outDir + output_prefix + '_initialization'
  230. harvard_img = nib.load(fharvard)
  231. harvard_data = harvard_img.get_fdata()
  233. subcort_ind = [LEFT_THAL_IDX,
  234. LEFT_CAUDATE_IDX,
  235. LEFT_PUT_IDX,
  236. LEFT_GLOB_IDX,
  237. LEFT_HIPPO_IDX,
  238. LEFT_AMYG_IDX,
  239. LEFT_ACCUM_IDX,
  240. LEFT_HIPPOAMYG_IDX,
  241. LEFT_STRIATUM_IDX,
  242. RIGHT_THAL_IDX,
  243. RIGHT_CAUDATE_IDX,
  244. RIGHT_PUT_IDX,
  245. RIGHT_GLOB_IDX,
  246. RIGHT_HIPPO_IDX,
  247. RIGHT_AMYG_IDX,
  248. RIGHT_ACCUM_IDX,
  249. RIGHT_HIPPOAMYG_IDX,
  250. RIGHT_STRIATUM_IDX,
  251. VENT_IDX]
  253. subcort_label = ['LEFT_THALAMUS',
  254. 'LEFT_CAUDATE',
  255. 'LEFT_PUTAMEN',
  256. 'LEFT_GLOBUS',
  257. 'LEFT_HIPPO',
  258. 'LEFT_AMYGDALA',
  259. 'LEFT_ACCUMBENS',
  260. 'LEFT_HIPPOAMYG',
  261. 'LEFT_STRIATUM',
  262. 'RIGHT_THALAMUS',
  263. 'RIGHT_CAUDATE',
  264. 'RIGHT_PUTAMEN',
  265. 'RIGHT_GLOBUS',
  266. 'RIGHT_HIPPO',
  267. 'RIGHT_AMYGDALA',
  268. 'RIGHT_ACCUMBENS',
  269. 'RIGHT_HIPPOAMYG',
  270. 'RIGHT_STRIATUM',
  271. 'VENTRICLES']
  273. for sind, slabel in zip(subcort_ind, subcort_label):
  274. tmplabel = np.zeros(harvard_data.shape[0:3])
  275. if sind == RIGHT_ACCUM_IDX or sind == LEFT_ACCUM_IDX:
  276. tmplabel[binary_erosion(harvard_data[:, :, :, sind] > 25)] = 1
  277. elif sind == RIGHT_HIPPOAMYG_IDX:
  278. tmplabel[binary_erosion(np.logical_or(
  279. harvard_data[:, :, :, RIGHT_HIPPO_IDX] > 35, harvard_data[:, :, :, RIGHT_AMYG_IDX] > 35))] = 1
  280. elif sind == LEFT_HIPPOAMYG_IDX:
  281. tmplabel[binary_erosion(np.logical_or(
  282. harvard_data[:, :, :, LEFT_HIPPO_IDX] > 35, harvard_data[:, :, :, LEFT_AMYG_IDX] > 35))] = 1
  283. elif sind == RIGHT_STRIATUM_IDX:
  284. tmplabel[binary_erosion(np.logical_or(np.logical_or(harvard_data[:, :, :, RIGHT_PUT_IDX] > 15,
  285. harvard_data[:, :, :, RIGHT_CAUDATE_IDX] > 15), harvard_data[:, :, :, RIGHT_ACCUM_IDX] > 15))] = 1
  286. elif sind == LEFT_STRIATUM_IDX:
  287. tmplabel[binary_erosion(np.logical_or(np.logical_or(harvard_data[:, :, :, LEFT_PUT_IDX] > 15,
  288. harvard_data[:, :, :, LEFT_CAUDATE_IDX] > 15), harvard_data[:, :, :, LEFT_ACCUM_IDX] > 15))] = 1
  289. elif sind == RIGHT_AMYG_IDX or sind == LEFT_AMYG_IDX:
  290. tmplabel[binary_erosion(harvard_data[:, :, :, sind] > 35)] = 1
  291. elif sind == RIGHT_HIPPO_IDX or sind == LEFT_HIPPO_IDX:
  292. tmplabel[binary_erosion(harvard_data[:, :, :, sind] > 35)] = 1
  293. elif sind == VENT_IDX:
  294. tmplabel[harvard_data[:, :, :, LEFT_VENT_IDX] > 25] = 1
  295. tmplabel[harvard_data[:, :, :, RIGHT_VENT_IDX] > 25] = 1
  296. else:
  297. tmplabel[binary_erosion(harvard_data[:, :, :, sind] > 50)] = 1
  299. # # remove large MD voxels in GM if argument provided
  300. if sind != VENT_IDX:
  301. if md_file:
  302. tmplabel[md_data > md_thresh] = 0
  303. if fa_file:
  304. tmplabel[fa_data > fa_thresh] = 0
  306. finit_tmplabel = finitlabels_prefix + '_' + slabel + fsl_ext()
  307. nib.save(nib.Nifti1Image(tmplabel, harvard_img.affine), finit_tmplabel)
  308. system('mirtk extract-connected-components ' +
  309. finit_tmplabel + ' ' + finit_tmplabel)
  311. # Extract surfaces
  312. finit_surf = finitlabels_prefix + '_' + slabel + '.vtk'
  313. system('mirtk extract-surface ' + finit_tmplabel +
  314. ' ' + finit_surf + ' -isovalue 0.5')
  315. system('mirtk extract-connected-points ' +
  316. finit_surf + ' ' + finit_surf)
  317. system('mirtk smooth-surface ' + finit_surf + ' ' +
  318. finit_surf + ' -iterations 50 -lambda 0.05')
  320. # Label Surfaces
  321. system('mirtk project-onto-surface ' + finit_surf + ' ' + finit_surf +
  322. ' -constant ' + str(sind) + ' -pointdata -name struct_label')
  324. # Combine left and right thalamus
  325. lh_thalamus = sutil.read_surf_vtk(
  326. finitlabels_prefix + '_LEFT_THALAMUS.vtk')
  327. rh_thalamus = sutil.read_surf_vtk(
  328. finitlabels_prefix + '_RIGHT_THALAMUS.vtk')
  329. appender = vtk.vtkAppendPolyData()
  330. appender.AddInputData(lh_thalamus)
  331. appender.AddInputData(rh_thalamus)
  332. appender.Update()
  333. thalamus_surf = appender.GetOutput()
  334. sutil.write_surf_vtk(thalamus_surf, finitlabels_prefix + '_THALAMUS.vtk')
  336. return
  338. def extract_surfaces_from_labels(flabels, label_list, outDir, fout_prefix):
  339. """
  340. Given a 3D nifti file containing a voxelwise labeling of structures, output a smoothed 3D surface representation of the labels
  342. Parameters
  343. ----------
  344. flabels: string
  345. filename for 3D nifti file containing labels
  346. label_list: string
  347. list of label integers indicating which labels to convert to surfaces (e.g. [1, 3])
  348. outDir: string
  349. directory to output the meshes
  350. fout_prefix: string
  351. prefix for output files
  353. Returns
  354. -------
  355. list of outputed vtk files
  356. """
  358. # Read in labels output nifti images
  359. label_img = nib.load(flabels)
  360. label_data = label_img.get_fdata()
  362. fout_list = []
  363. for label_ind in label_list:
  364. # Write tmp nifti
  365. ftmpNII = outDir + 'tmplabel' + str(label_ind) + fsl_ext()
  366. tmplabel = np.zeros(label_data.shape)
  367. tmplabel[label_data == label_ind] = 1
  368. nib.save(nib.Nifti1Image(tmplabel,
  369. label_img.affine), ftmpNII)
  371. # Extract surface around label
  372. fout_surf = fout_prefix + '_' + str(label_ind) + '.vtk'
  373. system('mirtk extract-surface ' + ftmpNII +
  374. ' ' + fout_surf + ' -isovalue 0.5')
  375. system('rm ' + ftmpNII)
  377. # Smooth surface for the label
  378. system('mirtk smooth-surface ' + fout_surf + ' ' +
  379. fout_surf + ' -iterations 50 -lambda 0.2')
  381. # Add label to surface mesh
  382. system('mirtk project-onto-surface ' + fout_surf + ' ' + fout_surf +
  383. ' -constant ' + str(label_ind) + ' -pointdata -name struct_label')
  384. fout_list = fout_list + [fout_surf]
  386. return fout_list
  389. def deform_subcortical_surfaces(fdwi, ffa, fmd, fharvard_native, segDir, initSegDir, subID, cpu_num=0):
  390. """
  391. Script to segment subcortical brain regions with 3D surface based deformation using DTI images/maps
  393. Parameters
  394. ----------
  395. fdwi: string
  396. filename for mean diffusion weigthed image (suggested b1000) used for globus pallidus segmentation
  397. ffa: string
  398. filename for fa map
  399. ffmd: string
  400. filename for md map
  401. fharvard_native: string
  402. FSL harvard probabilistic atlas transformed to native space
  403. segDir: string
  404. parent directory to store pipeline output
  405. initSegDir: string
  406. directory containing initial structure surfaces generated using initial_voxel_labels_from_harvard
  407. subID: string
  408. microbrain subject ID
  410. Optional Parameters
  411. -------------------
  412. cpu_num: integer
  413. number of cpu threads for MIRTK deform-mesh to use (default is mirtk default of all available threads)
  415. Returns
  416. -------
  417. none
  418. """
  420. min_edgelength = 0.6
  421. max_edgelength = 1.1
  423. min_hippo_edgelength = 0.6
  424. max_hippo_edgelength = 1.1
  426. curv_w = 8.0
  427. gcurv_w = 2.0
  429. step_size = 0.1
  430. step_num = 200
  432. averages = '4 2 1'
  434. initial_seg_prefix = '_initialization'
  435. seg_prefix = '_refined'
  437. if cpu_num > 0:
  438. cpu_str = ' -threads ' + str(cpu_num) + ' '
  439. else:
  440. cpu_str = ' '
  442. # Make output folders
  443. MDFA_Dir = segDir + 'md_plus_fa_maps/'
  444. if not path.exists(MDFA_Dir):
  445. system('mkdir ' + MDFA_Dir)
  447. probForceDir = segDir + 'probability_force_maps/'
  448. if not path.exists(probForceDir):
  449. system('mkdir ' + probForceDir)
  451. meshDir = segDir + 'mesh_output/'
  452. if not path.exists(meshDir):
  453. system('mkdir ' + meshDir)
  455. voxelDir = segDir + 'voxel_output/'
  456. if not path.exists(voxelDir):
  457. system('mkdir ' + voxelDir)
  459. atroposDir = segDir + 'atropos_hippoamyg_seg/'
  460. if not path.exists(atroposDir):
  461. system('mkdir ' + atroposDir)
  463. fsubcortseg_vtk = meshDir + subID + seg_prefix + '_subcortGM.vtk'
  464. if not path.exists(fsubcortseg_vtk):
  465. print('Subcortical Deformation Segmentation')
  467. # Make composite map of FA + csf probabilities. (This map defines the borders of the caudate and the thalamus)
  468. fmd_plusFA = MDFA_Dir + \
  469. path.basename(fmd.replace(fsl_ext(), '_plusFA' + fsl_ext()))
  470. system('fslmaths ' + fmd + ' -mul 1000 -add ' + ffa + ' ' + fmd_plusFA)
  471. MD_plusFA_data = nib.load(fmd_plusFA).get_fdata()
  473. fa_img = nib.load(ffa)
  475. harvard_img = nib.load(fharvard_native)
  476. harvard_data = harvard_img.get_fdata()
  478. dwi_img = nib.load(fdwi)
  479. dwi_data = dwi_img.get_fdata()
  481. # Probabilistic based force for subcortical structures
  482. glob_prob_force = np.zeros(dwi_data.shape)
  483. glob_pos_ind = [LEFT_WHITE_IDX, LEFT_PUT_IDX, LEFT_CAUDATE_IDX, LEFT_ACCUM_IDX,
  484. RIGHT_WHITE_IDX, RIGHT_PUT_IDX, RIGHT_CAUDATE_IDX, RIGHT_ACCUM_IDX]
  485. glob_neg_ind = [LEFT_GLOB_IDX, RIGHT_GLOB_IDX]
  486. for pos_ind in glob_pos_ind:
  487. glob_prob_force = glob_prob_force + harvard_data[:, :, :, pos_ind]
  489. for neg_ind in glob_neg_ind:
  490. glob_prob_force = glob_prob_force - harvard_data[:, :, :, neg_ind]
  492. fglob_prob_force = probForceDir + subID + '_glob_prob_force' + fsl_ext()
  493. nib.save(nib.Nifti1Image(glob_prob_force,
  494. fa_img.affine), fglob_prob_force)
  496. # Deform globus pallidus based on meanDWI and resticting movement into high FA regions
  497. fglobus_lh = initSegDir + subID + initial_seg_prefix + '_LEFT_GLOBUS.vtk'
  498. fglobus_lh_refined = meshDir + subID + seg_prefix + '_LEFT_GLOBUS.vtk'
  499. if not path.exists(fglobus_lh_refined):
  500. system('mirtk deform-mesh ' + fglobus_lh + ' ' + fglobus_lh_refined + ' -image ' + fdwi + ' -edge-distance 1.0 -edge-distance-averaging ' + averages + ' -edge-distance-smoothing 1 -edge-distance-median 1 -distance-image ' + fglob_prob_force + ' -distance 0.5 -distance-smoothing 1 -distance-averaging ' + averages + ' -distance-measure normal -optimizer EulerMethod -step ' + str(step_size) + ' -steps ' + str(
  501. step_num) + ' -epsilon 1e-6 -delta 0.001 -min-active 1% -reset-status -nointersection -fast-collision-test -min-width 0.01 -min-distance 0.01 -repulsion 4.0 -repulsion-distance 0.5 -repulsion-width 2.0 -curvature ' + str(curv_w) + ' -gauss-curvature ' + str(gcurv_w) + ' -edge-distance-type ClosestMaximum' + cpu_str + '-ascii -remesh 1 -min-edge-length ' + str(min_edgelength) + ' -max-edge-length ' + str(max_edgelength))
  502. sutil.surf_to_volume_mask(fdwi, fglobus_lh_refined, 1,
  503. voxelDir + subID + seg_prefix + '_LEFT_GLOBUS' + fsl_ext())
  505. fglobus_rh = initSegDir + subID + initial_seg_prefix + '_RIGHT_GLOBUS.vtk'
  506. fglobus_rh_refined = meshDir + subID + seg_prefix + '_RIGHT_GLOBUS.vtk'
  507. if not path.exists(fglobus_rh_refined):
  508. system('mirtk deform-mesh ' + fglobus_rh + ' ' + fglobus_rh_refined + ' -image ' + fdwi + ' -edge-distance 1.0 -edge-distance-averaging ' + averages + ' -edge-distance-smoothing 1 -edge-distance-median 1 -distance-image ' + fglob_prob_force + ' -distance 0.5 -distance-smoothing 1 -distance-averaging ' + averages + ' -distance-measure normal -optimizer EulerMethod -step ' + str(step_size) + ' -steps ' + str(step_num) +
  509. ' -epsilon 1e-6 -delta 0.001 -min-active 1% -reset-status -nointersection -fast-collision-test -min-width 0.01 -min-distance 0.01 -repulsion 4.0 -repulsion-distance 0.5 -repulsion-width 2.0 -curvature ' + str(curv_w) + ' -gauss-curvature ' + str(gcurv_w) + ' -edge-distance-type ClosestMaximum' + cpu_str + '-ascii -remesh 1 -min-edge-length ' + str(min_edgelength) + ' -max-edge-length ' + str(max_edgelength))
  510. sutil.surf_to_volume_mask(fdwi, fglobus_rh_refined, 1,
  511. voxelDir + subID + seg_prefix + '_RIGHT_GLOBUS' + fsl_ext())
  513. # Generate map for Striatum deformation
  514. fa_globus = np.zeros(MD_plusFA_data.shape)
  515. fa_globus[:] = MD_plusFA_data[:]
  517. lh_globus_data = nib.load(
  518. voxelDir + subID + seg_prefix + '_LEFT_GLOBUS' + fsl_ext()).get_fdata()
  519. rh_globus_data = nib.load(
  520. voxelDir + subID + seg_prefix + '_RIGHT_GLOBUS' + fsl_ext()).get_fdata()
  521. fa_globus[lh_globus_data == 1] = 2
  522. fa_globus[rh_globus_data == 1] = 2
  523. fmd_plusFA_globus = MDFA_Dir + \
  524. path.basename(fmd_plusFA.replace(fsl_ext(), '_globus' + fsl_ext()))
  525. if not path.exists(fmd_plusFA_globus):
  526. nib.save(nib.Nifti1Image(
  527. fa_globus, fa_img.affine), fmd_plusFA_globus)
  529. # Probabilistic atlas based force for striatum
  530. striatum_prob_force = np.zeros(dwi_data.shape)
  531. striatum_pos_ind = [LEFT_CORTEX_IDX, LEFT_WHITE_IDX,
  532. LEFT_GLOB_IDX, RIGHT_CORTEX_IDX, RIGHT_WHITE_IDX, RIGHT_GLOB_IDX]
  533. striatum_neg_ind = [LEFT_PUT_IDX, LEFT_CAUDATE_IDX, LEFT_ACCUM_IDX,
  534. RIGHT_PUT_IDX, RIGHT_CAUDATE_IDX, RIGHT_ACCUM_IDX]
  535. for pos_ind in striatum_pos_ind:
  536. striatum_prob_force = striatum_prob_force + \
  537. harvard_data[:, :, :, pos_ind]
  539. for neg_ind in striatum_neg_ind:
  540. striatum_prob_force = striatum_prob_force - \
  541. harvard_data[:, :, :, neg_ind]
  543. fstriatum_prob_force = probForceDir + subID + '_striatum_prob_force' + fsl_ext()
  544. if not path.exists(fstriatum_prob_force):
  545. nib.save(nib.Nifti1Image(striatum_prob_force,
  546. fa_img.affine), fstriatum_prob_force)
  548. fstriatum_lh = initSegDir + subID + initial_seg_prefix + '_LEFT_STRIATUM.vtk'
  549. fstriatum_lh_refined = meshDir + subID + seg_prefix + '_LEFT_STRIATUM.vtk'
  550. if not path.exists(fstriatum_lh_refined):
  551. system('mirtk deform-mesh ' + fstriatum_lh + ' ' + fstriatum_lh_refined + ' -image ' + fmd_plusFA_globus + ' -edge-distance 1.0 -edge-distance-averaging ' + averages + ' -edge-distance-smoothing 1 -edge-distance-median 1 -distance-image ' + fstriatum_prob_force + ' -distance 0.25 -distance-smoothing 1 -distance-averaging ' + averages + ' -distance-measure normal -optimizer EulerMethod -step ' + str(step_size) + ' -steps ' + str(step_num) +
  552. ' -epsilon 1e-6 -delta 0.001 -min-active 1% -reset-status -nointersection -fast-collision-test -min-width 0.01 -min-distance 0.01 -repulsion 4.0 -repulsion-distance 0.5 -repulsion-width 2.0 -curvature ' + str(curv_w) + ' -gauss-curvature ' + str(gcurv_w) + ' -edge-distance-type ClosestMaximum' + cpu_str + '-ascii -remesh 1 -min-edge-length ' + str(min_edgelength) + ' -max-edge-length ' + str(max_edgelength) + ' -edge-distance-min-intensity 0.3')
  553. sutil.surf_to_volume_mask(fdwi, fstriatum_lh_refined, 1,
  554. voxelDir + subID + seg_prefix + '_LEFT_STRIATUM' + fsl_ext())
  556. fstriatum_rh = initSegDir + subID + initial_seg_prefix + '_RIGHT_STRIATUM.vtk'
  557. fstriatum_rh_refined = meshDir + subID + seg_prefix + '_RIGHT_STRIATUM.vtk'
  558. if not path.exists(fstriatum_rh_refined):
  559. system('mirtk deform-mesh ' + fstriatum_rh + ' ' + fstriatum_rh_refined + ' -image ' + fmd_plusFA_globus + ' -edge-distance 1.0 -edge-distance-averaging ' + averages + ' -edge-distance-smoothing 1 -edge-distance-median 1 -distance-image ' + fstriatum_prob_force + ' -distance 0.25 -distance-smoothing 1 -distance-averaging ' + averages + ' -distance-measure normal -optimizer EulerMethod -step ' + str(step_size) + ' -steps ' + str(step_num) +
  560. ' -epsilon 1e-6 -delta 0.001 -min-active 1% -reset-status -nointersection -fast-collision-test -min-width 0.01 -min-distance 0.01 -repulsion 4.0 -repulsion-distance 0.5 -repulsion-width 2.0 -curvature ' + str(curv_w) + ' -gauss-curvature ' + str(gcurv_w) + ' -edge-distance-type ClosestMaximum' + cpu_str + '-ascii -remesh 1 -min-edge-length ' + str(min_edgelength) + ' -max-edge-length ' + str(max_edgelength) + ' -edge-distance-min-intensity 0.3')
  561. sutil.surf_to_volume_mask(fdwi, fstriatum_rh_refined, 1,
  562. voxelDir + subID + seg_prefix + '_RIGHT_STRIATUM' + fsl_ext())
  564. fa_md_striatum = np.zeros(MD_plusFA_data.shape)
  565. fa_md_striatum[:] = MD_plusFA_data[:]
  567. lh_striatum_data = nib.load(
  568. voxelDir + subID + seg_prefix + '_LEFT_STRIATUM' + fsl_ext()).get_fdata()
  569. rh_striatum_data = nib.load(
  570. voxelDir + subID + seg_prefix + '_RIGHT_STRIATUM' + fsl_ext()).get_fdata()
  571. fa_md_striatum[lh_striatum_data == 1] = 2
  572. fa_md_striatum[rh_striatum_data == 1] = 2
  573. fmd_plusFA_striatum = MDFA_Dir + \
  574. path.basename(fmd_plusFA.replace(
  575. fsl_ext(), '_striatum' + fsl_ext()))
  577. if not path.exists(fmd_plusFA_striatum):
  578. nib.save(nib.Nifti1Image(fa_md_striatum,
  579. fa_img.affine), fmd_plusFA_striatum)
  581. # Probabilistic atlas based force for thalamus
  582. thal_prob_force = np.zeros(dwi_data.shape)
  583. thal_pos_ind = [LEFT_CORTEX_IDX, LEFT_WHITE_IDX, LEFT_HIPPO_IDX,
  584. RIGHT_CORTEX_IDX, RIGHT_WHITE_IDX, RIGHT_HIPPO_IDX]
  585. thal_neg_ind = [LEFT_THAL_IDX, RIGHT_THAL_IDX]
  586. for pos_ind in thal_pos_ind:
  587. thal_prob_force = thal_prob_force + harvard_data[:, :, :, pos_ind]
  589. for neg_ind in thal_neg_ind:
  590. thal_prob_force = thal_prob_force - harvard_data[:, :, :, neg_ind]
  592. fthal_prob_force = probForceDir + subID + '_thalamus_prob_force' + fsl_ext()
  593. if not path.exists(fthal_prob_force):
  594. nib.save(nib.Nifti1Image(thal_prob_force,
  595. fa_img.affine), fthal_prob_force)
  597. fthalamus = initSegDir + subID + initial_seg_prefix + '_THALAMUS.vtk'
  598. fthalamus_refined = meshDir + subID + seg_prefix + '_THALAMUS.vtk'
  599. if not path.exists(fthalamus_refined):
  600. system('mirtk deform-mesh ' + fthalamus + ' ' + fthalamus_refined + ' -image ' + fmd_plusFA_striatum + ' -edge-distance 1.0 -edge-distance-averaging ' + averages + ' -edge-distance-smoothing 1 -edge-distance-median 1 -distance-image ' + fthal_prob_force + ' -distance 0.25 -distance-smoothing 1 -distance-averaging ' + averages + ' -distance-measure normal -optimizer EulerMethod -step ' + str(step_size) + ' -steps ' + str(step_num) +
  601. ' -epsilon 1e-6 -delta 0.001 -min-active 1% -reset-status -nointersection -fast-collision-test -min-width 0.01 -min-distance 0.01 -repulsion 4.0 -repulsion-distance 0.5 -repulsion-width 2.0 -curvature ' + str(curv_w) + ' -gauss-curvature ' + str(gcurv_w) + ' -edge-distance-type ClosestMaximum' + cpu_str + '-ascii -remesh 1 -min-edge-length ' + str(min_edgelength) + ' -max-edge-length ' + str(max_edgelength) + ' -edge-distance-min-intensity 1.25')
  603. thalamus_surf = sutil.read_surf_vtk(fthalamus_refined)
  604. [lh_thalamus, rh_thalamus] = sutil.split_surface_by_label(
  605. thalamus_surf, label=[LEFT_THAL_IDX, RIGHT_THAL_IDX], label_name='struct_label')
  607. fthalamus_lh_refined = meshDir + subID + seg_prefix + '_LEFT_THALAMUS.vtk'
  608. if not path.exists(fthalamus_lh_refined):
  609. sutil.write_surf_vtk(lh_thalamus, fthalamus_lh_refined)
  610. sutil.surf_to_volume_mask(fdwi, fthalamus_lh_refined, 1,
  611. voxelDir + subID + seg_prefix + '_LEFT_THALAMUS' + fsl_ext())
  613. fthalamus_rh_refined = meshDir + subID + seg_prefix + '_RIGHT_THALAMUS.vtk'
  614. if not path.exists(fthalamus_rh_refined):
  615. sutil.write_surf_vtk(rh_thalamus, fthalamus_rh_refined)
  616. sutil.surf_to_volume_mask(fdwi, fthalamus_rh_refined, 1,
  617. voxelDir + subID + seg_prefix + '_RIGHT_THALAMUS' + fsl_ext())
  619. # Hippocampus/Amygdala segmentation
  620. hipamyg_prob_force = np.zeros(dwi_data.shape)
  621. hipamyg_pos_ind = [LEFT_CORTEX_IDX, LEFT_WHITE_IDX, LEFT_THAL_IDX,
  622. RIGHT_CORTEX_IDX, RIGHT_WHITE_IDX, RIGHT_THAL_IDX, BRAIN_STEM_IDX]
  623. hipamyg_neg_ind = [LEFT_HIPPO_IDX, LEFT_AMYG_IDX,
  624. RIGHT_HIPPO_IDX, RIGHT_AMYG_IDX]
  625. for pos_ind in hipamyg_pos_ind:
  626. hipamyg_prob_force = hipamyg_prob_force + \
  627. harvard_data[:, :, :, pos_ind]
  629. for neg_ind in hipamyg_neg_ind:
  630. hipamyg_prob_force = hipamyg_prob_force - \
  631. harvard_data[:, :, :, neg_ind]
  633. fhipamyg_prob_force = probForceDir + subID + '_hipamyg_prob_force' + fsl_ext()
  634. nib.save(nib.Nifti1Image(hipamyg_prob_force,
  635. fa_img.affine), fhipamyg_prob_force)
  637. fhippoamyg_lh = initSegDir + subID + initial_seg_prefix + '_LEFT_HIPPOAMYG.vtk'
  638. fhippoamyg_lh_refined = atroposDir + subID + seg_prefix + '_LEFT_HIPPOAMYG.vtk'
  639. if not path.exists(fhippoamyg_lh_refined):
  640. system('mirtk deform-mesh ' + fhippoamyg_lh + ' ' + fhippoamyg_lh_refined + ' -image ' + fmd_plusFA + ' -edge-distance 1.0 -edge-distance-averaging ' + averages + ' -edge-distance-smoothing 1 -edge-distance-median 1 -distance-image ' + fhipamyg_prob_force + ' -distance 0.25 -distance-smoothing 1 -distance-averaging ' + averages + ' -distance-measure normal -optimizer EulerMethod -step ' + str(step_size) + ' -steps ' + str(step_num) +
  641. ' -epsilon 1e-6 -delta 0.001 -min-active 1% -reset-status -nointersection -fast-collision-test -min-width 0.01 -min-distance 0.01 -repulsion 4.0 -repulsion-distance 0.5 -repulsion-width 2.0 -curvature ' + str(curv_w) + ' -gauss-curvature ' + str(gcurv_w) + ' -edge-distance-type ClosestMaximum' + cpu_str + '-ascii -remesh 1 -min-edge-length ' + str(min_hippo_edgelength) + ' -max-edge-length ' + str(max_hippo_edgelength) + ' -edge-distance-min-intensity 1.0')
  642. sutil.surf_to_volume_mask(fdwi, fhippoamyg_lh_refined, 1,
  643. atroposDir + subID + seg_prefix + '_LEFT_HIPPOAMYG' + fsl_ext())
  645. fhippoamyg_rh = initSegDir + subID + initial_seg_prefix + '_RIGHT_HIPPOAMYG.vtk'
  646. fhippoamyg_rh_refined = atroposDir + subID + seg_prefix + '_RIGHT_HIPPOAMYG.vtk'
  647. if not path.exists(fhippoamyg_rh_refined):
  648. system('mirtk deform-mesh ' + fhippoamyg_rh + ' ' + fhippoamyg_rh_refined + ' -image ' + fmd_plusFA + ' -edge-distance 1.0 -edge-distance-averaging ' + averages + ' -edge-distance-smoothing 1 -edge-distance-median 1 -distance-image ' + fhipamyg_prob_force + ' -distance 0.25 -distance-smoothing 1 -distance-averaging ' + averages + ' -distance-measure normal -optimizer EulerMethod -step ' + str(step_size) + ' -steps ' + str(step_num) +
  649. ' -epsilon 1e-6 -delta 0.001 -min-active 1% -reset-status -nointersection -fast-collision-test -min-width 0.01 -min-distance 0.01 -repulsion 4.0 -repulsion-distance 0.5 -repulsion-width 2.0 -curvature ' + str(curv_w) + ' -gauss-curvature ' + str(gcurv_w) + ' -edge-distance-type ClosestMaximum' + cpu_str + '-ascii -remesh 1 -min-edge-length ' + str(min_hippo_edgelength) + ' -max-edge-length ' + str(max_hippo_edgelength) + ' -edge-distance-min-intensity 1.0')
  650. sutil.surf_to_volume_mask(fdwi, fhippoamyg_rh_refined, 1,
  651. atroposDir + subID + seg_prefix + '_RIGHT_HIPPOAMYG' + fsl_ext())
  653. # Run 3 channel tissue segmentation to separate amygdala from hippocampus
  654. # prepare prob maps and separate primary eigenvector into separate nifti file
  655. atropos_prefix_lh = atroposDir + subID + \
  656. initial_seg_prefix + '_INITAMYGHIPPO_LH'
  657. atropos_prefix_lh_label = atroposDir + subID + \
  658. initial_seg_prefix + '_INITAMYGHIPPO_LH_LABELS' + fsl_ext()
  659. atropos_prefix_lh_probout = atroposDir + subID + \
  660. initial_seg_prefix + '_INITAMYGHIPPO_LH_PROBOUT_'
  662. famyg_lh_prob = atropos_prefix_lh + '_01' + fsl_ext()
  663. if not path.exists(famyg_lh_prob):
  664. nib.save(nib.Nifti1Image(
  665. harvard_data[:, :, :, LEFT_AMYG_IDX]/100, harvard_img.affine), famyg_lh_prob)
  667. fhippo_lh_prob = atropos_prefix_lh + '_02' + fsl_ext()
  668. if not path.exists(fhippo_lh_prob):
  669. nib.save(nib.Nifti1Image(
  670. harvard_data[:, :, :, LEFT_HIPPO_IDX]/100, harvard_img.affine), fhippo_lh_prob)
  672. fcortex_lh_prob = atropos_prefix_lh + '_03' + fsl_ext()
  673. if not path.exists(fcortex_lh_prob):
  674. nib.save(nib.Nifti1Image(
  675. harvard_data[:, :, :, LEFT_CORTEX_IDX]/100, harvard_img.affine), fcortex_lh_prob)
  677. atropos_prefix_rh = atroposDir + subID + \
  678. initial_seg_prefix + '_INITAMYGHIPPO_RH'
  679. atropos_prefix_rh_label = atroposDir + subID + \
  680. initial_seg_prefix + '_INITAMYGHIPPO_RH_LABELS' + fsl_ext()
  681. atropos_prefix_rh_probout = atroposDir + subID + \
  682. initial_seg_prefix + '_INITAMYGHIPPO_RH_PROBOUT_'
  684. famyg_rh_prob = atropos_prefix_rh + '_01' + fsl_ext()
  685. if not path.exists(famyg_rh_prob):
  686. nib.save(nib.Nifti1Image(
  687. harvard_data[:, :, :, RIGHT_AMYG_IDX]/100, harvard_img.affine), famyg_rh_prob)
  689. fhippo_rh_prob = atropos_prefix_rh + '_02' + fsl_ext()
  690. if not path.exists(fhippo_rh_prob):
  691. nib.save(nib.Nifti1Image(
  692. harvard_data[:, :, :, RIGHT_HIPPO_IDX]/100, harvard_img.affine), fhippo_rh_prob)
  694. fcortex_rh_prob = atropos_prefix_rh + '_03' + fsl_ext()
  695. if not path.exists(fcortex_rh_prob):
  696. nib.save(nib.Nifti1Image(
  697. harvard_data[:, :, :, RIGHT_CORTEX_IDX]/100, harvard_img.affine), fcortex_rh_prob)
  699. PriorWeight = 0.3
  700. if not path.exists(atropos_prefix_lh_label):
  701. system('Atropos' +
  702. ' -a [' + fmd_plusFA + ']' +
  703. ' -x ' + atroposDir + subID + seg_prefix + '_LEFT_HIPPOAMYG' + fsl_ext() +
  704. ' -i PriorProbabilityImages[3, ' + atropos_prefix_lh + '_%02d' + fsl_ext() + ',' + str(PriorWeight) + ',0.0001]' +
  705. ' -m [0.3, 2x2x2] ' +
  706. ' --use-partial-volume-likelihoods false ' +
  707. ' -s 1x3 -s 1x2 ' +
  708. ' -o [' + atropos_prefix_lh_label + ',' + atropos_prefix_lh_probout + '%02d' + fsl_ext() + ']' +
  709. ' -k HistogramParzenWindows[1.0,32]' +
  710. ' -v 1')
  712. if not path.exists(atropos_prefix_rh_label):
  713. system('Atropos' +
  714. ' -a [' + fmd_plusFA + ']' +
  715. ' -x ' + atroposDir + subID + seg_prefix + '_RIGHT_HIPPOAMYG' + fsl_ext() +
  716. ' -i PriorProbabilityImages[3, ' + atropos_prefix_rh + '_%02d' + fsl_ext() + ',' + str(PriorWeight) + ',0.0001]' +
  717. ' -m [0.3, 2x2x2] ' +
  718. ' --use-partial-volume-likelihoods false ' +
  719. ' -s 1x3 -s 1x2 ' +
  720. ' -o [' + atropos_prefix_rh_label + ',' + atropos_prefix_rh_probout + '%02d' + fsl_ext() + ']' +
  721. ' -k HistogramParzenWindows[1.0,32]' +
  722. ' -v 1')
  724. lh_labels = nib.load(atropos_prefix_lh_label).get_fdata()
  725. rh_labels = nib.load(atropos_prefix_rh_label).get_fdata()
  727. # Generate new surfaces for hippocampus and amygdala
  728. fhippoamyg_lh_step1 = atroposDir + subID + \
  729. initial_seg_prefix + '_LEFT_HIPPOAMYG_STEP1'
  730. fsurf_list = extract_surfaces_from_labels(
  731. atropos_prefix_lh_label, [1, 2], atroposDir, fhippoamyg_lh_step1)
  732. famyg_lh_step1 = fsurf_list[0]
  733. fhippo_lh_step1 = fsurf_list[1]
  735. fhippoamyg_rh_step1 = atroposDir + subID + \
  736. initial_seg_prefix + '_RIGHT_HIPPOAMYG_STEP1'
  737. fsurf_list = extract_surfaces_from_labels(
  738. atropos_prefix_rh_label, [1, 2], atroposDir, fhippoamyg_rh_step1)
  739. famyg_rh_step1 = fsurf_list[0]
  740. fhippo_rh_step1 = fsurf_list[1]
  742. # Build amygdala prob force
  743. amyg_prob_force = np.zeros(dwi_data.shape)
  744. amyg_pos_ind = [LEFT_CORTEX_IDX, LEFT_WHITE_IDX, LEFT_THAL_IDX, LEFT_HIPPO_IDX, RIGHT_HIPPO_IDX,
  745. RIGHT_CORTEX_IDX, RIGHT_WHITE_IDX, RIGHT_THAL_IDX, BRAIN_STEM_IDX]
  746. amyg_neg_ind = [LEFT_AMYG_IDX, RIGHT_AMYG_IDX]
  747. for pos_ind in amyg_pos_ind:
  748. amyg_prob_force = amyg_prob_force + \
  749. harvard_data[:, :, :, pos_ind]
  751. for neg_ind in amyg_neg_ind:
  752. amyg_prob_force = amyg_prob_force - \
  753. harvard_data[:, :, :, neg_ind]
  755. famyg_prob_force = probForceDir + subID + '_amyg_prob_force' + fsl_ext()
  756. nib.save(nib.Nifti1Image(amyg_prob_force,
  757. fa_img.affine), famyg_prob_force)
  759. # Add hippocampus/cortex labels to md+fa map
  760. fa_md_hippocortex = np.zeros(MD_plusFA_data.shape)
  761. fa_md_hippocortex[:] = MD_plusFA_data[:]
  763. fa_md_hippocortex[lh_labels == 2] = 2
  764. fa_md_hippocortex[rh_labels == 2] = 2
  765. fa_md_hippocortex[lh_labels == 3] = 2
  766. fa_md_hippocortex[rh_labels == 3] = 2
  768. fmd_plusFA_hippocortex = MDFA_Dir + \
  769. path.basename(fmd_plusFA.replace(
  770. fsl_ext(), '_hippocortex' + fsl_ext()))
  771. if not path.exists(fmd_plusFA_hippocortex):
  772. nib.save(nib.Nifti1Image(fa_md_hippocortex,
  773. fa_img.affine), fmd_plusFA_hippocortex)
  775. famyg_lh_refined = meshDir + subID + seg_prefix + '_LEFT_AMYGDALA.vtk'
  776. if not path.exists(famyg_lh_refined):
  777. system('mirtk deform-mesh ' + famyg_lh_step1 + ' ' + famyg_lh_refined + ' -image ' + fmd_plusFA_hippocortex + ' -edge-distance 1.0 -edge-distance-averaging ' + averages + ' -edge-distance-smoothing 1 -edge-distance-median 1 -distance-image ' + famyg_prob_force + ' -distance 0.25 -distance-smoothing 1 -distance-averaging ' + averages + ' -distance-measure normal -optimizer EulerMethod -step ' + str(step_size) + ' -steps ' + str(step_num) +
  778. ' -epsilon 1e-6 -delta 0.001 -min-active 1% -reset-status -nointersection -fast-collision-test -min-width 0.01 -min-distance 0.01 -repulsion 4.0 -repulsion-distance 0.5 -repulsion-width 2.0 -curvature ' + str(curv_w) + ' -gauss-curvature ' + str(gcurv_w) + ' -edge-distance-type ClosestMaximum' + cpu_str + '-ascii -remesh 1 -min-edge-length ' + str(min_hippo_edgelength) + ' -max-edge-length ' + str(max_hippo_edgelength) + ' -edge-distance-min-intensity 1.0')
  779. system('mirtk project-onto-surface ' + famyg_lh_refined + ' ' + famyg_lh_refined +
  780. ' -constant ' + str(LEFT_AMYG_IDX) + ' -pointdata -name struct_label')
  781. sutil.surf_to_volume_mask(fdwi, famyg_lh_refined, 1,
  782. voxelDir + subID + seg_prefix + '_LEFT_AMYGDALA' + fsl_ext())
  784. famyg_rh_refined = meshDir + subID + seg_prefix + '_RIGHT_AMYGDALA.vtk'
  785. if not path.exists(famyg_rh_refined):
  786. system('mirtk deform-mesh ' + famyg_rh_step1 + ' ' + famyg_rh_refined + ' -image ' + fmd_plusFA_hippocortex + ' -edge-distance 1.0 -edge-distance-averaging ' + averages + ' -edge-distance-smoothing 1 -edge-distance-median 1 -distance-image ' + famyg_prob_force + ' -distance 0.25 -distance-smoothing 1 -distance-averaging ' + averages + ' -distance-measure normal -optimizer EulerMethod -step ' + str(step_size) + ' -steps ' + str(step_num) +
  787. ' -epsilon 1e-6 -delta 0.001 -min-active 1% -reset-status -nointersection -fast-collision-test -min-width 0.01 -min-distance 0.01 -repulsion 4.0 -repulsion-distance 0.5 -repulsion-width 2.0 -curvature ' + str(curv_w) + ' -gauss-curvature ' + str(gcurv_w) + ' -edge-distance-type ClosestMaximum' + cpu_str + '-ascii -remesh 1 -min-edge-length ' + str(min_hippo_edgelength) + ' -max-edge-length ' + str(max_hippo_edgelength) + ' -edge-distance-min-intensity 1.0')
  788. system('mirtk project-onto-surface ' + famyg_rh_refined + ' ' + famyg_rh_refined +
  789. ' -constant ' + str(RIGHT_AMYG_IDX) + ' -pointdata -name struct_label')
  790. sutil.surf_to_volume_mask(fdwi, famyg_rh_refined, 1,
  791. voxelDir + subID + seg_prefix + '_RIGHT_AMYGDALA' + fsl_ext())
  793. # Build amygdala prob force
  794. hippo_prob_force = np.zeros(dwi_data.shape)
  795. hippo_pos_ind = [LEFT_CORTEX_IDX, LEFT_WHITE_IDX, LEFT_THAL_IDX, LEFT_AMYG_IDX, RIGHT_AMYG_IDX,
  796. RIGHT_CORTEX_IDX, RIGHT_WHITE_IDX, RIGHT_THAL_IDX, BRAIN_STEM_IDX]
  797. hippo_neg_ind = [LEFT_HIPPO_IDX, RIGHT_HIPPO_IDX]
  798. for pos_ind in hippo_pos_ind:
  799. hippo_prob_force = hippo_prob_force + \
  800. harvard_data[:, :, :, pos_ind]
  802. for neg_ind in hippo_neg_ind:
  803. hippo_prob_force = hippo_prob_force - \
  804. harvard_data[:, :, :, neg_ind]
  806. fhippo_prob_force = probForceDir + subID + '_hippo_prob_force' + fsl_ext()
  807. nib.save(nib.Nifti1Image(hippo_prob_force,
  808. fa_img.affine), fhippo_prob_force)
  810. # Add hippocampus/cortex labels to md+fa map
  811. fa_md_amyg = np.zeros(MD_plusFA_data.shape)
  812. fa_md_amyg[:] = MD_plusFA_data[:]
  814. fa_md_amyg[lh_labels == 1] = 2
  815. fa_md_amyg[rh_labels == 1] = 2
  816. fa_md_amyg[lh_labels == 3] = 2
  817. fa_md_amyg[rh_labels == 3] = 2
  819. fmd_plusFA_amyg = MDFA_Dir + \
  820. path.basename(fmd_plusFA.replace(
  821. fsl_ext(), '_amyg' + fsl_ext()))
  822. if not path.exists(fmd_plusFA_amyg):
  823. nib.save(nib.Nifti1Image(fa_md_amyg,
  824. fa_img.affine), fmd_plusFA_amyg)
  826. fhippo_lh_refined = meshDir + subID + seg_prefix + '_LEFT_HIPPO.vtk'
  827. if not path.exists(fhippo_lh_refined):
  828. system('mirtk deform-mesh ' + fhippo_lh_step1 + ' ' + fhippo_lh_refined + ' -image ' + fmd_plusFA_amyg + ' -edge-distance 1.0 -edge-distance-averaging ' + averages + ' -edge-distance-smoothing 1 -edge-distance-median 1 -distance-image ' + fhippo_prob_force + ' -distance 0.25 -distance-smoothing 1 -distance-averaging ' + averages + ' -distance-measure normal -optimizer EulerMethod -step ' + str(step_size) + ' -steps ' + str(step_num) +
  829. ' -epsilon 1e-6 -delta 0.001 -min-active 1% -reset-status -nointersection -fast-collision-test -min-width 0.01 -min-distance 0.01 -repulsion 4.0 -repulsion-distance 0.5 -repulsion-width 2.0 -curvature ' + str(curv_w) + ' -gauss-curvature ' + str(gcurv_w) + ' -edge-distance-type ClosestMaximum' + cpu_str + '-ascii -remesh 1 -min-edge-length ' + str(min_hippo_edgelength) + ' -max-edge-length ' + str(max_hippo_edgelength) + ' -edge-distance-min-intensity 1.0')
  830. system('mirtk project-onto-surface ' + fhippo_lh_refined + ' ' + fhippo_lh_refined +
  831. ' -constant ' + str(LEFT_HIPPO_IDX) + ' -pointdata -name struct_label')
  832. sutil.surf_to_volume_mask(fdwi, fhippo_lh_refined, 1,
  833. voxelDir + subID + seg_prefix + '_LEFT_HIPPO' + fsl_ext())
  835. fhippo_rh_refined = meshDir + subID + seg_prefix + '_RIGHT_HIPPO.vtk'
  836. if not path.exists(fhippo_rh_refined):
  837. system('mirtk deform-mesh ' + fhippo_rh_step1 + ' ' + fhippo_rh_refined + ' -image ' + fmd_plusFA_amyg + ' -edge-distance 1.0 -edge-distance-averaging ' + averages + ' -edge-distance-smoothing 1 -edge-distance-median 1 -distance-image ' + fhippo_prob_force + ' -distance 0.25 -distance-smoothing 1 -distance-averaging ' + averages + ' -distance-measure normal -optimizer EulerMethod -step ' + str(step_size) + ' -steps ' + str(step_num) +
  838. ' -epsilon 1e-6 -delta 0.001 -min-active 1% -reset-status -nointersection -fast-collision-test -min-width 0.01 -min-distance 0.01 -repulsion 4.0 -repulsion-distance 0.5 -repulsion-width 2.0 -curvature ' + str(curv_w) + ' -gauss-curvature ' + str(gcurv_w) + ' -edge-distance-type ClosestMaximum' + cpu_str + '-ascii -remesh 1 -min-edge-length ' + str(min_hippo_edgelength) + ' -max-edge-length ' + str(max_hippo_edgelength) + ' -edge-distance-min-intensity 1.0')
  839. system('mirtk project-onto-surface ' + fhippo_rh_refined + ' ' + fhippo_rh_refined +
  840. ' -constant ' + str(RIGHT_HIPPO_IDX) + ' -pointdata -name struct_label')
  841. sutil.surf_to_volume_mask(fdwi, fhippo_rh_refined, 1,
  842. voxelDir + subID + seg_prefix + '_RIGHT_HIPPO' + fsl_ext())
  844. # Ventricle segmentation
  845. # Probabilistic based force for subcortical structures
  846. vent_prob_force = np.zeros(dwi_data.shape)
  847. vent_pos_ind = [LEFT_WHITE_IDX, LEFT_PUT_IDX, LEFT_CAUDATE_IDX, LEFT_ACCUM_IDX, LEFT_GLOB_IDX, LEFT_THAL_IDX,
  848. RIGHT_WHITE_IDX, RIGHT_PUT_IDX, RIGHT_CAUDATE_IDX, RIGHT_ACCUM_IDX, RIGHT_GLOB_IDX, RIGHT_THAL_IDX]
  850. vent_neg_ind = [LEFT_VENT_IDX, RIGHT_VENT_IDX]
  852. for pos_ind in vent_pos_ind:
  853. vent_prob_force = vent_prob_force + harvard_data[:, :, :, pos_ind]
  855. for neg_ind in vent_neg_ind:
  856. vent_prob_force = vent_prob_force - harvard_data[:, :, :, neg_ind]
  858. fvent_prob_force = probForceDir + subID + '_vent_prob_force' + fsl_ext()
  859. nib.save(nib.Nifti1Image(vent_prob_force,
  860. fa_img.affine), fvent_prob_force)
  862. # Deform ventricles on mean DWI
  863. fvent = initSegDir + subID + initial_seg_prefix + '_VENTRICLES.vtk'
  864. fvent_refined = meshDir + subID + seg_prefix + '_VENTRICLES.vtk'
  865. if not path.exists(fvent_refined):
  866. system('mirtk deform-mesh ' + fvent + ' ' + fvent_refined + ' -image ' + fdwi + ' -edge-distance 1.0 -edge-distance-averaging ' + averages + ' -edge-distance-smoothing 1 -edge-distance-median 1 -distance-image ' + fvent_prob_force + ' -distance 0.5 -distance-smoothing 1 -distance-averaging ' + averages + ' -distance-measure normal -optimizer EulerMethod -step ' + str(step_size) + ' -steps ' + str(
  867. step_num) + ' -epsilon 1e-6 -delta 0.001 -min-active 1% -reset-status -nointersection -fast-collision-test -min-width 0.01 -min-distance 0.01 -repulsion 4.0 -repulsion-distance 0.5 -repulsion-width 2.0 -curvature ' + str(curv_w) + ' -gauss-curvature ' + str(gcurv_w) + ' -edge-distance-type ClosestMaximum' + cpu_str + '-ascii -remesh 1 -min-edge-length ' + str(min_edgelength) + ' -max-edge-length ' + str(max_edgelength))
  868. system('mirtk project-onto-surface ' + fvent_refined + ' ' + fvent_refined +
  869. ' -constant ' + str(VENT_IDX) + ' -pointdata -name struct_label')
  870. sutil.surf_to_volume_mask(fdwi, fvent_refined, 1,
  871. voxelDir + subID + seg_prefix + '_VENTRICLES' + fsl_ext())
  873. # Append subcort surfs together into single vtk file
  874. fsubcortseg_vtk = meshDir + subID + seg_prefix + '_subcortGM.vtk'
  875. if not path.exists(fsubcortseg_vtk):
  876. appender = vtk.vtkAppendPolyData()
  877. appender.AddInputData(sutil.read_surf_vtk(fglobus_lh_refined))
  878. appender.AddInputData(sutil.read_surf_vtk(fglobus_rh_refined))
  879. appender.AddInputData(sutil.read_surf_vtk(fstriatum_lh_refined))
  880. appender.AddInputData(sutil.read_surf_vtk(fstriatum_rh_refined))
  881. appender.AddInputData(sutil.read_surf_vtk(fthalamus_lh_refined))
  882. appender.AddInputData(sutil.read_surf_vtk(fthalamus_rh_refined))
  883. appender.AddInputData(sutil.read_surf_vtk(famyg_rh_refined))
  884. appender.AddInputData(sutil.read_surf_vtk(famyg_lh_refined))
  885. appender.AddInputData(sutil.read_surf_vtk(fhippo_rh_refined))
  886. appender.AddInputData(sutil.read_surf_vtk(fhippo_lh_refined))
  887. appender.AddInputData(sutil.read_surf_vtk(fvent_refined))
  888. appender.Update()
  890. deepGM_surf = appender.GetOutput()
  891. sutil.write_surf_vtk(deepGM_surf, fsubcortseg_vtk)
  892. else:
  893. print('Subcortical Segmentation already performed: Skipping')
  895. return meshDir, voxelDir
  898. def segment(procDir, subID, preproc_suffix, shell_suffix, cpu_num=0):
  899. """
  900. Subcortical brain segmentation with 3D surface based deformation using DTI images/maps
  902. Parameters
  903. ----------
  904. procDir: string
  905. parent director containing subject
  906. subID: string
  907. the microbrain subject to be processed
  908. preproc_suffix: string
  909. suffix detailing what preprocessing has been performed
  910. shell_suffix: string
  911. suffix detailing what shells were used when running microbrain
  913. Optional Parameters
  914. -------------------
  915. cpu_num: integer
  916. number of threads to use for computationally intense tasks
  918. Returns
  919. -------
  920. none
  921. """
  923. subDir = procDir + '/' + subID
  924. segDir = subDir + '/subcortical_segmentation/'
  925. regDir = subDir + '/registration/'
  927. if preproc_suffix == '':
  928. suffix = '_' + shell_suffix
  929. fdwi = subDir + '/meanDWI/' + subID + '_mean_b' + \
  930. shell_suffix.split('b')[-1] + '_n4' + fsl_ext()
  931. else:
  932. suffix = '_' + preproc_suffix + '_' + shell_suffix
  933. fdwi = subDir + '/meanDWI/' + subID + '_' + \
  934. preproc_suffix + '_mean_b' + \
  935. shell_suffix.split('b')[-1] + '_n4' + fsl_ext()
  937. ffa = subDir + '/DTI_maps/' + subID + suffix + '_FA' + fsl_ext()
  938. fmd = subDir + '/DTI_maps/' + subID + suffix + '_MD' + fsl_ext()
  940. # Register probability maps
  941. ftemplate = get_fsl_standard_dir() + 'FSL_HCP1065_FA_1mm.nii.gz'
  943. fharvard = get_fsl_atlas_dir() + 'HarvardOxford/HarvardOxford-sub-prob-1mm.nii.gz'
  944. fharvard_native = register_probatlas_to_native(
  945. ffa, ftemplate, fharvard, regDir, cpu_num=cpu_num)
  947. # Parent Directory to store all output
  948. if not path.exists(segDir):
  949. system('mkdir ' + segDir)
  951. # Directory to store initial meshes
  952. initSegDir = segDir + 'structure_initialization/'
  953. if not path.exists(initSegDir):
  954. system('mkdir ' + initSegDir)
  956. initial_voxel_labels_from_harvard(
  957. fharvard_native, subID, initSegDir, md_file=fmd, fa_file=ffa)
  959. meshDir, voxelDir, = deform_subcortical_surfaces(
  960. fdwi, ffa, fmd, fharvard_native, segDir, initSegDir, subID, cpu_num=cpu_num)
  962. return meshDir, voxelDir

mbrain_segment.py at commit 47d9efb, under MIT · at the source

Overview

Authors: Graham Little1,2, J. Alejandro Acosta‐Franco1, Christian Beaulieu1
ORCID iDs: Graham Little, Christian Beaulieu
  1. Department of Radiology and Diagnostic Imaging & Biomedical Engineering University of Alberta Edmonton Alberta Canada
  2. Department of Computer Science Université de Sherbrooke Sherbrooke Quebec Canada
Institutions: Université de Sherbrooke (Canada); University of Alberta (Canada)
Journal: NMR in biomedicine, volume 39, issue 8, article e70353
Dates: received 13 November 2025; accepted 16 June 2026; published online 9 July 2026; in print August 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1002/nbm.70353 · PMID 42423344 · PMCID PMC13348013 · OpenAlex W4387453483
Open access: hybrid, a free copy (OpenAlex)
Preprint: doi:10.1101/2023.10.04.560912
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism)
Methods: Statistics, Machine learning, fMRI & imaging
MeSH: Aging*, Diffusion Tensor Imaging*, Gray Matter*, Adolescent, Adult, Aged, Algorithms, Automation, Female, Humans, Male, Middle Aged, Young Adult (* major topic)
Topic: Advanced Neuroimaging Techniques and Applications (Radiology, Nuclear Medicine and Imaging, Medicine), according to OpenAlex
Funding: University Hospital Foundation; Canadian Institutes of Health Research; Women's and Children's Health Research Institute; Canada Research Chairs; Natural Sciences and Engineering Research Council of Canada; Unifying Neuroscience and Artificial Intelligence in Quebec
Citations: not cited yet (Europe PMC); 66 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.

Repository

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

LittleBrainLab/microbrain

License: MIT
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 47d9efbc459ea480d28fa7f6cc5c4f594008009a, 22 January 2025
Languages: Python (19)
Size: 48 files, 19 scripts
Software Heritage: not archived
Found in: the text, “Surface‐Based Deep GM Segmentation on DTI”
Holds: README, environment (requirements.txt, setup.cfg, setup.py)
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (10 files), NiBabel (9 files), FSL (7 files), DIPY (6 files), ANTs (4 files), SciPy (2 files), dcm2niix (1 file), FreeSurfer (1 file), scikit-image (1 file), Connectome Workbench (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
20 files
At the source: github.com/LittleBrainLab/microbrain

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;
  • 19 scripts, each with its path and the digest of its content;
  • 6 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

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

Data

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

Data 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:

  • it says that the data are available on request

Read it in the paper: doi.org/10.1002/nbm.70353.

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 3, 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, 13 MeSH terms, 6 funders, 66 references.

Cite

This paper

Little, G., Acosta‐Franco, J. A., & Beaulieu, C. (2026). Automated Surface-Based Segmentation of Deep Gray Matter Regions Based on Diffusion Tensor Images Reveals Unique Age Trajectories Over the Healthy Lifespan. NMR in biomedicine, 39(8), e70353. https://doi.org/10.1002/nbm.70353

BibTeX

@article{little2026automated,
author = {Little, Graham and Acosta‐Franco, J. Alejandro and Beaulieu, Christian},
title = {{Automated Surface-Based Segmentation of Deep Gray Matter Regions Based on Diffusion Tensor Images Reveals Unique Age Trajectories Over the Healthy Lifespan}},
journal = {NMR in biomedicine},
year = {2026},
month = aug,
volume = {39},
number = {8},
pages = {e70353},
publisher = {Wiley},
issn = {0952-3480},
doi = {10.1002/nbm.70353},
url = {https://doi.org/10.1002/nbm.70353},
pmid = {42423344},
pmcid = {PMC13348013}
}

RIS

TY - JOUR
AU - Little, Graham
AU - Acosta‐Franco, J. Alejandro
AU - Beaulieu, Christian
TI - Automated Surface-Based Segmentation of Deep Gray Matter Regions Based on Diffusion Tensor Images Reveals Unique Age Trajectories Over the Healthy Lifespan
T2 - NMR in biomedicine
J2 - NMR Biomed
PY - 2026
DA - 2026/08/01
VL - 39
IS - 8
SP - e70353
SN - 0952-3480
PB - Wiley
DO - 10.1002/nbm.70353
UR - https://doi.org/10.1002/nbm.70353
LA - en
ER -

CSL-JSON

{
"id": "10.1002/nbm.70353",
"type": "article-journal",
"title": "Automated Surface-Based Segmentation of Deep Gray Matter Regions Based on Diffusion Tensor Images Reveals Unique Age Trajectories Over the Healthy Lifespan",
"container-title": "NMR in biomedicine",
"author": [
{
"family": "Little",
"given": "Graham"
},
{
"family": "Acosta‐Franco",
"given": "J. Alejandro"
},
{
"family": "Beaulieu",
"given": "Christian"
}
],
"container-title-short": "NMR Biomed",
"volume": "39",
"issue": "8",
"page": "e70353",
"DOI": "10.1002/nbm.70353",
"PMID": "42423344",
"PMCID": "PMC13348013",
"ISSN": "0952-3480",
"publisher": "Wiley",
"URL": "https://doi.org/10.1002/nbm.70353",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
1
]
]
}
}

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.70355
High-Resolution Diffusion Kurtosis Imaging of Hippocampus Subfields Across the Healthy Lifespan.
Journal: NMR in biomedicine
In common: structural MRI / diffusion, 11 references, author Christian Beaulieu
[2] doi:10.1371/journal.pbio.3003856 [code]
Aging and metabolism contribute separately to brain-body health.
Journal: PLoS biology
In common: Connectome Workbench, ANTs, FreeSurfer, 4 other tools, structural MRI / diffusion, 6 references
[3] doi:10.1162/imag.a.1183 [code]
Learning-based segmentation of diffusion-weighted MR images with arbitrary <i>q</i>-space samplings.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: scikit-image, NiBabel, SciPy, 1 other tool, structural MRI / diffusion, 9 references
[4] doi:10.7554/elife.108408 [code]
Frequency and laminar profile of feature-specific visual activity revealed by interleaved EEG-fMRI.
Journal: eLife
In common: DIPY, Connectome Workbench, ANTs, 6 other tools, 2 references
[5] doi:10.1038/s41467-026-71719-y [code]
Brain functional-structural gradient coupling reflects development, behavior and genetic influences.
Journal: Nature communications
In common: dcm2niix, DIPY, ANTs, 5 other tools, 2 references
[6] doi:10.1016/j.ynirp.2026.100360 [code]
CIVET-Chimp: An automated pipeline for MRI-based cortical surface extraction in chimpanzees.
Journal: Neuroimage. Reports
In common: Connectome Workbench, ANTs, FreeSurfer, 2 other tools, structural MRI / diffusion, 5 references
[7] doi:10.1016/j.crmeth.2026.101473 [code]
AmygdalaGo-BOLT for boundary-aware segmentation of the human amygdala.
Journal: Cell reports methods
In common: Connectome Workbench, ANTs, FreeSurfer, 5 other tools, structural MRI / diffusion, 2 references
[8] 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: DIPY, ANTs, FreeSurfer, 3 other tools, structural MRI / diffusion, 5 references
[9] doi:10.1038/s41598-026-54446-8 [code]
Deep learning-based Desikan-Killiany parcellation of the brain using diffusion MRI.
Journal: Scientific reports
In common: DIPY, FreeSurfer, scikit-image, 3 other tools, structural MRI / diffusion, 4 references
[10] doi:10.1038/s41467-026-71918-7 [code]
Developmental disinhibition gates language lateralization in childhood.
Journal: Nature communications
In common: dcm2niix, ANTs, FreeSurfer, 4 other tools, 3 references

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.