#!/bin/tcsh -f

if ("$1" == "-clean") then

	rm -f mask_SETR2.nii.gz
	rm -f SETR2RegParam_demean.1D
	rm -f mean_SETR2RegRCR.nii.gz
	rm -f SETR2RegRCRpc.nii.gz
	rm -f pos_beta_ds_SETR2RegRCRpc.nii.gz
	rm -f neg_beta_ds_SETR2RegRCRpc.nii.gz
	rm -f pos_ATSnum_ds_SETR2RegRCRpc.1D
	rm -f neg_ATSnum_ds_SETR2RegRCRpc.1D
	goto END

endif

if ( ! -f SETR2RegRCR.nii.gz || \
     ! -f SETR2RegParam.1D ) then
     echo "Missing DEMO files (SETR2RegRCR.nii.gz and SETR2RegParam.1D)."
     echo "They should have been in this archive. If you did not delete them"
     echo "by mistake, post a complaint on the AFNI"
     echo "message board"
     exit 1
endif

# create brain mask with 3dAutomask
if ( ! -f mask_SETR2.nii.gz) then 
	echo "Creating brain mask of SETR2RegRCR.nii.gz with 3dAutomask"
	3dAutomask -prefix mask_SETR2.nii.gz SETR2RegRCR.nii.gz
else
	echo "Reusing mask_SETR2.nii.gz"
endif

# demean realignment parameters
if ( ! -f SETR2RegParam_demean.1D) then
	echo "Removing mean of realignement parameters SETR2RegParam.1D"
	1d_tool.py -infile SETR2RegParam.1D -demean -write SETR2RegParam_demean.1D
else
	echo "Reusing SETR2RegParam_demean.1D"
endif

# normalize signal to signal percentage
if ( ! -f SETR2RegParam_demean.1D) then
	echo "Computing voxelwise mean volume mean_SETR2RegRCR.nii.gz"
	3dTstat -mean -prefix mean_SETR2RegRCR.nii.gz SETR2RegRCR.nii.gz
else
	echo "Reusing mean_SETR2RegRCR.nii.gz"
endif

if ( ! -f SETR2RegRCRpc.nii.gz) then
	echo "Computing percentage signal change of SETR2RegRCR.nii.gz"
	3dcalc -a SETR2RegRCR.nii.gz -b mean_SETR2RegRCR.nii.gz -expr '(a-b)/b' -prefix SETR2RegRCRpc.nii.gz
else
	echo "Reusing SETR2RegRCRpc.nii.gz"
endif

echo "NOW, 3dPFM would run but we will skip this step to save time."
echo "beta_ds_SETR2RegRCRpc.nii.gz and Z_Tstats_beta_ds_SETR2RegRCRpc.nii.gz are already computed."
echo "If you still want to run 3dPFM, execute the following command:"
# run 3dPFM
echo "3dPFM -input SETR2RegRCRpc.nii.gz -mask mask_SETR2.nii.gz -algorithm dantzig \ " 
echo "		-criteria bic -LHS SETR2RegParam_demean.1D -hrf SPMG1 -jobs 14 -outALL ds "

# This example considers the SPM canonical HRF for the deconvolution and the (demeaned) realignment parameters as LHS regressors. 
# The deconvolution is done with the dantzig selector algorithm where the regularization parameter is set according to the 
# Bayesian Information Criterion (BIC) in each voxel. Furthermore, the program will save all the possible output volumes 
# and use parallel computation using the snow R-package in 14 nodes (i.e. remove this option or change number of jobs to match your needs). 
# Computations will only be performed within non-zero voxels in mask mask_SETR2.nii.gz

# In this case, 3dPFM will creates 3D+time volumes with the results of the deconvolution, as well as the corresponding statistics 
# (see 3dPFM -help for more information). In order to explore the results, we recommend to see movie with the results of the 
# deconvolution as well as to create Activation Time Series that summarize the results in time (see below). 

