In [1]:

    

import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline


    

In [209]:

    

d1 = np.array([[0,0], [0,3], [4,0], [0,0]])

theta = np.pi * 2 / 3

Rot = np.matrix([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])
print(Rot)
translation = np.matrix([[5, 4]])

def Transform(p1, Rot, translation):
    p2 = Rot * p1.T 
    p2 = p2.T + translation
    return np.asarray(p2)
    
d2 = Transform(d1, Rot, translation) + \
    np.random.normal(0, 0.1, size=2*d1.shape[0]).reshape(d1.shape)

d2


    

    




[[-0.5       -0.8660254]
 [ 0.8660254 -0.5      ]]






    
Out[209]:





array([[ 5.10001217,  3.80883984],
       [ 2.51215817,  2.54176546],
       [ 2.97742202,  7.33203299],
       [ 4.90784786,  4.00401863]])

In [210]:

    

plt.figure()
plt.axis('equal')

plt.plot(d1[:,0], d1[:,1], 'bo')
plt.plot(d1[:,0], d1[:,1], 'b')
plt.fill(d1[:,0], d1[:,1], 'b')

plt.plot(d2[:,0], d2[:,1], 'g^')
plt.plot(d2[:,0], d2[:,1], 'g')
plt.fill(d2[:,0], d2[:,1], 'g')
plt.show()


    

In [211]:

    

d1 = np.array([[-2,2], [-1.5,2.4], [-1,2], [-0.5,2],[0,2.4], 
               [0.5,2], [1,1], [1.5,1], [1.5,0.5], [1,0.5], [1,1]])

d2 = Transform(d1, Rot, translation) + \
    np.random.normal(0, 0.1, size=2*d1.shape[0]).reshape(d1.shape)

def plotshapes(d1, d2, d3=None, d4=None):
    plt.figure()
    plt.axis('equal')
    plt.plot(d1[:,0], d1[:,1], 'bo')
    plt.plot(d2[:,0], d2[:,1], 'g^')
    plt.plot(d1[:,0], d1[:,1], 'b')
    plt.plot(d2[:,0], d2[:,1], 'g')
    if d3 is not None:
        plt.plot(d3[:,0], d3[:,1], 'co')
        plt.plot(d3[:,0], d3[:,1], 'c--')
    if d4 is not None:
        plt.plot(d4[:,0], d4[:,1], 'm^')
        plt.plot(d4[:,0], d4[:,1], 'm--')
    
    plt.show()
    
plotshapes(d1, d2)


    

In [203]:

    

def simTransform(pref, pcmp, showerror = False):
    err_before = np.mean(np.sum(np.square(pref - pcmp), axis=1))
    ref_mean = np.mean(pref, axis=0)
    prefcentered = np.asmatrix(pref) - np.asmatrix(ref_mean)
    
    cmp_mean = np.mean(pcmp, axis=0)
    pcmpcentered = np.asmatrix(pcmp) - np.asmatrix(cmp_mean)   
    
    Sxx = np.sum(np.square(pcmpcentered[:,0]))
    Syy = np.sum(np.square(pcmpcentered[:,1]))
    Sxxr = prefcentered[:,0].T * pcmpcentered[:,0] #(ref_x, x)
    Syyr = prefcentered[:,1].T * pcmpcentered[:,1] #(ref_y, y)
    Sxyr = prefcentered[:,1].T * pcmpcentered[:,0] #(ref_y, x)
    Syxr = prefcentered[:,0].T * pcmpcentered[:,1] #(ref_x, y)
    a = (Sxxr + Syyr)/(Sxx + Syy) #(Sxxr + Syyr) / (Sxx + Syy)
    b = (Sxyr - Syxr) / (Sxx + Syy)
    a = np.asscalar(a)
    b = np.asscalar(b)
    Rot = np.matrix([[a, -b],[b, a]])
    translation = -Rot * np.asmatrix(cmp_mean).T + np.asmatrix(ref_mean).T
    outx, outy = [], []
    res = Rot * np.asmatrix(pcmp).T + translation
    
    if showerror:
        err_after = np.mean(np.sum(np.square(pref - res.T), axis=1))
        print("Error before: %.4f    after: %.4f\n"%(err_before, err_after))
    return res.T


    

