Следующий код является реализацией DLT (прямое линейное преобразование).
У меня есть проблема в коде:
1. Я накапливаю ошибку, вызванную по шуму, чтобы получить рациональные результаты, без накопления шумовой тренд не выглядит правильным, и после накопления я не думаю, что кривая в любом случае верна, но я не могу понять, почему! Подскажите пожалуйста, в чем проблема в моей реализации гауссовского шума?
import numpy as np
import matplotlib.pyplot as plt
#------------------------------
np.random.seed(10)
def Normalization(nd, x):
x = np.asarray(x)
m, s = np.mean(x, axis=0), np.std(x)
if nd == 2:
Tr = np.array([[s, 0, m[0]], [0, s, m[1]], [0, 0, 1]])
else:
Tr = np.array([[s, 0, 0, m[0]], [0, s, 0, m[1]], [0, 0, s, m[2]], [0, 0, 0, 1]])
Tr = np.linalg.inv(Tr)
x = np.dot( Tr, np.concatenate( (x.T, np.ones((1,x.shape[0]))) ) )
x = x[0:nd, :].T
return Tr, x
def DLTcalib(nd, xyz, uv):
if (nd != 3):
raise ValueError('%dD DLT unsupported.' %(nd))
# Converting all variables to numpy array
xyz = np.asarray(xyz)
uv = np.asarray(uv)
n = xyz.shape[0]
# Validating the parameters:
if uv.shape[0] != n:
raise ValueError('Object (%d points) and image (%d points) have different number of points.' %(n, uv.shape[0]))
if (xyz.shape[1] != 3):
raise ValueError('Incorrect number of coordinates (%d) for %dD DLT (it should be %d).' %(xyz.shape[1],nd,nd))
if (n < 6):
raise ValueError('%dD DLT requires at least %d calibration points. Only %d points were entered.' %(nd, 2*nd, n))
Txyz, xyzn = Normalization(nd, xyz)
Tuv, uvn = Normalization(2, uv)
A = []
for i in range(n):
x, y, z = xyzn[i, 0], xyzn[i, 1], xyzn[i, 2]
u, v = uvn[i, 0], uvn[i, 1]
A.append( [x, y, z, 1, 0, 0, 0, 0, -u * x, -u * y, -u * z, -u] )
A.append( [0, 0, 0, 0, x, y, z, 1, -v * x, -v * y, -v * z, -v] )
# Convert A to array
A = np.asarray(A)
# Find the 11 parameters:
U, S, V = np.linalg.svd(A)
# The parameters are in the last line of Vh and normalize them
L = V[-1, :] / V[-1, -1]
H = L.reshape(3, nd + 1)
H = np.dot( np.dot( np.linalg.pinv(Tuv), H ), Txyz )
H = H / H[-1, -1]
L = H.flatten('C')
uv2 = np.dot( H, np.concatenate( (xyz.T, np.ones((1, xyz.shape[0]))) ) )
uv2 = uv2 / uv2[2, :]
err = np.sqrt( np.mean(np.sum( (uv2[0:2, :].T - uv)**2, 1)) )
return L, err
def DLT():
# Known 3D coordinates
xyz = [[-875, 0, 9.755], [442, 0, 9.755], [1921, 0, 9.755], [2951, 0.5, 9.755],
[-4132, 0.5, 23.618],[-876, 0, 23.618]]
# Known pixel coordinates
uv = [[76, 706], [702, 706], [1440, 706], [1867, 706], [264, 523], [625, 523]]
#------------------------------------------------------------------------------
nd = 3
errs = []
for sigma in range(0,100,1): # 30
# print('len__XYZ',len(xyz))
noise = np.random.normal(0,sigma/10,(6,2))
uv_noise = uv+noise
#print(type(noise),noise.shape)
P, err = DLTcalib(nd, xyz, uv_noise)
if(sigma):
err +=errs[sigma-1]
errs.append(err)
#-------------------------------
print('Matrix')
#print(P)
P = P.reshape(3,4)
H1 = P[:,:3]
H2 = P[:,3]
Q, R = np.linalg.qr(H1)
print('Q=',Q,'\n----------\nR',R)
K = np.linalg.inv(R)
K = K / K[-1,-1]
Rot = np.linalg.inv(Q)
C = np.dot(np.linalg.inv(H1),H2)
t = - np.dot(Rot, C)
print('\n----------\nK =',K,'\n----------\nRot=',Rot,'\n----------\nt=',C)
print('\nError = ',err)
#-------------------------------
with plt.style.context('Solarize_Light2'):
plt.figure(num='Error vs Noise')
plt.plot(errs,'r--',label='Error in pixels')
plt.legend()
# plt.xscale('symlog')
# plt.xlim(0., 10.0)
plt.xlabel(r'$\sigma$',fontsize=14)
plt.ylabel('Error in pixels',fontsize=14)
plt.show()
if __name__ == "__main__":
DLT()