diff --git a/mdevaluate/pbc.py b/mdevaluate/pbc.py
index 90165f7..fb6b50f 100644
--- a/mdevaluate/pbc.py
+++ b/mdevaluate/pbc.py
@@ -198,36 +198,28 @@ def nojump(frame, usecache=True):
             [m[:frame.step + 1, selection].sum(axis=0) for m in reader.nojump_matrixes]
         ).T) * frame.box.diagonal()
     return frame - delta
-
-
-def pbc_points(coordinates, box, thickness=0, index=False, inclusive=True, center=None):
-    """
-    Returns the points their first periodic images. Does not fold them back into the box.
-    Thickness 0 means all 27 boxes. Positive means the box+thickness. Negative values mean that less than the box is returned.
-    index=True also returns the indices with indices of images being their originals values.
-    inclusive=False returns only images, does not work with thickness <= 0
-    """
-    if center is None:
-        center = box/2
-    allcoordinates = np.copy(coordinates)
-    indices = np.tile(np.arange(len(coordinates)),(27))
-    for x in range(-1, 2, 1):
-            for y in range(-1, 2, 1):
-                for z in range(-1, 2, 1):
-                    vv = np.array([x, y, z], dtype=float)
-                    if not (vv == 0).all() :
-                        allcoordinates = np.concatenate((allcoordinates, coordinates + vv*box), axis=0)
     
+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.
+    Thickness 0 means all 27 boxes. Positive means the box+thickness. 
+    Negative values mean that less than the box is returned.
+    index=True also returns the indices with indices of images being their originals values.
+    """
+    if shear:
+        box[2,0] = box[2,0]%box[0,0]
+    grid = np.array([ [i, j, k] for k in [-1,0,1] for j in [1,0,-1] for i in [-1,0,1] ])
+    coordinatesPBC = np.concatenate([coordinates+v@box for v in grid], axis=0)
+    size = np.diag(box)
+    indices = np.tile(np.arange(len(coordinates)),(27))
     if thickness != 0:
-        mask = np.all(allcoordinates < center+box/2+thickness, axis=1)
-        allcoordinates = allcoordinates[mask]
+        mask = np.all(coordinatesPBC > -thickness, axis=1)
+        coordinatesPBC = coordinatesPBC[mask]
         indices = indices[mask]
-        mask = np.all(allcoordinates > center-box/2-thickness, axis=1)
-        allcoordinates = allcoordinates[mask]
+        mask = np.all(coordinatesPBC < size+thickness, axis=1)
+        coordinatesPBC = coordinatesPBC[mask]
         indices = indices[mask]
-    if not inclusive and thickness > 0:
-        allcoordinates = allcoordinates[len(coordinates):]
-        indices = indices[len(coordinates):]
     if index:
-        return (allcoordinates, indices)
-    return allcoordinates
+        return coordinatesPBC, indices
+    return coordinatesPBC