bend raster map

fitbspline.py

@@ -1,6 +1,7 @@
 #!/usr/bin/env python
 
 import sys
+import math
 from math import sqrt
 import numpy as np
 import matplotlib.pyplot as plt
@@ -9,6 +10,13 @@ from itertools import tee, izip
 from matplotlib.patches import Polygon
 from matplotlib.collections import PatchCollection
 import matplotlib
+from PIL import Image
+
+def y2lat(a):
+    return 180.0/math.pi*(2.0*math.atan(math.exp(a*math.pi/180.0))-math.pi/2.0)
+
+def lat2y(a):
+    return 180.0/math.pi*math.log(math.tan(math.pi/4.0+a*(math.pi/180.0)/2.0))
 
 def pairwise(iterable):
     "s -> (s0,s1), (s1,s2), (s2,s3), ..."
@@ -43,7 +51,7 @@ def getxing(p0, p1, p2, p3):
         # lines are parallel and never meet
         return None
     s = (vy*(p3[0]-p0[0])+vx*(p0[1]-p3[1]))/a
-    if s < 1.0:
+    if 0.0 < s < 1.0:
         return (p0[0]+s*ux, p0[1]+s*uy)
     else:
         return None
@@ -61,11 +69,10 @@ def ptInQuadrilateral(p, p0, p1, p2, p3):
     # it might be that the two normals cross at some point
     # in that case the two triangles are created differently
     cross = getxing(p0, p1, p2, p3)
-    #if cross:
+    if cross:
-    #    return ptInTriangle(p, p0, cross, p3) or ptInTriangle(p, p2, cross, p1)
+        return ptInTriangle(p, p0, cross, p3) or ptInTriangle(p, p2, cross, p1)
-    #else:
+    else:
-    #    return ptInTriangle(p, p0, p1, p2) or ptInTriangle(p, p2, p3, p0)
+        return ptInTriangle(p, p0, p1, p2) or ptInTriangle(p, p2, p3, p0)
-    return ptInTriangle(p, p0, p1, p2) or ptInTriangle(p, p2, p3, p0)
 
 def get_st(Ax,Ay,Bx,By,Cx,Cy,Dx,Dy,Xx,Xy):
     d = Bx-Ax-Cx+Dx
@@ -110,38 +117,17 @@ def get_st(Ax,Ay,Bx,By,Cx,Cy,Dx,Dy,Xx,Xy):
             if -0.0000000001 <= r22 <= 1.0000000001:
                 t.append(r22)
     if not s or not t:
-        return None
+        return [],[]
     if len(s) == 1 and len(t) == 2:
         s = [s[0],s[0]]
     if len(s) == 2 and len(t) == 1:
         t = [t[0],t[0]]
     return s, t
 
-def find_coeffs(pa, pb):
+def main(x,y,width,smoothing,subdiv):
-    matrix = []
-    for p1, p2 in zip(pa, pb):
-        matrix.append([p1[0], p1[1], 1, 0, 0, 0, -p2[0]*p1[0], -p2[0]*p1[1]])
-        matrix.append([0, 0, 0, p1[0], p1[1], 1, -p2[1]*p1[0], -p2[1]*p1[1]])
-
-    A = np.matrix(matrix, dtype=np.float)
-    B = np.array(pb).reshape(8)
-
-    #res = np.dot(np.linalg.inv(A), B)
-    res = np.dot(np.linalg.inv(A.T * A) * A.T, B)
-    return np.array(res).reshape(8)
-
-def main():
-    width = 2/5.0
     halfwidth = width/2.0
-    x = []
+    tck,u = interpolate.splprep([x,y],s=smoothing)
-    y = []
+    unew = np.linspace(0,1.0,subdiv+1)
-    with open(sys.argv[1]) as f:
-        for l in f:
-            a,b = l.split()
-            x.append(float(a))
-            y.append(float(b))
-    tck,u = interpolate.splprep([x,y],s=10)
-    unew = np.arange(0,1.1,0.1)
     out = interpolate.splev(unew,tck)
     heights = []
     offs = []
@@ -202,7 +188,6 @@ def main():
     qy.append(out[1][-1]+dy/dl)
     quads = []
     patches = []
-    #for (p0x,p0y,p1x,p1y),(p3x,p3y,p2x,p2y) in pairwise(zip(px,py,qx,qy)):
     for (p3x,p3y,p2x,p2y),(p0x,p0y,p1x,p1y) in pairwise(zip(px,py,qx,qy)):
         quads.append(((p0x,p0y),(p1x,p1y),(p2x,p2y),(p3x,p3y)))
         polygon = Polygon(((p0x,p0y),(p1x,p1y),(p2x,p2y),(p3x,p3y)), True)
@@ -215,77 +200,77 @@ def main():
             if ptInQuadrilateral(pt,p0,p1,p2,p3):
                 found.append(i)
         if found:
