Python计算机视觉第五章-多视图几何

本文主要是介绍Python计算机视觉第五章-多视图几何，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

目录

5.1　外极几何

5.1.1一个简单的数据集

5.1.2　用Matplotlib绘制三维数据

5.1.3　计算F：八点法

5.1.4　外极点和外极线

5.2　照相机和三维结构的计算

        5.2.1　三角剖分

5.2.2　由三维点计算照相机矩阵

5.2.3　由基础矩阵计算照相机矩阵

5.3　多视图重建

5.3.1　稳健估计基础矩阵

5.3.2　三维重建示例

5.3.3　多视图的扩展示例

5.4　立体图像

5.1　外极几何

        多视图几何是利用在不同视点所拍摄图像间的关系，来研究照相机之间或者特征之间关系的一门科学。图像的特征通常是兴趣点，本章使用的也是兴趣点特征。多视图几何中最重要的内容是双视图几何。

        如果有一个场景的两个视图以及视图中的对应图像点，那么根据照相机间的空间相对位置关系、照相机的性质以及三维场景点的位置，可以得到对这些图像点的一些几何关系约束。我们通过外极几何来描述这些几何关系。

5.1.1一个简单的数据集

        下面实验是关于使用简单的数据集来测试图像点、三维点和照相机参数矩阵的关系。实验将包括数据生成、相机投影、图像点的回归以及误差分析。

        实验代码：

import numpy as np
import matplotlib.pyplot as plt# 1. 生成三维点
np.random.seed(0)
num_points = 10
points_3D = np.random.rand(num_points, 3) * 10  # 生成10个三维点，范围在0到10之间# 2. 相机内参矩阵 (假设焦距 fx=fy=800, 主点 (cx, cy)=(320, 240) 在640x480的图像分辨率中)
focal_length = 800
cx, cy = 320, 240
camera_matrix = np.array([[focal_length, 0, cx],[0, focal_length, cy],[0, 0, 1]
])# 3. 相机外参矩阵 (旋转矩阵和平移向量)
rotation_matrix = np.eye(3)  # 单位旋转矩阵
translation_vector = np.array([0, 0, -10])  # 沿Z轴平移10个单位# 4. 计算图像点
def project_points(points_3D, camera_matrix, rotation_matrix, translation_vector):points_3D_homogeneous = np.hstack((points_3D, np.ones((points_3D.shape[0], 1))))extrinsic_matrix = np.hstack((rotation_matrix, translation_vector.reshape(-1, 1)))points_2D_homogeneous = camera_matrix @ (extrinsic_matrix @ points_3D_homogeneous.T)points_2D_homogeneous /= points_2D_homogeneous[2, :]  # 归一化return points_2D_homogeneous[:2, :].Tpoints_2D = project_points(points_3D, camera_matrix, rotation_matrix, translation_vector)# 5. 可视化图像点
plt.scatter(points_2D[:, 0], points_2D[:, 1], c='red', label='Projected Points')
plt.title('Projected 2D Points')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.gca().invert_yaxis()  # 图像坐标系Y轴反向
plt.show()# 6. 计算重投影误差
def compute_reprojection_error(points_3D, points_2D, camera_matrix, rotation_matrix, translation_vector):points_2D_projected = project_points(points_3D, camera_matrix, rotation_matrix, translation_vector)error = np.linalg.norm(points_2D - points_2D_projected, axis=1)return error.mean()reprojection_error = compute_reprojection_error(points_3D, points_2D, camera_matrix, rotation_matrix, translation_vector)
print(f'Reprojection Error: {reprojection_error:.4f}')

        分析：

1. 数据生成：我们生成了一个包含10个三维点的点云，并设置了一个相机内参矩阵。相机外参矩阵由单位旋转矩阵和平移向量组成，将三维点投影到图像平面上。

2. 图像点投影：通过相机内参矩阵和外参矩阵将三维点投影到二维平面上。实验结果显示了这些投影点在图像上的分布。

3. 重投影误差：计算了投影后的二维点和理论二维点之间的重投影误差。由于我们使用了完全正确的外参矩阵，重投影误差应该很小

结果：

• 图像点：图像点的分布反映了三维点在相机坐标系下的位置和投影方式。
• 重投影误差：在理论上，因为我们使用了准确的参数，重投影误差应该非常小。实际应用中，误差可以用于评估相机标定和模型的准确性。

5.1.2　用Matplotlib绘制三维数据

