Su'S Blog

等距圆柱投影法

等距圆柱投影（equidistant cylindrical projection）是一种简单的地图投影方法，在这种投影方法中：假设球面和圆柱面相切于赤道, 将球面上的经纬线投影到圆柱面上, 然后沿圆柱面的一条母线展开成平面的一种投影。

投影坐标换算

假设鱼眼图的大小Hd Wd，将其使用等距圆柱投影法投影到大小为Hs Ws的矩形图内。

在矩形图中建立经纬度笛卡尔坐标系统

代码示例：

ys, xs = np.indices((Hs, Ws), np.float32)
y_proj = Hs / 2.0 - ys
x_proj = xs - Ws / 2.0

以矩形图的中心点为坐标原点，原点右侧为x轴的正方向，上侧为y轴的正方向。x_proj 保存经度坐标，y_proj 保存纬度坐标。

建立鱼眼图到矩形图的坐标变换

鱼眼相机的成像过程分为两步：第一步，三维空间点线性地投影到一个球面上，它是一个虚拟的单位球面，它的球心与相机坐标系的原点重合；第二步，单位球面上的点投影到图像平面上形成鱼眼图，这个过程是非线性的。整个过程如上图所示。

在三维空间坐标系统中，假定坐标系的原点为单位球的球心。y轴穿过单位球的两个极点，单位球经度为0° 的经线与z轴的负半轴相交。已知单位球的半径为1，其球面上点Q的经纬度为(i, j),那么对应的三维坐标为:

由此可得，矩形图中经纬坐标为(m,n)，将其归一化到单位球中其经纬坐标为(m/Ws - 1/2,n/Hs - 1/2)，其三维空间为：

设鱼眼图的视野角范围为fov（取值范围 0 ~ 2 * Pi），则在三维空间坐标为：

化简后可得：

代码示例：

theta_alt = x_proj * fov / Ws
phi_alt = y_proj * np.pi / Hs
x = np.sin(theta_alt) * np.cos(phi_alt)
y = np.sin(phi_alt)
z = np.cos(theta_alt) * np.cos(phi_alt)

由上可知，单位球上的点在笛卡尔坐标下使用 (x,y,z) 表示，对应的球面极坐标表示为（单位球中 r = 1）：

如上图所示，theta为线段OM与z轴正方向的夹角，phi为z轴和点M的半平面与平面zox所构成的角。则由笛卡尔坐标点 (x,y,z) 计算theta和phi的公式如下：

代码示例：

theta = np.arctan2(y, x)
phi = np.arctan2(np.sqrt(x2 + y2), z)

由鱼眼相机的成像过程，在视野角为fov的情况下球面的半径 p = Hd * phi / fov ，则鱼眼图上的笛卡尔坐标点 (x_fish,y_fish) 的极坐标表示为：

将鱼眼图的坐标原点从图形中央移至左上角，则以左上角为坐标原点的坐标为：

代码示例：

p = Hd * phi / fov
y_fish = p * np.sin(theta)
x_fish = p * np.cos(theta)
ymap = Hd / 2.0 - y_fish
xmap = Wd / 2.0 + x_fish

由此，如果矩形图的坐标为 (x,y) ，则对应的鱼眼图的坐标为(xmap,ymap)。

代码示例：

disimg = cv2.remap(fishimg, xmap, ymap, cv2.INTER_LINEAR)

完整代码

    
    fov = fov * np.pi / 180.0
    ys, xs = np.indices((Hs, Ws), np.float32)
    y_proj = Hs / 2.0 - ys
    x_proj = xs - Ws / 2.0
    theta_alt = x_proj * fov / W
    phi_alt = y_proj * np.pi / H
    x = np.sin(theta_alt) * np.cos(phi_alt)
    y = np.sin(phi_alt)
    z = np.cos(theta_alt) * np.cos(phi_alt)
    theta = np.arctan2(y, x)    
    phi = np.arctan2(np.sqrt(x**2 + y**2), z)
    p = Hd * phi / fov
    y_fish = p * np.sin(theta)
    x_fish = p * np.cos(theta)
    ymap = Hd / 2.0 - y_fish
    xmap = Wd / 2.0 + x_fish

    fishimg = cv2.imread("path/example.jpg")
    disimg = cv2.remap(fishimg, xmap, ymap, cv2.INTER_LINEAR)

参考资料：

1、Equirectangular projection：https://en.wikipedia.org/wiki/Equirectangular_projection

2、苏柳.全景播放技术中的帧渲染[J].现代计算机(专业版),2015(23):39-41+54.