diff --git a/bedpostx b/bedpostx
index 8bbbde137b6dd1c3f3cf66a5be0186c550861943..9a20960e55ac95868c559ec4e1d7416bbaced8cd 100755
--- a/bedpostx
+++ b/bedpostx
@@ -18,10 +18,14 @@ fi
Usage() {
echo ""
- echo "Usage: bedpostX <subject directory> "
+ echo "Usage: bedpostX <subject directory> [options]"
echo ""
echo "expects to find bvals and bvecs in subject directory"
echo "expects to find data nodif_brain_mask nodif in subject directory"
+ echo "options:"
+ echo "-n (number of fibres per voxel, default 2)"
+ echo "-w (ARD weight, more weight means less fibres killed per voxel, default 1)"
+ echo "-b (burnin period, default 1000)"
echo ""
exit
}
@@ -42,6 +46,23 @@ make_absolute(){
[ "$1" = "" ] && Usage
+
+#parse option arguments
+nfibres=2
+fudge=1
+burnin=1000
+while [ ! -z "$2" ]
+do
+ case "$2" in
+ -n) nfibres=$3;shift;;
+ -w) fudge=$3;shift;;
+ -b) burnin=$3;shift;;
+ *) break;;
+ esac
+ shift
+done
+
+
subjdir=`make_absolute $1`
subjdir=`echo $subjdir | sed 's/\/$/$/g'`
@@ -82,7 +103,6 @@ mkdir -p ${subjdir}.bedpostX/logs/pid_${$}
mkdir -p ${subjdir}.bedpostX/xfms
mailto=`whoami`@fmrib.ox.ac.uk
-
echo Queuing preprocessing stages
preprocid=`qsub -M $mailto -V -q short.q \
-o ${subjdir}.bedpostX/logs -e ${subjdir}.bedpostX/logs \
@@ -90,10 +110,11 @@ preprocid=`qsub -M $mailto -V -q short.q \
awk '{print $3}'`
echo Queuing parallel processing stage
+echo I am using this file
nslices=`avwval ${subjdir}/data dim3`
bedpostid=`qsub -M $mailto -hold_jid $preprocid -V -q long.q -t 1-$nslices \
-o ${subjdir}.bedpostX/logs -e ${subjdir}.bedpostX/logs \
- ${FSLDIR}/bin/sge_bedpostX_single_slice.sh $subjdir | \
+ ${FSLDIR}/bin/sge_bedpostX_single_slice.sh $subjdir $nfibres $fudge $burnin | \
awk '{print $3}' | awk -F. '{print $1}'`
echo Queuing post processing stage
@@ -105,7 +126,7 @@ mergeid=`qsub -M $mailto -hold_jid $bedpostid -V -q short.q \
echo
echo Type ${subjdir}.bedpostX/monitor to show progress.
cat <<EOM > ${subjdir}.bedpostX/monitor
-#!/bin/sh
+# #!/bin/sh
finished=0
logdir=${subjdir}.bedpostX/logs
while [ \$finished -eq 0 ] ; do
@@ -126,7 +147,7 @@ chmod +x ${subjdir}.bedpostX/monitor
echo Type ${subjdir}.bedpostX/cancel to terminate the task.
cat <<EOC > ${subjdir}.bedpostX/cancel
-#!/bin/sh
+# #!/bin/sh
qdel $mergeid $bedpostid $preprocid
EOC
chmod +x ${subjdir}.bedpostX/cancel
diff --git a/bedpostx_single_slice.sh b/bedpostx_single_slice.sh
index 9a851a040b5a6022a25d48d930e914a386538e13..c6c96d389c3c8d823f5fd000146f130d8a88ddfc 100755
--- a/bedpostx_single_slice.sh
+++ b/bedpostx_single_slice.sh
@@ -14,5 +14,5 @@ ${FSLDIR}/bin/xfibres\
--mask=$subjdir/nodif_brain_mask_slice_$slicezp\
-b $subjdir/bvals -r $subjdir/bvecs\
--forcedir --logdir=$subjdir.bedpostX/diff_slices/data_slice_$slicezp\
- --nj=1000 --bi=600 --se=20 --upe=24 --nfibres=2 > $subjdir.bedpostX/logs/log$slicezp && echo Done
+ --fudge=$3 --nj=5000 --bi=$4 --se=100 --upe=24 --nfibres=$2 > $subjdir.bedpostX/logs/log$slicezp && echo Done
diff --git a/old_bedpostX b/old_bedpostX
index 4be5dae2bdbf0bd6c4eaecfb012ef565dae28202..7b499c54ce782e1ce8780703772b54c7c02e9a71 100644
--- a/old_bedpostX
+++ b/old_bedpostX
@@ -8,10 +8,14 @@ fi
Usage() {
echo ""
- echo "Usage: bedpostX <subject directory> "
+ echo "Usage: bedpostX <subject directory> [options]"
echo ""
echo "expects to find bvals and bvecs in subject directory"
echo "expects to find data nodif_brain_mask nodif in subject directory"
+ echo "options:"
+ echo "-n (number of fibres per voxel, default 2)"
+ echo "-w (ARD weight, more weight means less fibres killed per voxel, default 1)"
+ echo "-b (burnin period, default 1000)"
echo ""
exit
}
@@ -104,6 +108,22 @@ fi
Lock;
+#parse option arguments
+nfibres=2
+fudge=1
+burnin=1000
+while [ ! -z "$2" ]
+do
+ case "$2" in
+ -n) nfibres=$3;shift;;
+ -w) fudge=$3;shift;;
+ -b) burnin=$3;shift;;
+ *) break;;
+ esac
+ shift
+done
+
+
subjdir=`make_absolute $1`
subjdir=`echo $subjdir | sed 's/\/$/$/g'`
@@ -167,14 +187,14 @@ ${FSLDIR}/bin/avwslice ${subjdir}/nodif_brain_mask
if [ "x$FSLMACHINELIST" = "x" ] ; then
echo "processing data on local host"
- ${FSLDIR}/bin/bedpostX_proc $subjdir $nslices ${subjdir}.bedpostX/logs/pid_${$} &
+ ${FSLDIR}/bin/bedpostX_proc $subjdir $nslices -n $nfibres -w $fudge -b $burnin ${subjdir}.bedpostX/logs/pid_${$} &
else
echo "processing data on hosts: $FSLMACHINELIST"
for machine in $FSLMACHINELIST; do
echo "if [ -r /usr/local/etc/fslconf/fsl.sh ];then . /usr/local/etc/fslconf/fsl.sh;fi; if [ -r /etc/fslconf/fsl.sh ];then . /etc/fslconf/fsl.sh;fi; if [ -r \${HOME}/.fslconf/fsl.sh ]; then . \${HOME}/.fslconf/fsl.sh; fi; if [ x\${FSLDIR} != "x" ];then \${FSLDIR}/bin/bedpostX_proc $subjdir $nslices ${subjdir}.bedpostX/logs/pid_${$}; else echo FSLDIR not set in any default location on machine `hostname`;fi"| $FSLREMOTECALL $machine /bin/sh &
done
fi
-
+echo finished
b=0
finished=0;
@@ -212,6 +232,9 @@ while [ $fib -le $numfib ];do
fib=`echo "$fib +1"|bc`;
done
+${FSLDIR}/bin/avwmerge -z ${subjdir}.bedpostX/meanSignal `${FSLDIR}/bin/imglob -oneperimage ${subjdir}.bedpostX/diff_slices/data_slice_*/meanSignal`
+${FSLDIR}/bin/avwmerge -z ${subjdir}.bedpostX/stdSignal `${FSLDIR}/bin/imglob -oneperimage ${subjdir}.bedpostX/diff_slices/data_slice_*/stdSignal`
+
echo Removing intermediate files
if [ `imtest ${subjdir}.bedpostX/merged_th1samples` -eq 1 ];then
diff --git a/old_bedpostX_proc b/old_bedpostX_proc
index 4cdbd5f61af3a28ae1ffc54768b1091c5c65d28a..56d192ee61a62be7a3e1988839c42fdd8069de9e 100644
--- a/old_bedpostX_proc
+++ b/old_bedpostX_proc
@@ -6,21 +6,36 @@
Usage() {
echo ""
- echo "Usage: bedpostX_proc <subject_dir> <nslices> [piddir]"
+ echo "Usage: bedpostX_proc <subject_dir> <nslices> [-n nfibres -w ard_weight -b burnin_period] [piddir]"
echo ""
exit
}
[ "$2" = "" ] && Usage
-[ "$3" = "" ] || touch ${3}/`hostname`_fdt_${$}
+#parse option arguments
+nfibres=2
+fudge=1.5
+burnin=2000
+while [ ! -z "$3" ]
+do
+ case "$3" in
+ -n) nfibres=$4;shift;;
+ -w) fudge=$4;shift;;
+ -b) burnin=$4;shift;;
+ *) piddir=$3;;
+ esac
+ shift
+done
-subjdir=$1
-nslices=$2
+[ "$piddir" = "" ] || touch ${piddir}/`hostname`_fdt_${$}
+subjdir=$1
+nslices=$2
+
slice=0
while [ $slice -lt $nslices ];do
slicezp=`${FSLDIR}/bin/zeropad $slice 4`
@@ -28,13 +43,13 @@ while [ $slice -lt $nslices ];do
echo `hostname`_${$} > ${subjdir}.bedpostX/logs/.$slicezp
sleep 10
if [ `hostname`_${$} = `cat ${subjdir}.bedpostX/logs/.$slicezp | sed -n '1p'` ] ; then
- nice ${FSLDIR}/bin/xfibres --data=$subjdir/data_slice_$slicezp --mask=$subjdir/nodif_brain_mask_slice_$slicezp -b $subjdir/bvals -r $subjdir/bvecs --forcedir --logdir=$subjdir.bedpostX/diff_slices/data_slice_$slicezp --nj=1000 --bi=600 --bn=0 --se=20 --nfibres=2 --fudge=1.5> $subjdir.bedpostX/logs/log$slicezp
-#bi can be changed to 2000 if convergence looks dodgy
+ nice ${FSLDIR}/bin/xfibres --data=$subjdir/data_slice_$slicezp --mask=$subjdir/nodif_brain_mask_slice_$slicezp -b $subjdir/bvals -r $subjdir/bvecs --forcedir --logdir=$subjdir.bedpostX/diff_slices/data_slice_$slicezp --nj=1000 --bi=$burnin --bn=0 --se=20 --nfibres=$nfibres --fudge=$fudge > $subjdir.bedpostX/logs/log$slicezp
+
touch ${subjdir}.bedpostX/logs/.${slicezp}_finished
fi
fi
slice=`echo "$slice + 1" | bc`
done
-[ "$3" = "" ] || rm ${3}/`hostname`_fdt_${$}
-sleep 10
\ No newline at end of file
+[ "$piddir" = "" ] || rm ${piddir}/`hostname`_fdt_${$}
+sleep 10
diff --git a/xfibres.cc b/xfibres.cc
index 1f28291dc700a7deca74183b112cf46602b8f4f0..0baed21bf997712b0eee145e53f33692ec0f6d0d 100644
--- a/xfibres.cc
+++ b/xfibres.cc
@@ -123,7 +123,12 @@ class Samples{
Matrix m_dsamples;
Matrix m_S0samples;
Matrix m_lik_energy;
-
+
+// // storing signal
+// Matrix m_mean_sig;
+// Matrix m_std_sig;
+// Matrix m_sig2;
+
vector<Matrix> m_thsamples;
vector<Matrix> m_phsamples;
vector<Matrix> m_fsamples;
@@ -135,7 +140,7 @@ class Samples{
vector<Matrix> m_dyadic_vectors;
vector<RowVector> m_mean_fsamples;
vector<RowVector> m_mean_lamsamples;
-
+
float m_sum_d;
float m_sum_S0;
vector<SymmetricMatrix> m_dyad;
@@ -151,7 +156,7 @@ class Samples{
public:
- Samples( volume<int> vol2matrixkey,Matrix matrix2volkey,int nvoxels):
+ Samples( volume<int> vol2matrixkey,Matrix matrix2volkey,int nvoxels,int nmeasures):
opts(xfibresOptions::getInstance()),m_vol2matrixkey(vol2matrixkey),m_matrix2volkey(matrix2volkey){
m_beenhere=m_vol2matrixkey*0;
@@ -172,6 +177,13 @@ public:
m_S0samples=0;
m_lik_energy.ReSize(nsamples,nvoxels);
+ // m_mean_sig.ReSize(nmeasures,nvoxels);
+// m_mean_sig=0;
+// m_std_sig.ReSize(nmeasures,nvoxels);
+// m_std_sig=0;
+// m_sig2.ReSize(nmeasures,nvoxels);
+// m_sig2=0;
+
m_mean_dsamples.ReSize(nvoxels);
m_mean_dsamples=0;
m_mean_S0samples.ReSize(nvoxels);
@@ -222,6 +234,12 @@ public:
m_sum_f[f]+=mfib.fibres()[f].get_f();
m_sum_lam[f]+=mfib.fibres()[f].get_lam();
}
+
+// for(int i=1;i<=m_sig2.Nrows();i++){
+// float sig=mfib.get_signal()(i);
+// m_mean_sig(i,vox)+=sig;
+// m_sig2(i,vox)+=(sig*sig);
+// }
}
void finish_voxel(int vox){
@@ -314,6 +332,7 @@ public:
cart2sph(vec,th,ph);
if(f==0)
mfibre.addfibre(th,ph,m_mean_fsamples[f](voxbest),1,false);//no a.r.d. on first fibre
+ //mfibre.addfibre(th,ph,m_mean_fsamples[f](voxbest),1,true);//a.r.d. on first fibre
else
mfibre.addfibre(th,ph,m_mean_fsamples[f](voxbest),1,true);
@@ -347,7 +366,7 @@ public:
save_volume4D(tmp,logger.appendDir("S0samples"));
tmp.setmatrix(m_lik_energy,mask);
save_volume4D(tmp,logger.appendDir("lik_energy"));
-
+
//Sort the output based on mean_fsamples
//
vector<Matrix> sumf;
@@ -511,7 +530,8 @@ class xfibresVoxelManager{
m_multifibre.set_d(D);
m_multifibre.set_S0(S0);
if(opts.nfibres.value()>0){
- m_multifibre.addfibre(th,ph,f,1,false);//no a.r.d. on first fibre
+ //m_multifibre.addfibre(th,ph,f,1,false);//no a.r.d. on first fibre
+ m_multifibre.addfibre(th,ph,f,1,true);// a.r.d. on first fibre
for(int i=2; i<=opts.nfibres.value(); i++){
m_multifibre.addfibre();
}
@@ -566,6 +586,8 @@ class xfibresVoxelManager{
int main(int argc, char *argv[])
{
try{
+ cout<<"running SAAD'S version of xfibre"<<endl;
+
// Setup logging:
Log& logger = LogSingleton::getInstance();
xfibresOptions& opts = xfibresOptions::getInstance();
@@ -577,6 +599,7 @@ int main(int argc, char *argv[])
bvals=read_ascii_matrix(opts.bvalsfile.value());
bvecs=read_ascii_matrix(opts.bvecsfile.value());
if(bvecs.Nrows()>3) bvecs=bvecs.t();
+ if(bvals.Nrows()>1) bvals=bvals.t();
for(int i=1;i<=bvecs.Ncols();i++){
float tmpsum=sqrt(bvecs(1,i)*bvecs(1,i)+bvecs(2,i)*bvecs(2,i)+bvecs(3,i)*bvecs(3,i));
if(tmpsum!=0){
@@ -600,7 +623,7 @@ int main(int argc, char *argv[])
ColumnVector alpha, beta;
Amat=form_Amat(bvecs,bvals);
cart2sph(bvecs,alpha,beta);
- Samples samples(vol2matrixkey,matrix2volkey,datam.Ncols());
+ Samples samples(vol2matrixkey,matrix2volkey,datam.Ncols(),datam.Nrows());
for(int vox=1;vox<=datam.Ncols();vox++){
cout <<vox<<"/"<<datam.Ncols()<<endl;