eXTReMe Tracker
Oct 292012
 

As it turns out the if statement in bash isn’t as friendly towards floating point variables as it is for integer and string variables.

First of all the -gt and -lt switches don’t work.Even the unary operators >,<, etc. don’t work. The only way is to do a comparison between the floating point variables using bc and using the logical outout as a comparison string.

a=11.5; b=11.3
compare_result=`echo "$a >$b" | bc`
 if [ $compare_result -gt 0 ]; then
 echo $a
else
echo $b
 fi

 

Oct 242012
 

1. Convert the brain.mgz volume to analyze:

mri_convert $SUBJECTS_DIR/bert/mri/brain.mgz brain.img

2. Align brain.img with your functional data in spm (eg, func.img)

3. Create your spm statistical output (eg, t.img)

4. Create the registration file:

tkregister2 --s subjectname --mov meanswaf.img --regheader --reg  register.dat --surf orig

This will bring up the tkregister window with the orig volume. Hit the “Compare” button to see the functional. The green line will be the surface. Make sure the alignment is good (ie, the green line follows the bright intensity patterns in the functional). Hit the “Save” button to save the registration. This will create a file called register.dat. If you make modifications to the registration and then want to view or edit it later, re-run the above command WITHOUT –regheader. This requires a lot of manual tweaking to get the functional brain aligned with Bert brain.

5. View your functional data in the volume:

tkmedit subjectname orig.mgz -overlay spmT_0001.img -overlay-reg register.dat

6. View your functional data on the surface:

tksurfer subjectname lh inflated -overlay spmT_0001.img -overlay-reg register.dat

Adopted from spmPainting

Oct 012012
 

Batch analysis

#!/bin/bash

# run a batch analysis on one subject
#
# for info about installation, usage, and dependencies see doc/readme.txt
# or execute with no parameters for a usage message
#
# mjp rul3zzzzz
# 2008-01-25

#### globals ####

exec_dir=`echo $0 | sed "s|\(.*/\).*|\1|"`

# name of the config file to read params from
config_filename="UNSET"

# where all the analysis specific files are (build_model.m, etc)
analysis_dir="UNSET"

# the top level directory containing all the subject directories for this study
study_subjs_dir="UNSET"

# path to batch_analysis_analysis dir
batch_analysis_home=/software/batch_analysis

# analysis unspecific matlab files dir (spm jobman templates, etc.)
matlab_dir="$batch_analysis_home"/matlab
matlab_dir_override=0

# analysis unspecific shell scripts dir (extract_motion_parms.sh, etc.)
script_dir="$batch_analysis_home"/scripts
script_dir_override=0

# path to customized art
art_dir="$batch_analysis_home"/matlab/art
art_dir_override=0

# analysis package homes
spm_dir=/software/spm5_1782
fsl_dir=/etc/fsl/fsl.sh

# freesurfer path
freesurfer_home=$FREESURFER_HOME

# stages to start and stop analysis at if partial analysis desiredx
startstage="first"
stopstage="last"

# individual stages 

# be careful, just because you set something to 1 doesnt mean that the
# result is passed to the next stage. for example, setting only
# do_reslice=1 doesnt set do_normalize=0, and by default normalized
# ("w" prefix) files are passed along to the next analysis stage

do_slice_timing=0
do_fieldmap_correct=0
do_motion_correct=1
do_normalize=1
do_surf_project=0
do_outliers=1
do_view_funcruns=1
do_prepare_analysis=1
do_run_firstlevel=1
do_run_contrasts=1
do_surf_paint=0

# fieldmap correct options
fieldmap_defaults_file=
fieldmap_vdm_file=
fieldmap_total_epi_readout_time=

# motion correction options
do_reslice=0
use_moco_for_motion_estimate=0

# smoothing options
do_vol_smooth=1
do_surf_vol_smooth=0
do_surf_smooth=0

fwhm=

# outlier options
skip_art_if_outliers_file_exists=0 # if outliers exist dont overwrite
same_outliers_as_base_model=1 # whether the current model has same outliers as
                              # base (same_runs_as_base_model better be 1, too!)

# estimation options
TR=2
do_vol_estimate=1
do_surf_estimate=0

design_reporting=1
delete_old_analysis_without_asking=0 # delete old analyses automatically

same_runs_as_base_model=1 # whether the current model has same runnums as base
fir_model=0
do_bayes_estimation=0 # use the spm baysian estimation method 
no_design_matrix_temporal_filter=0 # skip temporal filtering of design matrix
no_ar_1_correct=0 # skip temporal correlation correction
use_temporal_deriv=1   # set to 0 in set_parm if do_slice_timing is set to 1
fsl_estimate=0 # whether to use fsl estimation in place of spm estimation
est_model_dir=            # look in a different dir than default for betas

# surf estimate options
dont_overwrite_resliced_big_niftis=0 # if sliced niftis exist dont overwrite

# surf paint options
paint_beta=0
paint_con=1
paint_t=1

## dont set these
prefix= # accumulated depending on processing steps

#### functions ####

# let em know how it is!
# TODO: add full descriptions of options (complete descriptions of -c and -m)
usage() {
    echo "$0 [options] subj"
    echo " options:"
    echo "   -c <config_filename>"
    echo "      name/value pairs, one to a line. any variable can be read from a config"
    echo "      file. for supported variables, see doc/example.cfg"

    echo "   -d <analysis_dir> "
    echo "      path to directory containing analysis specific matlab scripts"
    echo "       (build_model.m, build_contrasts.m, etc)"

    echo "   -D <study_subjs_dir> "
    echo "      path to directory containing the individual subject directories"

    echo "   -I <dicom_dir> "
    echo "      path to directory containing the dicom files for this subject"
    echo "      if this does not exist the most recent directory starting "
    echo "      with 'TrioTim' in the subject directory will be assumed "

    echo "   -b <batch_analysis_home> "
    echo "      path to directory containing the batch analysis distribution"

    echo "   -s <stage>: start at stage <stage>"
    echo "    stage can be:   "
    echo "       extract_motion"
    echo "       extract_volumes"
    echo "       slice_timing"
    echo "       fieldmap_correct"
    echo "       motion_correct"
    echo "       normalize"
    echo "       surf_project"
    echo "       smooth"
    echo "       outliers"
    echo "       view_funcruns"
    echo "       prepare_analysis"
    echo "       run_firstlevel"
    echo "       run_contrasts"
    echo "       surf_paint"

    echo "   -S <stage>: stop after stage <stage> (same options as start stage)"
    echo "   -R run numbers are different for this model than base model"
    echo "   -m <model>: create model_<model> instead of model (see doc/readme.txt)"
    echo "   -e [<estimated_model>]: use model_<estimated_model> to compute contrasts"
    echo "   -h, -?: print usage"
}

if [ ! "$1" ]; then
    usage
    exit
fi

# sets a parameter to a specified value
#
# in: paramter name and value pair
set_parm() {
    name="$1"
    value="$2"

    echo setting $name=\"$value\"

    # make the parm available
    export $name="$value"

    # check for special handling
    if [ "$name" == "batch_analysis_home" ]; then
	set_batch_analysis_home "$2"
    fi
    if [ "$name" == "matlab_dir" ]; then
	echo setting matlab_dir_override=1
	matlab_dir_override=1
    fi
    if [ "$name" == "script_dir" ]; then
	echo setting script_dir_override=1
	script_dir_override=1
    fi
    if [ "$name" == "art_dir" ]; then
	echo setting art_dir_override=1
	art_dir_override=1
    fi
    if [ "$name" == "est_model_dir" ]; then
	echo setting est_model_dir_override=1
	est_model_dir_override=1
    fi
    if [ "$name" == "slice_timing" -a "$value" == "1" ]; then
	echo "Using Slice Timing Correction; disabling temporal derivatives"
	use_temporal_deriv=0
    fi
}

