LucasKanade图像对齐算法实现不收敛？

import cv2 import numpy as np import matplotlib.pyplot as plt def calculate_steepest_descent(grad_x_warped, grad_y_warped, h): rows, columns = grad_x_warped.shape steepest_descent = np.zeros((rows, columns, 8)) warp_jacobian = np.zeros((2, 8)) # 2 x 8 because it's a homography, would be 2 x 6 if it was affine current_gradient = np.zeros((1, 2)) # Convert homography matrix into parameter array for better readability with the math functions later p = h.flatten() for y in range(rows): for x in range(columns): # Calculate Jacobian of the warp at each pixel, which contains the partial derivatives of the # warp parameters with respect to x and y coordinates, evaluated at the current value # of parameters common_denominator = (p[6]*x + p[7]*y + 1) warp_jacobian[0, 0] = (x) / common_denominator warp_jacobian[0, 1] = (y) / common_denominator warp_jacobian[0, 2] = (1) / common_denominator warp_jacobian[0, 3] = 0 warp_jacobian[0, 4] = 0 warp_jacobian[0, 5] = 0 warp_jacobian[0, 6] = (-(p[0]*(x**2) + p[1]*x*y + p[2]*x)) / (common_denominator ** 2) warp_jacobian[0, 7] = (-(p[1]*(y**2) + p[0]*x*y + p[2]*y)) / (common_denominator ** 2) warp_jacobian[1, 0] = 0 warp_jacobian[1, 1] = 0 warp_jacobian[1, 2] = 0 warp_jacobian[1, 3] = (x) / common_denominator warp_jacobian[1, 4] = (y) / common_denominator warp_jacobian[1, 5] = (1) / common_denominator warp_jacobian[1, 6] = (-(p[3]*(x**2) + p[4]*x*y + p[5]*x)) / (common_denominator ** 2) warp_jacobian[1, 7] = (-(p[4]*(y**2) + p[3]*x*y + p[5]*y)) / (common_denominator ** 2) # Get the x and y gradient intensity values corresponding to the current pixel location current_gradient[0, 0] = grad_x_warped[y, x] current_gradient[0, 1] = grad_y_warped[y, x] # Calculate full Jacobian (aka steepest descent image) at current pixel value steepest_descent[y, x, :] = np.dot(current_gradient, warp_jacobian) return steepest_descent def calculate_hessian(steepest_descent): rows, columns, channels = steepest_descent.shape hessian = np.zeros((channels, channels)) for y in range(rows): for x in range(columns): steepest_descent_single = steepest_descent[y, x, :][np.newaxis, :] steepest_descent_single_transpose = np.transpose(steepest_descent_single) hessian_current = np.dot(steepest_descent_single_transpose, steepest_descent_single) hessian += hessian_current return hessian def calculate_sd_param_updates(steepest_descent, img_error): rows, columns, channels = steepest_descent.shape sd_param_updates = np.zeros((8, 1)) for y in range(rows): for x in range(columns): steepest_descent_single = steepest_descent[y, x, :][np.newaxis, :] steepest_descent_single_transpose = np.transpose(steepest_descent_single) img_error_single = img_error[y, x] sd_param_updates += np.dot(steepest_descent_single_transpose, img_error_single) return sd_param_updates def calculate_final_param_updates(sd_param_updates, hessian): hessian_inverse = np.linalg.inv(hessian) final_param_updates = np.dot(hessian_inverse, sd_param_updates) return final_param_updates if __name__ == "__main__": # Load image reference = cv2.imread('test.png') reference = cv2.cvtColor(reference, cv2.COLOR_BGR2GRAY) # Generate template as small block from within reference image using homography # 'h' is the ground truth homography for warping reference image onto template image template_size = (100, 100) h = np.float32([[1, 0, -100],[0, 1, -100],[0, 0, 1]]) h_ground_truth = h.copy() template = cv2.warpPerspective(reference, h, template_size) # Convert template corner points to reference image coordinate plane template_corners = np.array([[0, 0],[0, 100],[100, 100],[100, 0]]) h_inverse = np.linalg.inv(h) reference_corners = cv2.perspectiveTransform(np.array([template_corners], dtype='float32'), h_inverse) # Small perturbation to ground truth homography h_mod = np.random.uniform(low=-1.0, high=1.0, size=(h.shape)) h_mod = np.array([[1, 1, 1],[1, 1, 1],[1, 1, 1]]) h_mod[0, 0] = h_mod[0, 0] * 0 h_mod[0, 1] = -h_mod[0, 1] * 0 h_mod[0, 2] = h_mod[0, 2] * 10 h_mod[1, 0] = h_mod[1, 0] * 0 h_mod[1, 1] = h_mod[1, 1] * 0 h_mod[1, 2] = h_mod[1, 2] * 10 h_mod[2, 0] = h_mod[2, 0] * 0 h_mod[2, 1] = h_mod[2, 1] * 0 h_mod[2, 2] = h_mod[2, 1] * 0 h = h + h_mod # Warp reference image to template image based on initial perturbed homography reference_transformed = cv2.warpPerspective(reference, h, template_size) # ############################## # Lucas-Kanade algorithm below # This is supposed to calculate the homography that undoes the small perturbation # and returns a homography as close as possible to the ground truth homography # ############################## # Precompute image gradients grad_x = cv2.Sobel(reference,cv2.CV_64F,1,0,ksize=1) grad_y = cv2.Sobel(reference,cv2.CV_64F,0,1,ksize=1) # Loop algorithm for given # of steps for i in range(1000): # Step 1 # Warp reference image onto coordinate frame of template reference_transformed = cv2.warpPerspective(reference, h, template_size) # Step 2 # Compute error image img_error = template - reference_transformed # fig_overlay = plt.figure() # ax1 = fig_overlay.add_subplot(1,3,1) # plt.imshow(img_warped) # ax2 = fig_overlay.add_subplot(1,3,2) # plt.imshow(template) # ax3 = fig_overlay.add_subplot(1,3,3) # plt.imshow(img_error) # plt.show() # Step 3 # Warp the gradients grad_x_warped = cv2.warpPerspective(grad_x, h, template_size) grad_y_warped = cv2.warpPerspective(grad_y, h, template_size) # Step 4 & 5 # Use Jacobian of warp to calculate steepest descent images steepest_descent = calculate_steepest_descent(grad_x_warped, grad_y_warped, h) # fig_overlay = plt.figure() # ax1 = fig_overlay.add_subplot(1,8,1) # plt.imshow(steepest_descent[:, :, 0]) # ax2 = fig_overlay.add_subplot(1,8,2) # plt.imshow(steepest_descent[:, :, 1]) # ax3 = fig_overlay.add_subplot(1,8,3) # plt.imshow(steepest_descent[:, :, 2]) # ax4 = fig_overlay.add_subplot(1,8,4) # plt.imshow(steepest_descent[:, :, 3]) # ax5 = fig_overlay.add_subplot(1,8,5) # plt.imshow(steepest_descent[:, :, 4]) # ax6 = fig_overlay.add_subplot(1,8,6) # plt.imshow(steepest_descent[:, :, 5]) # ax7 = fig_overlay.add_subplot(1,8,7) # plt.imshow(steepest_descent[:, :, 6]) # ax8 = fig_overlay.add_subplot(1,8,8) # plt.imshow(steepest_descent[:, :, 7]) # plt.show() # Step 6 # Compute Hessian matrix hessian = calculate_hessian(steepest_descent) # Step 7 # Compute steepest descent parameter updates by # dot producting error image with steepest descent images sd_param_updates = calculate_sd_param_updates(steepest_descent, img_error) # Step 8 # Compute final parameter updates final_param_updates = calculate_final_param_updates(sd_param_updates, hessian) # Step 9 # Update the parameters h = h.reshape(-1,1) h[:-1] += final_param_updates h = h.reshape(3,3) # Step 10 # Calculate norm of parameter updates final_param_update_norm = np.linalg.norm(final_param_updates) print("Final Param Norm: {}".format(final_param_update_norm)) reference_transformed = cv2.warpPerspective(reference, h, template_size) cv2.imwrite('warps/warp_{}.png'.format(i), reference_transformed) # Warp source image to destination based on homography reference_transformed = cv2.warpPerspective(reference, h, template_size) cv2.imwrite('final_warp.png', reference_transformed)

1条回答

网友

1楼 · 发布于 2024-04-27 17:40:12

你的评论不应该是完全的原因

但这可能是其中的一部分

解线性方程组时，不要求逆

hessian_inverse = np.linalg.inv(hessian)

再乘以它

^{pr2}$

这不仅是浪费，而且可能导致比求解线性方程组更大的数值不稳定性。在

而是使用方法solve。在

计算逆矩阵将重复对单位矩阵的每一列执行solve所需的一些操作。这些手术都不需要。在

相关问题更多 >

编程相关推荐

热门问题

热门文章