火星表面...

灿视学长

发布于 2021-05-28 15:02:46

1.7K0

发布于 2021-05-28 15:02:46

文章被收录于专栏：灿视学长灿视学长

各位于晏，亦菲，大家好。我是灿视。

今天，看到了于老师发表了一篇文章，《我为中国火星第一图做鱼眼矫正》。各位可以先去看看于老师的文章，于老师也很大方的开了自己的code。于老师代码写的很简洁，效果也很好。其中代码效果如下所示，展示了火星表面是什么样的：

但是，如果各位的 C++ 基础不是很好的话，可能跑不起来。又或者像我这种不太懂鱼眼曲面相片如何变换成平面图像的原理，就会有点痛苦。

因此，今天当我把模型train 起来之后，我就查了一下鱼眼相片矫正的资料，主要是参考论文 《基于双经度模型的鱼眼图像畸变矫正方法》。现在我相关重点知识分享给各位，希望给各位带来收获。

1. 鱼眼相片简介

鱼眼镜头是一种焦距较短，视角接近甚至超过180°的超广角镜头，因其外形酷似鱼眼而得名。其具有超大视角、信息量丰富、体积小且隐蔽性强等优点而被广泛应用于安全监测、视频会议、机器人导航、全景泊车、智能交通、医疗检测及机器视觉等领域。

然而由于其焦距短视场大的特点以及光学原理约束，导致鱼眼镜头图像桶形畸变较为严重，视觉效果差，很难满足正常需求。如何实现鱼眼图像畸变的精确矫正已成为国内外学者研究的热点。

2. 矫正方案简介

目前，鱼眼图像畸变矫正主要采用两种方式，即基于投影变换模型的矫正方法与和基于标定的鱼眼镜头畸变矫正方法。而又可以将鱼眼镜头投影成像模型分为4 类: 正交投影模型、等距投影模型、球面透视投影模型、等立体角投影模型。

其中基于球面透视投影主要是通过拟合多项式以优化目标函数，从而估计出矫正模型参数，推导出矫正后的图像，该方法研究较早，但计算复杂，实时性差。基于标定的矫正算法，主要是通过借助外部设备对鱼眼图像内外参数进行标定，通过真实坐标与鱼眼成像平面坐标之间坐标转换，实现鱼眼图像畸变矫正，该方法矫正精度高，但对实验设备精度要求较高。

3.双经度畸变矫正方法介绍

传统经纬矫正算法因其无需外部设备标定而被广泛使用，算法基本思想是通过建立鱼眼图像坐标点与半球面模型坐标点之间的映射关系，将球面点坐标转换为球面经纬度坐标，并采用等距投影或正交投影原理将球面上点经度映射为矫正图像水平坐标，纬度映射为矫正图像垂直坐标，从而实现鱼眼镜头的畸变矫正。

其传统经纬斤正算法示意图如图 1 所示, 图中球面上点

m 、 l

和

经度 / 纬度映射为平面图像横向 / 纵向坐标

m^{\prime} / l^{\prime}

和

k^{\prime}

, 其中相同经度上的点

和

, 斤正后具有相同的列坐标

和

k^{\prime}

,同一纬度上的点

和

斤正后具有相同的行坐标

'和

'通常在忽略镜头误差的情况下, 采用等距投影或正交投影建立鱼眼图像点与球面坐标点之间的映射关系, 会导致斤正后图像, 在

轴方向拱形琦变得到一定的斤正, 但

轴方向上拱形琦变很难得到有效矫正, 且靠近顶点位置拉伸更明显。

本文着重介绍基于球面的双经度畸变矫正模型，该模型解决传统经纬矫正算法竖直方向矫正效果好，水平方向仍存在拱形畸变的这一缺点。模型拟通过横向经度和纵向经度对球面进行分割，建立鱼眼图像坐标点与球面双经度坐标的对应关系，然后将纵向和横向经度值映射为平面坐标的横坐标和纵坐标，如图2 所示。

图中显示伴球面上同一横向经度点

m 、 n 、 l

, 经过斤正后映射为平面上的点

m^{\prime} n^{\prime} 、l^{\prime}

具有相同行坐标。如果采用传统经纬斤正算法, 由于

m 、 n 、 l

这3点在球面不同纬度线上, 其对应的映射目标图像上斤正点分别为

m^{\prime \prime} n^{\prime \prime} l^{\prime \prime}

点,可见 3 点对应水平位置基本保持不变, 导致传统经纬斤正算法对水平方向拱形琦变斤正作用下降;度矫正算法后矫正后3 点

m^{\prime}

n^{\prime}

l^{\prime}

位于同一行坐标，对水平方向拱形琦变起到较好的斤正效果, 可见本文所提算法能够同时对鱼眼图像横向与纵向琦变进行矫正。

4. 光学中心和球面半径

的计算

目前有资料证明，根据鱼眼镜头球面模型成像几何性质，空间直线投影为球面上大圆，而球面上大圆投射为鱼眼图像上椭圆，椭圆的长轴的长度与球面模型的直径长度相等。因此，本文将对空间直线在鱼眼图像上对应的畸变曲线进行采样，然后根据采样点拟合椭圆方程，并求得中心位置及长轴大小，从而确定光学中心与球面半径。