以下是使用Matplotlib绘制三维数据的实验代码：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D# 1. 生成三维数据
np.random.seed(0)
num_points = 100
x = np.random.uniform(-10, 10, num_points)
y = np.random.uniform(-10, 10, num_points)
z = np.random.uniform(-10, 10, num_points)# 2. 创建一个三维图
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')# 3. 绘制三维散点图
sc = ax.scatter(x, y, z, c='blue', marker='o')# 4. 设置图形标题和标签
ax.set_title('3D Scatter Plot of Random Points')
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
ax.set_zlabel('Z axis')# 5. 显示图形
plt.show()

        分析：

1. 数据生成：生成了100个随机三维点，每个点的x、y、z坐标在-10到10之间均匀分布。

2. 绘图：

   • 使用matplotlib的Axes3D来创建一个三维坐标系。
   • 使用scatter函数绘制三维散点图，点的颜色为蓝色，标记样式为圆形。

3. 图形解释：

   • 三维散点图：图中每个点表示一个三维坐标位置。通过旋转和缩放图形，可以观察数据的分布情况。
   • 坐标轴标签：X、Y、Z轴标签帮助我们理解三维数据在空间中的分布。

结果：

• 数据分布：在图中，点的分布是随机的，没有明显的聚集或规律。这种分布是由生成数据时使用的均匀分布所决定的。
• 图形直观性：三维散点图能够帮助我们直观地理解数据的空间结构，特别是在处理实际数据集时，这种可视化方式可以揭示数据的潜在结构和关系。

5.1.3　计算F：八点法

        八点法是通过对应点来计算基础矩阵的算法。

下面是使用Python计算基础矩阵的代码，基于八点法：

import numpy as np
import cv2def normalize_points(points):"""Normalize the points to improve numerical stability."""mean = np.mean(points, axis=0)std = np.std(points, axis=0)T = np.array([[1/std[0], 0, -mean[0]/std[0]],[0, 1/std[1], -mean[1]/std[1]],[0, 0, 1]])points_normalized = np.dot(T, np.vstack((points.T, np.ones(points.shape[0]))))return points_normalized[:2].T, Tdef compute_fundamental_matrix(pts1, pts2):"""Compute the fundamental matrix using the eight-point algorithm."""# Normalize pointspts1, T1 = normalize_points(pts1)pts2, T2 = normalize_points(pts2)# Construct the matrix AA = np.vstack([pts1[:, 0] * pts2[:, 0],pts1[:, 0] * pts2[:, 1],pts1[:, 0],pts1[:, 1] * pts2[:, 0],pts1[:, 1] * pts2[:, 1],pts1[:, 1],pts2[:, 0],pts2[:, 1],np.ones_like(pts2[:, 0])]).T# Compute the fundamental matrix F_, _, V = np.linalg.svd(A)F = V[-1].reshape(3, 3)# Enforce the rank-2 constraintU, S, Vt = np.linalg.svd(F)S[2] = 0F = np.dot(U, np.dot(np.diag(S), Vt))# Denormalize FF = np.dot(T2.T, np.dot(F, T1))F = F / F[2, 2]return F# Example usage with random points (replace with actual data)
pts1 = np.random.rand(8, 2) * 1000
pts2 = np.random.rand(8, 2) * 1000F = compute_fundamental_matrix(pts1, pts2)
print("Fundamental Matrix:\n", F)

分析：

1. 数据生成：示例中使用了随机生成的8个点对作为输入数据。实际应用中，应使用实际图像中的匹配点对。

2. 点归一化：点归一化有助于提高计算的数值稳定性。函数normalize_points对点进行均值和标准差归一化。

3. 基础矩阵计算：

   • 矩阵构造：构造了包含点对应关系的矩阵A。
   • SVD分解：通过SVD计算基础矩阵F，并应用rank-2约束使其符合要求。
   • 反归一化：将计算得到的基础矩阵F从归一化坐标系统转换回原始坐标系统。

结果：

• 基础矩阵：计算出的基础矩阵F用于描述两幅图像中的点对应关系。它应符合图像点对之间的几何约束。
• 准确性：计算的基础矩阵F的准确性依赖于输入点对的质量。在实际应用中，应选择精确的点匹配以获得更好的结果。

5.1.4　外极点和外极线

• 外极点 (Epipole): 在立体视觉中，外极点是一个图像中对应于另一个图像的相机中心的投影点。具体来说，若有两幅图像I1和I2，并且F是它们的基础矩阵，则I1中的外极点e1是I2中相机中心在I1中的对应点，反之亦然。

• 外极线 (Epipolar Line): 外极线是对应于一个图像中点的直线，它是另一幅图像中所有可能对应点的位置所在的直线。若e2是I2中的外极点，那么在I1中与I2中的一个点p2对应的外极线是通过e1并且与p2匹配的直线。

        下面是关于使用Python实现外极点和外极线的可视化。分析结果并理解外极点和外极线的几何关系的实验：

import numpy as np
import matplotlib.pyplot as plt
import cv2def compute_epipoles(F):"""Compute the epipoles from the fundamental matrix."""U, _, Vt = np.linalg.svd(F)e1 = Vt[-1]e2 = U[:, -1]e1 /= e1[2]  # Normalizee2 /= e2[2]  # Normalizereturn e1, e2def compute_epipolar_lines(F, points, epipole):"""Compute epipolar lines in the image for given points."""lines = np.dot(F.T, np.hstack((points, np.ones((points.shape[0], 1)))).T).Tlines = lines / np.sqrt(lines[:, 0]**2 + lines[:, 1]**2)[:, np.newaxis]return linesdef plot_epipolar_lines(img, lines, points):"""Plot epipolar lines on the image."""plt.imshow(img, cmap='gray')for line, point in zip(lines, points):x0, y0 = 0, -line[2] / line[1]x1, y1 = img.shape[1], -(line[2] + line[0] * img.shape[1]) / line[1]plt.plot([x0, x1], [y0, y1], 'r-')plt.plot(point[0], point[1], 'bo')plt.show()# Example usage with synthetic data
F = np.array([[1e-6, 0, -1e-3], [0, 1e-6, -1e-3], [1e-3, 1e-3, 1]])# Calculate epipoles
e1, e2 = compute_epipoles(F)
print("Epipole in image 1:", e1)
print("Epipole in image 2:", e2)# Synthetic points and image
points = np.array([[100, 200], [300, 400], [500, 600]])
img = np.ones((800, 800))# Compute epipolar lines
lines = compute_epipolar_lines(F, points, e1)# Plot
plot_epipolar_lines(img, lines, points)

        分析：

1. 外极点计算：

   • 使用compute_epipoles函数从基础矩阵F中计算两个图像的外极点。计算得到的外极点e1和e2应符合几何约束，即它们是另一个图像的相机中心在当前图像中的投影。

2. 外极线计算：

   • compute_epipolar_lines函数计算每个点在图像中的对应外极线。外极线是通过F和图像点的关系得到的。

3. 结果可视化：

   • 使用plot_epipolar_lines函数可视化图像中的外极线和点。外极线的绘制帮助理解点在另一幅图像中的匹配区域。

4. 准确性：

   • 外极点和外极线的准确性直接取决于基础矩阵的计算精度。在实际应用中，通常使用更准确的数据和方法（如RANSAC）来提高计算结果的质量。

结果：

5.2　照相机和三维结构的计算

        5.2.1　三角剖分

        三角剖分是将一个复杂的平面区域分解成若干个不重叠的三角形的过程，常用于计算机图形学、地理信息系统和数值分析中。它有助于简化复杂形状的处理，并为后续的计算提供支持。

        下面实验实现一个三角剖分算法，并通过该算法恢复出给定点集的三维位置：

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay
from mpl_toolkits.mplot3d import Axes3D# 生成随机二维点
np.random.seed(0)
points_2d = np.random.rand(30, 2)# 执行Delaunay三角剖分
tri = Delaunay(points_2d)# 可视化二维点和三角剖分结果
plt.figure()
plt.triplot(points_2d[:,0], points_2d[:,1], tri.simplices, lw=0.5)
plt.plot(points_2d[:,0], points_2d[:,1], 'o')
plt.title('2D Delaunay Triangulation')
plt.show()# 假设二维点的z坐标为0，恢复三维位置
points_3d = np.hstack((points_2d, np.zeros((points_2d.shape[0], 1))))# 可视化三维点（示例）
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points_3d[:,0], points_3d[:,1], points_3d[:,2])
ax.set_title('3D Points')
plt.show()# 打印三角剖分的顶点和三角形索引
print("Delaunay Triangulation:\n", tri.simplices)

分析：

1. 二维可视化：图中显示了随机生成的点及其Delaunay三角剖分的结果。各三角形表示了点集的局部结构。
2. 三维重建：在此实验中，我们假设二维点的z坐标为0，三维点只是二维点的扩展。因此，所有点都位于z=0的平面上。

结果：

        三角剖分有效地将点集分解为三角形，并在二维平面中提供了良好的可视化。通过这种剖分，可以在实际应用中进行更复杂的三维重建和分析。

5.2.2　由三维点计算照相机矩阵

        相机矩阵（Camera Matrix），在计算机视觉中，通常指的是相机内参矩阵（IntrinsicMatrix）和外参矩阵（Extrinsic Matrix）的组合。在这里，我们主要关注从三维点计算相机外参矩阵（包括旋转矩阵和平移向量）的过程。

import numpy as np
import cv2# 已知的三维点（世界坐标系中的点）
object_points = np.array([[0, 0, 0],     # 点1[1, 0, 0],     # 点2[1, 1, 0],     # 点3[0, 1, 0],     # 点4[0, 0, 1],     # 点5[1, 0, 1],     # 点6[1, 1, 1],     # 点7[0, 1, 1]      # 点8
], dtype=np.float32)# 对应的二维图像点
image_points = np.array([[100, 150],    # 点1[200, 150],    # 点2[200, 250],    # 点3[100, 250],    # 点4[120, 180],    # 点5[220, 180],    # 点6[220, 280],    # 点7[120, 280]     # 点8
], dtype=np.float32)# 相机内参矩阵（假设）
camera_matrix = np.array([[800, 0, 320],[0, 800, 240],[0, 0, 1]
], dtype=np.float32)# 相机的畸变系数（假设无畸变）
dist_coeffs = np.zeros((4, 1), dtype=np.float32)# 计算外参矩阵
success, rotation_vector, translation_vector = cv2.solvePnP(object_points, image_points, camera_matrix, dist_coeffs)# 计算旋转矩阵
rotation_matrix, _ = cv2.Rodrigues(rotation_vector)# 构造相机矩阵 [R | t]
camera_matrix_ext = np.hstack((rotation_matrix, translation_vector))print("旋转矩阵:\n", rotation_matrix)
print("平移向量:\n", translation_vector)
print("外参矩阵 [R | t]:\n", camera_matrix_ext)

分析：

1. 旋转矩阵 (Rotation Matrix): 表示相机坐标系相对于世界坐标系的旋转。
2. 平移向量 (Translation Vector): 表示相机坐标系相对于世界坐标系的平移。
3. 外参矩阵: 是一个 3x4 矩阵，结合了旋转矩阵和平移向量，表示了相机的外部参数。

结果：

• 旋转矩阵: 3x3 的矩阵，用于将三维点从世界坐标系转换到相机坐标系的方向。
• 平移向量: 一个 3x1 的向量，表示相机在世界坐标系中的位置。

5.2.3　由基础矩阵计算照相机矩阵

        在计算机视觉中，基础矩阵是描述两个相机视图中点之间几何关系的一个重要矩阵。给定两个视图的匹配点对，我们可以用基础矩阵来计算从一个视图到另一个视图的投影关系。基础矩阵 FF 是一个 3×33×3 的矩阵，它可以通过以下公式将一个点在第一个图像中的坐标映射到第二个图像中的坐标：

                                                                x′TFx=0

下面是一个简单的实验代码示例，展示如何计算基础矩阵，并从中估计相机矩阵。

import numpy as np
import cv2# 示例数据（替换为你自己的数据）
points1 = np.array([[100, 150],[200, 250],[300, 350],[400, 450]
], dtype=np.float32)points2 = np.array([[110, 160],[210, 260],[310, 360],[410, 460]
], dtype=np.float32)# 检查点集数量
if len(points1) < 8 or len(points2) < 8:print("错误：需要至少 8 对匹配点。")
else:# 计算基础矩阵F, mask = cv2.findFundamentalMat(points1, points2, method=cv2.FM_8POINT)# 检查基础矩阵 F 是否为 Noneif F is None:print("计算基础矩阵失败，返回值为 None。请检查输入点是否足够或匹配是否正确。")else:# 输出基础矩阵以进行调试print("基础矩阵 F:\n", F)print("基础矩阵 F 的维度:", F.shape)# 相机内参矩阵 (假设相机内参已知)K1 = np.array([[1000, 0, 320],[0, 1000, 240],[0, 0, 1]])K2 = K1  # 假设两个相机的内参矩阵相同# 计算本质矩阵E = K2.T @ F @ K1# 输出本质矩阵以进行调试print("本质矩阵 E:\n", E)print("本质矩阵 E 的维度:", E.shape)# 分解本质矩阵_, R, t, _ = cv2.recoverPose(E, points1, points2, K1)# 输出结果print("旋转矩阵 R:\n", R)print("平移向量 t:\n", t)

分析：

1. 基础矩阵 FF：计算出的基础矩阵 FF 描述了两个图像平面之间的几何关系。

2. 本质矩阵 EE：本质矩阵通过基础矩阵 FF 和相机内参矩阵 KK 计算得到，它包含了旋转和位移的信息。

3. 旋转矩阵 RR 和 平移向量 tt：这两个矩阵描述了从一个视图到另一个视图的相机姿态变化。旋转矩阵和位移向量共同定义了相机的相对位置和方向。

结果：

        通过计算基础矩阵并从中估计相机矩阵，我们能够理解和重建两个视图之间的空间关系。这对于立体视觉中的图像对齐、三维重建等任务具有重要意义。此实验演示了如何使用OpenCV来实现这些计算，并提供了一种基础的方法来处理实际的视觉数据。

5.3　多视图重建

        假设照相机已经标定，计算重建可以分为下面 4 个步骤：

import numpy as np
import cv2
import matplotlib.pyplot as plt# 生成一些示例数据（实际应用中你应该用自己的匹配点）
np.random.seed(0)
points1 = np.random.rand(100, 2) * 1000
points2 = points1 + (np.random.rand(100, 2) - 0.5) * 20  # 加入一些随机噪声
# 随机加入离群点
outliers = np.random.rand(20, 2) * 1000
points2[:20] = outliers# 计算基础矩阵并进行 RANSAC 稳健估计
F, mask = cv2.findFundamentalMat(points1, points2, method=cv2.FM_RANSAC)# 过滤出内点
inlier_points1 = points1[mask.ravel() == 1]
inlier_points2 = points2[mask.ravel() == 1]# 计算内点和离群点数量
num_inliers = np.sum(mask)
num_outliers = len(mask) - num_inliers# 绘图
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.scatter(points1[:, 0], points1[:, 1], c='r', label='All Points')
plt.scatter(points2[:, 0], points2[:, 1], c='b', label='All Points')
plt.title('Original Points')
plt.legend()plt.subplot(1, 2, 2)
plt.scatter(inlier_points1[:, 0], inlier_points1[:, 1], c='g', label='Inliers')
plt.scatter(inlier_points2[:, 0], inlier_points2[:, 1], c='y', label='Inliers')
plt.title('Inliers')
plt.legend()plt.show()# 输出结果
print(f"计算的基础矩阵:\n{F}")
print(f"内点数量: {num_inliers}")
print(f"离群点数量: {num_outliers}")

import numpy as np
import cv2
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D# 读取图像
img1 = cv2.imread('1.jpg', cv2.IMREAD_GRAYSCALE)
img2 = cv2.imread('1_blurred.jpg', cv2.IMREAD_GRAYSCALE)# 检测特征点和计算描述符
orb = cv2.ORB_create()
keypoints1, descriptors1 = orb.detectAndCompute(img1, None)
keypoints2, descriptors2 = orb.detectAndCompute(img2, None)# 匹配特征点
bf = cv2.BFMatcher(cv2.NORM_HAMMING, crossCheck=True)
matches = bf.match(descriptors1, descriptors2)# 提取匹配点
points1 = np.float32([keypoints1[m.queryIdx].pt for m in matches])
points2 = np.float32([keypoints2[m.trainIdx].pt for m in matches])# 相机内参数矩阵（示例值，请根据实际情况调整）
K = np.array([[1000, 0, img1.shape[1] / 2],[0, 1000, img1.shape[0] / 2],[0, 0, 1]])# 计算基础矩阵
F, _ = cv2.findFundamentalMat(points1, points2, cv2.FM_8POINT)# 从基础矩阵计算本质矩阵
E = K.T @ F @ K# 相机矩阵
_, R, t, _ = cv2.recoverPose(E, points1, points2, K)# 三角测量
points4D = cv2.triangulatePoints(np.eye(3, 4),np.hstack((R, t)),points1.T,points2.T
)
points3D = points4D[:3] / points4D[3]# 绘制三维点云
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points3D[0], points3D[1], points3D[2], s=1)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()# 输出三维点云的基本信息
print(f"重建的三维点数: {points3D.shape[1]}")

这篇关于Python计算机视觉第五章-多视图几何的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！