In [212]:

    

d3 = simTransform(pref=d1, pcmp=d2, showerror=True)
print(d3)
plotshapes(d1, d2, d4=d3)


    

    




Error before: 20.4283    after: 0.0127

[[-1.85291902  2.01625591]
 [-1.47396707  2.50134012]
 [-1.07379556  2.02663372]
 [-0.51026754  1.92930963]
 [-0.08665807  2.36771163]
 [ 0.42053322  2.03005122]
 [ 0.91639739  0.90297445]
 [ 1.50053308  1.05409419]
 [ 1.5455354   0.66182626]
 [ 1.00007022  0.42760748]
 [ 1.11453796  0.8821954 ]]






    

In [1]:

    

def simTransform1(ref_x, ref_y, x, y, reportRMSD = False):
    rmsd_before = np.sqrt(np.sum(np.square(ref_x - x) + np.square(ref_y - y)))
    err_before = np.mean(np.square(ref_x - x) + np.square(ref_y - y))
    ref_x = ref_x - np.mean(ref_x)
    ref_y = ref_y - np.mean(ref_y)
    x = x - np.mean(x)
    y = y - np.mean(y)
    Sxx = np.sum(np.square(x))
    Syy = np.sum(np.square(y))
    Sxxr = np.dot(ref_x, x)
    Syyr = np.dot(ref_y, y)
    Sxyr = np.dot(ref_y, x)
    Syxr = np.dot(ref_x, y)
    a = np.sum(np.array((ref_x, ref_y)).T * np.array((x,y)).T)/(Sxx + Syy) #(Sxxr + Syyr) / (Sxx + Syy)
    b = (Sxyr - Syxr) / (Sxx + Syy)
    #scale_factor = 1.0/np.sqrt(a**2 + b**2)
    Mat = np.matrix([[a, -b],[b, a]])
    outx, outy = [], []
    for i,(px,py) in enumerate(zip(x,y)):
        v = Mat * np.matrix([px, py]).T #np.multiply(Mat, np.array((px,py)).T)
        outx.append(v[0,0])
        outy.append(v[1,0])
        
    outx = np.array(outx)
    outy = np.array(outy)
    rmsd_after = np.sqrt(np.sum(np.square(ref_x - outx) + np.square(ref_y - outy)))
    if reportRMSD:
        err_after = np.mean(np.square(ref_x - outx) + np.square(ref_y - outy))
        print("RMSD before: %.4f  after: %.4f\n"%(rmsd_before, rmsd_after))
        print("Error before: %.4f  after: %.4f\n"%(err_before, err_after))
    return(outx, outy)


    

In [2]:

    

def simTransform2(ref_x, ref_y, x, y, reportRMSD = False):
    err_before = np.mean(np.square(ref_x - x) + np.square(ref_y - y))
    n = x.shape[0]
    A = np.vstack((x, y)).T
    B = np.vstack((ref_x, ref_y)).T
    A = np.matrix(A)
    B = np.matrix(B)
    rmsd_before = np.sqrt(np.sum(np.square(ref_x - x) + np.square(ref_y - y)))
    centroid_A = np.mean(A, axis=0)
    centroid_B = np.mean(B, axis=0)
    
    AA = A - np.tile(centroid_A, (n,1))
    BB = B - np.tile(centroid_B, (n,1))
    
    H = np.transpose(AA) * BB
    
    U, s, Vt = np.linalg.svd(H)

    R = Vt.T * U.T #np.matrix([[a, -b],[b, a]])
    if np.linalg.det(R) < 0:
        #sys.stderr.write("Warning: reflection detected!")
        Vt[1,:] *= -1
        R = Vt.T * U.T
        
    tr = -R * centroid_A.T + centroid_B.T
    A2 = R * A.T + np.tile(tr, (1,n))
    A2 = A2.T
    outx = np.asarray(A2[:,0]).flatten()
    outy = np.asarray(A2[:,1]).flatten()
    if reportRMSD:
        err_after = np.mean(np.square(ref_x - outx) + np.square(ref_y - outy))
        print("Error before: %.4f  after: %.4f\n"%(err_before, err_after))

    
    return (outx, outy)


    

In [ ]: