AFNI version info (afni -ver): Precompiled binary macos_13_ARM_clang: Jan 12 2025 (Version AFNI_24.3.10 'Elagabalus')
Hi, I’m using afni_proc.py for preprocessing for the first time. I specified -anat_has_skull yes, but I don’t see 3dSkullStrip or @SSwarper in the proc_subj script.
How does AFNI perform skull stripping in this case? Additionally, I noticed that auto_warp.py includes -skull_strip_input no. Shouldn’t this be yes if skull stripping is expected?
Thanks for your help!
afni_proc.py:
afni_proc.py \
-subj_id ${subj} \
-script proc_${subj} \
-out_dir ${dir_preproc} \
-blocks despike tshift align tlrc volreg blur mask scale regress \
-copy_anat ${dir_raw}/T1.${subj}.nii \
-anat_has_skull yes \
-anat_uniform_method unifize \
-anat_unif_GM yes \
-dsets ${dir_raw}/func.${subj}.r01.nii ${dir_raw}/func.${subj}.r02.nii ${dir_raw}/func.${subj}.r03.nii ${dir_raw}/func.${subj}.r04.nii ${dir_raw}/func.${subj}.r05.nii ${dir_raw}/func.${subj}.r06.nii \
-radial_correlate_blocks tcat volreg \
-blip_reverse_dset ${dir_raw}/dist_PA.${subj}.nii \
-tlrc_base MNI152_2009_template_SSW.nii.gz \
-tlrc_NL_warp \
-align_opts_aea \
-cost lpc+ZZ \
-giant_move \
-check_flip \
-volreg_align_to MIN_OUTLIER \
-volreg_align_e2a \
-volreg_tlrc_warp \
-blur_size 4.0 \
-mask_epi_anat yes \
-regress_motion_per_run \
-regress_censor_motion 0.3 \
-regress_censor_outliers 0.05 \
-regress_apply_mot_types demean deriv \
-regress_est_blur_epits \
-regress_est_blur_errts \
proc_subj:
#!/bin/tcsh -xef
echo "auto-generated by afni_proc.py, Fri Mar 7 13:52:41 2025"
echo "(version 7.82, December 6, 2024)"
echo "execution started: `date`"
# to execute via tcsh:
# tcsh -xef proc_D003 |& tee output.proc_D003
# to execute via bash:
# tcsh -xef proc_D003 2>&1 | tee output.proc_D003
# =========================== auto block: setup ============================
# script setup
# take note of the AFNI version
afni -ver
# check that the current AFNI version is recent enough
afni_history -check_date 7 Mar 2024
if ( $status ) then
echo "** this script requires newer AFNI binaries (than 7 Mar 2024)"
echo " (consider: @update.afni.binaries -defaults)"
exit
endif
# the user may specify a single subject to run with
if ( $#argv > 0 ) then
set subj = $argv[1]
else
set subj = D003
endif
# assign output directory name
set output_dir = /Users/andy/Desktop/Research/LDT/fMRI_data/preproc_afni/D003
# verify that the results directory does not yet exist
if ( -d $output_dir ) then
echo output dir "$subj.results" already exists
exit
endif
# set list of runs
set runs = (`count_afni -digits 2 1 6`)
# create results and stimuli directories
mkdir -p $output_dir
mkdir $output_dir/stimuli
# copy anatomy to results dir
3dcopy /Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/T1.D003.nii \
$output_dir/T1.D003
# copy template to results dir (for QC)
3dcopy /Users/andy/abin/MNI152_2009_template_SSW.nii.gz \
$output_dir/MNI152_2009_template_SSW.nii.gz
# copy any -ROI_import datasets as ROI_import_LABEL
3dcopy /Users/andy/abin/APQC_atlas_MNI_2009c_asym.nii.gz \
$output_dir/ROI_import_MNI_2009c_asym
# will extract automatic -blip_forward_dset in tcat block, below
# copy external -blip_reverse_dset dataset
3dTcat -prefix $output_dir/blip_reverse \
/Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/dist_PA.D003.nii
# ============================ 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 \
/Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/func.D003.r01.nii'[0..$]'
3dTcat -prefix $output_dir/pb00.$subj.r02.tcat \
/Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/func.D003.r02.nii'[0..$]'
3dTcat -prefix $output_dir/pb00.$subj.r03.tcat \
/Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/func.D003.r03.nii'[0..$]'
3dTcat -prefix $output_dir/pb00.$subj.r04.tcat \
/Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/func.D003.r04.nii'[0..$]'
3dTcat -prefix $output_dir/pb00.$subj.r05.tcat \
/Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/func.D003.r05.nii'[0..$]'
3dTcat -prefix $output_dir/pb00.$subj.r06.tcat \
/Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/func.D003.r06.nii'[0..$]'
# and make note of repetitions (TRs) per run
set tr_counts = ( 418 396 396 397 396 396 )
# -------------------------------------------------------
# enter the results directory (can begin processing data)
cd $output_dir
# -------------------------------------------------------
# extract initial volumes as automatic -blip_forward_dset
3dTcat -prefix blip_forward pb00.$subj.r01.tcat+orig"[0..29]"
# ---------------------------------------------------------
# QC: compute correlations with spherical ~averages
@radial_correlate -nfirst 0 -polort 3 -do_clean yes \
-rdir radcor.pb00.tcat \
pb00.$subj.r*.tcat+orig.HEAD
# ---------------------------------------------------------
# QC: look for columns of high variance
find_variance_lines.tcsh -polort 3 -nerode 2 \
-rdir vlines.pb00.tcat \
pb00.$subj.r*.tcat+orig.HEAD |& tee out.vlines.pb00.tcat.txt
# =============================== anat_unif ================================
# perform 'unifize' uniformity correction on anatomical dataset
3dUnifize -prefix T1.D003_unif -GM T1.D003+orig
# ========================== auto block: outcount ==========================
# QC: 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.05 of automask voxels are outliers
# - step() defines which TRs to remove via censoring
1deval -a outcount.r$run.1D -expr "1-step(a-0.05)" > 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
endif
end
# 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
# ---------------------------------------------------------
# check for potential 4095 saturation issues
foreach run ( $runs )
3dTto1D -method 4095_warn -input pb00.$subj.r$run.tcat+orig \
|& tee -a out.4095_all.txt
end
# make a file showing any resulting warnings
awk '/warning/ {print}' out.4095_all.txt | tee out.4095_warn.txt
# get run number and TR index for minimum outlier volume
set minindex = `3dTstat -argmin -prefix - outcount_rall.1D\'`
set ovals = ( `1d_tool.py -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
# ================================ despike =================================
# apply 3dDespike to each run
foreach run ( $runs )
3dDespike -NEW -nomask -prefix pb01.$subj.r$run.despike \
pb00.$subj.r$run.tcat+orig
end
# ================================= tshift =================================
# time shift data so all slice timing is the same
foreach run ( $runs )
3dTshift -tzero 0 -quintic -prefix pb02.$subj.r$run.tshift \
pb01.$subj.r$run.despike+orig
end
# ================================== blip ==================================
# compute blip up/down non-linear EPI distortion correction
# create median datasets from forward and reverse time series
3dTstat -median -prefix rm.blip.med.fwd blip_forward+orig
3dTstat -median -prefix rm.blip.med.rev blip_reverse+orig
# automask the median datasets
3dAutomask -apply_prefix rm.blip.med.masked.fwd rm.blip.med.fwd+orig
3dAutomask -apply_prefix rm.blip.med.masked.rev rm.blip.med.rev+orig
# compute the midpoint warp between the median datasets
3dQwarp -plusminus -pmNAMES Rev For \
-pblur 0.05 0.05 -blur -1 -1 \
-noweight -minpatch 9 \
-source rm.blip.med.masked.rev+orig \
-base rm.blip.med.masked.fwd+orig \
-prefix blip_warp
# warp median datasets (forward and each masked) for QC checks
# (and preserve obliquity)
3dNwarpApply -quintic -nwarp blip_warp_For_WARP+orig \
-source rm.blip.med.fwd+orig \
-prefix blip_med_for
3dNwarpApply -quintic -nwarp blip_warp_For_WARP+orig \
-source rm.blip.med.masked.fwd+orig \
-prefix blip_med_for_masked
3dNwarpApply -quintic -nwarp blip_warp_Rev_WARP+orig \
-source rm.blip.med.masked.rev+orig \
-prefix blip_med_rev_masked
3drefit -atrcopy blip_reverse+orig \
IJK_TO_DICOM_REAL \
blip_med_rev_masked+orig
# warp EPI time series data
foreach run ( $runs )
3dNwarpApply -quintic -nwarp blip_warp_For_WARP+orig \
-source pb02.$subj.r$run.tshift+orig \
-prefix pb03.$subj.r$run.blip
3drefit -atrcopy pb02.$subj.r$run.tshift+orig \
IJK_TO_DICOM_REAL \
pb03.$subj.r$run.blip+orig
end
# --------------------------------
# extract volreg registration base
3dbucket -prefix vr_base_min_outlier \
pb03.$subj.r$minoutrun.blip+orig"[$minouttr]"
# ================================= align ==================================
# for e2a: compute anat alignment transformation to EPI registration base
# (new anat will be intermediate, stripped, T1.D003_unif_ns+orig)
align_epi_anat.py -anat2epi -anat T1.D003_unif+orig \
-save_skullstrip -suffix _al_junk \
-epi vr_base_min_outlier+orig -epi_base 0 \
-epi_strip 3dAutomask \
-cost lpc+ZZ -giant_move -check_flip \
-volreg off -tshift off
# ================================== tlrc ==================================
# warp anatomy to standard space (non-linear warp)
auto_warp.py -base MNI152_2009_template_SSW.nii.gz -input \
T1.D003_unif_ns+orig \
-skull_strip_input no -unifize_input no
# move results up out of the awpy directory
# - NL-warped anat, affine warp, NL warp
# - use typical standard space name for anat
# - wildcard is a cheap way to go after any .gz
# - be sure NIFTI sform_code=2 means standard space
3dbucket -DAFNI_NIFTI_VIEW=tlrc \
-prefix T1.D003_unif_ns awpy/T1.D003_unif_ns.aw.nii*
mv awpy/anat.aff.Xat.1D .
mv awpy/anat.aff.qw_WARP.nii .
# ================================= volreg =================================
# align each dset to base volume, blip warp, to anat, warp to tlrc space
# (final warp input is same as blip input)
# verify that we have a +tlrc warp dataset
if ( ! -f T1.D003_unif_ns+tlrc.HEAD ) then
echo "** missing +tlrc warp dataset: T1.D003_unif_ns+tlrc.HEAD"
exit
endif
# register and warp
foreach run ( $runs )
# register each volume to the base image
3dvolreg -verbose -zpad 1 -base vr_base_min_outlier+orig \
-1Dfile dfile.r$run.1D -prefix rm.epi.volreg.r$run \
-cubic \
-1Dmatrix_save mat.r$run.vr.aff12.1D \
pb03.$subj.r$run.blip+orig
# create an all-1 dataset to mask the extents of the warp
3dcalc -overwrite -a pb03.$subj.r$run.blip+orig -expr 1 \
-prefix rm.epi.all1
# catenate blip/volreg/epi2anat/tlrc xforms
cat_matvec -ONELINE \
anat.aff.Xat.1D \
T1.D003_unif_al_junk_mat.aff12.1D -I \
mat.r$run.vr.aff12.1D > mat.r$run.warp.aff12.1D
# apply catenated xform: blip/volreg/epi2anat/tlrc/NLtlrc
# then apply non-linear standard-space warp
3dNwarpApply -master T1.D003_unif_ns+tlrc -dxyz 2.5 \
-source pb02.$subj.r$run.tshift+orig \
-nwarp "anat.aff.qw_WARP.nii mat.r$run.warp.aff12.1D \
blip_warp_For_WARP+orig" \
-prefix rm.epi.nomask.r$run
# warp the all-1 dataset for extents masking
3dNwarpApply -master T1.D003_unif_ns+tlrc -dxyz 2.5 \
-source rm.epi.all1+orig \
-nwarp "anat.aff.qw_WARP.nii mat.r$run.warp.aff12.1D" \
-interp cubic \
-ainterp 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
end
# 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 pb04.$subj.r$run.volreg
end
# warp the volreg base EPI dataset to make a final version
cat_matvec -ONELINE \
anat.aff.Xat.1D \
T1.D003_unif_al_junk_mat.aff12.1D -I > mat.basewarp.aff12.1D
3dNwarpApply -master T1.D003_unif_ns+tlrc -dxyz 2.5 \
-source vr_base_min_outlier+orig \
-nwarp "anat.aff.qw_WARP.nii mat.basewarp.aff12.1D" \
-prefix final_epi_vr_base_min_outlier
# create an anat_final dataset, aligned with stats
3dcopy T1.D003_unif_ns+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
# -----------------------------------------
# warp anat follower datasets (non-linear)
# warp follower dataset T1.D003_unif+orig
3dNwarpApply -source T1.D003_unif+orig \
-master anat_final.$subj+tlrc \
-ainterp wsinc5 -nwarp anat.aff.qw_WARP.nii anat.aff.Xat.1D\
-prefix anat_w_skull_warped
# ---------------------------------------------------------
# QC: compute correlations with spherical ~averages
@radial_correlate -nfirst 0 -polort 3 -do_clean yes \
-rdir radcor.pb04.volreg \
pb04.$subj.r*.volreg+tlrc.HEAD
# ================================== blur ==================================
# blur each volume of each run
foreach run ( $runs )
3dmerge -1blur_fwhm 4.0 -doall -prefix pb05.$subj.r$run.blur \
pb04.$subj.r$run.volreg+tlrc
end
# ================================== mask ==================================
# create 'full_mask' dataset (union mask)
foreach run ( $runs )
3dAutomask -prefix rm.mask_r$run pb05.$subj.r$run.blur+tlrc
end
# 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 T1.D003_unif_ns+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 tighter EPI mask by intersecting with anat mask
3dmask_tool -input full_mask.$subj+tlrc mask_anat.$subj+tlrc \
-inter -prefix mask_epi_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, MNI152_2009_template_SSW.nii.gz)
3dresample -master mask_epi_anat.$subj+tlrc -prefix ./rm.resam.group \
-input /Users/andy/abin/MNI152_2009_template_SSW.nii.gz'[0]'
# convert to binary group mask; fill gaps and holes
3dmask_tool -dilate_input 5 -5 -fill_holes -input rm.resam.group+tlrc \
-prefix mask_group
# note Dice coefficient of anat and template masks
3ddot -dodice mask_anat.$subj+tlrc mask_group+tlrc \
|& tee out.mask_at_dice.txt
# ================================= 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 pb05.$subj.r$run.blur+tlrc
3dcalc -a pb05.$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 pb06.$subj.r$run.scale
end
# ================================ regress =================================
# compute de-meaned motion parameters (for use in regression)
1d_tool.py -infile dfile_rall.1D -set_run_lengths $tr_counts \
-demean -write motion_demean.1D
# compute motion parameter derivatives (for use in regression)
1d_tool.py -infile dfile_rall.1D -set_run_lengths $tr_counts \
-derivative -demean -write motion_deriv.1D
# convert motion parameters for per-run regression
1d_tool.py -infile motion_demean.1D -set_run_lengths 418 396 396 397 396 396 \
-split_into_pad_runs mot_demean
1d_tool.py -infile motion_deriv.1D -set_run_lengths 418 396 396 397 396 396 \
-split_into_pad_runs mot_deriv
# create censor file motion_${subj}_censor.1D, for censoring motion
1d_tool.py -infile dfile_rall.1D -set_run_lengths $tr_counts \
-show_censor_count -censor_prev_TR \
-censor_motion 0.3 motion_${subj}
# combine multiple censor files
1deval -a motion_${subj}_censor.1D -b outcount_${subj}_censor.1D \
-expr "a*b" > censor_${subj}_combined_2.1D
# resample any -ROI_import dataset onto the EPI grid
# (and copy its labeltable)
3dresample -master final_epi_vr_base_min_outlier+tlrc -rmode NN \
-input ROI_import_MNI_2009c_asym+tlrc \
-prefix ROI_import_MNI_2009c_asym_resam
3drefit -copytables ROI_import_MNI_2009c_asym+tlrc \
ROI_import_MNI_2009c_asym_resam+tlrc
3drefit -cmap INT_CMAP ROI_import_MNI_2009c_asym_resam+tlrc
# note TRs that were not censored
# (apply from a text file, in case of a lot of censoring)
1d_tool.py -infile censor_${subj}_combined_2.1D \
-show_trs_uncensored space \
> out.keep_trs_rall.txt
set ktrs = "1dcat out.keep_trs_rall.txt"
# ------------------------------
# run the regression analysis
3dDeconvolve -input pb06.$subj.r*.scale+tlrc.HEAD \
-censor censor_${subj}_combined_2.1D \
-ortvec mot_demean.r01.1D mot_demean_r01 \
-ortvec mot_demean.r02.1D mot_demean_r02 \
-ortvec mot_demean.r03.1D mot_demean_r03 \
-ortvec mot_demean.r04.1D mot_demean_r04 \
-ortvec mot_demean.r05.1D mot_demean_r05 \
-ortvec mot_demean.r06.1D mot_demean_r06 \
-ortvec mot_deriv.r01.1D mot_deriv_r01 \
-ortvec mot_deriv.r02.1D mot_deriv_r02 \
-ortvec mot_deriv.r03.1D mot_deriv_r03 \
-ortvec mot_deriv.r04.1D mot_deriv_r04 \
-ortvec mot_deriv.r05.1D mot_deriv_r05 \
-ortvec mot_deriv.r06.1D mot_deriv_r06 \
-polort 3 -float \
-num_stimts 0 \
-fout -tout -x1D X.xmat.1D -xjpeg X.jpg \
-x1D_uncensored X.nocensor.xmat.1D \
-fitts fitts.$subj \
-errts errts.${subj} \
-x1D_stop \
-bucket stats.$subj
# -- use 3dTproject to project out regression matrix --
# (make errts like 3dDeconvolve, but more quickly)
3dTproject -polort 0 -input pb06.$subj.r*.scale+tlrc.HEAD \
-censor censor_${subj}_combined_2.1D -cenmode ZERO \
-ort X.nocensor.xmat.1D -prefix errts.${subj}.tproject
# if 3dDeconvolve fails, terminate the script
if ( $status != 0 ) then
echo '---------------------------------------'
echo '** 3dDeconvolve error, failing...'
echo ' (consider the file 3dDeconvolve.err)'
exit
endif
# display any large pairwise correlations from the X-matrix
1d_tool.py -show_cormat_warnings -infile X.xmat.1D |& tee out.cormat_warn.txt
# display degrees of freedom info from X-matrix
1d_tool.py -show_df_info -infile X.xmat.1D |& tee out.df_info.txt
# create an all_runs dataset to match the fitts, errts, etc.
3dTcat -prefix all_runs.$subj pb06.$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}.tproject+tlrc"[$ktrs]"
3dcalc -a rm.signal.all+tlrc \
-b rm.noise.all+tlrc \
-expr 'a/b' -prefix TSNR.$subj
# --------------------------------------------------
# compute TSNR stats for dset labels: brain, MNI_2009c_asym
compute_ROI_stats.tcsh \
-out_dir tsnr_stats_regress \
-stats_file tsnr_stats_regress/stats_auto_brain.txt \
-dset_ROI mask_epi_anat.$subj+tlrc \
-dset_data TSNR.$subj+tlrc \
-rset_label brain \
-rval_list 1
compute_ROI_stats.tcsh \
-out_dir tsnr_stats_regress \
-stats_file tsnr_stats_regress/stats_auto_MNI_2009c_asym.txt \
-dset_ROI ROI_import_MNI_2009c_asym_resam+tlrc \
-dset_data TSNR.$subj+tlrc \
-rset_label MNI_2009c_asym \
-rval_list ALL_LT
# ---------------------------------------------------
# compute and store GCOR (global correlation average)
# (sum of squares of global mean of unit errts)
3dTnorm -norm2 -prefix rm.errts.unit errts.${subj}.tproject+tlrc
3dmaskave -quiet -mask mask_epi_anat.$subj+tlrc \
rm.errts.unit+tlrc > mean.errts.unit.1D
3dTstat -sos -prefix - mean.errts.unit.1D\' > out.gcor.1D
echo "-- GCOR = `cat out.gcor.1D`"
# ---------------------------------------------------
# compute correlation volume
# (per voxel: correlation with masked brain average)
3dmaskave -quiet -mask mask_epi_anat.$subj+tlrc \
errts.${subj}.tproject+tlrc > mean.errts.1D
3dTcorr1D -prefix corr_brain errts.${subj}.tproject+tlrc mean.errts.1D
# --------------------------------------------------
# compute sum of baseline (all) regressors
3dTstat -sum -prefix sum_baseline.1D X.nocensor.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 = `1d_tool.py -infile X.xmat.1D -show_trs_uncensored encoded \
-show_trs_run $run`
if ( $trs == "" ) continue
3dFWHMx -detrend -mask mask_epi_anat.$subj+tlrc \
-ACF files_ACF/out.3dFWHMx.ACF.epits.r$run.1D \
all_runs.$subj+tlrc"[$trs]" >> blur.epits.1D
end
# 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 = `1d_tool.py -infile X.xmat.1D -show_trs_uncensored encoded \
-show_trs_run $run`
if ( $trs == "" ) continue
3dFWHMx -detrend -mask mask_epi_anat.$subj+tlrc \
-ACF files_ACF/out.3dFWHMx.ACF.errts.r$run.1D \
errts.${subj}.tproject+tlrc"[$trs]" >> blur.errts.1D
end
# 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
# ========================= auto block: QC_review ==========================
# generate quality control review scripts and HTML report
# generate a review script for the unprocessed EPI data
gen_epi_review.py -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)
# write AP uvars into a simple txt file
cat << EOF > out.ap_uvars.txt
mot_limit : 0.3
out_limit : 0.05
copy_anat : T1.D003+orig.HEAD
errts_dset : errts.${subj}.tproject+tlrc.HEAD
mask_dset : mask_epi_anat.$subj+tlrc.HEAD
template : MNI152_2009_template_SSW.nii.gz
ss_review_dset : out.ss_review.$subj.txt
max_4095_warn_dset : out.4095_warn.txt
vlines_tcat_dir : vlines.pb00.tcat
EOF
# and convert the txt format to JSON
cat out.ap_uvars.txt | afni_python_wrapper.py -eval "data_file_to_json()" \
> out.ap_uvars.json
# initialize gen_ss_review_scripts.py with out.ap_uvars.json
gen_ss_review_scripts.py -exit0 \
-init_uvars_json out.ap_uvars.json \
-write_uvars_json out.ss_review_uvars.json
# ========================== auto block: finalize ==========================
# remove temporary files
\rm -fr rm.* awpy
# --------------------------------------------------
# if the basic subject review script is here, run it
# (want this to be the last text output)
if ( -e @ss_review_basic ) then
./@ss_review_basic |& tee out.ss_review.$subj.txt
# generate html ss review pages
# (akin to static images from running @ss_review_driver)
apqc_make_tcsh.py -review_style pythonic -subj_dir . \
-uvar_json out.ss_review_uvars.json
apqc_make_html.py -qc_dir QC_$subj
echo "\nconsider running: \n"
echo " afni_open -b /Users/andy/Desktop/Research/LDT/fMRI_data/preproc_afni/D003/QC_$subj/index.html"
echo ""
endif
# return to parent directory (just in case...)
cd ..
echo "execution finished: `date`"
# ==========================================================================
# script generated by the command:
#
# afni_proc.py -subj_id D003 -script proc_D003 -out_dir \
# /Users/andy/Desktop/Research/LDT/fMRI_data/preproc_afni/D003 -blocks \
# despike tshift align tlrc volreg blur mask scale regress -copy_anat \
# /Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/T1.D003.nii \
# -anat_has_skull yes -anat_uniform_method unifize -anat_unif_GM yes \
# -dsets \
# /Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/func.D003.r01.nii \
# /Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/func.D003.r02.nii \
# /Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/func.D003.r03.nii \
# /Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/func.D003.r04.nii \
# /Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/func.D003.r05.nii \
# /Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/func.D003.r06.nii \
# -radial_correlate_blocks tcat volreg -blip_reverse_dset \
# /Users/andy/Desktop/Research/LDT/fMRI_data/raw/D003/dist_PA.D003.nii \
# -tlrc_base MNI152_2009_template_SSW.nii.gz -tlrc_NL_warp \
# -align_opts_aea -cost lpc+ZZ -giant_move -check_flip -volreg_align_to \
# MIN_OUTLIER -volreg_align_e2a -volreg_tlrc_warp -blur_size 4.0 \
# -mask_epi_anat yes -regress_motion_per_run -regress_censor_motion 0.3 \
# -regress_censor_outliers 0.05 -regress_apply_mot_types demean deriv \
# -regress_est_blur_epits -regress_est_blur_errts -html_review_style \
# pythonic -execute