I am fairly new to AFNI and am having trouble with my alignment between my functional and my anatomical scans.

I have used different cost functions (lpc+Z, lpc, lpa), different outlier limits (0.3, 0.25, 0.2), different motion limits (0.3, 0.25), different volume registrations (minimal outlier, 3rd TR), and different alignment functions (giant move, ginormous move), deobliqued vs non-deobliqued anatomicals, however, my alignment isn’t great. The stats functional (REML) seems to be going outside the final anatomical (activation outside anatomical).

The anatomical (brain.nii) has been motion corrected and manually edited by our collaborators. They said that “the T1 image they collected is a multi-echo image with multiple TEs and a RMS of these TEs. The multiple TEs give slightly different tissue contrast which can be used by freesurfer to better delineate the structure boundaries (although our contact is not sure if they did this or if they just used the RMS with manual corrections for each image).”
The final anatomical does not have a skull, but I am also in possession of the MPRAGE T1 scans that have not been skull stripped (we also ran with these MPRAGE T1 scans and the same issue is occurring).

I also deobliqued the functional scans so that they would be centered.
Before I ran my script, the deobliqued functional scans were near the raw anatomical scan (brain.nii) in alignment, however, the functional scans were a little off from each other.

Final Output: When I look at my Volume Registered (3Dvolreg) outputs, they are poorly aligned with the final_anat even before looking at the statistical functional maps (REML). I’ve included a picture of the run 1 3Dvolreg. The other runs are relatively similar - they are just not 100% aligned with each other. The stats REML does not align with the anat_final. The activation seems to be going outside the anat_final.

What can I do to get the best alignment possible? I have included my script below. Any advice would help!



set list of runs

set runs = (count -digits 2 1 6)

create results and stimuli directories

mkdir $output_dir
mkdir $output_dir/stimuli

copy stim files into stimulus directory


copy anatomy to results dir


============================ auto block: tcat ============================

apply 3dTcat to copy input dsets to results dir, while

removing the first 0 TRs

3dTcat -prefix $output_dir/pb00.subj.r01.tcat \ /space/kmsyn03/projects/MGH_alc/Stroop/individual/SUBJ/Afni/SUBJ_stroop_RUN01_deobl+orig'[0..]’
3dTcat -prefix $output_dir/pb00.subj.r02.tcat \ /space/kmsyn03/projects/MGH_alc/Stroop/individual/SUBJ/Afni/SUBJ_stroop_RUN02_deobl+orig'[0..]’
3dTcat -prefix $output_dir/pb00.subj.r03.tcat \ /space/kmsyn03/projects/MGH_alc/Stroop/individual/SUBJ/Afni/SUBJ_stroop_RUN03_deobl+orig'[0..]’
3dTcat -prefix $output_dir/pb00.subj.r04.tcat \ /space/kmsyn03/projects/MGH_alc/Stroop/individual/SUBJ/Afni/SUBJ_stroop_RUN04_deobl+orig'[0..]’
3dTcat -prefix $output_dir/pb00.subj.r05.tcat \ /space/kmsyn03/projects/MGH_alc/Stroop/individual/SUBJ/Afni/SUBJ_stroop_RUN05_deobl+orig'[0..]’
3dTcat -prefix $output_dir/pb00.subj.r06.tcat \ /space/kmsyn03/projects/MGH_alc/Stroop/individual/SUBJ/Afni/SUBJ_stroop_RUN06_deobl+orig'[0..]’

and make note of repetitions (TRs) per run

set tr_counts = ( 146 146 146 146 146 146 )


enter the results directory (can begin processing data)

cd $output_dir

========================== auto block: outcount ==========================

data check: compute outlier fraction for each volume

touch out.pre_ss_warn.txt
foreach run ( $runs )
3dToutcount -automask -fraction -polort 3 -legendre
pb00.$subj.r$run.tcat+orig > outcount.r$run.1D

# censor outlier TRs per run, ignoring the first 0 TRs
# - censor when more than 0.3 of automask voxels are outliers
# - step() defines which TRs to remove via censoring
1deval -a outcount.r$run.1D -expr "1-step(a-0.3)" > rm.out.cen.r$run.1D

# outliers at TR 0 might suggest pre-steady state TRs
if ( `1deval -a outcount.r$run.1D"{0}" -expr "step(a-0.4)"` ) then
    echo "** TR #0 outliers: possible pre-steady state TRs in run $run" \
        >> out.pre_ss_warn.txt


catenate outlier counts into a single time series

cat outcount.r*.1D > outcount_rall.1D

catenate outlier censor files into a single time series

cat rm.out.cen.r*.1D > outcount_${subj}_censor.1D

get run number and TR index for minimum outlier volume

set minindex = 3dTstat -argmin -prefix - outcount_rall.1D\'
set ovals = ( -set_run_lengths $tr_counts \ -index_to_run_tr $minindex )

save run and TR indices for extraction of vr_base_min_outlier

set minoutrun = $ovals[1]
set minouttr = $ovals[2]
echo “min outlier: run $minoutrun, TR $minouttr” | tee out.min_outlier.txt

================================= tshift =================================

time shift data so all slice timing is the same

foreach run ( $runs )
3dTshift -tzero 0 -quintic -prefix pb01.$subj.r$run.tshift


extract volreg registration base

3dbucket -prefix vr_base_min_outlier

================================= align ==================================

for e2a: compute anat alignment transformation to EPI registration base

(new anat will be current brain+orig) -anat2epi -anat brain+orig
-suffix _al_junk
-epi vr_base_min_outlier+orig -epi_base 0
-epi_strip 3dAutomask
-anat_has_skull no
-cost lpc+ZZ -giant_move
-volreg off -tshift off

================================== tlrc ==================================

warp anatomy to standard space

@auto_tlrc -base MNI_avg152T1+tlrc -input brain+orig -no_ss

store forward transformation matrix in a text file

cat_matvec brain+tlrc::WARP_DATA -I > warp.anat.Xat.1D

================================= volreg =================================

align each dset to base volume, align to anat, warp to tlrc space

verify that we have a +tlrc warp dataset

if ( ! -f brain+tlrc.HEAD ) then
echo “** missing +tlrc warp dataset: brain+tlrc.HEAD”

register and warp

foreach run ( $runs )
# register each volume to the base
3dvolreg -verbose -zpad 1 -base vr_base_min_outlier+orig
-1Dfile dfile.r$run.1D -prefix rm.epi.volreg.r$run
-1Dmatrix_save mat.r$run.vr.aff12.1D

# create an all-1 dataset to mask the extents of the warp
3dcalc -overwrite -a pb01.$subj.r$run.tshift+orig -expr 1   \
       -prefix rm.epi.all1

# catenate volreg/epi2anat/tlrc xforms
cat_matvec -ONELINE                                         \
           brain+tlrc::WARP_DATA -I                         \
           brain_al_junk_mat.aff12.1D -I                    \
           mat.r$run.vr.aff12.1D > mat.r$run.warp.aff12.1D

# apply catenated xform: volreg/epi2anat/tlrc
3dAllineate -base brain+tlrc                                \
            -input pb01.$subj.r$run.tshift+orig             \
            -1Dmatrix_apply mat.r$run.warp.aff12.1D         \
            -mast_dxyz 4                                    \
            -prefix rm.epi.nomask.r$run

# warp the all-1 dataset for extents masking
3dAllineate -base brain+tlrc                                \
            -input rm.epi.all1+orig                         \
            -1Dmatrix_apply mat.r$run.warp.aff12.1D         \
            -mast_dxyz 4 -final NN -quiet                   \
            -prefix rm.epi.1.r$run

# make an extents intersection mask of this run
3dTstat -min -prefix rm.epi.min.r$run rm.epi.1.r$run+tlrc


make a single file of registration params

cat dfile.r*.1D > dfile_rall.1D


create the extents mask: mask_epi_extents+tlrc

(this is a mask of voxels that have valid data at every TR)

3dMean -datum short -prefix rm.epi.mean rm.epi.min.r*.HEAD
3dcalc -a rm.epi.mean+tlrc -expr ‘step(a-0.999)’ -prefix mask_epi_extents

and apply the extents mask to the EPI data

(delete any time series with missing data)

foreach run ( $runs )
3dcalc -a rm.epi.nomask.r$run+tlrc -b mask_epi_extents+tlrc
-expr ‘a*b’ -prefix pb02.$subj.r$run.volreg

warp the volreg base EPI dataset to make a final version

cat_matvec -ONELINE
brain+tlrc::WARP_DATA -I
brain_al_junk_mat.aff12.1D -I > mat.basewarp.aff12.1D

3dAllineate -base brain+tlrc
-input vr_base_min_outlier+orig
-1Dmatrix_apply mat.basewarp.aff12.1D
-mast_dxyz 4
-prefix final_epi_vr_base_min_outlier

create an anat_final dataset, aligned with stats

3dcopy brain+tlrc anat_final.$subj

record final registration costs

3dAllineate -base final_epi_vr_base_min_outlier+tlrc -allcostX
-input anat_final.$subj+tlrc |& tee out.allcostX.txt

================================== blur ==================================

blur each volume of each run

foreach run ( $runs )
3dmerge -1blur_fwhm 8.0 -doall -prefix pb03.$subj.r$run.blur

================================== mask ==================================

create ‘full_mask’ dataset (union mask)

foreach run ( $runs )
3dAutomask -dilate 1 -prefix rm.mask_r$run pb03.$subj.r$run.blur+tlrc

create union of inputs, output type is byte

3dmask_tool -inputs rm.mask_r*+tlrc.HEAD -union -prefix full_mask.$subj

---- create subject anatomy mask, mask_anat.$subj+tlrc ----

(resampled from tlrc anat)

3dresample -master full_mask.$subj+tlrc -input brain+tlrc
-prefix rm.resam.anat

convert to binary anat mask; fill gaps and holes

3dmask_tool -dilate_input 5 -5 -fill_holes -input rm.resam.anat+tlrc
-prefix mask_anat.$subj

compute overlaps between anat and EPI masks

3dABoverlap -no_automask full_mask.$subj+tlrc mask_anat.$subj+tlrc
|& tee out.mask_ae_overlap.txt

note Dice coefficient of masks, as well

3ddot -dodice full_mask.$subj+tlrc mask_anat.$subj+tlrc
|& tee out.mask_ae_dice.txt

---- create group anatomy mask, mask_group+tlrc ----

(resampled from tlrc base anat, MNI_avg152T1+tlrc)

3dresample -master full_mask.$subj+tlrc -prefix ./

convert to binary group mask; fill gaps and holes

3dmask_tool -dilate_input 5 -5 -fill_holes -input
-prefix mask_group

================================= scale ==================================

scale each voxel time series to have a mean of 100

(be sure no negatives creep in)

(subject to a range of [0,200])

foreach run ( $runs )
3dTstat -prefix rm.mean_r$run pb03.$subj.r$run.blur+tlrc
3dcalc -a pb03.$subj.r$run.blur+tlrc -b rm.mean_r$run+tlrc
-c mask_epi_extents+tlrc
-expr ‘c * min(200, a/b*100)*step(a)*step(b)’
-prefix pb04.$subj.r$run.scale

================================ regress =================================

compute de-meaned motion parameters (for use in regression) -infile dfile_rall.1D -set_nruns 6
-demean -write motion_demean.1D

compute motion parameter derivatives (for use in regression) -infile dfile_rall.1D -set_nruns 6
-derivative -demean -write motion_deriv.1D

create censor file motion_${subj}_censor.1D, for censoring motion -infile dfile_rall.1D -set_nruns 6
-show_censor_count -censor_prev_TR
-censor_motion 0.25 motion_${subj}

combine multiple censor files

1deval -a motion_{subj}_censor.1D -b outcount_{subj}censor.1D
-expr “a*b” > censor

note TRs that were not censored

set ktrs = -infile censor_${subj}_combined_2.1D \ -show_trs_uncensored encoded


run the regression analysis

3dDeconvolve -input pb04.subj.r*.scale+tlrc.HEAD \ -censor censor_{subj}_combined_2.1D
-polort 3
-num_stimts 16
-stim_times 1 stimuli/SUBJ_stimuli1.txt ‘TENT(-4.4,13.2,8)’
-stim_label 1 stimuli1
-stim_times 2 stimuli/SUBJ_stimuli2.txt ‘TENT(-4.4,13.2,8)’
-stim_label 2 stimuli2
-stim_times 3 stimuli/SUBJ_stimuli3.txt ‘TENT(-4.4,13.2,8)’
-stim_label 3 stimuli3
-stim_times 4 stimuli/SUBJ_stimuli4.txt ‘TENT(-4.4,13.2,8)’
-stim_label 4 stimuli4
-stim_file 5 motion_demean.1D’[0]’ -stim_base 5 -stim_label 5 roll_01
-stim_file 6 motion_demean.1D’[1]’ -stim_base 6 -stim_label 6 pitch_01
-stim_file 7 motion_demean.1D’[2]’ -stim_base 7 -stim_label 7 yaw_01
-stim_file 8 motion_demean.1D’[3]’ -stim_base 8 -stim_label 8 dS_01
-stim_file 9 motion_demean.1D’[4]’ -stim_base 9 -stim_label 9 dL_01
-stim_file 10 motion_demean.1D’[5]’ -stim_base 10 -stim_label 10 dP_01
-stim_file 11 motion_deriv.1D’[0]’ -stim_base 11 -stim_label 11 roll_02
-stim_file 12 motion_deriv.1D’[1]’ -stim_base 12 -stim_label 12 pitch_02
-stim_file 13 motion_deriv.1D’[2]’ -stim_base 13 -stim_label 13 yaw_02
-stim_file 14 motion_deriv.1D’[3]’ -stim_base 14 -stim_label 14 dS_02
-stim_file 15 motion_deriv.1D’[4]’ -stim_base 15 -stim_label 15 dL_02
-stim_file 16 motion_deriv.1D’[5]’ -stim_base 16 -stim_label 16 dP_02
-iresp 1 iresp_stimuli1.$subj
-iresp 2 iresp_stimuli2.$subj
-iresp 3 iresp_stimuli3.$subj
-iresp 4 iresp_stimuli4.$subj
-gltsym ‘SYM: stimuli2 -stimuli1’
-glt_label 1 I-C
-fout -tout -x1D X.xmat.1D -xjpeg X.jpg
-x1D_uncensored X.nocensor.xmat.1D
-fitts fitts.subj \ -errts errts.{subj}
-bucket stats.$subj