假设椭圆方程：

A u^{2}+2 B u v+C v^{2}+D u+E v+F=0

其中，式子:

A,B,C,D,E,F

为椭圆方程参数，其中

B^{2} < AC

，

A,C > 0

，根据资料可以得到椭圆中心

(u_{0}, v_{0})

，长轴的长度

\begin{aligned} u_{0} &=\frac{C D-B E}{2\left(B^{2}-A C\right)} \\ \end{aligned}

v_{0}=\frac{A E-B D}{2\left(B^{2}-A C\right)}

\left.L=\sqrt{\left(\frac{C D^{2}-2 B D E+A E^{2}}{4\left(B^{2}-A C\right)}+F\right)\left(\frac{A+C+\sqrt{(A-C)^{2}+4 B^{2}}}{2\left(B^{2}-A C\right)}\right.}\right)

其中，具体的椭圆参数估计流程如下：

在鱼眼图像上选取真实空间为直线的畸变弧线，并在此弧线上进行采样，，采样点数不低于

个，并记录个点的行列坐标;

采用最小二乘法确定椭圆的系数

A、B、C、D、E、F

。

根据上面公式计算椭圆中心

(u_{0}, v_{0})

，长轴长度

大小。

重复上述步骤，继续采样至少

条弧线，计算出

(u_{0},v_{0}),L

的值。

对上述求取的各

(u_{0},v_{0})，L

求平均值，获取较准确的中心于半径值。

5. 双经度畸变矫正

在图

(3)

中，图

(a)

衣示斤正后的目标图像坐标图,

p^{\prime}

为目标图像上一点, 其坐标为

(i, j) ;

图

(b)

表示半球面双经度俯视图, 点

为与目标图像上点

p^{\prime}

相对应的逆向映射点;图

(c)

为图

(b)

球面模型侧视图,

点与图

(b) p

点相对应, 其坐标为

(x, y, z), p_{1}, p_{2}

分别为

点在

x o z

面和

y o z

面上的投影点,

\varphi 、

分别为

o p_{1}

与

轴正向的夹角和

o p_{2}

与

轴正向的夹角,设

点的坚直方向经度线经度与水平方向经度线经度分别为

\alpha 、 \beta

, 则

\alpha=\pi-\theta

\beta=\pi-\varphi

。

p^{\prime \prime}

点为

点在

x o y

面上的投影点; 图

(d)

表示原鱼眼图像坐标图,

p^{\prime \prime \prime}

对应图

。当投影模型为半球面时, 横向经度与纵向经度取值范围均为

0 \sim \pi

, 直接映射时目标图像太小。为保证目标图像与原图像大小相当, 故以

2 \times R

作为目标图像的行和列数目。因此目标图像上点

p^{\prime}(i, j)

与其对应球面上双经度坐标

(\alpha, \beta)

之间对应关系如下：

\frac{i}{2 R}=\frac{\alpha}{\pi} \Rightarrow \alpha=\left(\frac{\pi}{2 R}\right) \times i

\frac{j}{2 R}=\frac{\beta}{\pi} \Rightarrow \beta=\left(\frac{\pi}{2 R}\right) \times j

式中:

为球面模型的半径，

\alpha

表示双经度坐标中竖直方向经度线经度值，

\beta

为水平方向经度线经度值，

( i，j)

为目标图像上点

p^{\prime}

的坐标值。

由图3(c)可知，

点竖直方向

\alpha

为

op_{2}

与

负半轴的夹角，

点水平方向经度

\beta

为

op_{1}

与

负半轴的夹角，因此可得：

\begin{array}{l} \pi-\theta=\alpha(0 \leqslant \alpha \leqslant \pi) \\ \pi-\varphi=\beta(0 \leqslant \beta \leqslant \pi) \end{array}

继续进行化简，则：

\begin{array}{l} \theta=\pi-\left(\frac{\pi}{2 R}\right) \times i \\ \varphi=\pi-\left(\frac{\pi}{2 R}\right) \times j \end{array}

此外, 根据球面特征以及角度与坐标之间关系, 可建立角度

\varphi 、\theta

与

点坐标值

(x, y, z)

之间的对应关系, 从而得到目标图像点

p^{\prime}

坐标值

(i, j)

与

(x, y, z)

之间的关系式。

\begin{array}{l} x^{2}+y^{2}+z^{2}=R^{2}(z \geqslant 0) \\ \tan (\varphi)=\frac{z}{x}(0 \leqslant \varphi \leqslant \pi) \\ \tan (\theta)=\frac{z}{y}(0 \leqslant \theta \leqslant \pi) \end{array}

联立上面的式子，可得到

(z,y,z)

与

\varphi、\theta

的关系：