# set the batch_analysis home
# updates the matlab, scripts, and art dir if they are not overriden.
#
# in: path to the batch_analysis home directory
set_batch_analysis_home() {
    batch_analysis_home=$1

    if [ $matlab_dir_override == "0" ]; then
	matlab_dir="$batch_analysis_home"/matlab
    fi

    if [ $script_dir_override == "0" ]; then
	script_dir="$batch_analysis_home"/scripts
    fi

    if [ $art_dir_override == "0" ]; then
	art_dir="$batch_analysis_home"/matlab/art
    fi
}

# parse command line arguments
#
# in: $@
parse_args() {
   # get the options
    while getopts ":c:d:D:e:h:s:S:m:RI:h?" Option
    do
	case $Option in
	    c ) set_parm config_filename "$OPTARG"
		;;
	    d ) set_parm analysis_dir "$OPTARG"
		;;
	    D ) set_parm study_subjs_dir "$OPTARG"
		;;
	    I ) set_parm dicom_dir "$OPTARG"
		;;
	    b ) set_parm batch_analysis_home "$OPTARG"
		;;
	    s ) set_parm startstage "$OPTARG"
		;;
	    S ) set_parm stopstage "$OPTARG"
		;;
	    m ) set_parm model "$OPTARG"
		;;
	    e ) set_parm est_model_dir "$OPTARG"
		;;
	    R ) set_parm same_runs_as_base_model 0;;
	    * ) usage;
		exit;;
	esac
    done
    shift $(($OPTIND - 1))

    subj=$1
}

# read parameters from a config file
# TODO: handle comments better
#
# in: configuration filename
read_configfile() {
    filename=$1;

    set -- `cat $filename | sed "/#/d; /^\s*$/d;"`
    while [ "$1" ]; do
	varval=`env | grep "^$1="`

	if [ ! "$varval" ]; then
	    set_parm $1 $2
	else
	    echo warning: not overidding $varval
	fi
	shift 2
    done
}

# validate configuration, report errors
# checks that required directories have been set
# TODO: check that the required files are in the directories, too
validate_config() {
    # analysis directory
    if [ "$analysis_dir" == "UNSET" ]; then
	echo "ERROR! analysis directory has not been set"
	usage;
	exit 1;
    fi

    if [ ! -d "$analysis_dir" ]; then
	echo "ERROR! invalid analysis directory: $analysis_dir"
	usage;
	exit 1;
    fi

    # study subjects directory
    if [ "$study_subjs_dir" == "UNSET" ]; then
	echo "ERROR! study_subjs directory has not been set"
	usage;
	exit 1;
    fi

    if [ ! -d "$study_subjs_dir" ]; then
	echo "ERROR! invalid study_subjs directory: $study_subjs_dir"
	usage;
	exit 1;
    fi

    # subject name
    if [ "$subj" == "UNSET" ]; then
	echo "ERROR! subect name is unset"
	usage;
	exit 1;
    fi

    # check paths
    if [ ! `which fslmaths` ]; then
	source "$fsl_dir/fsl.sh"
    fi

    if [ ! `which mri_info` ]; then
	export SUBJECTS_DIR=$PWD/..
	source "$freesurfer_home/SetUpFreeSurfer.sh"
    fi
}

# dump the configuration to stdout
# TODO: do this for only our parms
dump_config() {
    env
}

# determine if a path is absolute
path_is_absolute() {
    if [ `echo $1 | sed "s/^\(.\).*/\1/"` == "/" ]; then
	return 1
    else
	return 0
    fi
}

# get the stem of the nifti files
get_niistem() {
    export niistem=`ls nii/[0-9]*.nii  | tail -n 1 | sed "s/-.*.nii/-/" | sed "s|nii/||"`
}

# map stage names to numbers
#
# in: name of desired stage number
get_stage_number() {
    name=$1

    case $name in
	first                        ) stage_num=1 ;;
	extract_motion               ) stage_num=1 ;;
	extract_volumes              ) stage_num=2 ;;
	slice_timing                 ) stage_num=3 ;;
	fieldmap_correct             ) stage_num=4 ;;
	motion_correct               ) stage_num=5 ;;
	normalize                    ) stage_num=6 ;;
	surf_project                 ) stage_num=7 ;;
	smooth                       ) stage_num=8 ;;
	mask                         ) stage_num=9 ;;
	outliers                     ) stage_num=10;;
	view_funcruns                ) stage_num=11;;
	prepare_analysis             ) stage_num=12;;
	run_firstlevel               ) stage_num=13;;
	run_contrasts                ) stage_num=14;;
	surf_paint                   ) stage_num=15;;
	last                         ) stage_num=15;;
	*                            ) stage_num=0;;
    esac
}

# map stage names to volume prefix addition
#
# in: name of desired stage number
update_prefix() {
    name=$1

    case $name in
	slice_timing     )
	    if [ "$do_slice_timing" == 1 ] ; then
		prefix=a$prefix
	    fi
	    ;;

	fieldmap_correct     )
	    if [ "$do_fieldmap_correct" == 1 ] ; then
		prefix=u$prefix
	    fi
	    ;;

	motion_correct   )
	    if [ "$do_reslice" == 1 -a "$do_normalize" == 0 ] ; then
		prefix=r$prefix
	    fi
	    ;;

	normalize        )
	    if [ "$do_normalize" == 1 ] ; then
		prefix=w$prefix
	    fi
	    ;;

	surf_project        )
	    if [ "$do_surf_project" == 1 ] ; then
		lhprefix=lh$prefix
		rhprefix=rh$prefix
	    fi
	    ;;

	smooth           )
	    if [ "$do_vol_smooth" == 1 ] ; then
		prefix=s$prefix
	    fi
	    if [ "$do_vol_surf_smooth" == 1 ] ; then
		prefix=S$prefix
	    fi
	    if [ "$do_surf_smooth" == 1 ] ; then
		lhprefix=S$lhprefix
		rhprefix=S$rhprefix
	    fi
	    ;;
	* ) ;;
    esac

    export prefix=$prefix
}

# get the motion from the dicom headers
#
# in:
#     subject directory name
#     dicom directory name
extract_motion() {
    subjdir="$1";
    dcmdir="$2";

    # estimate motion
    echo "extracting motion"
    mkdir -p "$subjdir"/motion
    cd "$dcmdir";
    "$script_dir"/extract_motion_parms.sh -S -t spm -d ../motion/
    cd ..

}

# convert the dicom directory into nifti volumes
#
# in: dicom directory name
extract_volumes() {
    dcmdir="$1"

    # extract nii files
    echo "extracting nii files from $dcmdir"
    mkdir -p nii
    "$script_dir"/dicomdir2nii.sh "$dcmdir" nii

}

# print available volumes and how many frames they have
print_avail_vols() {
    echo nothing
}

# slice-timing correction, if requested. (default is no)
slice_timing() {
    runnums="$1"

	#create matlab script to do slice timing correction via spm
    slice_timing_script=scripts/slice_timing.m
    echo "% auto-generated by $0" > $slice_timing_script
    echo "load "$matlab_dir"/slice_timing_job.mat" >> $slice_timing_script
    echo "addpath('scripts');" >> $slice_timing_script

    # build file list
    get_niistem
    echo "jobs{1}.temporal{1}.st.scans = {" >> $slice_timing_script
    for i in $runnums; do
	echo "'$PWD/nii/$niistem$i.nii'," >> $slice_timing_script
    done
    echo "};" >> $slice_timing_script

    if [ "$slice_timing_ref_slice" ]; then
	echo "jobs{1}.temporal{1}.st.refslice = $slice_timing_ref_slice;" >> $slice_timing_script
    fi

    if [ "$slice_timing_slice_order" ]; then
	echo "jobs{1}.temporal{1}.st.so = [$slice_timing_slice_order];" >> $slice_timing_script
    fi

    echo "spm_jobman('run',jobs);" >> $slice_timing_script

    "$script_dir"/run_matlab_script.sh "$slice_timing_script"
}

# set some filenames for fieldmap correcting
set_fieldmap_params() {

    get_niistem

    # paths to the images (and echo times if necessary)
    if [ ! "$fieldmap_mag_img" -a ! "$fieldmap_phs_img" ]; then
	# lotsa assumptions here ...
	for img in `ls $dicom_dir/*-1.dcm`; do
	    isfieldmap=`dicom_hdr $img | grep "ACQ Protocol Name" | grep "field_mapping"`
	    if [ "$isfieldmap" ]; then
		# check for mag img set (assume its first)
		if [ ! "$fieldmap_mag_img" ]; then
		    echo "found magnitude image: $img"

		    num=`echo $img | sed "s/.*[0-9]\+-\([0-9]\+\)-1.dcm/\1/"`
		    export fieldmap_mag_img=$PWD/nii/$niistem$num.nii;
		elif [ ! "$fieldmap_phs_img" ]; then
		    echo "found phase image: $img"

		    num=`echo $img | sed "s/.*[0-9]\+-\([0-9]\+\)-1.dcm/\1/"`
		    export fieldmap_phs_img=$PWD/nii/$niistem$num.nii;

		else
		    echo "ignoring extra fieldmap $img"
		fi
	    fi
	done
    fi

    # set echo times
    if [ ! "$fieldmap_short_echo" ]; then
	num=`echo $fieldmap_mag_img | sed "s/.*[0-9]\+-\([0-9]\+\).nii/\1/"`
	export fieldmap_short_echo=`dicom_hdr $dicom_dir/*-$num-1.dcm | grep "Echo Time" | sed "s|.*Echo Time//||"`
	echo "automatically set fieldmap_short_echo to $fieldmap_short_echo"
    fi

    if [ ! "$fieldmap_long_echo" ]; then
	num=`echo $fieldmap_phs_img | sed "s/.*[0-9]\+-\([0-9]\+\).nii/\1/"`
	export fieldmap_long_echo=`dicom_hdr $dicom_dir/*-$num-1.dcm | grep "Echo Time" | sed "s|.*Echo Time//||"`
	echo "automatically set fieldmap_long_echo to $fieldmap_long_echo"
    fi

    # set vdm filename
    if [ ! "$fieldmap_vdm_filename" ]; then
	num=`echo $fieldmap_phs_img | sed "s/.*[0-9]\+-\([0-9]\+\).nii/\1/"`
	export fieldmap_vdm_filename=$PWD/nii/vdm5_sc$niistem$num.nii
	echo "automatically set fieldmap_vdm_filename to $fieldmap_vdm_filename"
    fi

}

# spm fieldmap correction
#
# in: first functional run number
fieldmap_correct() {
    runnum="$1";
    get_niistem

    set_fieldmap_params

    # target epi image
    if [ ! "$fieldmap_epi_img" ]; then
	fieldmap_epi_img=$PWD/nii/$niistem$runnum.nii
    fi

    # make paths absolute
    path_is_absolute $fieldmap_mag_img
    if [ "$?" -ne 1 ]; then
	fieldmap_mag_img=$PWD/nii/$fieldmap_mag_img
    fi

    path_is_absolute $fieldmap_phs_img
    if [ "$?" -ne 1 ]; then
	fieldmap_phs_img=$PWD/nii/$fieldmap_phs_img
    fi

    path_is_absolute $fieldmap_epi_img
    if [ "$?" -ne 1 ]; then
	fieldmap_epi_img=$PWD/nii/$fieldmap_epi_img
    fi

    # error check
    if [ ! -f "$fieldmap_mag_img" -o ! -f "$fieldmap_phs_img" ]; then
	echo "magnitude or phase images not found! can't fieldmap correct"
	echo "fieldmap_mag_img: $fieldmap_mag_img"
	echo "fieldmap_phs_img: $fieldmap_phs_img"
	exit 0
    fi

    # readout time

    ## eventually
    #bandwidth=`dicom_hdr $dicom_dir/*$runnum-1.dcm | sed "s|.*Bandwidth//||"`;
    if [ ! "$fieldmap_total_epi_readout_time" ]; then
	echo "readout time not specified! can't fieldmap correct"
	return
    fi

    # defaults file
    if [ ! "$fieldmap_defaults_file" ]; then
	fieldmap_defaults_file="$matlab_dir/fieldmap_defaults.m"

	if [ ! "$fieldmap_short_echo" -o ! "$fieldmap_long_echo" ]; then
	    echo "both short and long echo times must be set! can't fieldmap correct"
	    return
	fi
    fi

    echo "$script_dir"/fieldmap_correct.sh \
	-m "$fieldmap_mag_img" \
	-p "$fieldmap_phs_img" \
	-e "$fieldmap_epi_img" \
	-f "$fieldmap_defaults_file" \
	-s "$fieldmap_short_echo" \
	-l "$fieldmap_long_echo" \
	-t "$fieldmap_total_epi_readout_time" \
	-d "$study_subjs_dir" \
	-D "$matlab_dir" \
	$subj

    "$script_dir"/fieldmap_correct.sh \
	-m "$fieldmap_mag_img" \
	-p "$fieldmap_phs_img" \
	-e "$fieldmap_epi_img" \
	-f "$fieldmap_defaults_file" \
	-s "$fieldmap_short_echo" \
	-l "$fieldmap_long_echo" \
	-t "$fieldmap_total_epi_readout_time" \
	-d "$study_subjs_dir" \
	-D "$matlab_dir" \
	$subj

}

# spm realign the functional runs
#
# in: functional run numbers
motion_correct() {
    runnums="$1";

    # set flags and prefix
    if [ "$do_reslice" == 1 ]; then
	reslice_flag=-R
    fi

    if [ "$do_fieldmap_correct" == 1 ]; then
	set_fieldmap_params

	reslice_flag=-R
	unwarp_flag="-U -V $fieldmap_vdm_filename"
    fi

    echo "$script_dir"/motion_correct.sh \
	-p "$prefix" \
	-r "$runnums" \
	-d "$study_subjs_dir" \
	$reslice_flag \
	$unwarp_flag \
	-D "$matlab_dir" \
	$subj

    "$script_dir"/motion_correct.sh \
	-p "$prefix" \
	-r "$runnums" \
	-d "$study_subjs_dir" \
	$reslice_flag \
	$unwarp_flag \
	-D "$matlab_dir" \
	$subj
}

