Tractography between cortical ROIs

Good afternoon, Paul, I hope you are well.

Last week I posted a question about making cortical ROI masks for tractography, which you helped me answer, thank you. I am trying to isolate specific WM tracts such as the uncinate fasciculus, SLF, ILF, IFOF etc. I am doing this by defining the cortical start and end points of these tracts and creating masks so that tractography is only done between these ROI start and end points. Here is how I'm doing it:

#!/bin/tcsh

#UNCINATE FASCICULUS MASK
3dresample                                                                   \
    -master     $BASE/$SUBJ/INTERMED/dt_DT.nii.gz                            \
    -inset      $BASE/$SUBJ/INTERMED/aparc+aseg.nii.gz                       \
    -rmode      NN                                                           \
    -prefix     $BASE/$SUBJ/INTERMED/aparc+aseg_dwi.nii.gz                   \
    -overwrite

#Make mask temporal (LHS)
3dcalc                                                                       \
    -a          $BASE/$SUBJ/INTERMED/aparc+aseg_dwi.nii.gz                   \
    -expr       'equals(a,1033)+equals(a,1015)+equals(a,1009)'               \
    -prefix     $BASE/$SUBJ/INTERMED/UF_L_temporal_mask.nii.gz               \
    -overwrite

#Make mask frontal (LHS)
3dcalc                                                                       \
    -a          $BASE/$SUBJ/INTERMED/aparc+aseg_dwi.nii.gz                   \
    -expr       'equals(a,1014)+equals(a,1012)'                              \
    -prefix     $BASE/$SUBJ/INTERMED/UF_L_frontal_mask.nii.gz                \
    -overwrite

#Inflate temporal (LHS)
3dROIMaker                                                                   \
    -inset    $BASE/$SUBJ/INTERMED/UF_L_temporal_mask.nii.gz                 \
    -refset   $BASE/$SUBJ/INTERMED/UF_L_temporal_mask.nii.gz                 \
    -thresh   0.5                                                            \
    -prefix   $BASE/$SUBJ/INTERMED/ROI_MAP_temporal                          \
    -volthr   1                                                              \
    -inflate  1

#Inflate frontal (LHS)
3dROIMaker                                                                   \
    -inset    $BASE/$SUBJ/INTERMED/UF_L_frontal_mask.nii.gz                  \
    -refset   $BASE/$SUBJ/INTERMED/UF_L_frontal_mask.nii.gz                  \
    -thresh   0.5                                                            \
    -prefix   $BASE/$SUBJ/INTERMED/ROI_MAP_frontal                           \
    -volthr   1                                                              \
    -inflate  1

#Combine temporal and frontal (LHS)
3dcalc                                                                       \
    -a       $BASE/$SUBJ/INTERMED/ROI_MAP_temporal_L_GMI+orig.BRIK.gz        \
    -b       $BASE/$SUBJ/INTERMED/ROI_MAP_frontal_L_GMI+orig.BRIK.gz         \
    -expr    'a*1 + b*2'                                                     \
    -prefix  $BASE/$SUBJ/INTERMED/ROI_network.nii.gz

### for presentation purposes
time 3dTrackID                                                               \
    -logic            AND                                                    \
    -mode             MINIP                                                  \
    -alg_Nseed_X      3                                                      \
    -alg_Nseed_Y      3                                                      \
    -alg_Nseed_Z      3                                                      \
    -dti_in           $BASE/$SUBJ/INTERMED/dt                                \
    -netrois          $BASE/$SUBJ/INTERMED/ROI_network.nii.gz                \
    -mini_num         5                                                      \
    -uncert           $BASE/$SUBJ/INTERMED/dt_UNC.nii.gz                     \
    -mask             $BASE/$SUBJ/INTERMED/DTI_mask.nii.gz                   \
    -dump_rois        AFNI                                                   \
    -nifti                                                                   \
    -no_indipair_out                                                         \
    -prefix           $BASE/$SUBJ/INTERMED/o.mPR_L                           \
    -do_trk_out                                                              \
    -overwrite        -echo_edu


time 3dTrackID                                                               \
    -mode             PROB                                                   \
    -netrois          $BASE/$SUBJ/INTERMED/ROI_network.nii.gz                \
    -uncert           $BASE/$SUBJ/INTERMED/dt_UNC.nii.gz                     \
    -dti_in           $BASE/$SUBJ/INTERMED/dt                                \
    -mask             $BASE/$SUBJ/INTERMED/DTI_mask.nii.gz                   \
    -algopt           /scratch/brtphi011/control_trial/DTI_scripts_cluster/ALGOPTS_PROB.dat \
    -dump_rois        AFNI                                                   \
    -nifti                                                                   \
    -no_indipair_out                                                         \
    -prefix           $BASE/$SUBJ/INTERMED/o.PR_L                            \
    -overwrite        -echo_edu

What I want to ask is whether this is the best way to be doing this? I am working with neonate data and have my cortical parcellations (aparc+aseg.nii.gz) from infant freesurfer. The isolated tracts I am getting look generally consistent with the literature but there is just no way for me to quality check my results.

Please could you let me know if you know of any ways I could improve my method?

Warm regards,
Philippa

Hi, Philippa-

Those steps look reasonable, just as long as some assumptions are true.

For example, for the 3dresample, that assumes that the aparc+aseg.nii.gz dataset is in the same space as the dt_DT.nii.gz dataset. That is, the T1w volume from which the aparc segmentation was estimated overlaps (very) well with the DTI dataset, like having been aligned already with 3dAllineate or something. Is that the case? I'll note that the AFNI program fat_proc_map_to_dti exists actually for specifically the purpose of aligning a T1w volume to a DTI b0 dataset, and being able to send along ROI maps with it. An example of that is here. An advantage of using the program to align the T1w volume to the DTI b0 is that it comes "pre set" to align volumes with opposite contrast, and it will produce QC images automatically to check (as that link shows). But it isn't necessary to use this---doing so would just combine the alignment and ROI resampling in a single step.

The 3dcalc commands are fine for selecting out ROIs. Note that each one is currently taking multiple ROIs and merging them into a single ROI, Like, ROIs #1033, 1015 and 1009 in the original aparc dataset become merged into a single ROI. Is that intended? That is certainly fine to do.

From a tracking point of view, you will eventually want your N ROIs of interest in the same volume, each labelled with a different integer in the volume. So, if you wanted to, you could do what you are doing now, with a single 3dcalc command and what appear to be your 2 ROIs of interest in a single volume from the start; there are many ways to select the ROIs and combine, and yours is correct, but I will choose a terser version where the ROI values are specified at the time of reading the file in:

# Make mask ALL: 
#  LHS_temporal = 1
#  LHS_frontal = 2
3dcalc                                                                       \
    -a          $BASE/$SUBJ/INTERMED/aparc+aseg_dwi.nii.gz"<1033,1015,1009>"                   \
    -b          $BASE/$SUBJ/INTERMED/aparc+aseg_dwi.nii.gz"<1014,1012>"                   \
    -expr       '1*bool(a) + 2*bool(b)'               \
    -prefix     $BASE/$SUBJ/INTERMED/UF_L_ALL_mask.nii.gz               \
    -overwrite

Then, only a single 3dROIMaker command is needed. There is a benefit to this in general: before inflation, the ROIs created in this way are guaranteed not to overlap. When you run 3dROIMaker, the program is built to also not allow them to overlap as they inflate. However, if you inflate first and then combine, there is a chance that the ROIs could overlap when combined with the later 3dcalc command, which would be awkward. I guess your ROIs are far apart, but in case you adapt this later for other ROIs or share the code with a colleague, it is better to avoid potential woe So, I think it is safer to keep all your non-overlapping ROIs in a file first, and then run 3dROIMaker. Continuing from the above 3dROIMaker command, you would then run:

3dROIMaker                                                                   \
    -inset    $BASE/$SUBJ/INTERMED/UF_L_ALL_mask.nii.gz                  \
    -refset   $BASE/$SUBJ/INTERMED/UF_L_ALL_mask.nii.gz                  \
    -thresh   0.5                                                            \
    -prefix   $BASE/$SUBJ/INTERMED/ROI_MAP_ALL                          \
    -volthr   1                                                              \
    -inflate  1

Note that during this step, you could also include the DTI-based FA map, with a threshold value, to constrain inflation into the WM; but since you are just inflating by 1 voxel, that might not be necessary.

At this point, you would no longer need the final 3dcalc command to combine ROIs (the "Combine temporal and frontal (LHS)" one). And you could use the $BASE/$SUBJ/INTERMED/ROI_MAP_ALL* dset for tracking.

