Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- """
- """
- import sys
- from constants import *
- from numpy import *
- from surveydata import *
- class SurveyMath():
- def __init__(self):
- self.initEOP()
- self.initImagePoints()
- self.EOP = array([], double)
- self.EOPStd = array([], double)
- self.sigma0 = array([], double)
- self.nPoints = len(self.pixs)
- self.nImages = len(self.pixs[0])
- #~ print self.pixs
- self.EOP.resize(self.nImages, 6)
- self.EOPStd.resize(self.nImages, 6)
- self.sigma0.resize(self.nImages, 1)
- for i in range(0, self.nImages):
- pass
- EOP, EOPStd, sigma0 = self.spaceResection(self.gps, self.pixs[:,i], affPara, f, ppxo, ppyo, self.EOPo_set, 0, 0, 0, 0, 0)
- self.EOP[i,:] = EOP
- self.EOPStd[i,:] = EOPStd
- self.sigma0[i,:] = sigma0
- #~ EOP, EOPStd, sigma0 = self.spaceResection(self.gps,\
- #~ self.pixs[:,1], array([-3.0656, 0.001681, 0.000000, 2.2979, 0.000000, -0.001680], double),\
- #~ f, ppxo, ppyo, self.EOPo_set, 0, 0, 0, 0, 0)
- #~ self.EOP[1,:] = EOP
- #~ self.EOPStd[1,:] = EOPStd
- #~ self.sigma0[1,:] = sigma0
- #~ self.
- #~ for i in range(0,l):
- #~ self.EOP[:,i], self.EOPStd[:,i], self.sigma0[:,i] = self.spaceResection(self.gps,\
- #~ self.pixs[:,i], affPara, f, ppxo, ppyo, self.EOPo_set, 0, 0, 0, 0, 0)
- #~ EOP, EOPStd, sigma0 = self.spaceResection(self.gps,\
- #~ self.pixs[:,i], affPara, f, ppxo, ppyo, self.EOPo_set, 0, 0, 0, 0, 0)
- #~ print EOP
- #~ self.EOP[i,:] = EOP
- #~ self.EOPStd.append(EOPStd)
- #~ self.sigma0.append(sigma0)
- #~ print self.EOP[:,1]
- #~ for i in range(0,l):
- #~
- for i in range(0, 1):
- self.spaceIntersection(self.pixs[i], affPara, self.EOP, IOP)
- #~ print self.pixs[i][0]
- #~ self.EOPRight, self.EOPStdRight, self.sigma0Right = self.spaceResection(self.gps, self.pixR, affParaRight, fRight, ppxoRight,
- #~ ppyoRight, self.EOPo_set_right, 0, 0, 0, 0, 0)
- #~ self.spaceIntersection(self.gps, self.pixL, self.pixR, self.EOPLeft, self.EOPRight, leftIOP, rightIOP, affParaLeft, affParaRight)
- #~ self.spaceIntersection(self.gps, self.pix
- #~ self.spaceIntersection(self.gps, self.pix, affPara, self.EOP, IOP)
- #~ print self.EOPo_set_left
- #~ print self.EOPo_set_right
- def initEOP(self):
- self.EOPo_set = array(p0_EOP, double)
- #~ self.EOPo_set_right = array(p0_right_EOP, double)
- #~ print self.EOPo_set_left[0:3]
- #~ print self.EOPo_set_right[0:3]
- self.EOPo_set[0:3] = self.EOPo_set[0:3] / 180 * pi
- #~ self.EOPo_set_right[0:3] = self.EOPo_set_right[0:3] / 180 * pi
- self.EOPo_set = self.EOPo_set[0:6]
- #~ self.EOPo_set_right = self.EOPo_set_right[0:6]
- def initImagePoints(self):
- fileImg1 = open('p0_leftCamA_imagePnts.txt', 'r')
- fileImg2 = open('p0_rightCam_imagePnts.txt', 'r')
- fileCtrl = open('petrieControl_NED.txt', 'r')
- pixL = []
- pixR = []
- gps = []
- i = 0
- self.pixs = zeros((31, 2, 2), double)
- while 1:
- line1 = fileImg1.readline()
- line2 = fileImg2.readline()
- #~ lineR = fileR.readline()
- lineCtrl = fileCtrl.readline()
- if not lineCtrl:
- break
- #~ pixL.append(array([[line1.split()[1], line1.split()[2]],[line2.split()[1], line2.split()[2]]], double))
- #~ pixR.append([lineR.split()[1], lineR.split()[2]])
- self.pixs[i] = array([[line1.split()[1], line1.split()[2]], [line2.split()[1], line2.split()[2]]], double)
- #~ pixL.append(array([line1.split()[1], line1.split()[2]], double))
- #~ pixR.append(array([line2.split()[1], line2.split()[2]], double))
- gps.append([lineCtrl.split()[1], lineCtrl.split()[2],
- lineCtrl.split()[3]])
- i += 1
- #~ self.pixs = array([pixL, pixR], double)
- #~ self.pixR = array(pixR, double)
- self.gps = array(gps, double)
- def spaceResection(self, gps, pix, affPara, f, ppx0, ppy0, EOPo, k1, k2, k3, p1, p2):
- photoGCP = array([], double)
- length = len(gps)
- photoGCP.resize(pix.shape)
- photoGCP[:,0:1] = (affPara[1]*pix[:,0:1])+(affPara[2]*pix[:,1:2]) + affPara[0]
- photoGCP[:,1:2] = (affPara[4]*pix[:,0:1])+(affPara[5]*pix[:,1:2]) + affPara[3]
- delta = EOPo.transpose()
- tolerance = ones(len(delta), double) * 0.001
- B2 = zeros((62,6), double)
- epsilon = zeros(2*length, double)
- Xnot = EOPo
- count = 0
- while count < 50:
- row = 1
- omega = Xnot[0]
- phi = Xnot[1]
- kappa = Xnot[2]
- X1 = Xnot[3]
- Y1 = Xnot[4]
- Z1 = Xnot[5]
- M = rotation3d(omega, phi, kappa)
- i = 0
- while i < length:
- r = (M[0,0]*(gps[i,0] - X1)) + (M[0,1]*(gps[i,1] - Y1)) + (M[0,2]*(gps[i,2] - Z1))
- s = (M[1,0]*(gps[i,0] - X1)) + (M[1,1]*(gps[i,1] - Y1)) + (M[1,2]*(gps[i,2] - Z1))
- q = (M[2,0]*(gps[i,0] - X1)) + (M[2,1]*(gps[i,1] - Y1)) + (M[2,2]*(gps[i,2] - Z1))
- B2[(2*(i+1)-1)-1,0] = (f/(q**2))*((r*((-M[2,2]*(gps[i,1] - Y1))+
- (M[2,1]*(gps[i,2] - Z1))))-(q*((-M[0,2]*(gps[i,1] - Y1))+
- (M[0,1]*(gps[i,2] - Z1)))))
- B2[(2*(i+1)-1)-1,1] = (f/(q**2))*((r*((cos(phi)*(gps[i,0] - X1))+
- (sin(omega)*sin(phi)*(gps[i,1] - Y1)) - (cos(omega)*sin(phi)*
- (gps[i,2] - Z1))))-(q*((-sin(phi)*cos(kappa)*(gps[i,0] - X1))+
- (sin(omega)*cos(phi)*cos(kappa)*(gps[i,1] - Y1)) - (cos(omega)*
- cos(phi)*cos(kappa)*(gps[i,2] - Z1)))))
- B2[(2*(i+1)-1)-1,2] = (-f/q)*((M[1,0]*(gps[i,0] - X1))+(M[1,1]*
- (gps[i,1] - Y1))+(M[1,2]*(gps[i,2] - Z1)))
- B2[(2*(i+1)-1)-1,3] = -(f/(q**2))*((r*M[2,0]) - (q*M[0,0]))
- B2[(2*(i+1)-1)-1,4] = -(f/(q**2))*((r*M[2,1]) - (q*M[0,1]))
- B2[(2*(i+1)-1)-1,5] = -(f/(q**2))*((r*M[2,2]) - (q*M[0,2]))
- B2[2*(i+1)-1,0] = (f/(q**2))* ((s*((-M[2,2]*(gps[i,1] - Y1))+(M[2,1]*
- (gps[i,2] - Z1))))-(q*((-M[1,2]*(gps[i,1] - Y1))+(M[1,1]*
- (gps[i,2] - Z1)))))
- B2[2*(i+1)-1,1] = (f/(q**2))*((s*((cos(phi)*(gps[i,0] - X1))+
- (sin(omega)*sin(phi)*(gps[i,1] - Y1)) - (cos(omega)*sin(phi)*
- (gps[i,2] - Z1))))-(q*((sin(phi)*sin(kappa)*(gps[i,0] - X1))-
- (sin(omega)*cos(phi)*sin(kappa)*(gps[i,1] - Y1)) + (cos(omega)*
- cos(phi)*sin(kappa)*(gps[i,2] - Z1)))))
- B2[2*(i+1)-1,2] = (f/q)*((M[0,0]*(gps[i,0] - X1))+(M[0,1]*(gps[i,1] - Y1))
- +(M[0,2]*(gps[i,2] - Z1)))
- B2[2*(i+1)-1,3] = -(f/(q**2))*((s*M[2,0]) - (q*M[1,0]))
- B2[2*(i+1)-1,4] = -(f/(q**2))*((s*M[2,1]) - (q*M[1,1]))
- B2[2*(i+1)-1,5] = -(f/(q**2))*((s*M[2,2]) - (q*M[1,2]))
- rDist2 = ((photoGCP[i,0]**2)+(photoGCP[i,1]**2));
- deltaX = photoGCP[i,0]*rDist2*(k1+rDist2*(k2+rDist2*k3)) + p1*(rDist2+(2*photoGCP[i,0]**2)) + 2*p2*photoGCP[i,0]*photoGCP[i,1]
- deltaY = photoGCP[i,1]*rDist2*(k1+rDist2*(k2+rDist2*k3)) + p1*(rDist2+ (2*photoGCP[i,1]**2)) + 2*p2*photoGCP[i,0]*photoGCP[i,1]
- epsilon[(2*(i+1)-1)-1] = photoGCP[i,0] - ppx0 + (f*(r/q)) + deltaX
- epsilon[2*(i+1)-1] = photoGCP[i,1] - ppy0 + (f*(s/q)) + deltaY
- i +=1
- delta = dot(linalg.inv(dot(B2.transpose(), B2)), dot(B2.transpose(), epsilon.transpose()))
- Xnot += delta
- count += 1
- EOP = Xnot
- sigma0 = dot(epsilon, epsilon.transpose())
- EOPstd = sqrt(sigma0*(diag(linalg.inv(dot(B2.transpose(), B2)))))
- return EOP, EOPstd, sigma0
- def spaceIntersection(self, pixs, affPara, EOP, IOP):
- X = array(EOP[:,3:4], double)
- Y = array(EOP[:,4:5], double)
- Z = array(EOP[:,5:6], double)
- M = array([], double)
- M.resize(self.nImages, 3, 3, 1)
- #~ for i in range(0,l):
- #~ M[i] = rotation3d(EOP[0],EOP[1],EOP[2])
- #~ M = rotation3d(EOP[0],EOP[1],EOP[2])
- f = IOP[0]
- ppxo = IOP[1]
- ppyo = IOP[2]
- k1 = IOP[3]
- k2 = IOP[4]
- k3 = IOP[5]
- p1 = IOP[6]
- p2 = IOP[7]
- photo = array([], double)
- photo.resize(self.nImages, self.nPoints, 2)
- #~ print photo.shape
- #~ for i in range(0,self.nImages):
- #~ print M[i,:]
- #~ print rotation3d(EOP[i,0:1], EOP[i,1:2], EOP[i,2:3])
- #~ print EOP[0,0:1]
- #~ print EOP[5,0:1]
- for i in range(0, self.nImages):
- M[i] = rotation3d(EOP[i,0:1], EOP[i,1:2], EOP[i,2:3])
- #~ print str(pixs[i,:,0:1]) + '\n'
- photo[i,:,0:1] = (affPara[1]*pixs[i,:,0:1])+(affPara[2]*pixs[i,:,1:2])+affPara[0]
- photo[i,:,1:2] = (affPara[4]*pixs[i,:,0:1])+(affPara[5]*pixs[i,:,1:2])+affPara[3]
- #~ print M[i]
- #~ M[i] =
- #~ photo[i,:,0:1] = (affPara[1]*pixs[i,:,0:1])+(affPara[2]*pixs[i,:,1:2]) + affPara[0]
- #~ print str(pixs[i,:][:]) + '\n'
- #~ print M
- X0 = array([32.5600, -12.7000, -1.2500], double)
- M.resize(self.nImages, 3, 3)
- B1 = zeros((2, 3), double)
- epsilon1 = zeros((self.nPoints, 2, 1))
- epsilon = zeros((self.nPoints, 4, 1))
- #~ for itera in range(0, 20):
- #~ for i in range(0, self.nImages):
- for itera in range(0, 20):
- for i in range(0, 1):
- Brow = 0
- r1 = (M[i,0,0]*(X0[0] - X[i])) + (M[i,0,1]*(X0[1] - Y[i])) + (M[i,0,2]*(X0[2] - Z[i]))
- s1 = (M[i,1,0]*(X0[0] - X[i])) + (M[i,1,1]*(X0[1] - Y[i])) + (M[i,1,2]*(X0[2] - Z[i]))
- q1 = (M[i,2,0]*(X0[0] - X[i])) + (M[i,2,1]*(X0[1] - Y[i])) + (M[i,2,2]*(X0[2] - Z[i]))
- #~ B1[2*(i+1)-1-1][0][0]
- B1[0][0] = (f/(q1**2))*((r1*M[i,2,0])-(q1*M[i,0,0]))
- B1[0][1] = (f/(q1**2))*((r1*M[i,2,1])-(q1*M[i,0,1]))
- B1[0][2] = (f/(q1**2))*((r1*M[i,2,2])-(q1*M[i,0,2]))
- B1[1][0] = (f/(q1**2))*((s1*M[i,2,0])-(q1*M[i,1,0]))
- B1[1][1] = (f/(q1**2))*((s1*M[i,2,1])-(q1*M[i,1,1]))
- B1[1][2] = (f/(q1**2))*((s1*M[i,2,2])-(q1*M[i,1,2]))
- #~ print B1
- if i == 0:
- B = B1
- rDist2 = ((photo[i,0:1,0]**2)+(photo[i,0:1,1]**2))
- deltaX = photo[i,0:1,0]*rDist2*(k1+rDist2*(k2+rDist2*k3)) + p1*(rDist2+ (2*photo[i,0:1,0]**2)) + 2*p2*photo[i,0:1,0]*photo[i,0:1,1]
- deltaY = photo[i,0:1,1]*rDist2*(k1+rDist2*(k2+rDist2*k3)) + p1*(rDist2+ (2*photo[i,0:1,1]**2)) + 2*p2*photo[i,0:1,0]*photo[i,0:1,1]
- epsilon1[i,2*Brow-1-1,0:1] = photo[i,0:1,0] - ppxo + (f*(r1/q1)) + deltaX #pnt(Brow,1) - ppx01 + (fLeft*(r1/q1));
- epsilon1[i,2*Brow-1,0:1] = photo[i,0:1,1] - ppyo + (f*(s1/q1)) + deltaY #pnt(Brow,2) - ppy01 + (fLeft*(s1/q1));
- Brow = Brow + 1
- #~ photo[:,0:1] = (affPara[1]*pix[:,0:1])+(affPara[2]*pix[:,1:2]) + affPara[0]
- #~ photo[:,1:2] = (affPara[4]*pix[:,0:1])+(affPara[5]*pix[:,1:2]) + affPara[3]
- #~
- #~ for i in range(0,20):
- #~ print i
- #~ def spaceIntersection(self, gps, pixL, pixR, EOPleft, EOPright, IOPleft, IOPright, affParaLeft, affParaRight):
- #~ X = EOPleft[3]
- #~ Y = EOPleft[4]
- #~ Z = EOPleft[5]
- #~
- #~ X_right = EOPright[3]
- #~ Y_right = EOPright[4]
- #~ Z_right = EOPright[5]
- #~
- #~ M = rotation3d(EOPleft[0],EOPleft[1],EOPleft[2])
- #~ M_right = rotation3d(EOPright[0],EOPright[1],EOPright[2])
- #~
- #~ fLeft = IOPleft[0]
- #~ ppxoLeft = IOPleft[1]
- #~ ppyoLeft = IOPleft[2]
- #~ k1Left = IOPleft[3]
- #~ k2Left = IOPleft[4]
- #~ k3Left = IOPleft[5]
- #~ p1Left = IOPleft[6]
- #~ p2Left = IOPleft[7]
- #~
- #~ fRight = IOPright[0]
- #~ ppxoRight = IOPright[1]
- #~ ppyoRight = IOPright[2]
- #~ k1Right = IOPright[3]
- #~ k2Right = IOPright[4]
- #~ k3Right = IOPright[5]
- #~ p1Right = IOPright[6]
- #~ p2Right = IOPright[7]
- #~
- #~ photoLeft = array([], double)
- #~ photoRight = array([], double)
- #~
- #~ length = len(gps)
- #~ photoLeft.resize(pixL.shape)
- #~ photoRight.resize(pixR.shape)
- #~
- #~ photoLeft[:,0:1] = (affParaLeft[1]*pixL[:,0:1])+(affParaLeft[2]*pixL[:,1:2]) + affParaLeft[0]
- #~ photoLeft[:,1:2] = (affParaLeft[4]*pixL[:,0:1])+(affParaLeft[5]*pixL[:,1:2]) + affParaLeft[3]
- #~
- #~ photoRight[:,0:1] = (affParaRight[1]*pixR[:,0:1])+(affParaRight[2]*pixR[:,1:2]) + affParaRight[0]
- #~ photoRight[:,1:2] = (affParaRight[4]*pixR[:,0:1])+(affParaRight[5]*pixR[:,1:2]) + affParaRight[3]
- #~
- #~ B1 = zeros((length,2, 3), double)
- #~ B2 = zeros((length,2, 3), double)
- #~ B = zeros((length, 4, 3), double)
- #~ epsilon1 = zeros((length, 2, 1))
- #~ epsilon2 = zeros((length, 2, 1))
- #~ epsilon = zeros((length, 4, 1))
- #~
- #~
- #~ i = 0
- #~ while i < 50:
- #~ Brow = 1
- #~ r1 = (M[0,0]*(gps[:,0:1] - X)) + (M[0,1]*(gps[:,1:2] - Y)) + (M[0,2]*(gps[:,2:3] - Z))
- #~ s1 = (M[1,0]*(gps[:,0:1] - X)) + (M[1,1]*(gps[:,1:2] - Y)) + (M[1,2]*(gps[:,2:3] - Z))
- #~ q1 = (M[2,0]*(gps[:,0:1] - X)) + (M[2,1]*(gps[:,1:2] - Y)) + (M[2,2]*(gps[:,2:3] - Z))
- #~
- #~ B1[:,2*Brow-1-1,0:1] = (fLeft/(q1**2))*((r1*M[2,0])-(q1*M[0,0]))
- #~ B1[:,2*Brow-1-1,1:2]= (fLeft/(q1**2))*((r1*M[2,1])-(q1*M[0,1]))
- #~ B1[:,2*Brow-1-1,2:3] = (fLeft/(q1**2))*((r1*M[2,2])-(q1*M[0,2]))
- #~
- #~ B1[:,2*Brow-1,0:1] = (fLeft/(q1**2))*((s1*M[2,0])-(q1*M[1,0]))
- #~ B1[:,2*Brow-1,1:2] = (fLeft/(q1**2))*((s1*M[2,1])-(q1*M[1,1]))
- #~ B1[:,2*Brow-1,2:3] = (fLeft/(q1**2))*((s1*M[2,2])-(q1*M[1,2]))
- #~
- #~
- #~ rDist2 = ((photoLeft[:,0:1]**2)+(photoLeft[:,1:2]**2))
- #~ deltaX = photoLeft[:,0:1]*rDist2*(k1Left+rDist2*(k2Left+rDist2*k3Left)) + p1Left*(rDist2+ (2*photoLeft[:,0:1]**2)) + 2*p2Left*photoLeft[:,0:1]*photoLeft[:,1:2]
- #~ deltaY = photoLeft[:,1:2]*rDist2*(k1Left+rDist2*(k2Left+rDist2*k3Left)) + p1Left*(rDist2+ (2*photoLeft[:,1:2]**2)) + 2*p2Left*photoLeft[:,0:1]*photoLeft[:,1:2]
- #~
- #~
- #~ epsilon1[:,2*Brow-1-1,0:1] = photoLeft[:,0:1] - ppxoLeft + (fLeft*(r1/q1)) + deltaX #pnt(Brow,1) - ppx01 + (fLeft*(r1/q1));
- #~ epsilon1[:,2*Brow-1,0:1] = photoLeft[:,1:2] - ppyoLeft + (fLeft*(s1/q1)) + deltaY #pnt(Brow,2) - ppy01 + (fLeft*(s1/q1));
- #~
- #~ r2 = (M_right[0,0]*(gps[:,0:1] - X_right)) + (M_right[0,1]*(gps[:,1:2] - Y_right)) + (M_right[0,2]*(gps[:,2:3] - Z_right))
- #~ s2 = (M_right[1,0]*(gps[:,0:1] - X_right)) + (M_right[1,1]*(gps[:,1:2] - Y_right)) + (M_right[1,2]*(gps[:,2:3] - Z_right))
- #~ q2 = (M_right[2,0]*(gps[:,0:1] - X_right)) + (M_right[2,1]*(gps[:,1:2] - Y_right)) + (M_right[2,2]*(gps[:,2:3] - Z_right))
- #~
- #~ B2[:,2*Brow-1-1,0:1] = (fRight/(q1**2))*((r2*M_right[2,0])-(q2*M_right[0,0]))
- #~ B2[:,2*Brow-1-1,1:2]= (fRight/(q1**2))*((r2*M_right[2,1])-(q2*M_right[0,1]))
- #~ B2[:,2*Brow-1-1,2:3] = (fRight/(q1**2))*((r2*M_right[2,2])-(q2*M_right[0,2]))
- #~
- #~ B2[:,2*Brow-1,0:1] = (fRight/(q2**2))*((s2*M_right[2,0])-(q2*M_right[1,0]))
- #~ B2[:,2*Brow-1,1:2] = (fRight/(q2**2))*((s2*M_right[2,1])-(q2*M_right[1,1]))
- #~ B2[:,2*Brow-1,2:3] = (fRight/(q2**2))*((s2*M_right[2,2])-(q2*M_right[1,2]))
- #~
- #~ rDist2 = ((photoRight[:,0:1]**2)+(photoRight[:,1:2]**2))
- #~ deltaX = photoRight[:,0:1]*rDist2*(k1Right+rDist2*(k2Right+rDist2*k3Right)) + p1Right*(rDist2+ (2*photoRight[:,0:1]**2)) + 2*p2Right*photoRight[:,0:1]*photoRight[:,1:2]
- #~ deltaY = photoRight[:,1:2]*rDist2*(k1Right+rDist2*(k2Right+rDist2*k3Right)) + p1Right*(rDist2+ (2*photoRight[:,1:2]**2)) + 2*p2Right*photoRight[:,0:1]*photoRight[:,1:2]
- #~
- #~ epsilon2[:,2*Brow-1-1,0:1] = photoRight[:,0:1] - ppyoRight + (fRight*(r2/q2)) + deltaX #pnt(Brow,1) - ppx01 + (fLeft*(r1/q1));
- #~ epsilon2[:,2*Brow-1,0:1] = photoRight[:,1:2] - ppyoRight + (fRight*(s2/q2)) + deltaY #pnt(Brow,2) - ppy01 + (fLeft*(s1/q1));
- #~
- #~
- #~ B[:,2:4,:] = B2
- #~
- #~ epsilon[:,0:2,:] = epsilon1
- #~ epsilon[:,2:4,:] = epsilon2
- #~ x = 0
- #~ while x < length:
- #~ b = B[x]
- #~ e = epsilon[x]
- #~ deltaX = dot(dot(linalg.inv(dot(b.T, b)), b.T), e)
- #~ gps[x] += deltaX
- #~
- #~ x +=1
- #~
- #~
- #~ print B.transpose()
- #~ deltaX = dot(dot(linalg.inv(dot(B.transpose(), B)), B.transpose()), epsilon)
- #~
- #~ i += 1
- surveyMath = SurveyMath()
Advertisement
Add Comment
Please, Sign In to add comment