if 3dDeconvolve fails, terminate the script

if ( $status != 0 ) then
echo ‘---------------------------------------’
echo ‘** 3dDeconvolve error, failing…’
echo ’ (consider the file 3dDeconvolve.err)’

display any large pairwise correlations from the X-matrix -show_cormat_warnings -infile X.xmat.1D |& tee out.cormat_warn.txt

look for odd timing in files for TENT functions -multi_timing stimuli/SUBJ_stimuli1.txt
-tr 2.2 -warn_tr_stats |& tee out.TENT_warn.txt

– execute the 3dREMLfit script, written by 3dDeconvolve –

tcsh -x stats.REML_cmd

if 3dREMLfit fails, terminate the script

if ( $status != 0 ) then
echo ‘---------------------------------------’
echo ‘** 3dREMLfit error, failing…’

create an all_runs dataset to match the fitts, errts, etc.

3dTcat -prefix all_runs.$subj pb04.$subj.r*.scale+tlrc.HEAD


create a temporal signal to noise ratio dataset

signal: if ‘scale’ block, mean should be 100

noise : compute standard deviation of errts

3dTstat -mean -prefix rm.signal.all all_runs.$subj+tlrc"[ktrs]" 3dTstat -stdev -prefix rm.noise.all errts.{subj}_REML+tlrc"[$ktrs]"
3dcalc -a rm.signal.all+tlrc
-b rm.noise.all+tlrc
-c full_mask.$subj+tlrc
-expr ‘c*a/b’ -prefix TSNR.$subj


compute and store GCOR (global correlation average)

(sum of squares of global mean of unit errts)

3dTnorm -norm2 -prefix rm.errts.unit errts.${subj}_REML+tlrc
3dmaskave -quiet -mask full_mask.$subj+tlrc rm.errts.unit+tlrc
> gmean.errts.unit.1D
3dTstat -sos -prefix - gmean.errts.unit.1D' > out.gcor.1D
echo “-- GCOR = cat out.gcor.1D


compute correlation volume

(per voxel: average correlation across masked brain)

(now just dot product with average unit time series)

3dcalc -a rm.errts.unit+tlrc -b gmean.errts.unit.1D -expr ‘a*b’ -prefix rm.DP
3dTstat -sum -prefix corr_brain rm.DP+tlrc


compute sum of non-baseline regressors from the X-matrix

(use to get list of regressor colums)

set reg_cols = -infile X.nocensor.xmat.1D -show_indices_interest
3dTstat -sum -prefix sum_ideal.1D X.nocensor.xmat.1D"[$reg_cols]"

also, create a stimulus-only X-matrix, for easy review

1dcat X.nocensor.xmat.1D"[$reg_cols]" > X.stim.xmat.1D

============================ blur estimation =============================

compute blur estimates

touch blur_est.$subj.1D # start with empty file

create directory for ACF curve files

mkdir files_ACF

– estimate blur for each run in epits –

touch blur.epits.1D

restrict to uncensored TRs, per run

foreach run ( $runs )
set trs = -infile X.xmat.1D -show_trs_uncensored encoded \ -show_trs_run $run
if ( $trs == “” ) continue
3dFWHMx -detrend -mask full_mask.$subj+tlrc
-ACF files_ACF/out.3dFWHMx.ACF.epits.r$run.1D
all_runs.$subj+tlrc"[$trs]" >> blur.epits.1D

compute average FWHM blur (from every other row) and append