# spm sptatial normalization
#
# in: functional run numbers
normalize() {
    runnums="$1";

    echo "$script_dir"/normalize_group_of_subjects.sh \
	 -t "$spm_dir"/templates/EPI.nii \
	 -d "$study_subjs_dir" \
	 -D "$matlab_dir" \
	 -B -C \
	 -r \"$1\" \
	 $subj
     "$script_dir"/normalize_group_of_subjects.sh \
	 -t "$spm_dir"/templates/EPI.nii \
	 -d "$study_subjs_dir" \
	 -D "$matlab_dir" \
	 -B -C \
	 -r "$1" \
	 $subj
}

# spm spatial smoothing
#
# in: list of functional run numbers
vol_smooth() {
    runnums="$1";

    if [ "$fwhm" ]; then
	fwhm_flag="-f $fwhm"
    fi

    echo "$script_dir"/volume_smooth.sh \
	-r "$runnums" \
	-p $prefix \
	-d "$study_subjs_dir" \
	-D "$matlab_dir" \
	$fwhm_flag \
	$subj

    "$script_dir"/volume_smooth.sh \
	-r "$runnums" \
	-p $prefix \
	-d "$study_subjs_dir" \
	-D "$matlab_dir" \
	$fwhm_flag \
	$subj
}

# unzip a file if it exists
#
# in
#  base filename 
gunzip_file() {
    if [ -e "$1.gz" ]; then
	echo gunzip -f "$1.gz"
	gunzip -f "$1.gz"
    fi
}

# make a mask of all ones same size as template volume
#
# in
#  $1 template volume
#  $2 filename
create_mask_of_ones() {
    echo fslmaths "$1" -Tmean -mul 0 -add 1 "$2";
    fslmaths "$1" -Tmean -mul 0 -add 1 "$2";
    gunzip_file "$2"
}

# register mean functional volume to surface anatomical
#
#
reg_subject() {
    meanfile=`ls nii/mean*.nii | sed "s|nii/||"`
    if [ ! "$meanfile" ]; then
	echo error, could not find meanfile in nii directory
	return 1
    fi

    meanname=`echo $meanfile | sed "s|.*/||; s|.nii||"`
    regfile=mri/transforms/"$meanname"_to_struc.xfm
    if [ ! -e "$regfile" ]; then
	echo "$script_dir"/register_func_surf.sh \
	    -s $subj \
	    -f $PWD/nii/$meanfile \
	    -d "$study_subjs_dir" \
	    -r "$regfile"

	"$script_dir"/register_func_surf.sh \
	    -s $subj \
	    -f $PWD/nii/$meanfile \
	    -d "$study_subjs_dir" \
	    -r "$regfile"
    fi

    export regfile=$regfile
    export meanfile=$meanfile
}

# spatial smoothing on the surface but keep in the volume
#
# in: list of functional run numbers
surf_vol_smooth() {
    runnums="$1";

    if [ ! "$regfile" ]; then
	reg_subject
    fi

    echo "$script_dir"/smooth_vol_on_surf.sh \
	-r "$runnums" \
	-R "$regfile" \
	-P $prefix \
	-d "$study_subjs_dir" \
	$subj
    "$script_dir"/smooth_vol_on_surf.sh \
	-r "$runnums" \
	-R "$regfile" \
	-P $prefix \
	-d "$study_subjs_dir" \
	$subj
}

# surface projection
#
# in: list of functional run numbers
surf_project() {
    runnums="$1"

    if [ ! "$regfile" ]; then
	reg_subject
    fi

   echo "$script_dir"/project_to_surface.sh \
       -r "$runnums" \
       -d "$study_subjs_dir" \
       -p $prefix \
       -P 0.5  \
       -R "$regfile" \
       $subj

   "$script_dir"/project_to_surface.sh \
       -r "$runnums" \
       -d "$study_subjs_dir" \
       -p $prefix \
       -P 0.5  \
       -R "$regfile" \
       $subj
}

# smooth data on the surface
#
# in: list of functional run numbers
surf_smooth() {
    runnums="$1"

    echo smooth_surf_on_surf.sh \
	-r "$runnums" \
       -d "$study_subjs_dir" \
       -L $lhprefix \
       -R $rhprefix \
	$subj
    smooth_surf_on_surf.sh \
	-r "$runnums" \
       -d "$study_subjs_dir" \
       -L $lhprefix \
       -R $rhprefix \
	$subj
}

# split big niftis into multiple slices
#
# in:
#  $1 list of files
slice_big_nifti() {
    files="$1"

    echo slicing big niftis: $files

    # create a matlab script to take care of the big nifti issue
    bignii_script=scripts/bignii.m

    echo "% auto-generated by $0" > $bignii_script
    echo "addpath $matlab_dir" >> $bignii_script

    for i in $files; do
	nii=`echo $i | sed "s/fs.nii$/nii/"`

	if [ -e "$nii" -a "$dont_overwrite_resliced_big_niftis" == 1 ]; then
	    echo "skipping $nii because im not overwritting resliced niftis"
	else
	    echo "reshape_big_nifti('$i','$nii');" >> $bignii_script
	fi
    done
    "$script_dir"/run_matlab_script.sh "$bignii_script"
}

# unsplit niftis that were previously split into multiple slices
#
# in:
#  $1 number of surface vertices
#  $2 list of files
unslice_big_nifti() {
    numverts="$1"
    files="$2"

    echo unslicing to big niftis: $files

    # create a matlab script to un-take-care of the big nifti issue
    bignii_script=scripts/unbignii.m

    echo "% auto-generated by $0" > $bignii_script
    echo "addpath $matlab_dir" >> $bignii_script

    for i in $files; do
	fsnii=`echo $i | sed "s/img$/fs.nii/"`
	echo "unreshape_big_nifti('$i',$numverts,'$fsnii');" >> $bignii_script
    done
    "$script_dir"/run_matlab_script.sh "$bignii_script"
}

# run art to determine outliers
#
# in: functional run numbers
art_batch() {
    runnums="$1";

    # get image names for each session

    # build file list
    get_niistem

    if [ "$use_moco_for_motion_estimate" == 1 ]; then
	motstem="spm_motion-"
    else
	motstem="rp_$niistem"
    fi

    ind=0;
    sess_str=""
    for i in $runnums; do
	let ind="ind+1"

	if [ "$use_moco_for_motion_estimate" == 1 ]; then
	    let motind="i+1"
	else
	    let motind="i"
	fi

	sess_str="$sess_str session $ind image $niistem$i.nii"
	mot_str="$mot_str session $ind motion $motstem$motind.txt"
    done

    # make the art session file
    art_sess_file=scripts/art_sess_file.asf
    echo "sessions: $ind" > $art_sess_file
    echo "global_mean: 1" >> $art_sess_file
    echo "drop_flag: 0" >> $art_sess_file
    echo "motion_file_type: 0" >> $art_sess_file
    echo "image_dir: $PWD/nii" >> $art_sess_file

    if [ "$use_moco_for_motion_estimate" == 1 ]; then
	echo "motion_dir: $PWD/motion" >> $art_sess_file
    else
	echo "motion_dir: $PWD/nii" >> $art_sess_file
    fi

    echo "motion_fname_from_image_fname: 0" >> $art_sess_file
    echo "end" >> $art_sess_file
    echo ""
    echo $sess_str >>  $art_sess_file
    echo $mot_str >>  $art_sess_file
    echo "end" >>  $art_sess_file

    art_script=scripts/art_batch.m

    echo "% auto-generated by $0" > $art_script
    echo "addpath " $art_dir >> $art_script
    echo "a=which('art');" >> $art_script
    echo "fprintf('art is at %s\n',a);" >> $art_script
    echo "art('sess_file','$art_sess_file');" >> $art_script
    echo "fprintf('make sure to save the outliers in the mat file:\n$PWD/nii/$outlierfile\npress enter to continue...'), pause;" >> $art_script
    "$script_dir"/run_matlab_script.sh -with-display "$art_script"
}

