diff --git a/poisson/continuous/shapes.py b/poisson/continuous/shapes.py
index 73bb5a71aa7533412ea94d71069534954421a6c1..52a676132efd660314512cf6890a5bccf5e890da 100644
--- a/poisson/continuous/shapes.py
+++ b/poisson/continuous/shapes.py
@@ -866,7 +866,6 @@ def rotate(x, angle, axis=None, center='self'):
         if None, will be the origin according to the dimension
         if self, will be the average of the points
         (center of regular polygon for 'inplace' rotation)
-
     '''
     x = np.asarray(x)
     center_parent = center
@@ -877,7 +876,7 @@ def rotate(x, angle, axis=None, center='self'):
     if len(x.shape) == 2:
         npt, ndim = x.shape
     elif len(x.shape) == 1:
-        npt, ndim = 1, x.shape[0]
+        npt, ndim = 0, x.shape[0]
         print(x.shape)
         x = x[None, :]
     else: raise ValueError
@@ -904,7 +903,60 @@ def rotate(x, angle, axis=None, center='self'):
             np.repeat([rot_mat], len(x), axis=0), x - center_mat)
            + center_mat)
 
-    if npt is 1:
+    if npt is 0:
+        return ret[0]
+    else:
+        return ret
+
+def reflect(x, axis):
+    '''returns reflections of point(s) x wrt to axis
+
+    Parameters:
+    -----------
+    x: array like
+        point or list of points in 2 or 3d
+    axis: array like
+        2 (3) points that generate the line (plan)
+        wrt which to reflect in 2d (3d)
+
+    Returns:
+    --------
+    ret: array like (shape of x)
+        the reflected point(s)
+    '''
+    x = np.asarray(x)
+    if len(x.shape) == 1:
+        ndim, npt = x.shape[0], 1
+        x = x[None, :]
+    else:
+        npt, ndim = x.shape
+    assert (ndim == 2 or ndim == 3)
+
+    assert isinstance(axis, (list, np.ndarray))
+    axis = np.asarray(axis)
+    assert axis.shape[1] == ndim
+    A = axis[0]
+    Amat = np.repeat([A], len(x), axis=0)
+
+    if ndim == 2:
+        u = (axis[1:] - A)[0]
+        v = np.array([u[1], -u[0]])
+        # u, v is othogonal basis
+        mirror = np.array([[1, 0],[0, -1]])
+        Mat = np.array([u, v])
+
+    elif ndim == 3:
+        u, C = axis[1:] - A
+        w = np.cross(u, C)
+        v = np.cross(w, u)
+        # u, v, w is othogonal basis
+        mirror = np.array([[1, 0, 0],[0, 1, 0], [0 ,0 ,-1]])
+        Mat = np.array([u, v, w])
+
+    Matinv = np.linalg.inv(Mat)
+    transf = Mat.T @ mirror @ Matinv.T
+    ret = np.einsum('kij, kj -> ki', np.repeat([transf], len(x), axis=0), (x - Amat)) + Amat
+    if npt is 0:
         return ret[0]
     else:
         return ret
@@ -1210,6 +1262,34 @@ def __test_plot_rotation_translation():
     plt.axes().set_aspect('equal', 'datalim')
     plt.show()
 
+def __test_reflection():
+    '''test if two reflexions is equal to identity'''
+    np.random.seed(112)
+
+    ndim = 2
+    axis = np.random.rand(ndim, ndim)
+    pt = np.random.rand(150, ndim)
+    assert(np.allclose(pt, reflect(reflect(pt, axis), axis)))
+
+    ndim = 3
+    axis = np.random.rand(ndim, ndim)
+    pt = np.random.rand(150, ndim)
+    assert(np.allclose(pt, reflect(reflect(pt, axis), axis)))
+
+    print('test reflexion passed')
+
+def __test_plot_reflection():
+    ndim = 2
+    axis = np.random.rand(ndim, ndim) - 0.5
+    pt = np.random.rand(10, ndim)
+    newpt = reflect(pt, axis)
+
+    plt.figure(figsize=(12,12))
+    plt.plot(*pt.T, marker='o')
+    plt.plot(*axis.T, lw=4)
+    plt.plot(*newpt.T, marker='^')
+    plt.axes().set_aspect('equal')
+    plt.show()
 
 #test_plot_ellipse()
 #test_plot_rectangle()