Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import matplotlib.pyplot as plt
- import numpy as np
- import random
- import cv2 as cv
- from scipy.interpolate import RectBivariateSpline
- def plot_image(img1, img2, color):
- fig = plt.figure()
- ax1 = fig.add_subplot(1, 2, 1)
- ax1.grid(False)
- ax2 = fig.add_subplot(1, 2, 2)
- ax2.grid(False)
- if (color == 'colored'):
- ax1.imshow(img1)
- ax2.imshow(img2)
- elif (color == "gray"):
- ax1.imshow(img1, cmap='gray', interpolation="bicubic")
- ax2.imshow(img2, cmap='gray', interpolation="bicubic")
- plt.show()
- def BGR2RGB(image):
- return cv.cvtColor(image,cv.COLOR_BGR2RGB)
- def read_image(name1, name2):
- img1 = cv.imread(name1)
- img2 = cv.imread(name2)
- RGB_image1 = BGR2RGB(img1)
- RGB_image2 = BGR2RGB(img2)
- # plot_image(RGB_image1, RGB_image2, 'colored')
- return img1, img2
- def calculate_h_matrix(points):
- n = len(points[0])
- b = np.zeros((n*2,1))
- #build b vector on the form : [x0,y0,x1,y1,x2,y2,x3,y3]
- for i in range(0,n):
- b[2*i] = points[1][i][0]
- b[2*i+1] = points[1][i][1]
- A = build_matrix_a(points,n)
- # print("A matrix:")
- # print(A)
- #solve for Ax=b
- h = np.linalg.lstsq(A, b)[0];
- #rounding h vector
- # h = np.around(h)
- # print(h)
- return h
- def build_matrix_a(points,n):
- A = np.zeros(((2 * n), 8))
- for i in range(n):
- x = points[0][i][0]
- y = points[0][i][1]
- x_dash = points[1][i][0]
- y_dash = points[1][i][1]
- # print(x,y,x_dash,y_dash)
- A[2*i] = np.array([x,y,1,0,0,0,-x*x_dash,-y*x_dash])
- A[2*i +1] = np.array([0,0,0,x,y,1,-x*y_dash,-y*y_dash])
- return A
- def h_to_matrix(h):
- H = np.zeros((3, 3))
- for i in range(0, 3):
- for j in range(0, 3):
- if 3 * i + j < 8:
- H[i][j] = h[3 * i + j]
- H[2][2] = 1
- return H
- def transform_point(H, point):
- output_point = np.dot(H, point)
- output_point[0] = output_point[0] / output_point[2]
- output_point[1] = output_point[1] / output_point[2]
- output_point[2] = 1
- # print(output_point)
- return output_point
- image1, image2 = read_image('img1.jpg', 'img2.jpg')
- '''
- fig = plt.figure()
- figA = fig.add_subplot(1,1,1)
- figA.imshow(image1)
- ptsa = plt.ginput(n=10, timeout=100)
- figA = fig.add_subplot(1,1,1)
- figA.imshow(image2)
- ptsb = plt.ginput(n=10, timeout=100)
- # figB = fig.add_subplot(1,2,2)
- # Display the image
- # lower use to flip the image
- # figA.imshow(image1)#,origin='lower')
- # figB.imshow(image2)#,origin='lower')
- plt.axis('image')
- plt.show()
- # n = number of points to read
- # pts = plt.ginput(n=3, timeout=10)
- print(ptsa);
- print(ptsb);
- '''
- def get_warp_dimension(image, H):
- height = image.shape[0]
- width = image.shape[1]
- boundaries = np.array([[0, 0], [height - 1, 0], [0, width - 1], [height - 1, width - 1]])
- i_s = np.zeros(4, dtype=int)
- j_s = np.zeros(4, dtype=int)
- for k in range(0, 4):
- i = boundaries[k][0]
- j = boundaries[k][1]
- mapped = transform_point(H,np.array([[j], [i], [1]]))
- # print(mapped.shape)
- j_s[k] = int(mapped[0][0])
- i_s[k] = int(mapped[1][0])
- #i -> col
- #j -> row
- min_mapped_i = np.min(i_s)
- min_mapped_j = np.min(j_s)
- max_mapped_i = np.max(i_s)
- max_mapped_j = np.max(j_s)
- new_height = (max_mapped_i - min_mapped_i + 1)
- new_width = (max_mapped_j - min_mapped_j + 1)
- height_shift = -min_mapped_i
- width_shift = -min_mapped_j
- return (new_height, new_width, height_shift, width_shift)
- def warp(image, H):
- height = image.shape[0]
- width = image.shape[1]
- (wraped_image_height, wraped_image_width, height_shift, width_shift) = get_warp_dimension(image, H)
- wraped_image = np.zeros((wraped_image_height, wraped_image_width, 3), dtype=np.uint8);
- for i in range(0, image.shape[0]):
- for j in range(0, image.shape[1]):
- mapped = transform_point(H,np.array([[j], [i], [1]]))
- mapped_j = int(mapped[0][0])
- mapped_i = int(mapped[1][0])
- wraped_image[mapped_i + height_shift][mapped_j + width_shift] = image[i][j]
- return wraped_image
- def inv_transform(point, H):
- H_inv = np.linalg.inv(H)
- result = np.dot(H_inv, point)
- w=result[2]
- result[0] = result[0]/w
- result[1] = result[1]/w
- result[2] = 1
- return result
- def remove_holes(img, warped_img, H, height_shift, width_shift):
- warped_height = warped_img.shape[0]
- warped_width = warped_img.shape[1]
- img_height = img.shape[0]
- img_width = img.shape[1]
- b, g, r = cv.split(img)
- bi = RectBivariateSpline(np.arange(img_height), np.arange(img_width), b)
- gi = RectBivariateSpline(np.arange(img_height), np.arange(img_width), g)
- ri = RectBivariateSpline(np.arange(img_height), np.arange(img_width), r)
- for i in range(0, warped_height):
- for j in range(0, warped_width):
- if int(warped_img[i][j][0]) == 0 and int(warped_img[i][j][1]) == 0 and int(warped_img[i][j][2]) == 0:
- inv_warped_img = inv_transform(np.array([[(j - width_shift)], [(i - height_shift)], [1]]), H)
- inv_i = int(inv_warped_img[1][0])
- inv_j = int(inv_warped_img[0][0])
- warped_img[i][j][0] = bi.ev(inv_i, inv_j)
- warped_img[i][j][1] = gi.ev(inv_i, inv_j)
- warped_img[i][j][2] = ri.ev(inv_i, inv_j)
- return warped_img
- #manual matched points
- points = np.zeros((2,10,2))
- 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])
- ,np.array([539.19,539.193,555.70,559.83,491.70, 270.80,322.419,336.87,322.41,332.74])))
- 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]),
- np.array([504.096,504.096,522.677,522.67,458.67,237.77,293.51,307.96, 299.70,314.16])))
- points[0] = before
- points[1] = after
- # print(points)
- h = calculate_h_matrix(points)
- # print(h)
- H = h_to_matrix(h)
- '''
- # print(get_warp_dimension(image1, H))
- warped_image_height, warped_image_width, height_shift, width_shift = get_warp_dimension(image1, H)
- warped_image = warp(image1, h_to_matrix(h))
- print(len(warped_image))
- print(len(warped_image[0]))
- fig = plt.figure()
- figA = fig.add_subplot(1,1,1)
- figA.imshow(warped_image)
- plt.show()
- # figA.imshow()
- #
- #
- # plt.imshow(warp(image1, h_to_matrix(h)))
- warped_img_without_holes = remove_holes(image1, warped_image, H, height_shift, width_shift)
- fig = plt.figure()
- figA = fig.add_subplot(1,1,1)
- figA.imshow(warped_img_without_holes)
- plt.show()
- ref_image_height = image2.shape[0]
- ref_image_width = image2.shape[1]
- print('Ref-image size ' + str(image2.shape))
- mergedImage_height = ref_image_height + height_shift
- if warped_image_height > mergedImage_height:
- mergedImage_height = warped_image_height
- mergedImage_width = ref_image_width + width_shift
- if warped_image_width > mergedImage_width:
- mergedImage_width = mergedImage_width
- mergedImage = np.zeros((mergedImage_height , mergedImage_width , 3), dtype=np.uint8);
- # sketch the reference image
- for i in range(0,ref_image_height):
- for j in range(0, ref_image_width):
- # if not( i+ shiftHeight < newHeight and j + shiftWidth < newWidth):
- mergedImage[i + height_shift][j + width_shift] = image2[i][j]
- fig = plt.figure()
- figA = fig.add_subplot(1, 1, 1)
- figA.imshow(mergedImage)
- plt.show()
- print(warped_img_without_holes.shape[0],warped_img_without_holes.shape[1])
- print('shape : \n')
- print(warped_img_without_holes.shape)
- # sketch the destination image
- print(mergedImage_height,mergedImage_width)
- for i in range(0,mergedImage_height-1):
- for j in range(0, mergedImage_width-1):
- if j>=warped_img_without_holes.shape[0] or i >= warped_img_without_holes.shape[0]:
- continue
- 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 ):
- mergedImage[i][j] = warped_img_without_holes[i][j]
- # im = Image.fromarray(mergedImage)
- # im.save("after_add_ref_image.jpg")
- #
- fig = plt.figure()
- figA = fig.add_subplot(1,1,1)
- figA.imshow(mergedImage)
- plt.show()
- '''
- def verifyH(imageA,imageB,h):
- fig = plt.figure()
- figA = fig.add_subplot(1, 2, 1)
- figB = fig.add_subplot(1, 2, 2)
- figA.imshow(image1) # ,origin='lower')
- figB.imshow(image2)
- ptsa = plt.ginput(n=4, timeout=10)
- print((ptsa))
- print(ptsa[0])
- 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)
- for i in range(4):
- tran_point = transform_point(H,np.array([ptsa[i][0],ptsa[i][1],1]))
- figB.scatter(tran_point[0], tran_point[1], c='b', s=40)
- plt.show()
- verifyH(image1,image2,H)
- #10 points
- # [(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)]
- # [(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