diff --git a/fsl/transform/nonlinear.py b/fsl/transform/nonlinear.py
index daf13adcf35f85995f1b309acf38dc2d6c145d37..77bc212b87d0fc71dfd3f4155fd5ed36a0d12f95 100644
--- a/fsl/transform/nonlinear.py
+++ b/fsl/transform/nonlinear.py
@@ -251,7 +251,9 @@ class DeformationField(NonLinearTransform):
         if from_ is None: from_ = self.refSpace
         if to    is None: to    = self.srcSpace
 
-        coords = np.asanyarray(coords)
+        coords   = np.asanyarray(coords)
+        outshape = coords.shape
+        coords   = coords.reshape((-1, 3))
 
         # We may need to pre-transform the
         # coordinates so they are in the
@@ -299,7 +301,7 @@ class DeformationField(NonLinearTransform):
         outcoords          = np.full(coords.shape, np.nan)
         outcoords[voxmask] = disps
 
-        return outcoords
+        return outcoords.reshape(outshape)
 
 
 class CoefficientField(NonLinearTransform):