# view each functional run in fslview
#
# in: functional run numbers
view_funcruns() {
    runnums="$1"

    get_niistem
    for i in $runnums; do
	fslview $PWD/nii/$niistem$i.nii;
    done
}

# create spm model
#
# in:
#  $1 functional run numbers
#  $2 model name or "" for none
#  $3 modelspace vol or surf
#  $4 hemisphere (if modelspace is surf)
build_model() {
    runnums="$1";
    modelname="$2";
    if [ "$3" == "surf" -o "$3" == "surface" ]; then
	modelspace=surf
	hemi="$4"
    else
	modelspace=vol
    fi

    modeldir="model"
    modelmfun="build_model"
    dir=$PWD/firstlevel;

    # create a matlab script
    model_script=scripts/model

    # look for a model type modifier
    if [ "$modelname" ]; then
	modeldir="$modeldir"_"$modelname"
	modelmfun="$modelmfun"_"$modelname"
	dir="$dir"_"$2"
	model_script="$model_script"_"$modelname"

	if [ "$same_outliers_as_base_model" -eq 1 ]; then
	    outlierfile=outliers.mat
	else
	    outlierfile=outliers_"$modelname".mat
	fi
    else
	outlierfile=outliers.mat
    fi

    # check if this is a surface analysis
    if [ "$modelspace" == "surf" ]; then
	dir="$dir"/$hemi
	model_script="$model_script"_$hemi
    fi

    model_script="$model_script".m

    mkdir -p "$modeldir"

    echo "addpath('"$analysis_dir"/');" > $model_script
    echo "dbstop if error;" >> $model_script

    # write nii file list based on model space
    echo "runs = {" >> $model_script
    niinum=`ls nii/[0-9]*.nii  | tail -n 1 | sed "s/-.*.nii/-/;s/nii.//"`
    if [ "$modelspace" == "surf" ]; then
	if [ "$hemi" == "lh" ]; then
	    hemiprefix=$lhprefix
	else
	    hemiprefix=$rhprefix
	fi

	niistem=`ls nii/$hemiprefix[0-9]*.nii  | tail -n 1 | sed "s/-.*.nii/-/"`

	explicit_mask=
	for i in $runnums; do
	    echo "'$PWD/nii/$niistem$i.nii'," >> $model_script

	    # build a mask of all ones same size as surface volume
	    if [ ! "$explicit_mask" ]; then
		explicit_mask=$PWD/nii/ones_mask_$hemi.nii
		create_mask_of_ones "$PWD/nii/$niistem$i.nii" "$explicit_mask"
	    fi

	    # determine if there are big niftis that need slicing
	    bn=`ls $PWD/nii/$hemiprefix*-$i.fs.nii`;
	    if [ "$bn"  ]; then
		slice_big_nifti "$bn"
	    fi

	done

    else # find volumes
	get_niistem
	niinum=`ls nii/[0-9]*.nii  | tail -n 1 | sed "s/-.*.nii/-/;s/nii.//"`
	for i in $runnums; do
	    echo "'$PWD/nii/$prefix$niistem$i.nii'," >> $model_script
	done
    fi

    echo "};" >> $model_script

    # build motionfile list
    if [ "$use_moco_for_motion_estimate" == 1 ]; then
	motstem="motion/spm_motion-"
    else
	motstem="nii/rp_$niinum"
    fi
    echo "motfiles = {" >> $model_script
    for i in $runnums; do
	if [ "$use_moco_for_motion_estimate" == 1 ]; then
	    let j="i+1"
	else
	    let j="i"
	fi
	echo "'$PWD/$motstem$j.txt'," >> $model_script
    done
    echo "};" >> $model_script

    # load outliers
    echo "" >> $model_script
    echo "out_idx = [];" >> $model_script
    echo "if(exist('$PWD/nii/$outlierfile','file'))" >> $model_script
    echo "  load('$PWD/nii/$outlierfile');" >> $model_script
    echo "end" >> $model_script

    # load model goby
    echo "load "$matlab_dir"/firstlevel_job.mat" >> $model_script

    # set tr
    echo "jobs{1}.stats{1}.fmri_spec.timing.RT = $TR;" >> $model_script

    # delete previous analysis if there and required
    if [ -e "$dir/SPM.mat" -a "$delete_old_analysis_without_asking" == 1 ]; then
	echo "WARNING: removing existing analysis in $dir!!!!"
	rm "$dir/SPM.mat"
    fi

    mkdir -p "$dir"
    echo "jobs{1}.stats{1}.fmri_spec(1).dir = {'$dir'};" >> $model_script

    # temporal derivatives
    if [ "$fir_model" == 1 ]; then
	echo "jobs{1}.stats{1}.fmri_spec(1).bases = struct();" >> $model_script
	echo "jobs{1}.stats{1}.fmri_spec(1).bases.fir.length = 20;" >> $model_script
	echo "jobs{1}.stats{1}.fmri_spec(1).bases.fir.order = 10;" >> $model_script
    elif [ "$use_temporal_deriv" == 1 ]; then
	echo "jobs{1}.stats{1}.fmri_spec(1).bases.hrf.derivs = [1 0];" >> $model_script
    else
	echo "jobs{1}.stats{1}.fmri_spec(1).bases.hrf.derivs = [0 0];" >> $model_script
    fi

    # build model
    echo " jobs{1}.stats{1}.fmri_spec(1).sess = $modelmfun('$subj',runs,motfiles,out_idx);"  >> $model_script

    # run
    echo "dbstop if error"  >> $model_script
    echo "spm_jobman('run',jobs);" >> $model_script
    if [ "$design_reporting" -eq "1" ]; then
	echo "pause;" >> $model_script
    fi

    # set mask (SPM sucks!!)
    if [ "$explicit_mask" ]; then
	echo "load $dir/SPM;" >> $model_script
	#echo "SPM.xM = [];" >> $model_script
	echo "SPM.xM.VM = spm_vol('$explicit_mask');" >> $model_script
	echo "SPM.xM.I = 0;" >> $model_script
	echo "SPM.xM.T = [];" >> $model_script
	echo "SPM.xM.TH = ones(size(SPM.xM.TH))*(-Inf);" >> $model_script
	echo "SPM.xM.xs = struct('Masking', 'explicit masking only');" >> $model_script
	echo "save $dir/SPM SPM;" >> $model_script
    fi

    # set default non-sphericity for surfaces to be diagonal (dirty)
    if [ "$no_ar_1_correct" -eq "1" ]; then
	echo "WARNING!!! not using spm built-in AR(1) correction!!!"
	echo "load $dir/SPM;" >> $model_script
	echo "SPM.xVi.Vi = {speye(size(SPM.xX.X,1))};" >> $model_script
	echo "SPM.xVi.V  = speye(size(SPM.xX.X,1));" >> $model_script
	echo "save $dir/SPM SPM;" >> $model_script
    fi

    "$script_dir"/run_matlab_script.sh "$model_script"

    # turn off temporal filtering of the design matrix (SPM sux!!)
    if [ "$no_design_matrix_temporal_filter" -eq "1" ]; then
	echo "WARNING!!! not temporally filtering the design matrix!!!"
	echo "load $dir/SPM;" >> $model_script
	echo "SPM.xY.K = SPM.xX.K;" >> $model_script
	echo "SPM.xX.K = 1;" >> $model_script
	echo "save $dir/SPM SPM;" >> $model_script
    fi
}

