Lucas-Kanade 图像对齐算法实现不收敛？

Question

我最近一直在尝试实现用于图像对齐的 Lucas-Kanade 算法，本文对此进行了详细介绍：https://www.ri.cmu.edu/pub_files/pub3/baker_simon_2004_1/baker_simon_2004_1.pdf

我已经设法实现了我链接的论文第 4 页中详述的算法，但损失似乎没有收敛。我一直在查看我的代码和数学，但似乎无法弄清楚我可能哪里出错了。

到目前为止，我尝试的是实现整个算法，并重新计算计算扭曲的雅可比行列式，以及对我的代码进行一般检查。

下面是我的代码，以及 Pastebin 上更易读的版本：https://pastebin.com/j28mUV65

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)

它应该只需要一个参考图像来进行测试。

预期的结果是算法收敛到与我在代码中计算的基本事实单应性相匹配的单应性，但损失似乎反而爆炸了，我最终得到一个完全不正确的单应性。

Answer 1

这应该是一条评论，因为我不确定是不是您的问题的全部原因

但它可能是其中的一部分

求解线性方程组不要计算逆

hessian_inverse = np.linalg.inv(hessian)

然后乘以它

final_param_updates = np.dot(hessian_inverse, sd_param_updates)

这既浪费又会导致比求解线性方程组通常具有更多的数值不稳定性。

改为使用方法solve。

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