And the 3dTrackID commands look good to me.

--pt

Oh, and just in case it is useful, a couple of visualization points/helper programs.

You can create DEC maps of your DTI fit results, which is a nice Red-Blue-Green colorized way to see how the tensors looks; see here for example of using fat_proc_decmap.

I noticed that you are running the mini-probabilistic tracking for visualization purposes, and I would do that, too. But I would also visualize the probability WM tract regions using fat_proc_connec_vis, as done here.

--pt

Hi Paul, thank you for your response.

I have not used 3dAllineate to align my T1w volume and my DWI dataset. Should I do this before I run Infant Freesurfer (which is the very first step of my whole pipeline)? Infant Freesurfer produces the aparc+aseg.nii.gz file from an mprage.nii.gz T1. Or can I just align the aparc+aseg.nii.gz file to the dt_DT.nii.gz at this stage?

Another thing I was unsure about is whether I need to run 3dAllineate before fat_proc_map_to_dti? Or does running fat_proc_map_to_dti replace 3dAllineate?

I ran fat_proc_map_to_dti without 3dAllineate:

fat_proc_map_to_dti \
        -source $BASE/$SUBJ/ANATOM/mprage.nii.gz  \
        -followers_NN $BASE/$SUBJ/INTERMED/aparc+aseg.nii.gz  \
        -base $BASE/$SUBJ/INTERMED/FILT_AP/FILT_AP.nii'[0]' \
        -prefix $BASE/$SUBJ/INTERMED/fp_map2dti/fp_dwi

The mprage.nii.gz file is what I input into Infant Freesurfer to get the aparc+aseg.nii.gz file. The FILT_AP.nii file is the raw DWI AP dataset after being visually "filtered" for dropout signals and motion artefacts. These volumes were left out. Should I be using the AP or the PA scan? Or should I use combined dataset, after it has been processed by TORTOISE DIFFPREP & DRBUDDI (although now it is called TORTOISEProcess with the new version 4 of TORTOISE)? My quality check files do not look good at all...

Warm regards,
Philippa

Hi, Philippa-

The alignment here would be run after Infant FreeSurfer. It will be important/necessary to have the anatomical parcellation ROIs get to the space where the DWI/DTI data is, to be able to use them as tracking targets.

For the alignment, you can use fat_proc_map_to_dti instead of 3dAllineate. In actuality, fat_proc_map_to_dti is a wrapper for 3dAllineate with different (hopefully useful) pieces added in.

Judging by the name of the QC image you shared, that is just the initial overlap of the "base" and "source" volumes in your fat_proc_map_to_dti command; that is, the overlap of the first AP volume and the anatomical mprage, respectively. Alignment would proceed from that point.

Now, a subtle point is that I think it looks like at least one of those volumes has obliquity present? I will guess it is the mprage; can you please verify that with:

3dinfo -is_oblique $BASE/$SUBJ/ANATOM/mprage.nii.gz  

(In the output: 1 -> yes, 0 -> no). Different software deal with obliquity differently. The "regular" FreeSurfer applies it, so that the datasets in the output directory are in the "scanner" coordinates. As a result, you should really use the anatomical in the (infant) FreeSurfer results as the "source" for fat_proc_map_to_dti, because that is a better grid+space match for the aparc dataset. In the regular FreeSurfer output, this would be named like *_SurfVol.*.

On the DWI side, to decide what "base" volume to use in fat_proc_map_to_dti, you will want to use a b=0 volume in the final DWI space, so in the output from TORTOISEProcess. That should be as distortion corrected as processing can make it. When using TORTOISE, I will guess you use a T2w volume as an anatomical reference? If so, that will probably be called "struct*" in the output, and it should be on the same space+grid as the final DWI data. If there is a fair amount of brightness inhomogeneity in the DWI b=0 volume, then using the struct*.nii* dset that is on the same grid+space, and which has a T2w contrast, would also be fine as the base.