\begin{array}{l} x=\frac{R}{\sqrt{\tan ^{2} \varphi+1+\left(\tan ^{2} \varphi\right) /\left(\tan ^{2} \theta\right)}} \\ y=\frac{R}{\sqrt{\tan ^{2} \theta+1+\left(\tan ^{2} \theta\right) /\left(\tan ^{2} \varphi\right)}} \\ z=\frac{R}{\sqrt{1+1 /\left(\tan ^{2} \varphi\right)+1 /\left(\tan ^{2} \theta\right)}} \end{array}

更进一步，可以得到

(x,y,z)

与

(i，j)

之间的关系：

\begin{array}{l} x=\frac{R}{\sqrt{\tan ^{2}\left(\pi-\frac{\pi}{2 R} \times j\right)+1+\tan ^{2}\left(\pi-\frac{\pi}{2 R} \times j\right) /\left(\tan ^{2}\left(\pi-\frac{\pi}{2 R} \times i\right)\right)}} \\ y=\frac{R}{\sqrt{\tan ^{2}\left(\pi-\frac{\pi}{2 R} \times i\right)+1+\tan ^{2}\left(\pi-\frac{\pi}{2 R} \times i\right) /\left(\tan ^{2}\left(\pi-\frac{\pi}{2 R} \times j\right)\right)}} \\ z=\frac{R}{\sqrt{1+1 / \tan ^{2}\left(\pi-\frac{\pi}{2 R} \times j\right)+1 /\left(\tan ^{2}\left(\pi-\frac{\pi}{2 R} \times i\right)\right)}} \end{array}

对于鱼眼图像与球面坐标之间的映射通常采用等距投影或者正交投影2种方式，暂时先不展开整理了。各位可以查阅下这篇论文～

6. 参考代码实例

import cv2
import numpy as np
import math
import time

# 鱼眼有效区域截取
def cut(img):
    img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    (_, thresh) = cv2.threshold(img_gray, 20, 255, cv2.THRESH_BINARY)
    contours, hierarchy = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    cnts = sorted(contours, key=cv2.contourArea, reverse=True)[0]
    x,y,w,h = cv2.boundingRect(cnts)
    r = max(w/ 2, h/ 2)
    # 提取有效区域
    img_valid = img[y:y+h, x:x+w]
    return img_valid, int(r)

# 鱼眼矫正
def undistort(src,r):
    # r： 半径， R: 直径
    R = 2*r
    # Pi: 圆周率
    Pi = np.pi
    # 存储映射结果
    dst = np.zeros((R, R, 3))
    src_h, src_w, _ = src.shape
    # 圆心
    x0, y0 = src_w//2, src_h//2

    for dst_y in range(0, R):

        theta = Pi - (Pi/R)*dst_y
        temp_theta = math.tan(theta)**2

        for dst_x in range(0, R):
            # 取坐标点 p[i][j]
            # 计算 sita 和 fi
            phi = Pi - (Pi/R)*dst_x
            temp_phi = math.tan(phi)**2

            tempu = r/(temp_phi+ 1 + temp_phi/temp_theta)**0.5
            tempv = r/(temp_theta + 1 + temp_theta/temp_phi)**0.5

            if (phi < Pi/2):
                u = x0 + tempu
            else:
                u = x0 - tempu

            if (theta < Pi/2):
                v = y0 + tempv
            else:
                v = y0 - tempv
      # 这里走的是最近邻插值
            if (u>=0 and v>=0 and u+0.5<src_w and v+0.5<src_h):
                dst[dst_y, dst_x, :] = src[int(v+0.5)][int(u+0.5)]

    return dst

if __name__ == "__main__":
    t = time.perf_counter()
    frame = cv2.imread('./imgs/chinese-1st-mars-image.jpeg')
    cut_img,R = cut(frame)
    result_img = undistort(cut_img, R+200) # 这里如果R不加200，
    cv2.imwrite('./imgs/res.jpg',result_img)
    print(time.perf_counter()-t)

该算法得到的效果图如下所示：

与于老师的算法基本上处理方法是类似的，都是针对像素点进行矫正，各位可以详细看看对比中的细节。

为天问一号点赞！

参考资料

https://codechina.csdn.net/mirrors/HLearning/fisheye/-/tree/master/02.double_longitude
https://mp.weixin.qq.com/s/xWFJPlcWUZLGxXD1AS425g

大家好，我是灿视。目前是位算法工程师 + 创业者 + 奶爸的时间管理者！

我曾在19，20年联合了各大厂面试官，连续推出两版《百面计算机视觉》，受到了广泛好评，帮助了数百位同学们斩获了BAT等大小厂算法Offer。现在，我们继续出发，持续更新最强算法面经。

我曾经花了4个月，跨专业从双非上岸华五软工硕士，也从不会编程到进入到百度与腾讯实习。

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2021-05-21，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

本文分享自灿视学长微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

编程算法

登录后参与评论

0 条评论

热度

火星表面...

火星表面...

1. 鱼眼相片简介

2. 矫正方案简介

3.双经度畸变矫正方法介绍

4. 光学中心和球面半径

5. 双经度畸变矫正

6. 参考代码实例

参考资料

社区

活动

资源

关于

腾讯云开发者

热门产品

热门推荐

更多推荐