s3.rest_surf_preproc

#!/bin/tcsh
# This script will complete all functional preprocessing and send data to surface for resting state
# Steps include: Slice time correction, Motion correction, Co-registration, Tissue extraction, Vol2Surf, Nuisance regression, Intensity normalization, Surface smoothing
#
#<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
# Script assumes that T1+orig and functional runs +orig are in the data directory
# Also assumes that @SUMA_Make_Spec_FS has been run for surface files
#
# Assumed directory structure:
# data_dir/
#	subject_ID/
#		bold/
#		   subID_T1_+orig
#		   subID_run_+orig
# SUMA_dir/
#	subject_ID/
#		SUMA/
#
#<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>

# Enter subjects
set subj 		= ()

# Set functional run identifier (i.e. TOM; BM; LL)      !! Do not enter anatomical !!
# same ones used to grab the correct raw dicom files and to name individual runs in s1.dataset
set runs 		= ()

# Set destination and study directory
set data_dir		= 
set SUMA_dir 		= 

# log file will contain errors and a history of runs processed
set log_name		= preproc.log
set log_dir		= 
set ow_log		= no			# do you want to overwrite previous log?
						# if no, script will append to the previous log file
# Set processing Options
set default_vol 	= yes			# yes 	-> default volreg aligns everything to the 1st volume of the 1st run
# if no:					# no 	-> allows you to specify what run and volume below
set volreg_run		= 			# name of run that contains the volume you want to register everything to (can be blank if default.vol == yes)
set volreg_vol 		= 			# volume number to register everything to (index begin at 0), (can be blank if default.vol == yes)

set resolution 		= 3 			# resolution (in mm) that you want the functional grid (used in co-registration)
set censor 		= (1)			# Motion censor limit (in mm) for nuisance regression, further motion parameters can be assessed with s4.motion
set band 		= (0.01 .1)		# Bandpass Filter
set fwhm 		= 6			# FWHM blur (in mm) for surface smoothing

# Set script options
set overwrite 		= yes			# yes 	-> delete conflicting files and create anew
						# no 	-> if conflicting file exists, proc step skipped 
set clean_output 	= no			# yes 	-> clean output, only essential final files will be saved
						# no 	-> all files will be saved, can be used for troubleshooting
set want_tsnr 		= no 			# yes 	-> will output a tsnr dataset based on the volreg step, and TSNR.$run.out which = mean TSNR text file
set skip_bad_motion	= no 			# yes 	-> calculate percent of censored TRs based on below parameters and will skip 
# if yes 						   if above the threshold
set bad_percent 	= 			# motion threshold, runs with percent TRs censored > bad_percent will be skipped

# known issues to be solved later:


#<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
#<><><><><><><><><><><><><><><> Do Not Change  <><><><><><><><><><><><><><><><><>
#<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>


# Set up log file
if (-d $log_dir) then
	set log = $log_dir/$log_name
else
	echo " "
	echo "Specified log directory does not exist, I will put log in:" ~
	echo " "
	set log = ~/$log_name
endif
if ($ow_log == yes && -f $log) \rm -f $log
set now = `date +"%a %b %d %Y %r"`
if (! -f $log) then
	touch $log
else
	echo " " >> $log
endif
echo "==================================================================" >> $log
echo "s3.rest_surf_preproc ran by $user on $now " >> $log
echo "==================================================================" >> $log
# Directory check
if (! -d $data_dir && ! -d $SUMA_dir) then
	echo " "
	echo "Data and/or SUMA directory does not exist, script failing..." | tee -a $log
	echo "It's a trap! - Admiral Ackbar" | tee -a $log
	echo " "
	exit
else
	cd $data_dir
endif
# prepare file name for clean output
if ($clean_output == yes) then
	set o = "rm."
else if ($clean_output == no) then
	set o = ""
