Skip to content

Instantly share code, notes, and snippets.

@olivierdalang olivierdalang/gist:ba97fc986ade4545068d Secret
Last active Aug 29, 2015

Embed
What would you like to do?
Learn more about clone URLs
Download ZIP
worldfile transform
Raw
gistfile1.py
""" 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
You can’t perform that action at this time.