Yet another subtle point is:

  • With adult data, the T2w volume or b=0 DWI volume have opposite contrast to the associated T1w volume (that is, ventricles are dark in one but bright in another, etc.); therefore, the default cost function in fat_proc_map_to_dti is "lpc", which is good at doing alignment under those circumstances.
  • With infant data, we should check if that is still the case, because depending on the infant age, the myelination and other factors can affect the relative tissue contrast of dsets.
  • So, make sure to check whether the "base" and "source" volumes you are using in fat_proc_map_to_dti have similar or opposite contrasts. If they are opposite, then try running the command as discussed above. If they are similar, then you would need to change the cost function, such as using -cost lpa as a first try; if that fails, and/or if there is a lot of brightness inhomogeneity, then maybe trying -cost nmi would be good.

Please let me know how that goes.

--pt

Thank you very much Paul, this has helped a lot!

Kind regards,
Philippa

Good morning, Paul, I hope you are well.

This is a slightly unrelated issue that I am encountering, but I thought I should send it in this series so that you understand my pipeline context

I am having an issue with 3dTrackID. I have aligned my T1 and DWI, and am trying to track the SLF.

BASE=/scratch/brtphi011/fasd_trial/pippa_DTI_subjects
SING_IMAGE=/opt/exp_soft/singularity-containers/freesurfer/freesurfer-8.1.0.sif

# List of subjects

SUBJECTS=(2003 2024 2036 2044 2056 2071 2078 2096 2102 2114 2125 2143 2150 2165 2171 2214 2228 2233 2242 2256 2263 2008 2029 2039 2046 2064 \
 2072 2081 2097 2103 2117 2132 2144 2151 2166 2176 2217 2229 2236 2243 2257 2264 2009 2031 2040 2050 2068 2074 2084 2099 2107 2123 2134 \
 2146 2158 2167 2212 2218 2230 2238 2252 2259 2020 2033 2041 2052 2070 2076 2090 2100 2109 2124 2137 2149 2159 2170 2213 2227 2231 2240 \
 2254 2260)

# ==========================
# Select subject for this array task
# ==========================
SUBJ=${SUBJECTS[$SLURM_ARRAY_TASK_ID]}
BASE_ANAT=$BASE/$SUBJ/ANATOM

echo "Processing subject $SUBJ"
##########################################################################################################################################
#SLF MASK

#Make mask ALL:
#start = 1
#end =2
3dcalc                  \
        -a $BASE/$SUBJ/INTERMED/fp_map2dti/fp_dwi_aparc+aseg.nii.gz"<1029,1031>"        \
        -b $BASE/$SUBJ/INTERMED/fp_map2dti/fp_dwi_aparc+aseg.nii.gz"<1027,2018,2020>"   \
        -expr '1*bool(a) + 2*bool(b)'   \
        -prefix $BASE/$SUBJ/INTERMED/SLF_MASK_LHS.nii.gz           \
        -overwrite

rm -r $BASE/$SUBJ/INTERMED/SLF_ROI_MASK_LHS_GMI+orig.BRIK.gz
rm -r $BASE/$SUBJ/INTERMED/SLF_ROI_MASK_LHS_GMI+orig.HEAD
rm -r $BASE/$SUBJ/INTERMED/SLF_ROI_MASK_LHS_GM+orig.BRIK.gz
rm -r $BASE/$SUBJ/INTERMED/SLF_ROI_MASK_LHS_GM+orig.HEAD
rm -r $BASE/$SUBJ/INTERMED/SLF_ROI_MASK_LHS_GMI.niml.lt
rm -r $BASE/$SUBJ/INTERMED/SLF_ROI_MASK_LHS_GM.niml.lt

3dROIMaker              \
        -inset $BASE/$SUBJ/INTERMED/SLF_MASK_LHS.nii.gz \
        -refset $BASE/$SUBJ/INTERMED/SLF_MASK_LHS.nii.gz        \
        -thresh 0.5     \
        -prefix $BASE/$SUBJ/INTERMED/SLF_ROI_MASK_LHS   \
        -volthr 1	\
        -inflate 3

### for presentation purposes
time 3dTrackID   \
         -logic AND \
         -mode MINIP	  \
         -alg_Nseed_X 6       \
         -alg_Nseed_Y 6      \
         -alg_Nseed_Z 6     \
         -dti_in $BASE/$SUBJ/INTERMED/DTI/dt                                            \
         -netrois $BASE/$SUBJ/INTERMED/SLF_ROI_MASK_LHS_GMI+orig.BRIK.gz \
         -mini_num 5 \
         -uncert $BASE/$SUBJ/INTERMED/DTI/dt_UNC.nii.gz   \
         -mask $BASE/$SUBJ/INTERMED/DTI/DTI_mask.nii.gz                                 \
         -dump_rois AFNI    \
         -nifti                            \
         -no_indipair_out                     \
         -alg_Thresh_FA 0.1 \                                               \
         -prefix /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/$SUBJ/INTERMED/o.mPR_SLF_LHS  \
         -overwrite	  \
         -echo_edu

But I keep getting the following error about an illegal file name for all of my subjects, and I cannot work out why it is giving this error. PROB tractography is not having any issues with the file name outputs.

Processing subject 2024
++ 3dcalc: AFNI version=AFNI_26.0.01 (Jan  9 2026) [64-bit]
++ Authored by: A cast of thousands
++ Output dataset /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/SLF_MASK_LHS.nii.gz
++ Value of threshold is: 0.500000
++ Each subrik of refset will be applied to corresponding inset brik.
++ GM map is done.
++ SKEL_STOP = 0
++ No refset labeltable for naming things.
++ GMI map is done.
++ Put labeltables on both the GM and GMI files.
++ ROI logic type is: AND
++ Tracking mode: MINIP
++ Number of ROIs in netw[0] = 2
++ Have labeltable for naming things.
++ Running with 5 mini-prob iterations
++ SEARCHING for vector files with prefix '/scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/DTI/dt*':
	FINDING:	 'V1'  'V2'  'V3' 
++ SEARCHING for scalar files with prefix '/scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/DTI/dt*':
	FINDING:	 not:'CHI'  not:'DT'  'FA'  'L1'  not:'L2'  not:'L3'  'MD'  'RD'  not:'UNC'  not:'cwts'  not:'debugbriks' 
++ Done with scalar search, found: 4 parameters
	--> so will have 15 output data matrices.
++ Tracking progress count: start ...
	[  1/5  ] -> 0.15 min
	[  2/5  ] -> 0.32 min
	[  3/5  ] -> 0.48 min
	[  4/5  ] -> 0.63 min
	[  5/5  ] -> 0.80 min
++ Done tracking, tidying up outputs...
++ From tracking, net[0] has 5810 tracks.
++ Writing output (RPI, same as your input): (null) ...
++ Number of pairwise connections (bundles) in netw[0] = 1
++ Not writing out *INDI* and *PAIR* maps.
++ Writing out individual files to: (null)/NET_ ...
** ERROR: Illegal filename for NIfTI output: (null)/NET_000_ROI_001__001.nii.gz
** ERROR: Illegal filename for NIfTI output: (null)/NET_000_ROI_001__002.nii.gz
** ERROR: Illegal filename for NIfTI output: (null)/NET_000_ROI_002__002.nii.gz

+++ Command Echo:
   3dTrackID -logic AND -mode MINIP -alg_Nseed_X 6 -alg_Nseed_Y 6 -alg_Nseed_Z 6 -dti_in /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/DTI/dt -netrois /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/SLF_ROI_MASK_LHS_GMI+orig.BRIK.gz -mini_num 5 -uncert /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/DTI/dt_UNC.nii.gz -mask /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/DTI/DTI_mask.nii.gz -dump_rois AFNI -nifti -no_indipair_out -alg_Thresh_FA 0.1   -prefix /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/o.mPR_SLF_LHS -overwrite 



real	0m48.451s
user	0m48.240s
sys	0m0.067s
++ Tracking mode: PROB
++ Number of ROIs in netw[0] = 2
++ Have labeltable for naming things.
++ SEARCHING for vector files with prefix '/scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/DTI/dt*':
	FINDING:	 'V1'  'V2'  'V3' 
++ SEARCHING for scalar files with prefix '/scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/DTI/dt*':
	FINDING:	 not:'CHI'  not:'DT'  'FA'  'L1'  not:'L2'  not:'L3'  'MD'  'RD'  not:'UNC'  not:'cwts'  not:'debugbriks' 
++ Done with scalar search, found: 4 parameters
	--> so will have 15 output data matrices.
++ Effective Monte iterations: 5000. Fraction threshold set: 0.10000
	--> Ntrack voxel threshold: 500.