-            if len(found) > 2:
+            if len(found) > 1:
-                print found
+                print "point found in two quads"
-            containingquad.append(found)
+                return None
+            containingquad.append(found[0])
         else:
-            print "can't find quad for point"
             containingquad.append(None)
-            #exit(1)
+    # check if the only points for which no quad could be found are in the
-    print containingquad
+    # beginning or in the end
-    trans = []
+    # find the first missing ones:
-    print width, height
+    for i,q in enumerate(containingquad):
-    srcquads = []
+        if q != None:
-    for off,h,srcquad in zip(offs,heights,quads):
+            break
-        #targetquad = ((0,height-off),(width,height-off),(width,height-off-h),(0,height-off-h))
+    # find the last missing ones
+    for j,q in izip(xrange(len(containingquad)-1, -1, -1), reversed(containingquad)):
+        if q != None:
+            break
+    # remove the first and last missing ones
+    if i != 0 or j != len(containingquad)-1:
+        containingquad = containingquad[i:j+1]
+        x = x[i:j+1]
+        y = y[i:j+1]
+    # check if there are any remaining missing ones:
+    if None in containingquad:
+        print "cannot find quad for point"
+        return None
+    for off,h in zip(offs,heights):
         targetquad = ((0,off+h),(width,off+h),(width,off),(0,off))
-        trans.append(find_coeffs(srcquad,targetquad))
         patches.append(Polygon(targetquad,True))
     tx = []
     ty = []
-    #targetquad = (0,height),(width,height),(width,0),(0,0)
-    #srcquad = (min(x),max(y)),(max(x),max(y)),(max(x),min(y)),(min(x),min(y)) 
-    #trans = find_coeffs(srcquad,targetquad)
-    #for (rx,ry) in zip(x,y):
-    #    a,b,c,d,e,f,g,h = trans
-    #    u = (a*rx + b*ry + c)/(g*rx + h*ry + 1)
-    #    v = (d*rx + e*ry + f)/(g*rx + h*ry + 1)
-    #    tx.append(u)
-    #    ty.append(v)
     assert len(containingquad) == len(x) == len(y)
-    assert len(out[0]) == len(out[1]) == len(px) == len(py) == len(qx) == len(qy) == len(quads)+1 == len(heights)+1 == len(offs)+1 == len(trans)+1
+    assert len(out[0]) == len(out[1]) == len(px) == len(py) == len(qx) == len(qy) == len(quads)+1 == len(heights)+1 == len(offs)+1
-    for (rx,ry),l in zip(zip(x,y),containingquad):
+    for (rx,ry),i in zip(zip(x,y),containingquad):
-        if not l:
+        if i == None:
             continue
-        for i in l[:1]:
+        (ax,ay),(bx,by),(cx,cy),(dx,dy) = quads[i]
-            (ax,ay),(bx,by),(cx,cy),(dx,dy) = quads[i]
+        s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,rx,ry)
-            s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,rx,ry)
+        # if more than one solution, take second 
-            if len(s) != 1 or len(t) != 1:
+        # TODO: investigate if this is always the right solution
-                print "fail"
-                exit(1)
-            #a,b,c,d,e,f,g,h = trans[i]
-            ##den = -a*e+a*h*ry+b*d-b*g*ry-d*h*rx+e*g*rx
-            ##tx.append((-b*f+b*ry+c*e-c*h*ry-e*rx+f*h*rx)/den)
-            ##ty.append((a*f-a*ry-c*d+c*g*ry+d*rx-f*g*rx)/den)
-            #u = (a*rx + b*ry + c)/(g*rx + h*ry + 1)
-            #v = (d*rx + e*ry + f)/(g*rx + h*ry + 1)
-            u = s[0]*width
-            v = offs[i]+t[0]*heights[i]
-            tx.append(u)
-            ty.append(v)
-    sx = []
-    sy = []
-    for ((x1,y1),(x2,y2)),((ax,ay),(bx,by),(cx,cy),(dx,dy)),off,h in zip(pairwise(zip(*out)),quads,offs,heights):
-        s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,x1,y1)
         if len(s) != 1 or len(t) != 1:
