Non-linear warp of epi to experimental T1 for multiple epi runs + concatenating affine + non-linear

AFNI version info (afni -ver): Version AFNI_24.1.02 'Publius Septimius Geta'

Hello!

I have 8 epi runs and a T1 anatomical which I am pre-processing using afni_proc.py (blocks tcat despike align volreg surf blur scale; code below).

I am aligning the epi to the anatomical. The alignment is quite good, but a little bit off around somatosensory/ motor areas, where we are investigating.

I want to apply non-linear warping to improve the warp of epi to anatomical, and I have achieved some improvement in alignment of the vr_base_min_outlier (from afni_proc.py) to our affine aligned anat (from afni_proc.py) with 3dQWarp (as suggested in this post Cannot get good alignment with align_epi_anat.py or 3dAllineate - #6 by ptaylor; code below).

This obviously only warps the single sub-brick in the vr_base_min_outlier. I want non-linear warping for ALL runs of my epi data. Questions:

  1. Do I calculate the non-linear warp once (on one epi sub-brick, as I have done), and then apply this warp to the rest of the epi data? Or should I be calculating the warp separately for each run? Or something more complicated which calculates the warp over time (all sub-bricks)?
  2. Can I/ should I concatenate this warp in with my volreg/post_vr_allin/epi2anat xforms so they only get applied once? If so, how? I had a look at how a non-linear warp from a tlrc block is concatenated with the affine alignment in afni_proc.py, but the output files are different as it uses @SSWarper (and I am using 3dQwarp) so I was unsure how to proceed.
  3. Should I be using @SSWarper for my non-linear alignment of epi to anat instead? And then copying the concatenation process in volreg?
  4. (maybe should have been question 1) Can you confirm I cannot use afni_proc.py to do this? I believe this is the case? If you can, can you let me know how to implement? I did some reading but got the impression it was not possible.

Thank you so much for your help,
I really didn’t want to get this wrong - it’s quite complicated!

Harriet

3dQwarp code

set anat = ${subj}_acq-UNIDEN_run-1_T1w_SS_al_junk
set epi = vr_base_min_outlier
3dresample \
                    -master ${epi}+orig \
                    -input ${anat}+orig \
                    -prefix ${anat}_resamp
3dQwarp \
            -source ${anat}${resamp_suffix}+orig.HEAD \
            -base ${epi}+orig \
            -prefix ${anat}_3dQ \
            -lpc -maxlev 0 -verb \
            -iwarp

3dNwarpApply \
         -nwarp ${anat}_3dQ_WARPINV+orig  \
         -source ${epi}+orig \
         -prefix ${epi}_3dQ_applied

Volreg code from proc script

# catenate volreg/post_vr_allin/epi2anat xforms
    cat_matvec -ONELINE                                                       \
               ${subj}_acq-UNIDEN_run-1_T1w_SS_al_junk_mat.aff12.1D -I \
               mat.vr_xrun_allin.r$run.aff12.1D                               \
               mat.r$run.vr.aff12.1D > mat.r$run.warp.aff12.1D

    # apply catenated xform: volreg/post_vr_allin/epi2anat
    3dAllineate -base ${subj}_acq-UNIDEN_run-1_T1w_SS+orig             \
                -input pb01.$subj.r$run.despike+orig                          \
                -1Dmatrix_apply mat.r$run.warp.aff12.1D                       \
                -mast_dxyz 0.75 -final wsinc5                                 \
                -prefix rm.epi.nomask.r$run

Needs to become more like this? From this post MNI normalization question with 3dNwarpApply - #6 by philippn

# catenate volreg/epi2anat/tlrc xforms
    cat_matvec -ONELINE                                                  \
               anatQQ.FT.aff12.1D                                        \
               anatSS.FT_al_junk_mat.aff12.1D -I                         \
               mat.r$run.vr.aff12.1D > mat.r$run.warp.aff12.1D

    # apply catenated xform: volreg/epi2anat/tlrc/NLtlrc
    # then apply non-linear standard-space warp
    3dNwarpApply -master anatQQ.FT+tlrc -dxyz 2.5                        \
                 -source pb01.$subj.r$run.tshift+orig                    \
                 -nwarp "anatQQ.FT_WARP.nii mat.r$run.warp.aff12.1D"     \
                 -prefix rm.epi.nomask.r$run

Other parts from my proc script that also need to be modified as they also pertain to combining/ applying warps????

    # warp the all-1 dataset for extents masking 
    3dAllineate -base ${subj}_acq-UNIDEN_run-1_T1w_SS+orig             \
                -input rm.epi.all1+orig                                       \
                -1Dmatrix_apply mat.r$run.warp.aff12.1D                       \
                -mast_dxyz 0.75 -final NN -quiet                              \
                -prefix rm.epi.1.r$run

Continued

# warp the volreg base EPI dataset to make a final version
cat_matvec -ONELINE                                                           \
    ${subj}_acq-UNIDEN_run-1_T1w_SS_al_junk_mat.aff12.1D -I  >         \
    mat.basewarp.aff12.1D

3dAllineate -base ${subj}_acq-UNIDEN_run-1_T1w_SS+orig                 \
            -input vr_base_min_outlier+orig                                   \
            -1Dmatrix_apply mat.basewarp.aff12.1D                             \
            -mast_dxyz 0.75 -final wsinc5                                     \
            -prefix final_epi_vr_base_min_outlier

Full afni_proc.py script and proc script if needed:

#!/bin/tcsh

# Chosen script for MND pre proc
# Edited to try for block and phase at once

set sub = sub-005
set ses = ses-01

set base_dir = /90days/uqhdemp1/Data/MND/bids

set extension = .nii.gz
set anat_raw_dir = $base_dir/${sub}/${ses}/anat
set anat_raw = ${sub}_${ses}_acq-UNIDEN_run-1_T1w
set anat_ss_dir = $base_dir/derivatives/hdbet/${sub}/${ses}/anat
set anat_ss = ${anat_raw}_SS

set dset_dir = $base_dir/${sub}/${ses}/func
set dset_base = ${sub}_${ses}_task

set surf_dir = $base_dir/derivatives/FastSurfer/output/${sub}_${ses}/SUMA
set surf_anat = ${sub}_${ses}_SurfVol.nii
set surf_spec = ${sub}_${ses}_lh.spec

# Using .8x3; Could also try 3.2 (x4); think it used to error when trying to use 2mm ( as blur needs to be bigger relative to voxel size, I think)
set blur_size = 2.4


afni_proc.py \
	-subj_id			${sub}_${ses}							\
	-dsets				$dset_dir/${dset_base}*${extension}				\
	-blocks				tcat despike align volreg surf blur scale			\
	-radial_correlate_blocks	tcat volreg							\
	-copy_anat			$anat_ss_dir/${anat_ss}${extension}				\
	-anat_has_skull			no								\
	-anat_follower			orig_anat_w_skull anat $anat_raw_dir/${anat_raw}${extension}	\
	-surf_anat			$surf_dir/$surf_anat						\
	-surf_spec			$surf_dir/$surf_spec						\
	-tcat_remove_first_trs    	0								\
	-align_opts_aea				-giant_move						\
						-partial_coverage					\
						-cost lpc+ZZ						\
	-align_unifize_epi		local								\
	-volreg_align_e2a										\
	-volreg_align_to		MIN_OUTLIER							\
	-volreg_post_vr_allin 		yes								\
	-volreg_pvra_base_index 	MIN_OUTLIER							\
	-volreg_interp 			-Fourier							\
	-volreg_warp_final_interp 	wsinc5								\
	-volreg_compute_tsnr		yes								\
	-blur_size			$blur_size							\
	-html_review_style		basic
## Start sub loop
foreach sub ( $subs )

    # Create subj var
    set subj = "${sub}_${ses}"

    ## Determine number of phase block runs
    # Edit to do if statement based on sub, and check if sub is in a text file list
    # set list of runs and make note of repetitions (TRs) per run
    # if () then
        set runs = (`count -digits 2 1 8`)
        set tr_counts = ( 230 230 230 230 195 195 195 195 )
    # else
    #    set runs = (`count -digits 2 1 7`)
    #    set tr_counts = ( 230 230 230 230 195 195 195 )
    # endif
    set last_run = $runs[$#runs]
    
    # Set blurs
    set blurs = ("0" "2.4")

    # assign output directory name
    set output_dir = $subj.results

    # verify that the results directory does not yet exist
    if ( -d $output_dir ) then
        echo output dir "$subj.results" already exists
        exit
    else
        # create results and stimuli directories
        mkdir -p $output_dir
        mkdir $output_dir/stimuli
    endif


    # =========================== auto block: setup ============================
    # script setup

    echo "auto-generated by afni_proc.py, Thu Sep 19 12:21:35 2024"
    echo "(version 7.74, April 8, 2024)"
    echo "execution started: `date`"

    # 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


    # copy anatomy to results dir
    3dcopy                                                                                                                 \
        /90days/uqhdemp1/Data/MND/bids/derivatives/hdbet/${sub}/${ses}/anat/${subj}_acq-UNIDEN_run-1_T1w_SS.nii.gz \
        $output_dir/${subj}_acq-UNIDEN_run-1_T1w_SS

    # copy anatomical follower datasets into the results dir
    3dcopy                                                                                            \
        /90days/uqhdemp1/Data/MND/bids/${sub}/${ses}/anat/${subj}_acq-UNIDEN_run-1_T1w.nii.gz \
        $output_dir/copy_af_orig_anat_w_skull

    # ============================ 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      \
        /90days/uqhdemp1/Data/MND/bids/${sub}/${ses}/func/${subj}_task-BLOCK1_bold.nii.gz'[0..$]'
    3dTcat -prefix $output_dir/pb00.$subj.r02.tcat      \
        /90days/uqhdemp1/Data/MND/bids/${sub}/${ses}/func/${subj}_task-BLOCK2_bold.nii.gz'[0..$]'
    3dTcat -prefix $output_dir/pb00.$subj.r03.tcat      \
        /90days/uqhdemp1/Data/MND/bids/${sub}/${ses}/func/${subj}_task-BLOCK3_bold.nii.gz'[0..$]'
    3dTcat -prefix $output_dir/pb00.$subj.r04.tcat      \
        /90days/uqhdemp1/Data/MND/bids/${sub}/${ses}/func/${subj}_task-BLOCK4_bold.nii.gz'[0..$]'
    3dTcat -prefix $output_dir/pb00.$subj.r05.tcat      \
        /90days/uqhdemp1/Data/MND/bids/${sub}/${ses}/func/${subj}_task-PHASE1_bold.nii.gz'[0..$]'
    3dTcat -prefix $output_dir/pb00.$subj.r06.tcat      \
        /90days/uqhdemp1/Data/MND/bids/${sub}/${ses}/func/${subj}_task-PHASE2_bold.nii.gz'[0..$]'
    3dTcat -prefix $output_dir/pb00.$subj.r07.tcat      \
        /90days/uqhdemp1/Data/MND/bids/${sub}/${ses}/func/${subj}_task-PHASE3_bold.nii.gz'[0..$]'

    if ($last_run == 08) then
        3dTcat -prefix $output_dir/pb00.$subj.r08.tcat      \
            /90days/uqhdemp1/Data/MND/bids/${sub}/${ses}/func/${subj}_task-PHASE4_bold.nii.gz'[0..$]'
    endif

    # enter the results directory (can begin processing data)
    cd $output_dir

    # ---------------------------------------------------------
    # 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

    # ========================== 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

        # 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

    # ---------------------------------------------------------
    # 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

    # --------------------------------
    # extract volreg registration base
    3dbucket -prefix vr_base_min_outlier                        \
        pb01.$subj.r$minoutrun.despike+orig"[$minouttr]"

    # ================================= align ==================================
    # for e2a: compute anat alignment transformation to EPI registration base
    # (new anat will be current ${subj}_acq-UNIDEN_run-1_T1w_SS+orig)
    # run (localized) uniformity correction on EPI base
    3dLocalUnifize -input vr_base_min_outlier+orig -prefix                        \
        vr_base_min_outlier_unif

    align_epi_anat.py -anat2epi -anat ${subj}_acq-UNIDEN_run-1_T1w_SS+orig \
        -suffix _al_junk                                                       \
        -epi vr_base_min_outlier_unif+orig -epi_base 0                         \
        -epi_strip 3dAutomask                                                  \
        -anat_has_skull no                                                     \
        -giant_move -partial_coverage -cost lpc+ZZ                             \
        -volreg off -tshift off

    # ================================= volreg =================================
    # align each dset to base volume, across runs, to anat

    # register and warp
    foreach run ( $runs )
        # extract MIN_OUTLIER index for current run
        set min_outlier_index = `3dTstat -argmin -prefix - outcount.r$run.1D\'`
        
        # extract volreg base for this run
        3dbucket -prefix vr_base_per_run_r$run                                    \
            pb01.$subj.r$run.despike+orig"[$min_outlier_index]"
        
        # and compute xforms for cross-run allin to vr_base
        3dAllineate -base vr_base_min_outlier+orig                                \
                    -source vr_base_per_run_r$run+orig                            \
                    -1Dfile vr_xrun_allin_dfile.m12.r$run.1D                      \
                    -1Dmatrix_save mat.vr_xrun_allin.r$run.aff12.1D               \
                    -autoweight -source_automask                                  \
                    -lpa -cubic

        # register each volume to the base image
        3dvolreg -verbose -zpad 1 -base vr_base_per_run_r$run+orig                \
                -1Dfile dfile.r$run.1D -prefix rm.epi.volreg.r$run               \
                -Fourier                                                         \
                -1Dmatrix_save mat.r$run.vr.aff12.1D                             \
                pb01.$subj.r$run.despike+orig

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

        # catenate volreg/post_vr_allin/epi2anat xforms
        cat_matvec -ONELINE                                                       \
                ${subj}_acq-UNIDEN_run-1_T1w_SS_al_junk_mat.aff12.1D -I \
                mat.vr_xrun_allin.r$run.aff12.1D                               \
                mat.r$run.vr.aff12.1D > mat.r$run.warp.aff12.1D

        # apply catenated xform: volreg/post_vr_allin/epi2anat
        3dAllineate -base ${subj}_acq-UNIDEN_run-1_T1w_SS+orig             \
                    -input pb01.$subj.r$run.despike+orig                          \
                    -1Dmatrix_apply mat.r$run.warp.aff12.1D                       \
                    -mast_dxyz 0.75 -final wsinc5                                 \
                    -prefix rm.epi.nomask.r$run

        # warp the all-1 dataset for extents masking 
        3dAllineate -base ${subj}_acq-UNIDEN_run-1_T1w_SS+orig             \
                    -input rm.epi.all1+orig                                       \
                    -1Dmatrix_apply mat.r$run.warp.aff12.1D                       \
                    -mast_dxyz 0.75 -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+orig
    end

    # make a single file of registration params
    cat dfile.r*.1D > dfile_rall.1D

    # compute motion magnitude time series: the Euclidean norm
    # (sqrt(sum squares)) of the motion parameter derivatives
    1d_tool.py -infile dfile_rall.1D                                              \
            -set_run_lengths $tr_counts                   \
            -derivative -collapse_cols euclidean_norm                          \
            -write motion_${subj}_enorm.1D

    # ----------------------------------------
    # create the extents mask: mask_epi_extents+orig
    # (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+orig -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+orig -b mask_epi_extents+orig               \
            -expr 'a*b' -prefix pb02.$subj.r$run.volreg
    end

    # warp the volreg base EPI dataset to make a final version
    cat_matvec -ONELINE                                                           \
        ${subj}_acq-UNIDEN_run-1_T1w_SS_al_junk_mat.aff12.1D -I  >         \
        mat.basewarp.aff12.1D

    3dAllineate -base ${subj}_acq-UNIDEN_run-1_T1w_SS+orig                 \
                -input vr_base_min_outlier+orig                                   \
                -1Dmatrix_apply mat.basewarp.aff12.1D                             \
                -mast_dxyz 0.75 -final wsinc5                                     \
                -prefix final_epi_vr_base_min_outlier

    # create an anat_final dataset, aligned with stats
    3dcopy ${subj}_acq-UNIDEN_run-1_T1w_SS+orig anat_final.$subj

    # record final registration costs
    3dAllineate -base final_epi_vr_base_min_outlier+orig -allcostX                \
                -input anat_final.$subj+orig |& tee out.allcostX.txt

    # --------------------------------------
    # create a TSNR dataset, just from run 1
    3dTstat -mean -prefix rm.signal.vreg.r01 pb02.$subj.r01.volreg+orig
    3dDetrend -polort 3 -prefix rm.noise.det -overwrite pb02.$subj.r01.volreg+orig
    3dTstat -stdev -prefix rm.noise.vreg.r01 rm.noise.det+orig
    3dcalc -a rm.signal.vreg.r01+orig                                             \
        -b rm.noise.vreg.r01+orig                                              \
        -c mask_epi_extents+orig                                               \
        -expr 'c*a/b' -prefix TSNR.vreg.r01.$subj

    # ---------------------------------------------------------
    # QC: compute correlations with spherical ~averages
    @radial_correlate -nfirst 0 -polort 3 -do_clean yes                           \
                    -rdir radcor.pb02.volreg                                    \
                    pb02.$subj.r*.volreg+orig.HEAD

    # ======================= surf (map data to surface) =======================
    # map EPI data to the surface domain

    # set directory variables
    set surface_dir =                                                          \
        /90days/uqhdemp1/Data/MND/bids/derivatives/FastSurfer/output/${subj}/SUMA

    # align the surface anatomy with the current experiment anatomy
    @SUMA_AlignToExperiment -exp_anat anat_final.$subj+orig                    \
                            -surf_anat $surface_dir/${subj}_SurfVol.nii \
                            -wd -strip_skull surf_anat                         \
                            -atlas_followers -overwrite_resp S                 \
                            -prefix ${subj}_SurfVol_Alnd_Exp 

    # map volume data to the surface of each hemisphere
    foreach hemi ( lh )
        foreach run ( $runs )
            3dVol2Surf -spec $surface_dir/${subj}_${hemi}.spec          \
                    -sv ${subj}_SurfVol_Alnd_Exp+orig                       \
                    -surf_A smoothwm                                        \
                    -surf_B pial                                            \
                    -f_index nodes                                          \
                    -f_steps 10                                             \
                    -map_func ave                                           \
                    -oob_value 0                                            \
                    -grid_parent pb02.$subj.r$run.volreg+orig               \
                    -out_niml pb03.$subj.$hemi.r$run.surf.niml.dset 
        end
    end

    # make local script for running suma, and make it executable
    echo suma -spec $surface_dir/${subj}_lh.spec                        \
            -sv ${subj}_SurfVol_Alnd_Exp+orig > run_suma
    chmod 755 run_suma

    # =========================== blur (on surface) ============================
    foreach hemi ( lh )
        foreach run ( $runs )
            foreach blur ( $blurs )
                if ( $blur == 0 ) then
                    echo $blur
                    echo "just copying and renaming"
                    # Copy file with next proc block name and blur 0
                    3dcopy \
                        pb03.$subj.$hemi.r$run.surf.niml.dset \
                        pb04.$subj.$hemi.r$run.blur${blur}.niml.dset
                else
                    echo $blur
                    # to save time, estimate blur parameters only once
                    if ( ! -f surf.smooth.params.1D ) then
                        SurfSmooth -spec $surface_dir/${subj}_${hemi}.spec    \
                                -surf_A smoothwm                                  \
                                -input pb03.$subj.$hemi.r$run.surf.niml.dset      \
                                -met HEAT_07                                      \
                                -target_fwhm $blur                                  \
                                -blurmaster pb03.$subj.$hemi.r$run.surf.niml.dset \
                                -detrend_master                                   \
                                -output pb04.$subj.$hemi.r$run.blur${blur}.niml.dset     \
                                | tee surf.smooth.params.1D
                        set bfile  = pb04.$subj.$hemi.r$run.blur${blur}.niml.dset.1D.smrec
                        set params = ( `1dcat $bfile | tail -n 1` )
                    else
                        # get blur params from smrec file
                        SurfSmooth -spec $surface_dir/${subj}_${hemi}.spec    \
                                -surf_A smoothwm                                  \
                                -input pb03.$subj.$hemi.r$run.surf.niml.dset      \
                                -met HEAT_07                                      \
                                -Niter $params[1]                                 \
                                -sigma $params[3]                                 \
                                -output pb04.$subj.$hemi.r$run.blur${blur}.niml.dset
                    endif
                endif
            end
        end
    end

    # ================================= 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 hemi ( lh )
        foreach run ( $runs )
            foreach blur ( $blurs )
                3dTstat -prefix rm.$hemi.mean_r$run.blur${blur}.niml.dset    \
                    pb04.$subj.$hemi.r$run.blur${blur}.niml.dset
                3dcalc -a pb04.$subj.$hemi.r$run.blur${blur}.niml.dset  \
                    -b rm.$hemi.mean_r$run.blur${blur}.niml.dset          \
                    -expr 'min(200, a/b*100)*step(a)*step(b)' \
                    -prefix pb05.$subj.$hemi.r$run.blur${blur}.scale.niml.dset
            end
        end
    end

    # ========================= 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

    # ========================== 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 ) 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 basic -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 $subj.results/QC_$subj/index.html"
        echo ""
    endif

    # return to parent directory (just in case...)
    cd ..

    echo "execution finished: `date`"

end
# end sub loop

# ==========================================================================
# script generated by the command:
#
# afni_proc.py -subj_id sub-005_ses-01 -dsets                                                                            \
#     /90days/uqhdemp1/Data/MND/bids/sub-005/ses-01/func/sub-005_ses-01_task-BLOCK1_bold.nii.gz                          \
#     /90days/uqhdemp1/Data/MND/bids/sub-005/ses-01/func/sub-005_ses-01_task-BLOCK2_bold.nii.gz                          \
#     /90days/uqhdemp1/Data/MND/bids/sub-005/ses-01/func/sub-005_ses-01_task-BLOCK3_bold.nii.gz                          \
#     /90days/uqhdemp1/Data/MND/bids/sub-005/ses-01/func/sub-005_ses-01_task-BLOCK4_bold.nii.gz                          \
#     /90days/uqhdemp1/Data/MND/bids/sub-005/ses-01/func/sub-005_ses-01_task-PHASE1_bold.nii.gz                          \
#     /90days/uqhdemp1/Data/MND/bids/sub-005/ses-01/func/sub-005_ses-01_task-PHASE2_bold.nii.gz                          \
#     /90days/uqhdemp1/Data/MND/bids/sub-005/ses-01/func/sub-005_ses-01_task-PHASE3_bold.nii.gz                          \
#     /90days/uqhdemp1/Data/MND/bids/sub-005/ses-01/func/sub-005_ses-01_task-PHASE4_bold.nii.gz                          \
#     -blocks tcat despike align volreg surf blur scale                                                                  \
#     -radial_correlate_blocks tcat volreg -copy_anat                                                                    \
#     /90days/uqhdemp1/Data/MND/bids/derivatives/hdbet/sub-005/ses-01/anat/sub-005_ses-01_acq-UNIDEN_run-1_T1w_SS.nii.gz \
#     -anat_has_skull no -anat_follower orig_anat_w_skull anat                                                           \
#     /90days/uqhdemp1/Data/MND/bids/sub-005/ses-01/anat/sub-005_ses-01_acq-UNIDEN_run-1_T1w.nii.gz                      \
#     -surf_anat                                                                                                         \
#     /90days/uqhdemp1/Data/MND/bids/derivatives/FastSurfer/output/sub-005_ses-01/SUMA/sub-005_ses-01_SurfVol.nii        \
#     -surf_spec                                                                                                         \
#     /90days/uqhdemp1/Data/MND/bids/derivatives/FastSurfer/output/sub-005_ses-01/SUMA/sub-005_ses-01_lh.spec            \
#     -tcat_remove_first_trs 0 -align_opts_aea -giant_move -partial_coverage                                             \
#     -cost lpc+ZZ -align_unifize_epi local -volreg_align_e2a                                                            \
#     -volreg_align_to MIN_OUTLIER -volreg_post_vr_allin yes                                                             \
#     -volreg_pvra_base_index MIN_OUTLIER -volreg_interp -Fourier                                                        \
#     -volreg_warp_final_interp wsinc5 -volreg_compute_tsnr yes -blur_size                                               \
#     2.4 -html_review_style basic

Howdy-

Just taking a step back to see if I understand. Is the issue that you have an extra, pre-calculated nonlinear warp to applied directly to the EPI that tries to undo some B0 inhomogeneity distortion, and you want to include that in your processing stream?

If that is the case, you can include that directly within your afni_proc.py command, as a somewhat recent addition:

-blip_warp_dset DSET    : specify extra options for 3dQwarp

        e.g. -blip_warp_dset epi_b0_WARP.nii.gz

    This option allows the user to pass a pre-computed distortion warp
    dataset, to replace the computation of a warp in the blip block.
    The most likely use is to first run epi_b0_correct.py for a b0
    distortion map computation, rather than the reverse phase encoding
    method that would be computed with afni_proc.py.

    When applying this option in afni_proc.py, instead of using options
    like:

        -blip_forward_dset DSET_FORWARD \
        -blip_reverse_dset DSET_REVERSE \
        -blip_opts_qw      OPTIONS ...  \

    use just this one option to pass the warp:

        -blip_warp_dset epi_b0_WARP.nii.gz \

    Please see 'epi_b0_correct.py -help' for more information.

Note: I think you will need to update your version of AFNI to have this available, because it went into the codebase on Aug 5, 2024. The next version built after that was 24.2.02. Just updating to current available AFNI will get that for you, too, in likely the simplest manner.

But if I am misunderstanding, please let me know.

-pt

1 Like

Hi pt,

Thank you very much for your quick response.

That is just about right – we want to improve the alignment of our EPI to anat (it is good, but could be a bit better in some key spots). We were trying to improve this alignment with 3dQwarp.

As you say, that misalignment is probably caused by b0 inhomogeneity distortion (we have no b0 maps or reverse phase scans to correct in other ways), so, I am keen to try the blip block in afni_proc.py with a pre-computed warp, as you suggest.

Thus, I should perform the following steps?

  1. Align EPI to anat (e.g., using align_epi_anat.py)
  2. Non-linear alignment of EPI to anat (using 3dQwarp)
  3. Retain the WARP file from 3dQwarp (mine is named: sub-005_ses-01_acq-UNIDEN_run-1_T1w_SS_al_junk_3dQ_WARPINV+orig), discard all others as not needed.
  4. Perform our actual pre-processing with afni_proc.py, using these options for blip block:
    -blip_warp_dset sub-005_ses-01_acq-UNIDEN_run-1_T1w_SS_al_junk_3dQ_WARPINV+orig

Just to make sure I didn’t misunderstand, the afni_proc.py blip_warp_dset option instructions say this: “The most likely use is to first run epi_b0_correct.py for a b0 distortion map computation”.

Did you mean that instead of steps 1-3 above, I should run epi_b0_correct.py to calculate the non-linear warp? (i.e., instead of 3dQwarp in step 2 above) I think I don’t have the right files to run this (e.g., frequency dset).

Thank you!

Harriet

Hi again,
I have been preparing my data to try the blip approach you suggested (I have not yet been able to run it as we are having some issues with the new afni version). However, having done the data preparation, I am wondering if the blip approach IS actually sensible for our needs.

As a reminder, I am using non-linear warping to improve the epi2anat alignment. This non-linear warp is not actually pre-calculated (before pre-processing) – which would have made it very easy to put in when starting our afni_proc.py script to pre-process from the start.

Rather, I have been calculating the non-linear warp after our tcat, despike, align and volreg blocks. This is because in the volreg block all our transformations (volreg/post_vr_allin/epi2anat xforms) are concatenated to give us the final_epi_vr_base_min_outlier file. I then use 3dQwarp to calc the warp between the final_epi_vr_base_min_outlier and our anatomical.

Thus, if I wanted to enter a -blip_warp_dset when I first run afni_proc.py, I first need to run most of our pre-processing to get the non-linear warp, then start again and feed that warp to afni_proc.py.
Question 1: Is that correct, or is there a way to do this without doing most of our pre-processing twice over?

I have been wondering again if I shouldn’t just run the pre-processing steps above as described (tcat, despike, align, volreg), use 3dQwarp to get the final_epi_vr_base_min_outlier to anat warp, then use 3dNwarpApply to apply this warp to the epi data for our 4 runs.

Question 2: But, should I be calculating non-linear warp for each run separately then applying it to just that run (rather than calculating the warp on the min_outlier and apply to all runs)?
Question 3: Is ok that I am calculating the non-linear warp on 1 sub-brick (final_epi_vr_base_min_outlier) but applying the warp to +time data (whole runs of epi)?
Question 4: Should I be trying to concatenate this non-linear warp with my other transformations (in the volreg step - volreg/post_vr_allin/epi2anat xforms)?

Thank you once again for your help!

H

Code

# Set warp direction: 1 = anat to epi; 2 = epi to anat
set warpdir = 2

	# set sub = sub-005
	set ses = ses-01
	set subj = ${sub}_${ses}

	## 3dQwarp page suggests we may want to always have epi as base - so partial coverage is less of an issue (see 3dQwarp help page)
	set anat = anat_final.$subj${exp_suf}
	set epi = final_epi_vr_base_min_outlier${exp_suf}
	
	# Resample base to source (so are same grid) - may not be needed (as may be done in align step)
		if ( $resamp == 1 ) then
			set resamp_suffix = _resamp	
				3dresample \
					-master ${epi}+orig \
					-input ${anat}+orig \
					-prefix ${anat}_resamp
		else
			set resamp_suffix = 
		endif
		
		# Non-linear warping
		3dQwarp \
			-source ${anat}${resamp_suffix}+orig.HEAD \
			-base ${epi}+orig \
			-prefix ${anat}_3dQ \
			-lpc -maxlev 0 -verb \
			-iwarp

		rm ${anat}_resamp+orig*
	endif

		if ( $warpdir == 2 ) then
			# Epi to anat --> use inverted warp (as always calculate warp with epi as base (in other script) to make partial volume less of an issue)
			set warp = ${anat}_3dQ_WARPINV+orig
		else
			# Anat to epi --> regular warp, see above
			set warp = ${anat}_3dQ_WARP+orig
		endif

		# Apply warp
		3dNwarpApply \
					-nwarp $warp  \
					-source ${epi}+orig \
					-prefix ${epi}_3dQ_applied
	endif

Code to apply warp to the epi runs + time

# Set warp direction: 1 = anat to epi; 2 = epi to anat
set warpdir = 2

# Set exp suf
set exp_suf = .fm

# Set runs
set runs = ( 01 02 03 04 )

	## Set up
	# set sub = sub-005
	set ses = ses-01
	set subj = ${sub}_${ses}

	## Set inputs - thing you want to warp should be like what you are warping to (i.e., SS if warp to epi - also no skull; AND pre-aligned)
	set anat = anat_final.$subj${exp_suf}

	## Set warp based on direction
	if ( $warpdir == 2 ) then
			# Epi to anat --> use inverted warp (as always calculate warp with epi as base (in other script) to make partial volume less of an issue)
			set warp = ${anat}_3dQ_WARPINV+orig
		else
			# Anat to epi --> regular warp, see above
			set warp = ${anat}_3dQ_WARP+orig
	endif

	foreach run ( $runs )
		set epi = pb02.${subj}.r${run}.volreg${exp_suf}

		# Apply warp
		3dNwarpApply \
					-nwarp ${anat}_3dQ_WARPINV+orig  \
					-source ${epi}+orig \
					-prefix ${epi}_3dQ_applied
	
	end