gistfile1.py

import numpy as np
import matplotlib.pyplot as plt

def pseudoinverse (A):
    U, S, V = np.linalg.svd(A)
    S = np.diag(1/S)
    S_z=np.zeros((V[1].size,U[1].size))
    S_z[:S[0].size,:S[1].size] = S
    return V.T.dot(S_z.dot(U.T))

fig = plt.figure(figsize=(18, 9))
subp1 = fig.add_subplot(121)

cloud = np.load("shape003.npy")
colorLetters = ('r', 'g', 'b', 'y')

pointsNumber = cloud.shape[1]
print(pointsNumber)

A = np.vstack((cloud[0,:,0][newaxis].T, cloud[1,:,1][newaxis].T))
A = np.hstack((A, np.vstack((np.ones((pointsNumber,1)), np.zeros((pointsNumber,1))))))
A = np.hstack((A, np.vstack((np.zeros((pointsNumber,1)), np.ones((pointsNumber,1))))))
b = np.vstack((cloud[0,:,1][newaxis].T, -cloud[1,:,0][newaxis].T))

A_inv=pseudoinverse(A.T.dot(A))
x = A_inv.dot(A.T.dot(b))

print(x)

for c, i in enumerate(cloud):
    subp1.scatter(i[:,0], i[:,1], c=colorLetters[c]) 
    
x1, x2 = 20,100

k1 = x[0]
k2 = -1/k1
b1 = x[1]
b2 = -x[2]/k1

plot([x1, x2], [k1*x1+b1, k1*x2+b1], color=colorLetters[0], linewidth=2)
plot([x1, x2], [k2*x1+b2, k2*x2+b2], color=colorLetters[1], linewidth=2)

subp1.set_xlim((-100, 250))
subp1.set_ylim((-250, 100))