vcc_calib

#!/bin/sh
set -e


if [ ! -e $TOOLBOX_PATH/bart ] ; then
	echo "\$TOOLBOX_PATH is not set correctly!" >&2
	exit 1
fi
export PATH=$TOOLBOX_PATH:$PATH
export BART_COMPAT_VERSION="v0.4.04"


CROP=0.
KRN=6
REG=24
LOW=20
THRESH=0.001

while getopts "c:k:r:l:t:" opt; do
	case $opt in
	l)
		LOW=$OPTARG
	;;
	c)
		CROP=$OPTARG
	;;
	k)
		KRN=$OPTARG
	;;
	r)
		REG=$OPTARG
	;;
	t)
		THRESH=$OPTARG
	;;
	\?)
		echo "Invalid option: -$OPTARG" >&2
	;;
	esac
done

shift $((OPTIND-1))


if [ $# -lt 2 ] ; then

        echo "Usage: $0 <kspace> <maps>"
        exit 1
fi


kspace=$(readlink -f "$1")
maps=$(readlink -f "$2")



#WORKDIR=$(mktemp -d)
# Mac: http://unix.stackexchange.com/questions/30091/fix-or-alternative-for-mktemp-in-os-x
WORKDIR=`mktemp -d 2>/dev/null || mktemp -d -t 'mytmpdir'`
trap 'rm -rf "$WORKDIR"' EXIT
cd $WORKDIR

# create virtual conjugate channels
bart flip 7 $kspace tmp1
bart circshift 0 1 tmp1 tmp2
bart circshift 1 1 tmp2 tmp1
bart circshift 2 1 tmp1 tmp2
bart conj tmp2 vcs


# calibration
bart join 3 $kspace vcs both


bart ecalib -k${KRN} -r${REG} -t${THRESH} -S -c${CROP} both bsens val
rm both.*



# phase centering - rotate phase up to sign
COILS=$(bart show -d3 bsens)
HCOILS=$(expr ${COILS} \/ 2)


bart extract 3 0 ${HCOILS} bsens maps1.coo
bart extract 3 ${HCOILS} ${COILS} bsens maps2.coo
bart fmac -s8 maps1.coo maps2.coo rel.coo
bart cpyphs rel.coo ph.coo
bart spow 0.5 ph.coo phsqrt.coo
bart fmac -C maps1.coo phsqrt.coo sens

rm bsens.*



# extract low resolution phase
bart caldir ${LOW} $kspace sensLR
bart cpyphs sensLR ph
rm sensLR.*

bart scale  -- 1.i ph phi
bart join 4 ph phi phc
rm ph.* phi.*


# align sign
bart fmac -s8 -C sens phc tmp1.coo
bart creal tmp1.coo tmp2.coo
bart cpyphs tmp2.coo tmp1.coo
bart fmac sens tmp1.coo $maps