高阶畸变校正

本文最后更新于：2023年4月15日 晚上

我们之前介绍了镜头畸变，记录了镜头畸变校正的常用方法，当对畸变精度要求更高时，传统的模型已经到达极限，本文提出一种新的畸变校正方案，以进一步提高畸变校正精度。

背景

对于超高精度的畸变校正要求，需要：

  1. 精确的标定板
  2. 精准的角点提取算法
  3. 合理、正确的畸变模型
  4. 准确地求解模型参数
  5. 减少畸变校正过程中的截断信息丢失
  6. 精确严苛的评估方法

等条件，任何一个步骤出问题都将使结果万劫不复。

相机为远心镜头，经过光学调整，确认存在畸变，但是没有除了流行的畸变模型外的模型可用。

这里记录我对畸变校正问题精确求解的探索过程，以及相关实验结果。

标定板

采用等间距的 bump 构成的板子，通过拟合每一个圆，提取圆心位置来定位角点。

角点提取

核心采用了 Halcon 的 fit_circle_contour_xld 算子或 fit_ellipse_contour_xld 算子

1
2
3
4
5
6
edges_sub_pix (ImageReduced1, Edges, 'canny', 1, 10, 40)

sort_contours_xld (Edges1, SortedContours, 'character', 'true', 'row')
fit_circle_contour_xld (SortedContours, 'algebraic', -1, 0, 0, 3, 2, CenterRowTuple, CenterColTuple, Radius, StartPhi, EndPhi, PointOrder)

fit_ellipse_contour_xld (SortedContours, 'fitzgibbon', -1, 0, 0, 200, 3, 2, Row1, Column1, Phi, Radius1, Radius2, StartPhi1, EndPhi1, PointOrder1)

使用拟合后的圆心作为角点坐标。

评估方法

我们要在尽可能不引入不确定的误差项的前提下，精确、公平地评估畸变校正的效果。

精度高的时候，所有的事情都是不确定的，我们能抓住的只有标定板角点之间是边长相等的正方形这一个信息。

据此我设计了评估畸变校正效果的方法：

  1. 获取现有畸变图像的角点 $m^*$
  2. 通过当前畸变校正方法 $D^{-1}$ 将这些点校正到当前认为的理想位置 $m’=D{-1}(m*)$
  3. 度量 $m’$ 中的相邻点距离，计算均值 $l_m$，选定一个固定的常数长度 $l_0$，计算倍率变化 $r = \frac{l_0}{l_m}$
  4. 修改 $m’$ 的点为 $m’‘=rm’$，这样 $m’'$ 的平均间隔为 $l_0$ 了
  5. 度量 $m’'$ 所有边长、对角线的标准差 $std_l,std_d$ ，标准差越小说明精度越高

畸变模型

传统模型

传统畸变校正需要考虑 相机成像的几何原理，计算相机的内参、外参以建立坐标系转换模型，一共需要求解 11 个参数，才能获取真实世界坐标到像素坐标的理想情况下（无畸变）的转换关系；
$$
s {m}=\mathbf{A}[\mathbf{R} \mid \mathbf{t}] {\mathbf{M}}
$$
如果我们虚拟一批理想的阵列作为初始位置 $ M_w$，经过上述模型可以得到理论上阵列角点在图像上的像素坐标位置 ${m}$；

我们从真实图像上获取的仅是畸变后的角点位置 $m^*$，畸变 $m$ 点的过程我们用 $D(m) $ 表示，有 $m^*=D(m) $；

就算传统模型是真正正确的，需要解最少 5 个参数；

$$ \begin{aligned} x_{\text {distorted }} & =x\left(1+k_{1} r^{2}+k_{2} r^{4}+k_{3} r^{6}\right)+2 p_{1} x y+p_{2}\left(r^{2}+2 x^{2}\right) \\ y_{\text {distorted }} & =y\left(1+k_{1} r^{2}+k_{2} r^{4}+k_{3} r{6}\right)+p_{1}\left(r{2}+2 y^{2}\right)+2 p_{2} x y\end{aligned} $$ 那么事实上的模型为： $$ m^*=D(\mathbf{A}[\mathbf{R} \mid \mathbf{t}] {\mathbf{M}}) $$