# estimate firstlevel model, just choose between fsl and spm
#
# in:
#  $1 model name or "" for none
#  $2 modelspace vol or surf
#  $3 hemisphere (if modelspace is surf)
run_firstlevel() {
    if [ "$fsl_estimate" -eq "1" ]; then
	# run spm first to trick it into thinking that the model has 
        # been estimated so that contrasts can be run
	echo "running dummy spm estimation......"
	run_firstlevel_spm "$1" "$2" "$3"

	echo "running reallife fsl estimation......"
	run_firstlevel_fsl "$1" "$2" "$3"
    else
	run_firstlevel_spm "$1" "$2" "$3"
    fi
}

# estimate firstlevel model using spm estimation
#
# in:
#  $1 model name or "" for none
#  $2 modelspace vol or surf
#  $3 hemisphere (if modelspace is surf)
run_firstlevel_spm() {
    modelname="$1";
    if [ "$2" == "surf" -o "$2" == "surface" ]; then
	modelspace=surf
	hemi="$3"
    else
	modelspace=vol
    fi

    firstleveldir=firstlevel

    # create a matlab script to do model estimation
    firstlevel_script=scripts/firstlevel

    # look for a model type modifier
    if [ "$modelname" ]; then
	firstleveldir="$firstleveldir"_"$modelname"
	firstlevel_script="$firstlevel_script"_"$modelname"
    fi

    # check if this is a surface analysis
    if [ "$modelspace" == "surf" ]; then
	firstleveldir="$firstleveldir"/$hemi
	firstlevel_script="$firstlevel_script"_$hemi
    fi

    firstlevel_script="$firstlevel_script".m

    echo "% auto-generated by $0" > $firstlevel_script

    if [ "$do_bayes_estimation" -eq "1" ]; then
	echo "load "$matlab_dir"/estimate_job_bayes.mat;" >> $firstlevel_script
    else
	echo "load "$matlab_dir"/estimate_job.mat;" >> $firstlevel_script
    fi
#    echo "addpath('scripts');" >> $firstlevel_script

    echo "spm_defaults;"  >> $firstlevel_script
    echo "jobs{1}.stats{1}.fmri_est(1).spmmat = {'$PWD/$firstleveldir/SPM.mat'};" >> $firstlevel_script

    echo "dbstop if error"  >> $firstlevel_script
    echo "spm_jobman('run',jobs);" >> $firstlevel_script
#    echo "keyboard;" >> $firstlevel_script
    "$script_dir"/run_matlab_script.sh "$firstlevel_script"

}

