diff --git a/fitbspline.py b/fitbspline.py index dbe3814..7bc128d 100644 --- a/fitbspline.py +++ b/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: - # return ptInTriangle(p, p0, cross, p3) or ptInTriangle(p, p2, cross, p1) - #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) + if cross: + return ptInTriangle(p, p0, cross, p3) or ptInTriangle(p, p2, cross, p1) + else: + 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): - 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 +def main(x,y,width,smoothing,subdiv): halfwidth = width/2.0 - x = [] - y = [] - 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) + tck,u = interpolate.splprep([x,y],s=smoothing) + unew = np.linspace(0,1.0,subdiv+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: - print found - containingquad.append(found) + if len(found) > 1: + print "point found in two quads" + return None + containingquad.append(found[0]) else: - print "can't find quad for point" containingquad.append(None) - #exit(1) - print containingquad - trans = [] - print width, height - srcquads = [] - for off,h,srcquad in zip(offs,heights,quads): - #targetquad = ((0,height-off),(width,height-off),(width,height-off-h),(0,height-off-h)) + # check if the only points for which no quad could be found are in the + # beginning or in the end + # find the first missing ones: + for i,q in enumerate(containingquad): + if q != None: + break + # 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 - for (rx,ry),l in zip(zip(x,y),containingquad): - if not l: + 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),i in zip(zip(x,y),containingquad): + if i == None: continue - for i in l[:1]: - (ax,ay),(bx,by),(cx,cy),(dx,dy) = quads[i] - s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,rx,ry) - if len(s) != 1 or len(t) != 1: - 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) + (ax,ay),(bx,by),(cx,cy),(dx,dy) = quads[i] + s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,rx,ry) + # if more than one solution, take second + # TODO: investigate if this is always the right solution if len(s) != 1 or len(t) != 1: - print "fail" - exit(1) - #u = (a*ax + b*ay + c)/(g*ax + h*ay + 1) - #v = (d*ax + e*ay + f)/(g*ax + h*ay + 1) - u = s[0]*width - v = off+t[0]*h - sx.append(u) - sy.append(v) - #u = (a*bx + b*by + c)/(g*bx + h*by + 1) - #v = (d*bx + e*by + f)/(g*bx + h*by + 1) - s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,x2,y2) - if len(s) != 1 or len(t) != 1: - print "fail" - exit(1) - u = s[0]*width - v = off+t[0]*h - sx.append(u) - sy.append(v) + s = s[1] + t = t[1] + else: + s = s[0] + t = t[0] + u = s*width + v = offs[i]+t*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: + # return None + # u = s[0]*width + # v = off+t[0]*h + # 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.plot(tx,ty) + plt.imshow(np.asarray(im1),extent=[0,width,0,height]) + 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