这里我们有的数据是 $m ^ *$ 和 $ M_w$，讲道理可以直接 16 个参数一起优化，但是现在常用的方法都是先解出内外参，再求畸变参数，这必然是一个来回震荡迭代的过程，存在诸多会产生误差的过程，包括：

  1. 成像模型本身的误差
  2. 畸变校正模型的误差
  3. 交错迭代程度的误差

因此我尝试调整模型，解决此类误差引入的问题。

正向模型

采用最优化的思路，构建一个模型，根据当前数据产生一个 Loss，通过优化这个 Loss 确定模型参数。

正向模型的含义是模型自变量为理想坐标，输出的结果为畸变坐标，我们常见的畸变校正模型大多都是正向模型，模型的正反基本上是定死的，很难反过来；

我首先认为理想模型到像素模型的本质是一个仿射变换，那么我可以用一个仿射变换直接等效相机内外参；

得到像素理想模型后直接添加畸变模型，获取到的点应该和我们提取的角点 $ m^*$ 接近，那么 Loss 就设置为二者之间的距离和。

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
import mtutils as mt
import numpy as np
from scipy.optimize import minimize
import cv2


def undistort(corner_info, img):
    base_ideal_point_list = list()
    real_point_list = list()
    
    line_list = corner_info['line_list']
    
    for Y in range(corner_info['H_num']):
        for X in range(corner_info['W_num']):
            real_point_list.append(line_list[Y][X])
            base_ideal_point_list.append([X, Y])

    base_ideal_point_array = np.ones([3, len(base_ideal_point_list)])
    base_ideal_point_array[:2, :] = np.array(base_ideal_point_list).T

    real_point_array = np.ones([3, len(base_ideal_point_list)])
    real_point_array[:2, :] = np.array(real_point_list).T

    def make_ideal_matrix(aff_para, base_ideal_point_array):
        dx, dy, theta, sx = aff_para
        movd_matrix = np.matrix([[1, 0, dx], [0, 1, dy], [0, 0, 1]])
        rotate_matrix = np.matrix([[np.cos(theta), - np.sin(theta), 0], [np.sin(theta), np.cos(theta), 0], [0, 0, 1]])
        scale_matrix = np.matrix([[sx, 0, 0], [0, sx, 0], [0, 0, 1]])

        affined_points = movd_matrix * rotate_matrix * scale_matrix * base_ideal_point_array
        return affined_points

    def make_ideal_matrix2(aff_para, base_ideal_point_array):
        a, b, c, d, e, f = aff_para
        affine_matrix = np.matrix([[a, b, c], [d, e, f], [0, 0, 1]])

        affined_points = affine_matrix * base_ideal_point_array
        return affined_points

    def affine_init_para_loss(aff_para):
        affined_points = make_ideal_matrix(aff_para, base_ideal_point_array)
        diff_matrix = affined_points - real_point_array
        loss = (np.array(diff_matrix) ** 2).sum()
        return loss

    def affine_init_para_loss2(aff_para):
        affined_points = make_ideal_matrix2(aff_para, base_ideal_point_array)
        diff_matrix = affined_points - real_point_array
        loss = (np.array(diff_matrix) ** 2).sum()
        return loss

    x0_affine = (140, 52, 0, 91)
    x0_affine2 = (91, 0, 140, 0, 91, 50)
    
    res_affine = minimize(affine_init_para_loss, x0_affine, method='Powell')
    res_affine2 = minimize(affine_init_para_loss2, x0_affine2, method='Powell')

    affined_points = make_ideal_matrix(res_affine.x, base_ideal_point_array)
    affined_points = np.array(affined_points).T[:, :2]


    down_scale = 1
    h, w = img.shape[:2]

    def distort_model(distort_para, affined_points):
        cx, cy, k1x, k2x, p1x, p2x, k0x, k1y, k2y, p1y, p2y, k0y = distort_para
        r_array = ((affined_points / down_scale - [cx, cy]) ** 2).sum(axis = 1) ** 0.5
        r_2 = r_array ** 2
        r_4 = r_array ** 1

        X = affined_points[:, 0] / down_scale
        Y = affined_points[:, 1] / down_scale
        X_2 = X ** 2
        Y_2 = Y ** 2
        XY = X * Y

        X_distorted = X * k0x *(1 + k1x * r_2 + k2x * r_4) + 2 * p1x * XY + p2x * (r_2 + 2 * X_2)
        Y_distorted = Y * k0y * (1 + k1y * r_2 + k2y * r_4) + p1y * (r_2 + 2 * Y_2) + 2 * p2y * XY

        distorted_points = np.vstack([X_distorted, Y_distorted]).T
        return distorted_points

    real_points = real_point_array.T[:, :2] / down_scale

    def distort_loss(total_para):
        a, b, c, d, e, f, cx, cy, k1x, k2x, p1x, p2x, k0x, k1y, k2y, p1y, p2y, k0y = total_para
        aff_para = a, b, c, d, e, f
        distort_para = cx, cy, k1x, k2x, p1x, p2x, k0x, k1y, k2y, p1y, p2y, k0y
        affined_points = make_ideal_matrix2(aff_para, base_ideal_point_array)
        affined_points = np.array(affined_points).T[:, :2]
        distort_points = distort_model(distort_para, affined_points)
        diff_matrix = (distort_points - real_points) * down_scale
        loss = (np.abs(np.array(diff_matrix))).sum() / len(diff_matrix)
        return loss

    x0_distort = 91, 0, 140, 0, 91, 50, 2000, 1500, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1
    res_distort = minimize(distort_loss, x0_distort, method='Powell')
    print(res_distort.fun)
    print()
    print(res_distort.x)

这个模型可以获得正向计算的过程，如果需要生成畸变映射矩阵，仅需将所有图像上的整数像素点作为自变量输入模型，对应输出的坐标就是畸变图像的位置，那么将这些畸变坐标保存下来就是畸变矩阵了；

不幸的是这个方法的效果不如传统模型，我怀疑是同时拟合理想模型和畸变校正参数太过于困难，因此修改了模型定义，提出了更加直接的反向模型。

反向模型

反向模型的含义是畸变校正模型输入的是畸变坐标，输出的是理想坐标，这种模型很少见，没有理论依据，索性我就放弃理论依据，自己胡搞一个模型。

$$ \begin{array}{c}r=\sqrt{\left(x-c_{x}\right)^{2}+\left(y-c_{y}\right)^{2}} \\ x^{\prime}=x\left(k_{1 x} r+k_{2 x} r^{2}\right)+p_{1 x} x y+p_{2 x} x^{2}+p_{3 x} y^{2}+p_{4 x} y+p_{5 x} x \\ y^{\prime}=y\left(k_{1 y} r+k_{2 y} r^{2}\right)+p_{1 y} x y+p_{2 y} x^{2}+p_{3 y} y^{2}+p_{4 y} x+p_{5 y} y\end{array} $$

其中自变量 $x,y$ 为畸变坐标，$x’,y’$ 为变换后的理想坐标，模型参数一共 16 个。

这样对于输入的畸变坐标点 $m^*$，就可以通过该模型获取对应的理想坐标点了，然后理想坐标点啊，就可以通过我们的评估系统得到误差，计算出一个 Loss，然后最优化这个 Loss，求解模型参数，这是端到端的过程，没有引入中间模型，消除了中间模型的拟合误差；

而且模型专注于反畸变，参数更多表达能力更强，可以拟合更多我们不知道的畸变因素。

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
import mtutils as mt
import numpy as np
from scipy.optimize import minimize
import cv2


def undistort2(corner_info, img):

    X_list = corner_info['cols']
    Y_list = corner_info['rows']

    h, w = img.shape[:2]

    corner_list_temp = list()
    for index in range(len(X_list)):
        corner_list_temp.append([X_list[index], Y_list[index]])

    analysis_info, data_line_list = eva_res_info(np.array(corner_list_temp), corner_info['H-Num'], corner_info['W-Num'])
    print(analysis_info)

    data_array = np.array(data_line_list)

    X = data_array[:,:, 0]
    Y = data_array[:,:, 1]

    def make_new_points(para, X=X, Y=Y):
        Temp_X = X - w / 2
        Temp_Y = Y - h / 2

        cx, cy, k1x, k2x, k1y, k2y, p1x, p2x, p1y, p2y, p3x, p3y, p4x, p4y, p5x, p5y= para
        r = ((Temp_X - cx) ** 2 + (Temp_Y - cy) ** 2) ** 0.5
        r_2 = r ** 2
        r_4 = r ** 4
        xy = Temp_X * Temp_Y
        x_2 = Temp_X ** 2
        y_2 = Temp_Y ** 2

        X_r = Temp_X * (k1x * r + k2x * r_2) + p1x * xy + p2x * x_2 + p3x * y_2 + p4x * Temp_Y + p5x * Temp_X
        Y_r = Temp_Y * (k1y * r + k2y * r_2) + p1y * xy + p2y * y_2 + p3y * x_2 + p4y * Temp_X + p5y * Temp_Y

        return X_r + w / 2, Y_r + h / 2

    def make_loss(para):
        X_r, Y_r = make_new_points(para)

        loss1 = np.abs(((X_r[:, 1:] - X_r[:,:-1]) ** 2 + (Y_r[:, 1:] - Y_r[:,:-1]) ** 2) ** 0.5 - StandardLength).sum()
        loss2 = np.abs(((X_r[1:, :] - X_r[:-1,:]) ** 2 + (Y_r[1:, :] - Y_r[:-1,:]) ** 2) ** 0.5 - StandardLength).sum()
        loss3 = np.abs(((X_r[1:, :-1] - X_r[:-1, 1:]) ** 2 + (Y_r[1:, :-1] - Y_r[:-1,1:]) ** 2) ** 0.5 - StandardLength * (2 ** 0.5)).sum()
        loss4 = np.abs(((X_r[1:, 1:] - X_r[:-1,:-1]) ** 2 + (Y_r[1:, 1:] - Y_r[:-1,:-1]) ** 2) ** 0.5 - StandardLength * (2 ** 0.5)).sum()

        Loss = loss1 + loss2 + loss3 + loss4
        return Loss

    x0 = -130, -100, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1
    print(make_loss(x0))

    res_distort = minimize(make_loss, x0, method='Powell', options={'maxiter': 1e5, 'disp': True})

    print(res_distort.fun)
    print()
    print(res_distort.x)

其中有个调参的小 trick，我发现将标准长度 StandardLength 设计为结果的均值，最终误差会显著减小。

生成 map

此时我已经获取了精度高于之前模型 6 倍左右的模型了，但当我需要将该反向模型导出为一个 map 用于畸变校正时发现，模型是反的，我该如何得到理想点为整数的畸变点坐标呢？

一个点一个点算当然可以解出来，但是太慢了，于是我使用了一种梯度下降法：

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
# make map

   e_thre = 1e-6

   X_dis_array = np.arange(w)
   Y_dis_array = np.arange(h)
   X_Mesh, Y_Mesh = np.meshgrid(X_dis_array, Y_dis_array)
   X_flatten = X_Mesh.flatten()
   Y_flatten = Y_Mesh.flatten()

   target_X_flatten = X_flatten.copy()
   target_Y_flatten = Y_flatten.copy()

   OK_X_array = np.zeros([0])
   OK_Y_array = np.zeros([0])
   OK_tX_array = np.zeros([0])
   OK_tY_array = np.zeros([0])

   with mt.tqdm(total=len(target_X_flatten)) as pbar:
       while(len(target_X_flatten) >0):

           NX_Mesh, NY_Mesh = make_new_points(res_distort.x, X_flatten, Y_flatten)
           ok_points_mark = ((NX_Mesh - target_X_flatten) ** 2 + (NY_Mesh - target_Y_flatten) ** 2) ** 0.5 <= e_thre

           Cur_OK_X_array = X_flatten[ok_points_mark]
           Cur_OK_Y_array = Y_flatten[ok_points_mark]
           Cur_OK_tX_array = target_X_flatten[ok_points_mark]
           Cur_OK_tY_array = target_Y_flatten[ok_points_mark]

           OK_X_array = np.concatenate([OK_X_array, Cur_OK_X_array])
           OK_Y_array = np.concatenate([OK_Y_array, Cur_OK_Y_array])
           OK_tX_array = np.concatenate([OK_tX_array, Cur_OK_tX_array])
           OK_tY_array = np.concatenate([OK_tY_array, Cur_OK_tY_array])

           X_flatten = X_flatten - (NX_Mesh - target_X_flatten) 
           Y_flatten = Y_flatten - (NY_Mesh - target_Y_flatten) 

           valid_mark = (X_flatten >= 0) * (X_flatten <= w-1) * (Y_flatten >= 0) * (Y_flatten <= h-1)

           next_mark = valid_mark * (~ok_points_mark)

           X_flatten = X_flatten[next_mark]
           Y_flatten = Y_flatten[next_mark]
           target_X_flatten = target_X_flatten[next_mark]
           target_Y_flatten = target_Y_flatten[next_mark]
           pbar.update((~next_mark).sum())

   map_x = np.zeros([h, w], dtype='float32')
   map_y = np.zeros([h, w], dtype='float32')

   map_x[OK_tY_array.astype('int32'), OK_tX_array.astype('int32')] = OK_X_array
   map_y[OK_tY_array.astype('int32'), OK_tX_array.astype('int32')] = OK_Y_array

原始代码有些乱，我解释一下思路：

  1. 生成理想点像素整数坐标阵列 $M_g$
  2. 初始畸变坐标位置等同于理想点阵列 $M_{d0}=M_g$
  3. 将畸变坐标 $M_{di}$ 带入模型，计算得到理想点坐标阵列 $M_{oi}$
  4. 度量残差 $E_i=abs (M_{di} - M_{oi})$
  5. 将残差小于阈值 $\Delta e$ 的点拿走，作为合理的畸变坐标
  6. 剩余的畸变坐标迭代：$M_{di+1}=M_{di}-E_i$
  7. 如果 $M_{di+1}$ 不为空，进入步骤 3

最后收集所有误差满足需求的畸变坐标，整理成 map 就是我们要的矩阵了。

$4096 \times 3072$ 的图像在工作机上计算时间 10s 以内。

叠加校正

经过上述过程获得的畸变校正模型已经优于传统的畸变校正模型结果了，但是我不经意间发现，在传统畸变校正生成的图上再做上述过程可以进一步提高精度。

过程都是工程化实现的东西这里就不详细写了，特别要注意的是此时在做畸变校正时需要校正两次，传统一次、反模型一次，这样很慢很不优雅，那么就需要 remap 叠加 的方法生成统一的 map。

所有上述过程中，可能产生误差的地方需要使用 float 类型，以保留原始数据信息。

实验结果

名称 含义
circle 原始数据
circle_vvd_undis 直接使用反向模型
undis_res_base 使用所有数据的传统模型
undis_res_opencv 使用单张图像的传统模型
undis_res_opencv_vvd_undis 在单张图像传统模型基础上做反向模型












参考资料



文章链接：
https://www.zywvvd.com/notes/study/camera-imaging/undistort-super/undistort-super/

Camera
#Camera
高阶畸变校正
https://www.zywvvd.com/notes/study/camera-imaging/undistort-super/undistort-super/
作者
Yiwei Zhang
发布于
2023年3月12日
许可协议