else
	echo " "
	echo "WARNING: clean_out option $clean_output not understood, setting to no" | tee -a $log
	echo " "
	set o = ""
endif
set nonomatch
# Subject loop
foreach s ($subj)
	echo " "
	echo "========================================================================="
	echo "========================================================================="
	echo " "
	echo "Current Subject: $s"
	echo " "
# Can I find the current subject?
	echo " " >> $log 
	echo $s >> $log
	echo "------" >> $log
	if (! -d $s/bold) then
		echo "Cannot find subject $s" | tee -a $log
		echo "Next subject..."
	else
		cd $s/bold
# Cycle through functional runs and grab runs that exist for current subject
#m = missing, t = temp, d = done (already), c = current
		set m_runs = ()
		set t_runs = ()
		set d_runs = () 
		set c_runs = ()
		foreach r ($runs)
			if (! -f ${s}_${r}+orig.HEAD) then
				echo "Missing run: $r"
				echo " "
				set m_runs = ($m_runs $r)
			else
				set t_runs = ($t_runs $r)
			endif
		end
		if ($#m_runs != 0) echo "Missing:" $m_runs >> $log
# remove conflicting files, take note of files already processed, and create a list of current files to process
		if ($overwrite == yes) then
			echo " "
			echo "Removing previous files"
			echo " "
			foreach r ($t_runs)
				\rm -f *"$r".pb* *"$r"*.mat* *"$r".MoPar* *mask* *rm* *T1_al* *T1_ns* *anat_final* *FullMask* *TSNR* *aff12* *mLV* *mWM* *SurfVol* *aseg* *censor* *CENSORTR* *A2E* *shft* *REML* *xmat* *.jpg*
			end
		endif
		foreach r ($t_runs)
			if (-f ${s}_$r.pb07.lh.smooth.niml.dset && -f ${s}_$r.pb07.rh.smooth.niml.dset) then
				set d_runs = ($d_runs $r)
			else
				set c_runs = ($c_runs $r)
			endif
		end
		if ($#d_runs != 0) then
			echo "Subject ${s}, runs: ${d_runs} already exist, I won't waste your time"
			echo "Already exist:" $d_runs >> $log
		endif
# Process current runs, if there are any
		if ($#c_runs != 0) then
# Volume based functional preprocessing
# Slice time correction
			echo " "
			echo "========================= Slice Time Correction ========================="
			echo " "
			foreach r ($c_runs)
				echo " "
				echo "Current run: $r"
				echo " "
				if (! -f ${o}${s}_$r.pb01.Tshift+orig.HEAD) then
					3dTshift 							\
						-prefix ${o}${s}_$r.pb01.Tshift 			\
						${s}_$r+orig
					if (! -f ${o}${s}_$r.pb01.Tshift+orig.HEAD) then
						echo "slice timing, run: $r failed..." >> $log
					endif
				else
					echo "${s}_$r.pb01.Tshift already exists"
				endif
			end
			echo " "
			echo "=========================== Volume Registration ==========================="
			echo " "
# Prepare proper volume for registration
			if ($default_vol == yes) then
				set v_run = ${o}${s}_$c_runs[1].pb01.Tshift+orig
				set v_vol = 0
			else if ($default_vol == no) then
				set v_run = ${o}${s}_$volreg_run.pb01.Tshift+orig
				set v_vol = $volreg_vol
			else
				echo " "
				echo "Default option $default_vol not understood, using defaults..." | tee -a $log
				echo " "
				set v_run = ${o}${s}_$c_runs[1].pb01.Tshift+orig
				set v_vol = 0
			endif
			echo "base: ${v_run}:${v_vol}" >> $log
# Align anatomy to epi base
			if (! -f ${s}_T1_ns+orig.HEAD) then
				align_epi_anat.py 										\
					-anat2epi -anat ${s}_T1+orig 								\
					-save_skullstrip 									\
					-suffix _al 										\
					-epi $v_run -align_centers yes								\
					-epi_strip 3dAutomask 									\
					-epi_base $v_vol -volreg off -tshift off
				if (! -f ${s}_T1_ns+orig.HEAD) then
					echo "align_epi_anat failed..." >> $log
				endif
			else
				echo " "
				echo "${s}_T1_ns already exists"
				echo " "
			endif
# Put everything into the same space 
			foreach r ($c_runs)
				echo " "		
				echo "Current run: $r"
				echo " "
				if (! -f ${o}${s}_$r.pb02.volreg+orig.HEAD) then
# Align all functional slices to same base and save transformations
					3dvolreg 										\
						-base $v_run"[$v_vol]" 								\
						-1Dfile ${s}_$r.MoPar.1D 							\
						-prefix rm.$r.volreg 								\
						-1Dmatrix_save ${o}${s}_$r.volreg.aff12.1D 					\
						-zpad 1 	 								\
						${o}${s}_$r.pb01.Tshift+orig
# Create an extent mask of base dataset
					3dcalc 											\
						-overwrite -a ${o}${s}_$r.pb01.Tshift+orig 					\
						-expr 1 -prefix rm.epi.all1
# Cat transformation parameters to apply to epi
					cat_matvec -ONELINE 									\
						${s}_T1_al_mat.aff12.1D -I 							\
						${o}${s}_$r.volreg.aff12.1D > ${o}${s}_$r.warp.mat.1D
# Apply transformation parameters to epi datasets (create temporary files until process is complete)
					3dAllineate 										\
						-base ${s}_T1_ns+orig 								\
						-input ${o}${s}_$r.pb01.Tshift+orig 						\
						-1Dmatrix_apply ${o}${s}_$r.warp.mat.1D 					\
						-mast_dxyz $resolution 								\
						-prefix rm.$r.nomask
# Apply transformations to mask as well
					3dAllineate 										\
						-base ${s}_T1_ns+orig 								\
						-input rm.epi.all1+orig 							\
						-1Dmatrix_apply ${o}${s}_$r.warp.mat.1D 					\
						-mast_dxyz $resolution -final NN -quiet 					\
						-prefix rm.$r.all1
# extents intersection
					3dTstat 										\
						-min -prefix rm.$r.min 								\
						rm.$r.all1+orig
				else
					echo " "
					echo "${s}_$r has already been aligned"
					echo " "
				endif
			end
			foreach r ($c_runs)
				if (! -f ${o}${s}_$r.pb02.volreg+orig.HEAD) then
					3dcalc 											\
						-a rm.$r.nomask+orig 								\
						-b rm.$r.all1+orig 								\
						-expr 'a*b' 									\
						-prefix ${o}${s}_$r.pb02.volreg
					if (! -f ${o}${s}_$r.pb02.volreg+orig.HEAD) then
						echo "volreg, run: $r failed..." >> $log
					endif
				else
					echo " "
					echo "${s}_$r.pb02.volreg+orig already exists"
					echo " "
				endif
			end
# test and skip bad motion if you want me to
			if ($skip_bad_motion == yes) then
				alias calc 'awk "BEGIN{ print \!* }" '
				set bm_runs = ()
				set gm_runs = ()
				set bm_perc = ()
				foreach r ($c_runs)
					set tot_tr = `1d_tool.py -infile ${s}_$r.MoPar.1D -show_tr_run_counts trs`
					set cen_tr = `1d_tool.py -infile ${s}_$r.MoPar.1D -quick_censor_count $censor`
					set perc = `calc int(($cen_tr / $tot_tr) * 100)`
					if ($perc >= $bad_percent) then
						set bm_runs = ($r $bm_runs)
						set bm_perc = ($perc $bm_perc)
					else
						set gm_runs = ($r $gm_runs)
					endif
				end
				if ($#bm_runs != 0) then
					echo "Skipped due to motion with ${censor}mm censor:" >> $log
					set i = 1
					set n = $#bm_perc
					while ($i <= $n)
						echo "  ${bm_runs[$i]}: ${bm_perc[$i]}% censored" >> $log
						@ i++
					end
				endif
				set c_runs = ($gm_runs)
				if ($#c_runs == 0) then
					echo " "
					echo "No runs survived motion threshold for ${s}, next..."
					echo " "
					continue
					continue
					continue
					continue
				endif
			endif
# Send functional data to surface
			echo " "
			echo "=================== Align Surface Volume to Experiment =================="
			echo " "
			if (! -f ${s}_SurfVol+orig.HEAD) then
				@SUMA_AlignToExperiment 									\
					-exp_anat ${s}_T1_ns+orig 								\
					-surf_anat $SUMA_dir/$s/SUMA/${s}_SurfVol+orig 						\
					-wd -align_centers 									\
					-prefix ${s}_SurfVol 									\
					-surf_anat_followers $SUMA_dir/$s/SUMA/aseg.nii
				if (! -f ${s}_SurfVol+orig.HEAD) then
					echo "AlignToExperiment failed..." >> $log
				endif
			else
				echo " "
				echo "${s}_SurfVol already exists"
				echo " "
			endif
# White matter and ventricle nuisance regressors
			echo " "
			echo "=================== Create and Re-sample Tissue Masks  =================="
			if (! -f ${s}_mWM_erode+orig.HEAD && ! -f ${s}_mLV_erode+orig.HEAD) then
				foreach t (mWM mLV)
					if ($t == mWM) then
						echo " "
						echo "White matter..."
						echo " "
							3dcalc 									\
								-a aseg_Alnd_Exp+orig 						\
								-expr 'equals(a,2)+equals(a,7)+equals(a,41)+equals(a,46)+equals(a,251)+equals(a,252)+equals(a,253)+equals(a,254)+equals(a,255)+equals(a,16)' \
								-prefix ${o}${s}_$t
							3dresample 								\
								-master ${o}${s}_$c_runs[1].pb02.volreg+orig 			\
								-rmode NN 							\
								-inset ${o}${s}_${t}+orig 					\
								-prefix ${o}${s}_${t}_re
						echo " "
						echo "I'm eroding..."
						echo " "
							3dcalc 									\
								-a ${o}${s}_${t}_re+orig 					\
								-prefix ${s}_${t}_erode 					\
								-b a+i -c a-i -d a+j -e a-j -f a+k -g a-k 			\
								-expr 'a*(1-amongst(0,a,b,c,d,e,f,g))'
						if (! -f ${s}_${t}_erode+orig.HEAD) then
							echo "white matter mask failed..." >> $log
						endif
					else
						echo " "
						echo "Ventricle..."
						echo " "
						3dcalc 										\
							-a aseg_Alnd_Exp+orig 							\
							-expr 'equals(a,4)+equals(a,43)' 					\
							-prefix ${o}${s}_$t
						echo " "
						echo "I'm eroding..."
						echo " "
							3dfractionize 								\
								-template ${o}${s}_$c_runs[1].pb02.volreg+orig 			\
								-input ${o}${s}_$t+orig 					\
								-prefix ${s}_${t}_erode 					\
								-vote -preserve -clip 0.8
						if (! -f ${s}_${t}_erode+orig.HEAD) then
							echo "ventricle mask failed..." >> $log
						endif
					endif
				end
			else
				echo " "
				echo "${s} tissue masks already exist"
				echo " "
			endif
			echo " "
			echo "==================== Extract Time Series from Masks  ===================="			
			foreach r ($c_runs)
				echo " "
				echo "Current run: $r"
				echo " "
				foreach t (mWM mLV)
					if ($t == mWM) then
						echo " "
						echo "White matter..."
						echo " "
						if (! -f ${o}${s}_${r}_$t.1D) then
							3dmaskave 								\
								-mask ${s}_${t}_erode+orig -q 					\
								${o}${s}_$r.pb02.volreg+orig 					\
								| 1d_tool.py -infile - -demean -write ${o}${s}_${r}_$t.1D
							if (! -f ${o}${s}_${r}_${t}.1D) then
								echo "time series from white matter, ${r} failed..." >> $log
							endif
						else
							echo " "
							echo "${s}_${r}_$t.1D already exists"
							echo " "
						endif
					else	
						echo " "		
						echo "Ventricle..."
						echo " "
						if (! -f ${o}${s}_${r}_$t.1D) then
							3dmaskave 								\
								-mask ${s}_${t}_erode+orig -q 					\
								${o}${s}_$r.pb02.volreg+orig 					\
								| 1d_tool.py -infile - -demean -write ${o}${s}_${r}_$t.1D
							if (! -f ${o}${s}_${r}_${t}.1D) then
								echo "time series from ventricle, ${r} failed..." >> $log
							endif
						else
							echo " "
							echo "${s}_${r}_$t.1D already exists"
							echo " "
						endif
					endif
				end
			end
# Send to surface
			echo " "
			echo "=========================== Volume to Surface ==========================="
			foreach r ($c_runs)
				echo " "
				echo "Current run: $r"
				echo " "
				foreach h (lh rh)
					if (! -f ${o}${s}_$r.pb03.$h.v2s.niml.dset) then
						3dVol2Surf 									\
							-spec $SUMA_dir/$s/SUMA/std.60.${s}_$h.spec 				\
							-surf_A smoothwm							\
							-surf_B pial 								\
							-sv ${s}_SurfVol+orig 							\
							-grid_parent ${o}${s}_$r.pb02.volreg+orig 				\
							-map_func ave -f_steps 10 -f_index nodes 				\
							-oob_index -1 -oob_value 0.0 						\
							-outcols_NSD_format 							\
							-out_niml ${o}${s}_$r.pb03.$h.v2s.niml.dset
						if (! -f ${o}${s}_$r.pb03.$h.v2s.niml.dset) then
							echo "vol2surf, ${r}:${h} failed..." >> $log
						endif
					else
						echo " "
						echo "${s}_$r.pb03.$h.v2s already exists"
						echo " "
					endif
				end
			end
# Intensity normalization
			echo " "
			echo "======================== Intensity Normalization ========================"
			foreach r ($c_runs)
				echo " "
				echo "Current run: $r"
				echo " "
				foreach h (lh rh)
					if (! -f ${o}${s}_$r.pb04.$h.scale.niml.dset) then
						3dTstat 									\
							-prefix ${o}${s}_$r.pb04.$h.mean.niml.dset				\
							${o}${s}_$r.pb03.$h.v2s.niml.dset
						3dcalc 										\
							-a ${o}${s}_$r.pb03.$h.v2s.niml.dset 					\
							-b ${o}${s}_$r.pb04.$h.mean.niml.dset					\
							-expr 'min(200, a/b*100)*step(a)*step(b)' 				\
							-prefix ${o}${s}_$r.pb04.$h.scale.niml.dset
						if (! -f ${o}${s}_$r.pb04.$h.scale.niml.dset) then
							echo "scaling, ${r}:${h} failed..." >> $log
						endif
					else
						echo " "
						echo "${s}_$r.pb04.$h.scale already exists"
						echo " "
					endif
				end
			end
# Bandpass filter
			echo " "
			echo "=========================== Bandpass Regressor =========================="
			foreach r ($c_runs)
				if (! -f ${o}${s}_$r.band.1D) then
					echo " "
					echo "Current run: $r"
					echo " "
					set tot_tr = `3dinfo -nv ${o}${s}_$r.pb04.lh.scale.niml.dset`
					set tr_len = `3dinfo -tr ${o}${s}_$r.pb04.lh.scale.niml.dset`
					1dBport -nodata $tot_tr $tr_len -band $band -invert -nozero > ${o}${s}_$r.band.1D
					if (! -f ${o}${s}_$r.band.1D) then
						echo "bandpass, ${r}:${h} failed..." >> $log
					endif
				else
					echo " "
					echo "${s}_$r.band.1D already exists"
					echo " "
				endif
			end
# Calculate motion regressors
			echo " "
			echo "=========================== Motion Regressors ==========================="
			foreach r ($c_runs)
			echo " "
			echo "Current run: $r"
			echo " "
				if (! -f ${o}${s}_$r.MoPar_derv.1D && ! -f ${o}${s}_$r.MoPar_demean.1D) then
					1d_tool.py 										\
						-infile ${s}_$r.MoPar.1D 							\
						-derivative 									\
						-write ${o}${s}_$r.MoPar_derv.1D
					1d_tool.py 										\
						-infile ${s}_$r.MoPar.1D 							\
						-demean 									\
						-write ${o}${s}_$r.MoPar_demean.1D
					if (! -f ${o}${s}_$r.MoPar_derv.1D && ! -f ${o}${s}_$r.MoPar_demean.1D) then
						echo "motion parameters, ${r} failed..." >> $log
					endif
				else
					echo " "
					echo "Demean and derivative files for ${s}_$r already exist"
					echo " "
				endif
				foreach c ($censor)
					if (! -f ${o}${s}_${r}.${c}mm_censor.1D) then 
						1d_tool.py 									\
							-infile ${s}_$r.MoPar.1D 						\
							-censor_prev_TR -show_censor_count					\
							-censor_motion $c ${o}${s}_$r.${c}mm
						if (! -f ${o}${s}_${r}.${c}mm_censor.1D) then
						echo "censoring at ${c}mm failed..." >> $log
						endif
					else
						echo " "
						echo "Censoring at ${c}mm failed..." >> $log
					endif
				end
			end
# Extract motion and tissue signal
			echo " "
			echo "==================== 3dDeconvolve Nuisance Regression ==================="
			foreach r ($c_runs)
				echo " "
				echo "Current run: $r"
				echo " "
				foreach c ($censor)
					foreach h (lh rh)
						if (! -f ${o}${s}_$r.pb05.$h.NRerrts.${c}mm.niml.dset) then
							3dDeconvolve 									\
								-input ${o}${s}_$r.pb04.$h.scale.niml.dset 				\
								-polort A -jobs 10							\
								-ortvec ${o}${s}_$r.band.1D bandpass					\
								-ortvec ${o}${s}_$r.MoPar_derv.1D MoPar_derv 				\
								-ortvec ${o}${s}_$r.MoPar_demean.1D MoPar_demean 			\
								-ortvec ${o}${s}_${r}_mWM.1D mWM					\
								-ortvec ${o}${s}_${r}_mLV.1D mLV					\
								-censor ${o}${s}_${r}.${c}mm_censor.1D 					\
								-tout -fout -x1D ${o}${s}_$r.$h.xmat.1D					\
								-x1D_uncensored ${o}${s}_$r.$h.nocen.xmat.1D 				\
								-errts ${o}${s}_$r.$h.NRerrts.${c}mm.niml.dset				\
								-x1D_stop 								\
								-bucket ${o}${s}_$r.$h.NRstats.${c}mm.niml.dset
							3dTproject 									\
								-polort 0 								\
								-input ${o}${s}_$r.pb04.$h.scale.niml.dset 				\
								-censor ${o}${s}_${r}.${c}mm_censor.1D -cenmode ZERO 			\
								-ort ${o}${s}_$r.$h.nocen.xmat.1D 					\
								-prefix ${o}${s}_$r.pb05.$h.NRerrts.${c}mm.niml.dset
							if (! -f ${o}${s}_$r.pb05.$h.NRerrts.${c}mm.niml.dset) then
								echo "nuisance regression, ${r}:${c}mm:${h} failed..." >> $log
							endif
						else
							echo " "
							echo "${s}_$r.pb05.$h.NRerrts.${c}mm already exists"
							echo " "
						endif
					end
				end
			end
# Spatial smoothing
			echo " "
			echo "=========================== Smooth on Surface ==========================="
			foreach r ($c_runs)
				echo " "
				echo "Current run: $r"
				echo " "
				foreach c ($censor)
					foreach h (lh rh)
						if (! -f ${s}_$r.pb06.$h.smooth.${c}mm.niml.dset) then
							SurfSmooth 									\
								-spec $SUMA_dir/$s/SUMA/std.60.${s}_$h.spec 				\
								-surf_A smoothwm							\
								-met HEAT_07 								\
								-target_fwhm $fwhm 							\
								-detrend_master 							\
								-blurmaster ${o}${s}_$r.pb05.$h.NRerrts.${c}mm.niml.dset 		\
								-input ${o}${s}_$r.pb05.$h.NRerrts.${c}mm.niml.dset 			\
								-output ${s}_$r.pb06.$h.smooth.${c}mm.niml.dset
							if (! -f ${s}_$r.pb06.$h.smooth.${c}mm.niml.dset) then
								echo "smoothing, ${r}:${c}mm:${h} failed..." >> $log
							endif
						else
							echo " "
							echo "${s}_$r.pb06.$h.smooth.${c}mm already exists"
							echo " "
						endif
					end
				end
			end
# tsnr dataset
			if ($want_tsnr == yes) then
			echo " "
			echo "================================ Calculating TSNR ==============================="
			foreach r ($c_runs)
				echo " "
				echo "Current run: $r"
				echo " "
				if (! -f rm.$r.mask+orig.HEAD) then
					3dAutomask 										\
						-prefix rm.$r.mask 								\
						${o}${s}_$r.pb02.volreg+orig
					if (! -f rm.$r.mask+orig.HEAD) then
						echo "brain mask, run: $r failed..." >> $log
					endif

				else
					echo "${s}_$r brain mask not needed"
				endif
			end
				foreach r ($c_runs)
					3dTstat 										\
						-cvarinv -prefix ${s}_$r.TSNR 							\
						${o}${s}_$r.pb02.volreg+orig.BRIK
					3dROIstats 										\
						-nzmean -mask rm.$r.mask+orig 							\
						${s}_$r.TSNR+orig.BRIK > TSNR.$r.out
					if (! -f ${s}_$r.TSNR+orig.HEAD) then
						echo "TSNR, run: $r failed..." >> $log
					endif
				end
			endif
# Print results of script to log
# g = good runs, b = bad runs
			set g_runs = ()
			set b_runs = ()
			foreach r ($c_runs)
				foreach c ($censor)
					if (-f ${s}_$r.pb06.lh.smooth.${c}mm.niml.dset && -f ${s}_$r.pb06.rh.smooth.${c}mm.niml.dset) then
						set g_runs = ($g_runs $r.${c}mm)
					else
						set b_runs = ($b_runs $r.${c}mm)
					endif
				end
			end
			if ($#b_runs != 0) echo "FAILED:" $b_runs >> $log
			if ($#g_runs != 0) then
				echo "Processed:" $g_runs >> $log
			else
				echo "Processed: none" >> $log
			endif
# Clean up
			echo " "
			echo "Cleaning up..."
			echo " "
			\rm -f rm.*
			if ($clean_output == yes) rm -f *T1_* *aseg* *A2E*
		else
			echo " "
			echo "No runs to process, DOH!"
			echo "Processed: none" >> $log
		endif
	endif
	cd $data_dir
end
unset nonomatch
echo " "
echo "========================================================================"
echo "Consider your resting state data preprocessed, may the force be with you"
echo "========================================================================"
echo " "
exit