We have been very interested in comparing slomoco to our normal preprocessing procedure.
Initial results indicate that our volume registration strategy is better (our max displacement values are about 1/2 of slomoco's values).  Slomoco, however, does 10x as well in subsequent despiking.  We presume the max displacement difference is because slomoco registers to the first volume and we choose the volume which has the lowest outlier profile and register to that.  Even with trimming of about 6 initial volumes...those early volumes are still not as good for volume registration.

We'd love to have the best of both world's, and so I added some code to slicemoco_newalgorithm.sh to get it to take a base volume from the commandline.  This seems to work and I get better looking maximum displacement profiles on early tests...however, about half the time when I run, I get a matlab crash dump (attached example).  I am on OSX Mavericks, Java 7, Matlab 2014B.  I have set the matlab java variable in my .cshrc, set DYLD Library paths in cshrc and matlab startup.m.  I've given matlab lots of java heap memory....I've tried running just one process at a time.  Sometimes everything works, and sometimes I get the matlab crash.  I can reprocess a run that failed initially and have it succeed. 

setenv DYLD_LIBRARY_PATH ""
setenv DYLD_VERSIONED_LIBRARY_PATH # ${AFNI_HOME} #:/opt/X11/lib:/opt/lib:/opt/local/lib
setenv DYLD_FALLBACK_LIBRARY_PATH ${AFNI_HOME}:/opt/local/lib:/usr/lib #:/opt/lib:/opt/X11/lib #/opt/local/lib
setenv MATLAB_JAVA /Library/Java/JavaVirtualMachines/jdk1.7.0_71.jdk/Contents/Home/jre

I think this alternative approach has real merit, but I don't understand the crash behavior.  If you can offer any insights I'd be very grateful.

-Dianne Patterson, Ph.D.
dkp@email.arizona.edu

=================
The slicemoco2.sh code is pasted below:
#! /bin/bash

function Usage () {
cat < slicemoco2.sh is a revision of
slicemoco_newalgorithm (distributed with PESTICA v2.1)
Algorithm: this script first runs slicewise in-plane (xy) 3DOF motion correction
then runs a slicewise 6DOF rigid-body correction for each slice
this script reads these motion parameters back in and regress on voxel timeseries
WARNING, make sure you have removed unsaturated images at start (such as first 4 vols)
You can test first volumes for spin saturation with: 3dToutcount | 1dplot -stdin -one
Is the first volume much higher than rest? If so, you may need to remove first several volumes first
If you don't know what this means, consult someone who does know, this is very important,
regression corrections (and analyses) perform poorly when the data has unsaturated volumes at the start
Simple method to remove 1st 4: 3dcalc -a "+orig[4..$]" -expr a -prefix .steadystate
Usage: slicemoco2.sh -d -b
-d = dataset: is the file prefix of
the 3D+time dataset to use and correct.
Note: this script will detect suffix for epi_filename
-s = slice_timing integer (default if not specified is = 1)
1=alternating interleaved ascending Siemens order (odd/even for odd#, even/odd for even#)
2=sequential ascending
3=alternating interleaved ascending GE/Phillips order (odd/even)
-r = perform in parallel with final PESTICA/IRF-RETROICOR regression correction
this assumes PESTICA estimation steps 1-6 have been run and exist in subdir _pestica/
this will pick up the final estimated IRFs and phase files and run final regression in parallel
-p = same, but assuming you used PMU data instead of PESTICA for the IRF-determination and correction
-b = base registration volume to be used by volreg (added by DP 2/22/2015)
EOF
exit 1
}
# default slice acq order is 1
SLOMOCO_SLICE_ORDER=1
phypest=0
phypmu=0
while getopts hd:s:rpb: opt; do
case "$opt" in
h)
Usage
exit 1
;;
d) # base 3D+time EPI dataset to use to perform corrections
epi=$OPTARG
;;
s) # slice timing integer, interleaved asc or sequential ascending supported
SLOMOCO_SLICE_ORDER=$OPTARG
;;
r)
phypest=1
;;
p)
phypmu=1
;;
b) # base registration volume for volreg
base_reg_num=$OPTARG
;;
:)
echo "option requires input"
exit 1
;;
esac
done
echo "================="
echo "slomoco2 knows volume ${base_reg_num} should be used for registration"
echo "================="
if [ $phypest -eq 1 ] ; then
if [ $phypmu -eq 1 ] ; then
echo "cannot set physio correction to both PMU and PESTICA, choose the appropriate one"
exit 1
fi
fi
# test if epi filename was set
if [ -z $epi ] ; then
echo "3D+time EPI dataset $epi must be given"
Usage
exit 1
fi
# test if base_reg_num was set
if [ -z $base_reg_num ] ; then
echo "base volume for registration is number ${base_reg_num}"
fi
# test for presence of input EPI file with one of the accepted formats and set suffix
if [ -f $epi.hdr ] ; then
suffix=".hdr"
elif [ -f $epi+orig.HEAD ] ; then
suffix="+orig.HEAD"
elif [ -f $epi.nii ] ; then
suffix=".nii"
elif [ -f $epi.nii.gz ] ; then
suffix=".nii.gz"
else
echo "3D+time EPI dataset $epi must exist, check filename (do not give suffix)"
echo "accepted formats: hdr +orig nii nii.gz"
Usage
echo ""
echo "***** $epi does not exist, exiting ..."
echo "accepted formats: hdr +orig nii nii.gz"
exit 1
fi
fullcommand="$0"
if [ -z $SLCMOCO_DIR ] ; then
echo "setting SLCMOCO_DIR to /MATLAB/slomoco_afni_07_2014"
SLCMOCO_DIR=/MATLAB/slomoco_afni_07_2014
#echo "setting SLCMOCO_DIR to directory where this script is located (not the CWD)"
#SLCMOCO_DIR='dirname $fullcommand'
# this command resides in the moco/ subdirectory of PESTICA_DIR
export SLCMOCO_DIR=$SLCMOCO_DIR
fi
# first test if we are running run_pestica.sh from the base SLCMOCO_DIR
homedir='pwd'
if [ $homedir == $SLCMOCO_DIR ] ; then
echo "you cannot run PESTICA from the downloaded/extracted SLCMOCO_DIR"
echo "please run this from the directory containing the data (or copy of the data)"
echo "that you want to correct. Exiting..."
exit 1
fi
# create PESTICA subdirectory if it does not exist
epi_pestica="$epi"_pestica
if [ ! -d $epi_pestica ] ; then
echo ""
echo "* Creating PESTICA Directory: $epi_pestica"
echo ""
mkdir $epi_pestica
echo "mkdir $epi_pestica" >> $epi_pestica/pestica_history.txt
echo "" >> $epi_pestica/pestica_history.txt
fi
echo "***** Using copied input in $epi_pestica/$epi+orig.HEAD as input timeseries"
# always copy input file into PESTICA subdirectory in AFNI format
if [ ! -f $epi_pestica/$epi+orig.HEAD ] ; then
echo "Copying: 3dcopy $epi$suffix $epi_pestica/$epi"
3dcopy $epi$suffix $epi_pestica/$epi
echo "3dcopy $epi$suffix $epi_pestica/$epi" >> $epi_pestica/pestica_history.txt
echo "" >> $epi_pestica/pestica_history.txt
fi
# write command line and SLCMOCO_DIR to history file
echo "'date'" >> $epi_pestica/pestica_history.txt
echo "PESTICA_afni-v2 command line: 'basename $fullcommand' $*" >> $epi_pestica/pestica_history.txt
echo "PESTICA env:
'env | grep PESTICA'" >> $epi_pestica/pestica_history.txt
echo "" >> $epi_pestica/pestica_history.txt
# do ALL work inside the pestica/ subdirectory
homedir='pwd'
if [ -f $epi_mask+orig.HEAD ] ; then
3dcopy $epi_mask+orig $epi_pestica/$epi_mask
fi
cd $epi_pestica
echo "cd $epi_pestica" >> pestica_history.txt
echo "" >> pestica_history.txt
epi_mask=${epi}.brain
#test for presence of mask file
if [ ! -f $epi_mask+orig.HEAD ] ; then
echo ""
echo "***** $epi_mask+orig.HEAD does not exist, creating mask"
echo "note, if you wish to use your own mask/brain file, kill this script"
echo "then 3dcopy your mask/brain to $epi_mask and re-run."
echo "PESTICA will use whatever resides in $epi_mask"
echo ""
echo "running 3dSkullStrip -input $epi+orig -prefix $epi_mask"
sleep 1
3dSkullStrip -input $epi+orig -prefix ___tmp_mask
# erode mask by one voxel
3dcalc -a ___tmp_mask+orig -prefix ___tmp_mask_ones+orig -expr 'step(a)'
3dcalc -a ___tmp_mask_ones+orig -prefix ___tmp_mask_ones_dil -b a+i -c a-i -d a+j -e a-j -f a+k -g a-k -expr 'amongst(1,a,b,c,d,e,f,g)'
3dcalc -a "$epi+orig[0]" -b ___tmp_mask_ones_dil+orig -prefix $epi_mask -expr 'a*step(b)'
rm ___tmp_mask*
echo ""
echo "done with skull-stripping - please check file and if not satisfied, I recommend running"
echo "3dSkullStrip with different parameters to attempt to get a satisfactory brain mask."
echo "Either way, this script looks in $epi_pestica/ for $epi_mask to use as your brain mask/strip"
echo " if 3dSkullStrip dies due to large size or other unknown reason, try: "
echo " 3dAutomask -dilate 1 -clfrac 0.4 -prefix $epi_pestica/$epi.brain $epi+orig"
sleep 1
echo "3dSkullStrip -input $epi+orig -prefix ___tmp_mask" >> pestica_history.txt
echo "3dcalc -a ___tmp_mask+orig -prefix ___tmp_mask_ones+orig -expr 'step(a)'" >> pestica_history.txt
echo "3dcalc -a ___tmp_mask_ones+orig -prefix ___tmp_mask_ones_dil -b a+i -c a-i -d a+j -e a-j -f a+k -g a-k -expr 'amongst(1,a,b,c,d,e,f,g)'" >> pestica_history.txt
echo "3dcalc -a "$epi+orig[0]" -b ___tmp_mask_ones_dil+orig -prefix $epi_mask -expr 'a*step(b)'" >> pestica_history.txt
echo "rm ___tmp_mask*" >> pestica_history.txt
echo "" >> pestica_history.txt
fi
echo "***** Using $epi_mask+orig.HEAD to mask out non-brain voxels"
if [ ! -f mocoafni.txt ] ; then
3dvolreg -prefix $epi.mocoafni.hdr -base ${base_reg_num} -zpad 8 -maxite 60 -x_thresh 0.005 -rot_thresh 0.008 -verbose -dfile $epi.mocoafni.txt -1Dfile $epi.mocoafni.1D -maxdisp1D $epi.mocoafni.maxdisp.1D -heptic $epi+orig
rm $epi.mocoafni.hdr $epi.mocoafni.img
cp $epi.mocoafni.txt mocoafni.txt
fi
if [ ! -d tempslmocoxy_afni_$epi ] ; then
echo "$SLCMOCO_DIR/run_correction_slicemocoxy_afni.sh -b $epi -p $epi.slicemocoxy_afni" >> pestica_history.txt
$SLCMOCO_DIR/run_correction_slicemocoxy_afni.sh -b $epi -p $epi.slicemocoxy_afni
fi
if [ ! -d tempslmoco_volslc_alg_vol_$epi.slicemocoxy_afni ] ; then
echo "$SLCMOCO_DIR/run_slicemoco_inside_fixed_vol.sh -b $epi.slicemocoxy_afni" >> pestica_history.txt
$SLCMOCO_DIR/run_slicemoco_inside_fixed_vol.sh -b $epi.slicemocoxy_afni
fi
echo ""
echo "Running Secondorder Motion Correction using SLOMOCO output"
echo ""
if [ $phypest -eq 1 ] ; then
echo matlab $MATLABLINE "<<<\"addpath $SLCMOCO_DIR; addpath $MATLAB_AFNI_DIR; addpath $PESTICA_DIR; addpath $PESTICA_DIR/matlab_retroicor; mocoparams=read_motion_newslicealg('tempslmoco_volslc_alg_vol_$epi.slicemocoxy_afni/motion.wholevol_zt'); load impulse_responses_secondpass.mat; slicemoco_newalgorithm_input_physio('$epi.slicemocoxy_afni+orig','$epi_mask+orig',mocoparams,cardph,respph,coeffs_card,coeffs_resp,$SLOMOCO_SLICE_ORDER); exit;\""
echo matlab $MATLABLINE "<<<\"addpath $SLCMOCO_DIR; addpath $MATLAB_AFNI_DIR; addpath $PESTICA_DIR; addpath $PESTICA_DIR/matlab_retroicor; mocoparams=read_motion_newslicealg('tempslmoco_volslc_alg_vol_$epi.slicemocoxy_afni/motion.wholevol_zt'); load impulse_responses_secondpass.mat; slicemoco_newalgorithm_input_physio('$epi.slicemocoxy_afni+orig','$epi_mask+orig',mocoparams,cardph,respph,coeffs_card,coeffs_resp,$SLOMOCO_SLICE_ORDER); exit;\"" >> pestica_history.txt
matlab $MATLABLINE <<<"addpath $SLCMOCO_DIR; addpath $MATLAB_AFNI_DIR; addpath $PESTICA_DIR; addpath $PESTICA_DIR/matlab_retroicor; mocoparams=read_motion_newslicealg('tempslmoco_volslc_alg_vol_$epi.slicemocoxy_afni/motion.wholevol_zt'); load impulse_responses_secondpass.mat; disp('Wait, script starting...'); slicemoco_newalgorithm_input_physio('$epi.slicemocoxy_afni+orig','$epi_mask+orig',mocoparams,cardph,respph,coeffs_card,coeffs_resp,$SLOMOCO_SLICE_ORDER); exit;"
corrstr="slicemocoxy_afni.zalg_moco2_irfret"
elif [ $phypmu -eq 1 ] ; then
echo matlab $MATLABLINE "<<<\"addpath $SLCMOCO_DIR; addpath $MATLAB_AFNI_DIR; addpath $PESTICA_DIR; addpath $PESTICA_DIR/matlab_retroicor; mocoparams=read_motion_newslicealg('tempslmoco_volslc_alg_vol_$epi.slicemocoxy_afni/motion.wholevol_zt'); load impulse_responses_pmu_secondpass.mat; slicemoco_newalgorithm_input_physio('$epi.slicemocoxy_afni+orig','$epi_mask+orig',mocoparams,cardph,respph,coeffs_card,coeffs_resp,$SLOMOCO_SLICE_ORDER); exit;\""
echo matlab $MATLABLINE "<<<\"addpath $SLCMOCO_DIR; addpath $MATLAB_AFNI_DIR; addpath $PESTICA_DIR; addpath $PESTICA_DIR/matlab_retroicor; mocoparams=read_motion_newslicealg('tempslmoco_volslc_alg_vol_$epi.slicemocoxy_afni/motion.wholevol_zt'); load impulse_responses_pmu_secondpass.mat; slicemoco_newalgorithm_input_physio('$epi.slicemocoxy_afni+orig','$epi_mask+orig',mocoparams,cardph,respph,coeffs_card,coeffs_resp,$SLOMOCO_SLICE_ORDER); exit;\"" >> pestica_history.txt
matlab $MATLABLINE <<<"addpath $SLCMOCO_DIR; addpath $MATLAB_AFNI_DIR; addpath $PESTICA_DIR; addpath $PESTICA_DIR/matlab_retroicor; mocoparams=read_motion_newslicealg('tempslmoco_volslc_alg_vol_$epi.slicemocoxy_afni/motion.wholevol_zt'); load impulse_responses_pmu_secondpass.mat; disp('Wait, script starting...'); slicemoco_newalgorithm_input_physio('$epi.slicemocoxy_afni+orig','$epi_mask+orig',mocoparams,cardph,respph,coeffs_card,coeffs_resp,$SLOMOCO_SLICE_ORDER); exit;"
corrstr="slicemocoxy_afni.zalg_moco2_irfret"
else
echo matlab $MATLABLINE "<<<\"addpath $SLCMOCO_DIR; mocoparams=read_motion_newslicealg('tempslmoco_volslc_alg_vol_$epi.slicemocoxy_afni/motion.wholevol_zt'); slicemoco_newalgorithm_input('$epi.slicemocoxy_afni+orig','$epi_mask+orig',mocoparams,$SLOMOCO_SLICE_ORDER); exit;\""
echo matlab $MATLABLINE "<<<\"addpath $SLCMOCO_DIR; mocoparams=read_motion_newslicealg('tempslmoco_volslc_alg_vol_$epi.slicemocoxy_afni/motion.wholevol_zt'); slicemoco_newalgorithm_input('$epi.slicemocoxy_afni+orig','$epi_mask+orig',mocoparams,$SLOMOCO_SLICE_ORDER); exit;\"" >> pestica_history.txt
echo ""
matlab $MATLABLINE <<<"addpath $SLCMOCO_DIR; mocoparams=read_motion_newslicealg('tempslmoco_volslc_alg_vol_$epi.slicemocoxy_afni/motion.wholevol_zt'); disp('Wait, script starting...'); slicemoco_newalgorithm_input('$epi.slicemocoxy_afni+orig','$epi_mask+orig',mocoparams,$SLOMOCO_SLICE_ORDER); disp('2nd order Moco Done!'); exit;"
corrstr="slicemocoxy_afni.zalg_moco2"
fi
matlab $MATLABLINE <<<"addpath $SLCMOCO_DIR; qa_slomoco('$epi.slicemocoxy_afni+orig','tempslmoco_volslc_alg_vol_$epi.slicemocoxy_afni/motion.wholevol_zt','tempslmocoxy_afni_$epi',$SLOMOCO_SLICE_ORDER); exit;"
cd $homedir
# change attributes to match original data
taxis_nums='3dAttribute TAXIS_NUMS $epi_pestica/$epi+orig'
taxis_floats='3dAttribute TAXIS_FLOATS $epi_pestica/$epi+orig'
taxis_offset='3dAttribute TAXIS_OFFSETS $epi_pestica/$epi+orig'
#change TR back to input dataset TR
#3drefit -TR $tr $epi+orig $epi.slicemocoxy_afni.zalg_moco2+orig
#copy over and save t-axis nums into new image registered dataset
3drefit -saveatr -atrint TAXIS_NUMS "$taxis_nums" $epi.$corrstr+orig
#copy over and save t-axis floats into new image
3drefit -saveatr -atrfloat TAXIS_FLOATS "$taxis_floats" $epi.$corrstr+orig
#copy over and save t-axis offsets into new image registered dataset
3drefit -saveatr -atrfloat TAXIS_OFFSETS "$taxis_offset" $epi.$corrstr+orig
3dNotes -h "slicewise_moco_inplane.sh $epi+orig" $epi.$corrstr+orig
echo "End of PESTICA script" >> $epi_pestica/pestica_history.txt
echo "'date'" >> $epi_pestica/pestica_history.txt
echo "" >> $epi_pestica/pestica_history.txt

Hi Dianne,

This is an excellent idea. The main reason I used the first volume was because I want to make this real-time. Thats not working yet, and if its hampering performance, then it was a bad idea! I think this is great, and I think it would be a good addition to SLOMOCO.

I looked at the crash dump, but I couldn't see anything that jumps out to what little I know about the inner workings of MATLAB (other than the interpreter). Its possible the interpreter in newer versions is less accommodating to text input from a batch file, and as of this week I now have a new version of MATLAB (2015) that I need to do a lot of work on. It is possible these are related, and if you're interested, let me know and I can get you any changes before the next version release goes out.

Erik