threaded error surface plotting
[imago.git] / gridf.py
1 import Image, ImageDraw, ImageFilter
2
3 from manual import lines as g_grid, l2ad, intersection, line as g_line
4 from intrsc import intersections_from_angl_dist
5 from linef import line_from_angl_dist
6 import pcf
7
8 class GridFittingFailedError(Exception):
9     pass
10
11 class MyGaussianBlur(ImageFilter.Filter):
12     name = "GaussianBlur"
13
14     def __init__(self, radius=2):
15         self.radius = radius
16     def filter(self, image):
17         return image.gaussian_blur(self.radius)
18
19 class V(object):
20     def __init__(self, x, y):
21         self.x = x
22         self.y = y
23     
24     def __add__(self, other):
25         return V(self.x + other.x, self.y + other.y)
26
27     def __sub__(self, other):
28         return V(self.x - other.x, self.y - other.y)
29
30     def __rmul__(self, other):
31         return V(other * self.x, other * self.y)
32
33     def __len__(self):
34         return 2;
35
36     def __getitem__(self, key):
37         if key == 0:
38             return self.x
39         elif key == 1:
40             return self.y
41         elif type(key) != int:
42             raise TypeError("V indices must be integers") 
43         else:
44             raise KeyError("V index ({}) out of range".format(key))
45
46     def __iter__(self):
47         yield self.x
48         yield self.y
49
50     @property
51     def normal(self):
52         return V(-self.y, self.x)
53
54 def projection(point, line, vector):
55     return V(*intersection(g_line(point, point + vector.normal), g_line(*line)))
56     
57 def error_surface(lines, a, b, c, d, hough, size, v1):
58     import matplotlib.pyplot as plt
59     from matplotlib import cm
60     import threading
61     import Queue
62     import time
63     import sys
64     import pickle
65
66     class Worker(threading.Thread):
67         def __init__(self, q_in, q_out, job):
68             threading.Thread.__init__(self)
69             self.q_in = q_in
70             self.q_out = q_out
71             self.job = job
72
73         def run(self):
74             while True:
75                 x = self.q_in.get()
76                 try:
77                     self.q_out.put((x, self.job(x)))
78                 except Exception:
79                     pass
80                 print x
81                 self.q_in.task_done()
82
83     X = []
84     Y = []
85     Z = []
86     s = 0.001
87     k = 200
88     for i in range(-k, k):
89         X.append(range(-k, k))
90         Y.append(2*k*[i])
91
92     job = lambda x: [distance(lines, get_grid(a + X[x][y] * s * v1, 
93                                               b + Y[x][y] * s * v1, 
94                                               c, d, hough, size),
95                                 size) for y in range(0,2 * k)]
96
97     q_in = Queue.Queue()
98     q_out = Queue.Queue()
99     for i in range(4):
100         t = Worker(q_in, q_out, job)
101         t.daemon = True
102         t.start()
103         
104     start = time.time()
105     for x in range(0, 2*k):
106         q_in.put(x)
107
108     q_in.join()
109
110     print time.time() - start
111
112     while True:
113         try:
114             Z.append(q_out.get_nowait())
115         except Queue.Empty:
116             break
117
118     Z.sort()
119     Z = [t for (x, t) in Z]
120
121     s_file = open('surface' + str(k), 'w')
122     pickle.dump((X, Y, Z), s_file)
123     s_file.close()
124     plt.imshow(Z, cmap=cm.gnuplot2, interpolation='bicubic', 
125                origin='upper', extent=(-k, k, -k, k), aspect='equal')
126     plt.colorbar()
127
128     plt.show()
129
130     sys.exit()
131
132 def find(lines, size, l1, l2, bounds, hough, do_something):
133     a, b, c, d = [V(*a) for a in bounds]
134     l1 = line_from_angl_dist(l1, size)
135     l2 = line_from_angl_dist(l2, size)
136     v1 = V(*l1[0]) - V(*l1[1])
137     v2 = V(*l2[0]) - V(*l2[1])
138     a = projection(a, l1, v1) 
139     b = projection(b, l1, v1) 
140     c = projection(c, l2, v2) 
141     d = projection(d, l2, v2) 
142
143     #error_surface(lines, a, b, c, d, hough, size, v1)
144
145     grid = get_grid(a, b, c, d, hough, size)
146     dist = distance(lines, grid, size)
147     print dist
148    
149     s = 0.02
150     while True:
151         ts1 = [(s, 0), (-s, 0), (s, s), (-s, -s), (-s, s), (s, -s), (0, s),  (0, -s)]
152         grids = [(get_grid(a + t[0] * v1, b + t[1] * v1, 
153                            c, d, hough, size), t) for t in ts1]
154         distances = [(distance(lines, grid, size), 
155                       grid, t) for grid, t in grids]
156         distances.sort(reverse=True)
157         if distances[0][0] > dist:
158             dist = distances[0][0]
159             grid = distances[0][1]
160             t = distances[0][2]
161             a, b = a + t[0] * v1, b + t[1] * v1
162             print dist
163             s *= 0.75
164         else: 
165            break
166
167     print "---"
168
169     s = 0.02
170     while True:
171         ts1 = [(s, 0), (-s, 0), (s, s), (-s, -s), (-s, s), (s, -s), (0, s),  (0, -s)]
172         grids = [(get_grid(a, b, 
173                            c + t[0] * v2, d + t[1] * v2, hough, size), t) for t in ts1]
174         distances = [(distance(lines, grid, size), 
175                       grid, t) for grid, t in grids]
176         distances.sort(reverse=True)
177         if distances[0][0] > dist:
178             dist = distances[0][0]
179             grid = distances[0][1]
180             t = distances[0][2]
181             c, d = c + t[0] * v2, d + t[1] * v2
182             print dist
183             s *= 0.75
184         else:
185             break
186
187     grid_lines = [[l2ad(l, size) for l in grid[0]], [l2ad(l, size) for l in grid[1]]]
188     return grid, grid_lines
189
190 def get_grid(a, b, c, d, hough, size):
191     l1 = hough.lines_from_list([a, b])
192     l2 = hough.lines_from_list([c, d])
193     c = intersections_from_angl_dist([l1, l2], size, get_all=True)
194     #TODO do something when a corner is outside the image
195     corners = (c[0] + c[1])
196     if len(corners) < 4:
197         print l1, l2, c
198         raise GridFittingFailedError
199     grid = g_grid(corners)
200     return grid
201
202 def distance(lines, grid, size):
203     im_l = Image.new('L', size)
204     dr_l = ImageDraw.Draw(im_l)
205     for line in sum(lines, []):
206         dr_l.line(line_from_angl_dist(line, size), width=1, fill=255)
207     im_l = im_l.filter(MyGaussianBlur(radius=15))
208     #GaussianBlur is undocumented class, may not work in future versions of PIL
209     im_g = Image.new('L', size)
210     dr_g = ImageDraw.Draw(im_g)
211     for line in grid[0] + grid[1]:
212         dr_g.line(line, width=1, fill=255)
213     #im_g = im_g.filter(MyGaussianBlur(radius=3))
214     #im_d, distance = combine(im_l, im_g)
215     distance = pcf.combine(im_l.tostring(), im_g.tostring())
216     return distance
217
218 def combine(bg, fg):
219     bg_l = bg.load()
220     fg_l = fg.load()
221     #res = Image.new('L', fg.size)
222     #res_l = res.load()
223
224     score = 0
225     area = 0
226
227     for x in xrange(fg.size[0]):
228         for y in xrange(fg.size[1]):
229             if fg_l[x, y]:
230                 #res_l[x, y] = bg_l[x, y] * fg_l[x, y]
231                 score +=  bg_l[x, y]
232                 area += 1
233
234     #return res, float(score)/area
235     return None, float(score)/area