Hi afni community.
I have a bizarre issue that I want to submit you. I was executing an afni proc py script, and I needed to import and intersect two masks.
afni_proc.py \
-subj_id $subjid \
-blocks despike ricor tshift align tlrc volreg mask blur scale regress \
-radial_correlate_blocks tcat volreg regress \
-copy_anat $subjid.anat.stripped.unifize.nii.gz \
-ROI_import Vi Vent_mask+tlrc \
-ROI_import WMi WM_mask+tlrc \
-anat_has_skull no \
-dsets $subjid.epi+orig.HEAD \
-blip_forward_dset $subjid.epi+orig.HEAD \
-blip_reverse_dset $subjid.pepolar+orig.HEAD \
-ricor_regs physio/$subjid.slibase.1D \
-ricor_regress_method per-run \
-align_unifize_epi local \
-align_opts_aea -cost lpc+ZZ \
-giant_move \
-tlrc_base MNI152_2009_template.nii.gz \
-tlrc_NL_warp \
-volreg_align_to MIN_OUTLIER \
-volreg_align_e2a \
-volreg_tlrc_warp \
-volreg_compute_tsnr yes \
-mask_epi_anat yes \
-mask_segment_anat yes \
-mask_segment_erode yes \
-mask_intersect Vf CSFe Vi \
-mask_intersect WMf WMe WMi \
-blur_size 3 \
-regress_censor_motion 0.2 \
-regress_censor_outliers 0.05 \
-regress_motion_per_run \
-regress_ROI Vf \
-regress_ROI WMf \
-regress_bandpass 0.01 0.1 \
-regress_make_corr_vols WMf Vf \
-regress_apply_mot_types demean deriv \
-regress_opts_3dD -jobs 4 \
-regress_3dD_stop \
-regress_reml_exec \
-regress_est_blur_epits \
-regress_est_blur_errts \
-html_review_style pythonic
Where Vent_mask+tlrc and WM_mask+tlrc are two masks that are on the MNI space (so fitted over the MNI2009 template), with 1x1x1mm, 256 256 256 grid. They don't need a warp, need just to be resampled with the epi grid before doing all the stuff.
The reason for I used -ROI_import is that it promises that:
ROI_import LABEL RSET : import a final grid ROI with the given labele.g. -ROI_import Glasser MNI_Glasser_HCP_v1.0.nii.gz
e.g. -ROI_import Benny my_habenula_rois.nii.gz
e.g. -ROI_import Benny path/to/ROIs/my_habenula_rois.nii.gz
Use this option to import an ROI dataset that is in the final space of
the EPI data. It will merely be resampled onto the final EPI grid
(not transformed).
o this might be based on the group template
o no warping will be done to this dataset
o this dataset WILL be resampled to match the final EPI
So I expected that everything went ok. However, after all the transformations (and 2h of execution) the script ended saying that "ROI_import_Vi does not exist".
Looking at the proc script things become even more bizarre... why it seems to try to use the Vi mask before resampling it? And why it forgets a "+tlrc" (which probably causes the not found error)?
Please send help!
Thnks
Ms
The proc (i removed some blocks lines like tcat and blip because there were too many lines):
#!/bin/tcsh -xef
echo "auto-generated by afni_proc.py, Wed Oct 9 16:58:18 2024"
echo "(version 7.79, Aug 30, 2024)"
echo "execution started: `date`"
# to execute via tcsh:
# tcsh -xef proc.ct10 |& tee output.proc.ct10
# to execute via bash:
# tcsh -xef proc.ct10 2>&1 | tee output.proc.ct10
# =========================== 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 = ct10
endif
# 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
endif
# set list of runs
set runs = (`count_afni -digits 2 1 1`)
# create results and stimuli directories
mkdir -p $output_dir
mkdir $output_dir/stimuli
# copy anatomy to results dir
3dcopy ct10.anat.stripped.unifize.nii.gz \
$output_dir/ct10.anat.stripped.unifize
# copy template to results dir (for QC)
3dcopy /home/martin/abin/MNI152_2009_template.nii.gz \
$output_dir/MNI152_2009_template.nii.gz
# copy any -ROI_import datasets as ROI_import_LABEL
3dcopy Vent_mask+tlrc $output_dir/ROI_import_Vi
3dcopy WM_mask+tlrc $output_dir/ROI_import_WMi
3dcopy /home/martin/abin/APQC_atlas_MNI_2009c_asym.nii.gz \
$output_dir/ROI_import_MNI_2009c_asym
# copy external -blip_forward_dset dataset
3dTcat -prefix $output_dir/blip_forward ct10.epi+orig
# copy external -blip_reverse_dset dataset
3dTcat -prefix $output_dir/blip_reverse ct10.pepolar+orig
# copy slice-based regressors for RETROICOR (rm first 0 TRs)
1dcat physio/ct10.slibase.1D > $output_dir/stimuli/ricor_orig_r01.1D
# ================================= align ==================================
# for e2a: compute anat alignment transformation to EPI registration base
# (new anat will be current ct10.anat.stripped.unifize+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 ct10.anat.stripped.unifize+orig \
-suffix _al_junk \
-epi vr_base_min_outlier_unif+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 (non-linear warp)
auto_warp.py -base MNI152_2009_template.nii.gz -input \
ct10.anat.stripped.unifize+orig \
-skull_strip_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 ct10.anat.stripped.unifize \
awpy/ct10.anat.stripped.unifize.aw.nii*
mv awpy/anat.un.aff.Xat.1D .
mv awpy/anat.un.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 ct10.anat.stripped.unifize+tlrc.HEAD ) then
echo "** missing +tlrc warp dataset: \
ct10.anat.stripped.unifize+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 \
pb04.$subj.r$run.blip+orig
# create an all-1 dataset to mask the extents of the warp
3dcalc -overwrite -a pb04.$subj.r$run.blip+orig -expr 1 \
-prefix rm.epi.all1
# catenate blip/volreg/epi2anat/tlrc xforms
cat_matvec -ONELINE \
anat.un.aff.Xat.1D \
ct10.anat.stripped.unifize_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 ct10.anat.stripped.unifize+tlrc -dxyz 1.5 \
-source pb03.$subj.r$run.tshift+orig \
-nwarp "anat.un.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 ct10.anat.stripped.unifize+tlrc -dxyz 1.5 \
-source rm.epi.all1+orig \
-nwarp "anat.un.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)
# (only 1 run, so just use 3dcopy to keep naming straight)
3dcopy rm.epi.min.r01+tlrc 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 pb05.$subj.r$run.volreg
end
# warp the volreg base EPI dataset to make a final version
cat_matvec -ONELINE \
anat.un.aff.Xat.1D \
ct10.anat.stripped.unifize_al_junk_mat.aff12.1D -I > \
mat.basewarp.aff12.1D
3dNwarpApply -master ct10.anat.stripped.unifize+tlrc -dxyz 1.5 \
-source vr_base_min_outlier+orig \
-nwarp "anat.un.aff.qw_WARP.nii mat.basewarp.aff12.1D" \
-prefix final_epi_vr_base_min_outlier
# create an anat_final dataset, aligned with stats
3dcopy ct10.anat.stripped.unifize+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
# --------------------------------------
# create a TSNR dataset, just from run 1
3dTstat -mean -prefix rm.signal.vreg.r01 pb05.$subj.r01.volreg+tlrc
3dDetrend -polort 4 -prefix rm.noise.det -overwrite pb05.$subj.r01.volreg+tlrc
3dTstat -stdev -prefix rm.noise.vreg.r01 rm.noise.det+tlrc
3dcalc -a rm.signal.vreg.r01+tlrc \
-b rm.noise.vreg.r01+tlrc \
-c mask_epi_extents+tlrc \
-expr 'c*a/b' -prefix TSNR.vreg.r01.$subj
# ---------------------------------------------------------
# QC: compute correlations with spherical ~averages
@radial_correlate -nfirst 0 -polort 4 -do_clean yes \
-rdir radcor.pb05.volreg \
pb05.$subj.r*.volreg+tlrc.HEAD
# ================================== mask ==================================
# create 'full_mask' dataset (union mask)
foreach run ( $runs )
3dAutomask -prefix rm.mask_r$run pb05.$subj.r$run.volreg+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 \
ct10.anat.stripped.unifize+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.nii.gz)
3dresample -master mask_epi_anat.$subj+tlrc -prefix ./rm.resam.group \
-input /home/martin/abin/MNI152_2009_template.nii.gz
# 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
# ---- segment anatomy into classes CSF/GM/WM ----
3dSeg -anat anat_final.$subj+tlrc -mask AUTO -classes 'CSF ; GM ; WM'
# copy resulting Classes dataset to current directory
3dcopy Segsy/Classes+tlrc .
# make individual ROI masks for regression (CSF GM WM and CSFe GMe WMe)
foreach class ( CSF GM WM )
# unitize and resample individual class mask from composite
3dmask_tool -input Segsy/Classes+tlrc"<$class>" \
-prefix rm.mask_${class}
3dresample -master pb05.$subj.r01.volreg+tlrc -rmode NN \
-input rm.mask_${class}+tlrc -prefix mask_${class}_resam
# also, generate eroded masks
3dmask_tool -input Segsy/Classes+tlrc"<$class>" -dilate_input -1 \
-prefix rm.mask_${class}e
3dresample -master pb05.$subj.r01.volreg+tlrc -rmode NN \
-input rm.mask_${class}e+tlrc -prefix mask_${class}e_resam
end
# create intersect mask Vf from masks CSFe and Vi
3dmask_tool -input mask_CSFe_resam+tlrc ROI_import_Vi \
-inter -prefix mask_inter_Vf
# create intersect mask WMf from masks WMe and WMi
3dmask_tool -input mask_WMe_resam+tlrc ROI_import_WMi \
-inter -prefix mask_inter_WMf
# ================================== blur ==================================
# blur each volume of each run
foreach run ( $runs )
3dmerge -1blur_fwhm 3 -doall -prefix pb06.$subj.r$run.blur \
pb05.$subj.r$run.volreg+tlrc
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 run ( $runs )
3dTstat -prefix rm.mean_r$run pb06.$subj.r$run.blur+tlrc
3dcalc -a pb06.$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 pb07.$subj.r$run.scale
end
# ================================ regress =================================
# compute de-meaned motion parameters (for use in regression)
1d_tool.py -infile dfile_rall.1D -set_nruns 1 \
-demean -write motion_demean.1D
# compute motion parameter derivatives (for use in regression)
1d_tool.py -infile dfile_rall.1D -set_nruns 1 \
-derivative -demean -write motion_deriv.1D
# convert motion parameters for per-run regression
1d_tool.py -infile motion_demean.1D -set_nruns 1 \
-split_into_pad_runs mot_demean
1d_tool.py -infile motion_deriv.1D -set_nruns 1 \
-split_into_pad_runs mot_deriv
# create censor file motion_${subj}_censor.1D, for censoring motion
1d_tool.py -infile dfile_rall.1D -set_nruns 1 \
-show_censor_count -censor_prev_TR \
-censor_motion 0.2 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
# create bandpass regressors (instead of using 3dBandpass, say)
1dBport -nodata 160 3.0 -band 0.01 0.1 -invert -nozero > bandpass_rall.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_Vi+tlrc \
-prefix ROI_import_Vi_resam
3drefit -copytables ROI_import_Vi+tlrc ROI_import_Vi_resam+tlrc
3drefit -cmap INT_CMAP ROI_import_Vi_resam+tlrc
3dresample -master final_epi_vr_base_min_outlier+tlrc -rmode NN \
-input ROI_import_WMi+tlrc \
-prefix ROI_import_WMi_resam
3drefit -copytables ROI_import_WMi+tlrc ROI_import_WMi_resam+tlrc
3drefit -cmap INT_CMAP ROI_import_WMi_resam+tlrc
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
# ------------------------------
# create 2 ROI regressors: Vf, WMf
# (get each ROI average time series and remove resulting mean)
foreach run ( $runs )
3dmaskave -quiet -mask mask_inter_Vf+tlrc \
pb05.$subj.r$run.volreg+tlrc \
| 1d_tool.py -infile - -demean -write rm.ROI.Vf.r$run.1D
3dmaskave -quiet -mask mask_inter_WMf+tlrc \
pb05.$subj.r$run.volreg+tlrc \
| 1d_tool.py -infile - -demean -write rm.ROI.WMf.r$run.1D
end
# and catenate the demeaned ROI averages across runs
cat rm.ROI.Vf.r*.1D > ROI.Vf_rall.1D
cat rm.ROI.WMf.r*.1D > ROI.WMf_rall.1D
# 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 pb07.$subj.r*.scale+tlrc.HEAD \
-censor censor_${subj}_combined_2.1D \
-ortvec bandpass_rall.1D bandpass \
-ortvec ROI.Vf_rall.1D ROI.Vf \
-ortvec ROI.WMf_rall.1D ROI.WMf \
-ortvec mot_demean.r01.1D mot_demean_r01 \
-ortvec mot_deriv.r01.1D mot_deriv_r01 \
-polort 4 \
-num_stimts 0 \
-jobs 4 \
-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
# 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
# -- 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...'
exit
endif
# create an all_runs dataset to match the fitts, errts, etc.
3dTcat -prefix all_runs.$subj pb07.$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 \
-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}_REML+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}_REML+tlrc > mean.errts.1D
3dTcorr1D -prefix corr_brain errts.${subj}_REML+tlrc mean.errts.1D
# compute 2 requested correlation volume(s)
# create correlation volume corr_af_WMf
3dcalc -a mask_inter_WMf+tlrc -b mask_epi_anat.$subj+tlrc -expr 'a*b' \
-prefix rm.fm.WMf
3dmaskave -q -mask rm.fm.WMf+tlrc \
errts.${subj}_REML+tlrc > mean.ROI.WMf.1D
3dTcorr1D -prefix corr_af_WMf errts.${subj}_REML+tlrc mean.ROI.WMf.1D
# create correlation volume corr_af_Vf
3dcalc -a mask_inter_Vf+tlrc -b mask_epi_anat.$subj+tlrc -expr 'a*b' \
-prefix rm.fm.Vf
3dmaskave -q -mask rm.fm.Vf+tlrc \
errts.${subj}_REML+tlrc > mean.ROI.Vf.1D
3dTcorr1D -prefix corr_af_Vf errts.${subj}_REML+tlrc mean.ROI.Vf.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 err_reml --
touch blur.err_reml.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.err_reml.r$run.1D \
errts.${subj}_REML+tlrc"[$trs]" >> blur.err_reml.1D
end
# 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 mask_epi_anat.$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}_REML+tlrc
# ---------------------------------------------------------
# QC: compute correlations with spherical ~averages
@radial_correlate -nfirst 0 -polort 4 -do_clean yes \
-rdir radcor.pb08.regress \
-mask mask_epi_anat.$subj+tlrc \
all_runs.$subj+tlrc.HEAD errts.${subj}_REML+tlrc.HEAD
# ========================= 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.2
out_limit : 0.05
copy_anat : ct10.anat.stripped.unifize+orig.HEAD
errts_dset : errts.${subj}_REML+tlrc.HEAD
mask_dset : mask_epi_anat.$subj+tlrc.HEAD
template : MNI152_2009_template.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 Segsy
# --------------------------------------------------
# if the basic subject review script is here, run it
# (want this to be the last text output)
if ( -e @ss_review_basic ) then
tcsh @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 $subj.results/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 ct10 -blocks despike ricor tshift align tlrc volreg \
# mask blur scale regress -radial_correlate_blocks tcat volreg regress \
# -copy_anat ct10.anat.stripped.unifize.nii.gz -ROI_import Vi \
# Vent_mask+tlrc -ROI_import WMi WM_mask+tlrc -anat_has_skull no -dsets \
# ct10.epi+orig.HEAD -blip_forward_dset ct10.epi+orig.HEAD \
# -blip_reverse_dset ct10.pepolar+orig.HEAD -ricor_regs \
# physio/ct10.slibase.1D -ricor_regress_method per-run -align_unifize_epi \
# local -align_opts_aea -cost lpc+ZZ -giant_move -tlrc_base \
# MNI152_2009_template.nii.gz -tlrc_NL_warp -volreg_align_to MIN_OUTLIER \
# -volreg_align_e2a -volreg_tlrc_warp -volreg_compute_tsnr yes \
# -mask_epi_anat yes -mask_segment_anat yes -mask_segment_erode yes \
# -mask_intersect Vf CSFe Vi -mask_intersect WMf WMe WMi -blur_size 3 \
# -regress_censor_motion 0.2 -regress_censor_outliers 0.05 \
# -regress_motion_per_run -regress_ROI Vf -regress_ROI WMf \
# -regress_bandpass 0.01 0.1 -regress_make_corr_vols WMf Vf \
# -regress_apply_mot_types demean deriv -regress_opts_3dD -jobs 4 \
# -regress_3dD_stop -regress_reml_exec -regress_est_blur_epits \
# -regress_est_blur_errts -html_review_style pythonic
Howdy, Ms-
I think you need to use -mask_import
instead of ROI_import
to be able to do the intersection.
# AP help on the mask_import option, including an example of mask intersection
-mask_import LABEL MSET : import a final grid mask with the given label
e.g. -mask_import Tvent template_ventricle_3mm+tlrc
* Note: -ROI_import basically makes -mask_import unnecessary.
Use this option to import a mask that is aligned with the final
EPI data _and_ is on the final grid (with -ROI_import, the ROI will
be resampled onto the final grid).
o this might be based on the group template
o this should already be resampled appropriately
o no warping or resampling will be done to this dataset
This mask can be applied via LABEL as other masks, using options
like: -regress_ROI, -regress_ROI_PC, -regress_make_corr_vols,
-regress_anaticor_label, -mask_intersect, -mask_union.
For example, one might import a ventricle mask from the template,
intersect it with the subject specific CSFe (eroded CSF) mask,
and possibly take the union with WMe (eroded white matter), before
using the result for principle component regression, as in:
-mask_import Tvent template_ventricle_3mm+tlrc \
-mask_intersect Svent CSFe Tvent \
-mask_union WM_vent Svent WMe \
-regress_ROI_PC WM_vent 3 \
--pt
rickr
October 10, 2024, 8:00pm
4
Hi Ms,
I think the problem has been fixed, but for a reminder to you and myself, I should have initially asked for the actual lines containing the error message. It took a while for me to realize the context of the error.
Do you build your own binaries or download a precompiled package? We can run a build tonight.
Thanks,
Hi Rick, Hi PTaylor
I'm sorry if my message has been too foggy. I try to explain more clearly.
When I executed the afni proc script reported above (so running the .proc file generated for the processing), the program terminated at some time reporting me an
ERROR: ROI_import_Vi not found
(something similar)
Since I didn't get what it was happening I analyzed the processing script above, and I found two problems.
E1: SYNTAX ERROR.
The script imports at line 58 the Vent_mask+tlrc
using the code: 3dcopy Vent_mask+tlrc $output_dir/ROI_import_Vi
. If I understood enough the AFNI behavior, this would mean that the mask is saved as ROI_import_Vi+tlrc
in the output folder.
So, when at line 269 (line numbers do match a real proc script because I deleted the tcat and blip) the program calls
3dmask_tool -input mask_CSFe_resam+tlrc ROI_import_Vi -inter -prefix mask_inter_Vf
The output is "ROI_import_Vi" not found. That is because it should be "+tlrc" here!
E2: CONCEPTUAL ERROR.
Even if you try to manually correct the above E1 and re-execute the code (I tried actually ), then the program will terminate at the same line 269 raising a Grid mismatch error
. In fact, looking at the same line, you see that the program is attempting to intersect the mask_CSFe_resam+tlrc
(that is on the epi grid) and the ROI_import_Vi(+tlrc!)
(that is on the MNI original grid!). While, for some unknown reason, the program resamples the ROI_import_Vi
only at line 334 with
3dresample -master final_epi_vr_base_min_outlier+tlrc -rmode NN -input ROI_import_Vi+tlrc -prefix ROI_import_Vi_resam
. But I think that it is this mask, the resamp version, that should be intersect with the CSFe mask at line 269! (Otherwise, why should one have used the option to import a "to-be-resampled" mask?)
Moral of the story: I overcome this by simply throwing away the -ROI_import
option and using instead -mask_import
(as PTaylor suggested me). Of course this meant I needed to resample the two masks to the epi grid outside and before the execution of afni proc, using a proxy grid that I exported manually.
But, I wanted to report these issues to you because it seems to me a pathological behavior of afni proc script, and maybe you could do something.
I downloaded afni using the standard procedure (Quick Setup) for Ubuntu22.04: 1.1.3. Linux, Ubuntu 22.04 — AFNI, SUMA and FATCAT: v24.3.02 . I think the afni version
that I used is the most update (I run the "update -defaults" at the beginning of this week so...)
Thank you for your support,
Ms
rickr
October 11, 2024, 3:28pm
6
Hi Ms,
You may have not seen my short message, but let me clarify the points either way:
I think that I fixed this problem yesterday, and we rebuilt the binaries last night.
It would be helpful for you to update your binaries and rerun the original command, to see if the fix is good for your command.
For reporting errors, I just meant that it is important to show commands (as you did) and the actual copy-and-pasted error messages (not done). Your interpretation is appreciated, but the entire messages are what we really want.
The question about whether you build your own or download was asked because if you build your own, that could have been done right then. But it doesn't matter, since you seem to download binaries.
Regarding -mask_import, that option indeed requires the mask to be on the same grid as the final output (and it does say so in the help). But we are moving away from -mask_import to the newer -ROI_import. So your original use was good.
Anyway, now that there are updated binaries, please feel free to update and try the original command.
-rick
Hi Rick.
I appreciate the support you are giving me and I apologize if my above message may look unkind or verbose. Sometimes my enthusiasm takes over. I will try to report my issues more concisely hereafter.
I updated the afni binaries, and ran again the script. It ended with the following error:
3dmask_tool -input mask_CSFe_resam+tlrc ROI_import_Vi+tlrc -inter -prefix mask_inter_Vf
++ processing 2 input dataset(s), NN=2...
++ padding all datasets by 0 (for dilations)
** FATAL ERROR: nvoxel mismatch
** Program compile date = Oct 10 2024
Thanks,
Ms
rickr
October 17, 2024, 1:35pm
8
Hi Ms,
Just so you know, I did replicate this a couple of days ago, but have not had time to fix it yet - it will happen. The problem is that ROI_import resampling is currently done in the regress block, while this intersect step is in the mask block. One of those will have to move.
Anyway, it will get fixed, hopefully soon. Sorry for the delay.