1. # -*- coding: utf-8 -*-
2. """
3. Created on Thu Jun  6 10:21:54 2019
4.  
5. @author: Student
6. """
7.  
8. import numpy as np
9. from matplotlib import pyplot as plt
10. from skimage import io
11. from scipy import ndimage as nd
12. from scipy import interpolate
13. from scipy import linalg
14. from skimage import color
15. import time
16.  
17.  
18. class LineBuilder:
19.     def __init__(self, line):
20.         self.line = line
21.         self.xs = list(line.get_xdata())
22.         self.ys = list(line.get_ydata())
23.         self.cid = line.figure.canvas.mpl_connect('button_press_event', self)
24.         self.cid = line.figure.canvas.mpl_connect('key_press_event', self)
25.         self.cid = line.figure.canvas.mpl_connect('scroll_event', self)
26.  
27.     def __call__(self, event):
28.         if(event.name=='key_press_event'):
29.             self.xs.append(int(self.xs[0]))
30.             self.ys.append(int(self.ys[0]))
31.             self.line.set_data(self.xs, self.ys)
32.             self.line.figure.canvas.draw()
33.             self.line.figure.canvas.mpl_disconnect(self.cid)
34.             print("xs = ", self.xs)
35.             print("ys = ", self.ys)
36.         if (event.name == 'button_press_event'):
37.             self.xs.append(int(event.xdata))
38.             self.ys.append(int(event.ydata))
39.             self.line.set_data(self.xs, self.ys)
40.             self.line.figure.canvas.draw()
41.             print("xs = ", self.xs)
42.             print("ys = ", self.ys)
43.         if(event.name=='scroll_event'):
44.             self.xs.pop()
45.             self.ys.pop()
46.             self.line.set_data(self.xs, self.ys)
47.             print("xs = ", self.xs)
48.             print("ys = ", self.ys)
49.        
50.  
51. def click_points(image):
52.     fig = plt.figure()
53.     ax = fig.add_subplot(111)            
54.     ax.imshow(image,cmap='gray')
55.     ax.set_title('click')
56.     line, = ax.plot([], [])
57.     linebuilder = LineBuilder(line)
58.     plt.show()
59.     plt.pause(7)
60.     xs=np.asarray(linebuilder.xs)
61.     ys=np.asarray(linebuilder.ys)
62.     return xs, ys
63.  
64.            
65. def reinterpolate_contours(xs, ys, dmin, dmax):
66.    
67.     n_xs = xs
68.     n_ys = ys
69.     grid_x, grid_y = np.meshgrid(n_xs, n_ys)
70.     a = n_ys.shape[0]
71.     b = n_xs.shape[0]
72.    
73.     for i in range(n_xs.shape[0]):
74.         if np.sqrt((n_xs[i+1] - n_ys[i+1])**2 + ((n_xs[i] - n_ys[i])**2)) > dmax:
75.             for j in range(grid_x.shape[1]):
76.                 I1 = b-n_xs[j]/b - grid_x[0,0] + n_xs[j]/b * grid_x[0,b]
77.                 I2 = b-n_xs[j]/b - grid_x[a,0] + n_xs[j]/b* grid_x[a,b]
78.                 I3 = a-n_ys[j]/a * I1 + n_ys[j]/a * I2
79.             grid_y[j], grid_x[j] = I3
80.    
81.     return n_xs, n_ys
82.  
83.  
84.  
85.  
86. def calculate_gradient_at_contour_points(image, sigma, xs, ys):
87.     fimage = nd.gaussian_filter(image, sigma)
88.     gx, gy = np.gradient(fimage)
89.     mag = (gx**2+gy**2)**0.5
90.     mag = (mag - mag.min()/(mag.max()-mag.min()))
91.     xs_gradient, ys_gradient = np.gradient(mag)
92.    
93.     return xs_gradient, ys_gradient
94.  
95. def fx(xs, ys):
96.        
97.     xs[ xs < 0 ] = 0.
98.     ys[ ys < 0 ] = 0.
99.  
100.     xs[ xs > image.shape[1]-1 ] = image.shape[1]-1
101.     ys[ ys > image.shape[0]-1 ] = image.shape[0]-1
102.  
103.     return xs_gradient[ (y.round().astype(int), x.round().astype(int)) ]
104.  
105. def fy(x, y):
106.        
107.     x[ x < 0 ] = 0.
108.     y[ y < 0 ] = 0.
109.    
110.     x[ x > image.shape[1]-1 ] = image.shape[1]-1
111.     y[ y > image.shape[0]-1 ] = image.shape[0]-1
112.            
113.     return ys_gradient[ (y.round().astype(int), x.round().astype(int)) ]
114.  
115.        
116.  
117.  
118. def generate_coefficient_matrix(N, alpha, beta):
119.     def create_matrix_A(a, b, N):
120.         row = np.r_[
121.         -2*a - 6*b,
122.         a + 4*b,
123.         -b,
124.         np.zeros(N-5),
125.         -b,
126.         a + 4*b
127.     ]
128.     A = np.zeros((N,N))
129.     for i in range(N):
130.         A[i] = np.roll(row, i)
131.         return A
132.  
133. def run():
134.    
135.  
136.     image = io.imread('xray.jpg')
137.     image = color.rgb2gray(image)
138.    
139.     plt.show()
140.     xs,ys = click_points(image)
141.     time.sleep(1)
142.     n_xs, n_ys = reinterpolate_contours(xs, ys, 5, 1)
143.     xs_gradient, ys_gradient = calculate_gradient_at_contour_points(image, 1, xs, ys)
144.    
145.  
146.     plt.figure(2)
147.     plt.imshow(image)
148.    
149. if __name__ == "__main__":
150.         run()