Alignment of two sets of points on the 2D plane

PointsAlignment.py

def pointAlignment(ref, move):
    """
    Mahd Sadjadi (2017)
    This is a function which takes two sets of
    points each with N points (both inputs are assumd to be
    numpy arrays of N*2 shape). The first input is 
    the reference grid and the second is the point set to 
    be aligned with the reference.
    This function shifts and rotate points to
    make the two sets aligned and returns the 
    result (with the shape N*2).
    """
    N = ref.shape[0] # number of atoms
    IM = np.ones(N) #array of ones
    XG,YG,ZG = ref.T # ref coordinates
    XM,YM,ZM = move.T # move coordinates
    XGA,YGA,ZGA = np.mean(ref,axis=0) # center of mass
    XMA,YMA,ZMA = np.mean(move,axis=0) # center of mass
    XGS = XG - XGA * IM # std xg
    YGS = YG - YGA * IM # std yg
    XMS = XM - XMA * IM
    YMS = YM - YMA * IM
    A = np.sum(XGS * XMS + YGS * YMS)/N
    B = np.sum(XGS * YMS - YGS * XMS)/N
    R = np.sqrt(A**2 + B**2)
    sin = B/R
    cos = A/R
    a = XGA * cos - YGA * sin - XMA
    b = XGA * sin + YGA * cos - YMA
    theta = np.arctan(B/A)
    
    deltax = (a * cos + b * sin)*IM - (1-cos)*XM + sin*YM
    deltay = (-a * sin + b * cos)*IM - sin*XM -(1-cos)*YM
    results = np.column_stack((XM+deltax,YM+deltay,ZM))
    return results