# estimate firstlevel model using spm estimation
#
# in:
#  $1 model name or "" for none
#  $2 modelspace vol or surf
#  $3 hemisphere (if modelspace is surf)
run_firstlevel_fsl() {
    modelname="$1";
    if [ "$2" == "surf" -o "$2" == "surface" ]; then
	modelspace=surf
	hemi="$3"
    else
	modelspace=vol
    fi

    # check that there is only one run
    morethanonerun=`echo $runnums | grep -c "[[:space:]]"`;
    if [ "$morethanonerun" == "1" ]; then
	echo "error: fsl model fit on more than a sinlge run is not supported yet"
	echo "runnums were $runnums"
	exit
    fi

    firstleveldir=firstlevel

    # create a matlab script to convert the spm model
    firstlevel_script=scripts/modelconvert

    # look for a model type modifier
    if [ "$modelname" ]; then
	firstleveldir="$firstleveldir"_"$modelname"
	firstlevel_script="$firstlevel_script"_"$modelname"
    fi

    spmfirstleveldir="$firstleveldir"
    firstleveldir="$firstleveldir"_fsl

    # check if this is a surface analysis
    if [ "$modelspace" == "surf" ]; then
	firstleveldir="$firstleveldir"/$hemi
	firstlevel_script="$firstlevel_script"_$hemi
    fi

    firstlevel_script="$firstlevel_script".m

    # remove old analysis
    if [ -e "$firstleveldir" ]; then
	rm -rf $firstleveldir
    fi

    echo mkdir -p "$PWD/$firstleveldir"
    mkdir -p "$PWD/$firstleveldir"

    # convert the spm design
    fsldesign=$PWD/"$firstleveldir"/design

    echo "% auto-generated by $0" > $firstlevel_script
    echo "addpath('"$matlab_dir"/');" > $firstlevel_script

    echo "spm2fslDesign('$PWD/$spmfirstleveldir/SPM.mat','$fsldesign');" >> $firstlevel_script

    echo "dbstop if error"  >> $firstlevel_script
#    echo "keyboard;" >> $firstlevel_script
    "$script_dir"/run_matlab_script.sh "$firstlevel_script"

    # run the estimation
    run_file="$fsldesign"_datafiles.txt
    echo "film_gls -rn $firstleveldir -sa -ms 5 `cat $run_file` design.mat 100.0"
    film_gls -rn "$firstleveldir" -sa -ms 5 `cat $run_file` "$fsldesign".mat 100.0

    # copy all results from the plus directory (grrrrr)
    mv "$firstleveldir"+/* "$firstleveldir"
    rmdir "$firstleveldir"+

    # copy the pes to betas
    echo "copying pes to  betas...."
    for pe in `ls $firstleveldir/pe*`; do
	num=`echo $pe | sed "s/.*pe\([0-9].*\).nii/\1/"`
	betaname=`printf "%s/beta_%04d.img" "$spmfirstleveldir" $num`
	echo mri_convert $pe $betaname
	mri_convert $pe $betaname
    done
}

# run firstlevel contrasts
#
# in:
#  $1 model name or "" for none
#  $2 estimated model dir (for betas) if different from default model dir
#  $3 modelspace vol or surf
#  $4 hemisphere (if modelspace is surf)
run_contrasts() {
    if [ "$fsl_estimate" -eq "1" ]; then
	run_contrasts_fsl "$1" "$2" "$3" "$4"
    else
	run_contrasts_spm "$1" "$2" "$3" "$4"
    fi
}

# run firstlevel contrasts using model estimated using spm
#
# in:
#  $1 model name or "" for none
#  $2 estimated model dir (for betas) if different from default model dir
#  $3 modelspace vol or surf
#  $4 hemisphere (if modelspace is surf)
run_contrasts_spm() {
    modelname="$1";
    estimated_model_dir="$2"

    if [ "$3" == "surf" -o "$3" == "surface" ]; then
	modelspace=surf
	hemi="$4"
    else
	modelspace=vol
    fi

    firstleveldir=firstlevel
    contrastmfun="build_contrasts"

    # create a matlab script to do constrasts
    contrasts_script=scripts/contrasts

    # look for a model type modifier
    if [ "$modelname" ]; then
	firstleveldir="$firstleveldir"_"$modelname"
	contrastmfun="$contrastmfun"_"$modelname"
	contrasts_script="$contrasts_script"_"$modelname"
    fi

    # look for a surface modelspace
    if [ "$modelspace" == "surf" ]; then
	firstleveldir="$firstleveldir"/$hemi
	contrastmfun="$contrastmfun"
	contrasts_script="$contrasts_script"_$hemi
    fi

    contrasts_script="$contrasts_script".m

    echo "% auto-generated by $0" > $contrasts_script
    echo "addpath('"$analysis_dir"/');" >> $contrasts_script

    echo "load $PWD/$firstleveldir/SPM;"  >> $contrasts_script
    echo "SPM.swd = '$PWD/$firstleveldir/';"  >> $contrasts_script
    echo "save $PWD/$firstleveldir/SPM.mat SPM;"  >> $contrasts_script

    echo "jobs{1}.stats{1}.con.spmmat = {'$PWD/$firstleveldir/SPM.mat'};"  >> $contrasts_script

    echo "jobs{1}.stats{1}.con.consess = $contrastmfun('$subj','$PWD/$firstleveldir/SPM.mat'$estimated_model_dir);"   >> $contrasts_script

    echo "dbstop if error"  >> $contrasts_script
#    echo "keyboard" >> $contrasts_script
    echo "spm_jobman('run',jobs);" >> $contrasts_script
    "$script_dir"/run_matlab_script.sh "$contrasts_script"
}

# run firstlevel fsl contrasts
#
# in:
#  $1 model name or "" for none
#  $2 estimated model dir (for betas) if different from default model dir
#  $3 modelspace vol or surf
#  $4 hemisphere (if modelspace is surf)
run_contrasts_fsl() {
    modelname="$1";
    estimated_model_dir="$2"

    if [ "$3" == "surf" -o "$3" == "surface" ]; then
	modelspace=surf
	hemi="$4"
    else
	modelspace=vol
    fi

    firstleveldir=firstlevel
    contrastmfun="build_contrasts"

    # create a matlab script to do constrasts
    contrasts_script=scripts/contrasts

    # look for a model type modifier
    if [ "$modelname" ]; then
	firstleveldir="$firstleveldir"_"$modelname"
	contrastmfun="$contrastmfun"_"$modelname"
	contrasts_script="$contrasts_script"_"$modelname"
    fi

    # look for a surface modelspace
    if [ "$modelspace" == "surf" ]; then
	firstleveldir="$firstleveldir"/$hemi
	contrastmfun="$contrastmfun"
	contrasts_script="$contrasts_script"_$hemi
    fi

    contrasts_script="$contrasts_script".m

    echo "% auto-generated by $0" > $contrasts_script
    echo "addpath('"$analysis_dir"/');" >> $contrasts_script

    echo "load $PWD/$firstleveldir/SPM;"  >> $contrasts_script
    echo "SPM.swd = '$PWD/$firstleveldir/';"  >> $contrasts_script
    echo "save $PWD/$firstleveldir/SPM.mat SPM;"  >> $contrasts_script

    echo "jobs{1}.stats{1}.con.spmmat = {'$PWD/$firstleveldir/SPM.mat'};"  >> $contrasts_script

    echo "jobs{1}.stats{1}.con.consess = $contrastmfun('$subj','$PWD/$firstleveldir/SPM.mat'$estimated_model_dir);"   >> $contrasts_script

    echo "dbstop if error"  >> $contrasts_script
#    echo "keyboard" >> $contrasts_script
    echo "spm_jobman('run',jobs);" >> $contrasts_script
    "$script_dir"/run_matlab_script.sh "$contrasts_script"
}

# convert contrasts for a model into paint files for surface visualization
#
# in:
#  $1 model name or "" for none
surf_paint() {
    modelname="$1";

    # look for a model type modifier
    firstleveldir=firstlevel
    if [ "$modelname" ]; then
	firstleveldir="$firstleveldir"_"$modelname"
    fi

    paint_script=scripts/paint.m

    # create surface viewable files of the contrasts
    echo "" > $paint_script
    echo "% create surface viewable files of the contrasts" >> $paint_script
    echo "addpath $matlab_dir" >> $paint_script

    for hemi in lh rh; do 

	numverts=`mris_info surf/$hemi.white 2> /dev/null | sed -n "/num vertices/p" | sed "s/.* //"`

	# build a list of files to unslice
	slifiles=
	if [ "$paint_beta" == 1 ]; then # beta images
            # zero nans in the beta images
	    for file in `ls $PWD/$firstleveldir/$hemi/beta_*.img`; do
		zeroed=`echo $file | sed s/beta_/beta_zeronan_/ | sed s/.img/.nii/`
		echo fslmaths $file -nan $zeroed
		fslmaths $file -nan $zeroed
		gunzip_file $zeroed
	    done
	    slifiles="$slifiles `ls $PWD/$firstleveldir/$hemi/beta_zeronan_????.nii`"
	fi
	if [ "$paint_con" == 1 ]; then # con images
            # zero nans in the con images
	    for file in `ls $PWD/$firstleveldir/$hemi/con_*.img`; do
		zeroed=`echo $file | sed s/con_/con_zeronan_/ | sed s/.img/.nii/`
		echo fslmaths $file -nan $zeroed
		fslmaths $file -nan $zeroed
		gunzip_file $zeroed
	    done
	    slifiles="$slifiles `ls $PWD/$firstleveldir/$hemi/con_zeronan_????.nii`"
	fi
	if [ "$paint_t" == 1 ]; then # t maps
	    slifiles="$slifiles `ls $PWD/$firstleveldir/$hemi/spmT_????.img`"
	fi

	if [ ! "$slifiles" ]; then
	    echo "no files to reslice"
	    return
	fi

	if [ ! "$numverts" ]; then
	    echo cant determine number of surface vertices.
	    echo maybe your freesurfer environment is not setup?
	    return 1
	fi

	for file in $slifiles; do
	    fsnii=`echo $file | sed "s/nii$\|img$/fs.nii/"`	    
	    echo "unreshape_big_nifti('$file',$numverts,'$fsnii');" >> $paint_script

	    #paint=`echo $fsnii | sed "s/.fs.nii/.mgh/"`
	    #echo "v = nifti('$file');"  >> $paint_script
	    #echo "save_mgh(v.dat(1:$numverts),'$paint',eye(4));" >> $paint_script

	done
    done
    "$script_dir"/run_matlab_script.sh "$paint_script"
}

#
#

#### main body of the run_subject.sh script #####

# configuration and validation

parse_args $@;

if [ "$config_filename" != "UNSET" ]; then
    read_configfile "$config_filename"
fi

validate_config;

# set up the start and stop stages
get_stage_number $startstage
start_stage_num=$stage_num
echo running from $startstage to $stopstage

get_stage_number $stopstage
stop_stage_num=$stage_num

#validate start and stop names
if [ "$start_stage_num" == "0" ]; then
    echo "error: invalid start stage name $startstage"
    usage
    exit 0;
fi
if [ "$stop_stage_num" == "0" ]; then
    echo "error: invalid stop stage name $stopstage"
    usage
    exit 0;
fi

dir=$PWD

if [ ! "$subj" ]; then
    usage
    exit
fi

subjdir="$study_subjs_dir/$subj"

if [ ! -d "$subjdir" ]; then
    echo "invalid subjects directory: $subjdir"
    usage
    exit
fi

if [ ! "$dicom_dir" ]; then
    # only the most recent mri data directory
    dicom_dir=`ls -dtr $subjdir/TrioTim* 2> /dev/null | tail -n 1`
fi

get_stage_number "extract_motion"
if [ "$start_stage_num" -le "$stage_num" ]; then
    echo "found dicom directory $dicom_dir"
    extract_motion $subjdir $dicom_dir
fi
if [ "$stop_stage_num" -le "$stage_num" ]; then
    exit;
fi

get_stage_number "extract_volumes"
if [ "$start_stage_num" -le "$stage_num" ]; then
    echo "found dicom directory $dicom_dir"
    extract_volumes $dicom_dir
fi
if [ "$stop_stage_num" -le "$stage_num" ]; then
    exit;
fi

mkdir -p scripts

# find the functional runs
if [ ! -e  scripts/func_runs.txt ]; then
    ls nii

    while [ 1 ] ; do
	echo "enter the run numbers of the functional runs (p to print again): "
	read runnums

	case runnums in
	    p ) ls nii;;
	    * ) echo "read $runnums"
		echo "enter 1 to accept, 0 to renter"
		read again
		if [ "$again" -eq 1 ]; then
		    break;
		fi
		;;
	esac
    done
    echo $runnums > scripts/func_runs.txt
else
    runnums=`cat scripts/func_runs.txt`
fi

get_stage_number "slice_timing"
if [ "$start_stage_num" -le "$stage_num" -a "$do_slice_timing" == 1 ]; then
    slice_timing "$runnums"
fi
if [ "$stop_stage_num" -le "$stage_num" ]; then
    exit;
fi
update_prefix "slice_timing"

get_stage_number "fieldmap_correct"
if [ "$start_stage_num" -le "$stage_num" -a "$do_fieldmap_correct" == 1 ]; then
    fieldmap_correct `echo $runnums s/ .*//`
fi
if [ "$stop_stage_num" -le "$stage_num" ]; then
    exit;
fi

get_stage_number "motion_correct"
if [ "$start_stage_num" -le "$stage_num" -a "$do_motion_correct" == 1 ]; then
    motion_correct "$runnums"
fi
if [ "$stop_stage_num" -le "$stage_num" ]; then
    exit;
fi
update_prefix "fieldmap_correct" # update the prefix here for fieldmap too coz
                                 # input to realign and reslice is normal runs
update_prefix "motion_correct"

get_stage_number "normalize"
if [ "$start_stage_num" -le "$stage_num" -a "$do_normalize" == 1 ]; then
    normalize "$runnums"
fi
if [ "$stop_stage_num" -le "$stage_num" ]; then
    exit;
fi
update_prefix "normalize"

get_stage_number "surf_project"
if [ "$start_stage_num" -le "$stage_num"  -a "$do_surf_project" == 1 ]; then
    surf_project "$runnums"
fi
if [ "$stop_stage_num" -le "$stage_num" ]; then
    exit;
fi
update_prefix "surf_project"

get_stage_number "smooth"
if [ "$start_stage_num" -le "$stage_num" ]; then
    # multiple smoothings are available
    if [ "$do_vol_smooth" == 1 ]; then
	vol_smooth "$runnums"
    fi

    if [ "$do_surf_vol_smooth" == 1 ]; then
	surf_vol_smooth "$runnums"
    fi

    if [ "$do_surf_smooth" == 1 ]; then
	surf_smooth "$runnums"
    fi
fi
if [ "$stop_stage_num" -le "$stage_num" ]; then
    exit;
fi
update_prefix "smooth"

echo prefix=$prefix

echo "Updating run nums for model: $model"
if [ "$model" -a "$same_runs_as_base_model" == 0 ]; then
    outlierfile=outliers_"$model".mat

# FIX: function this up. -twt
    if [ ! -e  scripts/func_runs_$model.txt ]; then
    	wc motion/*.txt

	while [ 1 ] ; do
	    echo "enter the numbers of the functional runs for this model (p to print again): "
	    read runnums

	    case runnums in
		p ) ls nii;;
		* ) echo "read $runnums"
		    echo "enter 1 to accept, 0 to renter"
		    read again
		    if [ "$again" -eq 1 ]; then
			break;
		    fi
		    ;;
	    esac
	done
	echo $runnums > scripts/func_runs_$model.txt
    else
	runnums=`cat scripts/func_runs_$model.txt`
    fi
elif [ "$same_runs_as_base_model" == 1 ]; then
    echo "assuming runnums in scripts/func_runs.txt is correct"
    outlierfile=outliers.mat
else
    outlierfile=outliers.mat
fi

echo "Run nums for model \"$model\": $runnums"
echo "Outlier file: $outlierfile"

get_stage_number "outliers"
if [ "$start_stage_num" -le "$stage_num" -a "$do_outliers" == 1 ]; then
    # view artifacts using art
    if [ -e nii/$outlierfile -a "$skip_art_if_outliers_file_exists" == 1 ]; then
	echo "skipping ART because skip_art_if_outliers_file_exists is 1 "
    else
	art_batch "$runnums"
    fi
fi
if [ "$stop_stage_num" -le "$stage_num" ]; then
    exit;
fi

get_stage_number "view_funcruns"
if [ "$start_stage_num" -le "$stage_num" -a "$do_view_funcruns" == 1 ]; then
    view_funcruns "$runnums"
fi
if [ "$stop_stage_num" -le "$stage_num" ]; then
    exit;
fi

get_stage_number "prepare_analysis"
if [ "$start_stage_num" -le "$stage_num"  -a "$do_prepare_analysis" == 1 ]; then
    echo preparing analysis with model=$model...

    if [ "$do_vol_estimate" == 1 ]; then
	build_model "$runnums" "$model" volume
    fi

    if [ "$do_surf_estimate" == 1 ]; then
	build_model "$runnums" "$model" surf lh
	build_model "$runnums" "$model" surf rh
    fi
fi
if [ "$stop_stage_num" -le "$stage_num" ]; then
    exit;
fi

# run first level
get_stage_number "run_firstlevel"
if [ "$start_stage_num" -le "$stage_num"  -a "$do_run_firstlevel" == 1 ]; then
    echo running first level with model=$model...

    if [ "$do_vol_estimate" == 1 ]; then
	run_firstlevel "$model" volume
    fi

    if [ "$do_surf_estimate" == 1 ]; then
	run_firstlevel "$model" surf lh
	run_firstlevel "$model" surf rh
    fi
fi
if [ "$stop_stage_num" -le "$stage_num" ]; then
    exit;
fi

# run contrasts
get_stage_number "run_contrasts"
if [ "$start_stage_num" -le "$stage_num"  -a "$do_run_contrasts" == 1 ]; then
    echo running contrasts with model=$model...

    if [ "$do_vol_estimate" == 1 ]; then
	run_contrasts "$model" "$est_model_dir" volume
    fi

    if [ "$do_surf_estimate" == 1 ]; then
	run_contrasts "$model" "$est_model_dir" surf lh
	run_contrasts "$model" "$est_model_dir" surf rh
    fi
fi
if [ "$stop_stage_num" -le "$stage_num" ]; then
    exit;
fi

# 
get_stage_number "surf_paint"
if [ "$start_stage_num" -le "$stage_num"  -a "$do_surf_paint" == 1 ]; then
    echo running surf_paint with model=$model...
    surf_paint "$model"
fi
if [ "$stop_stage_num" -le "$stage_num" ]; then
    exit;
fi