-            print "fail"
+            s = s[1]
-            exit(1)
+            t = t[1]
-        #u = (a*ax + b*ay + c)/(g*ax + h*ay + 1)
+        else:
-        #v = (d*ax + e*ay + f)/(g*ax + h*ay + 1)
+            s = s[0]
-        u = s[0]*width
+            t = t[0]
-        v = off+t[0]*h
+        u = s*width
-        sx.append(u)
+        v = offs[i]+t*heights[i]
-        sy.append(v)
+        tx.append(u)
-        #u = (a*bx + b*by + c)/(g*bx + h*by + 1)
+        ty.append(v)
-        #v = (d*bx + e*by + f)/(g*bx + h*by + 1)
+    #sx = []
-        s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,x2,y2)
+    #sy = []
-        if len(s) != 1 or len(t) != 1:
+    #for ((x1,y1),(x2,y2)),((ax,ay),(bx,by),(cx,cy),(dx,dy)),off,h in zip(pairwise(zip(*out)),quads,offs,heights):
-            print "fail"
+    #    s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,x1,y1)
-            exit(1)
+    #    if len(s) != 1 or len(t) != 1:
-        u = s[0]*width
+    #        return None
-        v = off+t[0]*h
+    #    u = s[0]*width
-        sx.append(u)
+    #    v = off+t[0]*h
-        sy.append(v)
+    #    sx.append(u)
+    #    sy.append(v)
+    #    s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,x2,y2)
+    #    if len(s) != 1 or len(t) != 1:
+    #        return None
+    #    u = s[0]*width
+    #    v = off+t[0]*h
+    #    sx.append(u)
+    #    sy.append(v)
+    im1 = Image.open("out.png")
+    im2 = Image.open("map.png")
+    bbox = [8.4320068359375,51.34090729023285,9.7119140625,53.556626004824615]
+    bbox[1] = lat2y(bbox[1])
+    bbox[3] = lat2y(bbox[3])
     colors = 100*np.random.rand(len(patches)/2)+100*np.random.rand(len(patches)/2)
     p = PatchCollection(patches, cmap=matplotlib.cm.jet, alpha=0.4)
     p.set_array(np.array(colors))
@@ -295,9 +280,54 @@ def main():
     fig, ax = plt.subplots()
     ax.add_collection(p)
     ax.set_aspect('equal')
-    plt.plot(x,y,out[0],out[1],px,py,qx,qy,sx,sy,tx,ty)
+    plt.imshow(np.asarray(im1),extent=[0,width,0,height])
-    #plt.plot(tx,ty)
+    plt.imshow(np.asarray(im2),extent=[bbox[0],bbox[2],bbox[1],bbox[3]])
+    plt.plot(x,y,out[0],out[1],px,py,qx,qy,tx,ty)
     plt.show()
+    iw,ih = im2.size
+    data = []
+    for i,(off,h,(p0,p1,p2,p3)) in enumerate(zip(offs,heights,quads)):
+        # first, account for the offset of the input image
+        p0 = p0[0]-bbox[0],p0[1]-bbox[1]
+        p1 = p1[0]-bbox[0],p1[1]-bbox[1]
+        p2 = p2[0]-bbox[0],p2[1]-bbox[1]
+        p3 = p3[0]-bbox[0],p3[1]-bbox[1]
+        # PIL expects coordinates in counter clockwise order
+        p1,p3 = p3,p1
+        #   x          lon
+        # ----- =     -----
+        #   w     bbox[2]-bbox[0]
+        # translate to pixel coordinates
+        p0 = (iw*p0[0])/(bbox[2]-bbox[0]),(ih*p0[1])/(bbox[3]-bbox[1])
+        p1 = (iw*p1[0])/(bbox[2]-bbox[0]),(ih*p1[1])/(bbox[3]-bbox[1])
+        p2 = (iw*p2[0])/(bbox[2]-bbox[0]),(ih*p2[1])/(bbox[3]-bbox[1])
+        p3 = (iw*p3[0])/(bbox[2]-bbox[0]),(ih*p3[1])/(bbox[3]-bbox[1])
+        # PIL starts coordinate system at the upper left corner, swap y coord
+        p0 = int(p0[0]),int(ih-p0[1])
+        p1 = int(p1[0]),int(ih-p1[1])
+        p2 = int(p2[0]),int(ih-p2[1])
+        p3 = int(p3[0]),int(ih-p3[1])
+        box=(0,int(ih*(height-off-h)/(bbox[3]-bbox[1])),
+                int(iw*width/(bbox[2]-bbox[0])),int(ih*(height-off)/(bbox[3]-bbox[1])))
+        quad=(p0[0],p0[1],p1[0],p1[1],p2[0],p2[1],p3[0],p3[1])
+        data.append((box,quad))
+    im_out = im2.transform((int(iw*width/(bbox[2]-bbox[0])),int(ih*height/(bbox[3]-bbox[1]))),Image.MESH,data,Image.BICUBIC)
+    im_out.save("out.png")
+    return True
 
 if __name__ == '__main__':
-    main()
+    x = []
+    y = []
+    with open(sys.argv[1]) as f:
+        for l in f:
+            a,b = l.split()
+            # apply mercator projection
+            b = lat2y(float(b))
+            x.append(float(a))
+            y.append(b)
+    width = 2.0/7.0
+    main(x,y,width,6,20)
+    #for smoothing in [1,2,4,8,12]:
+    #    for subdiv in range(10,30):
+    #        if main(x,y,width,smoothing,subdiv):
+    #            print width,smoothing,subdiv