diff --git a/bob/pad/face/test/test_block.py b/bob/pad/face/test/test_block.py
new file mode 100644
index 0000000000000000000000000000000000000000..6c23374e5e532230e497e3b25a1040a643229ac3
--- /dev/null
+++ b/bob/pad/face/test/test_block.py
@@ -0,0 +1,22 @@
+import numpy
+
+A_org = numpy.array(range(1, 17), "float64").reshape((4, 4))
+A_ans_0_3D = numpy.array(
+    [[[1, 2], [5, 6]], [[3, 4], [7, 8]], [[9, 10], [13, 14]], [[11, 12], [15, 16]]],
+    "float64",
+)
+A_ans_0_4D = numpy.array(
+    [[[[1, 2], [5, 6]], [[3, 4], [7, 8]]], [[[9, 10], [13, 14]], [[11, 12], [15, 16]]]],
+    "float64",
+)
+
+from bob.pad.face.utils.load_utils import block
+
+
+def test_block():
+    B = block(A_org, (2, 2), (0, 0))
+
+    assert (B == A_ans_0_4D).all()
+
+    B = block(A_org, (2, 2), (0, 0), flat=True)
+    assert (B == A_ans_0_3D).all()
diff --git a/bob/pad/face/utils/load_utils.py b/bob/pad/face/utils/load_utils.py
index 353df179604adc229c8aedbf39e67b241a67b2f6..04f82a8b76f320efcc905faab5e284da49743bfd 100644
--- a/bob/pad/face/utils/load_utils.py
+++ b/bob/pad/face/utils/load_utils.py
@@ -11,6 +11,46 @@ from imageio import get_reader
 from PIL import Image
 
 
+def block(X, block_size, block_overlap, flat=False):
+    """
+    Parameters
+    ----------
+    X : numpy.ndarray
+        The image to be split into blocks.
+    block_size : tuple
+        The size of the block.
+    block_overlap : tuple
+        The overlap of the block.
+
+    Returns
+    -------
+    numpy.ndarray
+        The image split into blocks.
+    """
+
+    assert len(block_size) == 2
+    assert len(block_overlap) == 2
+
+    size_ov_h = int(block_size[0] - block_overlap[0])
+    size_ov_w = int(block_size[1] - block_overlap[1])
+    n_blocks_h = int((X.shape[0] - block_overlap[0]) / size_ov_h)
+    n_blocks_w = int((X.shape[1] - block_overlap[1]) / size_ov_w)
+
+    blocks = numpy.zeros(shape=(n_blocks_h, n_blocks_w, size_ov_h, size_ov_w))
+    for h in range(n_blocks_h):
+        for w in range(n_blocks_w):
+
+            blocks[h, w, :, :] = X[
+                h * size_ov_h : h * size_ov_h + block_size[0],
+                w * size_ov_w : w * size_ov_w + block_size[1],
+            ]
+
+    if flat:
+        return blocks.reshape(n_blocks_h * n_blocks_w, blocks.shape[2], blocks.shape[3])
+
+    return blocks
+
+
 def scale(image, scaling_factor):
     """
     Scales and image using PIL
@@ -199,19 +239,21 @@ def blocks(data, block_size, block_overlap=(0, 0)):
     ValueError
         If data dimension is not between 2 and 4 (inclusive).
     """
+
     data = numpy.asarray(data)
     # if a gray scale image:
     if data.ndim == 2:
         output = block(data, block_size, block_overlap, flat=True)
     # if a color image:
     elif data.ndim == 3:
-        out_shape = list(data.shape[0:1]) + list(
-            block_output_shape(data[0], block_size, block_overlap, flat=True)
-        )
+        # out_shape = list(data.shape[0:1]) + list(
+        #    block_output_shape(data[0], block_size, block_overlap, flat=True)
+        # )
 
-        output = numpy.empty(out_shape, dtype=data.dtype)
+        # output = numpy.empty(out_shape, dtype=data.dtype)
+        output = []
         for i, img2d in enumerate(data):
-            block(img2d, block_size, block_overlap, output[i], flat=True)
+            output.append(block(img2d, block_size, block_overlap, flat=True))
         output = numpy.moveaxis(output, 0, 1)
     # if a color video:
     elif data.ndim == 4: