OSCR

Remote cortical degeneration related to structural connectivity following recent small subcortical infarcts.

Code ↔ Paper

7 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 7 matches
  1. [1] § Methods › Neuroimaging data processing ↔ cvdproc/pipelines/smri/fsl_anat_pipeline.py, lines 9 to 76 · score 0.68 · fsl_anat, Structural MRI, MNI space, tissue, pipeline, preprocessing · Report a mismatch
  2. [2] § Methods › Estimation of lesion-connected cortical region › Comparison between direct and indirect approaches ↔ cvdproc/pipelines/dmri/fdt/fdt_nipype.py, lines 1079 to 1160 · score 0.67 · medium connectivity, low connectivity, high connectivity, sulc, curv, thickness · Report a mismatch
  3. [3] § Methods › Estimation of lesion-connected cortical region › Direct estimation ↔ cvdproc/pipelines/dmri/dwi_pipeline.py, lines 1564 to 1630 · score 0.63 · Probabilistic tractography, seed mask, medium, native, threshold, space · Report a mismatch
  4. [4] § Methods › Neuroimaging data processing ↔ cvdproc/pipelines/smri/scn_pipeline.py, lines 10 to 70 · score 0.62 · structural covariance network, FreeSurfer, recon, pipeline, field · Report a mismatch
  5. [5] § Methods › cSVD score and vascular risk factor (VRF) score ↔ cvdproc/pipelines/smri/csvd_quantification/wmh/wmh_location_nipype.py, lines 40 to 114 · score 0.58 · periventricular WMH, deep WMH, Fazekas, cSVD · Report a mismatch
  6. [6] § Methods › Neuroimaging data processing ↔ cvdproc/pipelines/dmri/dwi_pipeline.py, lines 676 to 704 · score 0.57 · denoising, Synb0, MRtrix3, TOPUP, eddy, EPI · Report a mismatch
  7. [7] § Methods › Neuroimaging data processing ↔ cvdproc/pipelines/bash/synb0/run_synb0.sh, lines 64 to 143 · score 0.56 · Synb0 DISCO, EPI, FSL, undistorted, field, MRI · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 2,919 lines · 195 KB · MIT · 2 matches

  1. import os
  2. import pandas as pd
  3. import json
  4. import nibabel as nib
  5. import numpy as np
  6. import glob
  7. from cvdproc.config.paths import get_package_path
  8. from nipype import Node, Workflow, MapNode, Function
  9. from nipype.interfaces.utility import IdentityInterface, Merge, Select
  10. from cvdproc.bids_data.rename_bids_file import rename_bids_file
  11. # nipype built-in interfaces
  12. from nipype.interfaces.fsl import FLIRT, ExtractROI, ConvertXFM, MultiImageMaths
  13. from nipype.interfaces.freesurfer import MRIConvert
  14. from nipype.interfaces.mrtrix3 import Tractography, MRConvert, DWIPreproc, BuildConnectome, LabelConvert
  15. # custom interfaces
  16. from cvdproc.pipelines.common.files import FilterExisting, MergeFilename
  17. from cvdproc.pipelines.common.pad_dwi import PadDWI
  18. from cvdproc.pipelines.common.flip_b_table import FlipBTable
  19. from cvdproc.pipelines.common.register import Tkregister2fs2t1w, MRIConvertApplyWarp, MRIVol2Vol
  20. from cvdproc.pipelines.common.dwi_resample import ANTsResampleDWI
  21. from cvdproc.pipelines.common.image_calc import RemoveMaskRegion, CalculateScalarMaps
  22. from cvdproc.pipelines.dmri.fdt.fdt_nipype import B0AllAndAcqparam, IndexTxt, Topup, EddyCuda, OrderEddyOutputs, B0RefAndBrainMask, DTIFitBIDS, BedpostxGPUCustom, Probtrackx, ReadSingleValue, ExtractSurfaceParameters, ApplyGiiMaskToMgh, ApplyIndexTxtToMgh, DefineConnectionLevel
  23. from cvdproc.pipelines.dmri.synb0.synb0_nipype import Synb0
  24. from cvdproc.pipelines.dmri.psmd.psmd_nipype import PSMDCommandLine
  25. from cvdproc.pipelines.dmri.alps.alps_nipype import ALPS
  26. from cvdproc.pipelines.dmri.dsistudio.create_src_nipype import CreateSRC
  27. from cvdproc.pipelines.dmri.dsistudio.reconstruction_nipype import DSIstudioReconstruction, DSIstudioReconstructionOldVer, DSIstudioReconstruction2, DSIStudioBvecToFSL
  28. from cvdproc.pipelines.dmri.dsistudio.tracking_nipype import DSIstudioTracking
  29. from cvdproc.pipelines.dmri.dsistudio.tract_analysis import TractStatsInterface
  30. from cvdproc.pipelines.dmri.dsistudio.fixing_nipype import FixDSIStudioDWIAndBvec
  31. from cvdproc.pipelines.dmri.pved.pved_nipype import PVeD, PVeDAtropos, PVeDPrepareMNI
  32. from cvdproc.pipelines.dmri.mrtrix3.tcksample_nipype import TckSampleCommand, CalculateMeanTckSample
  33. from cvdproc.pipelines.dmri.mrtrix3.denoise_degibbs_nipype import MrtrixDenoise, MrtrixDegibbs
  34. from cvdproc.pipelines.dmri.mrtrix3.connectome_nipype import RemoveWMH, ConnectomePrepare
  35. from cvdproc.pipelines.dmri.mrtrix3.utils import RemoveTractRegion
  36. from cvdproc.pipelines.dmri.stats.dti_scalar_maps import GenerateNAWMMask, CalculateTDIWeightedScalars, TckSampleMultiScalarProfile, TckSampleMultiScalarBundle
  37. from cvdproc.pipelines.dmri.dipy.dipy_freewater_dti import FreeWaterTensor
  38. from cvdproc.pipelines.dmri.dipy.dipy_dki import DKIFit
  39. from cvdproc.pipelines.dmri.dipy.dipy_degibbs import DipyDegibbs
  40. from cvdproc.pipelines.dmri.freewater.single_shell_freewater import SingleShellFW
  41. from cvdproc.pipelines.dmri.freewater.markvcid_freewater import MarkVCIDFreeWater
  42. from cvdproc.pipelines.dmri.qsiprep.register import QSIPrepOrigToACPC
  43. from cvdproc.pipelines.dmri.visualization.vis import PlotTckOnSliceCmd
  44. from cvdproc.pipelines.smri.mirror.mirror_nipype import MirrorMask
  45. from cvdproc.pipelines.smri.fsl.fsl_anat_nipype import FSLANAT
  46. from cvdproc.pipelines.smri.freesurfer.synthstrip import SynthStrip
  47. from cvdproc.pipelines.smri.freesurfer.utils import MRIBinarize, MRIvol2surf, MergeRibbon
  48. class DWIPipeline:
  49. def __init__(self,
  50. subject: object,
  51. session: object,
  52. output_path: str,
  53. use_which_dwi: str = None,
  54. use_which_t1w: str = None,
  55. use_which_flair: str = None,
  56. use_freesurfer_longitudinal: bool = False,
  57. preprocess: bool = False,
  58. output_resolution: float = 2.0,
  59. degibbs: bool = True,
  60. flip_b_table_axis: list = [],
  61. preprocess_method: str = 'fdt',
  62. synb0: bool = False,
  63. use_which_reverse_b0: str = None,
  64. dti_fit: bool = False,
  65. dti_low_b: bool = True,
  66. dwi_t1w_register: bool = False,
  67. dki_fit: bool = False,
  68. dsistudio_dti: bool = False,
  69. dsistudio_gqi: bool = False,
  70. dsistudio_qsdr: bool = False,
  71. amico_noddi: bool = False,
  72. connectome: list = [],
  73. tractography: list = [],
  74. seed_mask: str = 'lesion_mask',
  75. use_which_mask: str = None,
  76. dtialps: bool = False,
  77. dtialps_register_method: int = 1,
  78. pved: bool = False,
  79. pved_method: str = 'v2',
  80. freewater: list = [],
  81. psmd: bool = False,
  82. psmd_exclude_seed_mask: bool = False,
  83. visual_pathway_analysis: bool = False,
  84. visual_pathway_analysis_flip_y: bool = False,
  85. calculate_dwi_metrics: bool = False,
  86. use_synthseg_wm_mask: bool = False,
  87. exclude_seed_mask: bool = True,
  88. exclude_wmh_mask: bool = False,
  89. keep_preproc_intermediate: bool = True,
  90. extract_from: str = None,
  91. **kwargs):
  92. """
  93. DWI pipeline
  94. Args:
  95. subject (object): Subject object
  96. session (object): Session object
  97. output_path (str): Output path
  98. use_which_dwi (str, optional): Which DWI file to use. Defaults to None.
  99. use_which_t1w (str, optional): Which T1w file to use. Defaults to None.
  100. use_which_flair (str, optional): Which FLAIR file to use. Defaults to None.
  101. use_freesurfer_longitudinal (bool, optional): Whether to use FreeSurfer longitudinal processing outpus. Defaults to False.
  102. preprocess (bool, optional): Whether to preprocess the DWI data. Defaults to False.
  103. output_resolution (float, optional): Output resolution in mm. Defaults to 2.0.
  104. degibbs (bool, optional): Whether to perform degibbsing. Defaults to True.
  105. flip_b_table_axis (list, optional): List of axes to flip the b-table. Defaults to []. [0], [1], [2] for x, y, z axes.
  106. preprocess_method (str, optional): Preprocessing method. 'fdt', 'mrtrix3', 'dsistudio_mc' or 'post_qsiprep'. Defaults to 'fdt'.
  107. synb0 (bool, optional): Whether to use Synb0 for TOPUP. Defaults to False.
  108. use_which_reverse_b0 (str, optional): Which reverse coded b0 image to use. Defaults to None.
  109. dti_fit (bool, optional): Whether to fit DTI model using FSL dtifit. Defaults to False.
  110. dti_low_b (bool, optional): Whether to only use b<=1500 volumes for DTI fitting. Defaults to True.
  111. dwi_t1w_register (bool, optional): Whether to register DWI to T1w space. Defaults to False.
  112. dki_fit (bool, optional): Whether to fit DKI model using Dipy. Defaults to False.
  113. dsistudio_dti (bool, optional): Whether to fit DTI model using DSI Studio. Defaults to False.
  114. dsistudio_gqi (bool, optional): Whether to fit GQI model using DSI Studio. Defaults to False.
  115. dsistudio_qsdr (bool, optional): Whether to fit QSDR model using DSI Studio. Defaults to False.
  116. amico_noddi (bool, optional): Whether to fit NODDI model using AMICO. Defaults to False.
  117. connectome (list, optional): Whether to create a connectome (Global). Defaults to []. Subsets of ['fdt', 'mrtrix3', 'dsistudio'].
  118. tractography (list, optional): Tractography method (ROI-based). Defaults to []. Subsets of ['fdt', 'mrtrix3', 'dsistudio'].
  119. seed_mask (str, optional): Folder name for seed mask. Expected to find <bids_dir>/derivatives/<seed_mask>/sub-<id>/ses-<id>/*<use_which_mask>*.nii.gz. Defaults to 'lesion_mask'. Must be in the T1w space.
  120. use_which_mask (str, optional): Which mask file to use (in the processed DWI space). Defaults to None.
  121. dtialps (bool, optional): Whether to perform DTI-ALPS analysis. Defaults to False.
  122. dtialps_register_method (int, optional): Registration method for DTI-ALPS. 1 for FA + FLIRT, 2 for FA + synthmorph.
  123. pved (bool, optional): Whether to perform PVeD analysis. Defaults to False.
  124. pved_method (str, optional): PVeD implementation. ``v1`` uses the legacy MATLAB/QSDR workflow; ``v2`` uses the faster Atropos workflow with MNI-space DTI maps. Defaults to ``v2``.
  125. freewater (list, optional): List of freewater estimation methods. Defaults to []. Subsets of ['single_shell_freewater', 'markvcid_freewater', 'dti_freewater'].
  126. psmd (bool, optional): Whether to perform PSMD analysis. Defaults to False.
  127. psmd_exclude_seed_mask (bool, optional): Whether to exclude seed mask in PSMD calculation. Defaults to False.
  128. visual_pathway_analysis (bool, optional): Whether to perform visual pathway analysis. Defaults to False. (An ongoing project! Experimental feature.)
  129. visual_pathway_analysis_flip_y (bool, optional): Whether to fix header (flip y) in visual pathway analysis. Defaults to False. (An ongoing project! Experimental feature.)
  130. calculate_dwi_metrics (bool, optional): Whether to calculate DWI scalar maps. Defaults to False.
  131. use_synthseg_wm_mask (bool, optional): Whether to use SynthSeg WM mask (instead of Freesurfer results). Defaults to False.
  132. exclude_seed_mask (bool, optional): Whether to exclude seed mask in DWI metrics calculation (For NAWM). Defaults to True.
  133. exclude_wmh_mask (bool, optional): Whether to exclude WMH mask in DWI metrics calculation (For NAWM). Defaults to False.
  134. keep_preproc_intermediate (bool, optional): Whether to keep the preproc_intermediate directory after preprocessing. Defaults to True. Set False to clean up after the final DWI, gradients, and brain mask are ready; workflow caches and QSIPrep outputs are not removed.
  135. extract_from (str, optional): Folder name to extract results from.
  136. """
  137. self.subject = subject
  138. self.session = session
  139. self.output_path = os.path.abspath(output_path)
  140. self.use_which_dwi = use_which_dwi # which dwi file to use
  141. self.use_which_t1w = use_which_t1w # which t1w file to use
  142. self.use_which_flair = use_which_flair # which flair file to use
  143. self.use_freesurfer_longitudinal = use_freesurfer_longitudinal # whether to use FreeSurfer longitudinal processing outputs
  144. # Preprocessing
  145. self.preprocess = preprocess # whether to preprocess
  146. self.output_resolution = output_resolution # output resolution in mm
  147. self.degibbs = degibbs # whether to perform degibbsing
  148. self.flip_b_table_axis = flip_b_table_axis # axes to flip the b-table
  149. self.preprocess_method = preprocess_method # preprocessing method
  150. self.synb0 = synb0 # whether to use synthetic b0 image for TOPUP
  151. self.use_which_reverse_b0 = use_which_reverse_b0 # whether to use reverse b0 image for TOPUP. if synb0 is True, this will be ignored
  152. # Renconstruction
  153. self.dti_fit = dti_fit # whether to fit DTI model using FSL FDT
  154. self.dti_low_b = dti_low_b # whether to only use b<=1500 volumes for DTI fitting
  155. self.dwi_t1w_register = dwi_t1w_register # whether to register DWI to T1w space
  156. self.dki_fit = dki_fit # whether to fit DKI model using Dipy
  157. self.dsistudio_dti = dsistudio_dti # whether to fit DTI model using DSI Studio
  158. self.dsistudio_gqi = dsistudio_gqi # whether to fit GQI model using DSI Studio
  159. self.dsistudio_qsdr = dsistudio_qsdr # whether to fit QSDR model using DSI Studio
  160. self.amico_noddi = amico_noddi # whether to fit NODDI model using AMICO
  161. # Global connectome and connectivity analysis
  162. self.connectome = connectome # whether to create a connectome
  163. # ROI-based tractography
  164. self.tractography = tractography
  165. self.seed_mask = seed_mask # seed mask for probabilistic tractography
  166. self.use_which_mask = use_which_mask # which mask to use (in T1w space)
  167. # DTI-ALPS
  168. self.dtialps = dtialps
  169. self.dtialps_register_method = dtialps_register_method # 1 for FA + FLIRT, 2 for T1 + FNIRT, 3 for T1 + Synthmorph
  170. # PVeD
  171. self.pved = pved
  172. self.pved_method = str(pved_method).lower()
  173. if self.pved_method not in {'v1', 'v2'}:
  174. raise ValueError(f"Unsupported PVeD method: {pved_method}. Expected 'v1' or 'v2'.")
  175. # Freewater
  176. self.freewater = freewater
  177. # PSMD
  178. self.psmd = psmd
  179. self.psmd_exclude_seed_mask = psmd_exclude_seed_mask
  180. # Visual Pathway Analysis (an ongoing project!)
  181. self.visual_pathway_analysis = visual_pathway_analysis
  182. self.visual_pathway_analysis_flip_y = visual_pathway_analysis_flip_y
  183. # Calculate scalar maps
  184. self.calculate_dwi_metrics = calculate_dwi_metrics # whether to calculate DWI metrics
  185. self.use_synthseg_wm_mask = use_synthseg_wm_mask # whether to use SynthSeg WM mask
  186. self.exclude_seed_mask = exclude_seed_mask # whether to exclude seed mask in DWI metrics calculation
  187. self.exclude_wmh_mask = exclude_wmh_mask # whether to exclude WMH mask in DWI metrics calculation
  188. self.keep_preproc_intermediate = keep_preproc_intermediate # whether to keep preproc_intermediate directory
  189. # extract results
  190. self.extract_from = extract_from
  191. #### Not user configurable ####
  192. base_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
  193. self.script_path2 = os.path.join(base_dir, 'bash', 'fdt_fs_processing.sh') # FreeSurfer processing
  194. self.psmd_skeleton_mask = os.path.join(base_dir, 'external', 'psmd', 'skeleton_mask_2019.nii.gz')
  195. def check_data_requirements(self):
  196. """
  197. check if the required data is available
  198. :return: bool
  199. """
  200. return self.session.get_dwi_files() is not None and self.session.get_t1w_files() is not None
  201. def create_workflow(self):
  202. # Find Image Files
  203. os.makedirs(self.output_path, exist_ok=True)
  204. dwi_files = self.session.get_dwi_files()
  205. phase_encoding_direction_dict = {
  206. "i": "1 0 0",
  207. "i-": "-1 0 0",
  208. "j": "0 1 0",
  209. "j-": "0 -1 0",
  210. "k": "0 0 1",
  211. "k-": "0 0 -1"
  212. }
  213. if self.use_which_dwi:
  214. dwi_files = [file for file in dwi_files if self.use_which_dwi in file['path']]
  215. if len(dwi_files) != 4:
  216. raise FileNotFoundError(f"No specific DWI file found for {self.use_which_dwi} or more than one found or not in standard BIDS format.")
  217. else:
  218. print("No specific DWI file selected. Using the first complete set.")
  219. grouped_files = {}
  220. for file in dwi_files:
  221. base_name = file['path'].split('_acq-')[0]
  222. if base_name not in grouped_files:
  223. grouped_files[base_name] = []
  224. grouped_files[base_name].append(file)
  225. first_complete_group = None
  226. for group in grouped_files.values():
  227. file_types = {file['type'] for file in group}
  228. if {'dwi', 'bval', 'bvec', 'json'}.issubset(file_types):
  229. first_complete_group = group
  230. break
  231. if not first_complete_group:
  232. raise FileNotFoundError("No complete set of DWI files found.")
  233. dwi_files = first_complete_group
  234. dwi_files_dict = {file['type']: file['path'] for file in dwi_files}
  235. print(f"[DWI Pipeline] Using DWI files: {dwi_files_dict}")
  236. dwi_image = dwi_files_dict['dwi']
  237. dwi_json = dwi_files_dict['json']
  238. dwi_bval = dwi_files_dict['bval']
  239. dwi_bvec = dwi_files_dict['bvec']
  240. # Get the total readout time and phase encoding direction from the json file
  241. with open(dwi_json, 'r') as f:
  242. dwi_json_data = json.load(f)
  243. try:
  244. total_readout_time = dwi_json_data['TotalReadoutTime']
  245. # as a str
  246. total_readout_time = str(total_readout_time)
  247. phase_encoding_direction = dwi_json_data['PhaseEncodingDirection']
  248. except KeyError:
  249. print("[Warn] Seems TotalReadoutTime or PhaseEncodingDirection not found in the json file.")
  250. total_readout_time = '0.05'
  251. phase_encoding_direction = 'j-' # default phase encoding direction
  252. print("[Warn] Set TotalReadoutTime to 0.05. But make sure the value is the same in all acquisitions.")
  253. print("[Warn] Set PhaseEncodingDirection to j-. It might not be correct")
  254. # Find Reverse b0 Image
  255. if not self.synb0 and self.use_which_reverse_b0 is not None:
  256. # TODO: Use the 'IntendedFor' field in the json file to find the reverse b0 image
  257. fmap_files = self.session.get_fmap_files()
  258. reverse_b0_files = [file for file in fmap_files if self.use_which_reverse_b0 in file['path']]
  259. if len(reverse_b0_files) != 4:
  260. raise FileNotFoundError(f"No specific reverse b0 file found for {self.use_which_reverse_b0} or more than one found.")
  261. reverse_b0_files_dict = {file['type']: file['path'] for file in reverse_b0_files}
  262. print(f"Using reverse b0 files: {reverse_b0_files_dict}")
  263. reverse_b0_image = reverse_b0_files_dict['dwi']
  264. reverse_b0_json = reverse_b0_files_dict['json']
  265. reverse_b0_bval = reverse_b0_files_dict['bval']
  266. reverse_b0_bvec = reverse_b0_files_dict['bvec']
  267. with open(reverse_b0_json, 'r') as f:
  268. reverse_b0_json_data = json.load(f)
  269. try:
  270. reverse_b0_total_readout_time = str(reverse_b0_json_data['TotalReadoutTime'])
  271. reverse_b0_phase_encoding_direction = reverse_b0_json_data['PhaseEncodingDirection']
  272. except KeyError:
  273. print("Seems TotalReadoutTime or PhaseEncodingDirection not found in the json file.")
  274. if not reverse_b0_total_readout_time:
  275. reverse_b0_total_readout_time = '0.05'
  276. print("Set TotalReadoutTime to 0.05. But make sure the value is the same in all acquisitions.")
  277. else:
  278. reverse_b0_image = None
  279. # get T1w image
  280. t1w_files = self.session.get_t1w_files()
  281. if self.use_which_t1w:
  282. t1w_files = [f for f in t1w_files if self.use_which_t1w in f]
  283. if len(t1w_files) != 1:
  284. raise FileNotFoundError(f"No specific T1w file found for {self.use_which_t1w} or more than one found.")
  285. t1w_file = t1w_files[0]
  286. print(f"[DWI Pipeline] Using T1w: {t1w_file}")
  287. else:
  288. t1w_files = [t1w_files[0]]
  289. t1w_file = t1w_files[0]
  290. print(f"No specific T1w file selected. Using the first one: {t1w_file}.")
  291. fs_output = self.session.freesurfer_dir
  292. # if fs_output exists
  293. if fs_output is None:
  294. fs_subjects_dir = ''
  295. fs_subject_id = ''
  296. fs_output_process = False
  297. print("[DWI Pipeline] No FreeSurfer output found. Skipping related processing.")
  298. else:
  299. fs_subjects_dir = os.path.dirname(fs_output)
  300. fs_subject_id = os.path.basename(fs_output)
  301. # automatically do related anat preprocessing
  302. fs_output_process = True
  303. print(f"[DWI Pipeline] FreeSurfer output found: {fs_output}. Will do related processing.")
  304. # if use longitudinal outputs, check if the longitudinal output exists (fs_subjects_dir/ses-<session_id>.long.sub-<subject_id>)
  305. if self.use_freesurfer_longitudinal:
  306. session_ids = self.subject.sessions_id
  307. fs_long_output_dir = os.path.join(fs_subjects_dir, f'ses-{self.session.session_id}.long.sub-{self.subject.subject_id}')
  308. if not os.path.exists(fs_long_output_dir):
  309. raise FileNotFoundError(f"FreeSurfer longitudinal output not found: {fs_long_output_dir}. Please check the directory or set use_freesurfer_longitudinal to False.")
  310. else:
  311. print(f"[DWI Pipeline] FreeSurfer longitudinal output found: {fs_long_output_dir}. Related to surface metrics extraction in tractography.")
  312. # for session_id in session_ids:
  313. # fs_long_output_dir = os.path.join(fs_subjects_dir, f'ses-{session_id}.long.sub-{self.subject.subject_id}')
  314. # if not os.path.exists(fs_long_output_dir):
  315. # raise FileNotFoundError(f"FreeSurfer longitudinal output not found: {fs_long_output_dir}. Please check the directory or set use_freesurfer_longitudinal to False.")
  316. # else:
  317. # print(f"[DWI Pipeline] FreeSurfer longitudinal output found: {fs_long_output_dir}. Related to surface metrics extraction in tractography.")
  318. if self.pved and self.pved_method == 'v2':
  319. pved_t1w_to_mni_warp = os.path.join(self.subject.bids_dir, 'derivatives', 'xfm', f'sub-{self.subject.subject_id}', f'ses-{self.session.session_id}', f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_from-T1w_to-MNI152NLin6ASym_warp.nii.gz')
  320. if os.path.isfile(pved_t1w_to_mni_warp):
  321. self.dti_fit = True
  322. self.dwi_t1w_register = True
  323. else:
  324. print(f"[DWI Pipeline] PVeD method='v2' requires an existing T1w-to-MNI152NLin6Asym warp, but it was not found: {pved_t1w_to_mni_warp}")
  325. print("[DWI Pipeline] Falling back to PVeD method='v1'.")
  326. self.pved_method = 'v1'
  327. if self.pved and self.pved_method == 'v1':
  328. # must run QSDR reconstruction first
  329. self.dsistudio_qsdr = True
  330. if self.visual_pathway_analysis and self.preprocess_method != 'post_qsiprep':
  331. self.dsistudio_gqi = True
  332. if self.dwi_t1w_register or self.dtialps or self.psmd:
  333. # must run DTI fitting first
  334. self.dti_fit = True
  335. # Decide one seed mask to use
  336. lesion_mask_dir = self.session._find_output(self.seed_mask)
  337. seed_mask = [file for file in os.listdir(lesion_mask_dir) if f'{self.use_which_mask}' in file]
  338. if seed_mask is not None and len(seed_mask) == 1:
  339. print(f"[DWI Pipeline] Using seed mask: {seed_mask[0]}.")
  340. seed_mask = os.path.join(lesion_mask_dir, seed_mask[0])
  341. elif seed_mask is not None and len(seed_mask) > 1:
  342. #print(f"Using the first mask found: {seed_mask[0]}.")
  343. seed_mask = os.path.join(lesion_mask_dir, seed_mask[0])
  344. print(f"[DWI Pipeline] Using seed mask: {seed_mask}.")
  345. else:
  346. seed_mask = ''
  347. print("[DWI Pipeline] No seed mask found.")
  348. #raise FileNotFoundError(f"No seed mask found in {lesion_mask_dir}. Please check the directory.")
  349. # Whether need mrtrix3 preproc dwi (nifti to mif)
  350. if 'mrtrix3' in self.connectome or 'mrtrix3' in self.tractography:
  351. mrtrix3_preproc = True
  352. else:
  353. mrtrix3_preproc = False
  354. # Whether need fdt bedpostx
  355. if 'fdt' in self.connectome or 'fdt' in self.tractography:
  356. fdt_bedpostx = True
  357. else:
  358. fdt_bedpostx = False
  359. t1w_to_mni_warp = os.path.join(self.subject.bids_dir, 'derivatives', 'xfm', f'sub-{self.subject.subject_id}', f'ses-{self.session.session_id}', f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_from-T1w_to-MNI152NLin6ASym_warp.nii.gz')
  360. if not os.path.exists(t1w_to_mni_warp):
  361. t1w_to_mni_warp = ""
  362. mni_to_t1w_warp = os.path.join(self.subject.bids_dir, 'derivatives', 'xfm', f'sub-{self.subject.subject_id}', f'ses-{self.session.session_id}', f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_from-MNI152NLin6ASym_to-T1w_warp.nii.gz')
  363. if not os.path.exists(mni_to_t1w_warp):
  364. mni_to_t1w_warp = ""
  365. # Whether need to exclude WMH mask in DWI metrics calculation
  366. wmh_mask_file = None
  367. wmh_probmap_file = None
  368. if self.exclude_wmh_mask or self.visual_pathway_analysis:
  369. wmh_output_dir = self.session._find_output('wmh_quantification')
  370. if wmh_output_dir is None:
  371. raise FileNotFoundError("WMH quantification output not found. Cannot exclude WMH mask.")
  372. # sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-T1w_label-WMH_desc-*_mask.nii.gz
  373. wmh_mask_files = [file for file in os.listdir(wmh_output_dir) if f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-T1w_label-WMH" in file and file.endswith('_mask.nii.gz')]
  374. if len(wmh_mask_files) == 1:
  375. print(f"[DWI Pipeline] Exclude WMH mask: {wmh_mask_files[0]}.")
  376. wmh_mask_file = os.path.join(wmh_output_dir, wmh_mask_files[0])
  377. # sub-XXX_ses-XXX_space-T1w_label-WMH_desc-XXX_probmap.nii.gz
  378. wmh_probmap_files = [file for file in os.listdir(wmh_output_dir) if f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-T1w_label-WMH_desc-" in file and file.endswith('_probmap.nii.gz')]
  379. if len(wmh_probmap_files) == 1:
  380. print(f"[DWI Pipeline] WMH probmap: {wmh_probmap_files[0]}.")
  381. wmh_probmap_file = os.path.join(wmh_output_dir, wmh_probmap_files[0])
  382. else:
  383. raise FileNotFoundError("[DWI Pipeline] WMH probmap in T1w space not found or multiple files found. Cannot exclude WMH probmap.")
  384. elif len(wmh_mask_files) == 2:
  385. mask_excluded_files = [file for file in wmh_mask_files if "MaskExcluded" in file]
  386. if len(mask_excluded_files) == 1:
  387. print(f"[DWI Pipeline] Exclude WMH mask: {mask_excluded_files[0]}.")
  388. wmh_mask_file = os.path.join(wmh_output_dir, mask_excluded_files[0])
  389. else:
  390. raise FileNotFoundError("[DWI Pipeline] WMH MaskExcluded file not found or multiple files found. Cannot exclude WMH mask.")
  391. wmh_probmap_files = [
  392. file for file in os.listdir(wmh_output_dir)
  393. if f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-T1w_label-WMH_desc-" in file
  394. and file.endswith("_probmap.nii.gz")
  395. ]
  396. if len(wmh_probmap_files) == 1:
  397. print(f"[DWI Pipeline] WMH probmap: {wmh_probmap_files[0]}.")
  398. wmh_probmap_file = os.path.join(wmh_output_dir, wmh_probmap_files[0])
  399. else:
  400. raise FileNotFoundError("[DWI Pipeline] WMH probmap in T1w space not found or multiple files found. Cannot exclude WMH probmap.")
  401. else:
  402. raise FileNotFoundError("[DWI Pipeline] WMH mask in T1w space not found. Cannot exclude WMH mask.")
  403. # Check chpseg segmentation output
  404. if self.session.anat_seg_dir is not None:
  405. # ./chpseg/T1w_chp_mask.nii.gz
  406. chpseg_mask_file = os.path.join(self.session.anat_seg_dir, 'chpseg', 'T1w', 'T1w_chp_mask.nii.gz')
  407. if os.path.exists(chpseg_mask_file) and self.dwi_t1w_register:
  408. print(f"[DWI Pipeline] CHP segmentation mask found: {chpseg_mask_file}. Will do related processing.")
  409. chpseg_process = True
  410. else:
  411. chpseg_process = False
  412. else:
  413. chpseg_process = False
  414. # Space entity ('preprocdwi' for our custom pipeline, 'ACPC' for QSIPrep)
  415. if self.preprocess_method == 'post_qsiprep':
  416. space_entity = 'ACPC'
  417. else:
  418. space_entity = 'preprocdwi'
  419. # ============================================
  420. # Main Workflow
  421. # ============================================
  422. dwi_workflow = Workflow(name='dwi_workflow')
  423. inputnode = Node(IdentityInterface(fields=['bids_dir', 't1w_file',
  424. 'dwi_file', 'bval_file',
  425. 'bvec_file', 'json_file',
  426. 'output_path',
  427. 'phase_encoding_number', 'total_readout_time',
  428. 'fs_output', 'seed_mask_dtispace', 'seed_mask_t1wspace',
  429. 'fs_subjects_dir', 'fs_subject_id',
  430. 'fs_processing_dir', 'qsiprep_out_bvec',
  431. 'mni_to_t1w_warp', 't1w_to_mni_warp']), name="inputnode")
  432. inputnode.inputs.bids_dir = self.subject.bids_dir
  433. inputnode.inputs.t1w_file = t1w_file
  434. inputnode.inputs.dwi_file = dwi_image
  435. inputnode.inputs.bval_file = dwi_bval
  436. inputnode.inputs.bvec_file = dwi_bvec
  437. inputnode.inputs.json_file = dwi_json
  438. inputnode.inputs.output_path = self.output_path
  439. inputnode.inputs.phase_encoding_number = phase_encoding_direction_dict[phase_encoding_direction]
  440. inputnode.inputs.total_readout_time = total_readout_time
  441. inputnode.inputs.fs_output = fs_output
  442. inputnode.inputs.fs_subjects_dir = fs_subjects_dir
  443. inputnode.inputs.fs_subject_id = fs_subject_id
  444. inputnode.inputs.seed_mask_t1wspace = seed_mask
  445. inputnode.inputs.mni_to_t1w_warp = mni_to_t1w_warp
  446. inputnode.inputs.t1w_to_mni_warp = t1w_to_mni_warp
  447. # ============================================
  448. # Preprocessing (DWI)
  449. # ============================================
  450. # Node to store the b0_all and acqparam for TOPUP
  451. b0_all_node = Node(IdentityInterface(fields=['b0_all', 'acqparam']), name='b0_all')
  452. # Node to store preprocessed DWI data
  453. preproc_dwi_node = Node(IdentityInterface(fields=['preproc_dwi', 'bvec', 'bval', 'dwi_mask',
  454. 'b0', 'preproc_t1w']) # These two used for QSIPrep output
  455. , name='preproc_dwi')
  456. if self.preprocess:
  457. # Directories for intermediate results
  458. preproc_intermediate_dir = os.path.join(self.output_path, 'preproc_intermediate')
  459. os.makedirs(preproc_intermediate_dir, exist_ok=True)
  460. print("[DWI Pipeline] Output resolution set to {} mm.".format(self.output_resolution))
  461. print("[DWI Pipeline] Flip b-table axis: {}.".format(self.flip_b_table_axis))
  462. # Padding: make sure the slices in the z direction are even (TOPUP requirement)
  463. # If no need to pad, the output will be the same as input
  464. padding_node = Node(PadDWI(), name='padding')
  465. dwi_workflow.connect(inputnode, 'dwi_file', padding_node, 'in_dwi')
  466. dwi_workflow.connect(inputnode, 'bvec_file', padding_node, 'in_bvec')
  467. dwi_workflow.connect(inputnode, 'bval_file', padding_node, 'in_bval')
  468. dwi_workflow.connect(inputnode, 'json_file', padding_node, 'in_json')
  469. padding_node.inputs.out_file = os.path.join(preproc_intermediate_dir, rename_bids_file(dwi_image, {"desc": "cropped"}, 'dwi', '.nii.gz'))
  470. # Denoise using MRtrix3
  471. denoise_node = Node(MrtrixDenoise(), name='denoise')
  472. dwi_workflow.connect(padding_node, 'out_file', denoise_node, 'dwi_img')
  473. denoise_node.inputs.output_dwi = os.path.join(preproc_intermediate_dir, 'dwi_denoised.nii.gz')
  474. # # Degibbs using MRtrix3
  475. # degibbs_node = Node(MrtrixDegibbs(), name='degibbs')
  476. # dwi_workflow.connect(denoise_node, 'output_dwi_img', degibbs_node, 'dwi_img')
  477. # dwi_workflow.connect(denoise_node, 'output_dwi_bvec', degibbs_node, 'dwi_bvec')
  478. # dwi_workflow.connect(denoise_node, 'output_dwi_bval', degibbs_node, 'dwi_bval')
  479. # degibbs_node.inputs.output_dir = preproc_intermediate_dir
  480. # Degibbs using DIPY
  481. if self.degibbs:
  482. degibbs_node = Node(DipyDegibbs(), name='degibbs')
  483. dwi_workflow.connect(denoise_node, 'output_dwi_img', degibbs_node, 'dwi_file')
  484. degibbs_node.inputs.output_dwi = os.path.join(preproc_intermediate_dir, 'dwi_degibbs.nii.gz')
  485. else:
  486. print("[DWI Pipeline] Skip degibbs.")
  487. degibbs_node = Node(IdentityInterface(fields=['output_dwi_img', 'output_dwi_bvec', 'output_dwi_bval']), name='pesudo_degibbs')
  488. dwi_workflow.connect(denoise_node, 'output_dwi_img', degibbs_node, 'output_dwi_img')
  489. dwi_workflow.connect(inputnode, 'bvec_file', degibbs_node, 'output_dwi_bvec')
  490. dwi_workflow.connect(inputnode, 'bval_file', degibbs_node, 'output_dwi_bval')
  491. # Get b0 image for TOPUP
  492. # Need a b0_all image and acqparam.txt file for FSL TOPUP
  493. if self.synb0:
  494. print("[DWI Pipeline] Use Synb0 for TOPUP.")
  495. synb0_node = Node(Synb0(), name='synb0')
  496. dwi_workflow.connect([
  497. (inputnode, synb0_node, [('t1w_file', 't1w_img'),
  498. ('json_file', 'dwi_json')]),
  499. (degibbs_node, synb0_node, [('output_dwi_img', 'dwi_img')])
  500. ])
  501. synb0_node.inputs.fmap_output_dir = os.path.join(self.subject.bids_dir, f"sub-{self.subject.subject_id}", f"ses-{self.session.session_id}", 'fmap')
  502. synb0_node.inputs.output_path_synb0 = os.path.join(self.subject.bids_dir, 'derivatives', 'synb0', f'sub-{self.subject.subject_id}', f"ses-{self.session.session_id}")
  503. dwi_workflow.connect(synb0_node, 'acqparam', b0_all_node, 'acqparam')
  504. dwi_workflow.connect(synb0_node, 'b0_all', b0_all_node, 'b0_all')
  505. elif reverse_b0_image is not None:
  506. print("[DWI Pipeline] Use reverse b0 image for TOPUP.")
  507. create_b0_acqparam_node = Node(B0AllAndAcqparam(), name='create_b0_acqparam')
  508. dwi_workflow.connect([
  509. (inputnode, create_b0_acqparam_node, [('phase_encoding_number', 'phase_encoding_number'),
  510. ('total_readout_time', 'total_readout_time')]),
  511. (degibbs_node, create_b0_acqparam_node, [('output_dwi_img', 'dwi_img'),
  512. ('output_dwi_bval', 'dwi_bval')])
  513. ])
  514. create_b0_acqparam_node.inputs.output_path = preproc_intermediate_dir
  515. create_b0_acqparam_node.inputs.reverse_dwi_img = reverse_b0_image
  516. create_b0_acqparam_node.inputs.reverse_dwi_bval = reverse_b0_bval
  517. create_b0_acqparam_node.inputs.reverse_phase_encoding_number = phase_encoding_direction_dict[reverse_b0_phase_encoding_direction]
  518. create_b0_acqparam_node.inputs.reverse_total_readout_time = reverse_b0_total_readout_time
  519. dwi_workflow.connect(create_b0_acqparam_node, 'b0_all', b0_all_node, 'b0_all')
  520. dwi_workflow.connect(create_b0_acqparam_node, 'acqparam', b0_all_node, 'acqparam')
  521. else:
  522. if self.preprocess_method != "dsistudio_mc":
  523. raise FileNotFoundError("No reverse b0 image found for TOPUP. Please check the directory or set synb0 to True to use synthetic b0 image.")
  524. print("[DWI Pipeline] No reverse b0 image available.")
  525. # Main preprocessing workflow, according to the preprocess method
  526. if self.preprocess_method == 'fdt':
  527. print("[DWI Pipeline] Use FSL FDT for DWI preprocessing.")
  528. topup_node = Node(Topup(), name='topup')
  529. dwi_workflow.connect([
  530. (b0_all_node, topup_node, [('b0_all', 'b0_all_file'),
  531. ('acqparam', 'acqparam_file')])
  532. ])
  533. topup_node.inputs.output_basename = os.path.join(preproc_intermediate_dir, 'topup_results')
  534. topup_node.inputs.output_b0_basename = os.path.join(preproc_intermediate_dir, 'hifi_b0')
  535. topup_node.inputs.config_file = 'b02b0.cnf'
  536. extract_b0_node = Node(ExtractROI(), name='extract_b0')
  537. dwi_workflow.connect(topup_node, 'b0_image', extract_b0_node, 'in_file')
  538. extract_b0_node.inputs.roi_file = os.path.join(preproc_intermediate_dir, 'dwi_b0.nii.gz')
  539. extract_b0_node.inputs.t_min = 0
  540. extract_b0_node.inputs.t_size = 1
  541. create_dwi_mask_node = Node(SynthStrip(), name='create_dwi_mask')
  542. dwi_workflow.connect([
  543. (extract_b0_node, create_dwi_mask_node, [('roi_file', 'image')])
  544. ])
  545. dwi_mask_filepath = os.path.join(preproc_intermediate_dir, 'dwi_brain_mask.nii.gz')
  546. create_dwi_mask_node.inputs.mask_file = dwi_mask_filepath
  547. create_index_node = Node(IndexTxt(), name='create_index')
  548. dwi_workflow.connect(inputnode, 'bval_file', create_index_node, 'bval_file')
  549. create_index_node.inputs.output_dir = preproc_intermediate_dir
  550. eddy_node = Node(EddyCuda(), name='eddy')
  551. dwi_workflow.connect([
  552. (degibbs_node, eddy_node, [('output_dwi_img', 'dwi_file')]),
  553. (inputnode, eddy_node, [('bval_file', 'bval_file'), ('bvec_file', 'bvec_file')]),
  554. (b0_all_node, eddy_node, [('acqparam', 'acqparam_file')]),
  555. (create_index_node, eddy_node, [('index_file', 'index_file')]),
  556. (create_dwi_mask_node, eddy_node, [('mask_file', 'mask_file')]),
  557. (topup_node, eddy_node, [('topup_basename', 'topup_basename')]),
  558. ])
  559. eddy_node.inputs.output_basename = os.path.join(preproc_intermediate_dir, 'eddy_corrected_data')
  560. order_eddy_outputs_node = Node(OrderEddyOutputs(), name='order_eddy_outputs')
  561. dwi_workflow.connect(eddy_node, 'output_filename', order_eddy_outputs_node, 'eddy_output_filename')
  562. dwi_workflow.connect(eddy_node, 'bvals', order_eddy_outputs_node, 'bval')
  563. order_eddy_outputs_node.inputs.output_resolution = self.output_resolution
  564. dwi_workflow.connect(eddy_node, 'eddy_output_dir', order_eddy_outputs_node, 'eddy_output_dir')
  565. dwi_workflow.connect(inputnode, 'output_path', order_eddy_outputs_node, 'new_output_dir')
  566. order_eddy_outputs_node.inputs.new_output_filename = rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '')
  567. flip_b_table_node = Node(FlipBTable(), name='flip_b_table')
  568. dwi_workflow.connect(order_eddy_outputs_node, 'ordered_bvec', flip_b_table_node, 'in_bvec')
  569. flip_b_table_node.inputs.flip_axis = self.flip_b_table_axis
  570. flip_b_table_node.inputs.out_bvec = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.bvec'))
  571. create_b0ref_and_brainmask_node = Node(B0RefAndBrainMask(), name='create_b0ref_and_brainmask')
  572. dwi_workflow.connect(order_eddy_outputs_node, 'ordered_dwi', create_b0ref_and_brainmask_node, 'input_dwi')
  573. dwi_workflow.connect(order_eddy_outputs_node, 'ordered_bval', create_b0ref_and_brainmask_node, 'input_bval')
  574. create_b0ref_and_brainmask_node.inputs.output_dir = self.output_path
  575. create_b0ref_and_brainmask_node.inputs.output_b0_filename = rename_bids_file(dwi_image, {"space": space_entity}, 'dwiref', '.nii.gz')
  576. create_b0ref_and_brainmask_node.inputs.output_b0_mask_filename = rename_bids_file(dwi_image, {"space": space_entity, "desc": "brain"}, 'mask', '.nii.gz')
  577. dwi_workflow.connect(order_eddy_outputs_node, 'ordered_dwi', preproc_dwi_node, 'preproc_dwi')
  578. dwi_workflow.connect(flip_b_table_node, 'out_bvec', preproc_dwi_node, 'bvec')
  579. dwi_workflow.connect(order_eddy_outputs_node, 'ordered_bval', preproc_dwi_node, 'bval')
  580. dwi_workflow.connect(create_b0ref_and_brainmask_node, 'b0_brain_mask', preproc_dwi_node, 'dwi_mask')
  581. dwi_workflow.connect(create_b0ref_and_brainmask_node, 'b0_image', preproc_dwi_node, 'b0')
  582. preproc_dwi_filename = rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.nii.gz')
  583. elif self.preprocess_method == 'mrtrix3':
  584. print("[DWI Pipeline] Use MRtrix3 for DWI preprocessing.")
  585. # Get b0_all.mif
  586. convert_b0_all_mif_node = Node(MRConvert(), name='convert_b0_all_mif')
  587. dwi_workflow.connect(b0_all_node, 'b0_all', convert_b0_all_mif_node, 'in_file')
  588. convert_b0_all_mif_node.inputs.out_file = os.path.join(preproc_intermediate_dir, 'b0_all.mif')
  589. convert_b0_all_mif_node.inputs.args = '-force'
  590. #convert_raw_dwi_mif_node.inputs.args = f'-import_pe_table {os.path.join(self.output_path, "acqparam.txt")}'
  591. # Convert denoised and degibbsed DWI to mif
  592. convert_denoised_degibbs_mif_node = Node(MRConvert(), name='convert_denoised_degibbs_mif')
  593. dwi_workflow.connect(degibbs_node, 'output_dwi_img', convert_denoised_degibbs_mif_node, 'in_file')
  594. dwi_workflow.connect(inputnode, 'bval_file', convert_denoised_degibbs_mif_node, 'in_bval')
  595. dwi_workflow.connect(inputnode, 'bvec_file', convert_denoised_degibbs_mif_node, 'in_bvec')
  596. convert_denoised_degibbs_mif_node.inputs.out_file = os.path.join(preproc_intermediate_dir, 'dwi_denoise_degibbs.mif')
  597. convert_denoised_degibbs_mif_node.inputs.args = '-force'
  598. # TOPUP and Eddy
  599. # I don't know how to integrate TOPUP and Eddy in mrtrix3 if use synb0 :(
  600. # preprocess
  601. dwi_preprocess_node = Node(DWIPreproc(), name='dwi_preprocess')
  602. dwi_workflow.connect(convert_denoised_degibbs_mif_node, 'out_file', dwi_preprocess_node, 'in_file')
  603. dwi_workflow.connect(convert_b0_all_mif_node, 'out_file', dwi_preprocess_node, 'in_epi')
  604. dwi_preprocess_node.inputs.out_file = os.path.join(preproc_intermediate_dir, 'dwi_preproc.mif')
  605. dwi_preprocess_node.inputs.rpe_options = 'pair'
  606. dwi_preprocess_node.inputs.pe_dir = phase_encoding_direction
  607. # alternative of preprocess
  608. # https://community.mrtrix.org/t/synb0-for-dwifslpreproc-how/6386/2
  609. # TODO
  610. # convert to nifti
  611. convert_preproc_dwi_node = Node(MRConvert(), name='convert_preproc_dwi')
  612. dwi_workflow.connect(dwi_preprocess_node, 'out_file', convert_preproc_dwi_node, 'in_file')
  613. convert_preproc_dwi_node.inputs.out_file = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.nii.gz'))
  614. convert_preproc_dwi_node.inputs.out_bvec = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.bvec'))
  615. convert_preproc_dwi_node.inputs.out_bval = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.bval'))
  616. convert_preproc_dwi_node.inputs.args = '-force'
  617. # get the mask
  618. extract_b0_node = Node(ExtractROI(), name='extract_b0')
  619. dwi_workflow.connect(convert_preproc_dwi_node, 'out_file', extract_b0_node, 'in_file')
  620. extract_b0_node.inputs.roi_file = os.path.join(preproc_intermediate_dir, 'dwi_b0.nii.gz')
  621. extract_b0_node.inputs.t_min = 0
  622. extract_b0_node.inputs.t_size = 1
  623. create_dwi_mask_node = Node(SynthStrip(), name='create_dwi_mask')
  624. dwi_workflow.connect(extract_b0_node, 'roi_file', create_dwi_mask_node, 'image')
  625. create_dwi_mask_node.inputs.mask_file = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity, "desc": "brain"}, 'mask', '.nii.gz'))
  626. dwi_workflow.connect(convert_preproc_dwi_node, 'out_file', preproc_dwi_node, 'preproc_dwi')
  627. dwi_workflow.connect(convert_preproc_dwi_node, 'out_bvec', preproc_dwi_node, 'bvec')
  628. dwi_workflow.connect(convert_preproc_dwi_node, 'out_bval', preproc_dwi_node, 'bval')
  629. dwi_workflow.connect(create_dwi_mask_node, 'mask_file', preproc_dwi_node, 'dwi_mask')
  630. preproc_dwi_filename = rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.nii.gz')
  631. elif self.preprocess_method == 'dsistudio_mc':
  632. print("[DWI Pipeline] Use DSI Studio for motion correction. This is designed for DSI acquisitions without reverse phase encoding.")
  633. create_src_preprocess = Node(CreateSRC(), name='create_src_preprocess')
  634. dwi_workflow.connect(degibbs_node, 'output_dwi_img', create_src_preprocess, 'source')
  635. dwi_workflow.connect(inputnode, 'bval_file', create_src_preprocess, 'bval')
  636. dwi_workflow.connect(inputnode, 'bvec_file', create_src_preprocess, 'bvec')
  637. create_src_preprocess.inputs.output = os.path.join(preproc_intermediate_dir, 'dwi_temp.src.gz')
  638. motion_correction_node = Node(DSIstudioReconstruction2(), name='dsistudio_hmc')
  639. dwi_workflow.connect(create_src_preprocess, 'out_file', motion_correction_node, 'source')
  640. motion_correction_node.inputs.thread_count = 12
  641. motion_correction_node.inputs.motion_correction = 1
  642. motion_correction_node.inputs.check_btable = 0
  643. motion_correction_node.inputs.make_isotropic = self.output_resolution
  644. motion_correction_node.inputs.save_nii = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.nii.gz'))
  645. fix_bvec_node = Node(DSIStudioBvecToFSL(), name='bvec_to_fsl')
  646. dwi_workflow.connect(motion_correction_node, 'out_bvec', fix_bvec_node, 'in_bvec')
  647. fix_bvec_node.inputs.out_bvec = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.bvec'))
  648. dwi_bvec_fix_node = Node(FixDSIStudioDWIAndBvec(), name='dwi_bvec_fix')
  649. dwi_workflow.connect(motion_correction_node, 'out_nii', dwi_bvec_fix_node, 'in_file')
  650. dwi_workflow.connect(fix_bvec_node, 'out_bvec', dwi_bvec_fix_node, 'in_bvec')
  651. dwi_bvec_fix_node.inputs.out_file = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.nii.gz'))
  652. dwi_bvec_fix_node.inputs.out_bvec = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.bvec'))
  653. dwi_workflow.connect(inputnode, 'dwi_file', dwi_bvec_fix_node, 'reference_file')
  654. dwi_resample_node = Node(ANTsResampleDWI(), name='dwi_resample_iso')
  655. dwi_workflow.connect(dwi_bvec_fix_node, 'out_file', dwi_resample_node, 'in_file')
  656. dwi_resample_node.inputs.out_file = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.nii.gz'))
  657. dwi_resample_node.inputs.voxel_size = self.output_resolution
  658. dwi_resample_node.inputs.interpolation = 4
  659. dwi_resample_node.inputs.bspline_order = 3
  660. dwi_resample_node.inputs.pixeltype = 6
  661. create_b0ref_and_brainmask_node = Node(B0RefAndBrainMask(), name='create_b0ref_and_brainmask')
  662. dwi_workflow.connect(dwi_resample_node, 'out_file', create_b0ref_and_brainmask_node, 'input_dwi')
  663. dwi_workflow.connect(motion_correction_node, 'out_bval', create_b0ref_and_brainmask_node, 'input_bval')
  664. create_b0ref_and_brainmask_node.inputs.output_dir = self.output_path
  665. create_b0ref_and_brainmask_node.inputs.output_b0_filename = rename_bids_file(dwi_image, {"space": space_entity}, 'dwiref', '.nii.gz')
  666. create_b0ref_and_brainmask_node.inputs.output_b0_mask_filename = rename_bids_file(dwi_image, {"space": space_entity, "desc": "brain"}, 'mask', '.nii.gz')
  667. dwi_workflow.connect(dwi_resample_node, 'out_file', preproc_dwi_node, 'preproc_dwi')
  668. dwi_workflow.connect(dwi_bvec_fix_node, 'out_bvec', preproc_dwi_node, 'bvec')
  669. dwi_workflow.connect(motion_correction_node, 'out_bval', preproc_dwi_node, 'bval')
  670. dwi_workflow.connect(create_b0ref_and_brainmask_node, 'b0_brain_mask', preproc_dwi_node, 'dwi_mask')
  671. dwi_workflow.connect(create_b0ref_and_brainmask_node, 'b0_image', preproc_dwi_node, 'b0')
  672. preproc_dwi_filename = rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.nii.gz')
  673. elif not self.preprocess and self.preprocess_method == 'post_qsiprep':
  674. # search for preprocessed dwi in qsiprep output
  675. qsiprep_dwi_path = os.path.join(self.subject.bids_dir, 'derivatives', 'qsiprep', f"sub-{self.subject.subject_id}", f"ses-{self.session.session_id}", 'dwi')
  676. # *space-ACPC_desc-preproc_dwi.nii.gz *space-ACPC_desc-preproc_dwi.bvec *space-ACPC_desc-preproc_dwi.bval *space-ACPC_desc-brain_mask.nii.gz
  677. target_dwi = None
  678. for file in os.listdir(qsiprep_dwi_path):
  679. if file.endswith('_desc-preproc_dwi.nii.gz'):
  680. target_dwi = os.path.join(qsiprep_dwi_path, file)
  681. break
  682. if target_dwi is None:
  683. raise FileNotFoundError(f"No preprocessed DWI file found in qsiprep output: {qsiprep_dwi_path}")
  684. target_bvec = target_dwi.replace('_dwi.nii.gz', '_dwi.bvec')
  685. target_bval = target_dwi.replace('_dwi.nii.gz', '_dwi.bval')
  686. target_mask = target_dwi.replace('_desc-preproc_dwi.nii.gz', '_desc-brain_mask.nii.gz')
  687. if not os.path.exists(target_bvec) or not os.path.exists(target_bval) or not os.path.exists(target_mask):
  688. raise FileNotFoundError(f"bvec, bval or brain mask file not found in qsiprep output: {qsiprep_dwi_path}")
  689. # search for preprocessed t1w in qsiprep output
  690. # situation 1: qsiprep/sub-XXX/anat/sub-XXX_space-ACPC_desc-preproc_T1w.nii.gz (subject-wise)
  691. # situation 2: qsiprep/sub-XXX/ses-YYY/anat/sub-XXX_ses-YYY_space-ACPC_desc-preproc_T1w.nii.gz (session-wise)
  692. qsiprep_t1w_path1 = os.path.join(self.subject.bids_dir, 'derivatives', 'qsiprep', f"sub-{self.subject.subject_id}", 'anat', f"sub-{self.subject.subject_id}_space-ACPC_desc-preproc_T1w.nii.gz")
  693. qsiprep_t1w_path2 = os.path.join(self.subject.bids_dir, 'derivatives', 'qsiprep', f"sub-{self.subject.subject_id}", f"ses-{self.session.session_id}", 'anat', f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-ACPC_desc-preproc_T1w.nii.gz")
  694. if os.path.exists(qsiprep_t1w_path1):
  695. target_t1w = qsiprep_t1w_path1
  696. elif os.path.exists(qsiprep_t1w_path2):
  697. target_t1w = qsiprep_t1w_path2
  698. else:
  699. raise FileNotFoundError(f"No preprocessed T1w file found in qsiprep output: {qsiprep_t1w_path1} or {qsiprep_t1w_path2}")
  700. # check QSIPrep outputs
  701. print(f"[DWI Pipeline] You selected to use QSIPrep preprocessed DWI and T1w images.")
  702. print(f"[DWI Pipeline] Found preprocessed DWI: {target_dwi}")
  703. print(f"[DWI Pipeline] Found preprocessed T1w: {target_t1w}")
  704. # Flip bvec
  705. print("[DWI Pipeline] Flip b-table axis: {}.".format(self.flip_b_table_axis))
  706. inputnode.inputs.qsiprep_out_bvec = target_bvec
  707. flip_b_table_node = Node(FlipBTable(), name='flip_b_table_qsiprep')
  708. dwi_workflow.connect([(inputnode, flip_b_table_node, [('qsiprep_out_bvec', 'in_bvec')])])
  709. flip_b_table_node.inputs.flip_axis = self.flip_b_table_axis
  710. flip_b_table_node.inputs.out_bvec = target_bvec
  711. preproc_dwi_node.inputs.preproc_dwi = target_dwi
  712. #preproc_dwi_node.inputs.bvec = target_bvec
  713. dwi_workflow.connect([(flip_b_table_node, preproc_dwi_node, [('out_bvec', 'bvec')])])
  714. preproc_dwi_node.inputs.bval = target_bval
  715. preproc_dwi_node.inputs.dwi_mask = target_mask
  716. preproc_dwi_node.inputs.b0 = target_dwi.replace('_desc-preproc_dwi.nii.gz', '_dwiref.nii.gz')
  717. preproc_dwi_node.inputs.preproc_t1w = target_t1w
  718. preproc_dwi_filename = rename_bids_file(dwi_image, {"space": "ACPC", "desc": "preproc"}, 'dwi', '.nii.gz')
  719. else:
  720. print("[DWI Pipeline] Assuming DWI has been preprocessed and arching to output directory.")
  721. preproc_dwi_node.inputs.preproc_dwi = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.nii.gz'))
  722. preproc_dwi_node.inputs.bvec = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.bvec'))
  723. preproc_dwi_node.inputs.bval = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.bval'))
  724. preproc_dwi_node.inputs.dwi_mask = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity, "desc": "brain"}, 'mask', '.nii.gz'))
  725. preproc_dwi_node.inputs.b0 = os.path.join(self.output_path, rename_bids_file(dwi_image, {"space": space_entity}, 'dwiref', '.nii.gz'))
  726. preproc_dwi_filename = rename_bids_file(dwi_image, {"space": space_entity, "desc": "preproc"}, 'dwi', '.nii.gz')
  727. # Wait for all final outputs, including the mask whose creation can still
  728. # read intermediate files after the final DWI has been written.
  729. if not self.keep_preproc_intermediate and self.preprocess_method != 'post_qsiprep':
  730. def _cleanup_intermediate(preproc_dwi, bvec, bval, dwi_mask, intermediate_dir):
  731. import os
  732. import shutil
  733. if os.path.isdir(intermediate_dir):
  734. target = os.path.realpath(intermediate_dir)
  735. for path in (preproc_dwi, bvec, bval, dwi_mask):
  736. if not os.path.isfile(path):
  737. raise FileNotFoundError(f"Cannot clean preprocessing intermediates: missing final output {path}")
  738. if os.path.commonpath([target, os.path.realpath(path)]) == target:
  739. raise ValueError(f"Cannot remove intermediate directory containing final output {path}")
  740. shutil.rmtree(intermediate_dir)
  741. print(f"[DWI Pipeline] Removed preprocessing intermediates: {intermediate_dir}")
  742. return preproc_dwi
  743. cleanup_fields = ['preproc_dwi', 'bvec', 'bval', 'dwi_mask']
  744. cleanup_node = Node(Function(input_names=cleanup_fields + ['intermediate_dir'],
  745. output_names=['dummy'],
  746. function=_cleanup_intermediate),
  747. name='cleanup_preproc_intermediate', always_run=True)
  748. cleanup_node.inputs.intermediate_dir = os.path.join(self.output_path, 'preproc_intermediate')
  749. for field in cleanup_fields:
  750. dwi_workflow.connect(preproc_dwi_node, field, cleanup_node, field)
  751. # MRtrix3 conversion (nifti to mif): preprocessed DWI and brain mask
  752. if mrtrix3_preproc:
  753. mrtrix3_output_dir = os.path.join(self.output_path, 'mrtrix3')
  754. os.makedirs(mrtrix3_output_dir, exist_ok=True)
  755. preproc_dwi_mif_node = Node(IdentityInterface(fields=['preproc_dwi', 'dwi_mask']), name='preproc_dwi_mif_node')
  756. preproc_dwi_mif = os.path.join(mrtrix3_output_dir, rename_bids_file(preproc_dwi_filename, {}, 'dwi', '.mif'))
  757. if not os.path.exists(preproc_dwi_mif):
  758. mri_convert_node = Node(MRConvert(), name='dwi_mif_convert')
  759. mri_convert_node.inputs.out_file = preproc_dwi_mif
  760. mri_convert_node.inputs.args = '-force'
  761. dwi_workflow.connect(preproc_dwi_node, 'preproc_dwi', mri_convert_node, 'in_file')
  762. dwi_workflow.connect(preproc_dwi_node, 'bvec', mri_convert_node, 'in_bvec')
  763. dwi_workflow.connect(preproc_dwi_node, 'bval', mri_convert_node, 'in_bval')
  764. dwi_workflow.connect(mri_convert_node, 'out_file', preproc_dwi_mif_node, 'preproc_dwi')
  765. else:
  766. preproc_dwi_mif_node.inputs.preproc_dwi = preproc_dwi_mif
  767. brain_mask_mif = os.path.join(mrtrix3_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': 'brain'}, 'mask', '.mif'))
  768. if not os.path.exists(brain_mask_mif):
  769. mask_mif_convert = Node(MRConvert(), name='mask_mif_convert')
  770. mask_mif_convert.inputs.out_file = brain_mask_mif
  771. mask_mif_convert.inputs.args = '-force'
  772. dwi_workflow.connect(preproc_dwi_node, 'dwi_mask', mask_mif_convert, 'in_file')
  773. dwi_workflow.connect(mask_mif_convert, 'out_file', preproc_dwi_mif_node, 'dwi_mask')
  774. else:
  775. preproc_dwi_mif_node.inputs.dwi_mask = brain_mask_mif
  776. # ===========================================
  777. # DWI scalar maps storage node
  778. dwi_scalarmaps_output_node = Node(IdentityInterface(fields=[
  779. 'fa_img', 'md_img', 'markvcid2_fw_img', 'ad', 'rd', # Tensor model
  780. 'odi_img', 'icvf_img', 'isovf_img', # NODDI model
  781. 'gqi_gfa_img', 'gqi_iso_img', 'gqi_qa_img', # GQI model
  782. 'chidia_img', # QSM!
  783. 'mk_img', 'ak_img', 'rk_img', # DKI model
  784. ]), name='dwi_scalarmaps_output')
  785. # default all to ''
  786. dwi_scalarmaps_output_node.inputs.fa_img = ''
  787. dwi_scalarmaps_output_node.inputs.md_img = ''
  788. dwi_scalarmaps_output_node.inputs.markvcid2_fw_img = ''
  789. dwi_scalarmaps_output_node.inputs.ad = ''
  790. dwi_scalarmaps_output_node.inputs.rd = ''
  791. dwi_scalarmaps_output_node.inputs.odi_img = ''
  792. dwi_scalarmaps_output_node.inputs.icvf_img = ''
  793. dwi_scalarmaps_output_node.inputs.isovf_img = ''
  794. dwi_scalarmaps_output_node.inputs.gqi_gfa_img = ''
  795. dwi_scalarmaps_output_node.inputs.gqi_iso_img = ''
  796. dwi_scalarmaps_output_node.inputs.gqi_qa_img = ''
  797. dwi_scalarmaps_output_node.inputs.chidia_img = ''
  798. dwi_scalarmaps_output_node.inputs.mk_img = ''
  799. dwi_scalarmaps_output_node.inputs.ak_img = ''
  800. dwi_scalarmaps_output_node.inputs.rk_img = ''
  801. dwi_scalarmaps_output_node.inputs.chidia_img = ''
  802. # ===========================================
  803. # ===========================================
  804. # Reconstruction (FSL FDT dtifit)
  805. # ===========================================
  806. dti_fit_output_node = Node(
  807. IdentityInterface(fields=['fa_img', 'md_img', 'tensor_img']),
  808. name='dti_fit_output'
  809. )
  810. if self.dti_fit:
  811. print("[DWI Pipeline] DTI fitting using FSL dtifit: True")
  812. dtifit_output_dir = os.path.join(self.output_path, 'dtifit')
  813. os.makedirs(dtifit_output_dir, exist_ok=True)
  814. target_fa_output = os.path.join(dtifit_output_dir, rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "tensor", "param": "fa"}, "dwimap", ".nii.gz"))
  815. target_ad_output = os.path.join(dtifit_output_dir, rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "tensor", "param": "ad"}, "dwimap", ".nii.gz"))
  816. if os.path.exists(target_ad_output):
  817. print(f"[DWI Pipeline] DTI AD image already exists: {target_ad_output}. Skipping DTI fit.")
  818. dti_fit_output_node.inputs.fa_img = target_fa_output
  819. dti_fit_output_node.inputs.md_img = os.path.join(dtifit_output_dir,rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "tensor", "param": "md"}, "dwimap", ".nii.gz"))
  820. dti_fit_output_node.inputs.ad_img = os.path.join(dtifit_output_dir, rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "tensor", "param": "ad"}, "dwimap", ".nii.gz"))
  821. dti_fit_output_node.inputs.rd_img = os.path.join(dtifit_output_dir, rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "tensor", "param": "rd"}, "dwimap", ".nii.gz"))
  822. dti_fit_output_node.inputs.tensor_img = os.path.join(dtifit_output_dir, rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "tensor", "param": "tensor"}, "dwimap", ".nii.gz"))
  823. else:
  824. dti_fit_node = Node(DTIFitBIDS(), name='dti_fit')
  825. dwi_workflow.connect(preproc_dwi_node, 'preproc_dwi', dti_fit_node, 'dwi_file')
  826. dwi_workflow.connect(preproc_dwi_node, 'bvec', dti_fit_node, 'bvec_file')
  827. dwi_workflow.connect(preproc_dwi_node, 'bval', dti_fit_node, 'bval_file')
  828. dwi_workflow.connect(preproc_dwi_node, 'dwi_mask', dti_fit_node, 'mask_file')
  829. dti_fit_node.inputs.output_basename = os.path.join(dtifit_output_dir, 'dti')
  830. dti_fit_node.inputs.bids_rename = True
  831. dti_fit_node.inputs.overwrite = True
  832. dti_fit_node.inputs.dti_low_b = self.dti_low_b
  833. dwi_workflow.connect(dti_fit_node, 'dti_fa', dti_fit_output_node, 'fa_img')
  834. dwi_workflow.connect(dti_fit_node, 'dti_md', dti_fit_output_node, 'md_img')
  835. dwi_workflow.connect(dti_fit_node, 'dti_ad', dti_fit_output_node, 'ad_img')
  836. dwi_workflow.connect(dti_fit_node, 'dti_rd', dti_fit_output_node, 'rd_img')
  837. dwi_workflow.connect(dti_fit_node, 'dti_tensor', dti_fit_output_node, 'tensor_img')
  838. # Connect DTI outputs to dwi_scalarmaps_output_node
  839. dwi_workflow.connect(dti_fit_output_node, 'fa_img', dwi_scalarmaps_output_node, 'fa_img')
  840. dwi_workflow.connect(dti_fit_output_node, 'md_img', dwi_scalarmaps_output_node, 'md_img')
  841. dwi_workflow.connect(dti_fit_output_node, 'ad_img', dwi_scalarmaps_output_node, 'ad_img')
  842. dwi_workflow.connect(dti_fit_output_node, 'rd_img', dwi_scalarmaps_output_node, 'rd_img')
  843. # ===========================================
  844. # Reconstruction (Dipy DKI fit)
  845. # ===========================================
  846. dki_fit_output_node = Node(
  847. IdentityInterface(fields=['mk_img', 'ak_img', 'rk_img']),
  848. name='dki_fit_output'
  849. )
  850. if self.dki_fit:
  851. print("[DWI Pipeline] DKI fitting using Dipy: True")
  852. dki_output_dir = os.path.join(self.output_path, 'dkifit')
  853. os.makedirs(dki_output_dir, exist_ok=True)
  854. target_mk_output = os.path.join(dki_output_dir, rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "dki", "param": "mk"}, "dwimap", ".nii.gz"))
  855. if os.path.exists(target_mk_output):
  856. print(f"[DWI Pipeline] DKI MK image already exists: {target_mk_output}. Skipping DKI fit.")
  857. dki_fit_output_node.inputs.mk_img = target_mk_output
  858. dki_fit_output_node.inputs.ak_img = os.path.join(dki_output_dir, rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "dki", "param": "ak"}, "dwimap", ".nii.gz"))
  859. dki_fit_output_node.inputs.rk_img = os.path.join(dki_output_dir, rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "dki", "param": "rk"}, "dwimap", ".nii.gz"))
  860. else:
  861. dki_fit_node = Node(DKIFit(), name='dki_fit')
  862. dwi_workflow.connect(preproc_dwi_node, 'preproc_dwi', dki_fit_node, 'dwi_file')
  863. dwi_workflow.connect(preproc_dwi_node, 'bvec', dki_fit_node, 'bvec_file')
  864. dwi_workflow.connect(preproc_dwi_node, 'bval', dki_fit_node, 'bval_file')
  865. dwi_workflow.connect(preproc_dwi_node, 'dwi_mask', dki_fit_node, 'mask_file')
  866. dki_fit_node.inputs.output_dir = dki_output_dir
  867. dki_fit_node.inputs.bids_rename = True
  868. dki_fit_node.inputs.overwrite = True
  869. dki_fit_node.inputs.num_threads = 8
  870. dwi_workflow.connect(dki_fit_node, 'mk_file', dki_fit_output_node, 'mk_img')
  871. dwi_workflow.connect(dki_fit_node, 'ak_file', dki_fit_output_node, 'ak_img')
  872. dwi_workflow.connect(dki_fit_node, 'rk_file', dki_fit_output_node, 'rk_img')
  873. # Connect DKI outputs to dwi_scalarmaps_output_node
  874. dwi_workflow.connect(dki_fit_output_node, 'mk_img', dwi_scalarmaps_output_node, 'mk_img')
  875. dwi_workflow.connect(dki_fit_output_node, 'ak_img', dwi_scalarmaps_output_node, 'ak_img')
  876. dwi_workflow.connect(dki_fit_output_node, 'rk_img', dwi_scalarmaps_output_node, 'rk_img')
  877. # ============================================
  878. # DWI and T1w registration (And other anat preprocessing)
  879. # ============================================
  880. # We put it here because we prefer to use FA image for registration (FLIRT), while QSIPrep uses b0 image !
  881. nawm_dseg = None
  882. if self.dwi_t1w_register:
  883. print("[DWI Pipeline] Register DWI to T1w image: True")
  884. anat_output_dir = os.path.join(self.output_path, 'anat')
  885. os.makedirs(anat_output_dir, exist_ok=True)
  886. xfm_dir = os.path.join(self.subject.bids_dir, 'derivatives', 'xfm', f"sub-{self.subject.subject_id}", f"ses-{self.session.session_id}")
  887. if self.session.xfm_dir is None:
  888. os.makedirs(os.path.join(self.subject.bids_dir, 'derivatives', 'xfm', f"sub-{self.subject.subject_id}", f"ses-{self.session.session_id}"), exist_ok=True)
  889. if self.preprocess_method == 'post_qsiprep':
  890. qsiprep_out_matrix = os.path.join(xfm_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_from-ACPC_to-T1w_xfm.mat")
  891. if os.path.exists(qsiprep_out_matrix):
  892. print(f"[DWI Pipeline] QSIPrep transform matrix already exists: {qsiprep_out_matrix}. Skipping.")
  893. dwi_to_t1w_reg_node = Node(IdentityInterface(fields=['out_matrix_file', 'out_file']), name='dwi_to_t1w_reg_qsiprep')
  894. dwi_to_t1w_reg_node.inputs.out_matrix_file = qsiprep_out_matrix
  895. dwi_to_t1w_reg_node.inputs.out_file = ''
  896. else:
  897. dwi_to_t1w_reg_node = Node(QSIPrepOrigToACPC(), name='dwi_to_t1w_reg_qsiprep')
  898. dwi_workflow.connect(inputnode, 't1w_file', dwi_to_t1w_reg_node, 'orig_t1w_file')
  899. dwi_workflow.connect(dti_fit_output_node, 'fa_img', dwi_to_t1w_reg_node, 'dwimap_file')
  900. dwi_to_t1w_reg_node.inputs.subject_id = self.subject.subject_id
  901. dwi_to_t1w_reg_node.inputs.session_id = self.session.session_id
  902. dwi_to_t1w_reg_node.inputs.output_dir = xfm_dir
  903. dwi_workflow.connect(preproc_dwi_node, 'preproc_t1w', dwi_to_t1w_reg_node, 'preprocess_t1w_file')
  904. else:
  905. dwi_to_t1w_reg_node = Node(FLIRT(), name='dwi_to_t1w_reg')
  906. t1w_stripped = os.path.join(xfm_dir, rename_bids_file(t1w_file, {'desc': 'brain', 'space': None}, 'T1w', '.nii.gz'))
  907. dwi_workflow.connect(dti_fit_output_node, 'fa_img', dwi_to_t1w_reg_node, 'in_file')
  908. if os.path.exists(t1w_stripped):
  909. dwi_to_t1w_reg_node.inputs.reference = t1w_stripped
  910. else:
  911. t1w_strip_node = Node(SynthStrip(), name='t1w_strip_for_dwi_reg')
  912. dwi_workflow.connect(inputnode, 't1w_file', t1w_strip_node, 'image')
  913. t1w_strip_node.inputs.out_file = t1w_stripped
  914. dwi_workflow.connect(t1w_strip_node, 'out_file', dwi_to_t1w_reg_node, 'reference')
  915. dwi_to_t1w_reg_node.inputs.out_file = os.path.join(xfm_dir, rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "tensor", "param": "fa", 'space': 'T1w'}, 'dwimap', '.nii.gz'))
  916. dwi_to_t1w_reg_node.inputs.out_matrix_file = os.path.join(xfm_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_from-{space_entity}_to-T1w_xfm.mat")
  917. # need to invert the matrix for later use
  918. invert_dwi_to_t1w_reg_node = Node(ConvertXFM(), name='invert_dwi_to_t1w_reg')
  919. dwi_workflow.connect(dwi_to_t1w_reg_node, 'out_matrix_file', invert_dwi_to_t1w_reg_node, 'in_file')
  920. invert_dwi_to_t1w_reg_node.inputs.out_file = os.path.join(xfm_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_from-T1w_to-{space_entity}_xfm.mat")
  921. invert_dwi_to_t1w_reg_node.inputs.invert_xfm = True
  922. # If seed_mask_t1wspace is provided, we need to transform it to DWI space
  923. if seed_mask != '':
  924. seed_mask_to_dwi_node = Node(FLIRT(), name='seed_mask_to_dwi')
  925. dwi_workflow.connect(inputnode, 'seed_mask_t1wspace', seed_mask_to_dwi_node, 'in_file')
  926. dwi_workflow.connect(invert_dwi_to_t1w_reg_node, 'out_file', seed_mask_to_dwi_node, 'in_matrix_file')
  927. dwi_workflow.connect(preproc_dwi_node, 'b0', seed_mask_to_dwi_node, 'reference')
  928. seed_mask_to_dwi_node.inputs.interp = 'nearestneighbour'
  929. seed_mask_to_dwi_node.inputs.apply_xfm = True
  930. seed_mask_to_dwi_node.inputs.out_file = os.path.join(os.path.dirname(seed_mask), rename_bids_file(seed_mask, {'space': space_entity}, 'mask', '.nii.gz'))
  931. if chpseg_process:
  932. chpseg_mask_to_dwi_node = Node(FLIRT(), name='chpseg_mask_to_dwi')
  933. chpseg_mask_to_dwi_node.inputs.in_file = chpseg_mask_file
  934. dwi_workflow.connect(invert_dwi_to_t1w_reg_node, 'out_file', chpseg_mask_to_dwi_node, 'in_matrix_file')
  935. dwi_workflow.connect(preproc_dwi_node, 'b0', chpseg_mask_to_dwi_node, 'reference')
  936. chpseg_mask_to_dwi_node.inputs.interp = 'nearestneighbour'
  937. chpseg_mask_to_dwi_node.inputs.apply_xfm = True
  938. chpseg_mask_to_dwi_node.inputs.out_file = os.path.join(anat_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-{space_entity}_desc-chpseg_dseg.nii.gz")
  939. if fs_output_process:
  940. # Freesurfer to T1w registration
  941. fs_to_t1w_xfm_node = Node(Tkregister2fs2t1w(), name='fs_to_t1w_xfm')
  942. dwi_workflow.connect(inputnode, 'fs_subjects_dir', fs_to_t1w_xfm_node, 'fs_subjects_dir')
  943. dwi_workflow.connect(inputnode, 'fs_subject_id', fs_to_t1w_xfm_node, 'fs_subject_id')
  944. fs_to_t1w_xfm_node.inputs.output_matrix = os.path.join(xfm_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_from-fs_to-T1w_xfm.mat")
  945. fs_to_t1w_xfm_node.inputs.output_inverse_matrix = os.path.join(xfm_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_from-T1w_to-fs_xfm.mat")
  946. concat_fs_to_dwi_xfm_node = Node(ConvertXFM(), name='concat_fs_to_dwi_xfm')
  947. dwi_workflow.connect(fs_to_t1w_xfm_node, 'output_matrix', concat_fs_to_dwi_xfm_node, 'in_file')
  948. concat_fs_to_dwi_xfm_node.inputs.concat_xfm = True
  949. dwi_workflow.connect(invert_dwi_to_t1w_reg_node, 'out_file', concat_fs_to_dwi_xfm_node, 'in_file2')
  950. concat_fs_to_dwi_xfm_node.inputs.out_file = os.path.join(xfm_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_from-fs_to-{space_entity}_xfm.mat")
  951. # Seed mask to fs space registration
  952. if seed_mask != '':
  953. seed_mask_to_fs_node = Node(MRIVol2Vol(), name='seed_mask_to_fs')
  954. seed_mask_to_fs_node.inputs.moving_image = seed_mask # in T1w space
  955. dwi_workflow.connect(fs_to_t1w_xfm_node, 'output_inverse_matrix', seed_mask_to_fs_node, 'fsl_matrix')
  956. seed_mask_to_fs_node.inputs.target_image = os.path.join(fs_output, 'mri', 'orig.mgz') # fs native space
  957. seed_mask_to_fs_node.inputs.interp = 'nearestneighbour'
  958. seed_mask_to_fs_node.inputs.output_image = os.path.join(os.path.dirname(seed_mask), rename_bids_file(seed_mask, {'space': 'fs'}, 'mask', '.nii.gz'))
  959. # Refine freesurfer aparc+aseg parcellation
  960. def get_aseg_file(fs_subjects_dir, fs_subject_id):
  961. import os
  962. return os.path.join(fs_subjects_dir, fs_subject_id, 'mri', 'aparc+aseg.mgz')
  963. get_aseg_node = Node(Function(input_names=['fs_subjects_dir', 'fs_subject_id'],
  964. output_names=['aseg_file'],
  965. function=get_aseg_file), name='get_aseg_file')
  966. dwi_workflow.connect(inputnode, 'fs_subjects_dir', get_aseg_node, 'fs_subjects_dir')
  967. dwi_workflow.connect(inputnode, 'fs_subject_id', get_aseg_node, 'fs_subject_id')
  968. aseg_mgz_to_nifti_node = Node(MRIConvert(), name='aseg_mgz_to_nifti')
  969. dwi_workflow.connect(get_aseg_node, 'aseg_file', aseg_mgz_to_nifti_node, 'in_file')
  970. aseg_mgz_to_nifti_node.inputs.out_file = os.path.join(anat_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-fs_aparcaseg.nii.gz")
  971. # convert to DWI space
  972. fs_aparcaseg_to_dwi_node = Node(FLIRT(), name='fs_aparcaseg_to_dwi')
  973. dwi_workflow.connect(aseg_mgz_to_nifti_node, 'out_file', fs_aparcaseg_to_dwi_node, 'in_file')
  974. dwi_workflow.connect(concat_fs_to_dwi_xfm_node, 'out_file', fs_aparcaseg_to_dwi_node, 'in_matrix_file')
  975. dwi_workflow.connect(preproc_dwi_node, 'b0', fs_aparcaseg_to_dwi_node, 'reference')
  976. fs_aparcaseg_to_dwi_node.inputs.interp = 'nearestneighbour'
  977. fs_aparcaseg_to_dwi_node.inputs.apply_xfm = True
  978. fs_aparcaseg_to_dwi_node.inputs.out_file = os.path.join(anat_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-{space_entity}_aparcaseg.nii.gz")
  979. # Remove WMH from aseg
  980. refine_aseg_node = Node(RemoveWMH(), name='refine_aseg')
  981. dwi_workflow.connect(fs_aparcaseg_to_dwi_node, 'out_file', refine_aseg_node, 'in_aseg')
  982. refine_aseg_node.inputs.out_aseg = os.path.join(anat_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-{space_entity}_desc-WMHremoved_aparcaseg.nii.gz")
  983. create_cortical_GM_mask_node = Node(MergeRibbon(), name='create_cortical_GM_mask')
  984. dwi_workflow.connect(inputnode, 'fs_subjects_dir', create_cortical_GM_mask_node, 'subjects_dir')
  985. dwi_workflow.connect(inputnode, 'fs_subject_id', create_cortical_GM_mask_node, 'subject_id')
  986. create_cortical_GM_mask_node.inputs.output_gm_mask = os.path.join(anat_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-fs_label-corticalGM_mask.nii.gz")
  987. cortical_GM_to_dwi_node = Node(FLIRT(), name='cortical_GM_to_dwi')
  988. dwi_workflow.connect(create_cortical_GM_mask_node, 'output_gm_mask', cortical_GM_to_dwi_node, 'in_file')
  989. dwi_workflow.connect(concat_fs_to_dwi_xfm_node, 'out_file', cortical_GM_to_dwi_node, 'in_matrix_file')
  990. dwi_workflow.connect(preproc_dwi_node, 'b0', cortical_GM_to_dwi_node, 'reference')
  991. cortical_GM_to_dwi_node.inputs.interp = "nearestneighbour"
  992. cortical_GM_to_dwi_node.inputs.apply_xfm = True
  993. cortical_GM_to_dwi_node.inputs.out_file = os.path.join(
  994. anat_output_dir, rename_bids_file(create_cortical_GM_mask_node.inputs.output_gm_mask, {"space": space_entity}, "mask", ".nii.gz")
  995. )
  996. # convert wmparc.mgz to DWI space
  997. def get_wmparc_file(fs_subjects_dir, fs_subject_id):
  998. import os
  999. return os.path.join(fs_subjects_dir, fs_subject_id, 'mri', 'wmparc.mgz')
  1000. get_wmparc_node = Node(Function(input_names=['fs_subjects_dir', 'fs_subject_id'],
  1001. output_names=['wmparc_file'],
  1002. function=get_wmparc_file), name='get_wmparc_file')
  1003. dwi_workflow.connect(inputnode, 'fs_subjects_dir', get_wmparc_node, 'fs_subjects_dir')
  1004. dwi_workflow.connect(inputnode, 'fs_subject_id', get_wmparc_node, 'fs_subject_id')
  1005. wmparc_mgz_to_nifti_node = Node(MRIConvert(), name='wmparc_mgz_to_nifti')
  1006. dwi_workflow.connect(get_wmparc_node, 'wmparc_file', wmparc_mgz_to_nifti_node, 'in_file')
  1007. wmparc_mgz_to_nifti_node.inputs.out_file = os.path.join(anat_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-fs_wmparc.nii.gz")
  1008. # convert to DWI space
  1009. fs_wmparc_to_dwi_node = Node(FLIRT(), name='fs_wmparc_to_dwi')
  1010. dwi_workflow.connect(wmparc_mgz_to_nifti_node, 'out_file', fs_wmparc_to_dwi_node, 'in_file')
  1011. dwi_workflow.connect(concat_fs_to_dwi_xfm_node, 'out_file', fs_wmparc_to_dwi_node, 'in_matrix_file')
  1012. dwi_workflow.connect(preproc_dwi_node, 'b0', fs_wmparc_to_dwi_node, 'reference')
  1013. fs_wmparc_to_dwi_node.inputs.interp = 'nearestneighbour'
  1014. fs_wmparc_to_dwi_node.inputs.apply_xfm = True
  1015. fs_wmparc_to_dwi_node.inputs.out_file = os.path.join(anat_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-{space_entity}_wmparc.nii.gz")
  1016. # exclude seed mask from wmparc (DWI space) if exists
  1017. final_wmparc_node = Node(IdentityInterface(fields=['final_wmparc']), name='final_wmparc')
  1018. if seed_mask != '':
  1019. remove_seed_from_wmparc_node = Node(RemoveMaskRegion(), name='remove_seed_from_wmparc')
  1020. remove_seed_from_wmparc_node.inputs.output_image = os.path.join(
  1021. anat_output_dir,
  1022. f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-{space_entity}_desc-seedRemoved_wmparc.nii.gz"
  1023. )
  1024. dwi_workflow.connect(fs_wmparc_to_dwi_node, 'out_file', remove_seed_from_wmparc_node, 'input_image')
  1025. dwi_workflow.connect(seed_mask_to_dwi_node, 'out_file', remove_seed_from_wmparc_node, 'mask_image')
  1026. dwi_workflow.connect(remove_seed_from_wmparc_node, 'output_image', final_wmparc_node, 'final_wmparc')
  1027. else:
  1028. dwi_workflow.connect(fs_wmparc_to_dwi_node, 'out_file', final_wmparc_node, 'final_wmparc')
  1029. # Load CHIDIA (QSM) image if available, independent of visual_pathway_analysis
  1030. chidia_in_t1w_file = os.path.join(
  1031. self.subject.bids_dir, 'derivatives', 'qsm_pipeline', f"sub-{self.subject.subject_id}", f"ses-{self.session.session_id}", "QSM_registered",
  1032. f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-T1w_Chidia.nii.gz",
  1033. )
  1034. if os.path.exists(chidia_in_t1w_file):
  1035. chidia_to_dwi_node = Node(FLIRT(), name='chidia_to_dwi')
  1036. chidia_to_dwi_node.inputs.in_file = chidia_in_t1w_file
  1037. dwi_workflow.connect(invert_dwi_to_t1w_reg_node, 'out_file', chidia_to_dwi_node, 'in_matrix_file')
  1038. dwi_workflow.connect(preproc_dwi_node, 'b0', chidia_to_dwi_node, 'reference')
  1039. chidia_to_dwi_node.inputs.interp = 'trilinear'
  1040. chidia_to_dwi_node.inputs.apply_xfm = True
  1041. chidia_to_dwi_node.inputs.out_file = os.path.join(anat_output_dir, rename_bids_file(chidia_in_t1w_file, {'space': space_entity}, "Chidia", ".nii.gz"))
  1042. dwi_workflow.connect(chidia_to_dwi_node, 'out_file', dwi_scalarmaps_output_node, 'chidia_img')
  1043. if self.visual_pathway_analysis:
  1044. if wmh_probmap_file is not None:
  1045. wmhprobmap_to_dwi_node = Node(FLIRT(), name='wmhprobmap_to_dwi')
  1046. wmhprobmap_to_dwi_node.inputs.in_file = wmh_probmap_file
  1047. dwi_workflow.connect(invert_dwi_to_t1w_reg_node, 'out_file', wmhprobmap_to_dwi_node, 'in_matrix_file')
  1048. dwi_workflow.connect(preproc_dwi_node, 'b0', wmhprobmap_to_dwi_node, 'reference')
  1049. wmhprobmap_to_dwi_node.inputs.interp = 'trilinear'
  1050. wmhprobmap_to_dwi_node.inputs.apply_xfm = True
  1051. wmhprobmap_to_dwi_node.inputs.out_file = os.path.join(anat_output_dir, rename_bids_file(wmh_probmap_file, {'space': space_entity}, "probmap", ".nii.gz"))
  1052. # NAWM mask creation
  1053. base_name = f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-{space_entity}"
  1054. wm_source_available = False
  1055. wm_source_node = None
  1056. wm_source_output = None
  1057. if self.use_synthseg_wm_mask:
  1058. synthseg_dir = os.path.join(
  1059. self.subject.bids_dir, "derivatives", "anat_seg",
  1060. f"sub-{self.subject.subject_id}", f"ses-{self.session.session_id}", "synthseg"
  1061. )
  1062. matches = glob.glob(os.path.join(synthseg_dir, "*synthseg.nii.gz"))
  1063. if len(matches) != 1:
  1064. raise RuntimeError(f"Expected exactly 1 SynthSeg output, found {len(matches)}: {matches}")
  1065. synthseg_target_output = matches[0]
  1066. synthseg_to_dwi_node = Node(FLIRT(), name="synthseg_to_dwi")
  1067. synthseg_to_dwi_node.inputs.in_file = synthseg_target_output
  1068. synthseg_to_dwi_node.inputs.interp = "nearestneighbour"
  1069. synthseg_to_dwi_node.inputs.apply_xfm = True
  1070. synthseg_to_dwi_node.inputs.out_file = os.path.join(
  1071. anat_output_dir,
  1072. rename_bids_file(synthseg_target_output, {"space": space_entity}, "synthseg", ".nii.gz")
  1073. )
  1074. dwi_workflow.connect(invert_dwi_to_t1w_reg_node, "out_file", synthseg_to_dwi_node, "in_matrix_file")
  1075. dwi_workflow.connect(preproc_dwi_node, "b0", synthseg_to_dwi_node, "reference")
  1076. wm_source_available = True
  1077. wm_source_type = "synthseg"
  1078. wm_source_node = synthseg_to_dwi_node
  1079. wm_source_output = "out_file"
  1080. elif fs_output_process:
  1081. wm_source_available = True
  1082. wm_source_type = "freesurfer"
  1083. wm_source_node = refine_aseg_node
  1084. wm_source_output = "out_aseg"
  1085. else:
  1086. wm_source_type = None
  1087. print("[DWI Pipeline] No valid WM source found. Skip NAWM mask creation.")
  1088. if wm_source_available:
  1089. create_WM_masks_node = MapNode(
  1090. MRIBinarize(), name="create_WM_masks", iterfield=["output_volume", "match"]
  1091. )
  1092. dwi_workflow.connect(wm_source_node, wm_source_output, create_WM_masks_node, "input_volume")
  1093. if wm_source_type == "synthseg":
  1094. create_WM_masks_node.inputs.match = [
  1095. [2, 41],
  1096. [2, 41, 7, 46],
  1097. [2], # Left-Cerebral-WM only
  1098. [41], # Right-Cerebral-WM only
  1099. ]
  1100. else:
  1101. create_WM_masks_node.inputs.match = [
  1102. [2, 41, 77, 251, 252, 253, 254, 255],
  1103. [2, 41, 77, 251, 252, 253, 254, 255, 7, 46],
  1104. [2], # Left-Cerebral-WM only
  1105. [41], # Right-Cerebral-WM only
  1106. ]
  1107. create_WM_masks_node.inputs.output_volume = [
  1108. os.path.join(anat_output_dir, f"{base_name}_label-cerebralWM_mask.nii.gz"),
  1109. os.path.join(anat_output_dir, f"{base_name}_label-WM_mask.nii.gz"),
  1110. os.path.join(anat_output_dir, f"{base_name}_label-rawNAWMLeft_mask.nii.gz"),
  1111. os.path.join(anat_output_dir, f"{base_name}_label-rawNAWMRight_mask.nii.gz"),
  1112. ]
  1113. def collect_exclusion_masks(seed_mask, wmh_mask, use_seed, use_wmh):
  1114. import os
  1115. masks = []
  1116. if use_seed and seed_mask and os.path.exists(seed_mask):
  1117. masks.append(seed_mask)
  1118. if use_wmh and wmh_mask and os.path.exists(wmh_mask):
  1119. masks.append(wmh_mask)
  1120. return masks
  1121. collect_excl_node = Node(Function(
  1122. input_names=["seed_mask", "wmh_mask", "use_seed", "use_wmh"],
  1123. output_names=["masks"],
  1124. function=collect_exclusion_masks
  1125. ), name="collect_exclusion_masks")
  1126. collect_excl_node.inputs.use_seed = bool(self.exclude_seed_mask and seed_mask != "")
  1127. collect_excl_node.inputs.use_wmh = bool(self.exclude_wmh_mask)
  1128. if self.exclude_seed_mask and seed_mask != "":
  1129. dwi_workflow.connect(seed_mask_to_dwi_node, "out_file", collect_excl_node, "seed_mask")
  1130. else:
  1131. collect_excl_node.inputs.seed_mask = ""
  1132. if self.exclude_wmh_mask or self.visual_pathway_analysis:
  1133. wmh_mask_to_dwi_node = Node(FLIRT(), name="wmh_mask_to_dwi")
  1134. wmh_mask_to_dwi_node.inputs.in_file = wmh_mask_file
  1135. wmh_mask_to_dwi_node.inputs.interp = "nearestneighbour"
  1136. wmh_mask_to_dwi_node.inputs.apply_xfm = True
  1137. wmh_mask_to_dwi_node.inputs.out_file = os.path.join(
  1138. anat_output_dir,
  1139. rename_bids_file(wmh_mask_file, {"space": space_entity}, "mask", ".nii.gz")
  1140. )
  1141. dwi_workflow.connect(invert_dwi_to_t1w_reg_node, "out_file", wmh_mask_to_dwi_node, "in_matrix_file")
  1142. dwi_workflow.connect(preproc_dwi_node, "b0", wmh_mask_to_dwi_node, "reference")
  1143. dwi_workflow.connect(wmh_mask_to_dwi_node, "out_file", collect_excl_node, "wmh_mask")
  1144. else:
  1145. collect_excl_node.inputs.wmh_mask = ""
  1146. def merge_exclusion_masks(mask_files, ref_file, out_file):
  1147. import os
  1148. import subprocess
  1149. mask_files = [] if mask_files is None else [m for m in mask_files if m]
  1150. os.makedirs(os.path.dirname(out_file), exist_ok=True)
  1151. if len(mask_files) == 0:
  1152. subprocess.check_call(["fslmaths", ref_file, "-mul", "0", out_file])
  1153. return out_file
  1154. cmd = ["fslmaths", mask_files[0]]
  1155. for mask_file in mask_files[1:]:
  1156. cmd += ["-add", mask_file]
  1157. cmd += ["-thr", "0.5", "-bin", out_file]
  1158. subprocess.check_call(cmd)
  1159. return out_file
  1160. merge_masks_to_exclude_node = Node(Function(
  1161. input_names=["mask_files", "ref_file", "out_file"],
  1162. output_names=["out_file"],
  1163. function=merge_exclusion_masks
  1164. ), name="merge_masks_to_exclude")
  1165. merge_masks_to_exclude_node.inputs.out_file = os.path.join(anat_output_dir, f"{base_name}_desc-exclusionMerged_mask.nii.gz")
  1166. dwi_workflow.connect(collect_excl_node, "masks", merge_masks_to_exclude_node, "mask_files")
  1167. dwi_workflow.connect(preproc_dwi_node, "b0", merge_masks_to_exclude_node, "ref_file")
  1168. # Left NAWM
  1169. get_raw_nawm_left_mask_node = Node(Select(index=2), name="select_raw_nawm_left_mask")
  1170. dwi_workflow.connect(create_WM_masks_node, "output_volume", get_raw_nawm_left_mask_node, "inlist")
  1171. get_final_nawm_left_node = Node(GenerateNAWMMask(), name="get_final_nawm_left_mask")
  1172. get_final_nawm_left_node.inputs.erode_mm = 2
  1173. get_final_nawm_left_node.inputs.output_mask = os.path.join(anat_output_dir, f"{base_name}_label-NAWMLeft_mask.nii.gz")
  1174. dwi_workflow.connect(get_raw_nawm_left_mask_node, "out", get_final_nawm_left_node, "wm_mask")
  1175. dwi_workflow.connect(merge_masks_to_exclude_node, "out_file", get_final_nawm_left_node, "exclude_mask")
  1176. # Right NAWM
  1177. get_raw_nawm_right_mask_node = Node(Select(index=3), name="select_raw_nawm_right_mask")
  1178. dwi_workflow.connect(create_WM_masks_node, "output_volume", get_raw_nawm_right_mask_node, "inlist")
  1179. get_final_nawm_right_node = Node(GenerateNAWMMask(), name="get_final_nawm_right_mask")
  1180. get_final_nawm_right_node.inputs.erode_mm = 2
  1181. get_final_nawm_right_node.inputs.output_mask = os.path.join(anat_output_dir, f"{base_name}_label-NAWMRight_mask.nii.gz")
  1182. dwi_workflow.connect(get_raw_nawm_right_mask_node, "out", get_final_nawm_right_node, "wm_mask")
  1183. dwi_workflow.connect(merge_masks_to_exclude_node, "out_file", get_final_nawm_right_node, "exclude_mask")
  1184. # Merge left and right into a single NAWM dseg (1=left, 2=right)
  1185. def merge_nawm_dseg(left_mask, right_mask, out_file):
  1186. import os, nibabel, numpy as np
  1187. os.makedirs(os.path.dirname(out_file), exist_ok=True)
  1188. left_img = nibabel.load(left_mask)
  1189. left_data = left_img.get_fdata().astype(bool)
  1190. right_data = nibabel.load(right_mask).get_fdata().astype(bool)
  1191. dseg_data = np.zeros(left_data.shape, dtype=np.uint8)
  1192. dseg_data[left_data] = 1
  1193. dseg_data[right_data] = 2
  1194. dseg_img = nibabel.Nifti1Image(dseg_data, left_img.affine, left_img.header)
  1195. dseg_img.header.set_data_dtype(np.uint8)
  1196. nibabel.save(dseg_img, out_file)
  1197. return out_file
  1198. merge_nawm_dseg_node = Node(Function(
  1199. input_names=["left_mask", "right_mask", "out_file"],
  1200. output_names=["out_file"],
  1201. function=merge_nawm_dseg
  1202. ), name="merge_nawm_dseg")
  1203. nawm_dseg = os.path.join(anat_output_dir, f"{base_name}_label-NAWM_dseg.nii.gz")
  1204. merge_nawm_dseg_node.inputs.out_file = nawm_dseg
  1205. dwi_workflow.connect(get_final_nawm_left_node, "nawm_mask", merge_nawm_dseg_node, "left_mask")
  1206. dwi_workflow.connect(get_final_nawm_right_node, "nawm_mask", merge_nawm_dseg_node, "right_mask")
  1207. else:
  1208. nawm_dseg = None
  1209. # ============================================
  1210. # DSI Studio Reconstruction
  1211. # ============================================
  1212. if self.dsistudio_gqi or self.dsistudio_qsdr:
  1213. # first need to convert to .src.gz format
  1214. os.makedirs(os.path.join(self.output_path, 'dsistudio'), exist_ok=True)
  1215. create_src_node = Node(CreateSRC(), name='create_src')
  1216. dwi_workflow.connect(preproc_dwi_node, 'preproc_dwi', create_src_node, 'source')
  1217. create_src_node.inputs.output = os.path.join(self.output_path, 'dsistudio', rename_bids_file(preproc_dwi_filename, {"desc": None, "desc": "preproc"}, 'dwi', '.src.gz'))
  1218. dti_out = os.path.join(self.output_path, 'dsistudio', rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "dti"}, 'dwimap', '.fib.gz'))
  1219. gqi_out = os.path.join(self.output_path, 'dsistudio', rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "gqi"}, 'dwimap', '.fib.gz'))
  1220. qsdr_out = os.path.join(self.output_path, 'dsistudio', rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "qsdr"}, 'dwimap', '.fib.gz'))
  1221. if self.dsistudio_dti:
  1222. dti_reconstruction_node = Node(DSIstudioReconstruction(), name='dti_reconstruction')
  1223. dwi_workflow.connect(create_src_node, 'out_file', dti_reconstruction_node, 'source')
  1224. dti_reconstruction_node.inputs.method = 1
  1225. dti_reconstruction_node.inputs.output = dti_out
  1226. dti_reconstruction_node.inputs.other_output = 'fa,ad,rd,md,tensor'
  1227. #dti_reconstruction_node.inputs.check_btable = 1
  1228. #dwi_workflow.connect(preproc_dwi_node, 'dwi_mask', dti_reconstruction_node, 'mask')
  1229. if self.dsistudio_gqi:
  1230. gqi_reconstruction_node = Node(DSIstudioReconstruction(), name='gqi_reconstruction')
  1231. dwi_workflow.connect(create_src_node, 'out_file', gqi_reconstruction_node, 'source')
  1232. gqi_reconstruction_node.inputs.method = 4
  1233. gqi_reconstruction_node.inputs.other_output = 'fa,ad,rd,md,iso,rdi,nrdi,tensor'
  1234. gqi_reconstruction_node.inputs.output = gqi_out
  1235. gqi_reconstruction_node.inputs.check_btable = 1
  1236. #dwi_workflow.connect(preproc_dwi_node, 'dwi_mask', gqi_reconstruction_node, 'mask')
  1237. if self.dsistudio_qsdr:
  1238. qsdr_reconstruction_node = Node(DSIstudioReconstructionOldVer(), name='qsdr_reconstruction')
  1239. dwi_workflow.connect(create_src_node, 'out_file', qsdr_reconstruction_node, 'source')
  1240. qsdr_reconstruction_node.inputs.method = 7
  1241. qsdr_reconstruction_node.inputs.qsdr_reso = 2.0
  1242. qsdr_reconstruction_node.inputs.other_output = 'fa,ad,rd,md,iso,rdi,nrdi,tensor'
  1243. qsdr_reconstruction_node.inputs.output = qsdr_out
  1244. qsdr_reconstruction_node.inputs.thread_count = 8
  1245. qsdr_reconstruction_node.inputs.check_btable = 1
  1246. #dwi_workflow.connect(preproc_dwi_node, 'dwi_mask', qsdr_reconstruction_node, 'mask')
  1247. if self.preprocess_method == 'post_qsiprep':
  1248. qsirecon_output_dir = self.session._find_output('qsirecon-DSIStudio')
  1249. gqi_fib = os.path.join(qsirecon_output_dir, 'dwi', rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "gqi"}, 'dwimap', '.fib.gz'))
  1250. if not os.path.exists(gqi_fib):
  1251. raise FileNotFoundError(f"GQI .fib.gz file not found in {qsirecon_output_dir}. Please run QSIRecon with DSIStudio reconstruction first.")
  1252. dwi_scalarmaps_output_node.inputs.gqi_gfa_img = os.path.join(qsirecon_output_dir, 'dwi', rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "gqi", "param": "gfa"}, 'dwimap', '.nii.gz'))
  1253. dwi_scalarmaps_output_node.inputs.gqi_iso_img = os.path.join(qsirecon_output_dir, 'dwi', rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "gqi", "param": "iso"}, 'dwimap', '.nii.gz'))
  1254. dwi_scalarmaps_output_node.inputs.gqi_qa_img = os.path.join(qsirecon_output_dir, 'dwi', rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "gqi", "param": "qa"}, 'dwimap', '.nii.gz'))
  1255. # ===========================================
  1256. # AMICO NODDI Reconstruction
  1257. # ===========================================
  1258. if self.amico_noddi:
  1259. from cvdproc.pipelines.dmri.amico.amico_nipype import AmicoNoddi
  1260. amico_output_dir = os.path.join(self.output_path, 'NODDI')
  1261. amico_noddi_node = Node(AmicoNoddi(), name='amico_noddi')
  1262. dwi_workflow.connect(preproc_dwi_node, 'preproc_dwi', amico_noddi_node, 'dwi')
  1263. dwi_workflow.connect(preproc_dwi_node, 'bval', amico_noddi_node, 'bval')
  1264. dwi_workflow.connect(preproc_dwi_node, 'bvec', amico_noddi_node, 'bvec')
  1265. dwi_workflow.connect(preproc_dwi_node, 'dwi_mask', amico_noddi_node, 'mask')
  1266. amico_noddi_node.inputs.output_dir = amico_output_dir
  1267. amico_noddi_node.inputs.direction_filename = rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "noddi", 'param': 'direction'}, 'dwimap', '.nii.gz')
  1268. amico_noddi_node.inputs.icvf_filename = rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "noddi", 'param': 'icvf'}, 'dwimap', '.nii.gz')
  1269. amico_noddi_node.inputs.isovf_filename = rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "noddi", 'param': 'isovf'}, 'dwimap', '.nii.gz')
  1270. amico_noddi_node.inputs.od_filename = rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "noddi", 'param': 'odi'}, 'dwimap', '.nii.gz')
  1271. amico_noddi_node.inputs.modulated_icvf_filename = rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "noddi", 'param': 'ModulatedICVF'}, 'dwimap', '.nii.gz')
  1272. amico_noddi_node.inputs.modulated_od_filename = rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "noddi", 'param': 'ModulatedODI'}, 'dwimap', '.nii.gz')
  1273. amico_noddi_node.inputs.config_filename = rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "noddi"}, 'config', '.pickle')
  1274. # Connect outputs to dwi_scalarmaps_output_node
  1275. dwi_workflow.connect(amico_noddi_node, 'icvf', dwi_scalarmaps_output_node, 'icvf_img')
  1276. dwi_workflow.connect(amico_noddi_node, 'isovf', dwi_scalarmaps_output_node, 'isovf_img')
  1277. dwi_workflow.connect(amico_noddi_node, 'od', dwi_scalarmaps_output_node, 'odi_img')
  1278. else:
  1279. # if these files already exist, also set in dwi_scalarmaps_output_node
  1280. amico_output_dir = os.path.join(self.output_path, 'NODDI')
  1281. icvf_img = os.path.join(amico_output_dir, rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "noddi", 'param': 'icvf'}, 'dwimap', '.nii.gz'))
  1282. isovf_img = os.path.join(amico_output_dir, rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "noddi", 'param': 'isovf'}, 'dwimap', '.nii.gz'))
  1283. od_img = os.path.join(amico_output_dir, rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "noddi", 'param': 'odi'}, 'dwimap', '.nii.gz'))
  1284. if os.path.exists(icvf_img) and os.path.exists(isovf_img) and os.path.exists(od_img):
  1285. dwi_scalarmaps_output_node.inputs.icvf_img = icvf_img
  1286. dwi_scalarmaps_output_node.inputs.isovf_img = isovf_img
  1287. dwi_scalarmaps_output_node.inputs.odi_img = od_img
  1288. print('[DWI Pipeline] NODDI scalar maps already exist, using them.')
  1289. # ===========================================
  1290. # FSL FDT bedpostx
  1291. # ===========================================
  1292. if fdt_bedpostx:
  1293. bedpostx_output_dir = os.path.join(self.output_path, 'bedpostx')
  1294. bedpostx_node = Node(BedpostxGPUCustom(), name="bedpostx")
  1295. dwi_workflow.connect(preproc_dwi_node, "preproc_dwi", bedpostx_node, "dwi_img")
  1296. dwi_workflow.connect(preproc_dwi_node, "bvec", bedpostx_node, "bvec")
  1297. dwi_workflow.connect(preproc_dwi_node, "bval", bedpostx_node, "bval")
  1298. dwi_workflow.connect(preproc_dwi_node, "dwi_mask", bedpostx_node, "mask")
  1299. bedpostx_node.inputs.out_dir = bedpostx_output_dir
  1300. # ===========================================
  1301. # Connectome
  1302. # ===========================================
  1303. if self.connectome:
  1304. connectome_output_dir = os.path.join(self.output_path, 'connectome')
  1305. os.makedirs(connectome_output_dir, exist_ok=True)
  1306. # Prepare connectome node
  1307. prepare_connectome_node = Node(ConnectomePrepare(), name='prepare_connectome')
  1308. dwi_workflow.connect(preproc_dwi_mif_node, 'preproc_dwi', prepare_connectome_node, 'preproc_dwi_mif')
  1309. dwi_workflow.connect(preproc_dwi_mif_node, 'dwi_mask', prepare_connectome_node, 'dwi_mask_mif')
  1310. prepare_connectome_node.inputs.output_dir = connectome_output_dir
  1311. dwi_workflow.connect(refine_aseg_node, 'out_aseg', prepare_connectome_node, 'aseg')
  1312. prepare_connectome_node.inputs.wm_response = rename_bids_file(preproc_dwi_filename, {'param': 'fod', 'label': 'WM', 'desc': None}, 'dwimap', '.txt')
  1313. prepare_connectome_node.inputs.wm_fod = rename_bids_file(preproc_dwi_filename, {'param': 'fod', 'label': 'WM', 'desc': None}, 'dwimap', '.mif')
  1314. prepare_connectome_node.inputs.wm_fod_norm = rename_bids_file(preproc_dwi_filename, {'param': 'fod', 'label': 'WM', 'desc': 'norm'}, 'dwimap', '.mif')
  1315. prepare_connectome_node.inputs.gm_response = rename_bids_file(preproc_dwi_filename, {'param': 'fod', 'label': 'GM', 'desc': None}, 'dwimap', '.txt')
  1316. prepare_connectome_node.inputs.gm_fod = rename_bids_file(preproc_dwi_filename, {'param': 'fod', 'label': 'GM', 'desc': None}, 'dwimap', '.mif')
  1317. prepare_connectome_node.inputs.gm_fod_norm = rename_bids_file(preproc_dwi_filename, {'param': 'fod', 'label': 'GM', 'desc': 'norm'}, 'dwimap', '.mif')
  1318. prepare_connectome_node.inputs.csf_response = rename_bids_file(preproc_dwi_filename, {'param': 'fod', 'label': 'CSF', 'desc': None}, 'dwimap', '.txt')
  1319. prepare_connectome_node.inputs.csf_fod = rename_bids_file(preproc_dwi_filename, {'param': 'fod', 'label': 'CSF', 'desc': None}, 'dwimap', '.mif')
  1320. prepare_connectome_node.inputs.csf_fod_norm = rename_bids_file(preproc_dwi_filename, {'param': 'fod', 'label': 'CSF', 'desc': 'norm'}, 'dwimap', '.mif')
  1321. prepare_connectome_node.inputs.sift_mu = rename_bids_file(preproc_dwi_filename, {'model': 'sift2', 'desc': None}, 'mu', '.txt')
  1322. prepare_connectome_node.inputs.five_tt_dwi = rename_bids_file(preproc_dwi_filename, {'desc': None}, '5tt', '.mif')
  1323. prepare_connectome_node.inputs.gmwmSeed_dwi = rename_bids_file(preproc_dwi_filename, {'desc': None}, 'gmwmSeed', '.mif')
  1324. prepare_connectome_node.inputs.streamlines = rename_bids_file(preproc_dwi_filename, {'desc': None}, 'streamlines', '.tck')
  1325. prepare_connectome_node.inputs.sift_weights = rename_bids_file(preproc_dwi_filename, {'model': 'sift2', 'desc': None}, 'streamlineweights', '.txt')
  1326. # Prepare parcellation
  1327. fs_parcellation_config = Node(Merge(2), name='merge_fs_parcellation_config')
  1328. fs_parcellation_config.inputs.in1 = get_package_path('data', 'labelconvert', 'fs_default.txt')
  1329. fs_parcellation_config.inputs.in2 = get_package_path('data', 'labelconvert', 'custom', 'fs_aparc.txt')
  1330. prepare_fs_parcellation_node = MapNode(LabelConvert(), name='prepare_parcellation', iterfield=['in_config', 'out_file'])
  1331. dwi_workflow.connect(refine_aseg_node, 'out_aseg', prepare_fs_parcellation_node, 'in_file')
  1332. prepare_fs_parcellation_node.inputs.in_lut = get_package_path('data', 'labelconvert_in', 'FreeSurferColorLUT.txt')
  1333. prepare_fs_parcellation_node.inputs.args = '-force'
  1334. dwi_workflow.connect(fs_parcellation_config, 'out', prepare_fs_parcellation_node, 'in_config')
  1335. prepare_fs_parcellation_node.inputs.out_file = [
  1336. os.path.join(connectome_output_dir, rename_bids_file(preproc_dwi_filename, {'seg': 'aparcaseg', 'desc': None}, 'parcellation', '.mif')),
  1337. os.path.join(connectome_output_dir, rename_bids_file(preproc_dwi_filename, {'seg': 'aparc', 'desc': None}, 'parcellation', '.mif'))
  1338. ]
  1339. # Connectome construction
  1340. construct_connectome_node = MapNode(BuildConnectome(), name='construct_connectome', iterfield=['in_parc', 'out_file'])
  1341. dwi_workflow.connect(prepare_connectome_node, 'global_streamlines', construct_connectome_node, 'in_file')
  1342. construct_connectome_node.inputs.out_file = [
  1343. os.path.join(connectome_output_dir, rename_bids_file(preproc_dwi_filename, {'seg': 'aparcaseg', 'desc': None}, 'connectivity', '.csv')),
  1344. os.path.join(connectome_output_dir, rename_bids_file(preproc_dwi_filename, {'seg': 'aparc', 'desc': None}, 'connectivity', '.csv'))
  1345. ]
  1346. dwi_workflow.connect(prepare_fs_parcellation_node, 'out_file', construct_connectome_node, 'in_parc')
  1347. dwi_workflow.connect(prepare_connectome_node, 'sift_weights', construct_connectome_node, 'in_weights')
  1348. construct_connectome_node.inputs.args = '-zero_diagonal -symmetric -scale_invnodevol'
  1349. # ===========================================
  1350. # Tractography: Designed for lesion-based probabilistic tractography
  1351. # ===========================================
  1352. tractography_output_dir = os.path.join(self.output_path, 'tractography')
  1353. if 'fdt' in self.tractography:
  1354. fdt_tractography_output_dir = os.path.join(tractography_output_dir, 'fdt')
  1355. os.makedirs(fdt_tractography_output_dir, exist_ok=True)
  1356. uncorrected_tractography_output_dir = os.path.join(fdt_tractography_output_dir, 'PathLengthUncorrected')
  1357. seed_mask_tractography = Node(Probtrackx(), name='seed_mask_tractography')
  1358. dwi_workflow.connect(bedpostx_node, 'source_for_probtrackx', seed_mask_tractography, 'source')
  1359. dwi_workflow.connect(bedpostx_node, 'mask_for_probtrackx', seed_mask_tractography, 'dwi_mask')
  1360. dwi_workflow.connect(seed_mask_to_fs_node, 'output_image', seed_mask_tractography, 'seed_mask') # in fs space
  1361. dwi_workflow.connect(concat_fs_to_dwi_xfm_node, 'out_file', seed_mask_tractography, 'seed_to_dwi_xfm')
  1362. seed_mask_tractography.inputs.seed_ref = os.path.join(fs_output, 'mri', 'orig.mgz') # fs native space
  1363. dwi_workflow.connect(create_cortical_GM_mask_node, 'output_gm_mask', seed_mask_tractography, 'waypoints') # in fs space
  1364. seed_mask_tractography.inputs.path_length_correction = False
  1365. seed_mask_tractography.inputs.output_dir = uncorrected_tractography_output_dir
  1366. seed_mask_tractography.inputs.nsamples = 10000
  1367. seed_mask_tractography.inputs.args = '-l --forcedir --opd --ompl'
  1368. waytotal_node = Node(ReadSingleValue(), name='waytotal_node')
  1369. dwi_workflow.connect(seed_mask_tractography, 'waytotal', waytotal_node, 'input_file')
  1370. lh_fdtpaths_to_surf = Node(MRIvol2surf(), name='lh_fdtpaths_to_surf')
  1371. dwi_workflow.connect(seed_mask_tractography, 'fdt_paths', lh_fdtpaths_to_surf, 'volume')
  1372. lh_fdtpaths_to_surf.inputs.subjects_dir = fs_subjects_dir
  1373. lh_fdtpaths_to_surf.inputs.regheader = fs_subject_id
  1374. lh_fdtpaths_to_surf.inputs.hemi = 'lh'
  1375. lh_fdtpaths_to_surf.inputs.output_surf = os.path.join(uncorrected_tractography_output_dir, f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_space-fsaverage_TDI.mgh')
  1376. lh_fdtpaths_to_surf.inputs.proj_frac = 0.5
  1377. lh_fdtpaths_to_surf.inputs.target = 'fsaverage'
  1378. mask_lh_fdtpaths_surf = Node(ApplyIndexTxtToMgh(), name='mask_lh_fdtpaths_surf')
  1379. dwi_workflow.connect(lh_fdtpaths_to_surf, 'output_surf', mask_lh_fdtpaths_surf, 'measure_mgh')
  1380. mask_lh_fdtpaths_surf.inputs.index_txt = get_package_path('data', 'standard', 'fsaverage', 'lh.aparc.label_medial_wall.txt')
  1381. mask_lh_fdtpaths_surf.inputs.output_mgh = os.path.join(uncorrected_tractography_output_dir, f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_space-fsaverage_desc-NoMedialWall_TDI.mgh')
  1382. rh_fdtpaths_to_surf = Node(MRIvol2surf(), name='rh_fdtpaths_to_surf')
  1383. dwi_workflow.connect(seed_mask_tractography, 'fdt_paths', rh_fdtpaths_to_surf, 'volume')
  1384. rh_fdtpaths_to_surf.inputs.subjects_dir = fs_subjects_dir
  1385. rh_fdtpaths_to_surf.inputs.regheader = fs_subject_id
  1386. rh_fdtpaths_to_surf.inputs.hemi = 'rh'
  1387. rh_fdtpaths_to_surf.inputs.output_surf = os.path.join(uncorrected_tractography_output_dir, f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_space-fsaverage_TDI.mgh')
  1388. rh_fdtpaths_to_surf.inputs.proj_frac = 0.5
  1389. rh_fdtpaths_to_surf.inputs.target = 'fsaverage'
  1390. mask_rh_fdtpaths_surf = Node(ApplyIndexTxtToMgh(), name='mask_rh_fdtpaths_surf')
  1391. dwi_workflow.connect(rh_fdtpaths_to_surf, 'output_surf', mask_rh_fdtpaths_surf, 'measure_mgh')
  1392. mask_rh_fdtpaths_surf.inputs.index_txt = get_package_path('data', 'standard', 'fsaverage', 'rh.aparc.label_medial_wall.txt')
  1393. mask_rh_fdtpaths_surf.inputs.output_mgh = os.path.join(uncorrected_tractography_output_dir, f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_space-fsaverage_desc-NoMedialWall_TDI.mgh')
  1394. define_connection_level_node = Node(DefineConnectionLevel(), name='define_connection_levels')
  1395. dwi_workflow.connect(mask_lh_fdtpaths_surf, 'output_mgh', define_connection_level_node, 'lh_mgh')
  1396. dwi_workflow.connect(mask_rh_fdtpaths_surf, 'output_mgh', define_connection_level_node, 'rh_mgh')
  1397. define_connection_level_node.inputs.low_threshold = 3.8e-5
  1398. dwi_workflow.connect(waytotal_node, 'value', define_connection_level_node, 'divisor')
  1399. define_connection_level_node.inputs.output_files = {
  1400. 'lh_unconn': os.path.join(uncorrected_tractography_output_dir, f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_space-fsaverage_desc-UnConn_mask.mgh'),
  1401. 'rh_unconn': os.path.join(uncorrected_tractography_output_dir, f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_space-fsaverage_desc-UnConn_mask.mgh'),
  1402. 'lh_lowconn': os.path.join(uncorrected_tractography_output_dir, f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_space-fsaverage_desc-LowConn_mask.mgh'),
  1403. 'rh_lowconn': os.path.join(uncorrected_tractography_output_dir, f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_space-fsaverage_desc-LowConn_mask.mgh'),
  1404. 'lh_medconn': os.path.join(uncorrected_tractography_output_dir, f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_space-fsaverage_desc-MediumConn_mask.mgh'),
  1405. 'rh_medconn': os.path.join(uncorrected_tractography_output_dir, f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_space-fsaverage_desc-MediumConn_mask.mgh'),
  1406. 'lh_highconn': os.path.join(uncorrected_tractography_output_dir, f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_space-fsaverage_desc-HighConn_mask.mgh'),
  1407. 'rh_highconn': os.path.join(uncorrected_tractography_output_dir, f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_space-fsaverage_desc-HighConn_mask.mgh'),
  1408. }
  1409. # remask unconn regions
  1410. mask_lh_unconn_node = Node(ApplyGiiMaskToMgh(), name='mask_lh_unconn_regions')
  1411. dwi_workflow.connect(define_connection_level_node, 'lh_unconn_mask', mask_lh_unconn_node, 'measure_mgh')
  1412. mask_lh_unconn_node.inputs.mask_gii = get_package_path('data', 'standard', 'fsaverage', 'lh.thickness.shape.gii')
  1413. mask_lh_unconn_node.inputs.output_mgh = os.path.join(uncorrected_tractography_output_dir, f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_space-fsaverage_desc-UnConn_mask.mgh')
  1414. mask_rh_unconn_node = Node(ApplyGiiMaskToMgh(), name='mask_rh_unconn_regions')
  1415. dwi_workflow.connect(define_connection_level_node, 'rh_unconn_mask', mask_rh_unconn_node, 'measure_mgh')
  1416. mask_rh_unconn_node.inputs.mask_gii = get_package_path('data', 'standard', 'fsaverage', 'rh.thickness.shape.gii')
  1417. mask_rh_unconn_node.inputs.output_mgh = os.path.join(uncorrected_tractography_output_dir, f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_space-fsaverage_desc-UnConn_mask.mgh')
  1418. # extract surface parameters
  1419. uncorr_extract_surf_params = Node(ExtractSurfaceParameters(), name='uncorr_extract_surf_params')
  1420. uncorr_extract_surf_params.inputs.fs_subjects_dir = fs_subjects_dir
  1421. uncorr_extract_surf_params.inputs.sessions = self.subject.sessions_id
  1422. if self.use_freesurfer_longitudinal:
  1423. uncorr_extract_surf_params.inputs.sessions = self.subject.sessions_id + [f'long.sub-{self.subject.subject_id}']
  1424. dwi_workflow.connect(mask_lh_unconn_node, 'output_mgh', uncorr_extract_surf_params, 'lh_unconn_mask')
  1425. dwi_workflow.connect(mask_rh_unconn_node, 'output_mgh', uncorr_extract_surf_params, 'rh_unconn_mask')
  1426. dwi_workflow.connect(define_connection_level_node, 'lh_lowconn_mask', uncorr_extract_surf_params, 'lh_low_conn_mask')
  1427. dwi_workflow.connect(define_connection_level_node, 'rh_lowconn_mask', uncorr_extract_surf_params, 'rh_low_conn_mask')
  1428. dwi_workflow.connect(define_connection_level_node, 'lh_medconn_mask', uncorr_extract_surf_params, 'lh_medium_conn_mask')
  1429. dwi_workflow.connect(define_connection_level_node, 'rh_medconn_mask', uncorr_extract_surf_params, 'rh_medium_conn_mask')
  1430. dwi_workflow.connect(define_connection_level_node, 'lh_highconn_mask', uncorr_extract_surf_params, 'lh_high_conn_mask')
  1431. dwi_workflow.connect(define_connection_level_node, 'rh_highconn_mask', uncorr_extract_surf_params, 'rh_high_conn_mask')
  1432. uncorr_extract_surf_params.inputs.output_dir = uncorrected_tractography_output_dir
  1433. uncorr_extract_surf_params.inputs.csv_file_name = f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_desc-ConnectivityLevels_surfpars.csv'
  1434. if 'mrtrix3' in self.tractography:
  1435. seed_based_track_node = Node(IdentityInterface(fields=['seed_based_track']), name='seed_based_track')
  1436. tckgen_output_dir = os.path.join(tractography_output_dir, 'mrtrix3')
  1437. tckgen_node = Node(Tractography(), name='tckgen')
  1438. os.makedirs(tckgen_output_dir, exist_ok=True)
  1439. dwi_workflow.connect(preproc_dwi_mif_node, 'preproc_dwi', tckgen_node, 'in_file')
  1440. dwi_workflow.connect(preproc_dwi_node, 'dwi_mask', tckgen_node, 'roi_mask')
  1441. dwi_workflow.connect(seed_mask_to_dwi_node, 'out_file', tckgen_node, 'seed_image')
  1442. tckgen_node.inputs.out_file = os.path.join(tckgen_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': 'ROIseed'}, 'streamlines', '.tck'))
  1443. tckgen_node.inputs.algorithm = 'Tensor_Prob'
  1444. tckgen_node.inputs.select = 10000
  1445. tckgen_node.inputs.nthreads = 4
  1446. tckgen_node.inputs.args = '-force'
  1447. dwi_workflow.connect(tckgen_node, 'out_file', seed_based_track_node, 'seed_based_track')
  1448. if nawm_dseg is not None:
  1449. nawm_without_tract = os.path.join(tckgen_output_dir, rename_bids_file(nawm_dseg, {'desc': 'ConnTractRemoved'}, 'dseg', '.nii.gz'))
  1450. remove_tract_region_node = Node(RemoveTractRegion(), name='remove_tract_region')
  1451. dwi_workflow.connect(seed_based_track_node, 'seed_based_track', remove_tract_region_node, 'tck_file')
  1452. dwi_workflow.connect(merge_nawm_dseg_node, 'out_file', remove_tract_region_node, 'wm_mask')
  1453. remove_tract_region_node.inputs.out_tract_mask = os.path.join(tckgen_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': 'ROIseedTract'}, 'mask', '.nii.gz'))
  1454. remove_tract_region_node.inputs.out_tdi_norm = os.path.join(tckgen_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': 'ROIseedTract'}, 'tdi', '.nii.gz'))
  1455. remove_tract_region_node.inputs.out_wm_mask = nawm_without_tract
  1456. # ===========================================
  1457. # DTI-ALPS
  1458. # ===========================================
  1459. #if self.dtialps and self.preprocess_method != 'post_qsiprep':
  1460. if self.dtialps:
  1461. print("[DWI Pipeline] DTI-ALPS analysis: True")
  1462. dti_alps_node = Node(ALPS(), name='dti_alps')
  1463. dwi_workflow.connect(dti_fit_output_node, 'fa_img', dti_alps_node, 'fa_img')
  1464. dti_alps_node.inputs.output_dir = os.path.join(self.output_path, 'DTI-ALPS')
  1465. dti_alps_node.inputs.alps_dir = get_package_path("pipelines", "external", "alps")
  1466. # self.dtialps_register_method: 1 for 'flirt', 2 for 'synthmorph'
  1467. dti_alps_node.inputs.register_method = 'flirt' if self.dtialps_register_method == 1 else 'synthmorph' if self.dtialps_register_method == 2 else 'ants'
  1468. dwi_workflow.connect(dti_fit_output_node, 'tensor_img', dti_alps_node, 'tensor_img')
  1469. # BIDS output filenames in MNI space
  1470. mni_entities = {"space": "MNI152NLin6ASym", "desc": None, "model": "tensor"}
  1471. dti_alps_node.inputs.fa_out = os.path.join(dti_alps_node.inputs.output_dir, rename_bids_file(preproc_dwi_filename, {**mni_entities, "param": "fa"}, "dwimap", ".nii.gz"))
  1472. dti_alps_node.inputs.xx_out = os.path.join(dti_alps_node.inputs.output_dir, rename_bids_file(preproc_dwi_filename, {**mni_entities, "param": "xx"}, "dwimap", ".nii.gz"))
  1473. dti_alps_node.inputs.yy_out = os.path.join(dti_alps_node.inputs.output_dir, rename_bids_file(preproc_dwi_filename, {**mni_entities, "param": "yy"}, "dwimap", ".nii.gz"))
  1474. dti_alps_node.inputs.zz_out = os.path.join(dti_alps_node.inputs.output_dir, rename_bids_file(preproc_dwi_filename, {**mni_entities, "param": "zz"}, "dwimap", ".nii.gz"))
  1475. # If warp already exists
  1476. target_warp = os.path.join(self.subject.bids_dir, 'derivatives', 'xfm', f'sub-{self.subject.subject_id}', f'ses-{self.session.session_id}', f'sub-{self.subject.subject_id}_ses-{self.session.session_id}_from-T1w_to-MNI152NLin6ASym_warp.nii.gz')
  1477. if self.dtialps_register_method == 2 and os.path.exists(target_warp) and self.dwi_t1w_register:
  1478. dwi_workflow.connect(inputnode, 't1w_file', dti_alps_node, 't1_img')
  1479. dwi_workflow.connect(dwi_to_t1w_reg_node, 'out_matrix_file', dti_alps_node, 'fa_to_t1w_affine')
  1480. dti_alps_node.inputs.t1_to_mni_warp = target_warp
  1481. # ===========================================
  1482. # PVeD
  1483. # ===========================================
  1484. if self.pved:
  1485. print(f"[DWI Pipeline] PVeD analysis: True (method={self.pved_method})")
  1486. if self.pved_method == 'v1':
  1487. pved_node = Node(PVeD(), name='pved_v1')
  1488. dwi_workflow.connect(qsdr_reconstruction_node, 'out_file', pved_node, 'qsdr_fib_file')
  1489. pved_node.inputs.output_dir = os.path.join(self.output_path, 'PVeD')
  1490. pved_node.inputs.script_path = get_package_path("pipelines", "matlab", "pved", "pved_single.m")
  1491. pved_node.inputs.spm_path = get_package_path("data", "matlab_toolbox", "spm12")
  1492. pved_node.inputs.pved_path = get_package_path("data", "matlab_toolbox", "EstPVeD", 'MATLAB', 'pkg_pved_est')
  1493. else:
  1494. pved_prepare_node = Node(PVeDPrepareMNI(), name='pved_prepare_mni')
  1495. dwi_workflow.connect(dti_fit_output_node, 'fa_img', pved_prepare_node, 'fa_file')
  1496. dwi_workflow.connect(dti_fit_output_node, 'md_img', pved_prepare_node, 'md_file')
  1497. dwi_workflow.connect(dti_fit_output_node, 'tensor_img', pved_prepare_node, 'tensor_file')
  1498. dwi_workflow.connect(dwi_to_t1w_reg_node, 'out_matrix_file', pved_prepare_node, 'dwi_to_t1w_affine')
  1499. dwi_workflow.connect(inputnode, 't1w_file', pved_prepare_node, 't1w_reference')
  1500. pved_prepare_node.inputs.t1w_to_mni_warp = pved_t1w_to_mni_warp
  1501. pved_prepare_node.inputs.output_dir = os.path.join(self.output_path, 'PVeD')
  1502. pved_prepare_node.inputs.subject_id = f"sub-{self.subject.subject_id}"
  1503. pved_prepare_node.inputs.session_id = f"ses-{self.session.session_id}"
  1504. pved_node = Node(PVeDAtropos(), name='pved_v2')
  1505. dwi_workflow.connect(pved_prepare_node, 'fa_file', pved_node, 'fa_file')
  1506. dwi_workflow.connect(pved_prepare_node, 'md_file', pved_node, 'md_file')
  1507. dwi_workflow.connect(pved_prepare_node, 'tensor_file', pved_node, 'tensor_file')
  1508. pved_node.inputs.output_dir = os.path.join(self.output_path, 'PVeD')
  1509. pved_node.inputs.subject_id = f"sub-{self.subject.subject_id}"
  1510. pved_node.inputs.session_id = f"ses-{self.session.session_id}"
  1511. # ===========================================
  1512. # Free Water
  1513. # ===========================================
  1514. freewater_node = Node(IdentityInterface(fields=['single_shell_fw_img', 'markvcid_fw_img', 'dti_fw_img']), name='freewater_node')
  1515. if 'single_shell_freewater' in self.freewater:
  1516. single_shell_freewater_node = Node(SingleShellFW(), name='single_shell_freewater')
  1517. dwi_workflow.connect(inputnode, 'bval_file', single_shell_freewater_node, 'fbval')
  1518. dwi_workflow.connect(preproc_dwi_node, 'preproc_dwi', single_shell_freewater_node, 'fdwi')
  1519. dwi_workflow.connect(preproc_dwi_node, 'bvec', single_shell_freewater_node, 'fbvec')
  1520. dwi_workflow.connect(preproc_dwi_node, 'dwi_mask', single_shell_freewater_node, 'mask_file')
  1521. fw_output_path = os.path.join(self.output_path, 'freewater', 'single_shell_freewater')
  1522. single_shell_freewater_node.inputs.working_directory = fw_output_path
  1523. single_shell_freewater_node.inputs.output_directory = fw_output_path
  1524. single_shell_freewater_node.inputs.crop_shells = False
  1525. dwi_workflow.connect(single_shell_freewater_node, 'output_fw', freewater_node, 'fw_img')
  1526. else:
  1527. freewater_node.inputs.fw_img = os.path.join(self.output_path, 'single_shell_freewater', 'freewater.nii.gz')
  1528. if 'markvcid_freewater' in self.freewater:
  1529. markvcid_freewater_node = Node(MarkVCIDFreeWater(), name='markvcid_freewater')
  1530. dwi_workflow.connect(preproc_dwi_node, 'preproc_dwi', markvcid_freewater_node, 'in_dwi')
  1531. dwi_workflow.connect(preproc_dwi_node, 'bvec', markvcid_freewater_node, 'in_dwi_bvec')
  1532. dwi_workflow.connect(preproc_dwi_node, 'bval', markvcid_freewater_node, 'in_dwi_bval')
  1533. dwi_workflow.connect(preproc_dwi_node, 'dwi_mask', markvcid_freewater_node, 'in_dwi_mask')
  1534. fw_output_path = os.path.join(self.output_path, 'freewater', 'markvcid_freewater')
  1535. markvcid_freewater_node.inputs.output_dir = fw_output_path
  1536. markvcid_freewater_node.inputs.script_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), '..', 'external', 'MarkVCID2', 'scripts_FW_CONSORTIUM'))
  1537. markvcid_freewater_node.inputs.output_fw_img = os.path.join(os.path.dirname(fw_output_path), rename_bids_file(preproc_dwi_filename, {'model': 'tensor', 'param': 'freewater', 'desc': 'MarkVCID2'}, 'dwimap', '.nii.gz'))
  1538. dwi_workflow.connect(markvcid_freewater_node, 'out_fw', freewater_node, 'markvcid_fw_img')
  1539. dwi_workflow.connect(markvcid_freewater_node, 'out_fw', dwi_scalarmaps_output_node, 'markvcid2_fw_img')
  1540. elif os.path.exists(os.path.join(self.output_path, 'freewater', rename_bids_file(preproc_dwi_filename, {'model': 'tensor', 'param': 'freewater', 'desc': 'MarkVCID2'}, 'dwimap', '.nii.gz'))):
  1541. dwi_scalarmaps_output_node.inputs.markvcid2_fw_img = os.path.join(self.output_path, 'freewater', rename_bids_file(preproc_dwi_filename, {'model': 'tensor', 'param': 'freewater', 'desc': 'MarkVCID2'}, 'dwimap', '.nii.gz'))
  1542. if 'dti_freewater' in self.freewater:
  1543. dti_freewater_node = Node(FreeWaterTensor(), name='dti_freewater')
  1544. dwi_workflow.connect(preproc_dwi_node, 'preproc_dwi', dti_freewater_node, 'dwi_file')
  1545. dwi_workflow.connect(preproc_dwi_node, 'bvec', dti_freewater_node, 'bvec_file')
  1546. dwi_workflow.connect(preproc_dwi_node, 'bval', dti_freewater_node, 'bval_file')
  1547. dwi_workflow.connect(preproc_dwi_node, 'dwi_mask', dti_freewater_node, 'mask_file')
  1548. fw_output_path = os.path.join(self.output_path, 'freewater', 'tensor_model_freewater')
  1549. dti_freewater_node.inputs.output_dir = fw_output_path
  1550. dwi_workflow.connect(dti_freewater_node, 'freewater_file', freewater_node, 'dti_fw_img')
  1551. else:
  1552. freewater_node.inputs.dti_fw_img = os.path.join(self.output_path, 'tensor_model_freewater', 'dti_freewater.nii.gz')
  1553. # ===========================================
  1554. # PSMD
  1555. # ===========================================
  1556. if self.psmd:
  1557. print("[DWI Pipeline] PSMD analysis: True (using pipeline dtifit FA/MD)")
  1558. psmd_output_dir = os.path.join(self.output_path, 'psmd')
  1559. psmd_node = Node(PSMDCommandLine(), name="psmd_node")
  1560. dwi_workflow.connect(dti_fit_output_node, 'fa_img', psmd_node, 'fa_image')
  1561. dwi_workflow.connect(dti_fit_output_node, 'md_img', psmd_node, 'md_image')
  1562. psmd_node.inputs.mask_file = self.psmd_skeleton_mask
  1563. if self.psmd_exclude_seed_mask:
  1564. dwi_workflow.connect(seed_mask_to_dwi_node, 'out_file', psmd_node, 'lesion_mask')
  1565. psmd_node.inputs.output_dir = psmd_output_dir
  1566. ###########################
  1567. # Visual Pathway Analysis #
  1568. ###########################
  1569. if self.visual_pathway_analysis:
  1570. vpa_output_dir = os.path.join(self.output_path, 'visual_pathway_analysis')
  1571. os.makedirs(vpa_output_dir, exist_ok=True)
  1572. from cvdproc.pipelines.dmri.vp_project.tt_process_nipype import TTGZToTCK, TTGZToTDI, FlipTT
  1573. gather_gqi_fib_node = Node(IdentityInterface(fields=['gqi_fib']), name='gqi_fib_gather')
  1574. if self.preprocess_method == 'post_qsiprep':
  1575. qsirecon_output_dir = self.session._find_output('qsirecon-DSIStudio')
  1576. gqi_fib = os.path.join(qsirecon_output_dir, 'dwi', rename_bids_file(preproc_dwi_filename, {"desc": None, "model": "gqi"}, 'dwimap', '.fib.gz'))
  1577. if not os.path.exists(gqi_fib):
  1578. raise FileNotFoundError(f"GQI .fib.gz file not found in {qsirecon_output_dir}. Please run QSIRecon with DSIStudio reconstruction first.")
  1579. gather_gqi_fib_node.inputs.gqi_fib = gqi_fib
  1580. fix_header = False
  1581. else:
  1582. fix_header = True
  1583. print("[DWI Pipeline] Fix header in VP analysis: True")
  1584. if self.dsistudio_dti:
  1585. dwi_workflow.connect(dti_reconstruction_node, 'out_file', gather_gqi_fib_node, 'gqi_fib')
  1586. else:
  1587. dwi_workflow.connect(gqi_reconstruction_node, 'out_file', gather_gqi_fib_node, 'gqi_fib')
  1588. vp_flip_y = self.visual_pathway_analysis_flip_y
  1589. # 1. Prepare ROIs
  1590. from cvdproc.pipelines.dmri.vp_project.prepare_roi_nipype import PrepareVPROI
  1591. prepare_vp_roi_node = Node(PrepareVPROI(), name='prepare_vp_roi')
  1592. dwi_workflow.connect(inputnode, 'fs_output', prepare_vp_roi_node, 'fs_output_dir')
  1593. prepare_vp_roi_node.inputs.subject_id = self.subject.subject_id
  1594. prepare_vp_roi_node.inputs.session_id = self.session.session_id
  1595. dwi_workflow.connect(concat_fs_to_dwi_xfm_node, 'out_file', prepare_vp_roi_node, 'fs_to_dwi_xfm')
  1596. dwi_workflow.connect(preproc_dwi_node, 'b0', prepare_vp_roi_node, 'dwi_ref')
  1597. prepare_vp_roi_node.inputs.output_dir = anat_output_dir
  1598. prepare_vp_roi_node.inputs.space_entity = space_entity
  1599. # ROI for Meyer's Loop
  1600. roi4meyersloop_to_t1w_node = MapNode(MRIConvertApplyWarp(), name='roi4meyersloop_to_t1w', iterfield=['input_image', 'output_image'])
  1601. dwi_workflow.connect(inputnode, 'mni_to_t1w_warp', roi4meyersloop_to_t1w_node, 'warp_image')
  1602. roi4meyersloop_to_t1w_node.inputs.input_image = [
  1603. get_package_path('data', 'standard', 'MNI152', 'meyers_loop_roi', 'meyers_loop_roi_L.nii.gz'),
  1604. get_package_path('data', 'standard', 'MNI152', 'meyers_loop_roi', 'meyers_loop_roi_R.nii.gz')
  1605. ]
  1606. roi4meyersloop_to_t1w_node.inputs.output_image = [
  1607. os.path.join(anat_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_space-T1w_desc-roi4meyersloop_mask.nii.gz"),
  1608. os.path.join(anat_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_space-T1w_desc-roi4meyersloop_mask.nii.gz")
  1609. ]
  1610. roi4meyersloop_to_t1w_node.inputs.interp = 'nearest'
  1611. roi4meyersloop_to_dwi_node = MapNode(FLIRT(), name='roi4meyersloop_to_dwi', iterfield=['in_file', 'out_file'])
  1612. dwi_workflow.connect(invert_dwi_to_t1w_reg_node, 'out_file', roi4meyersloop_to_dwi_node, 'in_matrix_file')
  1613. dwi_workflow.connect(roi4meyersloop_to_t1w_node, 'output_image', roi4meyersloop_to_dwi_node, 'in_file')
  1614. roi4meyersloop_to_dwi_node.inputs.out_file = [
  1615. os.path.join(anat_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_space-{space_entity}_desc-roi4meyersloop_mask.nii.gz"),
  1616. os.path.join(anat_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_space-{space_entity}_desc-roi4meyersloop_mask.nii.gz")
  1617. ]
  1618. dwi_workflow.connect(preproc_dwi_node, 'b0', roi4meyersloop_to_dwi_node, 'reference')
  1619. roi4meyersloop_to_dwi_node.inputs.interp = "nearestneighbour"
  1620. roi4meyersloop_to_dwi_node.inputs.apply_xfm = True
  1621. # select left and right roi separately
  1622. fetch_roi4meyersloop_L = Node(Select(index=0), name="fetch_roi4meyersloop_L")
  1623. fetch_roi4meyersloop_R = Node(Select(index=1), name="fetch_roi4meyersloop_R")
  1624. dwi_workflow.connect(roi4meyersloop_to_dwi_node, "out_file", fetch_roi4meyersloop_L, "inlist")
  1625. dwi_workflow.connect(roi4meyersloop_to_dwi_node, "out_file", fetch_roi4meyersloop_R, "inlist")
  1626. # 2. ROI Tractography
  1627. raw_track_output_dir = os.path.join(vpa_output_dir, 'raw_tracts')
  1628. os.makedirs(raw_track_output_dir, exist_ok=True)
  1629. lh_or_track_node = Node(DSIstudioTracking(), name='lh_or_track')
  1630. dwi_workflow.connect(gather_gqi_fib_node, 'gqi_fib', lh_or_track_node, 'source')
  1631. lh_or_track_node.inputs.output = os.path.join(raw_track_output_dir, 'lh_OR.tt.gz')
  1632. lh_or_track_node.inputs.thread_count = 8
  1633. lh_or_track_node.inputs.tract_count = 2500
  1634. #lh_or_track_node.inputs.turning_angle = 45
  1635. lh_or_track_node.inputs.max_length = 200
  1636. lh_or_track_node.inputs.args = '--tip_iteration=0 --method=0'
  1637. dwi_workflow.connect(prepare_vp_roi_node, 'lh_lgn_dilx_roi', lh_or_track_node, 'seed')
  1638. dwi_workflow.connect(prepare_vp_roi_node, 'lh_v1_ext2_roi', lh_or_track_node, 'end')
  1639. #dwi_workflow.connect(fetch_roi4meyersloop_L, 'out', lh_or_track_node, 'ter')
  1640. rh_or_track_node = Node(DSIstudioTracking(), name='rh_or_track')
  1641. dwi_workflow.connect(gather_gqi_fib_node, 'gqi_fib', rh_or_track_node, 'source')
  1642. rh_or_track_node.inputs.output = os.path.join(raw_track_output_dir, 'rh_OR.tt.gz')
  1643. rh_or_track_node.inputs.thread_count = 8
  1644. rh_or_track_node.inputs.tract_count = 2500
  1645. #rh_or_track_node.inputs.turning_angle = 45
  1646. rh_or_track_node.inputs.max_length = 200
  1647. rh_or_track_node.inputs.args = '--tip_iteration=0 --method=0'
  1648. dwi_workflow.connect(prepare_vp_roi_node, 'rh_lgn_dilx_roi', rh_or_track_node, 'seed')
  1649. dwi_workflow.connect(prepare_vp_roi_node, 'rh_v1_ext2_roi', rh_or_track_node, 'end')
  1650. #dwi_workflow.connect(fetch_roi4meyersloop_R, 'out', rh_or_track_node, 'ter')
  1651. lh_ot_track_node = Node(DSIstudioTracking(), name='lh_ot_track')
  1652. dwi_workflow.connect(gather_gqi_fib_node, 'gqi_fib', lh_ot_track_node, 'source')
  1653. lh_ot_track_node.inputs.output = os.path.join(raw_track_output_dir, 'lh_OT.tt.gz')
  1654. lh_ot_track_node.inputs.thread_count = 8
  1655. lh_ot_track_node.inputs.tract_count = 2500
  1656. lh_ot_track_node.inputs.max_length = 60
  1657. lh_ot_track_node.inputs.args = '--tip_iteration=2 --method=1'
  1658. dwi_workflow.connect(prepare_vp_roi_node, 'optic_chiasm_dil1_roi', lh_ot_track_node, 'seed')
  1659. dwi_workflow.connect(prepare_vp_roi_node, 'lh_lgn_dil1_roi', lh_ot_track_node, 'end')
  1660. rh_ot_track_node = Node(DSIstudioTracking(), name='rh_ot_track')
  1661. dwi_workflow.connect(gather_gqi_fib_node, 'gqi_fib', rh_ot_track_node, 'source')
  1662. rh_ot_track_node.inputs.output = os.path.join(raw_track_output_dir, 'rh_OT.tt.gz')
  1663. rh_ot_track_node.inputs.thread_count = 8
  1664. rh_ot_track_node.inputs.tract_count = 2500
  1665. rh_ot_track_node.inputs.max_length = 60
  1666. rh_ot_track_node.inputs.args = '--tip_iteration=2 --method=1'
  1667. dwi_workflow.connect(prepare_vp_roi_node, 'optic_chiasm_dil1_roi', rh_ot_track_node, 'seed')
  1668. dwi_workflow.connect(prepare_vp_roi_node, 'rh_lgn_dil1_roi', rh_ot_track_node, 'end')
  1669. # Experimental: Meyer's Loop
  1670. lr_ml_track_node = Node(DSIstudioTracking(), name='lr_meyersloop_track')
  1671. dwi_workflow.connect(gather_gqi_fib_node, 'gqi_fib', lr_ml_track_node, 'source')
  1672. lr_ml_track_node.inputs.output = os.path.join(raw_track_output_dir, 'lr_MeyersLoop.tt.gz')
  1673. lr_ml_track_node.inputs.thread_count = 8
  1674. lr_ml_track_node.inputs.tract_count = 2500
  1675. lr_ml_track_node.inputs.turning_angle = 45
  1676. lr_ml_track_node.inputs.args = '--tip_iteration=2 --method=1'
  1677. dwi_workflow.connect(fetch_roi4meyersloop_L, 'out', lr_ml_track_node, 'seed')
  1678. dwi_workflow.connect(prepare_vp_roi_node, 'lh_v1_ext2_roi', lr_ml_track_node, 'end')
  1679. rh_ml_track_node = Node(DSIstudioTracking(), name='rh_meyersloop_track')
  1680. dwi_workflow.connect(gather_gqi_fib_node, 'gqi_fib', rh_ml_track_node, 'source')
  1681. rh_ml_track_node.inputs.output = os.path.join(raw_track_output_dir, 'rh_MeyersLoop.tt.gz')
  1682. rh_ml_track_node.inputs.thread_count = 8
  1683. rh_ml_track_node.inputs.tract_count = 2500
  1684. rh_ml_track_node.inputs.turning_angle = 45
  1685. rh_ml_track_node.inputs.args = '--tip_iteration=2 --method=1'
  1686. dwi_workflow.connect(fetch_roi4meyersloop_R, 'out', rh_ml_track_node, 'seed')
  1687. dwi_workflow.connect(prepare_vp_roi_node, 'rh_v1_ext2_roi', rh_ml_track_node, 'end')
  1688. # 3. Refine tracts with ROIs
  1689. from cvdproc.pipelines.dmri.vp_project.refine_vp_nipype import RefineVP
  1690. vp_refine_node = Node(RefineVP(), name='vp_refine')
  1691. # LPS -> RAS if fix_header
  1692. if fix_header:
  1693. gather_raw_track_node = Node(Merge(6), name='gather_raw_track')
  1694. dwi_workflow.connect(lh_or_track_node, 'output', gather_raw_track_node, 'in1')
  1695. dwi_workflow.connect(lh_ot_track_node, 'output', gather_raw_track_node, 'in2')
  1696. dwi_workflow.connect(rh_or_track_node, 'output', gather_raw_track_node, 'in3')
  1697. dwi_workflow.connect(rh_ot_track_node, 'output', gather_raw_track_node, 'in4')
  1698. dwi_workflow.connect(lr_ml_track_node, 'output', gather_raw_track_node, 'in5')
  1699. dwi_workflow.connect(rh_ml_track_node, 'output', gather_raw_track_node, 'in6')
  1700. transform_tt_space_node = MapNode(FlipTT(), name='transform_tt_space', iterfield=['in_tt', 'out_tt'])
  1701. dwi_workflow.connect(gather_raw_track_node, 'out', transform_tt_space_node, 'in_tt')
  1702. transform_tt_space_node.inputs.out_tt = [
  1703. os.path.join(raw_track_output_dir, 'lh_OR_flipy.tt.gz'),
  1704. os.path.join(raw_track_output_dir, 'lh_OT_flipy.tt.gz'),
  1705. os.path.join(raw_track_output_dir, 'rh_OR_flipy.tt.gz'),
  1706. os.path.join(raw_track_output_dir, 'rh_OT_flipy.tt.gz'),
  1707. os.path.join(raw_track_output_dir, 'lr_ML_flipy.tt.gz'),
  1708. os.path.join(raw_track_output_dir, 'rh_ML_flipy.tt.gz'),
  1709. ]
  1710. transform_tt_space_node.inputs.flip_x_axis = False
  1711. transform_tt_space_node.inputs.flip_y_axis = True
  1712. transform_tt_space_node.inputs.flip_z_axis = False
  1713. select_lh_or_node = Node(Select(index=[0]), name="select_lh_or_ras")
  1714. select_lh_ot_node = Node(Select(index=[1]), name="select_lh_ot_ras")
  1715. select_rh_or_node = Node(Select(index=[2]), name="select_rh_or_ras")
  1716. select_rh_ot_node = Node(Select(index=[3]), name="select_rh_ot_ras")
  1717. select_lh_ml_node = Node(Select(index=[4]), name="select_lh_ml_ras")
  1718. select_rh_ml_node = Node(Select(index=[5]), name="select_rh_ml_ras")
  1719. dwi_workflow.connect(transform_tt_space_node, "out_tt", select_lh_or_node, "inlist")
  1720. dwi_workflow.connect(transform_tt_space_node, "out_tt", select_lh_ot_node, "inlist")
  1721. dwi_workflow.connect(transform_tt_space_node, "out_tt", select_rh_or_node, "inlist")
  1722. dwi_workflow.connect(transform_tt_space_node, "out_tt", select_rh_ot_node, "inlist")
  1723. dwi_workflow.connect(transform_tt_space_node, "out_tt", select_lh_ml_node, "inlist")
  1724. dwi_workflow.connect(transform_tt_space_node, "out_tt", select_rh_ml_node, "inlist")
  1725. dwi_workflow.connect(select_lh_or_node, "out", vp_refine_node, "lh_or")
  1726. dwi_workflow.connect(select_lh_ot_node, "out", vp_refine_node, "lh_ot")
  1727. dwi_workflow.connect(select_rh_or_node, "out", vp_refine_node, "rh_or")
  1728. dwi_workflow.connect(select_rh_ot_node, "out", vp_refine_node, "rh_ot")
  1729. dwi_workflow.connect(select_lh_ml_node, "out", vp_refine_node, "lh_ml")
  1730. dwi_workflow.connect(select_rh_ml_node, "out", vp_refine_node, "rh_ml")
  1731. else:
  1732. dwi_workflow.connect(lh_or_track_node, 'output', vp_refine_node, 'lh_or')
  1733. dwi_workflow.connect(lh_ot_track_node, 'output', vp_refine_node, 'lh_ot')
  1734. dwi_workflow.connect(rh_or_track_node, 'output', vp_refine_node, 'rh_or')
  1735. dwi_workflow.connect(rh_ot_track_node, 'output', vp_refine_node, 'rh_ot')
  1736. dwi_workflow.connect(lr_ml_track_node, 'output', vp_refine_node, 'lh_ml')
  1737. dwi_workflow.connect(rh_ml_track_node, 'output', vp_refine_node, 'rh_ml')
  1738. dwi_workflow.connect(prepare_vp_roi_node, 'optic_chiasm_dil1_roi', vp_refine_node, 'cho_roi')
  1739. dwi_workflow.connect(prepare_vp_roi_node, 'lh_lgn_dil1_roi', vp_refine_node, 'lh_lgn_roi')
  1740. dwi_workflow.connect(prepare_vp_roi_node, 'lh_lgn_dilx_roi', vp_refine_node, 'lh_lgn_dia_x_roi')
  1741. dwi_workflow.connect(prepare_vp_roi_node, 'lh_lgn_extendpart_roi', vp_refine_node, 'lh_lgn_extendpart_roi')
  1742. dwi_workflow.connect(prepare_vp_roi_node, 'rh_lgn_dil1_roi', vp_refine_node, 'rh_lgn_roi')
  1743. dwi_workflow.connect(prepare_vp_roi_node, 'rh_lgn_dilx_roi', vp_refine_node, 'rh_lgn_dia_x_roi')
  1744. dwi_workflow.connect(prepare_vp_roi_node, 'rh_lgn_extendpart_roi', vp_refine_node, 'rh_lgn_extendpart_roi')
  1745. dwi_workflow.connect(prepare_vp_roi_node, 'lh_v1_ext2_roi', vp_refine_node, 'lh_v1_roi')
  1746. dwi_workflow.connect(prepare_vp_roi_node, 'rh_v1_ext2_roi', vp_refine_node, 'rh_v1_roi')
  1747. dwi_workflow.connect(fetch_roi4meyersloop_L, 'out', vp_refine_node, 'lh_meyersloop_roi')
  1748. dwi_workflow.connect(fetch_roi4meyersloop_R, 'out', vp_refine_node, 'rh_meyersloop_roi')
  1749. vp_refine_node.inputs.output_dir = vpa_output_dir
  1750. vp_refine_node.inputs.output_lh_ot = rename_bids_file(preproc_dwi_filename, {'bundle': 'OT', 'hemi': 'L', 'desc': 'voxelspace'}, 'streamlines', '.tt.gz')
  1751. vp_refine_node.inputs.output_rh_ot = rename_bids_file(preproc_dwi_filename, {'bundle': 'OT', 'hemi': 'R', 'desc': 'voxelspace'}, 'streamlines', '.tt.gz')
  1752. vp_refine_node.inputs.output_lh_or = rename_bids_file(preproc_dwi_filename, {'bundle': 'OR', 'hemi': 'L', 'desc': 'voxelspace'}, 'streamlines', '.tt.gz')
  1753. vp_refine_node.inputs.output_rh_or = rename_bids_file(preproc_dwi_filename, {'bundle': 'OR', 'hemi': 'R', 'desc': 'voxelspace'}, 'streamlines', '.tt.gz')
  1754. vp_refine_node.inputs.output_lh_ml = rename_bids_file(preproc_dwi_filename, {'bundle': 'ML', 'hemi': 'L', 'desc': 'voxelspace'}, 'streamlines', '.tt.gz')
  1755. vp_refine_node.inputs.output_rh_ml = rename_bids_file(preproc_dwi_filename, {'bundle': 'ML', 'hemi': 'R', 'desc': 'voxelspace'}, 'streamlines', '.tt.gz')
  1756. if fix_header:
  1757. gather_fliped_refined_tracts = Node(Merge(6), name='gather_fliped_refined_tracts')
  1758. dwi_workflow.connect(vp_refine_node, 'refined_lh_ot', gather_fliped_refined_tracts, 'in1')
  1759. dwi_workflow.connect(vp_refine_node, 'refined_rh_ot', gather_fliped_refined_tracts, 'in2')
  1760. dwi_workflow.connect(vp_refine_node, 'refined_lh_or', gather_fliped_refined_tracts, 'in3')
  1761. dwi_workflow.connect(vp_refine_node, 'refined_rh_or', gather_fliped_refined_tracts, 'in4')
  1762. dwi_workflow.connect(vp_refine_node, 'refined_lh_ml', gather_fliped_refined_tracts, 'in5')
  1763. dwi_workflow.connect(vp_refine_node, 'refined_rh_ml', gather_fliped_refined_tracts, 'in6')
  1764. re_flip_node = MapNode(FlipTT(), name='re_flip_tt', iterfield=['in_tt', 'out_tt'])
  1765. dwi_workflow.connect(gather_fliped_refined_tracts, 'out', re_flip_node, 'in_tt')
  1766. re_flip_node.inputs.out_tt = [
  1767. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'bundle': 'OT', 'hemi': 'L', 'desc': 'voxelspace'}, 'streamlines', '.tt.gz')),
  1768. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'bundle': 'OT', 'hemi': 'R', 'desc': 'voxelspace'}, 'streamlines', '.tt.gz')),
  1769. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'bundle': 'OR', 'hemi': 'L', 'desc': 'voxelspace'}, 'streamlines', '.tt.gz')),
  1770. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'bundle': 'OR', 'hemi': 'R', 'desc': 'voxelspace'}, 'streamlines', '.tt.gz')),
  1771. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'bundle': 'ML', 'hemi': 'L', 'desc': 'voxelspace'}, 'streamlines', '.tt.gz')),
  1772. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'bundle': 'ML', 'hemi': 'R', 'desc': 'voxelspace'}, 'streamlines', '.tt.gz')),
  1773. ]
  1774. re_flip_node.inputs.flip_x_axis = False
  1775. re_flip_node.inputs.flip_y_axis = True
  1776. re_flip_node.inputs.flip_z_axis = False
  1777. gather_refined_tracts = Node(IdentityInterface(fields=['out']), name='gather_refined_tracts')
  1778. dwi_workflow.connect(re_flip_node, 'out_tt', gather_refined_tracts, 'out')
  1779. else:
  1780. gather_refined_tracts = Node(Merge(6), name='gather_refined_tracts')
  1781. dwi_workflow.connect(vp_refine_node, 'refined_lh_ot', gather_refined_tracts, 'in1')
  1782. dwi_workflow.connect(vp_refine_node, 'refined_rh_ot', gather_refined_tracts, 'in2')
  1783. dwi_workflow.connect(vp_refine_node, 'refined_lh_or', gather_refined_tracts, 'in3')
  1784. dwi_workflow.connect(vp_refine_node, 'refined_rh_or', gather_refined_tracts, 'in4')
  1785. dwi_workflow.connect(vp_refine_node, 'refined_lh_ml', gather_refined_tracts, 'in5')
  1786. dwi_workflow.connect(vp_refine_node, 'refined_rh_ml', gather_refined_tracts, 'in6')
  1787. # tract stats
  1788. tract_stats_node = MapNode(TractStatsInterface(), name='tract_stats', iterfield=['tract', 'output_txt'])
  1789. dwi_workflow.connect(gather_refined_tracts, 'out', tract_stats_node, 'tract')
  1790. dwi_workflow.connect(gather_gqi_fib_node, 'gqi_fib', tract_stats_node, 'source')
  1791. tract_stats_node.inputs.output_txt = [
  1792. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'OT', 'hemi': 'L'}, 'stats', '.txt')),
  1793. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'OT', 'hemi': 'R'}, 'stats', '.txt')),
  1794. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'OR', 'hemi': 'L'}, 'stats', '.txt')),
  1795. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'OR', 'hemi': 'R'}, 'stats', '.txt')),
  1796. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'ML', 'hemi': 'L'}, 'stats', '.txt')),
  1797. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'ML', 'hemi': 'R'}, 'stats', '.txt'))
  1798. ]
  1799. tract_stats_node.inputs.export = 'stat'
  1800. # tt.gz to tck
  1801. tt_to_tck_node = MapNode(TTGZToTCK(), name='tt_to_tck', iterfield=['tt_gz_path', 'out_tck'])
  1802. tt_to_tck_node.synchronize = True
  1803. dwi_workflow.connect(gather_refined_tracts, 'out', tt_to_tck_node, 'tt_gz_path')
  1804. dwi_workflow.connect(preproc_dwi_node, 'b0', tt_to_tck_node, 'ref_nii')
  1805. tt_to_tck_node.inputs.out_tck = [
  1806. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'OT', 'hemi': 'L'}, 'streamlines', '.tck')),
  1807. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'OT', 'hemi': 'R'}, 'streamlines', '.tck')),
  1808. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'OR', 'hemi': 'L'}, 'streamlines', '.tck')),
  1809. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'OR', 'hemi': 'R'}, 'streamlines', '.tck')),
  1810. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'ML', 'hemi': 'L'}, 'streamlines', '.tck')),
  1811. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'ML', 'hemi': 'R'}, 'streamlines', '.tck'))
  1812. ]
  1813. tt_to_tck_node.inputs.flip_y_axis = vp_flip_y
  1814. # tt.gz to TDI image
  1815. tt_to_tdi_node = MapNode(TTGZToTDI(), name='tt_to_tdi', iterfield=['tt_gz_path', 'out_tdi'])
  1816. tt_to_tdi_node.synchronize = True
  1817. dwi_workflow.connect(gather_refined_tracts, 'out', tt_to_tdi_node, 'tt_gz_path')
  1818. dwi_workflow.connect(preproc_dwi_node, 'b0', tt_to_tdi_node, 'ref_nii')
  1819. tt_to_tdi_node.inputs.out_tdi = [
  1820. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'OT', 'hemi': 'L'}, 'tdi', '.nii.gz')),
  1821. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'OT', 'hemi': 'R'}, 'tdi', '.nii.gz')),
  1822. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'OR', 'hemi': 'L'}, 'tdi', '.nii.gz')),
  1823. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'OR', 'hemi': 'R'}, 'tdi', '.nii.gz')),
  1824. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'ML', 'hemi': 'L'}, 'tdi', '.nii.gz')),
  1825. os.path.join(vpa_output_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': 'ML', 'hemi': 'R'}, 'tdi', '.nii.gz'))
  1826. ]
  1827. tt_to_tdi_node.inputs.flip_y_axis = vp_flip_y
  1828. # gather OR and OT tck files for QC
  1829. select_or_ot = Node(Select(index=[0, 1, 2, 3]), name="select_or_ot")
  1830. dwi_workflow.connect(tt_to_tck_node, "out_tck", select_or_ot, "inlist")
  1831. or_ot_qc_node = Node(PlotTckOnSliceCmd(), name="or_ot_qc")
  1832. or_ot_qc_node.inputs.script_path = get_package_path("pipelines", "dmri", "visualization", "plot_tck_cli.py")
  1833. dwi_workflow.connect(select_or_ot, "out", or_ot_qc_node, "tck_files")
  1834. dwi_workflow.connect(preproc_dwi_node, "b0", or_ot_qc_node, "ref_nii")
  1835. or_ot_qc_node.inputs.out_png = os.path.join(vpa_output_dir, "QC.png")
  1836. # Detect existing visual pathway tracking outputs for value extraction
  1837. # (allows re-extracting values from existing tracts without re-running tracking)
  1838. vpa_check_dir = os.path.join(self.output_path, 'visual_pathway_analysis')
  1839. vp_tck_files = [
  1840. os.path.join(vpa_check_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': b, 'hemi': h}, 'streamlines', '.tck'))
  1841. for h, b in [('L', 'OT'), ('R', 'OT'), ('L', 'OR'), ('R', 'OR'), ('L', 'ML'), ('R', 'ML')]
  1842. ]
  1843. vp_tdi_files = [
  1844. os.path.join(vpa_check_dir, rename_bids_file(preproc_dwi_filename, {'desc': None, 'bundle': b, 'hemi': h}, 'tdi', '.nii.gz'))
  1845. for h, b in [('L', 'OT'), ('R', 'OT'), ('L', 'OR'), ('R', 'OR'), ('L', 'ML'), ('R', 'ML')]
  1846. ]
  1847. vp_tck_exist = all(os.path.exists(f) for f in vp_tck_files)
  1848. vp_tdi_exist = all(os.path.exists(f) for f in vp_tdi_files)
  1849. if self.visual_pathway_analysis:
  1850. vp_extract = True
  1851. elif vp_tck_exist and vp_tdi_exist:
  1852. print("[DWI Pipeline] Visual pathway tracking outputs found. Enabling value extraction from existing tracts.")
  1853. vp_extract = True
  1854. else:
  1855. vp_extract = False
  1856. if vp_extract and not self.visual_pathway_analysis:
  1857. # Use existing tracking outputs for value extraction
  1858. tt_to_tck_node = Node(IdentityInterface(fields=['out_tck']), name='tt_to_tck_stub')
  1859. tt_to_tck_node.inputs.out_tck = vp_tck_files
  1860. tt_to_tdi_node = Node(IdentityInterface(fields=['out_tdi']), name='tt_to_tdi_stub')
  1861. tt_to_tdi_node.inputs.out_tdi = vp_tdi_files
  1862. # Find existing WMH probmap in DWI space
  1863. anat_dir = os.path.join(self.output_path, 'anat')
  1864. wmh_probmap_dwi_glob = os.path.join(anat_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-{space_entity}_*probmap.nii.gz")
  1865. wmh_probmap_dwi_candidates = glob.glob(wmh_probmap_dwi_glob)
  1866. wmh_probmap_dwi_file = wmh_probmap_dwi_candidates[0] if len(wmh_probmap_dwi_candidates) == 1 else ''
  1867. wmhprobmap_to_dwi_node = Node(IdentityInterface(fields=['out_file']), name='wmhprobmap_stub')
  1868. wmhprobmap_to_dwi_node.inputs.out_file = wmh_probmap_dwi_file
  1869. #########################
  1870. # Calculate DWI metrics #
  1871. #########################
  1872. if self.calculate_dwi_metrics:
  1873. dwi_metrics_output_dir = os.path.join(self.output_path, 'dwi_metrics_stats')
  1874. os.makedirs(dwi_metrics_output_dir, exist_ok=True)
  1875. dwi_metrics_node = Node(Merge(15), name='dwi_metrics_node')
  1876. dwi_workflow.connect(dwi_scalarmaps_output_node, 'fa_img', dwi_metrics_node, 'in1')
  1877. dwi_workflow.connect(dwi_scalarmaps_output_node, 'md_img', dwi_metrics_node, 'in2')
  1878. dwi_workflow.connect(dwi_scalarmaps_output_node, 'markvcid2_fw_img', dwi_metrics_node, 'in3')
  1879. dwi_workflow.connect(dwi_scalarmaps_output_node, 'ad_img', dwi_metrics_node, 'in4')
  1880. dwi_workflow.connect(dwi_scalarmaps_output_node, 'rd_img', dwi_metrics_node, 'in5')
  1881. dwi_workflow.connect(dwi_scalarmaps_output_node, 'odi_img', dwi_metrics_node, 'in6')
  1882. dwi_workflow.connect(dwi_scalarmaps_output_node, 'icvf_img', dwi_metrics_node, 'in7')
  1883. dwi_workflow.connect(dwi_scalarmaps_output_node, 'isovf_img', dwi_metrics_node, 'in8')
  1884. dwi_workflow.connect(dwi_scalarmaps_output_node, 'gqi_gfa_img', dwi_metrics_node, 'in9')
  1885. dwi_workflow.connect(dwi_scalarmaps_output_node, 'gqi_iso_img', dwi_metrics_node, 'in10')
  1886. dwi_workflow.connect(dwi_scalarmaps_output_node, 'gqi_qa_img', dwi_metrics_node, 'in11')
  1887. dwi_workflow.connect(dwi_scalarmaps_output_node, 'chidia_img', dwi_metrics_node, 'in12')
  1888. dwi_workflow.connect(dwi_scalarmaps_output_node, 'mk_img', dwi_metrics_node, 'in13')
  1889. dwi_workflow.connect(dwi_scalarmaps_output_node, 'ak_img', dwi_metrics_node, 'in14')
  1890. dwi_workflow.connect(dwi_scalarmaps_output_node, 'rk_img', dwi_metrics_node, 'in15')
  1891. exist_dwi_metrics_node = Node(FilterExisting(), name='exist_dwi_metrics')
  1892. dwi_workflow.connect(dwi_metrics_node, 'out', exist_dwi_metrics_node, 'input_file_list')
  1893. if vp_extract:
  1894. # Add WMH probmap to scalar maps
  1895. vp_dwi_metrics_node = Node(Merge(2), name='vp_dwi_metrics_node')
  1896. dwi_workflow.connect(dwi_metrics_node, 'out', vp_dwi_metrics_node, 'in1')
  1897. dwi_workflow.connect(wmhprobmap_to_dwi_node, 'out_file', vp_dwi_metrics_node, 'in2')
  1898. # -------------------------------------- #
  1899. # Calculate DWI metrics using .tck files #
  1900. # -------------------------------------- #
  1901. if 'mrtrix3' in self.tractography:
  1902. seed_tract_scalar_maps_node = Node(TckSampleMultiScalarBundle(), name='seed_tract_scalar_maps')
  1903. dwi_workflow.connect(seed_based_track_node, 'seed_based_track', seed_tract_scalar_maps_node, 'tck_file')
  1904. dwi_workflow.connect(dwi_metrics_node, 'out', seed_tract_scalar_maps_node, 'scalar_files')
  1905. seed_tract_scalar_maps_node.inputs.scalar_names = ["FA", "MD", "FW (MarkVCID2)", "AD", "RD", "ODI", "ICVF", "ISOVF", "GQI_GFA", "GQI_ISO", "GQI_QA", "CHIDIA", "MK", "AK", "RK"]
  1906. seed_tract_scalar_maps_node.inputs.stat_tck = 'mean'
  1907. seed_tract_scalar_maps_node.inputs.output_csv = os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_label-ROIseedTract_desc-mean_dwimap.csv")
  1908. if vp_extract:
  1909. vp_tract_scalar_maps_node = MapNode(TckSampleMultiScalarProfile(), name='vp_tract_scalar_maps', iterfield=['tck_file', 'output_csv'])
  1910. dwi_workflow.connect(tt_to_tck_node, 'out_tck', vp_tract_scalar_maps_node, 'tck_file')
  1911. dwi_workflow.connect(vp_dwi_metrics_node, 'out', vp_tract_scalar_maps_node, 'scalar_files')
  1912. vp_tract_scalar_maps_node.inputs.scalar_names = ["FA", "MD", "FW (MarkVCID2)", "AD", "RD", "ODI", "ICVF", "ISOVF", "GQI_GFA", "GQI_ISO", "GQI_QA", "CHIDIA", "WMHprobmap"]
  1913. vp_tract_scalar_maps_node.inputs.output_csv = [
  1914. os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_label-OT_desc-alongtract_dwimap.csv"),
  1915. os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_label-OT_desc-alongtract_dwimap.csv"),
  1916. os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_label-OR_desc-alongtract_dwimap.csv"),
  1917. os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_label-OR_desc-alongtract_dwimap.csv"),
  1918. os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_label-ML_desc-alongtract_dwimap.csv"),
  1919. os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_label-ML_desc-alongtract_dwimap.csv")
  1920. ]
  1921. # ------------------------------------------------ #
  1922. # Calculate DWI metrics using volumetric masks #
  1923. # ------------------------------------------------ #
  1924. # 1. Seed mask
  1925. if seed_mask != '' and seed_mask is not None:
  1926. calc_seedmask_node = Node(CalculateScalarMaps(), name="scalar_maps_for_seedmask")
  1927. calc_seedmask_node.inputs.colnames = ["FA", "MD", "FW (MarkVCID2)", "AD", "RD", "ODI", "ICVF", "ISOVF", "GQI_GFA", "GQI_ISO", "GQI_QA", "CHIDIA", "MK", "AK", "RK"]
  1928. calc_seedmask_node.inputs.roi_label = 1
  1929. calc_seedmask_node.inputs.output_csv = os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_label-seedmask_desc-mean_dwimap.csv")
  1930. dwi_workflow.connect(seed_mask_to_dwi_node, 'out_file', calc_seedmask_node, 'mask_file')
  1931. dwi_workflow.connect(dwi_metrics_node, 'out', calc_seedmask_node, 'data_files')
  1932. calc_seedmask_node.inputs.ignore_background = False
  1933. # 2. WMH mask
  1934. if wmh_mask_file != '' and wmh_mask_file is not None:
  1935. scalar_maps_for_wmhmask_node = Node(CalculateScalarMaps(), name='scalar_maps_for_wmhmask')
  1936. scalar_maps_for_wmhmask_node.inputs.colnames = ["FA", "MD", "FW (MarkVCID2)", "AD", "RD", "ODI", "ICVF", "ISOVF", "GQI_GFA", "GQI_ISO", "GQI_QA", "CHIDIA", "MK", "AK", "RK"]
  1937. scalar_maps_for_wmhmask_node.inputs.roi_label = 1
  1938. scalar_maps_for_wmhmask_node.inputs.output_csv = os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_label-WMH_desc-mean_dwimap.csv")
  1939. dwi_workflow.connect(wmh_mask_to_dwi_node, 'out_file', scalar_maps_for_wmhmask_node, 'mask_file')
  1940. dwi_workflow.connect(dwi_metrics_node, 'out', scalar_maps_for_wmhmask_node, 'data_files')
  1941. scalar_maps_for_wmhmask_node.inputs.ignore_background = False
  1942. # 3. NAWM dseg (1=left, 2=right): per-ROI and combined
  1943. if self.dwi_t1w_register and nawm_dseg is not None:
  1944. scalar_maps_for_nawm_perroi_node = Node(CalculateScalarMaps(), name='scalar_maps_for_nawm_perroi')
  1945. scalar_maps_for_nawm_perroi_node.inputs.colnames = ["FA", "MD", "FW (MarkVCID2)", "AD", "RD", "ODI", "ICVF", "ISOVF", "GQI_GFA", "GQI_ISO", "GQI_QA", "CHIDIA", "MK", "AK", "RK"]
  1946. scalar_maps_for_nawm_perroi_node.inputs.roi_labels = [1, 2]
  1947. scalar_maps_for_nawm_perroi_node.inputs.statistic = "mean"
  1948. scalar_maps_for_nawm_perroi_node.inputs.output_csv = os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_label-NAWM_desc-perROI_dwimap.csv")
  1949. dwi_workflow.connect(merge_nawm_dseg_node, 'out_file', scalar_maps_for_nawm_perroi_node, 'mask_file')
  1950. dwi_workflow.connect(dwi_metrics_node, 'out', scalar_maps_for_nawm_perroi_node, 'data_files')
  1951. scalar_maps_for_nawm_perroi_node.inputs.ignore_background = False
  1952. scalar_maps_for_nawm_combined_node = Node(CalculateScalarMaps(), name='scalar_maps_for_nawm_combined')
  1953. scalar_maps_for_nawm_combined_node.inputs.colnames = ["FA", "MD", "FW (MarkVCID2)", "AD", "RD", "ODI", "ICVF", "ISOVF", "GQI_GFA", "GQI_ISO", "GQI_QA", "CHIDIA", "MK", "AK", "RK"]
  1954. scalar_maps_for_nawm_combined_node.inputs.roi_labels = [1, 2]
  1955. scalar_maps_for_nawm_combined_node.inputs.combine_rois = True
  1956. scalar_maps_for_nawm_combined_node.inputs.statistic = "mean"
  1957. scalar_maps_for_nawm_combined_node.inputs.output_csv = os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_label-NAWM_desc-mean_dwimap.csv")
  1958. dwi_workflow.connect(merge_nawm_dseg_node, 'out_file', scalar_maps_for_nawm_combined_node, 'mask_file')
  1959. dwi_workflow.connect(dwi_metrics_node, 'out', scalar_maps_for_nawm_combined_node, 'data_files')
  1960. scalar_maps_for_nawm_combined_node.inputs.ignore_background = False
  1961. # 4. NAWM dseg without tract mask (combined)
  1962. if self.dwi_t1w_register and nawm_dseg is not None and 'mrtrix3' in self.tractography:
  1963. scalar_maps_for_nawm_without_tract_node = Node(CalculateScalarMaps(), name='scalar_maps_for_nawm_without_tract')
  1964. scalar_maps_for_nawm_without_tract_node.inputs.colnames = ["FA", "MD", "FW (MarkVCID2)", "AD", "RD", "ODI", "ICVF", "ISOVF", "GQI_GFA", "GQI_ISO", "GQI_QA", "CHIDIA", "MK", "AK", "RK"]
  1965. scalar_maps_for_nawm_without_tract_node.inputs.roi_labels = [1, 2]
  1966. scalar_maps_for_nawm_without_tract_node.inputs.combine_rois = True
  1967. scalar_maps_for_nawm_without_tract_node.inputs.statistic = "mean"
  1968. scalar_maps_for_nawm_without_tract_node.inputs.output_csv = os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_label-NAWMwithoutTract_desc-mean_dwimap.csv")
  1969. dwi_workflow.connect(remove_tract_region_node, 'out_wm_mask', scalar_maps_for_nawm_without_tract_node, 'mask_file')
  1970. dwi_workflow.connect(dwi_metrics_node, 'out', scalar_maps_for_nawm_without_tract_node, 'data_files')
  1971. scalar_maps_for_nawm_without_tract_node.inputs.ignore_background = False
  1972. if self.dwi_t1w_register and fs_output_process:
  1973. # 5. wmparc from fs
  1974. scalar_maps_for_wmparc_node = Node(CalculateScalarMaps(), name='scalar_maps_for_wmparc')
  1975. scalar_maps_for_wmparc_node.inputs.colnames = ["FA", "MD", "FW (MarkVCID2)", "AD", "RD", "ODI", "ICVF", "ISOVF", "GQI_GFA", "GQI_ISO", "GQI_QA", "CHIDIA", "MK", "AK", "RK"]
  1976. scalar_maps_for_wmparc_node.inputs.output_csv = os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_label-wmparc_desc-median_dwimap.csv")
  1977. scalar_maps_for_wmparc_node.inputs.statistic = "median" # From Dr.Maria's code
  1978. dwi_workflow.connect(final_wmparc_node, 'final_wmparc', scalar_maps_for_wmparc_node, 'mask_file')
  1979. dwi_workflow.connect(dwi_metrics_node, 'out', scalar_maps_for_wmparc_node, 'data_files')
  1980. scalar_maps_for_wmparc_node.inputs.ignore_background = False
  1981. # 6. aparc+aseg from fs
  1982. scalar_maps_for_aparc_aseg_node = Node(CalculateScalarMaps(), name='scalar_maps_for_aparc_aseg')
  1983. scalar_maps_for_aparc_aseg_node.inputs.colnames = ["FA", "MD", "FW (MarkVCID2)", "AD", "RD", "ODI", "ICVF", "ISOVF", "GQI_GFA", "GQI_ISO", "GQI_QA", "CHIDIA", "MK", "AK", "RK"]
  1984. scalar_maps_for_aparc_aseg_node.inputs.output_csv = os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_label-aparcaseg_desc-mean_dwimap.csv")
  1985. scalar_maps_for_aparc_aseg_node.inputs.statistic = "mean"
  1986. dwi_workflow.connect(fs_aparcaseg_to_dwi_node, 'out_file', scalar_maps_for_aparc_aseg_node, 'mask_file')
  1987. dwi_workflow.connect(dwi_metrics_node, 'out', scalar_maps_for_aparc_aseg_node, 'data_files')
  1988. scalar_maps_for_aparc_aseg_node.inputs.ignore_background = False
  1989. # 7. chpseg
  1990. if chpseg_process:
  1991. scalar_maps_for_chpseg_node = Node(CalculateScalarMaps(), name='scalar_maps_for_chpseg')
  1992. scalar_maps_for_chpseg_node.inputs.colnames = ["FA", "MD", "FW (MarkVCID2)", "AD", "RD", "ODI", "ICVF", "ISOVF", "GQI_GFA", "GQI_ISO", "GQI_QA", "CHIDIA", "MK", "AK", "RK"]
  1993. scalar_maps_for_chpseg_node.inputs.output_csv = os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_label-chpseg_desc-mean_dwimap.csv")
  1994. scalar_maps_for_chpseg_node.inputs.statistic = "mean"
  1995. dwi_workflow.connect(chpseg_mask_to_dwi_node, 'out_file', scalar_maps_for_chpseg_node, 'mask_file')
  1996. dwi_workflow.connect(dwi_metrics_node, 'out', scalar_maps_for_chpseg_node, 'data_files')
  1997. scalar_maps_for_chpseg_node.inputs.ignore_background = False
  1998. # 8. JHU atlas
  1999. if self.dwi_t1w_register and mni_to_t1w_warp != "":
  2000. jhu_atlas_to_t1w_node = Node(MRIConvertApplyWarp(), name='jhu_atlas_to_t1w')
  2001. dwi_workflow.connect(inputnode, 'mni_to_t1w_warp', jhu_atlas_to_t1w_node, 'warp_image')
  2002. jhu_atlas_to_t1w_node.inputs.input_image = get_package_path('data', 'standard', 'JHU', 'JHU-ICBM-labels-1mm.nii.gz')
  2003. jhu_atlas_to_t1w_node.inputs.output_image = os.path.join(anat_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-T1w_seg-JHU_dseg.nii.gz")
  2004. jhu_atlas_to_t1w_node.inputs.interp = 'nearest'
  2005. jhu_atlas_to_dwi_node = Node(FLIRT(), name='jhu_atlas_to_dwi')
  2006. dwi_workflow.connect(jhu_atlas_to_t1w_node, 'output_image', jhu_atlas_to_dwi_node, 'in_file')
  2007. dwi_workflow.connect(invert_dwi_to_t1w_reg_node, 'out_file', jhu_atlas_to_dwi_node, 'in_matrix_file')
  2008. dwi_workflow.connect(preproc_dwi_node, 'b0', jhu_atlas_to_dwi_node, 'reference')
  2009. jhu_atlas_to_dwi_node.inputs.interp = 'nearestneighbour'
  2010. jhu_atlas_to_dwi_node.inputs.apply_xfm = True
  2011. jhu_atlas_to_dwi_node.inputs.out_file = os.path.join(anat_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_space-{space_entity}_seg-JHU_dseg.nii.gz")
  2012. scalar_maps_for_jhu_node = Node(CalculateScalarMaps(), name='scalar_maps_for_jhu')
  2013. scalar_maps_for_jhu_node.inputs.colnames = ["FA", "MD", "FW (MarkVCID2)", "AD", "RD", "ODI", "ICVF", "ISOVF", "GQI_GFA", "GQI_ISO", "GQI_QA", "CHIDIA", "MK", "AK", "RK"]
  2014. scalar_maps_for_jhu_node.inputs.output_csv = os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_label-JHU_desc-mean_dwimap.csv")
  2015. scalar_maps_for_jhu_node.inputs.statistic = 'mean'
  2016. dwi_workflow.connect(jhu_atlas_to_dwi_node, 'out_file', scalar_maps_for_jhu_node, 'mask_file')
  2017. dwi_workflow.connect(dwi_metrics_node, 'out', scalar_maps_for_jhu_node, 'data_files')
  2018. scalar_maps_for_jhu_node.inputs.ignore_background = True
  2019. # visual pathway analysis: TDI image
  2020. if vp_extract:
  2021. scalar_maps_for_vptdi_node = MapNode(CalculateTDIWeightedScalars(), name='scalar_maps_for_vptdi', iterfield=['weight_file', 'output_csv'])
  2022. scalar_maps_for_vptdi_node.inputs.colnames = ["FA", "MD", "FW (MarkVCID2)", "AD", "RD", "ODI", "ICVF", "ISOVF", "GQI_GFA", "GQI_ISO", "GQI_QA", "CHIDIA", "WMHprobmap"]
  2023. dwi_workflow.connect(vp_dwi_metrics_node, 'out', scalar_maps_for_vptdi_node, 'data_files')
  2024. dwi_workflow.connect(tt_to_tdi_node, 'out_tdi', scalar_maps_for_vptdi_node, 'weight_file')
  2025. scalar_maps_for_vptdi_node.inputs.output_csv = [
  2026. os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_label-OT_desc-TDIweighted_dwimap.csv"),
  2027. os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_label-OT_desc-TDIweighted_dwimap.csv"),
  2028. os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_label-OR_desc-TDIweighted_dwimap.csv"),
  2029. os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_label-OR_desc-TDIweighted_dwimap.csv"),
  2030. os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-L_label-ML_desc-TDIweighted_dwimap.csv"),
  2031. os.path.join(dwi_metrics_output_dir, f"sub-{self.subject.subject_id}_ses-{self.session.session_id}_hemi-R_label-ML_desc-TDIweighted_dwimap.csv")
  2032. ]
  2033. scalar_maps_for_vptdi_node.inputs.background_value = -1
  2034. return dwi_workflow
  2035. def extract_results(self):
  2036. def _sanitize_metric_name(name):
  2037. """
  2038. Convert metric names to safe column names.
  2039. Example:
  2040. 'FW (MarkVCID2)' -> 'markvcid_FW'
  2041. """
  2042. mapping = {
  2043. "FW (MarkVCID2)": "markvcid_FW"
  2044. }
  2045. if name in mapping:
  2046. return mapping[name]
  2047. safe_name = str(name).strip()
  2048. safe_name = safe_name.replace(" ", "_")
  2049. safe_name = safe_name.replace("(", "")
  2050. safe_name = safe_name.replace(")", "")
  2051. safe_name = safe_name.replace("-", "_")
  2052. safe_name = safe_name.replace("/", "_")
  2053. return safe_name
  2054. def _read_single_row_csv_as_dict(csv_path):
  2055. """
  2056. Read a single-row CSV and return a dict using sanitized column names.
  2057. Missing file or empty file -> {}
  2058. """
  2059. if csv_path is None or not os.path.exists(csv_path):
  2060. return {}
  2061. df = pd.read_csv(csv_path)
  2062. if df.empty or df.shape[0] < 1:
  2063. return {}
  2064. row0 = df.iloc[0]
  2065. out = {}
  2066. for col in df.columns:
  2067. safe_col = _sanitize_metric_name(col)
  2068. out[safe_col] = row0[col]
  2069. return out
  2070. def _read_alongtract_metrics_wide(csv_path, n_points=100):
  2071. """
  2072. Read along-tract metrics CSV (expected 100 rows, one row per tract point)
  2073. and flatten it into a single-row dict in wide format:
  2074. FA_001 ... FA_100
  2075. MD_001 ... MD_100
  2076. ...
  2077. Missing file -> empty dict
  2078. Missing rows -> pad with np.nan
  2079. Extra rows -> truncate to n_points
  2080. """
  2081. if csv_path is None or not os.path.exists(csv_path):
  2082. return {}
  2083. df = pd.read_csv(csv_path)
  2084. if df.empty:
  2085. return {}
  2086. df = df.reset_index(drop=True)
  2087. if df.shape[0] < n_points:
  2088. pad_df = pd.DataFrame(np.nan, index=range(n_points - df.shape[0]), columns=df.columns)
  2089. df = pd.concat([df, pad_df], ignore_index=True)
  2090. elif df.shape[0] > n_points:
  2091. df = df.iloc[:n_points, :].copy()
  2092. out = {}
  2093. for col in df.columns:
  2094. prefix = _sanitize_metric_name(col)
  2095. for i in range(n_points):
  2096. out[f"{prefix}_{i+1:03d}"] = df.iloc[i][col]
  2097. return out
  2098. def _read_wmparc_metrics_wide(csv_path, wmparc_dict):
  2099. if csv_path is None or not os.path.exists(csv_path):
  2100. return {}
  2101. df = pd.read_csv(csv_path)
  2102. if df.empty or "roi_label" not in df.columns:
  2103. return {}
  2104. metric_cols = [c for c in df.columns if c != "roi_label"]
  2105. out = {}
  2106. for metric in metric_cols:
  2107. metric_name = _sanitize_metric_name(metric)
  2108. for _, row in df.iterrows():
  2109. roi_label = int(row["roi_label"])
  2110. if roi_label not in wmparc_dict:
  2111. continue
  2112. roi_name = wmparc_dict[roi_label]
  2113. out[f"{metric_name}_{roi_name}"] = row[metric]
  2114. return out
  2115. os.makedirs(self.output_path, exist_ok=True)
  2116. dwi_output_path = self.extract_from
  2117. # DTI-ALPS
  2118. alps_columns = ["Subject", "Session", "ALPS_L", "ALPS_R", "ALPS_mean", "G-ALPS_L", "G-ALPS_R", "G-ALPS_avg", "G-ALPS_sum"]
  2119. alps_results_df = pd.DataFrame(columns=alps_columns)
  2120. # PSMD
  2121. psmd_columns = ["Subject", "Session", "PSMD", "PSMD_Left", "PSMD_Right"]
  2122. psmd_results_df = pd.DataFrame(columns=psmd_columns)
  2123. # PVeD
  2124. pved_columns = ["Subject", "Session", "PVeD", "PVeD_L", "PVeD_R", "QA_index"]
  2125. pved_results_df = pd.DataFrame(columns=pved_columns)
  2126. # Track-based DWI metrics
  2127. track_dwi_metrics_columns = ["Subject", "Session", "Track_DTI_FA", "Track_DTI_MD", "Track_SS_FW", "Track_DTI_FW"]
  2128. track_dwi_metrics_df = pd.DataFrame(columns=track_dwi_metrics_columns)
  2129. # Mask-based mean dwimap results, grouped by mask type
  2130. mean_dwimap_seedmask_df = pd.DataFrame()
  2131. mean_dwimap_WMH_df = pd.DataFrame()
  2132. mean_dwimap_NAWM_df = pd.DataFrame()
  2133. mean_dwimap_NAWM_perroi_df = pd.DataFrame()
  2134. median_dwimap_wmparc_df = pd.DataFrame()
  2135. mean_dwimap_aparcaseg_df = pd.DataFrame()
  2136. mean_dwimap_chpseg_df = pd.DataFrame()
  2137. mean_dwimap_jhu_df = pd.DataFrame()
  2138. # Surface parameters
  2139. surface_parameters_df = pd.DataFrame()
  2140. mirror_surface_parameters_df = pd.DataFrame()
  2141. # Along-tract DWI metrics
  2142. alongtract_L_OR_df = pd.DataFrame()
  2143. alongtract_L_OT_df = pd.DataFrame()
  2144. alongtract_R_OR_df = pd.DataFrame()
  2145. alongtract_R_OT_df = pd.DataFrame()
  2146. # TDIweighted dwimap results
  2147. tdiweighted_L_OR_df = pd.DataFrame()
  2148. tdiweighted_L_OT_df = pd.DataFrame()
  2149. tdiweighted_R_OR_df = pd.DataFrame()
  2150. tdiweighted_R_OT_df = pd.DataFrame()
  2151. # Parse FreeSurferColorLUT.txt for label lookups
  2152. import re
  2153. import xml.etree.ElementTree as ET
  2154. fs_lut = {}
  2155. lut_path = get_package_path('data', 'labelconvert_in', 'FreeSurferColorLUT.txt')
  2156. with open(lut_path, 'r') as f:
  2157. for raw_line in f:
  2158. line = raw_line.strip()
  2159. if not line or line.startswith('#'):
  2160. continue
  2161. m = re.match(r'^(\d+)\s+(.+?)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s*$', line)
  2162. if m:
  2163. fs_lut[int(m.group(1))] = m.group(2).strip()
  2164. # Build wmparc dict from FS LUT (WM labels in 3000-4999 range)
  2165. wmparc_dict = {k: v for k, v in fs_lut.items() if 3000 <= k <= 4999}
  2166. # Parse JHU-labels.xml
  2167. jhu_lut = {}
  2168. xml_path = get_package_path('data', 'standard', 'JHU', 'JHU-labels.xml')
  2169. tree = ET.parse(xml_path)
  2170. for label_elem in tree.getroot().findall('.//label'):
  2171. jhu_lut[int(label_elem.get('index'))] = label_elem.text.strip()
  2172. print(f"Reading results from {dwi_output_path}...")
  2173. for subject_folder in os.listdir(dwi_output_path):
  2174. if not subject_folder.startswith("sub-"):
  2175. continue
  2176. subject_id = subject_folder.split("-", 1)[1]
  2177. subject_folder_path = os.path.join(dwi_output_path, subject_folder)
  2178. if not os.path.isdir(subject_folder_path):
  2179. continue
  2180. for session_folder in os.listdir(subject_folder_path):
  2181. if not session_folder.startswith("ses-"):
  2182. continue
  2183. session_id = session_folder.split("-", 1)[1]
  2184. session_path = os.path.join(subject_folder_path, session_folder)
  2185. if not os.path.isdir(session_path):
  2186. continue
  2187. dwi_metrics_stats_dir = os.path.join(session_path, "dwi_metrics_stats")
  2188. dti_alps_csv = os.path.join(session_path, "DTI-ALPS", "alps.stat", "alps.csv")
  2189. psmd_txt = os.path.join(session_path, "psmd", "psmd_out.txt")
  2190. pved_csv = os.path.join(session_path, "PVeD", "PVeD_metrics.csv")
  2191. tract_fa_csv = os.path.join(dwi_metrics_stats_dir, "track_in_dti_FA_mean.csv")
  2192. tract_md_csv = os.path.join(dwi_metrics_stats_dir, "track_in_dti_MD_mean.csv")
  2193. tract_ss_fw_csv = os.path.join(dwi_metrics_stats_dir, "track_in_freewater_mean.csv")
  2194. tract_dti_fw_csv = os.path.join(dwi_metrics_stats_dir, "track_in_dti_freewater_mean.csv")
  2195. seedmask_mean_csv = os.path.join(
  2196. dwi_metrics_stats_dir,
  2197. f"sub-{subject_id}_ses-{session_id}_label-seedmask_desc-mean_dwimap.csv"
  2198. )
  2199. wmh_mean_csv = os.path.join(
  2200. dwi_metrics_stats_dir,
  2201. f"sub-{subject_id}_ses-{session_id}_label-WMH_desc-mean_dwimap.csv"
  2202. )
  2203. nawm_mean_csv = os.path.join(
  2204. dwi_metrics_stats_dir,
  2205. f"sub-{subject_id}_ses-{session_id}_label-NAWM_desc-mean_dwimap.csv"
  2206. )
  2207. nawm_perroi_csv = os.path.join(
  2208. dwi_metrics_stats_dir,
  2209. f"sub-{subject_id}_ses-{session_id}_label-NAWM_desc-perROI_dwimap.csv"
  2210. )
  2211. wmparc_median_csv = os.path.join(
  2212. dwi_metrics_stats_dir,
  2213. f"sub-{subject_id}_ses-{session_id}_label-wmparc_desc-median_dwimap.csv"
  2214. )
  2215. aparcaseg_mean_csv = os.path.join(
  2216. dwi_metrics_stats_dir,
  2217. f"sub-{subject_id}_ses-{session_id}_label-aparcaseg_desc-mean_dwimap.csv"
  2218. )
  2219. chpseg_mean_csv = os.path.join(
  2220. dwi_metrics_stats_dir,
  2221. f"sub-{subject_id}_ses-{session_id}_label-chpseg_desc-mean_dwimap.csv"
  2222. )
  2223. jhu_mean_csv = os.path.join(
  2224. dwi_metrics_stats_dir,
  2225. f"sub-{subject_id}_ses-{session_id}_label-JHU_desc-mean_dwimap.csv"
  2226. )
  2227. surface_csv = os.path.join(session_path, "surface_parameters.csv")
  2228. mirror_surface_csv = os.path.join(session_path, "mirror_surface_parameters.csv")
  2229. alongtract_L_OR_csv = os.path.join(
  2230. dwi_metrics_stats_dir,
  2231. f"sub-{subject_id}_ses-{session_id}_hemi-L_label-OR_desc-alongtract_dwimap.csv"
  2232. )
  2233. alongtract_L_OT_csv = os.path.join(
  2234. dwi_metrics_stats_dir,
  2235. f"sub-{subject_id}_ses-{session_id}_hemi-L_label-OT_desc-alongtract_dwimap.csv"
  2236. )
  2237. alongtract_R_OR_csv = os.path.join(
  2238. dwi_metrics_stats_dir,
  2239. f"sub-{subject_id}_ses-{session_id}_hemi-R_label-OR_desc-alongtract_dwimap.csv"
  2240. )
  2241. alongtract_R_OT_csv = os.path.join(
  2242. dwi_metrics_stats_dir,
  2243. f"sub-{subject_id}_ses-{session_id}_hemi-R_label-OT_desc-alongtract_dwimap.csv"
  2244. )
  2245. tdiweighted_L_OR_csv = os.path.join(
  2246. dwi_metrics_stats_dir,
  2247. f"sub-{subject_id}_ses-{session_id}_hemi-L_label-OR_desc-TDIweighted_dwimap.csv"
  2248. )
  2249. tdiweighted_L_OT_csv = os.path.join(
  2250. dwi_metrics_stats_dir,
  2251. f"sub-{subject_id}_ses-{session_id}_hemi-L_label-OT_desc-TDIweighted_dwimap.csv"
  2252. )
  2253. tdiweighted_R_OR_csv = os.path.join(
  2254. dwi_metrics_stats_dir,
  2255. f"sub-{subject_id}_ses-{session_id}_hemi-R_label-OR_desc-TDIweighted_dwimap.csv"
  2256. )
  2257. tdiweighted_R_OT_csv = os.path.join(
  2258. dwi_metrics_stats_dir,
  2259. f"sub-{subject_id}_ses-{session_id}_hemi-R_label-OT_desc-TDIweighted_dwimap.csv"
  2260. )
  2261. # DTI-ALPS
  2262. if os.path.exists(dti_alps_csv):
  2263. alps_df = pd.read_csv(dti_alps_csv)
  2264. if {"alps_L", "alps_R", "alps"}.issubset(alps_df.columns):
  2265. new_data = pd.DataFrame([{
  2266. "Subject": f"sub-{subject_id}",
  2267. "Session": f"ses-{session_id}",
  2268. "ALPS_L": alps_df["alps_L"].values[0],
  2269. "ALPS_R": alps_df["alps_R"].values[0],
  2270. "ALPS_mean": alps_df["alps"].values[0],
  2271. "G-ALPS_L": alps_df["g_alps_L"].values[0] if "g_alps_L" in alps_df.columns else np.nan,
  2272. "G-ALPS_R": alps_df["g_alps_R"].values[0] if "g_alps_R" in alps_df.columns else np.nan,
  2273. "G-ALPS_avg": alps_df["g_alps_avg"].values[0] if "g_alps_avg" in alps_df.columns else np.nan,
  2274. "G-ALPS_sum": alps_df["g_alps_sum"].values[0] if "g_alps_sum" in alps_df.columns else np.nan
  2275. }])
  2276. alps_results_df = pd.concat([alps_results_df, new_data], ignore_index=True)
  2277. # PSMD
  2278. if os.path.exists(psmd_txt):
  2279. with open(psmd_txt, "r") as f:
  2280. line = f.readline().strip()
  2281. parts = line.split()
  2282. if len(parts) == 1:
  2283. new_data = pd.DataFrame([{
  2284. "Subject": f"sub-{subject_id}",
  2285. "Session": f"ses-{session_id}",
  2286. "PSMD": float(parts[0]),
  2287. "PSMD_Left": np.nan,
  2288. "PSMD_Right": np.nan
  2289. }])
  2290. elif len(parts) == 2:
  2291. new_data = pd.DataFrame([{
  2292. "Subject": f"sub-{subject_id}",
  2293. "Session": f"ses-{session_id}",
  2294. "PSMD": np.nan,
  2295. "PSMD_Left": float(parts[0]),
  2296. "PSMD_Right": float(parts[1])
  2297. }])
  2298. else:
  2299. new_data = pd.DataFrame([{
  2300. "Subject": f"sub-{subject_id}",
  2301. "Session": f"ses-{session_id}",
  2302. "PSMD": np.nan,
  2303. "PSMD_Left": np.nan,
  2304. "PSMD_Right": np.nan
  2305. }])
  2306. psmd_results_df = pd.concat([psmd_results_df, new_data], ignore_index=True)
  2307. # PVeD
  2308. if os.path.exists(pved_csv):
  2309. pved_df = pd.read_csv(pved_csv)
  2310. if {"PVeD_total", "PVeD_L", "PVeD_R"}.issubset(pved_df.columns):
  2311. qa_value = pved_df["QA_index"].values[0] if "QA_index" in pved_df.columns else np.nan
  2312. new_data = pd.DataFrame([{
  2313. "Subject": f"sub-{subject_id}",
  2314. "Session": f"ses-{session_id}",
  2315. "PVeD": pved_df["PVeD_total"].values[0],
  2316. "PVeD_L": pved_df["PVeD_L"].values[0],
  2317. "PVeD_R": pved_df["PVeD_R"].values[0],
  2318. "QA_index": qa_value
  2319. }])
  2320. pved_results_df = pd.concat([pved_results_df, new_data], ignore_index=True)
  2321. # Track-based DWI metrics
  2322. track_csvs = [tract_fa_csv, tract_md_csv, tract_ss_fw_csv, tract_dti_fw_csv]
  2323. if any(os.path.exists(csv) for csv in track_csvs):
  2324. new_data = pd.DataFrame([{
  2325. "Subject": f"sub-{subject_id}",
  2326. "Session": f"ses-{session_id}",
  2327. "Track_DTI_FA": pd.read_csv(tract_fa_csv)["mean"].values[0] if os.path.exists(tract_fa_csv) else None,
  2328. "Track_DTI_MD": pd.read_csv(tract_md_csv)["mean"].values[0] if os.path.exists(tract_md_csv) else None,
  2329. "Track_SS_FW": pd.read_csv(tract_ss_fw_csv)["mean"].values[0] if os.path.exists(tract_ss_fw_csv) else None,
  2330. "Track_DTI_FW": pd.read_csv(tract_dti_fw_csv)["mean"].values[0] if os.path.exists(tract_dti_fw_csv) else None
  2331. }])
  2332. track_dwi_metrics_df = pd.concat([track_dwi_metrics_df, new_data], ignore_index=True)
  2333. # Mean dwimap results grouped by mask type
  2334. mean_specs = [
  2335. ("seedmask", seedmask_mean_csv, mean_dwimap_seedmask_df),
  2336. ("WMH", wmh_mean_csv, mean_dwimap_WMH_df),
  2337. ("NAWM", nawm_mean_csv, mean_dwimap_NAWM_df),
  2338. ]
  2339. wmparc_metric_dict = _read_wmparc_metrics_wide(wmparc_median_csv, wmparc_dict)
  2340. if wmparc_metric_dict:
  2341. row = {
  2342. "Subject": f"sub-{subject_id}",
  2343. "Session": f"ses-{session_id}",
  2344. }
  2345. row.update(wmparc_metric_dict)
  2346. median_dwimap_wmparc_df = pd.concat(
  2347. [median_dwimap_wmparc_df, pd.DataFrame([row])],
  2348. ignore_index=True,
  2349. )
  2350. # aparc+aseg metrics
  2351. aparcaseg_metric_dict = _read_wmparc_metrics_wide(aparcaseg_mean_csv, fs_lut)
  2352. if aparcaseg_metric_dict:
  2353. row = {"Subject": f"sub-{subject_id}", "Session": f"ses-{session_id}"}
  2354. row.update(aparcaseg_metric_dict)
  2355. mean_dwimap_aparcaseg_df = pd.concat([mean_dwimap_aparcaseg_df, pd.DataFrame([row])], ignore_index=True)
  2356. # chpseg metrics (1=left, 2=right)
  2357. chpseg_dict = {1: "chp_L", 2: "chp_R"}
  2358. chpseg_metric_dict = _read_wmparc_metrics_wide(chpseg_mean_csv, chpseg_dict)
  2359. if chpseg_metric_dict:
  2360. row = {"Subject": f"sub-{subject_id}", "Session": f"ses-{session_id}"}
  2361. row.update(chpseg_metric_dict)
  2362. mean_dwimap_chpseg_df = pd.concat([mean_dwimap_chpseg_df, pd.DataFrame([row])], ignore_index=True)
  2363. # JHU atlas metrics
  2364. jhu_metric_dict = _read_wmparc_metrics_wide(jhu_mean_csv, jhu_lut)
  2365. if jhu_metric_dict:
  2366. row = {"Subject": f"sub-{subject_id}", "Session": f"ses-{session_id}"}
  2367. row.update(jhu_metric_dict)
  2368. mean_dwimap_jhu_df = pd.concat([mean_dwimap_jhu_df, pd.DataFrame([row])], ignore_index=True)
  2369. updated_mean_dfs = {}
  2370. for mask_label, csv_path, target_df in mean_specs:
  2371. metric_dict = _read_single_row_csv_as_dict(csv_path)
  2372. if metric_dict:
  2373. row = {
  2374. "DWI_Pipeline_ID": f"sub-{subject_id}_ses-{session_id}",
  2375. "Subject": f"sub-{subject_id}",
  2376. "Session": f"ses-{session_id}",
  2377. "Mask": mask_label,
  2378. }
  2379. row.update(metric_dict)
  2380. target_df = pd.concat([target_df, pd.DataFrame([row])], ignore_index=True)
  2381. updated_mean_dfs[mask_label] = target_df
  2382. mean_dwimap_seedmask_df = updated_mean_dfs["seedmask"]
  2383. mean_dwimap_WMH_df = updated_mean_dfs["WMH"]
  2384. mean_dwimap_NAWM_df = updated_mean_dfs["NAWM"]
  2385. # NAWM per-ROI metrics (2 rows: roi_label=1 left, roi_label=2 right)
  2386. nawm_perroi_dict = {1: "NAWM_Left", 2: "NAWM_Right"}
  2387. nawm_perroi_metric_dict = _read_wmparc_metrics_wide(nawm_perroi_csv, nawm_perroi_dict)
  2388. if nawm_perroi_metric_dict:
  2389. row = {"Subject": f"sub-{subject_id}", "Session": f"ses-{session_id}"}
  2390. row.update(nawm_perroi_metric_dict)
  2391. mean_dwimap_NAWM_perroi_df = pd.concat([mean_dwimap_NAWM_perroi_df, pd.DataFrame([row])], ignore_index=True)
  2392. # Surface parameters
  2393. if os.path.exists(surface_csv):
  2394. df = pd.read_csv(surface_csv)
  2395. df.insert(0, "DWI_Pipeline_ID", f"sub-{subject_id}_ses-{session_id}")
  2396. df.insert(1, "Subject", f"sub-{subject_id}")
  2397. df.insert(2, "Session", f"ses-{session_id}")
  2398. surface_parameters_df = pd.concat([surface_parameters_df, df], ignore_index=True)
  2399. if os.path.exists(mirror_surface_csv):
  2400. df_mirror = pd.read_csv(mirror_surface_csv)
  2401. df_mirror.insert(0, "DWI_Pipeline_ID", f"sub-{subject_id}_ses-{session_id}")
  2402. df_mirror.insert(1, "Subject", f"sub-{subject_id}")
  2403. df_mirror.insert(2, "Session", f"ses-{session_id}")
  2404. mirror_surface_parameters_df = pd.concat([mirror_surface_parameters_df, df_mirror], ignore_index=True)
  2405. # Visual Pathway: Along-tract DWI metrics
  2406. alongtract_specs = [
  2407. ("L", "OR", alongtract_L_OR_csv, alongtract_L_OR_df),
  2408. ("L", "OT", alongtract_L_OT_csv, alongtract_L_OT_df),
  2409. ("R", "OR", alongtract_R_OR_csv, alongtract_R_OR_df),
  2410. ("R", "OT", alongtract_R_OT_csv, alongtract_R_OT_df),
  2411. ]
  2412. updated_alongtract_dfs = {}
  2413. for hemi, tract_label, csv_path, target_df in alongtract_specs:
  2414. wide_metrics = _read_alongtract_metrics_wide(csv_path, n_points=100)
  2415. if wide_metrics:
  2416. row = {
  2417. "DWI_Pipeline_ID": f"sub-{subject_id}_ses-{session_id}",
  2418. "Subject": f"sub-{subject_id}",
  2419. "Session": f"ses-{session_id}",
  2420. "Hemisphere": hemi,
  2421. "Tract": tract_label,
  2422. }
  2423. row.update(wide_metrics)
  2424. target_df = pd.concat([target_df, pd.DataFrame([row])], ignore_index=True)
  2425. updated_alongtract_dfs[(hemi, tract_label)] = target_df
  2426. alongtract_L_OR_df = updated_alongtract_dfs[("L", "OR")]
  2427. alongtract_L_OT_df = updated_alongtract_dfs[("L", "OT")]
  2428. alongtract_R_OR_df = updated_alongtract_dfs[("R", "OR")]
  2429. alongtract_R_OT_df = updated_alongtract_dfs[("R", "OT")]
  2430. # Visual Pathway: TDIweighted dwimap results
  2431. tdiweighted_specs = [
  2432. ("L", "OR", tdiweighted_L_OR_csv, tdiweighted_L_OR_df),
  2433. ("L", "OT", tdiweighted_L_OT_csv, tdiweighted_L_OT_df),
  2434. ("R", "OR", tdiweighted_R_OR_csv, tdiweighted_R_OR_df),
  2435. ("R", "OT", tdiweighted_R_OT_csv, tdiweighted_R_OT_df),
  2436. ]
  2437. updated_tdiweighted_dfs = {}
  2438. for hemi, tract_label, csv_path, target_df in tdiweighted_specs:
  2439. metric_dict = _read_single_row_csv_as_dict(csv_path)
  2440. if metric_dict:
  2441. row = {
  2442. "DWI_Pipeline_ID": f"sub-{subject_id}_ses-{session_id}",
  2443. "Subject": f"sub-{subject_id}",
  2444. "Session": f"ses-{session_id}",
  2445. "Hemisphere": hemi,
  2446. "Tract": tract_label,
  2447. }
  2448. row.update(metric_dict)
  2449. target_df = pd.concat([target_df, pd.DataFrame([row])], ignore_index=True)
  2450. updated_tdiweighted_dfs[(hemi, tract_label)] = target_df
  2451. tdiweighted_L_OR_df = updated_tdiweighted_dfs[("L", "OR")]
  2452. tdiweighted_L_OT_df = updated_tdiweighted_dfs[("L", "OT")]
  2453. tdiweighted_R_OR_df = updated_tdiweighted_dfs[("R", "OR")]
  2454. tdiweighted_R_OT_df = updated_tdiweighted_dfs[("R", "OT")]
  2455. # Save results
  2456. alps_output_excel = os.path.join(self.output_path, "dwi_alps_summary.xlsx")
  2457. psmd_output_excel = os.path.join(self.output_path, "dwi_psmd_summary.xlsx")
  2458. pved_output_excel = os.path.join(self.output_path, "dwi_pved_summary.xlsx")
  2459. track_dwi_metrics_output_excel = os.path.join(self.output_path, "dwi_track_metrics_summary.xlsx")
  2460. mean_dwimap_seedmask_output_xlsx = os.path.join(self.output_path, "dwi_seedmask_summary.xlsx")
  2461. mean_dwimap_WMH_output_xlsx = os.path.join(self.output_path, "dwi_WMH_summary.xlsx")
  2462. mean_dwimap_NAWM_output_xlsx = os.path.join(self.output_path, "dwi_NAWM_summary.xlsx")
  2463. mean_dwimap_NAWM_perroi_output_xlsx = os.path.join(self.output_path, "dwi_NAWM_perROI_summary.xlsx")
  2464. median_dwimap_wmparc_output_xlsx = os.path.join(self.output_path, "dwi_wmparc_summary.xlsx")
  2465. surface_output_excel = os.path.join(self.output_path, "dwi_surface_parameters_summary.xlsx")
  2466. mirror_output_excel = os.path.join(self.output_path, "dwi_mirror_surface_parameters_summary.xlsx")
  2467. alongtract_L_OR_output_xlsx = os.path.join(self.output_path, "dwi_alongtract_hemi-L_label-OR_summary.xlsx")
  2468. alongtract_L_OT_output_xlsx = os.path.join(self.output_path, "dwi_alongtract_hemi-L_label-OT_summary.xlsx")
  2469. alongtract_R_OR_output_xlsx = os.path.join(self.output_path, "dwi_alongtract_hemi-R_label-OR_summary.xlsx")
  2470. alongtract_R_OT_output_xlsx = os.path.join(self.output_path, "dwi_alongtract_hemi-R_label-OT_summary.xlsx")
  2471. tdiweighted_L_OR_output_xlsx = os.path.join(self.output_path, "dwi_TDIweighted_hemi-L_label-OR_summary.xlsx")
  2472. tdiweighted_L_OT_output_xlsx = os.path.join(self.output_path, "dwi_TDIweighted_hemi-L_label-OT_summary.xlsx")
  2473. tdiweighted_R_OR_output_xlsx = os.path.join(self.output_path, "dwi_TDIweighted_hemi-R_label-OR_summary.xlsx")
  2474. tdiweighted_R_OT_output_xlsx = os.path.join(self.output_path, "dwi_TDIweighted_hemi-R_label-OT_summary.xlsx")
  2475. if not alps_results_df.empty:
  2476. alps_results_df.to_excel(alps_output_excel, header=True, index=False)
  2477. print(f"DTI-ALPS results saved to {alps_output_excel}")
  2478. else:
  2479. print("No DTI-ALPS results found.")
  2480. if not psmd_results_df.empty:
  2481. psmd_results_df.to_excel(psmd_output_excel, header=True, index=False)
  2482. print(f"PSMD results saved to {psmd_output_excel}")
  2483. else:
  2484. print("No PSMD results found.")
  2485. if not pved_results_df.empty:
  2486. pved_results_df.to_excel(pved_output_excel, header=True, index=False)
  2487. print(f"PVeD results saved to {pved_output_excel}")
  2488. else:
  2489. print("No PVeD results found.")
  2490. if not track_dwi_metrics_df.empty:
  2491. track_dwi_metrics_df.to_excel(track_dwi_metrics_output_excel, header=True, index=False)
  2492. print(f"Track-based DWI metrics results saved to {track_dwi_metrics_output_excel}")
  2493. else:
  2494. print("No track-based DWI metrics results found.")
  2495. if not mean_dwimap_seedmask_df.empty:
  2496. mean_dwimap_seedmask_df.to_excel(mean_dwimap_seedmask_output_xlsx, index=False)
  2497. print(f"Mean dwimap seedmask results saved to {mean_dwimap_seedmask_output_xlsx}")
  2498. else:
  2499. print("No mean dwimap seedmask results found.")
  2500. if not mean_dwimap_WMH_df.empty:
  2501. mean_dwimap_WMH_df.to_excel(mean_dwimap_WMH_output_xlsx, index=False)
  2502. print(f"Mean dwimap WMH results saved to {mean_dwimap_WMH_output_xlsx}")
  2503. else:
  2504. print("No mean dwimap WMH results found.")
  2505. if not mean_dwimap_NAWM_df.empty:
  2506. mean_dwimap_NAWM_df.to_excel(mean_dwimap_NAWM_output_xlsx, index=False)
  2507. print(f"Mean dwimap NAWM results saved to {mean_dwimap_NAWM_output_xlsx}")
  2508. else:
  2509. print("No mean dwimap NAWM results found.")
  2510. if not mean_dwimap_NAWM_perroi_df.empty:
  2511. mean_dwimap_NAWM_perroi_df.to_excel(mean_dwimap_NAWM_perroi_output_xlsx, index=False)
  2512. print(f"Mean dwimap NAWM per-ROI results saved to {mean_dwimap_NAWM_perroi_output_xlsx}")
  2513. else:
  2514. print("No mean dwimap NAWM per-ROI results found.")
  2515. if not median_dwimap_wmparc_df.empty:
  2516. median_dwimap_wmparc_df.to_excel(median_dwimap_wmparc_output_xlsx, index=False)
  2517. print(f"Median dwimap wmparc results saved to {median_dwimap_wmparc_output_xlsx}")
  2518. else:
  2519. print("No median dwimap wmparc results found.")
  2520. if not mean_dwimap_aparcaseg_df.empty:
  2521. mean_dwimap_aparcaseg_df.to_excel(os.path.join(self.output_path, "dwi_aparcaseg_summary.xlsx"), index=False)
  2522. print(f"Mean dwimap aparcaseg results saved to {os.path.join(self.output_path, 'dwi_aparcaseg_summary.xlsx')}")
  2523. else:
  2524. print("No mean dwimap aparcaseg resul

dwi_pipeline.py at commit 93e022c, under MIT · at the source

Overview

Authors: Youjie Wang1, Jingyu Cui1, Yuying Yan1, Tang Yang1, Yue Yuan1, Rumei Lei1, Rongfeng Luo1, Bo Wu1, Shuai Jiang1
ORCID iDs: Youjie Wang
  1. Department of Neurology, West China Hospital of Sichuan University, No. 37 Guo Xue Xiang, Chengdu 610041, China
Journal: NeuroImage. Clinical, volume 50, article 103987
Dates: received 20 November 2025; accepted 16 March 2026; published online 20 March 2026; in print 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1016/j.nicl.2026.103987 · PMID 41875823 · PMCID PMC13049619 · OpenAlex W7139949947
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: structural MRI / diffusion (modality), human (organism), stroke (population), clinical / translational (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Connectivity, fMRI & imaging, Preprocessing
Keywords: Secondary degeneration, Recent small subcortical infarcts, Diffusion tensor imaging
MeSH: Cerebral Cortex*, Cerebral Infarction*, Aged, Connectome, Diffusion Magnetic Resonance Imaging, Female, Humans, Longitudinal Studies, Magnetic Resonance Imaging, Male, Middle Aged (* major topic)
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: West China Hospital, Sichuan University; Sichuan University (2024HXBH041, 2023HXBH007, U24A20690); National Natural Science Foundation of China (82271328, 82371322, 82301661); Sichuan Provincial Science and Technology Support Program (2024YFFK0314)
Citations: not cited yet (Europe PMC); 43 references in the paper

Abstract

Background: Secondary cortical degeneration caused by the remote effects of subcortical infarcts has been implicated in long-term outcomes after acute ischemic stroke. However, this process remains insufficiently studied in recent small subcortical infarcts (RSSI). We aimed to verify RSSI-induced cortical damage, determine whether it can be captured by neuroimaging markers, and explore its association with clinical outcomes.

Methods: RSSI patients with longitudinal Magnetic Resonance Imaging (MRI) were included. Cortical degeneration was assessed using linear mixed-effects models, incorporating a direct approach based on individual diffusion weighted imaging and an indirect approach using the normative connectome from the Human Connectome Project (HCP). Principal component analysis (PCA) was employed to extract features of cortical alterations. The resulting component scores were used in general linear models to assess associations with neuroimaging markers and clinical outcomes.

Results: A total of 76 RSSI patients were analyzed. RSSI was found to induce progressive cortical thinning and volume loss in structurally connected regions. PCA identified a component reflecting parenchymal atrophy associated with diffusion-based markers of white matter integrity, as well as the presence of track/cap signs. Moreover, faster cortical degeneration in lesion-connected regions was significantly associated with a greater increase in Hamilton Anxiety Rating Scale (HAMA) scores (β = -2.38, 95% CI = -4.30 – -0.47, p = 0.017).

Conclusions: RSSI induces secondary cortical damage through structurally connected fiber tracts, which is detectable by neuroimaging markers of white matter integrity. These regional cortical alterations may be relevant to post-stroke outcomes and require validation in larger longitudinal studies.

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

Repositories

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

kjamison/nemo

License: MIT
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: d661f3fddba70dff1eb93f0c57b4f726f09723f9, 23 July 2026
Languages: Python (19), Shell (5), JavaScript (1)
Size: 264 files, 25 scripts
Software Heritage: archived
Found in: the text, “Indirect estimation”
Holds: README, license file, environment (requirements.txt)
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (18 files), SciPy (17 files), NiBabel (15 files), Matplotlib (4 files), Nilearn (4 files), FSL (2 files), MRtrix3 (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
27 files

LuuuXG/cvdproc

License: MIT
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 93e022c07bcf075aaf1d727dd6ca2261b1a69c57, 11 September 2026
Languages: Python (308), Shell (81), JavaScript (36), R (27), MATLAB (24), Jupyter (4)
Size: 771 files, 480 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, license file, environment (pyproject.toml, requirements.txt, setup.py), tests, documentation
Not found: CITATION.cff, continuous integration
Tools: NumPy (155 files), NiBabel (126 files), Nipype (119 files), pandas (55 files), FSL (54 files), FreeSurfer (46 files), SciPy (41 files), Matplotlib (33 files), PyTorch (14 files), DIPY (11 files), scikit-image (11 files), MRtrix3 (10 files), ANTs (9 files), Image Processing Toolbox (7 files), MONAI (6 files), Pillow (6 files), neuromaps (5 files), Tools for NIfTI and ANALYZE image (MATLAB) (5 files), ggplot2 (3 files), igraph (3 files), pydicom (3 files), SimpleITK (3 files), tidyverse (3 files), ggseg (2 files), imageio (2 files), Nilearn (2 files), reshape2 (2 files), TensorFlow (2 files), AFNI (1 file), Brain Connectivity Toolbox (1 file), Dcm2Bids (1 file), fMRIPrep (1 file), h5py (1 file), patchwork (1 file), PyBIDS (1 file), QSIPrep (1 file), scikit-learn (1 file), SPM (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
482 files

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

Tracing map

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

What the map holds:

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

Data will be made available on request.

Mr. Wang had full access to all the data in the study and took full responsibility for the integrity of the data and the accuracy of the analyses. The neuroimaging analysis code is available at https://github.com/LuuuXG/cvdproc. The source data is available from the corresponding author upon reasonable request and subject to institutional and ethical regulations.

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

Versions

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

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 9 authors, 3 keywords, 11 MeSH terms, 4 funders, 43 references.

Cite

This paper

Wang, Y., Cui, J., Yan, Y., Yang, T., Yuan, Y., Lei, R., Luo, R., Wu, B., & Jiang, S. (2026). Remote cortical degeneration related to structural connectivity following recent small subcortical infarcts. NeuroImage. Clinical, 50, 103987. https://doi.org/10.1016/j.nicl.2026.103987

BibTeX

@article{wang2026remote,
author = {Wang, Youjie and Cui, Jingyu and Yan, Yuying and Yang, Tang and Yuan, Yue and Lei, Rumei and Luo, Rongfeng and Wu, Bo and Jiang, Shuai},
title = {{Remote cortical degeneration related to structural connectivity following recent small subcortical infarcts}},
journal = {NeuroImage. Clinical},
year = {2026},
month = mar,
volume = {50},
pages = {103987},
publisher = {Elsevier},
issn = {2213-1582},
doi = {10.1016/j.nicl.2026.103987},
url = {https://doi.org/10.1016/j.nicl.2026.103987},
pmid = {41875823},
pmcid = {PMC13049619}
}

RIS

TY - JOUR
AU - Wang, Youjie
AU - Cui, Jingyu
AU - Yan, Yuying
AU - Yang, Tang
AU - Yuan, Yue
AU - Lei, Rumei
AU - Luo, Rongfeng
AU - Wu, Bo
AU - Jiang, Shuai
TI - Remote cortical degeneration related to structural connectivity following recent small subcortical infarcts
T2 - NeuroImage. Clinical
J2 - Neuroimage Clin
PY - 2026
DA - 2026/03/20
VL - 50
SP - 103987
SN - 2213-1582
PB - Elsevier
DO - 10.1016/j.nicl.2026.103987
UR - https://doi.org/10.1016/j.nicl.2026.103987
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.nicl.2026.103987",
"type": "article-journal",
"title": "Remote cortical degeneration related to structural connectivity following recent small subcortical infarcts",
"container-title": "NeuroImage. Clinical",
"author": [
{
"family": "Wang",
"given": "Youjie"
},
{
"family": "Cui",
"given": "Jingyu"
},
{
"family": "Yan",
"given": "Yuying"
},
{
"family": "Yang",
"given": "Tang"
},
{
"family": "Yuan",
"given": "Yue"
},
{
"family": "Lei",
"given": "Rumei"
},
{
"family": "Luo",
"given": "Rongfeng"
},
{
"family": "Wu",
"given": "Bo"
},
{
"family": "Jiang",
"given": "Shuai"
}
],
"container-title-short": "Neuroimage Clin",
"volume": "50",
"page": "103987",
"DOI": "10.1016/j.nicl.2026.103987",
"PMID": "41875823",
"PMCID": "PMC13049619",
"ISSN": "2213-1582",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.nicl.2026.103987",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
20
]
]
}
}

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.3390/bs15060746 [code]
Brain White Matter Alterations in Young Adults with Childhood Emotional Neglect Experience.
Journal: Behavioral sciences (Basel, Switzerland)
In common: QSIPrep, PyBIDS, MRtrix3, 20 other tools, structural MRI / diffusion
[2] doi:10.1038/s41592-025-02599-1 [code]
DeepPrep: an accelerated, scalable and robust pipeline for neuroimaging preprocessing empowered by deep learning.
Journal: Nature methods
In common: fMRIPrep, neuromaps, PyBIDS, 20 other tools, structural MRI / diffusion, 1 reference
[3] doi:10.1038/s41467-026-68714-8 [code]
Two axes of white matter development.
Journal: Nature communications
In common: QSIPrep, fMRIPrep, ggseg, 18 other tools, structural MRI / diffusion, 1 reference
[4] doi:10.1002/alz.70196 [code]
Cognitive implications and associated transcriptomic signatures of distinct regional iron depositions in cerebral small vessel disease.
Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association
In common: ggseg, neuromaps, Nipype, 14 other tools, stroke, structural MRI / diffusion, 3 references, author Youjie Wang
[5] doi:10.1038/s41598-023-32924-7 [code]
Reproducible Tract Profiles 2 (RTP2) suite, from diffusion MRI acquisition to clinical practice and research.
Journal: Scientific reports
In common: fMRIPrep, PyBIDS, MRtrix3, 16 other tools, structural MRI / diffusion, clinical / translational, 3 references
[6] doi:10.1038/s41398-025-03614-x [code]
Sex differences in response to violence: role of salience network expansion and connectivity on depression.
Journal: Translational psychiatry
In common: fMRIPrep, neuromaps, PyBIDS, 18 other tools, 1 reference
[7] doi:10.1038/s42003-025-08573-z [code]
Profiling brain morphology for autism spectrum disorder with two cross-culture large-scale consortia.
Journal: Communications biology
In common: ggseg, MONAI, SimpleITK, 20 other tools, structural MRI / diffusion
[8] doi:10.1016/j.neuroimage.2023.120109 [code]
Benchmarking explanation methods for mental state decoding with deep learning models.
Journal: NeuroImage
In common: fMRIPrep, PyBIDS, DIPY, 18 other tools, 1 reference
[9] doi:10.1093/brain/awae315 [code]
Enhancing cognitive performance prediction by white matter hyperintensity connectivity assessment.
Journal: Brain : a journal of neurology
In common: neuromaps, Nipype, pydicom, 15 other tools, stroke, structural MRI / diffusion, 4 references
[10] doi:10.1038/s41586-025-09708-2 [code]
A probabilistic histological atlas of the human brain for MRI segmentation.
Journal: Nature
In common: MRtrix3, DIPY, Nipype, 18 other tools, structural MRI / diffusion, 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.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1016/j.nicl.2026.103987.svg)](https://openscicode.org/badge/doi/10.1016/j.nicl.2026.103987)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1016/j.nicl.2026.103987.svg
   :target: https://openscicode.org/badge/doi/10.1016/j.nicl.2026.103987
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1016/j.nicl.2026.103987"><img src="https://openscicode.org/badge/doi/10.1016/j.nicl.2026.103987.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge 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: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

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