++ Tracking progress count: start ...
	[  10% ] -> 0.48 min
	[  20% ] -> 0.95 min
	[  30% ] -> 1.43 min
	[  40% ] -> 1.92 min
	[  50% ] -> 2.38 min
	[  60% ] -> 2.87 min
	[  70% ] -> 3.35 min
	[  80% ] -> 3.83 min
	[  90% ] -> 4.30 min
	[ 100% ] -> 4.78 min
++ Writing output (RPI, same as your input): /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/o.PR_SLF_LHS ...
++ Number of pairwise connections (bundles) in netw[0] = 1
++ Not writing out *INDI* and *PAIR* maps.
++ Writing out individual files to: /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/o.PR_SLF_LHS/NET_ ...

+++ Command Echo:
   3dTrackID -mode PROB -netrois /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/SLF_ROI_MASK_LHS_GMI+orig.BRIK.gz -uncert /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/DTI/dt_UNC.nii.gz -dti_in /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/DTI/dt -mask /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/DTI/DTI_mask.nii.gz -algopt /scratch/brtphi011/fasd_trial/DTI_scripts_pippa/ALGOPTS_PROB.dat -dump_rois AFNI -nifti -no_indipair_out -prefix /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/2024/INTERMED/o.PR_SLF_LHS -overwrite 



real	4m47.623s
user	4m46.870s
sys	0m0.055s

When tracking the uncinate fasciculus, I do not get this error, with the exact same code block for 3dTrackID.

Please let me know how I should proceed.
Thanks so much,
Warm regards,
Philippa

Hi, Philippa-

Hmm, OK, this is a subtle issue about scripting. When I saw that error, I thought, "Oh, there must not be '-prefix ..' option used." But then I saw you using one in your command. But then when I copied your command and tried running it, just replacing your filenames with ones in the FATCAT demo, I got the same error!

... and then I noticed this line where the FA threshold is given:

         -alg_Thresh_FA 0.1 \                                               \

The subtle, scripty thing here is that there are two backslashes there, not just one, to be the continuation of line character. Since there is whitespace to the right of the first one, that essentially ends the continuation-of-line-ness for the command. Hence, 3dTrackID does not see any options after that point, which includes -prefix ... So, indeed, for a subtle scripty reason, the command does not have a prefix to write out with, and when it tries to write, it is very unhappy, hence the errors. So, if you fix that line by just having a single \ with no whitespace to the right, your command should be in good shape.

(NB: the 3dTrackID program should actually whine much earlier on if no -prefix .. is provided, and so I will update it to do so, which would have helped make this issue clearer, too.)

If I can take the liberty of mentioning two other small considerations:

  • The total number of seeds per voxel for tracking is very large. Using this combo of options:

             -alg_Nseed_X 6       \
            -alg_Nseed_Y 6      \
             -alg_Nseed_Z 6     \
    

    .... means there are 6x6x6 = 216 seeds per voxel being run. I would probably start with something much smaller, say, using 2 in each of those cases, and having 2x2x2=8 seeds per voxel should still be plenty. Maybe it isn't that big a deal with your current use case, because you are tracking a sparse network, but for larger numbers of ROIs or larger mask or smaller voxels that could lead to huge output.

  • In lines 31-36, you are removing individual files with the Linux command rm. Certainly fine. But I would not include the -r option here---that means to make the removal recursive. Here, you are just 'recursing' through individual file names. Removing the -r won't change that behavior, but it will make your code safer against deleting a ton of things (like the full contents of a directory, including any subdirectories) if a scripting mistake were made. So, this isn't wrong, it just might prevent future woe if a typo slips in anywhere.

  • On an extremely minor note: in most cases, you have defined input and output files paths using variables, like in $BASE/$SUBJ/INTERMED. that is great, and what I would do, too. However, when making the path for outputting 3dTrackID's results, those variables are abandoned for a longer form, hardwired path: /scratch/brtphi011/fasd_trial/pippa_DTI_subjects/$SUBJ/INTERMED. If it were me, I think I would aim to keep the same style (and shorter pathname), and use that option with: -prefix $BASE/$SUBJ/INTERMED/o.mPR_SLF_LHS.

    • and actually, since that trio of path segments is used so often, I might be tempted to define a new path to the subject directory to stand in place of all of them, like: SDIR=$BASE/$SUBJ/INTERMED, and then you could just use $SDIR as a briefer form in paths.

But things look generally good there!

--pt