ثبت تصویر image registration با استفاده از OpenCV
فرض کنید که یک تصویر هوایی از یک ناحیه شهر دارید و بعد از مدتی یک تصویر دیگر هوایی دریافت می‌کنید که بخش‌هایی از تصویر قبلی در آن موجود است. حالا می‌خواهید دو تصویر را به صورت اتوماتیک طوری روی هم قرار دهید که قسمت‌های متناظر روی هم بی‌افتند. تصویر جدید و قدیم را باهم داشته باشید. یا فرض کنید از شهر یک سری عکس هوایی دریافت کرده اید که این عکس ها با هم اشتراک دارند و شما می‌خواهید آن‌ها را اتوماتیک به هم بچسبایند. یا فرض کنید از آسمان چندین عکس تهیه کرده‌اید و می‌خواهید آن‌ها را کنار هم قرار دهید و یک تصویر واحد بسازید. به این کار ثبت تصویر یا image registration می‌گویند. در این پس می‌خواهیم این کار را با استفاده از پایتون و OpenCV انجام دهیم.
راه‌های مختلفی برای این کار هست که امروز یکی از راه‌های ساده برای این کار را بررسی می‌کنیم و کدش را ارائه می‌دهیم. روش‌های بسیار متنوعی برای این کار وجود دارد و برای کابردهای مختلف روش‌های مختلفی وجود دارد. روش کلی برای انجام این کار وجود ندارد و باید برای کارهای مختلف روش های مختلفی ایجاد کرد تا بهترین نتیجه حاصل شود.
در پست بازشناسی یک شی در تصویر با استفاده از OpenCV در مورد نقاط کلیدی صحبت کردیم. یکی از بخش‌ها اصلی روش ارائه شده در این پست برای تبت تصویراستفاده از نقاط کلیدی و توصیفگرهای ORB است. تابع extract_orb برای تصویر یک سری نقطه کلیدی استخراج می‌کند و برای هر نقطه کلیدی یک توصیفگر ایجاد می‌کند.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 | def extract_orb(img): """ Extract orb feature from image :param img: input image :return: """ # convert to gray if len(img.shape) > 2: img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) else: img_gray = img.copy() # orb keypoint detector and descriptor orb_obj = cv2.ORB_create(nfeatures=1500) # extract orb keypoints and descriptors kp, desc = orb_obj.detectAndCompute(img_gray, None) # keypoints coordinate in 2D image kp = np.array([p.pt for p in kp]).T return kp, desc |
تابع match_orb توصیفگرهای نقاط کلیدی دو تصویر مختلف را می‌گیرد و با استفاده از KNN دو نقطه کلیدی نزدیک به هر نقطه کلیدی را پیدا می‌کند. در صورتی که فاصله توصیفگر نقطه کلیدی تا نقطه کلیدی اول به طرز قابل قبولی کمتر از فاصله نقطه کلیدی تا نقطه کلیدی دوم باشد، نشان می‌دهند که نقطه اول یک نقطه کلیدی متناظر خوب برای نقطه کلیدی است. در غیر این صورت از آن صرف نظر می‌شود و گفته می‌شود تطبیقی برای آن نقطه پیدا نشده‌است. خروجی این تابع اندیس نقاط کلیدی منطبق شده است.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | def match_orb(descriptor_source, descriptor_target): """ math keypoints of two images :param descriptor_source: :param descriptor_target: :return: """ # Matcher bf = cv2.BFMatcher() matches = bf.knnMatch(descriptor_source, descriptor_target, k=2) # keep good matches good_matches = np.array([], dtype=np.int32).reshape((0, 2)) matches_num = len(matches) for i in range(matches_num): if matches[i][0].distance <= 0.8 * matches[i][1].distance: temp = np.array([matches[i][0].queryIdx, matches[i][0].trainIdx]) good_matches = np.vstack((good_matches, temp)) return good_matches |
فرض کنید که نقاط A را داریم و نقاط B نقاط متناظر با نقاط A در یک دستگاه مختصات دیگر هستند که با استفاده از نگاشت X[A]+Y = B به دست می‌آیند. حالا فرض کنید که نقاط A و نقاط B را داده‌اند و می‌خواهیم این نگاشت خطی را پیدا کنیم. برای این کار از تابع estimate_affine_transform استفاده می‌کنیم.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | def estimate_affine_transform(src_points, trg_points): """ Estimate affine transform by getting a set of points and their transformed :param src_points: :param trg_points: :return: """ num_of_points = src_points.shape[1] M = np.zeros((2 * num_of_points, 6)) for i in range(num_of_points): temp = [[src_points[0, i], src_points[1, i], 0, 0, 1, 0], [0, 0, src_points[0, i], src_points[1, i], 0, 1]] M[2 * i: 2 * i + 2, :] = np.array(temp) b = trg_points.T.reshape((2 * num_of_points, 1)) theta = np.linalg.lstsq(M, b)[0] X_trans = theta[:4].reshape((2, 2)) Y_trans = theta[4:] return X_trans, Y_trans |
فرض کنید که نقاط A را در یک دستگاه مختصات داریم و می‌دانیم که اگر آن‌ها را با نگاشت خطی X[A] + Y به یک دستگاه مختصات دیگر ببریم، نقطه متناظر نقاط A، نقاط B می‌شود. حالا فرض کنید که A و B را داریم ولی نگاشت خطی X[ ]+Y را نداریم. به ما یک نگاشت خطی P[ ] + Q می‌دهند و می‌خواهیم بدانیم که چقدر این نگاشت خطی دقیق است و عمل نگاشت را درست انجام می‌دهد. برای این کار نقاط A را با استفاده از PA+Q به مختصات جدید می‌بریم و مجذور مربعات فاصله نقاط را با نقاط متناظرشان که B باشد می‌سنجیم. هر چه این مقدار کمتر باشد نشان می‌دهد نگاشت خطی داده شده بهتر است. این کاری است که تابع mse_error انجام می‌دهد.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | def mse_error(X_trans, Y_trans, src_points, trg_points): """ :param X_trans: transform parameter X[]+Y :param Y_trans: transform parameter X[]+Y :param src_points: source points :param trg_points: traget points :return: """ # transform source points by affine transfor X[]+Y src_map_on_trg = np.dot(X_trans, src_points) + Y_trans # calculate MSE error diff_square = np.power(src_map_on_trg - trg_points, 2) mse = np.sqrt(np.sum(diff_square, axis=0)) # return mse error return mse |
همه نقاط متناظری که با استفاده از تطبیق نقاط کلیدی با استفاده از توصیفگر‌هایشان به دست آمد، تقاط متناظر خوبی نیستند و اشتباهاتی در آن‌ها موجود است. تابع RANSAC هر بار تعداد از نقاط را انتخاب می‌کند. سپس با استفاده از آن نقاط و تابع estimate_affine_transform یک نگاشت خطی پیدا می‌کند و سپس با استفاده mse_error دقت نگاشت به دست آمده برای سایر نقاط را می‌سنجد. اینگونه بهترین نگاشت را پیدا می‌کند و تناظر‌های خوب و بد را از هم تشخیص می‌دهد.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 | def ransac(points_src, points_target): """ :param points_src: :param points_target: :return: """ # in each iteration select k point randomly K = 3 # if mse error in less than threshold consider it mse_threshold = 1.0 # maximum number of iterations number_of_iteration = 3000 inliers_num_final = 0 X_trans_final = None Y_trans_final = None inliers = None # for number of iteration repeat for i in range(number_of_iteration): # select K points randomly k_selected = np.random.randint(0, points_src.shape[1], (K, 1)) # find affine transform by K selected points X_trans, Y_trans = estimate_affine_transform(src_points=points_src[:, k_selected], trg_points=points_target[:, k_selected]) # calculate mse error of affine transform mse = mse_error(X_trans=X_trans, Y_trans=Y_trans, src_points=points_src, trg_points=points_target) if not (mse is None): inliers_tmp = np.where(mse < mse_threshold) inliers_num_tmp = len(inliers_tmp[0]) # keep best affine transform until now if inliers_num_tmp > inliers_num_final: inliers_num_final = inliers_num_tmp inliers = inliers_tmp X_trans_final = X_trans Y_trans_final = Y_trans else: pass # return affine transform X[]+Y and inliers return X_trans_final, Y_trans_final, inliers |
تابع affine_matrix نقاط کلیدی متناظر که توسط match_orb پیدا شده‌اند را می‌گیرد. سپس با استفاده از ransac تناظر‌های خوب از میان آن‌ها را پیدا می‌کند و تناظر‌های بد را کنار می‌گذارد. سپس با استفاده از تناظر‌های خوب پیدا شده توسط Ransac و تابع estimate_affine_transform، ماتریس نگاشت خطی M را پیدا می‌کند.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 | def affine_matrix(points_in_src_img, points_in_trg_img, good_matches): """ Calculate affine transform :param points_in_src_img: :param points_in_trg_img: :param good_matches: :return: """ # keypoints in image 1 that are in matched points points_in_src_img = points_in_src_img[:, good_matches[:, 0]] # keypoints in image 2 that are in matched points points_in_trg_img = points_in_trg_img[:, good_matches[:, 1]] # calculate inliers with RANSAC algorithm _, _, inliers = ransac(points_in_src_img, points_in_trg_img) # keypoints in image 1 that are in inliers points points_in_src_img = points_in_src_img[:, inliers[0]] # keypoints in image 2 that are in inliers points points_in_trg_img = points_in_trg_img[:, inliers[0]] # estimate affine transform X[]+Y by inliers points in first and second image X_trans, Y_trans = estimate_affine_transform(points_in_src_img, points_in_trg_img) M = np.hstack((X_trans, Y_trans)) return M |
تابع two_image_registraion_without_padding، دو تصویر را می‌گیرد و با استفاده از تابع‌هایی که گفته شده کار ثبت تصویر یا image registration را انجام می‌دهد. فقط مشکلی که این تابع دارد این است که بعد از انتقال تصویر اول به دستگاه مختصات تصویر دوم، ممکن است بخشی از تصویر اول بیرون بزند. این مشکل در تابع two_images_registraion_padding حل شده است. تابع multiple_image_registration مانند تابع two_images_registraion_padding است با این تفاوت که بیش از دو تصویر می‌گیرد.
کد کلی برنامه درانتهای مطلب آورده شده‌است. در ادامه به بعضی از نمونه‌های خروجی این برنامه می‌پردازیم:
نقشه تهران:
تصاویر ورودی:
خروجی الگوریتم:
تصویر هوایی تهران:
خروجی برنامه:
مقبره حکیم فردوسی:
خروجی برنامه:
راه شیری:
خروجی برنامه:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 | import numpy as np import cv2 def extract_orb(img): """ Extract orb feature from image :param img: input image :return: """ # convert to gray if len(img.shape) > 2: img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) else: img_gray = img.copy() # orb keypoint detector and descriptor orb_obj = cv2.ORB_create(nfeatures=1500) # extract orb keypoints and descriptors kp, desc = orb_obj.detectAndCompute(img_gray, None) # keypoints coordinate in 2D image kp = np.array([p.pt for p in kp]).T return kp, desc def match_orb(descriptor_source, descriptor_target): """ math keypoints of two images :param descriptor_source: :param descriptor_target: :return: """ # Matcher bf = cv2.BFMatcher() matches = bf.knnMatch(descriptor_source, descriptor_target, k=2) # keep good matches good_matches = np.array([], dtype=np.int32).reshape((0, 2)) matches_num = len(matches) for i in range(matches_num): if matches[i][0].distance <= 0.8 * matches[i][1].distance: temp = np.array([matches[i][0].queryIdx, matches[i][0].trainIdx]) good_matches = np.vstack((good_matches, temp)) return good_matches def estimate_affine_transform(src_points, trg_points): """ Estimate affine transform by getting a set of points and their transformed :param src_points: :param trg_points: :return: """ num_of_points = src_points.shape[1] M = np.zeros((2 * num_of_points, 6)) for i in range(num_of_points): temp = [[src_points[0, i], src_points[1, i], 0, 0, 1, 0], [0, 0, src_points[0, i], src_points[1, i], 0, 1]] M[2 * i: 2 * i + 2, :] = np.array(temp) b = trg_points.T.reshape((2 * num_of_points, 1)) theta = np.linalg.lstsq(M, b)[0] X_trans = theta[:4].reshape((2, 2)) Y_trans = theta[4:] return X_trans, Y_trans def mse_error(X_trans, Y_trans, src_points, trg_points): """ :param X_trans: transform parameter X[]+Y :param Y_trans: transform parameter X[]+Y :param src_points: source points :param trg_points: traget points :return: """ # transform source points by affine transfor X[]+Y src_map_on_trg = np.dot(X_trans, src_points) + Y_trans # calculate MSE error diff_square = np.power(src_map_on_trg - trg_points, 2) mse = np.sqrt(np.sum(diff_square, axis=0)) # return mse error return mse def ransac(points_src, points_target): """ :param points_src: :param points_target: :return: """ # in each iteration select k point randomly K = 3 # if mse error in less than threshold consider it mse_threshold = 1.0 # maximum number of iterations number_of_iteration = 3000 inliers_num_final = 0 X_trans_final = None Y_trans_final = None inliers = None # for number of iteration repeat for i in range(number_of_iteration): # select K points randomly k_selected = np.random.randint(0, points_src.shape[1], (K, 1)) # find affine transform by K selected points X_trans, Y_trans = estimate_affine_transform(src_points=points_src[:, k_selected], trg_points=points_target[:, k_selected]) # calculate mse error of affine transform mse = mse_error(X_trans=X_trans, Y_trans=Y_trans, src_points=points_src, trg_points=points_target) if not (mse is None): inliers_tmp = np.where(mse < mse_threshold) inliers_num_tmp = len(inliers_tmp[0]) # keep best affine transform until now if inliers_num_tmp > inliers_num_final: inliers_num_final = inliers_num_tmp inliers = inliers_tmp X_trans_final = X_trans Y_trans_final = Y_trans else: pass # return affine transform X[]+Y and inliers return X_trans_final, Y_trans_final, inliers def affine_matrix(points_in_src_img, points_in_trg_img, good_matches): """ Calculate affine transform :param points_in_src_img: :param points_in_trg_img: :param good_matches: :return: """ # keypoints in image 1 that are in matched points points_in_src_img = points_in_src_img[:, good_matches[:, 0]] # keypoints in image 2 that are in matched points points_in_trg_img = points_in_trg_img[:, good_matches[:, 1]] # calculate inliers with RANSAC algorithm _, _, inliers = ransac(points_in_src_img, points_in_trg_img) # keypoints in image 1 that are in inliers points points_in_src_img = points_in_src_img[:, inliers[0]] # keypoints in image 2 that are in inliers points points_in_trg_img = points_in_trg_img[:, inliers[0]] # estimate affine transform X[]+Y by inliers points in first and second image X_trans, Y_trans = estimate_affine_transform(points_in_src_img, points_in_trg_img) M = np.hstack((X_trans, Y_trans)) return M def two_image_registraion_without_padding(img_source, img_target): """ Gets two image an do image registration :param img_source: :param img_target: :return: """ # extract orb features for both images keypoint_source, descriptor_source = extract_orb(img_source) keypoint_target, descriptor_target = extract_orb(img_target) # math keypoints of two image matched_points = match_orb(descriptor_source, descriptor_target) # calculate affine matrix M = affine_matrix(keypoint_source, keypoint_target, matched_points) rows, cols, _ = img_target.shape # transform source image by H transformed_src_img = cv2.warpAffine(src=img_source, M=M, dsize=(cols, rows)) # find where will be corners of source image new_corners = np.dot(a=M, b=np.array([[0, img_source.shape[1]-1, 0, img_source.shape[1]-1], [0, 0, img_source.shape[0]-1, img_source.shape[0]-1], [1, 1, 1, 1]])) # merge two images merged_image = cv2.addWeighted(src1=transformed_src_img, alpha=0.5, src2=img_target, beta=0.5, gamma=0) return merged_image, new_corners def two_images_registraion_padding(img_source, img_target): """ Does image registration with padding. no part of imaged will be lost :param img_source: :param img_target: :return: """ # do image registration without padding merged_image_without_padding, new_corners = two_image_registraion_without_padding(img_source=img_source, img_target=img_target) # find who two padd x_min = int(np.floor(new_corners[0, :].min())) x_max = int(np.ceil(new_corners[0, :].max())) y_min = int(np.floor(new_corners[1, :].min())) y_max = int(np.ceil(new_corners[1, :].max())) left, right, top, bottom = 0, 0, 0, 0 if x_min < 0: left = -x_min if x_max > img_target.shape[1]: right = x_max-img_target.shape[1] if y_min < 0: top = -y_min if y_max > img_target.shape[0]: bottom = y_max - img_target.shape[0] # pad target image padded_img_target = cv2.copyMakeBorder(src=img_target, top=top, bottom=bottom, left=left, right=right, borderType=cv2.BORDER_CONSTANT) # do image registration merged_image_with_padding, _ = two_image_registraion_without_padding(img_source=img_source, img_target=padded_img_target) return merged_image_with_padding def multiple_image_registration(list_imgs): """ gets two or more than two images and do image registration by all of them :param list_imgs: :return: """ if len(list_imgs) < 2: raise ValueError('Need at least 2 images!') for i in range(1, len(list_imgs)): if i == 1: merged_image = two_images_registraion_padding(img_source=list_imgs[0], img_target=list_imgs[1]) else: merged_image = two_images_registraion_padding(img_source=merged_image, img_target=list_imgs[i]) return merged_image if __name__ == "__main__": np.random.seed(1) # ================================================ # test image registration with two image # read images img_source = cv2.imread("tehran-1.png") img_target = cv2.imread("tehran-2.png") # do image registration merged_image = two_images_registraion_padding(img_source=img_source, img_target=img_target) # save image cv2.imwrite(filename='result_of_registration_two_images_tehran_map.png', img=merged_image) # cv2.imshow('img_source', img_source) # cv2.imshow('img_target', img_target) cv2.imshow('result_of_registration_two_images_tehran_map', merged_image) # ================================================ # test image registration with two image # read images img_source = cv2.imread("ferdosi-1.png") img_target = cv2.imread("ferdosi-2.png") # do image registration merged_image = two_images_registraion_padding(img_source=img_source, img_target=img_target) # save image cv2.imwrite(filename='result_of_registration_two_images_hakim_ferdosi.png', img=merged_image) # cv2.imshow('img_source', img_source) # cv2.imshow('img_target', img_target) cv2.imshow('result_of_registration_two_images_hakim_ferdosi', merged_image) # ================================================ # test image registration with several images img1 = cv2.imread(filename='teh1.png') img2 = cv2.imread(filename='teh2.png') img3 = cv2.imread(filename='teh3.png') # do image registration with several image merged_image = multiple_image_registration(list_imgs=[img1, img2, img3]) cv2.imwrite(filename='result_of_registration_several_images_tehran.png', img=merged_image) cv2.imshow('result_of_registration_several_images_tehran', merged_image) # ================================================ # test image registration with several images img1 = cv2.imread(filename='milkyway1.png') img2 = cv2.imread(filename='milkyway2.png') img3 = cv2.imread(filename='milkyway3.png') # do image registration with several image merged_image = multiple_image_registration(list_imgs=[img1, img2, img3]) cv2.imwrite(filename='result_of_registration_several_images_milkyway.png', img=merged_image) cv2.imshow('result_of_registration_several_images_milkyway', merged_image) cv2.waitKey(0) cv2.destroyAllWindows() |
سلام وقت بخیر
مطلب خیلی جالبی بود یادمه ترم ۲ ارشد، استاد پردازش تصویرمون یه فوق تمرین داده بود ۳ تا تصویر ماهواره‌ای ای بهمون داد که این ۳ تصویر بهم بچسبونید هرکاری کردم نشد!
ممنونم استفاده بردیم