diff --git a/oxford_asl.in b/oxford_asl.in
index 6ffa012..47d537f 100755
--- a/oxford_asl.in
+++ b/oxford_asl.in
@@ -646,7 +646,7 @@ ApplyWarps() {
         Log "   - Moco: $moco_xfm"
     	appremat="--premat=$moco_xfm"
     fi
-    applywarp -i $tempdir/asldata_orig -r $tempdir/meanasl -o $tempdir/asldata $appremat -w $corrwarp --rel $APPLYWARP_OPTS
+    applywarp -i $tempdir/asldata_orig -r $tempdir/nativeref -o $tempdir/asldata $appremat -w $corrwarp --rel $APPLYWARP_OPTS
 
     # Correct for quantitative signal magnitude as volume has been locally scaled by the Jacobian
     fslmaths $tempdir/asldata -mul $jacobian $tempdir/asldata
@@ -663,14 +663,15 @@ ApplyWarps() {
             Log "   - Calib->ASL: $calib_xfm"
             appremat="--premat=$calib_xfm"
         fi
-        applywarp -i $tempdir/calib_orig -r $tempdir/meanasl -o $tempdir/calib $appremat -w $corrwarp --rel $APPLYWARP_OPTS
+
+        applywarp -i $tempdir/calib_orig -r $tempdir/nativeref -o $tempdir/calib $appremat -w $corrwarp --rel $APPLYWARP_OPTS
         fslmaths $tempdir/calib -mul $jacobian $tempdir/calib
         if [ ! -z $cref ]; then
-            applywarp -i $tempdir/cref_orig -r $tempdir/meanasl -o $tempdir/cref $appremat -w $corrwarp --rel $APPLYWARP_OPTS
+            applywarp -i $tempdir/cref_orig -r $tempdir/nativeref -o $tempdir/cref $appremat -w $corrwarp --rel $APPLYWARP_OPTS
             fslmaths $tempdir/cref -mul $jacobian $tempdir/cref
         fi
         if [ ! -z $cblip ]; then
-            applywarp -i $tempdir/cblip_orig -r $tempdir/meanasl -o $tempdir/cblip $appremat -w $corrwarp --rel $APPLYWARP_OPTS
+            applywarp -i $tempdir/cblip_orig -r $tempdir/nativeref -o $tempdir/cblip $appremat -w $corrwarp --rel $APPLYWARP_OPTS
             fslmaths $tempdir/cblip -mul $jacobian $tempdir/cblip
         fi
     fi
@@ -1194,6 +1195,14 @@ imcp $infile $tempdir/asldata # this is the MAIN data that we will reflect any c
 # take a copy that will not be subject to any subsequent corrections
 imcp $tempdir/asldata $tempdir/asldata_orig
 
+# Take the mean of the ASL data to use as a reference image in applywarp. 
+# This is required because in FSL6, applywarp duplicates the output for every
+# volume of the reference image, so it is imperative that it is a 3D image!
+# Note that in applywarp the only function of the reference image is to define
+# the output space, so this never needs to be changed as our 'native' space is
+# always the same space as the input ASL data
+fslmaths $tempdir/asldata -Tmean $tempdir/nativeref
+
 ### Motion Correction (main)
 # note motion correction within calibration data is done above
 if [ ! -z $moco ]; then
@@ -1214,7 +1223,7 @@ if [ ! -z $moco ]; then
 	done
     # Concatenate the motion matrices for use with applywarp
     cat $tempdir/asldata.mat/MAT* > $tempdir/asldata.cat
-    applywarp --ref=$tempdir/asldata --in=$tempdir/asldata_orig --out=$tempdir/asldata --premat=$tempdir/asldata.cat $APPLYWARP_OPTS
+    applywarp --ref=$tempdir/nativeref --in=$tempdir/asldata_orig --out=$tempdir/asldata --premat=$tempdir/asldata.cat $APPLYWARP_OPTS
 	# Convert all calibration images to align with asldata
 	flirt -in $tempdir/calib -out $tempdir/calib -ref $tempdir/asldata -init $tempdir/calib2asl.mat -applyxfm $APPLYXFM_OPTS
     if [ ! -z $cref ]; then
@@ -1244,7 +1253,7 @@ if [ ! -z $gdcwarp ]; then
 
     # Convert the supplied warp image to extract the jacbbian (as parts) 
     # (meanasl here is just a reference). Using relative warp convention
-    convertwarp -r $tempdir/meanasl -o $tempdir/asldist_userwarp -w $gdcwarp --rel -j $tempdir/distcorr/jacobian_parts
+    convertwarp -r $tempdir/nativeref -o $tempdir/asldist_userwarp -w $gdcwarp --rel -j $tempdir/distcorr/jacobian_parts
     
     # Calculation of the jacobian for the warp - method suggested in:
     # https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=FSL;d3fee1e5.0908
@@ -1475,7 +1484,7 @@ if [ ! -z $fasthasrun ]; then
 
     if [ ! -z $tempdir/biasfield_struct ]; then
 	# transform the bias field and invert to use for sensitivity correction in calibration
-	applywarp --ref=$tempdir/asldata --in=$tempdir/biasfield_struct --out=$tempdir/biasfield --premat=$tempdir/struct2asl.mat $APPLYWARP_OPTS
+	applywarp --ref=$tempdir/nativeref --in=$tempdir/biasfield_struct --out=$tempdir/biasfield --premat=$tempdir/struct2asl.mat $APPLYWARP_OPTS
     fi
 
     if [ ! -z $senscorr ]; then
@@ -1618,7 +1627,7 @@ if [ -e $tempdir/fmap.* ]; then
     Registration $distbase "-m $tempdir/mask --tissseg $tempdir/tissseg --imat $tempdir/asl2struct.mat --finalonly" distcorr
     
     # Generate the correction warp in ASL space
-    convertwarp -r $tempdir/meanasl -o $tempdir/asldist_warp -w $tempdir/distcorr/asl2struct_warp.nii.gz --postmat=$tempdir/distcorr/struct2asl.mat --rel -j $tempdir/distcorr/jacobian_parts
+    convertwarp -r $tempdir/nativeref -o $tempdir/asldist_warp -w $tempdir/distcorr/asl2struct_warp.nii.gz --postmat=$tempdir/distcorr/struct2asl.mat --rel -j $tempdir/distcorr/jacobian_parts
 
     # If we ALSO have a user supplied GDC warp then we need to merge that with the one we have found here
     # use convert warp to combine the warps. Note that we cannot do this in a single step becase we need
@@ -1626,7 +1635,7 @@ if [ -e $tempdir/fmap.* ]; then
     # the user warp 
     if [ ! -z $gdcwarp ]; then
 	    Log "Adding user-supplied GDC warp to distortion correction warp"
-        convertwarp -r $tempdir/meanasl -o $tempdir/asldist_warp -w $tempdir/asldist_warp --warp2=$tempdir/asldist_userwarp --rel -j $tempdir/distcorr/jacobian_parts
+        convertwarp -r $tempdir/nativeref -o $tempdir/asldist_warp -w $tempdir/asldist_warp --warp2=$tempdir/asldist_userwarp --rel -j $tempdir/distcorr/jacobian_parts
     fi
 
     # Calculation of the jacobian for the warp - method suggested in:
@@ -1982,9 +1991,9 @@ if [ ! -z $fasthasrun ] && [ -z $pvgm ]; then
     convert_xfm -omat $tempdir/struct2asl.mat -inverse $tempdir/asl2struct.mat
     
     # Gray matter - assume this will be PVE 1
-    applywarp --ref=$tempdir/asldata --in=$tempdir/pvgm_struct --out=$tempdir/pvgm_inasl --premat=$tempdir/struct2asl.mat $APPLYWARP_OPTS
+    applywarp --ref=$tempdir/nativeref --in=$tempdir/pvgm_struct --out=$tempdir/pvgm_inasl --premat=$tempdir/struct2asl.mat $APPLYWARP_OPTS
     # white matter  - assume this will be PVE 2
-    applywarp --ref=$tempdir/asldata --in=$tempdir/pvwm_struct --out=$tempdir/pvwm_inasl --premat=$tempdir/struct2asl.mat $APPLYWARP_OPTS
+    applywarp --ref=$tempdir/nativeref --in=$tempdir/pvwm_struct --out=$tempdir/pvwm_inasl --premat=$tempdir/struct2asl.mat $APPLYWARP_OPTS
     # threshold (upper and lower) the PVE to avoid artefacts of spline interpolation and also ignore very low PVE that could cause numerical issues.
     fslmaths $tempdir/pvgm_inasl -thr 0.1 -min 1 $tempdir/pvgm_inasl
     fslmaths $tempdir/pvwm_inasl -thr 0.1 -min 1 $tempdir/pvwm_inasl
@@ -2138,7 +2147,7 @@ if [ ! -z $pvcorr ]; then
         if [ ! -z $senscorr ]; then
 	      # Struct->ASL registration may have changed, so re-transform the bias field and invert to use for sensitivity correction
           Log "Re-registering bias field and calculating sensitivity map"
-          applywarp --ref=$tempdir/asldata --in=$tempdir/biasfield_struct --out=$outdir/native_space/biasfield --premat=$tempdir/struct2asl.mat $APPLYWARP_OPTS
+          applywarp --ref=$tempdir/nativeref --in=$tempdir/biasfield_struct --out=$outdir/native_space/biasfield --premat=$tempdir/struct2asl.mat $APPLYWARP_OPTS
 	      fslmaths $outdir/native_space/biasfield -recip $outdir/native_space/sensitivity
         fi
       fi