worldfile transform
""" Notes : | |
pointsPairs: array of points pairs, where the first point is in pixel coordinates, the second in world coordinates | |
lockRot: whether rotation should be locked | |
lockRatio: whether pixel ratio should be locked to 1/1 | |
Variables A,B,C,D,E,F described here : http://en.wikipedia.org/wiki/World_file | |
""" | |
def _leastSquareTransform(self, pointsPairs , lockRot, lockRatio): | |
"""Computes least squares transform for overdetermined system""" | |
if not lockRot and not lockRatio: | |
# We do a least square fit | |
# We need more than 3 points to be in a overdetermined situation | |
assert(len(pointsPairs)>3) | |
""" | |
# Definition of one point | |
pWldXn = A*pPixXn + B*pPixYn + C | |
pWldYn = D*pPixXn + E*pPixYn + F | |
# Definition of one point with errors | |
pWldXn = A*pPixXn + B*pPixYn + C + eXn | |
pWldYn = D*pPixXn + E*pPixYn + F + eYn | |
# Definition of one point's error | |
eXn = (pWldXn - A*pPixXn - B*pPixYn - C) | |
eYn = (pWldYn - D*pPixXn - E*pPixYn - F) | |
# We want the sum of errors squared to converge towards 0 | |
1) DER_A( SUM( eXn**2 )) = 0 | |
2) DER_B( SUM( eXn**2 )) = 0 | |
3) DER_C( SUM( eXn**2 )) = 0 | |
4) DER_D( SUM( eYn**2 )) = 0 | |
5) DER_E( SUM( eYn**2 )) = 0 | |
6) DER_F( SUM( eYn**2 )) = 0 | |
# Terms : | |
o = pWldXn | |
p = pWldYn | |
x = pPixXn | |
y = pPixYn | |
eXn = (o - A*x - B*y - C) | |
eYn = (p - D*x - E*y - F) | |
eXn^2 = C^2 + 2*y*B*C + 2*x*A*C - 2*o*C + y^2*B^2 + 2*x*y*A*B - 2*o*y*B + x^2*A^2 - 2*o*x*A + o^2 | |
eYn^2 = F^2 + 2*y*E*F + 2*x*D*F - 2*p*F + y^2*E^2 + 2*x*y*D*E - 2*p*y*E + x^2*D^2 - 2*p*x*D + p^2 | |
# Sums : | |
o = SUM(pWldXn) | |
p = SUM(pWldYn) | |
x = SUM(pPixXn) | |
y = SUM(pPixYn) | |
u = SUM(pPixXn^2) # x^2 | |
v = SUM(pPixYn^2) # y^2 | |
w = SUM(pPixXn*pPixYn) # x*y | |
j = SUM(pWldXn*pPixXn) # o*x | |
k = SUM(pWldXn*pPixYn) # o*y | |
l = SUM(pWldYn*pPixXn) # p*x | |
m = SUM(pWldyn*pPixYn) # p*y | |
n = SUM(1) | |
SUM(eXn^2) = n*C^2 + 2*y*B*C + 2*x*A*C - 2*o*C + v*B^2 + 2*w*A*B - 2*k*B + u*A^2 - 2*j*A + n*o^2 | |
SUM(eYn^2) = n*F^2 + 2*y*E*F + 2*x*D*F - 2*p*F + v*E^2 + 2*w*D*E - 2*m*E + u*D^2 - 2*l*D + n*p^2 | |
1) 2*x*C+2*w*B+2*u*A-2*j = 0 | |
2) 2*y*C+2*v*B+2*w*A-2*k = 0 | |
3) 2*n*C+2*y*B+2*x*A-2*o = 0 | |
4) 2*x*F+2*w*E+2*u*D-2*l = 0 | |
5) 2*y*F+2*v*E+2*w*D-2*m = 0 | |
6) 2*n*F+2*y*E+2*x*D-2*p = 0 | |
System1 (A;B;C) : 2*x*C+2*w*B+2*u*A-2*j = 0, 2*y*C+2*v*B+2*w*A-2*k = 0, 2*n*C+2*y*B+2*x*A-2*o = 0 | |
System2 (D;E;F) : 2*x*F+2*w*E+2*u*D-2*l = 0, 2*y*F+2*v*E+2*w*D-2*m = 0, 2*n*F+2*y*E+2*x*D-2*p = 0 | |
# Results | |
A = -(j*(n*v-y^2)+w*(o*y-k*n)+x*(k*y-o*v))/(u*(y^2-n*v)-2*w*x*y+v*x^2+n*w^2) | |
B = (u*(o*y-k*n)+x*(-j*y-o*w)+k*x^2+j*n*w)/(u*(y^2-n*v)-2*w*x*y+v*x^2+n*w^2) | |
C = -(u*(o*v-k*y)+j*w*y+(k*w-j*v)*x-o*w^2)/(u*(y^2-n*v)-2*w*x*y+v*x^2+n*w^2) | |
D = -(l*(n*v-y^2)+w*(p*y-m*n)+x*(m*y-p*v))/(u*(y^2-n*v)-2*w*x*y+v*x^2+n*w^2) | |
E = (u*(p*y-m*n)+x*(-l*y-p*w)+m*x^2+l*n*w)/(u*(y^2-n*v)-2*w*x*y+v*x^2+n*w^2) | |
F = -(u*(p*v-m*y)+l*w*y+(m*w-l*v)*x-p*w^2)/(u*(y^2-n*v)-2*w*x*y+v*x^2+n*w^2) | |
""" | |
o,p,x,y,u,v,w,j,k,l,m,n = 0,0,0,0,0,0,0,0,0,0,0,0 | |
for pair in pointsPairs: | |
o += pair[1].x() | |
p += pair[1].y() | |
x += pair[0].x() | |
y += pair[0].y() | |
u += pair[0].x()**2 | |
v += pair[0].y()**2 | |
w += pair[0].x()*pair[0].y() | |
j += pair[1].x()*pair[0].x() # o*x | |
k += pair[1].x()*pair[0].y() # o*y | |
l += pair[1].y()*pair[0].x() # p*x | |
m += pair[1].y()*pair[0].y() # p*y | |
n += 1 | |
self.a = -(j*(n*v-y**2)+w*(o*y-k*n)+x*(k*y-o*v))/(u*(y**2-n*v)-2*w*x*y+v*x**2+n*w**2) | |
self.b = (u*(o*y-k*n)+x*(-j*y-o*w)+k*x**2+j*n*w)/(u*(y**2-n*v)-2*w*x*y+v*x**2+n*w**2) | |
self.c = -(u*(o*v-k*y)+j*w*y+(k*w-j*v)*x-o*w**2)/(u*(y**2-n*v)-2*w*x*y+v*x**2+n*w**2) | |
self.d = -(l*(n*v-y**2)+w*(p*y-m*n)+x*(m*y-p*v))/(u*(y**2-n*v)-2*w*x*y+v*x**2+n*w**2) | |
self.e = (u*(p*y-m*n)+x*(-l*y-p*w)+m*x**2+l*n*w)/(u*(y**2-n*v)-2*w*x*y+v*x**2+n*w**2) | |
self.f = -(u*(p*v-m*y)+l*w*y+(m*w-l*v)*x-p*w**2)/(u*(y**2-n*v)-2*w*x*y+v*x**2+n*w**2) | |
elif lockRot and not lockRatio: | |
# same developement, but lockRot means we have no rotation, aka B=0 and D=0, so we remove all terms containing those | |
# We do a least square fit | |
# We need more than 2 points to be in a overdetermined situation | |
assert(len(pointsPairs)>2) | |
""" | |
System1 (A;C) : 2*x*C+2*u*A-2*j = 0,2*n*C+2*x*A-2*o = 0 | |
System2 (E;F) : 2*y*F+2*v*E-2*m = 0,2*n*F+2*y*E-2*p = 0 | |
# Results | |
A = (o*x-j*n)/(x^2-n*u) | |
B = 0 | |
C = (j*x-o*u)/(x^2-n*u) | |
D = 0 | |
E = (p*y-m*n)/(y^2-n*v) | |
F = (m*y-p*v)/(y^2-n*v) | |
""" | |
o,p,x,y,u,v,w,j,k,l,m,n = 0,0,0,0,0,0,0,0,0,0,0,0 | |
for pair in pointsPairs: | |
o += pair[1].x() | |
p += pair[1].y() | |
x += pair[0].x() | |
y += pair[0].y() | |
u += pair[0].x()**2 | |
v += pair[0].y()**2 | |
j += pair[1].x()*pair[0].x() # o*x | |
m += pair[1].y()*pair[0].y() # p*y | |
n += 1 | |
self.a = (o*x-j*n)/(x**2-n*u) | |
self.b = 0 | |
self.c = (j*x-o*u)/(x**2-n*u) | |
self.d = 0 | |
self.e = (p*y-m*n)/(y**2-n*v) | |
self.f = (m*y-p*v)/(y**2-n*v) | |
elif not lockRot and lockRatio: | |
# same developement, but lockRatio means we have square pixels, aka A=-E and B=D, so we replace those terms | |
# We do a least square fit | |
# We need more than 2 points to be in a overdetermined situation | |
assert(len(pointsPairs)>2) | |
""" | |
System1 (A;C;D;F) : 2*x*C+2*w*D+2*u*A-2*j = 0, 2*n*C+2*y*D+2*x*A-2*o = 0, 2*x*F+2*w*(-A)+2*u*D-2*l = 0, 2*n*F+2*y*(-A)+2*x*D-2*p = 0 | |
# Results | |
A = (x^2*(-p*y-j*n)+x*(l*n*y+n*p*w-n*o*u)+o*x^3-l*n^2*w+j*n^2*u)/(x^2*(y^2-2*n*u)-2*n*w*x*y+x^4+n^2*w^2+n^2*u^2) | |
C = (x*(j*y^2+u*(p*y-j*n)+w*(l*n-o*y))-j*n*w*y-l*n*u*y+j*x^3+(-p*w-o*u)*x^2+n*o*w^2+n*o*u^2)/(x^2*(y^2-2*n*u)-2*n*w*x*y+x^4+n^2*w^2+n^2*u^2) | |
D = (x^2*(o*y-l*n)+x*(-j*n*y-n*o*w-n*p*u)+p*x^3+j*n^2*w+l*n^2*u)/(x^2*(y^2-2*n*u)-2*n*w*x*y+x^4+n^2*w^2+n^2*u^2) | |
F = (x*(l*y^2+w*(-p*y-j*n)+u*(-o*y-l*n))-l*n*w*y+j*n*u*y+l*x^3+(o*w-p*u)*x^2+n*p*w^2+n*p*u^2)/(x^2*(y^2-2*n*u)-2*n*w*x*y+x^4+n^2*w^2+n^2*u^2) | |
""" | |
o,p,x,y,u,v,w,j,k,l,m,n = 0,0,0,0,0,0,0,0,0,0,0,0 | |
for pair in pointsPairs: | |
o += pair[1].x() | |
p += pair[1].y() | |
x += pair[0].x() | |
y += pair[0].y() | |
u += pair[0].x()**2 | |
w += pair[0].x()*pair[0].y() | |
j += pair[1].x()*pair[0].x() # o*x | |
l += pair[1].y()*pair[0].x() # p*x | |
n += 1 | |
self.a = (x**2*(-p*y-j*n)+x*(l*n*y+n*p*w-n*o*u)+o*x**3-l*n**2*w+j*n**2*u)/(x**2*(y**2-2*n*u)-2*n*w*x*y+x**4+n**2*w**2+n**2*u**2) | |
self.c = (x*(j*y**2+u*(p*y-j*n)+w*(l*n-o*y))-j*n*w*y-l*n*u*y+j*x**3+(-p*w-o*u)*x**2+n*o*w**2+n*o*u**2)/(x**2*(y**2-2*n*u)-2*n*w*x*y+x**4+n**2*w**2+n**2*u**2) | |
self.d = (x**2*(o*y-l*n)+x*(-j*n*y-n*o*w-n*p*u)+p*x**3+j*n**2*w+l*n**2*u)/(x**2*(y**2-2*n*u)-2*n*w*x*y+x**4+n**2*w**2+n**2*u**2) | |
self.f = (x*(l*y**2+w*(-p*y-j*n)+u*(-o*y-l*n))-l*n*w*y+j*n*u*y+l*x**3+(o*w-p*u)*x**2+n*p*w**2+n*p*u**2)/(x**2*(y**2-2*n*u)-2*n*w*x*y+x**4+n**2*w**2+n**2*u**2) | |
self.b = self.d | |
self.e = -self.a | |
else:#elif lockRot and lockRatio | |
# same developement, but lockRot and lockRatio means A=-E and B=D = 0, so we replace those terms | |
# We need more than 2 points to be in a overdetermined situation | |
assert(len(pointsPairs)>2) | |
""" | |
System1 (A;C;F) : 2*x*C+2*u*A-2*j = 0, 2*n*C+2*x*A-2*o = 0, 2*n*F+2*y*(-A)-2*p = 0 | |
# Results | |
A = (o*x-j*n)/(x^2-n*u) | |
C = (j*x-o*u)/(x^2-n*u) | |
F = (o*x*y-j*n*y+p*x^2-n*p*u)/(n*x^2-n^2*u) | |
""" | |
o,p,x,y,u,v,w,j,k,l,m,n = 0,0,0,0,0,0,0,0,0,0,0,0 | |
for pair in pointsPairs: | |
o += pair[1].x() | |
p += pair[1].y() | |
x += pair[0].x() | |
y += pair[0].y() | |
u += pair[0].x()**2 | |
j += pair[1].x()*pair[0].x() # o*x | |
n += 1 | |
self.a = (o*x-j*n)/(x**2-n*u) | |
self.b = 0 | |
self.c = (j*x-o*u)/(x**2-n*u) | |
self.d = 0 | |
self.e = -self.a | |
self.f = (o*x*y-j*n*y+p*x**2-n*p*u)/(n*x**2-n**2*u) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment