diff --git a/mdevaluate/pbc.py b/mdevaluate/pbc.py
index 7f0690c..99ee9f0 100644
--- a/mdevaluate/pbc.py
+++ b/mdevaluate/pbc.py
@@ -241,9 +241,7 @@ def nojump(frame, usecache=True):
     return frame - delta
 
 
-def pbc_points(
-    coordinates, box, thickness=0, index=False, shear=False, extra_image=False
-):
+def pbc_points(coordinates, box, thickness=0, index=False, shear=False):
     """
     Returns the points their first periodic images. Does not fold
     them back into the box.
@@ -253,22 +251,14 @@ def pbc_points(
     originals values.
     """
     if shear:
-        box[2, 0] = box[2, 0] % box[0, 0]
-    if shear or extra_image:
-        grid = np.array(
-            [
-                [i, j, k]
-                for k in [-2, -1, 0, 1, 2]
-                for j in [2, 1, 0, -1, -2]
-                for i in [-2, -1, 0, 1, 2]
-            ]
-        )
-        indices = np.tile(np.arange(len(coordinates)), 125)
-    else:
-        grid = np.array(
-            [[i, j, k] for k in [-1, 0, 1] for j in [1, 0, -1] for i in [-1, 0, 1]]
-        )
-        indices = np.tile(np.arange(len(coordinates)), 27)
+        box[2, 0] = box[2, 0] % box[2, 2]
+        if box[2, 0] > box[2, 2]:
+            box[2, 0] = box[2, 0] - box[2, 2]
+
+    grid = np.array(
+        [[i, j, k] for k in [-1, 0, 1] for j in [1, 0, -1] for i in [-1, 0, 1]]
+    )
+    indices = np.tile(np.arange(len(coordinates)), 27)
     coordinates_pbc = np.concatenate([coordinates + v @ box for v in grid], axis=0)
     size = np.diag(box)