set blurs = ( 3dTstat -mean -prefix - blur.epits.1D'{0..$(2)}'\' )
echo average epits FWHM blurs: $blurs
echo “$blurs # epits FWHM blur estimates” >> blur_est.$subj.1D

compute average ACF blur (from every other row) and append

set blurs = ( 3dTstat -mean -prefix - blur.epits.1D'{1..$(2)}'\' )
echo average epits ACF blurs: $blurs
echo “$blurs # epits ACF blur estimates” >> blur_est.$subj.1D

– estimate blur for each run in errts –

touch blur.errts.1D

restrict to uncensored TRs, per run

foreach run ( $runs )
set trs = -infile X.xmat.1D -show_trs_uncensored encoded \ -show_trs_run $run
if ( $trs == “” ) continue
3dFWHMx -detrend -mask full_mask.$subj+tlrc
-ACF files_ACF/out.3dFWHMx.ACF.errts.r$run.1D
errts.${subj}+tlrc"[$trs]" >> blur.errts.1D

compute average FWHM blur (from every other row) and append

set blurs = ( 3dTstat -mean -prefix - blur.errts.1D'{0..$(2)}'\' )
echo average errts FWHM blurs: $blurs
echo “$blurs # errts FWHM blur estimates” >> blur_est.$subj.1D

compute average ACF blur (from every other row) and append

set blurs = ( 3dTstat -mean -prefix - blur.errts.1D'{1..$(2)}'\' )
echo average errts ACF blurs: $blurs
echo “$blurs # errts ACF blur estimates” >> blur_est.$subj.1D

– estimate blur for each run in err_reml –

touch blur.err_reml.1D

restrict to uncensored TRs, per run

foreach run ( $runs )
set trs = -infile X.xmat.1D -show_trs_uncensored encoded \ -show_trs_run $run
if ( $trs == “” ) continue
3dFWHMx -detrend -mask full_mask.$subj+tlrc
-ACF files_ACF/out.3dFWHMx.ACF.err_reml.r$run.1D
errts.${subj}_REML+tlrc"[$trs]" >> blur.err_reml.1D

compute average FWHM blur (from every other row) and append

set blurs = ( 3dTstat -mean -prefix - blur.err_reml.1D'{0..$(2)}'\' )
echo average err_reml FWHM blurs: $blurs
echo “$blurs # err_reml FWHM blur estimates” >> blur_est.$subj.1D

compute average ACF blur (from every other row) and append

set blurs = ( 3dTstat -mean -prefix - blur.err_reml.1D'{1..$(2)}'\' )
echo average err_reml ACF blurs: $blurs
echo “$blurs # err_reml ACF blur estimates” >> blur_est.$subj.1D

add 3dClustSim results as attributes to any stats dset

mkdir files_ClustSim

run Monte Carlo simulations using method ‘ACF’

set params = ( grep ACF blur_est.$subj.1D | tail -n 1 )
3dClustSim -both -mask full_mask.$subj+tlrc -acf $params[1-3]
-cmd 3dClustSim.ACF.cmd -prefix files_ClustSim/ClustSim.ACF

run 3drefit to attach 3dClustSim results to stats

set cmd = ( cat 3dClustSim.ACF.cmd )
$cmd stats.subj+tlrc stats.{subj}_REML+tlrc

================== auto block: generate review scripts ===================

generate a review script for the unprocessed EPI data -script @epi_review.$subj
-dsets pb00.$subj.r*.tcat+orig.HEAD

generate scripts to review single subject results

(try with defaults, but do not allow bad exit status) -mot_limit 0.25 -out_limit 0.3 -exit0

========================== auto block: finalize ==========================

remove temporary files

\rm -f rm.*

if the basic subject review script is here, run it

(want this to be the last text output)

if ( -e @ss_review_basic ) ./@ss_review_basic |& tee out.ss_review.$subj.txt

return to parent directory

cd …

echo “execution finished: date

Hi, Becca-

Hmmm, OK. One thing to note-- we usually don’t deoblique beforehand, and that gets taken care of in the script. Note that “deobliquing” can mean two things: removing the obliquity information, so we pretend that the oblique coordinates were the scanner coords to start with (done via “3drefit -deoblique …”); or, applying hte obliquity information, so the dset gets resampled and put into the scanner coords (“3dWarp -deoblique …”).

MIN_OUTLIER should generally be the best way to choose a reference volume for motion correction. Outlier limits shouldn’t affect alignment.

I don’t know about the anatomical vol without seeing it… So I can’t comment on cost function usage (for standard, adult human EPI and T1w anatomical, lpc+ZZ or lpc is the cost function of choice-- the former is a bit more robust to weird things happening).

In general, including the command itself, rather than the script it creates, is a bit more concise and easier for us to comment on.

If you are willing to upload the data and your command, I can take a look at it and try running it; I will PM you instructions for that.