# To see a movie with the results of 3dPFM, we recommend to proceed as follows. Set the UNDERLAY volume to SETR2RegRCR.nii.gz 
# and, for instance, set the OVERLAY volume to beta_ds_SETR2RegRCR.nii.gz. Open a graph to see the time series. 
# Set the AFNI environment variable AFNI_SLAVE_FUNCTIME = YES, and AFNI_SLAVE_THROLAY=YES This will change the overlay to 
# the same sub-brick index as the functional dataset. Set the Overlay threshold to 0 (or any relevant value 
# if you want to set a threshold in the amplitude of the beta coefficients to be seen in overlay). In the graph time series, 
# press the left/right-arrow keys to move the cursor. Enjoy the movie. 

echo "Setting environment variables AFNI_SLAVE_FUNCTIME YES and AFNI_SLAVE_THROLAY YES"
# In tcsh environment
setenv AFNI_SLAVE_FUNCTIME YES
setenv AFNI_SLAVE_THROLAY YES

echo "Running AFNI to see results of 3dPFM"
afni \
-com "SWITCH_UNDERLAY A SETR2RegRCR.nii.gz" \
-com "SWITCH_OVERLAY A beta_ds_SETR2RegRCRpc.nii.gz" \
-com 'OPEN_WINDOW A.axialimage geom=800x800 mont=5x4:1 ifrac=1' \
-com 'OPEN_WINDOW A.sagittalimage geom=600x600 mont=7x6:1 ifrac=1' \
-com 'SET_PBAR_ALL A.+99 1.0 Spectrum:red_to_blue+gap' \
-com 'SET_THRESHNEW A 0' \
-com 'SET_IJK A 35 67 12' \
-com 'SET_VIEW A.orig' \
-com 'SET_XHAIRS A.SINGLE' \
-com 'OPEN_WINDOW A.axialgraph -keypress=v'

echo "Computing Activation Time Series, counting the number of voxels"
echo "with positive or negative beta coefficients at each time point"
if ( ! -f pos_beta_ds_SETR2RegRCRpc.nii.gz) then
	echo "Computing volume of positive events"
	3dcalc -a beta_ds_SETR2RegRCRpc.nii.gz -expr 'ispositive(a)' -prefix pos_beta_ds_SETR2RegRCRpc.nii.gz
else
	echo "Reusing pos_beta_ds_SETR2RegRCRpc.nii.gz"
endif

if ( ! -f pos_ATSnum_ds_SETR2RegRCRpc.1D) then
	echo "Computing positive Activation Time Series"
	3dmaskave -q -mask mask_SETR2.nii.gz -sum pos_beta_ds_SETR2RegRCRpc.nii.gz > pos_ATSnum_ds_SETR2RegRCRpc.1D
else
	echo "Reusing pos_ATSnum_ds_SETR2RegRCRpc.1D"
endif

if ( ! -f neg_beta_ds_SETR2RegRCRpc.nii.gz) then
	echo "Computing volume of negative events"
	3dcalc -a beta_ds_SETR2RegRCRpc.nii.gz -expr 'isnegative(a)' -prefix neg_beta_ds_SETR2RegRCRpc.nii.gz
else
	echo "Reusing neg_beta_ds_SETR2RegRCRpc.nii.gz"
endif

if ( ! -f neg_ATSnum_ds_SETR2RegRCRpc.1D) then
	echo "Computing negative Activation Time Series"
	3dmaskave -q -mask mask_SETR2.nii.gz -sum neg_beta_ds_SETR2RegRCRpc.nii.gz > neg_ATSnum_ds_SETR2RegRCRpc.1D
else
	echo "Reusing neg_ATSnum_ds_SETR2RegRCRpc.1D"
endif	
	

echo "Plotting positive and negative Activation Time Series (with different scales)"
1dplot -sepscl pos_ATSnum_ds_SETR2RegRCRpc.1D neg_ATSnum_ds_SETR2RegRCRpc.1D

END: