如何根据点对应计算球形相机位置？

Question

我在等距矩形图像中标记了 4 个点。[红点]

我还在俯视图像中标记了 4 个对应点 [红点]

如何计算摄像机在俯视图像上的位置？

到目前为止，我看到有 4 条射线 ( R1, R2, R3, R4) 从未知的相机中心延伸C = (Cx, Cy, Cz)穿过等距柱状图像中的点，并终止于俯视图像的像素坐标 ( P1, P2, P3, P4)。因此 4 个向量方程的形式为：

[Cx, Cy, Cz] + [Rx, Ry, Rz]*t = [x, y, 0] 

对于每封信件。所以

C + R1*t1 = P1 = [x1, y1, 0]
C + R2*t2 = P2 = [x2, y2, 0]
C + R3*t3 = P3 = [x3, y3, 0]
C + R4*t4 = P4 = [x4, y4, 0]

那么 7 个未知数和 12 个方程？这是我的尝试，但似乎没有给出合理的答案：

import numpy as np

def equi2sphere(x, y):
    width = 2000
    height = 1000
    theta = 2 * np.pi * x  / width - np.pi
    phi = np.pi * y / height
    return theta, phi

HEIGHT = 1000
MAP_HEIGHT = 788
#
# HEIGHT = 0
# MAP_HEIGHT = 0

# Point in equirectangular image, bottom left = (0, 0)
xs = [1190, 1325, 1178, 1333]
ys = [HEIGHT - 730,   HEIGHT - 730,  HEIGHT - 756,  HEIGHT - 760]

# import cv2
# img = cv2.imread('equirectangular.jpg')
# for x, y in zip(xs, ys):
#     img = cv2.circle(img, (x, y), 15, (255, 0, 0), -1)
# cv2.imwrite("debug_equirectangular.png", img)

# Corresponding points in overhead map, bottom left = (0, 0)
px = [269, 382, 269, 383]
py = [778, 778, 736, 737]

# import cv2
# img = cv2.imread('map.png')
# for x, y in zip(px, py):
#     img = cv2.circle(img, (x, y), 15, (255, 0, 0), -1)
# cv2.imwrite("debug_map.png", img)

As = []
bs = []
for i in range(4):

    x, y = xs[i], ys[i]

    theta, phi = equi2sphere(x, y)

    # convert to spherical
    p = 1
    sx = p * np.sin(phi) * np.cos(theta)
    sy = p * np.sin(phi) * np.sin(theta)
    sz = p * np.cos(phi)

    print(x, y, '->', np.degrees(theta), np.degrees(phi), '->', round(sx, 2), round(sy, 2), round(sz, 2))

    block = np.array([
        [1, 0, 0, sx],
        [0, 1, 0, sy],
        [1, 0, 1, sz],
    ])

    y = np.array([px[i], py[i], 0])

    As.append(block)
    bs.append(y)



A = np.vstack(As)
b = np.hstack(bs).T
solution = np.linalg.lstsq(A, b)
Cx, Cy, Cz, t = solution[0]

import cv2
img = cv2.imread('map_overhead.png')

for i in range(4):

    x, y = xs[i], ys[i]

    theta, phi = equi2sphere(x, y)

    # convert to spherical
    p = 1
    sx = p * np.sin(phi) * np.cos(theta)
    sy = p * np.sin(phi) * np.sin(theta)
    sz = p * np.cos(phi)

    pixel_x = Cx + sx * t
    pixel_y = Cy + sy * t
    pixel_z = Cz + sz * t
    print(pixel_x, pixel_y, pixel_z)
    img = cv2.circle(img, (int(pixel_x), img.shape[0] - int(pixel_y)), 15, (255,255, 0), -1)
    img = cv2.circle(img, (int(Cx), img.shape[0] - int(Cy)), 15, (0,255, 0), -1)


cv2.imwrite("solution.png", img)


# print(A.dot(solution[0]))
# print(b)

生成的相机位置（绿色）和投影点（青色）

编辑：修复的一个错误是 PI/4 中等距柱状图像中的经度偏移，它修复了旋转问题，但比例仍然以某种方式关闭。

Answer 1

编辑：使用 MAP 图片宽度/长度进行球面转换可以为相机中心提供更好的结果。积分位置还是有点乱。地图上有更好的相机中心解决方案：，点有些扁平化

我冒昧地重写了一些代码，使用变量和颜色添加点标识（在您的原始代码中，某些点在各个列表中的顺序不同）。如果想要处理更多数据点，这是更好的选择。是的，我选择了一个字典用于调试目的，但是 N 个点的列表确实是更好的选择，前提是它们在不同的投影之间正确索引配对。

我还调整了坐标以匹配我的图片。以及 x,y 变量的用法/命名以供我理解。

它仍然是不正确的，但是每个找到的位置之间存在某种一致性。

可能的原因

OpenCV 图像将 [0,0] 放在左上角。下面的代码与点坐标的约定一致，但我没有更改任何数学公式。

也许某些公式存在错误或不一致。您可能需要再次检查您的约定：符号、[0,0] 位置等。

我在公式中没有看到任何与相机位置和高度相关的输入，这可能是错误的来源。

您可以查看这个执行等距柱状投影的项目：https://github.com/NitishMutha/equirectangular-toolbox

from typing import Dict

import cv2
import numpy as np


def equi2sphere(x, y, width, height):
    theta = (2 * np.pi * x / width) - np.pi
    phi = (np.pi * y / height)
    return theta, phi


WIDTH = 805
HEIGHT = 374  # using stackoverflow PNG
MAP_WIDTH = 662
MAP_HEIGHT = 1056  # using stackoverflow PNG

BLUE = (255, 0, 0)
GREEN = (0, 255, 0)
RED = (0, 0, 255)
CYAN = (255, 255, 0)
points_colors = [BLUE, GREEN, RED, CYAN]

TOP_LEFT = "TOP_LEFT"
TOP_RIGHT = "TOP_RIGHT"
BOTTOM_LEFT = "BOTTOM_LEFT"
BOTTOM_RIGHT = "BOTTOM_RIGHT"


class Point:
    def __init__(self, x, y, color):
        self.x = x
        self.y = y
        self.c = color

    @property
    def coords(self):
        return (self.x, self.y)


# coords using GIMP which uses upperleft [0,0]
POINTS_ON_SPHERICAL_MAP: Dict[str, Point] = {TOP_LEFT    : Point(480, 263, BLUE),
                                             TOP_RIGHT   : Point(532, 265, GREEN),
                                             BOTTOM_LEFT : Point(473, 274, RED),
                                             BOTTOM_RIGHT: Point(535, 275, CYAN),
                                             }
# xs = [480, 532, 473, 535, ]
# ys = [263, 265, 274, 275, ]

img = cv2.imread('equirectangular.png')
for p in POINTS_ON_SPHERICAL_MAP.values():
    img = cv2.circle(img, p.coords, 5, p.c, -1)
cv2.imwrite("debug_equirectangular.png", img)

# coords using GIMP which uses upperleft [0,0]
# px = [269, 382, 269, 383]
# py = [278, 278, 320, 319]
POINTS_ON_OVERHEAD_MAP: Dict[str, Point] = {TOP_LEFT    : Point(269, 278, BLUE),
                                            TOP_RIGHT   : Point(382, 278, GREEN),
                                            BOTTOM_LEFT : Point(269, 320, RED),
                                            BOTTOM_RIGHT: Point(383, 319, CYAN),
                                            }

img = cv2.imread('map.png')

for p in POINTS_ON_OVERHEAD_MAP.values():
    img = cv2.circle(img, p.coords, 5, p.c, -1)
cv2.imwrite("debug_map.png", img)

As = []
bs = []
for point_location in [TOP_LEFT, TOP_RIGHT, BOTTOM_LEFT, BOTTOM_RIGHT]:
    x_spherical, y_spherical = POINTS_ON_SPHERICAL_MAP[point_location].coords
    theta, phi = equi2sphere(x=x_spherical, y=y_spherical, width=MAP_WIDTH, height=MAP_HEIGHT) # using the overhead map data for conversions
    # convert to spherical
    p = 1
    sx = p * np.sin(phi) * np.cos(theta)
    sy = p * np.sin(phi) * np.sin(theta)
    sz = p * np.cos(phi)
    print(f"{x_spherical}, {y_spherical} -> {np.degrees(theta):+.3f}, {np.degrees(phi):+.3f} -> {sx:+.3f}, {sy:+.3f}, {sz:+.3f}")
    block = np.array([[1, 0, 0, sx],
                      [0, 1, 0, sy],
                      [1, 0, 1, sz], ])
    x_map, y_map = POINTS_ON_OVERHEAD_MAP[point_location].coords
    vector = np.array([x_map, y_map, 0])
    As.append(block)
    bs.append(vector)

A = np.vstack(As)
b = np.hstack(bs).T
solution = np.linalg.lstsq(A, b)
Cx, Cy, Cz, t = solution[0]

img = cv2.imread("debug_map.png")

for point_location in [TOP_LEFT, TOP_RIGHT, BOTTOM_LEFT, BOTTOM_RIGHT]:
    x_spherical, y_spherical = POINTS_ON_SPHERICAL_MAP[point_location].coords
    theta, phi = equi2sphere(x=x_spherical, y=y_spherical, width=MAP_WIDTH, height=MAP_HEIGHT) # using the overhead map data for conversions
    # convert to spherical
    p = 1
    sx = p * np.sin(phi) * np.cos(theta)
    sy = p * np.sin(phi) * np.sin(theta)
    sz = p * np.cos(phi)
    pixel_x = Cx + sx * t
    pixel_y = Cy + sy * t
    pixel_z = Cz + sz * t
    print(f"{pixel_x:+0.0f}, {pixel_y:+0.0f}, {pixel_z:+0.0f}")
    img = cv2.circle(img, (int(pixel_x), int(pixel_y)), 5, POINTS_ON_SPHERICAL_MAP[point_location].c, -1)

img = cv2.circle(img, (int(Cx), int(Cy)), 4, (200, 200, 127), 3)

cv2.imwrite("solution.png", img)

与我最初的解决方案进行映射： 调试图： 等距矩形图像： 调试等距矩形：