Advertisement
Guest User

Untitled

a guest
Apr 26th, 2019
81
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 9.68 KB | None | 0 0
  1. import matplotlib.pyplot as plt
  2. import numpy as np
  3. import random
  4. import cv2 as cv
  5. from scipy.interpolate import RectBivariateSpline
  6.  
  7. def plot_image(img1, img2, color):
  8. fig = plt.figure()
  9. ax1 = fig.add_subplot(1, 2, 1)
  10. ax1.grid(False)
  11.  
  12. ax2 = fig.add_subplot(1, 2, 2)
  13. ax2.grid(False)
  14.  
  15. if (color == 'colored'):
  16. ax1.imshow(img1)
  17. ax2.imshow(img2)
  18. elif (color == "gray"):
  19. ax1.imshow(img1, cmap='gray', interpolation="bicubic")
  20. ax2.imshow(img2, cmap='gray', interpolation="bicubic")
  21.  
  22. plt.show()
  23.  
  24. def BGR2RGB(image):
  25. return cv.cvtColor(image,cv.COLOR_BGR2RGB)
  26.  
  27.  
  28. def read_image(name1, name2):
  29. img1 = cv.imread(name1)
  30. img2 = cv.imread(name2)
  31.  
  32. RGB_image1 = BGR2RGB(img1)
  33. RGB_image2 = BGR2RGB(img2)
  34.  
  35. # plot_image(RGB_image1, RGB_image2, 'colored')
  36.  
  37. return img1, img2
  38.  
  39. def calculate_h_matrix(points):
  40. n = len(points[0])
  41. b = np.zeros((n*2,1))
  42. #build b vector on the form : [x0,y0,x1,y1,x2,y2,x3,y3]
  43. for i in range(0,n):
  44. b[2*i] = points[1][i][0]
  45. b[2*i+1] = points[1][i][1]
  46. A = build_matrix_a(points,n)
  47. # print("A matrix:")
  48. # print(A)
  49. #solve for Ax=b
  50. h = np.linalg.lstsq(A, b)[0];
  51. #rounding h vector
  52. # h = np.around(h)
  53. # print(h)
  54. return h
  55.  
  56. def build_matrix_a(points,n):
  57. A = np.zeros(((2 * n), 8))
  58. for i in range(n):
  59. x = points[0][i][0]
  60. y = points[0][i][1]
  61. x_dash = points[1][i][0]
  62. y_dash = points[1][i][1]
  63. # print(x,y,x_dash,y_dash)
  64. A[2*i] = np.array([x,y,1,0,0,0,-x*x_dash,-y*x_dash])
  65. A[2*i +1] = np.array([0,0,0,x,y,1,-x*y_dash,-y*y_dash])
  66. return A
  67.  
  68.  
  69. def h_to_matrix(h):
  70. H = np.zeros((3, 3))
  71. for i in range(0, 3):
  72. for j in range(0, 3):
  73. if 3 * i + j < 8:
  74. H[i][j] = h[3 * i + j]
  75.  
  76. H[2][2] = 1
  77. return H
  78.  
  79.  
  80. def transform_point(H, point):
  81. output_point = np.dot(H, point)
  82. output_point[0] = output_point[0] / output_point[2]
  83. output_point[1] = output_point[1] / output_point[2]
  84. output_point[2] = 1
  85. # print(output_point)
  86. return output_point
  87.  
  88.  
  89. image1, image2 = read_image('img1.jpg', 'img2.jpg')
  90.  
  91. '''
  92. fig = plt.figure()
  93.  
  94.  
  95. figA = fig.add_subplot(1,1,1)
  96. figA.imshow(image1)
  97. ptsa = plt.ginput(n=10, timeout=100)
  98.  
  99. figA = fig.add_subplot(1,1,1)
  100. figA.imshow(image2)
  101. ptsb = plt.ginput(n=10, timeout=100)
  102.  
  103. # figB = fig.add_subplot(1,2,2)
  104.  
  105. # Display the image
  106. # lower use to flip the image
  107. # figA.imshow(image1)#,origin='lower')
  108. # figB.imshow(image2)#,origin='lower')
  109. plt.axis('image')
  110. plt.show()
  111. # n = number of points to read
  112. # pts = plt.ginput(n=3, timeout=10)
  113. print(ptsa);
  114. print(ptsb);
  115. '''
  116.  
  117. def get_warp_dimension(image, H):
  118. height = image.shape[0]
  119. width = image.shape[1]
  120.  
  121. boundaries = np.array([[0, 0], [height - 1, 0], [0, width - 1], [height - 1, width - 1]])
  122. i_s = np.zeros(4, dtype=int)
  123. j_s = np.zeros(4, dtype=int)
  124.  
  125. for k in range(0, 4):
  126. i = boundaries[k][0]
  127. j = boundaries[k][1]
  128. mapped = transform_point(H,np.array([[j], [i], [1]]))
  129. # print(mapped.shape)
  130. j_s[k] = int(mapped[0][0])
  131. i_s[k] = int(mapped[1][0])
  132.  
  133. #i -> col
  134. #j -> row
  135. min_mapped_i = np.min(i_s)
  136. min_mapped_j = np.min(j_s)
  137. max_mapped_i = np.max(i_s)
  138. max_mapped_j = np.max(j_s)
  139.  
  140. new_height = (max_mapped_i - min_mapped_i + 1)
  141. new_width = (max_mapped_j - min_mapped_j + 1)
  142.  
  143. height_shift = -min_mapped_i
  144. width_shift = -min_mapped_j
  145.  
  146. return (new_height, new_width, height_shift, width_shift)
  147.  
  148.  
  149. def warp(image, H):
  150. height = image.shape[0]
  151. width = image.shape[1]
  152.  
  153. (wraped_image_height, wraped_image_width, height_shift, width_shift) = get_warp_dimension(image, H)
  154. wraped_image = np.zeros((wraped_image_height, wraped_image_width, 3), dtype=np.uint8);
  155.  
  156. for i in range(0, image.shape[0]):
  157. for j in range(0, image.shape[1]):
  158. mapped = transform_point(H,np.array([[j], [i], [1]]))
  159. mapped_j = int(mapped[0][0])
  160. mapped_i = int(mapped[1][0])
  161. wraped_image[mapped_i + height_shift][mapped_j + width_shift] = image[i][j]
  162.  
  163. return wraped_image
  164.  
  165. def inv_transform(point, H):
  166. H_inv = np.linalg.inv(H)
  167. result = np.dot(H_inv, point)
  168. w=result[2]
  169. result[0] = result[0]/w
  170. result[1] = result[1]/w
  171. result[2] = 1
  172.  
  173. return result
  174.  
  175. def remove_holes(img, warped_img, H, height_shift, width_shift):
  176. warped_height = warped_img.shape[0]
  177. warped_width = warped_img.shape[1]
  178.  
  179. img_height = img.shape[0]
  180. img_width = img.shape[1]
  181.  
  182. b, g, r = cv.split(img)
  183. bi = RectBivariateSpline(np.arange(img_height), np.arange(img_width), b)
  184. gi = RectBivariateSpline(np.arange(img_height), np.arange(img_width), g)
  185. ri = RectBivariateSpline(np.arange(img_height), np.arange(img_width), r)
  186.  
  187. for i in range(0, warped_height):
  188. for j in range(0, warped_width):
  189. if int(warped_img[i][j][0]) == 0 and int(warped_img[i][j][1]) == 0 and int(warped_img[i][j][2]) == 0:
  190. inv_warped_img = inv_transform(np.array([[(j - width_shift)], [(i - height_shift)], [1]]), H)
  191. inv_i = int(inv_warped_img[1][0])
  192. inv_j = int(inv_warped_img[0][0])
  193.  
  194. warped_img[i][j][0] = bi.ev(inv_i, inv_j)
  195. warped_img[i][j][1] = gi.ev(inv_i, inv_j)
  196. warped_img[i][j][2] = ri.ev(inv_i, inv_j)
  197.  
  198. return warped_img
  199.  
  200.  
  201. #manual matched points
  202. points = np.zeros((2,10,2))
  203. before = np.column_stack((np.array([565.17,639.5,567.24,635.37,765.43,730.33,812.91,866.59,893.43,980.14])
  204. ,np.array([539.19,539.193,555.70,559.83,491.70, 270.80,322.419,336.87,322.41,332.74])))
  205. after = np.column_stack((np.array([102.72,177.04,100.66,174.98,305.04,288.53,362.85,422.72,441.30,521.822]),
  206. np.array([504.096,504.096,522.677,522.67,458.67,237.77,293.51,307.96, 299.70,314.16])))
  207. points[0] = before
  208. points[1] = after
  209.  
  210.  
  211.  
  212.  
  213. # print(points)
  214.  
  215. h = calculate_h_matrix(points)
  216. # print(h)
  217.  
  218. H = h_to_matrix(h)
  219. '''
  220. # print(get_warp_dimension(image1, H))
  221. warped_image_height, warped_image_width, height_shift, width_shift = get_warp_dimension(image1, H)
  222.  
  223. warped_image = warp(image1, h_to_matrix(h))
  224. print(len(warped_image))
  225. print(len(warped_image[0]))
  226.  
  227. fig = plt.figure()
  228. figA = fig.add_subplot(1,1,1)
  229. figA.imshow(warped_image)
  230. plt.show()
  231. # figA.imshow()
  232. #
  233. #
  234. # plt.imshow(warp(image1, h_to_matrix(h)))
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242. warped_img_without_holes = remove_holes(image1, warped_image, H, height_shift, width_shift)
  243. fig = plt.figure()
  244. figA = fig.add_subplot(1,1,1)
  245. figA.imshow(warped_img_without_holes)
  246.  
  247. plt.show()
  248.  
  249.  
  250.  
  251.  
  252. ref_image_height = image2.shape[0]
  253. ref_image_width = image2.shape[1]
  254.  
  255. print('Ref-image size ' + str(image2.shape))
  256.  
  257.  
  258. mergedImage_height = ref_image_height + height_shift
  259. if warped_image_height > mergedImage_height:
  260. mergedImage_height = warped_image_height
  261.  
  262. mergedImage_width = ref_image_width + width_shift
  263. if warped_image_width > mergedImage_width:
  264. mergedImage_width = mergedImage_width
  265.  
  266. mergedImage = np.zeros((mergedImage_height , mergedImage_width , 3), dtype=np.uint8);
  267.  
  268. # sketch the reference image
  269. for i in range(0,ref_image_height):
  270. for j in range(0, ref_image_width):
  271. # if not( i+ shiftHeight < newHeight and j + shiftWidth < newWidth):
  272. mergedImage[i + height_shift][j + width_shift] = image2[i][j]
  273.  
  274. fig = plt.figure()
  275. figA = fig.add_subplot(1, 1, 1)
  276. figA.imshow(mergedImage)
  277.  
  278. plt.show()
  279.  
  280. print(warped_img_without_holes.shape[0],warped_img_without_holes.shape[1])
  281. print('shape : \n')
  282. print(warped_img_without_holes.shape)
  283.  
  284. # sketch the destination image
  285. print(mergedImage_height,mergedImage_width)
  286. for i in range(0,mergedImage_height-1):
  287. for j in range(0, mergedImage_width-1):
  288. if j>=warped_img_without_holes.shape[0] or i >= warped_img_without_holes.shape[0]:
  289. continue
  290. if not(int(warped_img_without_holes[i][j][0]) == 0 and int(warped_img_without_holes[i][j][0]) == 0 and int(warped_img_without_holes[i][j][2]) == 0 ):
  291. mergedImage[i][j] = warped_img_without_holes[i][j]
  292.  
  293. # im = Image.fromarray(mergedImage)
  294. # im.save("after_add_ref_image.jpg")
  295. #
  296.  
  297. fig = plt.figure()
  298. figA = fig.add_subplot(1,1,1)
  299. figA.imshow(mergedImage)
  300.  
  301. plt.show()
  302.  
  303.  
  304. '''
  305.  
  306. def verifyH(imageA,imageB,h):
  307.  
  308. fig = plt.figure()
  309. figA = fig.add_subplot(1, 2, 1)
  310. figB = fig.add_subplot(1, 2, 2)
  311. figA.imshow(image1) # ,origin='lower')
  312. figB.imshow(image2)
  313. ptsa = plt.ginput(n=4, timeout=10)
  314. print((ptsa))
  315. print(ptsa[0])
  316.  
  317. figA.scatter(x=[ptsa[0][0],ptsa[1][0],ptsa[2][0],ptsa[3][0]],y=[ptsa[0][1],ptsa[1][1],ptsa[2][1],ptsa[3][1]],c='r', s=40)
  318. for i in range(4):
  319. tran_point = transform_point(H,np.array([ptsa[i][0],ptsa[i][1],1]))
  320. figB.scatter(tran_point[0], tran_point[1], c='b', s=40)
  321. plt.show()
  322.  
  323. verifyH(image1,image2,H)
  324. #10 points
  325.  
  326. # [(565.1774193548388, 539.1935483870968), (639.5, 539.1935483870968), (567.241935483871, 555.7096774193549), (635.3709677419355, 559.8387096774193), (765.4354838709678, 491.7096774193548), (730.3387096774194, 270.8064516129033), (812.9193548387096, 322.41935483870975), (866.5967741935484, 336.8709677419355), (893.4354838709677, 322.41935483870975), (980.1451612903226, 332.741935483871)]
  327. # [(102.7258064516129, 504.09677419354836), (177.04838709677418, 504.09677419354836), (100.66129032258064, 522.6774193548387), (174.98387096774192, 522.6774193548387), (305.04838709677415, 458.6774193548387), (288.5322580645161, 237.7741935483872), (362.85483870967744, 293.51612903225805), (422.7258064516129, 307.9677419354839), (441.3064516129033, 299.7096774193549), (521.8225806451613, 314.1612903225807)]
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement