/gist:ba97fc986ade4545068d Secret
Last active
August 29, 2015 14:20
Star
You must be signed in to star a gist
worldfile transform
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" 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