diff --git a/examples/my_nodes/dp_node_ants_apply.m b/examples/my_nodes/dp_node_ants_apply.m new file mode 100644 index 0000000..42cf348 --- /dev/null +++ b/examples/my_nodes/dp_node_ants_apply.m @@ -0,0 +1,45 @@ +classdef dp_node_ants_apply < dp_node + % Wrapper around antsApplyTransforms onto a given reference grid + % input.nii_fn : image to be transformed (e.g. atlas labels in MNI) + % input.ref_fn : reference grid / output space (e.g. trimmed FA) + % input.affine_fn, input.invwarp_fn : from dp_node_ants_reg_quick + % output.nii_fn : resampled image in ref_fn space + + properties + options = ''; % e.g. ' -n NearestNeighbor' for labels + use_inverse = 1; % default = inverse chain + end + + methods + + function obj = dp_node_ants_apply(options) + if (nargin>0), obj.options = options; end + obj.output_test = {'nii_fn'}; + end + + function output = i2o(obj, input) + output.nii_fn = dp.new_fn(input.op, input.nii_fn, '_ants'); + output.ref_fn = input.ref_fn; + end + + function output = execute(obj, input, output) + msf_mkdir(fileparts(output.nii_fn)); + if obj.use_inverse + % apply inverse chain: template→subject (fixed → moving) + t = sprintf('-t "[%s,1]" -t "%s" ', input.affine_fn, input.invwarp_fn); + else + % apply forward chain: subject→template (moving → fixed) + t = sprintf('-t "%s" -t "%s" ', input.warp_fn, input.affine_fn); + end + + if isempty(obj.options), obj.options = ' -n NearestNeighbor'; end + + % cmd using composed transform chain "t" (forward or inverse, depending on use_inverse) + cmd = sprintf('antsApplyTransforms -d 3 -i "%s" -r "%s" %s -o "%s"%s', ... + input.nii_fn, input.ref_fn, t, output.nii_fn, obj.options); + + + system(cmd); + end + end +end diff --git a/examples/my_nodes/dp_node_ants_reg_quick.m b/examples/my_nodes/dp_node_ants_reg_quick.m new file mode 100644 index 0000000..c30f7be --- /dev/null +++ b/examples/my_nodes/dp_node_ants_reg_quick.m @@ -0,0 +1,64 @@ +classdef dp_node_ants_reg_quick < dp_node + % antsRegistrationSyN.sh (moving = input.nii_fn, fixed = input.target_fn) + % Outputs affine + forward/inverse warps in + + properties + options = ''; % include threads and preset here, e.g. ' ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=16 -t s' + end + + methods + + function obj = dp_node_ants_reg_quick(options) + if (nargin>0), obj.options = options; end + obj.output_test = {'affine_fn','warp_fn','invwarp_fn','warped_fn'}; + + end + + function output = i2o(obj, input) + tmp = msf_fn_append(input.nii_fn, '_mni_'); + prefix = msf_fn_new_path(input.op, tmp); + [pdir,pfx,~] = msf_fileparts(prefix); + prefix = fullfile(pdir, pfx); % ANTs prefix base (no ext) + + output.prefix = prefix; + output.target_fn = input.target_fn; + output.affine_fn = [prefix '0GenericAffine.mat']; + output.warp_fn = [prefix '1Warp.nii.gz']; + output.invwarp_fn = [prefix '1InverseWarp.nii.gz']; + output.warped_fn = [prefix 'Warped.nii.gz']; % QC, remains in op + end + + + function output = execute(obj, input, output) + % set threads and build the command + thr = feature('numcores'); + envp = sprintf('ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=%d', thr); + + % if user didn’t provide flags, default to symmetric registration (-t s) + if isempty(strtrim(obj.options)) + flags = ' -t s'; + else + flags = [' ' strtrim(obj.options)]; + end + + msf_mkdir(fileparts(output.warped_fn)); + + + % full SyN variant + % cmd = sprintf('%s antsRegistrationSyN.sh -d 3 -f "%s" -m "%s" -o "%s"%s', ... + % envp, input.target_fn, input.nii_fn, output.prefix, flags); + + % quick SyN variant + cmd = sprintf('%s antsRegistrationSyNQuick.sh -d 3 -f "%s" -m "%s" -o "%s"%s', ... + envp, input.target_fn, input.nii_fn, output.prefix, flags); + + + % Always save ANTs output to log in the same folder + log_fn = fullfile(fileparts(output.warped_fn), 'ants_registration.log'); + cmd = [cmd, sprintf(' > "%s" 2>&1', log_fn)]; + + system(cmd); + end + + end +end diff --git a/examples/my_nodes/dp_node_atlas_crop_to_mask.m b/examples/my_nodes/dp_node_atlas_crop_to_mask.m new file mode 100644 index 0000000..9e6bc2b --- /dev/null +++ b/examples/my_nodes/dp_node_atlas_crop_to_mask.m @@ -0,0 +1,53 @@ +classdef dp_node_atlas_crop_to_mask < dp_node + % Crop an atlas/label map with a provided mask (same grid). + % + % Inputs (from upstream): + % nii_fn - atlas / label image to crop (e.g. warped JHU labels) + % mask_fn - mask in the same space (e.g. trimmed FA brain mask) + % op - optional output directory + % + % Output: + % nii_fn - cropped atlas, saved as /_crop.nii.gz + + properties + suffix = '_crop'; + end + + methods + function obj = dp_node_atlas_crop_to_mask(suffix) + if nargin > 0 + obj.suffix = suffix; + end + obj.output_test = {'nii_fn'}; + end + + function output = i2o(obj, input) + % Decide output folder + if isfield(input, 'op') && ~isempty(input.op) + op = input.op; + else + [op, ~, ~] = msf_fileparts(input.nii_fn); + end + + % Use dp.new_fn to append suffix before extension + output.nii_fn = dp.new_fn(op, input.nii_fn, obj.suffix); + end + + function output = execute(obj, input, output) + % Read atlas/labels and mask + [A, Ah] = mdm_nii_read(input.nii_fn); + B = mdm_nii_read(input.mask_fn); + + % Basic safety + if ~isequal(size(A), size(B)) + error('dp_node_atlas_crop_to_mask: atlas/mask size mismatch'); + end + + msf_mkdir(fileparts(output.nii_fn)); + + % Apply crop + mdm_nii_write(A .* (B > 0), output.nii_fn, Ah); + + end + end +end diff --git a/examples/my_nodes/dp_node_segm_jhu.m b/examples/my_nodes/dp_node_segm_jhu.m new file mode 100644 index 0000000..80cdfab --- /dev/null +++ b/examples/my_nodes/dp_node_segm_jhu.m @@ -0,0 +1,71 @@ +classdef dp_node_segm_jhu < dp_node_segm + + methods (Hidden) + function [labels, ids] = segm_info(obj) + + % /usr/local/fsl/data/atlases/JHU-labels.xml + + txt = { ... + '0: Unclassified', ... % Unclassified + '1: MCP', ... % Middle cerebellar peduncle + '2: PCT', ... % Pontine crossing tract (a part of MCP) + '3: GCC', ... % Genu of corpus callosum + '4: CCB', ... % Body of corpus callosum + '5: SCC', ... % Splenium of corpus callosum + '6: FX', ... % Fornix (column and body of fornix) + '7: CST_R', ... % Corticospinal tract R + '8: CST_L', ... % Corticospinal tract L + '9: ML_R', ... % Medial lemniscus R + '10: ML_L', ... % Medial lemniscus L + '11: ICP_R', ... % Inferior cerebellar peduncle R + '12: ICP_L', ... % Inferior cerebellar peduncle L + '13: SCP_R', ... % Superior cerebellar peduncle R + '14: SCP_L', ... % Superior cerebellar peduncle L + '15: CP_R', ... % Cerebral peduncle R + '16: CP_L', ... % Cerebral peduncle L + '17: ALIC_R', ... % Anterior limb of internal capsule R + '18: ALIC_L', ... % Anterior limb of internal capsule L + '19: PLIC_R', ... % Posterior limb of internal capsule R + '20: PLIC_L', ... % Posterior limb of internal capsule L + '21: RLIC_R', ... % Retrolenticular part of internal capsule R + '22: RLIC_L', ... % Retrolenticular part of internal capsule L + '23: ACR_R', ... % Anterior corona radiata R + '24: ACR_L', ... % Anterior corona radiata L + '25: SCR_R', ... % Superior corona radiata R + '26: SCR_L', ... % Superior corona radiata L + '27: PCR_R', ... % Posterior corona radiata R + '28: PCR_L', ... % Posterior corona radiata L + '29: PTR_R', ... % Posterior thalamic radiation (include optic radiation) R + '30: PTR_L', ... % Posterior thalamic radiation (include optic radiation) L + '31: SS_R', ... % Sagittal stratum (incl ILF and IFOF) R + '32: SS_L', ... % Sagittal stratum (incl ILF and IFOF) L + '33: EC_R', ... % External capsule R + '34: EC_L', ... % External capsule L + '35: CGC_R', ... % Cingulum (cingulate gyrus) R + '36: CGC_L', ... % Cingulum (cingulate gyrus) L + '37: CGH_R', ... % Cingulum (hippocampus) R + '38: CGH_L', ... % Cingulum (hippocampus) L + '39: FST_R', ... % Fornix (crus) / Stria terminalis R + '40: FST_L', ... % Fornix (crus) / Stria terminalis L + '41: SLF_R', ... % Superior longitudinal fasciculus R + '42: SLF_L', ... % Superior longitudinal fasciculus L + '43: SFOF_R', ... % Superior fronto-occipital fasciculus R + '44: SFOF_L', ... % Superior fronto-occipital fasciculus L + '45: UF_R', ... % Uncinate fasciculus R + '46: UF_L', ... % Uncinate fasciculus L + '47: TAP_R', ... % Tapetum R + '48: TAP_L' ... % Tapetum L + }; + + + + + a = @(x) strsplit(x, ' '); + b = @(x,i) x{i}; + c = @(x) str2num(x(1:(end-1))); + + labels = cellfun(@(x) b(a(x),2), txt, 'UniformOutput', false); + ids = cellfun(@(x) c(b(a(x),1)), txt, 'UniformOutput', false); + end + end +end diff --git a/examples/my_nodes/dp_node_segm_xtractseg.m b/examples/my_nodes/dp_node_segm_xtractseg.m new file mode 100644 index 0000000..5e071eb --- /dev/null +++ b/examples/my_nodes/dp_node_segm_xtractseg.m @@ -0,0 +1,165 @@ +classdef dp_node_segm_xtractseg < dp_node_segm + + % inputs + % dmri_fn + % bval_fn + % bvec_fn + % mask_fn + + % to do: flow management, e.g. deleting files and more + + methods + + function obj = dp_node_segm_xtractseg() + obj.output_test = {'labels_fn'}; + end + + function output = i2o(obj, input) + + output.op = input.op; + + f = @(x) fullfile(output.op, x); + % output.mask_fn = f('nodif_brain_mask.nii.gz'); + % output.dmri_fn = f('dmri.nii.gz'); + % output.bval_fn = f('x.bvals'); + % output.bvec_fn = f('x.bvecs'); + + output.mask_fn = f('nodif_brain_mask_MNI.nii.gz'); + output.dmri_fn = f('Diffusion_MNI.nii.gz'); + output.fa_fn = f('FA_MNI.nii.gz'); + output.bval_fn = f('x.bvals'); + output.bvec_fn = f('x.bvecs'); + + + + output.labels_fn = f('bundles.nii.gz'); + + output.tmp.bp = msf_tmp_path(); + output.tmp.do_delete = 1; + + end + + function output = execute(obj, input, output) + + % prepare + msf_mkdir(output.op); + + % put brain mask and dmri into the right place + copyfile(input.dmri_fn, fullfile(output.op, 'dmri.nii.gz')); + copyfile(input.bval_fn, output.bval_fn); + copyfile(input.bvec_fn, output.bvec_fn); + + % Preliminary implementation + if (~isempty(input.mask_fn)) + copyfile(input.mask_fn, output.mask_fn); + end + + + % convert relative to absolute path + [~,tmp] = fileattrib(output.op); + op = tmp.Name; + + % cmd = sprintf(... + % 'docker run -v "%s":/data -t "%s" TractSeg -i "%s" -o /data %s --bvals "%s" --bvecs "%s" %s', ... + % op, ... + % 'wasserth/tractseg_container:master', ... % docker name + % 'data/dmri.nii.gz', ... + % '--raw_diffusion_input', ... % options + % 'data/x.bvals', ... + % 'data/x.bvecs', ... + % '--preprocess'); % register to mni space + + cmd = sprintf(... + 'docker run -v "%s":/data -t "%s" TractSeg -i "%s" -o /data %s --bvals "%s" --bvecs "%s" %s', ... + op, ... + 'wasserth/tractseg_container:master', ... + 'data/dmri.nii.gz', ... + '--raw_diffusion_input', ... + 'data/x.bvals', ... + 'data/x.bvecs', ... + '--preprocess --tract_definition xtract'); + + + obj.syscmd(cmd); + + + % wrap up the files in a 4D volume + h = mdm_nii_read_header(output.dmri_fn); + labels = obj.segm_labels(); + R = zeros(h.dim(2), h.dim(3), h.dim(4), numel(labels)); + for c = 1:numel(labels) + tmp = mdm_nii_read(fullfile(output.op, ... + 'bundle_segmentations_MNI', ... + cat(2, labels{c}, '.nii.gz'))); + R(:,:,:,c) = tmp; + end + mdm_nii_write(uint8(R), output.labels_fn, h); + + end + + end + + methods (Hidden) + + function [labels, ids] = segm_info(obj) + + % dp_node_segm_xtractseg.segm_info() — XTRACT bundle names + + txt = { ... + '0: ac', ... % Anterior commissure + '1: af_l', ... % Arcuate fasciculus (left) + '2: af_r', ... % Arcuate fasciculus (right) + '3: ar_l', ... % Acoustic radiation (left) + '4: ar_r', ... % Acoustic radiation (right) + '5: atr_l', ... % Anterior thalamic radiation (left) + '6: atr_r', ... % Anterior thalamic radiation (right) + '7: cbd_l', ... % Cingulum bundle – dorsal (left) + '8: cbd_r', ... % Cingulum bundle – dorsal (right) + '9: cbp_l', ... % Cingulum bundle – perigenual (left) + '10: cbp_r', ... % Cingulum bundle – perigenual (right) + '11: cbt_l', ... % Cingulum bundle – temporal (left) + '12: cbt_r', ... % Cingulum bundle – temporal (right) + '13: cst_l', ... % Corticospinal tract (left) + '14: cst_r', ... % Corticospinal tract (right) + '15: fa_l', ... % Frontal aslant tract (left) + '16: fa_r', ... % Frontal aslant tract (right) + '17: fma', ... % Forceps major (splenium of corpus callosum) + '18: fmi', ... % Forceps minor (genu of corpus callosum) + '19: fx_l', ... % Fornix (left) + '20: fx_r', ... % Fornix (right) + '21: ifo_l', ... % Inferior fronto-occipital fasciculus (left) + '22: ifo_r', ... % Inferior fronto-occipital fasciculus (right) + '23: ilf_l', ... % Inferior longitudinal fasciculus (left) + '24: ilf_r', ... % Inferior longitudinal fasciculus (right) + '25: mcp', ... % Middle cerebellar peduncle + '26: mdlf_l', ... % Middle longitudinal fasciculus (left) + '27: mdlf_r', ... % Middle longitudinal fasciculus (right) + '28: or_l', ... % Optic radiation (left) + '29: or_r', ... % Optic radiation (right) + '30: slf1_l', ... % Superior longitudinal fasciculus I (left) + '31: slf1_r', ... % Superior longitudinal fasciculus I (right) + '32: slf2_l', ... % Superior longitudinal fasciculus II (left) + '33: slf2_r', ... % Superior longitudinal fasciculus II (right) + '34: slf3_l', ... % Superior longitudinal fasciculus III (left) + '35: slf3_r', ... % Superior longitudinal fasciculus III (right) + '36: str_l', ... % Superior thalamic radiation (left) + '37: str_r', ... % Superior thalamic radiation (right) + '38: uf_l', ... % Uncinate fasciculus (left) + '39: uf_r', ... % Uncinate fasciculus (right) + '40: vof_l', ... % Vertical occipital fasciculus (left) + '41: vof_r' ... % Vertical occipital fasciculus (right) + }; + + + a = @(x) strsplit(x, ' '); + b = @(x,i) x{i}; + c = @(x) 1 + str2num(x(1:(end-1))); + + labels = cellfun(@(x) b(a(x),2), txt, 'UniformOutput', false); + ids = cellfun(@(x) c(b(a(x),1)), txt, 'UniformOutput', false); + + + end + end +end + diff --git a/examples/my_nodes/my_mask_trim.m b/examples/my_nodes/my_mask_trim.m new file mode 100644 index 0000000..1b331b0 --- /dev/null +++ b/examples/my_nodes/my_mask_trim.m @@ -0,0 +1,53 @@ +classdef my_mask_trim < dp_node + + methods + + function output = i2o(obj, input) + + output = input; + output.mask_fn = dp.new_fn(input.op, input.mask_fn, '_trim'); + + + output.fa_fn = dp.new_fn(input.op, input.fa_fn, '_trim'); + output.md_fn = dp.new_fn(input.op, input.md_fn, '_trim'); + output.rd_fn = dp.new_fn(input.op, input.rd_fn, '_trim'); + output.ad_fn = dp.new_fn(input.op, input.ad_fn, '_trim'); + output.s0_fn = dp.new_fn(input.op, input.s0_fn, '_trim'); + + + end + + function output = execute(obj, input, output) + + % Load, trim, save the mask + [M,h] = mdm_nii_read(input.mask_fn); + + f = @(x,a,b) mio_smooth_4d(double(x), a) > b; + M = f(f(f(M, 0.7, 0.3), 0.75, 0.8), 0.7, 0.85); + M = f(M, 0.7, 0.8); + + mdm_nii_write(double(M), output.mask_fn, h); + + % Trim parameter maps + + s = {'fa_fn', 'md_fn', 'ad_fn', 'rd_fn', 's0_fn'}; + + for c = 1:numel(s) + + [I,h] = mdm_nii_read(input.(s{c})); + IM = medfilt3(I); + + if (c == 1) + ind = mio_smooth_4d(double(I - IM), 0.5) > 0.15; + end + + I(ind) = IM(ind); + + mdm_nii_write(I.*M, output.(s{c}), h); + end + end + + + end + +end \ No newline at end of file diff --git a/examples/my_nodes/my_primary_node.m b/examples/my_nodes/my_primary_node.m new file mode 100644 index 0000000..cae0f5e --- /dev/null +++ b/examples/my_nodes/my_primary_node.m @@ -0,0 +1,44 @@ +classdef my_primary_node < dp_node_primary + + properties + bp; + pattern = '*'; + end + + methods + + function obj = my_primary_node(bp, pattern) + obj.bp = bp; + obj.pattern = pattern; + end + + function outputs = get_iterable(obj) + + d = dir(fullfile(obj.bp, obj.pattern)); + + outputs = {}; + + % find folders + for c = 1:numel(d) + + if (d(c).name(1) == '.') + continue; + end + + d2 = dir(fullfile(obj.bp, d(c).name, '2*')); + + if (numel(d2) > 1) + 1; + end + + for c2 = 1:numel(d2) + + output.bp = obj.bp; + output.id = fullfile(d(c).name, d2(c2).name); + + outputs{end+1} = output; + end + end + end + end +end \ No newline at end of file diff --git a/examples/my_nodes/split_rois.m b/examples/my_nodes/split_rois.m new file mode 100644 index 0000000..48fe6ab --- /dev/null +++ b/examples/my_nodes/split_rois.m @@ -0,0 +1,95 @@ +classdef split_rois < dp_node_roi_from_label + + properties + roi_list + end + + methods + + + function obj = split_rois(name, segm_node, roi_list) + obj = obj@dp_node_roi_from_label(name,segm_node); + obj.roi_list = roi_list; + + slab_suffix = {'_I','_II','_III'}; + + % find ROI ids from the segm node's vocabulary + base_ids = []; + for k = 1:numel(roi_list) + hit = find(strcmpi(obj.roi_names, roi_list{k}), 1); + if ~isempty(hit) + % extract numeric value if roi_ids is a cell array + if iscell(obj.roi_ids) + base_ids(end+1) = obj.roi_ids{hit}; + else + base_ids(end+1) = obj.roi_ids(hit); + end + else + warning('split_rois: tract %s not found', roi_list{k}); + end + end + + % assign three new IDs per base tract + next_id = max([obj.roi_ids{:}]) + 1; % unwrap cells if needed + + + for i = 1:numel(base_ids) + base_name = obj.roi_names{i}; % or your roi_list{i} + for s = 1:3 + obj.roi_ids{end+1} = next_id; + obj.roi_names{end+1} = [roi_list{i} slab_suffix{s}]; + next_id = next_id + 1; + end + end + + end + + function [R,h] = roi_get_volume(obj, output, ~, c_roi) + + rid = obj.roi_ids; + if iscell(rid), ridv = [rid{:}]; else, ridv = rid; end + ix = find(ridv == c_roi, 1); + + name = obj.roi_names{ix}; + + % --- is this one of our synthetic slabs? (… I / II / III) --- + if endsWith(name,'_I') || endsWith(name,'_II') || endsWith(name,'_III') + + % which slab is requested? + if endsWith(name,'_I'), slab = 1; + elseif endsWith(name,'_II'), slab = 2; + else slab = 3; + end + + % base tract name without the suffix + base_name = regexprep(name, '(_I|_II|_III)$', ''); + + % look up the base tract's ROI id by exact (case-insensitive) name + ix_base = find(strcmpi(obj.roi_names, base_name), 1); + + if iscell(rid), base_id = rid{ix_base}; else, base_id = rid(ix_base); end + + % fetch the base mask via the parent implementation + [R,h] = roi_get_volume@dp_node_roi_from_label(obj, output, 1, base_id); + + switch slab + case 1 + R(:, :, 44:end) = 0; + case 2 + R(:, :, 1:44) = 0; + R(:, :, 66:end) = 0; + case 3 + R(:, :, 1:66) = 0; + end + + else + [R,h] = roi_get_volume@dp_node_roi_from_label(obj, output, 1, c_roi); + end + + + end + + end + + +end \ No newline at end of file diff --git a/examples/step2_pipe_tractseg_xtract_juelich.m b/examples/step2_pipe_tractseg_xtract_juelich.m new file mode 100644 index 0000000..d9aaeca --- /dev/null +++ b/examples/step2_pipe_tractseg_xtract_juelich.m @@ -0,0 +1,211 @@ +% Root to NIfTI data +base_path = '/path/to/your/data/nii'; + +% Subjects → add dmri/xps/op +id1 = my_primary_node(base_path, '*'); +% id1 = dp_node_io_filter_by_number(1).connect(id1); % keep only first n subjects + +id2 = dp_node_io_append({... + {'dmri_fn', @(x) fullfile(x.bp, x.id, 'DWI', 'DTI_dn_mc.nii.gz')}, ... + {'xps_fn', @(x) fullfile(x.bp, x.id, 'DWI', 'DTI_dn_mc_xps.mat')}, ... + {'op', @(x) fullfile(x.bp, x.id, 'tractseg')}}).connect(id1); + +% b0 → mask +m0 = dp_node_dmri_subsample_b0().connect(id2); +m1 = dp_node_segm_hd_bet().connect(m0, 'm1'); + +% xps→bval/bvec, gradcheck +t0 = dp_node_dmri_io_xps_to_bval_bvec().connect(id2, 't0'); +t1 = dp_node_mrtrix_dwigradcheck().connect(t0); + +% merge mask + grads, alias mask +t2 = dp_node_io_merge({t1, m1}); t2.do_prefix = 0; +t3 = dp_node_io('mask_fn', 'm1_mask_fn').connect(t2); + +% TractSeg (MNI outputs: Diffusion_MNI, FA_MNI, nodif_brain_mask_MNI, bundles.nii.gz) +t4 = dp_node_segm_tractseg().connect(t3); + +% XTRACT definitions in a separate output folder +xt_op = dp_node_io('op', @(x) fullfile(x.bp, x.id, 'tractseg_xtract')).connect(t3); +xt1 = dp_node_segm_xtractseg().connect(xt_op); + +% xt1.run('report'); +% xt1.run('execute'); +% xt1.run('execute', struct('do_overwrite',1)); + +% DTI (LLS) on TractSeg’s MNI-space diffusion data (b_delta = 1) into subfolder /dti_lls +dti1 = dp_node_dmri_xps_from_bval_bvec(1).connect(t4); +dti2 = dp_node_io('op', @(x) fullfile(x.bp, x.id, 'tractseg', 'dti_lls')).connect(dti1); +dti3 = dp_node_dmri_dti().connect(dti2); + +% Refine brain mask and trim FA/MD/AD/RD/S0 maps using smoothing + median filtering +dti_trim = my_mask_trim().connect(dti3); + + +% Register AD/RD output filenames for dti_lls (not used downstream in this script) +% dti3_append = dp_node_io_append({ +% {'ad_fn', @(x) fullfile(x.op, 'dti_lls_ad.nii.gz')}, ... +% {'rd_fn', @(x) fullfile(x.op, 'dti_lls_rd.nii.gz')} +% }).connect(dti3); + + +% Merge trimmed DTI maps with TractSeg outputs (keep original field names) +mdt = dp_node_io_merge({dti_trim, t4}); mdt.do_prefix = 0; + +% Expose only labels + FA/MD/AD/RD for TractSeg ROI stats (used by t5, t5_split) +troi = dp_node_io_rename({ + {'labels_fn','dp_node_segm_tractseg_labels_fn'}, ... + {'fa_fn','fa_fn'}, ... + {'md_fn','md_fn'}, ... + {'ad_fn','ad_fn'}, ... + {'rd_fn','rd_fn'} ... +}).connect(mdt); + +% Merge trimmed DTI maps with XTRACT-based TractSeg outputs (keep original field names) +mdtx = dp_node_io_merge({dti_trim, xt1}); mdtx.do_prefix = 0; + +% Expose only labels + FA/MD/AD/RD for XTRACT ROI stats (used by xt2, xt1_split) +xtroi = dp_node_io_rename({ + {'labels_fn','dp_node_segm_xtractseg_labels_fn'}, ... + {'fa_fn','fa_fn'}, ... + {'md_fn','md_fn'}, ... + {'ad_fn','ad_fn'}, ... + {'rd_fn','rd_fn'} ... +}).connect(mdtx); + + +% TractSeg bundle ROI stats → CSV (mean + voxel count) +t5 = dp_node_roi_from_label('tractseg', t4).connect(troi); +t_csv = dp_node_csv('tractseg_export','../report/csv', {'mean','n'}).connect(t5); + +% TractSeg: split selected bundles into 3 z-slabs (I/II/III) and export ROI stats +ts_split_list = { ... + 'CST_left','CST_right', ... + 'SCP_left','SCP_right', ... + 'ICP_left','ICP_right', ... + 'CC_1','CC_2','CC_3','CC_4','CC_5','CC_6','CC_7', ... + 'MCP', ... + 'IFO_left','IFO_right', ... + 'ILF_left','ILF_right', ... + 'CG_left','CG_right', ... + 'UF_left','UF_right'}; + +% TractSeg: split selected bundles into 3 z-slabs (I/II/III) and export ROI stats → CSV (mean + voxel count) +t5_split = split_rois('tractseg_slab', t4, ts_split_list).connect(troi); +t6_split = dp_node_csv('tractseg_export_with_slabs','../report/csv', {'mean','n'}).connect(t5_split); + +% XTRACT-based TractSeg bundle ROI stats → CSV (mean + voxel count) +xt2 = dp_node_roi_from_label('tractseg_xtract', xt1).connect(xtroi); +xt3 = dp_node_csv('tractseg_xtract_export','../report/csv', {'mean','n'}).connect(xt2); + +% XTRACT: split CST (left/right) into 3 z-slabs and export ROI stats → CSV (mean + voxel count) +xt1_split = split_rois('tractseg_xtract_slab', xt1, {'cst_l','cst_r'}).connect(xtroi); +xt_csv = dp_node_csv('tractseg_xtract_export_with_slabs','../report/csv', {'mean','n'}).connect(xt1_split); + + +% ------------- JHU ROIs using ANTs (quick registration) ---------------------- +% +% TractSeg places dMRI in a semi-MNI rigid space (t4.*_MNI). +% DTI maps (dti3) are computed on that same grid. +% To align the JHU atlas with this subject grid: +% (A) Register subject FA (semi-MNI) → JHU FA template using ANTs quick mode +% (B) Apply the inverse warp to bring JHU labels → subject trimmed FA grid +% (C) Crop warped labels by the subject mask +% +% Result: JHU atlas labels and trimmed DTI maps are voxel-aligned on the semi-MNI (TractSeg) grid. + + +% JHU atlas paths +jhu_fa_fn = '/usr/local/fsl/data/atlases/JHU/JHU-ICBM-FA-1mm.nii.gz'; +jhu_labels_fn = '/usr/local/fsl/data/atlases/JHU/JHU-ICBM-labels-1mm.nii.gz'; + +% ----------- ANTs quick registration branch -------- +amni0 = dp_node_io('op', @(x) fullfile(x.bp, x.id, 'tractseg_jhu_ants_quick')).connect(dti_trim); +% amni0.run('report'); + +% ANTs: subject trimmed FA (moving) -> JHU FA (fixed) +amni1 = dp_node_io('nii_fn', 'fa_fn').connect(amni0); % moving: trimmed FA +% amni1.run('report'); + +amni2 = dp_node_io('target_fn', @(x) jhu_fa_fn).connect(amni1); % fixed: JHU FA +% amni2.run('report'); + +amni3 = dp_node_ants_reg_quick().connect(amni2); +% amni3.run('report'); + +% Attach trimmed FA grid as reference for antsApply +amni4 = dp_node_io_merge({amni3, dti_trim}); +amni4.do_prefix = 0; + +% Use trimmed FA as ref_fn +amni5 = dp_node_io_append({{'ref_fn','my_mask_trim_fa_fn'}}).connect(amni4); + +% Warp full JHU labels to subject trimmed FA space (using inverse chain) +amni6 = dp_node_io('nii_fn', @(x) jhu_labels_fn).connect(amni5); +amni7 = dp_node_ants_apply(' -n NearestNeighbor').connect(amni6); +% amni7.nii_fn: JHU labels (quick) on trimmed FA grid + +% Crop warped labels by subject mask on same grid +amni_merge1 = dp_node_io_merge({amni7, dti_trim}); amni_merge1.do_prefix = 0; +amni_merge2 = dp_node_io_append({{'mask_fn','my_mask_trim_fa_fn'}}).connect(amni_merge1); +jhu_crop = dp_node_atlas_crop_to_mask().connect(amni_merge2); + + +% QC +if (0) + qc1 = dp_node_io_merge({dti_trim, amni3, jhu_crop}); qc1.do_prefix = 0; + + % subject trimmed FA + labels warped back from MNI to subject space + if (0) + qc2 = dp_node_io_rename({ + {'fa_fn', 'fa_fn'}, ... + {'labels_fn', 'dp_node_atlas_crop_to_mask_nii_fn'}, ... + {'mask_fn', 'mask_fn'} + }).connect(qc1); + end + + % subject FA warped into JHU (MNI) space + labels in MNI + if (0) + qc2 = dp_node_io_rename({ + {'fa_fn', 'dp_node_ants_reg_quick_warped_fn'}, ... + {'labels_fn', @(x) jhu_labels_fn} + }).connect(qc1); + end + + % JHU template FA + JHU atlas labels (reference anatomy in MNI space) + if (0) + qc2 = dp_node_io_rename({ + {'fa_fn', @(x) jhu_fa_fn}, ... + {'labels_fn', @(x) jhu_labels_fn} + }).connect(qc1); + end + + % qc2.run('mgui'); + + qc3 = dp_node_roi_from_label('tmp', dp_node_segm_jhu()).connect(qc2); + qc3.run('mgui'); +end + + +% ****** Merge for ROI stats ****** + +mda = dp_node_io_merge({dti_trim, jhu_crop}); +mda.do_prefix = 0; + +aroi = dp_node_io_rename({ + {'labels_fn','dp_node_atlas_crop_to_mask_nii_fn'}, ... + {'fa_fn','fa_fn'}, ... + {'md_fn','md_fn'}, ... + {'ad_fn','ad_fn'}, ... + {'rd_fn','rd_fn'} +}).connect(mda); + + +jhu_LUT = dp_node_segm_jhu(); %LUT +jroi = dp_node_roi_from_label('jhu_ants_quick', jhu_LUT).connect(aroi); +j_cvs = dp_node_csv('jhu_ants_export','../report/csv', {'mean','n'}).connect(jroi); + + + + +